森林是全球陆地生态系统中的最大有机碳库^[1], 森林生态系统贮存了全球陆地植被碳库的77%, 全球土壤碳库的39%^[2], 在维护区域生态环境和全球碳平衡中起着巨大的作用。准确估算区域森林碳储量及其分布, 对于评价森林生态系统的功能具有重要意义。但是, 由于方法和数据的缺乏, 不同学者估测的森林碳储量相差较大, 导致森林碳汇功能评价具有较大的不确定性^[3]。

日益发展的遥感技术具有快速、准确、对森林无破坏性并能进行宏观监测的优势, 使得遥感成为获取森林地上生物量的主要途径。目前基于遥感数据计算森林碳储量的方法主要有：回归估计法^[4-5], 神经网络法^[6-7], 遥感数据与过程模型融合的方法^[8]和空间仿真模拟^[9-10]等方法。与其余方法相比, 神经网络与空间仿真方法具有较高的估算精度^[11-12]。人工神经网络 (ANN) 是通过模拟人脑神经系统建立起来的一类模型, 在建模时不需要给出具体的数学函数, 可以一次性引入多个解释变量, 并同时输出多个估测量, 适合复杂的非线性模型的模拟, 被广泛地应用于工学、天文学、生物学等领域。近年来, 神经网络也常被用来估算森林生物量。与其他模型相比神经网络方法虽然可以提高森林的预测精度^[13-14], 但是在模拟预测过程中易陷入局部最优, 从而导致了部分区域模拟结果与实际值不相符^[15]。而且在模拟过程中, 由于隐含层个数和单元数的确定缺乏理论指导, 学习与记忆具有不确定性^[16], 虽然使用相同方法、相同数据但是却得出不同的结果, 使研究结果不能重复。

空间仿真模拟采用与地统计学方法相结合的随机算法模拟森林碳分布, 通过对局部森林碳分布特征量的分析, 得到局部森林碳分布函数, 再用蒙特卡洛方法实现对局部的估计^[10]。它不像克里格方法, 追求的是特定点位某个属性的局部最优估值^[17], 也不像回归方法那样只考虑保证总体平均数的估计精度。空间仿真模拟追求的是森林碳储量模拟的真实性, 尽可能地接近真实的空间分布^[9]。由于不同的方法对于同一地区碳储量的估算结果是不同的, 因此对同一地区使用不同方法进行碳储量估算, 可以对比各方法的优劣, 为森林碳估算提供指导。本文以龙泉市为研究对象, 采用BP神经网络与空间仿真模拟方法对其地上部分碳储量和碳分布进行仿真, 并对两种方法的估算能力进行对比分析。

1 研究区与数据 1.1 研究区概况

龙泉市 (118°43′—119°26′E, 27°42′—28°21′N) 地处浙江省西南部, 隶属丽水市, 东西长68.9 km, 南北宽70.8 km, 总面积3059 km²。

龙泉市在地貌上属于浙南山地, 地形复杂、海拔高低悬殊, 因此气候基本呈垂直变化分布, 光、温、水地域差异明显。该市地处亚热带季风气候区, 温暖湿润, 四季分明, 雨水充沛, 光热较优, 适宜各种林木生长, 植物资源丰富。全市共有高等植物1800余种, 其中木本植物1105种 (含种下分类群), 占全省3/4以上。

在全国森林资源经营管理分区方案中, 龙泉市属于南方山地丘陵区中的南方低山丘陵亚区, 是重要的集体林区。全市林业用地面积265633 hm², 占总面积87.17%, 森林蓄积量1 455.9万m³, 森林覆盖率84.2%, 乔木林年总生长量为101.9万m³, 生长率8.53%(2008年)。森林类型主要有常绿落叶阔叶混交林、针叶阔叶混交林、常绿阔叶林、黄山松林、马尾松林、杉木林、毛竹林以及山地矮林、灌丛等类型。

1.2 研究数据预处理 1.2.1 地面样地数据

