结合地统计学与Landsat 8影像的乔木林地上碳储量估算

邱新彩，郑冬梅，王海宾，安天宇，许等平，郝月兰，彭道黎

（1.北京林业大学 林学院，北京 100083；2.国家林业和草原局 调查规划设计院，北京 100714；3.国家林业和草原局 林产工业规划设计院，北京 100010）

森林是全球陆地生态系统中最大的有机碳库，在维护生态环境和碳平衡中发挥巨大作用[1]。森林碳储量是研究大气碳收支的重要参数，也是评价国家固碳潜力是否有效的数据支持[2]。准确评估森林碳储量，不仅可以推动碳循环和碳汇研究工作的进展，而且可为森林资源监测和经营管理、可持续发展和生态环境建设提供决策服务[3]。传统的森林碳储量估算方法，如采用森林资源一类、二类调查数据或生态样地的测树资料来估算碳储量，存在费时费力、调查间隔期长、适用范围小等缺点，遥感数据的全覆盖、多时相、多空间分辨率等诸多优势特征，为碳储量估算提供可行路径，与传统估算方法相比更具优势[4]。它可以与地面实测样地数据相结合，建立碳储量遥感估算方程，为快速、准确、实时地估算多尺度森林碳储量提供了有效手段[5-7]。

目前国内外关于森林碳储量研究主要集中在国家或区域尺度上，普遍通过直接或间接测定森林植被或优势树种的生物量，再乘以对应的碳元素含量推算而得[8-9]。Landsat 8 作为新型的Landsat 陆地资源卫星数据，提供了光谱信息丰富、时间分辨率较高、获取速度快的遥感信息源，在林业监测中获得了较好的效果[10-12]，但应用Landsat 8 影像估算森林生物量和碳储量的研究还不多。森林生物量和碳储量的估算方法多以参数和非参数方法为主，其中参数方法以回归分析为主[13-14]，非参数方法包括随机森林、人工神经网络、K 最近邻等的研究和应用较多[15-20]。但是这些方法大多都没有考虑空间异质性对森林生物量和碳储量的影响，缺乏明确的地理空间信息，被认为是非空间方法。地统计学因其兼顾了空间结构特征并具备完好的空间制图功能，被广泛地应用在地质学、生态学等领域。近年来，国内外学者采用不同的地统计学方法对森林变量进行空间插值及验证分析，在森林蓄积量和生物量研究中取得了较好的估算效果[21-23]。贺鹏[24]、王海宾[25]等采用了协同克里格插值法对森林蓄积量、森林生物量进行了空间插值研究，结果均显示了较好的适用性，获得了较高的预测精度。在文章中，作者均未采用遥感数据结合地统计学方法对蓄积量和生物量进行适用性研究，但作者对后续采用遥感数据结合地统计学方法估算变量进行了展望。因此，本文借鉴前人的研究展望，将遥感数据与协同克里格插值法相结合来估算乔木林地上碳储量，试验并讨论这种方法的可行性和适用性。此外，为与非空间估算方法对比，本文将常用的稳健估计方法作为建模方法参与模型构建，并对两种估算方法的精度进行对比分析。

1 研究区概况

浙江省（180°01′～123°10′E，27°06′～31°11′N）地处中国东南沿海，全省陆地面积1.02×105km2，地势从西南向东北呈阶梯状下降趋势，西南以山地为主，平均海拔800 m；中部多为丘陵，海拔在100～500 m；东北部是低平的冲积平原，海拔在10 m 以下。浙江省属亚热带季风性气候，四季分明，雨热同季，年平均气温15～18 ℃，无霜期230～270 d，年平均日照时数1 710～2 100 h，年 均 降 水 量980～2 000 mm。土壤类型多样，主要有红壤、黄壤、粗骨土、水稻土、潮土、盐土等类型。常绿阔叶林为该省的地带性植被。浙江省森林资源丰富，森林面积为6.04×106hm2，占林地面积的91.70%，森林覆盖率达59.43%，森林蓄积达2.81×108m3，占活立木总蓄积的89.58%。

2 数据准备及处理

2.1 专项调查数据

2012年，通过外业调查获取专项调查数据，依据林分年龄（分为幼龄林、中龄林、近成过熟林）和优势树种作为控制，设置和调查了96 个样地，尽量使其具有代表性，样地为0.067 hm2的正方形样地，主要包括针叶林、阔叶林、混交林等不同森林类型，专项调查样地分布见图1。样地属性因子包括树种、胸径（起测胸径为5 cm）、平均高、郁闭度、坡向、坡度、海拔、株数、年龄、乔木林地上生物量和乔木林地下生物量等因子。本研究用到的因子为乔木林地上生物量因子，样地调查数据的生物量统计结果见表1。

