Estimation and Uncertainty Analysis of Aboveground Carbon Storage of Pinus densata based on Random Forests and Monte Carlo

研究区位于云南省迪庆州香格里拉市，99°20′～100°19′ E、26°52′～28°52′ N，总面积为11 613 km²。地形起伏大，海拔高差达到4 042 m，森林覆盖率达74.99％，森林资源丰富。主要优势树种为云杉( Picea asperata Mast.) 、云南松( Pinus yunnanensis Franch.) 、高山松( Pinus densata Mast.) 、高山栎( Quercus semicarpifolia Smith.)等^[18]。高山松是香格里拉市的优势树种之一，占全市面积的16.18％。纯林主要分布在2 600～3 500 m海拔的向阳山坡和河流两岸。

使用USGS下载的Landsat8-OLI影像数据，选取3景2015年云量较少的影像，具体信息如表1所示。影像使用ENVI5.3软件进行预处理。通过辐射定标、大气校正，消除地形和气溶胶对地物反射的影响。最后使用坡度匹配法进行地形校正，两次校正后，北坡和南坡的地形阴影平均值接近，阴坡部分的反射率得到阳坡部分反射率的值补偿。

研究区共60个样地，样地调查时间为2015年11月至2016年3月。样地大小为30 m × 30 m，每块样地间隔3 km以上。详细记录了森林树种、胸径（DBH）、树高、样地中心坐标、海拔、坡度、坡位等属性^[19]。使用二元材积表计算每个样地中林木的蓄积量。

(1)高山松单株生物量模型^[20]如下：

其中，W为单株生物量，D为胸径，H为树高。

(2)高山松的碳储量计算公式为：

C_t为高山松碳储量，f_c为含碳系数，含碳系数0.513 1^[21]。外业调查数据见表2。

以碳储量为因变量，遥感影像因子为自变量。根据各变量与碳储量的相关性选择最优因子，构建基于遥感特征因子的碳储量估测模型。为研究纹理特征对碳储量估测的影响，设置像元大小为5 × 5、9 × 9的窗口，以提高模型的性能^[18,22]。

为了提高模型的精度，本研究提取的遥感因子主要包括7个单波段因子、5个植被指数、8个波段组合因子、3个信息增强因子、9个纹理因子、3个地形因子、1个植被生长因子。具体变量信息如表3所示。

随机森林（RF）是一种基于决策树的机器学习方法，其特点是采用带替换的随机抽样方法从样本数据中提取一组n个数据作为训练集，根据这些数据构建分类回归树。RF是目前最有效的非参数回归模型之一。与参数回归方法相比，该算法不需要检验变量的正态性和独立性等假设。在一定程度上可以避免过拟合现象的出现，对异常值的处理效果较好，可处理高维度数据^[23]。能够评估各个特征在模型中的重要性。节点的分裂属性通常用基尼系数（Gini）或信息熵来评估，基尼系数可用于计算RF回归中每个模型特征因子的重要性。将变量重要性评分用V来表示，假设有m个特征x₁, x₂, x₃, ···, x_c,现在要计算出每个特征x_j的Gini指数评分V_j^Gini,亦即第j个特征在RF所有决策树中节点分裂不纯度的平均改变量。计算公式为：

其中k为类别数，P_mk是k在m个节点中的占比。

特征x_j在节点m的重要性，即节点m分支前后的Gini指数变化量如下：

Gini_L和Gini_R分别表示第m个节点前后的Gini指数。每个特征的重要性可以通过对所有特征的基尼指数进行归一化来获得。

蒙特卡洛模拟（MC）的基本思想是反复模拟一个随机事件的发生，并根据随机事件发生的频率估计其概率特征。由于区域尺度碳储量估测模型中的不确定性来源复杂且难以测量，因此MC方法在解决这个问题方面具有显著优势。反复模拟碳储量模型的建立和估测，得到碳储量估测值和误差的概率分布^[24-25]。碳储量估测值的误差用均方根误差（RMSE）和相对均方根误差（rRMSE）表示。具体方法如下：

(1)使用遥感模型对所有样地的碳储量进行无偏估计 $ {\widehat{g}}_{i} $ 。

(2)假设残差服从 $ {\epsilon }～\mathrm{N}\left(0,{\sigma }^{2}\right) $ 的正态分布，采用一个符合 $ \epsilon \beta～\mathrm{N}\left(0,{\sigma ̂}^{2}\right) $ 的随机生成的数组来模拟新的残值 $ \epsilon \beta $ ，其中 $ {\sigma ̂}^{2} $ 是碳储量 $ {\widehat{g}}_{i} $ 的方差。