研究区地面调查数据为2009年龙泉市森林资源连续清查样地调查数据, 全市共有按系统抽样方法布设的固定样地102个, 样地间距6 km×4 km, 样地形状为正方形, 样地面积0.08 hm²。本研究采用的有效样地个数为99个, 样地数据特征值见表 1。由于现有的森林生物量模型有限, 本研究将树种分为杉木、马尾松、硬阔和软阔4个树种组, 根据已发表的生物量模型^[18-20]进行样地内单株地上部分生物量计算。如果样地含有毛竹, 毛竹生物量根据文献[21]的毛竹单株生物量模型计算。单木碳储量由单木生物量乘以碳储量转换系数0.5^[22]得到, 最后累加求得样地碳储量。结合样地碳储量与遥感因子利用神经网络和空间仿真方法来估计研究区森林地上部分碳储量及其碳密度。

1.2.2 遥感数据

遥感数据选用2009年龙泉市全境Landsat TM影像数据, 由119/40和119/41两景组成。并对其进行了几何校正和辐射校正, 总误差小于一个像元。对于用于建模的遥感因子, 通过ArcGIS提取样地所对应的遥感图像6个波段灰度值、相关植被指数以及其它波段组合等遥感变量。通过分析比较各样地碳密度与对应遥感变量之间的相关性, 选取其中与碳密度相关性较大的遥感变量进行碳储量的估算。主要遥感变量与森林碳密度之间的关系见表 2。

2 研究方法 2.1 BP神经网络仿真建模

BP神经网络是1种多层前馈神经网络, 该网络的主要特点是信号前向传递, 误差反向传递。主要由3部分组成：输入层、隐含层、输出层。隐含层可以分为一层或多层, 一个包含两层隐含层的BP神经网络的拓扑结构如图 1所示^[23]。

BP网络模型处理信息的基本原理是：输入信号通过隐层节点作用于输出层节点, 经过非线形变换, 产生输出信号, 网络训练的每个样本包括输入向量P和期望输出量Y, 网络输出值A与期望输出值Y之间的偏差, 通过调整输入节点与隐层节点的联接强度取值和隐层节点与输出节点之间的联接强度以及阈值, 使误差沿梯度方向下降, 经过反复学习训练, 确定与最小误差相对应的网络参数 (权值和阈值), 训练即告停止。此时经过训练的神经网络即能对类似样本的输入信息, 自行处理输出误差最小的经过非线形转换的信息。

设一个三层BP神经网络, 输入节点x_i；隐含节点y_i, 输出节点z_k, 输入节点与隐含层节点的网络权值为w_ij, 隐含层节点与输出节点的网络权值为w_ij, 当输出节点的期望输出为t_l时, BP模型的计算公式如下

隐含层节点的输出

(1)

输出节点的计算输出为

(2)

输出节点的误差计算公式为：

(3)

2.2 序列高斯协同仿真

为了得到研究区碳密度分布图, 采用基于地统计学的序列高斯协同仿真模拟来进行仿真。该方法通过将研究区划分成块的方法进行仿真并假设每个单元的估计值是一个随机函数在该位置的随机变量Z(u) 的实现, 其概率分布假定为正态分布, 并由该点周围的样地数据确定。地统计学中, 半方差函数用来描述随机函数空间关系, 相交的半方差函数可以度量两个随机函数相互的空间相关关系。设变量Z为森林碳储量, 则Z(u) 为定义在二维空间u处的随机函数, 其半方差函数γ_ZZ(h)、空间协方差C_ZZ(h) 通过关于距离h的方程式进行计算：

(4)

(5)

式中，α是变程范围内第α个样本；相距h的两个样本称为一个样本对；N为变程范围内样本对的数量；为了区分样本对中的两个样本数据, 分别称为头和尾, m_-h和m_+h分别为若干个样本对的尾和头数据的平均值。

采用序列高斯协同仿真模拟时, 需要由一个统计平均数和方差来确定特定的概率密度函数, 而这个统计平均数和方差可以用同位协同简单克里格估计来获得。

(6)

(7)

式中，u为一个被估计的位置; z^sck(u) 为位置u的森林碳估计值; z(u_α) 为位置u局部范围内第α个样地的森林碳实测值; n(u) 为在搜索范围内所获得的样地数; y(u) 为在像元u处的遥感影像数据; m_z、m_y分别为地面样地数据和遥感影像数据的均值; λ_α^sck(u) 和λ_y^sck(u) 分别为样地数据和影像数据的权重; C_zz(0) 为森林碳地面样地数据的方差; C_zy(0) 为森林碳和遥感影像数据的协方差。当h=(u_α-u) 时, C_zy(h) 为森林碳与遥感影像数据的交叉协方差函数。