图1 研究区内调查样地分布Fig.1 Distribution of survey samples in the study area

表1 调查样地生物量统计结果Table 1 Statistic results of biomass for the survey samples in the study area

依据获取的样地蓄积量，乘以生物量换算因子（Biomass expansion factor，BEF）得到乔木林地上生物量，为便于计算，将乔木林地上生物量转换为每公顷生物量（t·hm-2），每公顷生物量乘以相对应的含碳率即可得到乔木林地上碳储量。蓄积量直接采用已计算好的样地蓄积量数据。BEF采用沈楚楚[24]提供的浙江省4 大树种（组）的BEF，分别为杉木林（0.745 3 t·m-3）、马尾松林（0.883 9 t·m-3）、软阔类（1.657 2 t·m-3）、硬阔类（1.070 5 t·m-3）。将其他树种分别归并到软阔类和硬阔类树种（组）并采用相应的BEF 进行生物量计算。本文依据《土地利用、土地利用变化与林业碳汇计量监测技术指南》确定树种含碳率，指南中的树种含碳率可覆盖研究区内的各个优势树种（组）。

在获得乔木林地上碳储量数据后，采用二倍标准差法剔除异常值，结果剩余95 个专项调查数据，将其作为因变量用于本文Landsat 8 影像构建乔木林地上碳储量模型。

2.2 辅助数据

本研究采用的辅助数据为覆盖研究区的数字高程模型（Digital elevation model，DEM）数据，比例尺为1:50 000，用于提供地形等因子。

2.3 遥感数据及预处理

本研究选用一景Landsat 8 OLI（Operational Land Image，陆地成像仪传感器）影像作为遥感数据源，覆盖整个研究区，空间分辨率为30 米，影像轨道号为119/040，采集时间为2014年10月26日。所采用波段包括波段序号为2～7 的6 个波段，应用ENVI5.3 软件对Landsat 8 影像进行预处理，包括辐射定标、大气校正。

3 研究方法

3.1 建模变量设置

本研究选用植被指数（9 个）、纹理波段（48个）、主成分变换因子（3 个）、缨帽变换因子（3个）和地形因子（3 个）共66 个特征作为自变量参与模型构建。按照地面实测样地位置，提取各个样地对应的自变量因子。表2为本研究采用的9个植被指数。

纹理是一个非常复杂的空间属性，易受到传感器观测角度以及观测对象所处环境的不同而产生较大的变化，此外纹理特征的度量易受到窗口大小、方向以及步长的影响。为探究Landsat 8 多光谱波段纹理特征的应用潜力，本研究采用灰度共生矩阵方法分别提取Landsat8 多光谱波段的8种纹理特征，8 种纹理特征的计算见表3。

表2 研究中用到的植被指数†Table 2 Vegetation indices used in the study

表3 研究中用到的纹理特征†Table 3 Texture features used in the study

影像的多光谱波段间常存在较高的相关性，图像变换可以将图像信息压缩到少数几个有效波段，去除冗余信息。目前较为常用的有主成分分析法和缨帽变换法，在以往的研究中，主成分分析法的前3 个波段和缨帽变换的前3 个波段应用较为广泛[10,27]。本研究采用主成分分析获取的前3 个波段和缨帽变换的前3 个波段（亮度Brightness、绿度Greenness 和湿度Wetness）共计6 个波段变换因子作为自变量，通过分析确定相关自变量参与碳储量模型的构建。此外，应用DEM 数据提取对应的海拔、坡度、坡向3 个因子，参与模型构建。

3.2 建模变量优选

采用皮尔森相关性分析法对自变量进行筛选，提取出极显著相关的因子，进一步采用方差膨胀因子（Variance inflation factor，VIF）法筛选自变量，将筛选的自变量用于模型构建。VIF 是指解释变量之间存在多重共线性时的方差与不存在多重共线性时的方差之比，为容忍度的倒数，VIF 越大，显示共线性越严重[28]。经验判断方法表明：当0＜VIF＜10，不存在多重共线性；当10≤VIF＜100，存在较强的多重共线性；当VIF≥100，存在严重多重共线性。依据相关的文献资料[29-30]，本研究设置0＜VIF＜10，容忍度＞0.1 来进行控制，对自变量做进一步提取。

3.3 稳健估计