(3)将碳储量的估测值 $ {\widehat{g}}_{i} $ 与新的残差 $ \epsilon \beta $ 相结合，得到了新样地碳储量 $ {g}'_{i} $ 。计算方差 $ Var\left({g}'_{i}\right) $ 的公式为：

式（5）中，n为样地数。

(4)重复第2步和第3步，直到 $ Var\left({g}'_{i}\right) $ 趋于稳定，并且在m次模拟后，计算每个样地的碳储量均值μ、均方根误差RMSE和相对均方根误差rRMSE；计算公式为：

由于波段反射率、植被指数、样地数据与碳储量之间具有较高的相关性，因此根据各样地碳储量的差异，选取Stock volume、B53、B73、VFC、ND53进行相关性分析。所选特征因子之间的相关性如图1所示。右上角的图为成对特征散点图，沿对角线分布的为密度直方图，左下角为成对特征的核密度图。密度直方图显示了碳储量和蓄积量之间的正相关关系，核密度图和散点图显示出碳储量与蓄积量的相关性较强，因此把蓄积量作为建模因子。根据贝叶斯相关矩阵图可以看出5个特征因子与碳储量的相关性较强，表明波段反射率、植被指数、样地数据与碳储量的相关性较强，所以选取Stock volume、B53、B73、VFC、ND53 等5个因子构建碳储量模型。

Figure 1.  Schematic diagram of Bayesian correlation matrix of each characteristic factor

计算所选多光谱特征的各种属性，包括单波段光谱数据的均值、变异系数、标准差、方差、植被指数、波段比值、DEM数据、样地数据和基于GLCM的纹理特征。通过RFE去除具有高相关性的变量后^[26]，共保留了17个自变量。模型精度不一定会随着自变量数量的增加而增加，建模中包含过多的变量因素反而会增加模型的复杂度；因此，应剔除对模型精度没有显著贡献或降低模型精度的某些变量。只选择有代表性的变量进行建模。当使用10个因子建模（VFC、蓄积量、B53、B73、ND53、坡向、R9B4CR、R9B5VA、R5B3SM、R5B3SK）时，模型达到最佳精度（RMSE=5.539 t·hm⁻²,MAE=4.319 t·hm⁻²）。所有潜在变量的个数与模型精度的关系如图2所示。

Figure 2.  Relationship between the number of variables and the accuracy of the model

本研究中的模型是使用RF方法构建的。为量化不同类型变量对模型的拟合度和不确定性的影响，选取10个自变量，根据属性类型组合成4种不同的方案。模型1选择VFC、蓄积量、B53、B73、ND53进行建模；模型2选择坡向、VFC、蓄积量、B53、B73、ND53进行建模；模型3选择VFC、蓄积量、B53、B73、ND53、R9B4CR、R9B5VA、R5B3SM、R5B3SK进行建模；模型4选择坡向、VFC、蓄积量、B53、B73、ND53、R9B4CR、R9B5VA、R5B3SM、R5B3SK进行建模。最后，用MC来计算每个模型的不确定性。结果如表4所示。

结果显示，基于光谱和样地数据的RF-MC模型精度最低，R²为0.549，RMSE为11.289 t·hm⁻²，MAE为8.552 t·hm⁻²。分别引入DEM和纹理特征后，模型精度有所提高。其中引入纹理特征的模型精度提高较大，说明纹理特征可以有效提升模型精度。同时引入DEM和纹理特征的RF-MC模型精度最优，R²为0.892，RMSE为5.539 t·hm⁻²，MAE为4.319 t·hm⁻²。该模型比其他相关研究的拟合程度高，如温小荣等人对建德市杉木生物量的估测R²仅为0.8^[27]。

由表4可知，模型的不确定性随R²的增大而降低，模型的残差变异会影响预测结果，从而影响预测值的方差。因此残差值的增大会导致模型R²降低，模型的预测结果受到限制，不确定性随之增加。4种模型的不确定性随着DEM和纹理特征的引入逐渐降低，rRMSE分别为38.13%，31.97%，23.76%，18.70%。随着DEM和纹理信息的同时引入，原本仅基于光谱和样地数据模型的不确定性减少了19.43％，说明DEM和纹理特征对降低高山松地上碳储量反演模型的不确定性有一定作用。

4种模型中不同变量的重要性如图3所示。

Figure 3.  Relative importance of different variables

Gini指数用于计算RF回归中每个模型特征的重要性。每个因子对模型的重要性如图3所示，第1、2种模型中蓄积量占比最多；添加纹理特征后，模型3和模型4中纹理因子的贡献率超过蓄积量，但蓄积量仍然在各模型中占比很大，说明把蓄积量引入作为建模因子有很大的发展潜力。纹理特征作为建模贡献率最高的因子可有效提高模型的精度。

为了量化不同类型变量对模型精度的影响，本研究将变量组合成4种模型，分别进行建模和预测。预测结果的散点图如图4所示。

Figure 4.  Prediction results of different models

根据图4，模型精度与拟合度呈正相关。基于光谱和样地数据的RF-MC模型，模型拟合度R²为0.549，逐步引入DEM数据和GLCM纹理特征后，R²最终增加到0.892。随着不同类型变量的引入，R²逐渐增大，模型的不确定性得到有效降低。表明RF-MC方法在估测地上碳储量方面效果显著，建模因子的选择对模型的精度及不确定性有显著影响。

遥感影像提取的数据只能提供森林水平方向的光谱信息。胸径和树高两个因子反映树木的垂直结构。这些因子与碳储量估测密切相关，减少数据饱和问题^[28]。Wang^[29]使用iPad Pro LiDAR传感器成功估算了树木的DBH。Duan^[30]使用机载LiDAR和DEM数据基于树冠分割和地形偏差引起的修正方法提取树高。本研究中，引入蓄积量作为建模因子来提供森林的垂直结构信息，不受数据饱和度的影响，贡献率高（图3）。因此，采用水平和垂直值相结合的方式获取森林结构数据，以提高植被复杂地区的预测能力。如Sinha^[31]利用光学和合成孔径雷达（SAR）数据成功估测印度热带落叶林的碳储量。

相关研究表明温度、二氧化碳浓度和降水对碳储量估测具有复杂的影响。单独升高温度会增强植物蒸腾作用，导致光合作用受到抑制和碳储量减少。上述3个因素的综合作用显著改善植物的生长^[32-33]。气候变化通过改变树木物候和生长季节长度对树木的生长产生影响^[34]。气候因子被认为是一个气候周期的平均值，在时间上与碳储量估测数据不同步。而提高建模精度的关键是解决气候因子和碳储量估测数据的同步问题。地形和气候因素对植物地上生长的影响较大，不同地理特征和气候条件的地区反演碳储量需根据实际情况适当地选择和设置可变因素和模型参数。

研究表明，高分辨率数据能够提供丰富的纹理信息，添加纹理信息可有效提高碳储量预测模型的精度^[35]。由于高分辨率数据并未完全覆盖研究区，则选择空间分辨率为30 m的Landsat8-OLI遥感影像数据，该数据估测碳储量模型拟合度R²=0.892。综上所述，具有丰富纹理信息的高分辨率数据在碳储量估测中的应用有待进一步探索。

研究使用RF-MC方法的优势在于，可以对碳储量的不确定性进行稳定估测^[36]。传统的模型分析方法被引入MC，经多次模拟，模型参数及其协方差矩阵的可变性逐渐减小^[37]，降低估测结果的残差水平。随着模拟次数增加，碳储量估测值的不确定性趋于稳定^[38]。傅煜的研究结果表明，蒙特卡洛的不确定性值达到平稳状态的运算时间随建模样本量的增大而减少^[11]。McRoberts^[39]在美国明尼苏达州东北部设立研究区并采用蒙特卡洛方法分析了模型预测单木材积的不确定性，并指出利用模型计算大面积森林蓄积量引起的不确定性较小。改进MC方法是对输入参数进行有效抽样，减少输入参数组合的数量，以提高计算效率，如增大计算区域，样本量增加能提高估测精度降低不确定性，MC结果达到稳定所需要的运算时间逐渐缩短，可以在一定程度上提高工作效率。MC可以有效降低模型的不确定性，为研究碳储量和类似应用提供了显著优势^[40]。碳储量预测值在上下端波动较大，中间较小，目前的技术和设备无法完全消除碳储量预测中的数据饱和度问题^[41-42]。为降低饱和度对碳储量预测的影响，可以增加信息量更大的纹理特征、垂直结构信息、气候因子等不受数据饱和度影响的指标，有助于提高碳储量估测的精度，获得更好的预测结果。

本研究以高山松为研究对象，采用Landsat8-OLI遥感影像数据、DEM和样地数据，建立不同特征组合的RF回归模型，同时使用MC方法计算每个模型的不确定性。得出以下结论：①引入地形因子碳储量估测模型的精度得到提升，R²提高0.121、RMSE降低1.824 t·hm⁻²、MAE降低0.854 t·hm⁻²，不确定性降低6.16％。②纹理特征在原始图像亮度的基础上增强了空间信息的识别程度，可以显著提高森林地上碳储量的估测精度。③合理选择建模变量可以提高模型的精度，降低模型的不确定性。④RF-MC方法可以为区域尺度上的碳储量估测及其不确定性分析提供稳定可靠的结果。