由平均数和方差确定的密度函数f[Z(u)]可以用下式表示：

(8)

概率密度函数积分得到其条件累积分布, 并假定这个分布符合正态分布, 该分布中随机抽取一个数作为待估位置u处的模拟实现。本文选取与碳储量相关性最高的TM5参与仿真模拟。在仿真过程中, 样本半方差可以用下式模拟：

(9)

式中, c₀为块金值, c₁为结构参数, a为变程, h为距离。

2.3 模型精度验证

为了检验两种模型的精度, 本研究选取70个样地数据 (建模样本) 用来建模, 剩余29个样地数据 (检验样本) 用来检验模型精度, 两种方法两组数据保持一致。两组数据的统计见表 3。

模型确定后, 以决定系数 (R²)、均方根误差 (RMSE) 以及相对均方根误差 (RRMSE) 对模型拟合精度进行评价。

3 结果与分析 3.1 模型的确定 3.1.1 BP神经网络建模

Funahashi^[24]指出单隐含层BP神经网络模型能够以任意精度逼近任意函数。鉴于本研究输入输出因子相对简单, 选取单隐含层进行模型构建。利用matlab2010b神经网络工具箱及相关程序, 以样地所在位置对应的6个波段信息及相关性较高的MSAVI、RSR、TM5+TM7经归一化处理后作为输入层, 通过多次训练并调整神经网络目标误差以及隐含层神经元个数, 从而选取较为理想的模型来估算研究区森林碳储量, 最后对神经网络输出结果进行反归一化得到森林碳储量的预测值。用前面提到的检验样本数据集检验样地位置的估计值与实测值的吻合程度。计算可知, 检验样本预测值与实测值的相关性为0.67, 决定系数 (R²) 为0.45, 均方根误差 (RMSE) 为20.42 Mg/hm², 相对均方根误差 (RRMSE) 为0.63。其中建模数据与检验数据的预测值基本统计见表 4。

3.1.2 半方差函数拟合

在进行半方差函数拟合时, 通过反复调节块金值、基台值以及变程3个参数来选择较好的组合。当块金值为0.35, 基台值为0.65, 变程为7560时达到最优值。此时半方差函数为：

式中, h为距离, r^sph(h) 为标准半方差。仿真完成后, 提取检验样本所对应的像元值对仿真结果进行精度验证。计算可知, 预测值与实测值的相关性为0.68, 决定系数 (R²) 为0.47, 均方根误差 (RMSE) 为20.04 Mg/hm², 相对均方根误差 (RRMSE)0.63。建模数据与检验数据的预测值基本统计见表 5。

3.2 模拟结果

从统计特征来看, 基于BP神经网络模拟的森林平均碳密度为36.10 Mg/hm², 碳密度最大值为130.92 Mg/hm², 最小值为-107.33 Mg/hm², 碳总量为11042990 Mg, 标准差为20.09 Mg/hm²。基于空间仿真模拟的森林平均碳密度为37.23 Mg/hm², 碳密度最大值为116.78 Mg/hm², 最小值为0.18 Mg/hm², 碳总量为11388657 Mg, 标准差为19.35 Mg/hm²。

在森林碳分布方面, 由于神经网络方法估测结果含有部分负值, 与实际值不相符, 而且由于负值的存在使得森林碳分布表现不明显。因此在制作森林碳分布图时, 本研究将小于0的值赋值为0。图 2为两种方法估测的碳密度分布图。从碳分布来看, 基于空间仿真模拟的森林碳储量较高的区域主要分布在研究区南部, 西北部以及北部。基于神经网络模拟的森林碳储量较高的区域主要分布在研究区东南部和西北部, 整体上呈现出南北高中部低的趋势。与神经网络方法相比, 基于空间仿真方法估算的森林碳储量在分布上相对分散, 未表现出明显的南北多中间少的趋势。两者相差较大的地区在东南部, 基于神经网络估测的碳密度要高于空间仿真方法。

3.3 结果分析