稳健估计是实际情况偏离假设模型的不敏感性估计，主要目的是在建模中寻求样本数据分布的主体，不迁就异常样地，采用一定的抗差函数削弱其使用比例，以此来孤立异常样地和抵制其不利影响，从而使估测模型具有一定的“稳定性”和“抗干扰性”[31-32]。稳健估计实质上是求非线性函数的极值解，采用迭代函数法求解，具体步骤及原理详见文献[32]。

3.4 协同克里格插值

协同克里格法（Co-kriging，COK）在变量估测上具有较好的预测精度[24-25,33]，它是利用两个或两个以上的变量，将其中一个作为主变量，剩余变量作为辅助变量，将主变量的自相关性和主辅变量的交互相关性结合起来进行无偏最优估计[30]，详细的理论介绍可参考文献[24]和[25]。协同克里格插值的公式如下：

式中，Z*(x0)为待估点处的蓄积量估测值；λ1i和λ2i分别是主变量Z1和辅变量Z2实测值的权重，且n和m分别是参与估测x0点的蓄积量和其它辅助变量Z2的实测值数目。

3.5 模型验证

考虑到覆盖单景Landsat8 影像的专项调查样地数量较少的情况，本文采用留一交叉验证法进行评价。通过对比模型预测的乔木林地上碳储量与实测碳储量值，比较不同建模方法的预测精度，以评价不同模型的估算性能。本研究采用的评价指标有：决定系数（R2）、均方根误差（Root mean square error，RMSE）、平均绝对偏差（Mean absolute error，MAE）和总预报偏差的相对误差（RE），R2越大，模型的拟合精度越好，RMSE、MAE 和RE 值越小，表示模型预测精度越高。4 个评价指标计算公式如下：

式中：yi为碳储量实测值，为碳储量预测值，为碳储量实测值平均值，n为样本个数。

4 结果与分析

4.1 描述性统计分析

使用SPSS 20.0 统计软件，对单位面积乔木林地上的碳储量进行了描述性统计，并对其进行K-S检验，结果显示乔木林地上碳储量数据服从正态分布。表4是单位面积乔木林地上碳储量的描述性统计和正态分布性检验结果。

表4 乔木林地上碳储量描述性统计结果Table 4 Descriptive statistics of carbon storage in the arbor forest

4.2 自变量筛选结果

应用SPSS 20.0 软件，对提取的专项调查数据的乔木林地上碳储量和66 个自变量因子进行皮尔森相关性分析，筛选出与碳储量极显著相关（P＜0.01）的自变量因子共计22 个，见表5，包括NDVI、RVI、RVI54、RVI64、SAVI、NLI、ARVI、PCA2、Greenness、海拔、B1_Mean、B2_Mean、B3_Mean、B3_Entropy、B3_Second、B3_Correlation、B6_Mean、B6_Variance、B6_Contrast、B6_Diss、B6_Entropy、B6_Second。

对提取出的极显著相关的22 个自变量因子进一步优选，生成基于VIF 方法筛选的自变量集，筛选后的7 个自变量因子见表6，包括RVI、NLI、海拔、B1_Mean、B3_Correlation、B6_Mean、B6_Variance。

4.3 稳健估计

将筛选的自变量集构建乔木林地上碳储量模型，应用MATLAB2014a 计算各个模型的参数，最后得到稳健估计法构建的乔木林地上碳储量估算模型，其数学表达式为C=0.000 1×RVI+2.100 6×NLI+0.000 1×海拔-0.000 2×B1_Mean+0.000 5×B3_Mean-0.000 4×B6_Mean-0.000 5×B6_Variance，C 为乔木林地上碳储量，其他因子为对应的自变量因子。

表5 自变量因子与乔木林地上碳储量的Pearson 相关系数†Table 5 Pearson correlation coefficients between independent variables and arbor forest above ground carbon storage

表6 基于VIF 法筛选的自变量因子†Table 6 Independent factors filter based on the VIF method

4.4 协同克里格插值

基于VIF 法筛选的自变量，提取出与乔木林地上碳储量相关性最高的1 个自变量（RVI）和相关性最低的2 个自变量（B6_Mean、B3_Correlation），采用协同克里格插值对乔木林地上碳储量进行估算。为取得更好地估算效果，本研究对ArcGIS 10.1 地统计分析模块中的11 种变异模型进行了计算和分析，根据变异函数的理论和评价方法[25]，得到最优变异函数的拟合结果，变异函数拟合结果见表7。由表7可知，在协同克里格插值中，采用提取的3 个变量参与插值后，块金值为0，块金值/基台值小于25%，表明系统具有很强的相关性，依据变异函数的评价标准可知，在使用RVI、B6_Mean、B3_Correlation 时，K-贝塞耳模型可获得最优的拟合结果。为获得较高的预测精度，本研究选择VIF 中的RVI、B6_Mean、B3_Correlation 3 个变量拟合的最优变异函数（K-贝塞耳模型）作为最优变异函数，最后采用协同克里格法对乔木林地上碳储量进行插值，并对插值效果进行评价。

表7 基于VIF 提取的3 个自变量拟合的最优方差函数模型及参数Table 7 Optimal variance function model and parameters fitted by three different variables extracted by VIF

4.5 模型评价及对比

应用2 种方法对乔木林地上的碳储量进行估算，并采用4 种评价指标对预测结果进行比较，结果见表8。由表8可知，协同克里格插值法的决定系数（R2=0.45）高于稳健估计（R2=0.42），表示协同克里格插值的拟合精度较好。协同克里格插值法的均方根误差、平均绝对偏差和总预报偏差的相对误差（RMSE=9.88 t·hm-2，MAE=7.75 t·hm-2，RE=0.24%）均明显低于稳健估计法（RMSE=10.15 t·hm-2，MAE=8.03 t·hm-2，RE=0.27%），表示协同克里格模型预测精度更高。由以上分析可知，协同克里格插值法对乔木林地上碳储量的估算性能更强，具有较高的空间模拟预测精度。

表8 2 种建模方法的预测精度比较Table 8 Comparison of prediction accuracy of two models

5 结论与讨论

本研究基于Landsat 8 影像和DEM 数据提取植被指数、纹理特征、主成分波段、缨帽变换以及地形因子，结合研究区专项调查样地数据和辅助数据，采用皮尔森相关系数和VIF法提取自变量，基于筛选的自变量，采用稳健估计和协同克里格插值法构建乔木林地上碳储量模型，并对模型精度进行对比分析。结果表明协同克里格插值法构建模型的精度要优于稳健估计。

在自变量筛选上，采用了两步法对自变量进行了筛选，以期为建模提供更好的自变量因子。皮尔森相关系数法是用来反映因变量和自变量间的相关程度，通过剔除不相关的自变量因子，提取出与因变量相关性较高的自变量因子用于建模。但在实际的应用中，采用遥感影像衍生出的自变量因子间往往存在多重共线性问题，如何减小这些问题对模型稳定性的影响，也是变量筛选的关键因素之一。方差膨胀因子法在这方面显示了其优势，通过考察给定的解释变量被方程中其他所有解释变量所解释的程度，以此来判断其是否存在多重共线性，通过一定的标准剔除影响较大的自变量因子，减小多重共线性影响。因此，本文将两种方法结合使用，既兼顾了因变量与自变量间的相关性，又兼顾了自变量间的共线性问题，在一定程度上可提高所构建模型的稳定性和适用性。

研究结果显示协同克里格插值法的估算效果要优于稳健估计法，这是因为协同克里格插值法将主变量的自相关性和主辅变量的交互相关性结合起来，对乔木林地上碳储量进行无偏最优估计，可获得较好的估算精度。此外，本文将协同克里格插值和遥感影像相结合的方式，考虑了遥感自变量的空间连续特征，与前人只考虑不具有空间连续性特征的自变量进行插值相比，更具有优势，可以说是对文献[22]和[23]的探索式研究。

采用协同克里格插值法构建乔木林地上碳储量模型时，本研究选择了相关性最高的RVI 自变量和相关性最低的B6_Mean、B3_Correlation 两个自变量，这是因为在采用与乔木林地上碳储量相关性最好的3 个自变量插值时，由于变量间存在的共线性问题，导致插值精度较低，在试验不同3 个自变量插值时，也显示了同样的效果。当选择RVI、B6_Mean、B3_Correlation 进行插值时，可获得最高的精度。因此，在采用协同克里格插值法时，选取相关性较高且共线性较弱的自变量因子对因变量进行插值具有较好的估测效果和适用性。

在自变量因子上，地形因子与乔木林地上碳储量显示了较好的相关性，在经过两步筛选后，海拔因子仍被优选出来作为建模变量，说明海拔因子能够反映乔木林地上碳储量的空间垂直分布特征。但本研究在采用协同克里格插值法构建模型时，并未将海拔因子参与建模。本研究采用的DEM 数据比例尺为1:50 000，属于较小比例尺数据，在一定程度上均化了地理空间特征信息，未能充分挖掘出地理空间特征信息对建模精度的贡献。因此，如何进一步挖掘地形因子信息、采用适宜比例尺DEM 数据作为建模自变量，充分利用其特有的空间特征信息，提高所构建模型的稳定性和精度等问题有待于进一步研究。