为了对比两种方法的估测精度, 可将两种方法的估测结果与抽样统计估计结果对比, 结果见表 6。表 6显示, 两种方法估测的碳总量要高于系统抽样估计结果, 但神经网络方法估测的碳总量更接近抽样估计结果。与抽样估计结果相比, 基于空间仿真方法估计的碳总量要高9.40%；基于BP神经网络估计的碳总量要高8.82%。两种方法相比, BP神经网络方法估算的平均碳密度要低1.13 Mg/hm², 最大碳密度高14.41 Mg/hm², 最小碳密度低107.51 Mg/hm²。两种方法最小值相差较大是因为神经网络方法在建设用地、水体等区域碳密度估算存在过度拟合的现象, 使得估测值为负值。

判断估测精度的另一种方法, 是将估计结果与样地实测值进行比较。将全部样地实测数据与碳分布图叠加得到森林碳密度对比图 (图 3)。从图 3可以看出, 两种方法均可以在一定程度体现森林碳密度的空间分布格局。但神经网络方法在东南部森林碳密度估测值要高于样地实测值, 与实测样地差距较大。这可能是由神经网络算法本身的固有缺陷造成, 本研究采用BP算法来训练神经网络, 该算法在预测过程中易陷入局部最优, 从而导致了部分区域模拟结果偏高^[15]。

为了更准确地体现两种方法与样地实测值的一致性, 提取全部样地所在像元的预测值, 通过对比分析从数值上来体现两种方法与实测值的差异。从两种方法估计结果中分别提取与地面样地位置对应像元的碳密度估计结果, 以地面样地数据为真值, 分别就两种方法的估计结果进行比较。图 4为两种方法碳密度估测值与实测值的关系图。图中显示, 基于空间仿真估测的结果与样地实测数据的拟合程度较好, 相关性达到0.75, R²= 0.56。而基于神经网络估测的结果与样地实测数据的相关性为0.63, R²= 0.39。图 5为两种方法森林碳密度预测误差曲线图, 显示空间仿真方法的估算结果与实际值更相符, 最小绝对差值为0.07 Mg/hm², 最大绝对差值为52.78 Mg/hm², 均方根误差为14.23 Mg/hm²。基于神经网络方法的最小绝对差值为0.29 Mg/hm², 最大绝对差值为63.01 Mg/hm², 均方根误差为22.59 Mg/hm²。由此可见, 空间仿真方法对于样地位置的估算结果与实测值更接近, 更能表现森林碳密度的分布。

以上分析仅对样地实测数据与预测数据进行了对比分析, 不能从全部像元值上体现两种方法估测结果的差异。表 6显示, BP神经网络含有部分负值, 但碳总量却与空间仿真方法相差不多。这可能与两种方法估测值的频率分布有关, 频率分布曲线和累积频率分布曲线在揭示各值频率分布上具有明显优势。图 6和图 7分别为两种方法预测值的频率分布曲线和累积频率分布曲线。由图 6和图 7可见, 两种方法频率较高的值都集中在20—40 Mg/hm²的值域范围内, 累积频率都为40%。在此值域范围内神经网络平均值为29.43 Mg/hm², 空间仿真方平均值为29.73 Mg/hm²。与空间仿真方法相比, BP神经网络在40—70 Mg/hm²之间有一个明显的突起, 在此值域范围内累积频率百分比为31.5%, 均值为53.13 Mg/hm², 而空间仿真方法频率百分比为34.5%, 均值为51.09 Mg/hm², BP神经网络在此值域的频率分布低于空间仿真方法但均值高于空间仿真方法。两种方法在大于70 Mg/hm²的值相对较少, 累积频率百分比分别为5.88%和7.2%, 神经网络方法也低于空间仿真方法。由此可见, 神经网络方法在估算森林碳密度时, 低值的估计要多于空间仿真方法, 再加上神经网络方法含有部分负值 (所占比例为2.23%), 使得神经网络方法估算得的碳总量要低于空间仿真方法, 也使得神经网络估计的碳总量更接近抽样估计值。但空间仿真方法值域分布更为合理, 极值与样地数据更接近, 而且估计结果没有负值。

为了体现两种方法在碳密度分布上的差异, 将两种方法估测结果进行差值运算, 图 8为空间仿真结果减去BP神经网络模拟结果得到的碳密度差值图, 该图由两种方法对应像元位置的森林碳密度估计值相减得到。由图 8可知, 空间仿真方法在森林碳密度较高的区域要高于神经网络方法, 在研究区东南部却低于神经网络方法。通过前面分析可知, 空间仿真方法与样地实测值更接近, 更能表示研究区真实的碳分布。若以空间仿真方法为标准, 说明BP神经网络方法低估了森林碳密度较高区域的值, 高估了部分森林碳密度较低的值。在像元尺度上对差值图进行统计分析可得, 两种方法估测结果最大差值为89.68 Mg/hm², 最小差值为-74.39 Mg/hm², 均值为1.25 Mg/hm², 标准差18.09。说明两种方法总体上差距不大, 但正负方向上的极值差距较大。由差值频率分布图和累积频率分布图 (图 9) 可知, 两种方法估算结果的差值基本以0为对称轴呈左右分布, 值域在-40—40 Mg/hm²的范围内, 而那些差值较大的点, 只是个别像元的影响。

4 结论与讨论

本研究基于BP神经网络方法和空间仿真模拟法对龙泉市森林碳储量及碳分布进行了估算, 并对两种方法进行了对比分析。通过对比分析可知空间仿真模拟和神经网络方法在碳总量上与抽样数据估计结果相差不多, BP神经网络方法估算的森林碳总量更接近抽样估计结果, 但在值域分布上空间仿真方法更合理。两种方法在森林碳分布上相差较大, 空间仿真方法与实测值更接近, 更能反映真实的碳分布。这主要是因为空间仿真模拟是基于图像的地统计条件模拟技术, 以区域化变量理论为基础, 通过量化随机函数的空间关系, 得到其条件累积分布函数进而取得区域碳密度。它追求的是森林碳储量分布的真实性, 尽可能的接近碳分布的真实情况。而BP神经网络在进行网络训练时追求的是误差最小, 并没有考虑到变量间的空间关系, 在反应局部特征上可能会有一定的差异。

利用神经网络模型模拟森林碳储量分布已有相关学者进行了研究, 但是, 目前的研究内容主要是生物量建模^{[6, 25]}, 关于碳分布的研究区域主要为林场、湿地等小区域^[26-28], 对于县级以上碳分布还少有研究^[15]。小区域森林碳估算, 由于研究区用地类型单一, 植被类型简单, 测结果较为精确, 如翟晓江等^[29]利用神经网络方法对黄龙山林区森林生物量进行可估算预测精度达到87.49%, 相关性达到0.738；李丹丹等^[26]利用BP神经网络建立了旺业甸林场森林生物量模型仿真检验结果的平均相对误差为15.7%, 相对系数达0.8022。大区域森林碳估算时, 土地覆盖与土地利用类型的多样性以及植被的复杂性会对估算精度有一定的影响, 没有小区域森林碳估算精度高, 汪少华等^[15]利用BP神经网络方法估算了临安市森林碳储量, 检验样本与实测值得相关性为0.61, 决定系数为0.37；陈蜀蓉等^[30]利用Erf-BP神经网络对缙云县公益林森林生物量进行了估算, 预测数据的决定系数为0.513。此外, Cutler等^[31]利用多源遥感数据估算了3个不同研究区的森林地上部分生物量。研究表明, 只使用单一研究区数据进行生物量估算, 预测值与实测值的相关性在0.79以上。若将3个研究区数据结合使用, 相关性仅为0.55, 因此不同研究区间的差异也会对生物量估测精度产生一定的影响。

与空间仿真方法相比, 神经网络方法仿真得出的森林碳密度较高的区域较多, 部分地区碳密度预测值小于0。这可能由2方面的原因造成：第一, 神经网络本身的缺陷, 使用误差反向传播算法来训练神经网络, 容易使训练结果陷入局部最优, 从而导致了部分区域模拟结果与实际值不相符。第二, 本文所使用的有效样地数仅为99个, 研究指出如果森林生物量建模野外实测数据的分布和代表性不足, 对应用神经网络估算森林碳密度的估算精度有一定影响^[32], 再加上神经网络外延性较差^[33], 因此对于缺少样地信息的河流等区域的预测能力较差, 使预测值与实际值差别较大。