基于地统计学和多源遥感数据的森林碳密度估算*

严恩萍 赵运林 林 辉 莫登奎 王广兴,3

(1.中南林业科技大学 林业遥感大数据与生态安全湖南省重点实验室 长沙 410004； 2.中南林业科技大学林学院 长沙 410004；3.南伊利诺伊大学地理系 卡本代尔 629012)

基于地统计学和多源遥感数据的森林碳密度估算*

严恩萍1,2赵运林1,2林 辉1,2莫登奎1,2王广兴1,2,3

(1.中南林业科技大学 林业遥感大数据与生态安全湖南省重点实验室 长沙 410004； 2.中南林业科技大学林学院 长沙 410004；3.南伊利诺伊大学地理系 卡本代尔 629012)

【目的】 基于遥感影像空间分辨率和地面样地大小不一致的现象，采用地统计学和多源遥感数据进行森林碳密度估算，为MODIS数据在区域森林碳密度估算领域的应用提供参考。【方法】 以湖南省攸县为试验区，首先利用基于块的序列高斯协同模拟算法，将25.8 m×25.8 m的样地数据分别上推到250 m × 250 m、500 m × 500 m和1 000 m×1 000 m； 然后将上推后的样地数据分别与MOD13Q1、MOD09A1、MOD15A2数据结合，利用序列高斯协同模拟算法开展区域森林碳密度估算研究； 最后将最优结果用于湖南省森林碳密度估算。【结果】 Landsat5和MODIS数据与森林碳密度的敏感因子具有高度相似性，排在前3位的分别为1/TM3、1/TM2、1/TM1和1/Band1、1/Band4、1/Band3； 与植被指数产品MOD13Q1和MOD15A2相比，多光谱数据Landsat5和MOD09A1在攸县森林碳密度估算方面显示出巨大潜力，估算精度分别为82.02%和75.64%； 基于MOD09A1的序列高斯协同模拟算法具有很好的适用性，可用于湖南省森林碳密度的空间模拟，估算精度为74.07%。【结论】 采用基于块的序列高斯协同模拟算法，可以实现由地面样地到不同空间分辨率MODIS像元之间的转换； 由于空间分辨率的限制，MOD09A1数据在刻画空间细节方面不如Landsat5精细。该研究方法适用于地面调查样地大小和遥感影像空间分辨率不一致的区域森林碳密度估算。

林业遥感； 森林资源清查； 多源遥感； 基于块的序列高斯协同模拟； 森林碳密度

作为陆地生态系统的主体，森林贮存了陆地生态系统有机碳的76%～98%，在全球碳循环中发挥着重要作用(刘畅等, 2014; Yanetal., 2015)。森林碳蕴含着丰富的地表结构和功能信息，近年来，随着全球气候变化和“温室效应”的加剧，森林碳密度的空间分布和动态变化研究日益受到重视(严恩萍等, 2015; Zhaoetal., 2016)。遥感技术的快速发展和国家森林资源清查体系的持续推进，基于遥感影像和地面样地数据的估算方法日趋成熟，如回归分析、神经网络、KNN和空间模拟等(Tianetal., 2012; Wangetal., 2004a; 2009; Luetal., 2012; 戚玉娇等, 2015)，但受调查成本的限制，目前森林资源调查样地大小主要分布于10 m × 10 m～50 m × 50 m之间，然而区域森林碳密度制图采用的遥感影像空间分辨率通常在100 m × 100 m～1 000 m × 1 000 m。如湖南省森林资源清查的样地大小为25.8 m × 25.8 m，但是Ladsat5和MODIS数据的像元大小分别为30 m × 30 m和250 m × 250 m、500 m × 500 m、1 000 m × 1 000 m，存在明显的遥感影像空间分辨率和地面样地大小不一致的现象。因此，需要建立一种尺度转换关系，将小尺度的地面样地同大尺度的遥感像元结合起来。

目前，常用的尺度上推方法有最邻近像元法、窗口平均法和块克里格法。最邻近像元法假设较大的块由N个较小的像元组成，将离块中心最近的像元值分配给块，该方法虽然简单，但是效果不好，无法真实反映较小单元数据间的细微变化； 窗口平均法是一种简单且使用广泛的尺度上推方法，即假设窗口内所有像元具有相同的权重，将窗口内像元的平均值赋值给较大的窗口，但忽略了样本数据间的空间自相关性(Wangetal., 2004b)； 块克里格法是通过计算块内较小像元的克里金估计(冉有华等, 2009)，然后计算所有子像元平均值的一种地统计学方法，林业上常用的是序列高斯协同模拟算法。如张茂震等(2009； 2014)将森林资源清查样地数据与Landsat5影像结合，采用序列高斯协同模拟算法绘制了临安市和仙居县森林碳的空间分布，发现模拟结果与地面样地数据具有较好的空间一致性。沈希等(2011)以临安市为研究区，比较了序列高斯协同模拟和一元二次非线性回归算法刻画森林碳分布的真实程度，发现前者模拟结果更接近地面样地估算结果。但这些研究多单纯采用Landsat5数据，缺少利用序列高斯协同模拟算法结合多源遥感数据开展区域森林碳密度估算的研究，特别是结合MODIS这种多空间分辨率的大尺度数据。

鉴于此，本研究首先以森林资源清查数据和Landsat5影像为信息源，利用基于块的序列高斯协同模拟算法，将25.8 m × 25.8 m的样地数据分别上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m， 然后将MODIS数据和上推后的样地数据结合开展森林碳密度估算，最后将最优结果用于湖南省森林碳密度的估算研究，以期为MODIS数据在区域森林碳密度估算领域的应用提供参考。

1 材料与方法

1.1 研究区概况

研究区位于湖南省东部的攸县(图1)，113°09′09″—113°51′30″E，26°46′34″—27°26′30″N。境内四季分明，雨水充足，土壤肥沃，属中亚热带季风湿润气候常绿阔叶林带，年均气温17.8 ℃，无霜期292天，年降水量1 410 mm左右。全县现有林地166 866.7 hm2，其中森林蓄积量和森林覆盖率分别达到311.87 m3和57.24%。

图1 研究区位置Fig.1 Geographic position of study area

1.2 数据来源与预处理

1.2.1 固定样地数据 现有森林碳密度估算一般通过直接或间接测定获取森林生物量，再乘以生物量中的含碳率推算而得。本文以攸县2009年森林资源清查数据为基础，利用生物量回归方程分树种(组)计算生物量，再乘以相应树种的含碳率系数，获取研究区固定样地碳密度数据。固定样地是以4 km × 8 km抽样间隔设置的正方形样地，样地面积0.067 hm2。样地调查因子包括样地类别、地类、龄组、优势树种、树高、胸径等。

采用表1列出的回归模型(李海奎, 2010)，分树种(组)计算样地生物量，混交林生物量分别按比例[6杉(Cunninghamialanceolata)4马(Pinusmassoniana)、5软阔5硬阔、3.6杉2.4马2软阔2硬阔]计算，经济林、灌木林分别按平均生物量23.7 t·hm-2和19.76 t·hm-2计算(吴丹等， 2011)，上述3类树种含碳率均取0.5。对于无明确生物量回归模型的树种，采用近似树种参数替代。

将研究区固定样地和遥感影像叠加，剔除个别云覆盖样地，最后保留78块样地作为研究样本。通过随机抽样将样本数据分为建模样本和验证样本，抽样比例分别为2/3和1/3； 结合湖南省森林资源分布情况，采用分层抽样选取1/2的固定样地作为测算样本，剔除少量云覆盖样地，最后保留2 892块样地作为研究样本，其中建模样本占3/4，验证样本占1/4。1.2.2 遥感数据 采用Landsat5和MODIS 2种遥感影像。Landsat5数据空间分辨率为30 m × 30 m，含云量均低于1%，影像质量较好，用于攸县森林碳密度遥感反演的接收时间为2009年8月21日(轨道号123/41)，用于湖南省森林碳密度反演的Landsat5数据接收时间依次为2009年8月21日(轨道号123/40、123/41、123/42)、2009年8月28日(轨道号124/40、124/41、124/42)和2009年7月18日(轨道号125/40、125/41、125/42)。数据预处理包括辐射定标、大气校正和几何校正3种操作； 数据变换包括倒数运算、比值运算(两波段、三波段和四波段组合间的比值运算)、植被指数运算(NDVI、EVI、ARVI和SAVI 4种指数)、主成分变换和纹理变换(均值、角二阶矩、对比度、相关、相异、熵、逆差距和方差8种纹理因子)，参与运算的波段包括TM1、TM2、TM3、TM4、TM5和TM7共6个原始波段，共计生成88种基于Landsat5的遥感变量。

①WS表示树干生物量，WP表示树皮生物量，WB表示树枝生物量，WL表示树叶生物量，WT表示地上部分总生物量，WR表示地下部分生物量，D表示样地平均胸径，H表示样地平均树高。WS: stem biomass;WP: bark biomass;WB: branch biomass;WL: leaf biomass;WT: lotal aboveground biomass;D: tree diameter at breast height;H: tree height.

MODIS数据包括MOD13Q1、MOD09A1和MOD15A2共3种产品，空间分辨率分别为250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m，接收时间均为2009年8月21日(轨道号h28v06，用于湖南省森林碳密度反演的MODIS数据轨道号包括h27v05、h27v06和h28v06)，3种产品的预处理包括投影转换、影像镶嵌、影像裁剪以及图像修复4个步骤，均通过遥感软件ENVI 5.0实现。其中图像修复采用IDL语言编程实现(崔丽华等， 2009； Xiangetal., 2013)，包括条带去除和去云处理2部分。同Landsat5，MOD09A1数据最后选取Band3、Band4、Band1、Band2、Band6、Band7共6个波段分别进行倒数运算、比值运算、植被指数运算、主成分变换和纹理变换，共计生成88种基于MOD09A1的遥感变量； MOD13Q1和MOD15A2均包括2种植被指数产品，分别为NDVI、EVI和FPAR、LAI，遥感变量均包括原始植被指数和16种纹理因子。

为筛选合适的遥感因子，研究运用ArcGIS 10.2软件提取地面样地所在位置的遥感因子，利用SPSS 20.0软件分析遥感变量和样地森林碳密度值间的Pearson相关性，保留相关性最高的遥感因子参与后续模拟。

1.2.3 土地利用数据 根据研究区实际情况，参照《土地利用现状分类标准》(GB/T 21010—2015)，将研究区土地利用/覆盖分为森林、耕地、水体、草地、建设用地和其他用地6大类。在ENVI 5.0遥感软件支持下，利用先验知识建立训练样本，对处理后的Landsat5遥感影像进行监督分类，得到初始分类结果，然后结合森林资源连续清查数据和森林资源分布图，对监督分类结果进行修正，得到覆盖湖南省和攸县的2009年土地利用/覆盖数据。经检验，总体精度分别为85.1%和89.5%，Kappa系数均超过0.79，精度较高，满足研究的需要。最后提取分类结果中的森林信息，将其作为掩膜数据，用于后续序列高斯协同模拟和森林碳密度制图。

1.3 研究方法

1.3.1 序列高斯协同模拟 序列高斯协同模拟(sequential Gaussian co-simulation, SGCS)是以地面样地和遥感影像为基础，通过结合变异函数，利用随机模拟算法估计未知参数的方法。该算法假设研究区由n个等大小的像元构成，每个像元的估计值从已有估计值和周围样地数据确定的条件累积分布中通过随机抽样获得。该条件分布由一个统计平均数和方差确定，统计平均数和方差可以基于已有的地面样地数据和遥感影像通过点位协同简单克里格估计实现。具体计算公式如下：

(1)

(2)

式中：zsck(u)代表协同简单克里格算法像元位置u处森林碳密度的预测值；δ2(sck)(u)代表协同简单克里格算法的预测方差；z(uα)代表样本数据，α=1,2,…,n(u)；n(u)代表给定搜索范围内获得的样本数量；y(u)代表像元u处的光谱变量；λsckα、λscky分别代表协同简单克里格算法中样本数据和影像数据的权重；mz、my分别代表地面样本数据和遥感影像数据的均值；Czz(0)代表地面样本数据的方差；Czy(0)代表森林碳密度与遥感影像光谱变量的协方差； 且当h=uα-u时，Czy(h)代表森林碳密度预测值与光谱变量的交叉协方差函数。

当输出的森林碳密度分布图像元与固定样地大小一致时，SGCS模拟过程为： 1) 用随机抽样方法设置一种遍历每个像元的顺序； 2) 从分布中随机抽取一个像元位置值u，采用协同简单克里格算法分别估计像元位置u处的预测值和方差，由此预测值和方差确定一个条件累积分布函数； 3) 从分布中随机抽取一个值，将其作为随机变量在像元位置u的实现。重复步骤1～3，直至所有像元都有估计值，进而得到一张覆盖整个研究区的森林碳密度分布图。将该过程执行M次，就可得到M张分布图，通过求和取平均值获取研究区森林碳密度的均值分布图和方差分布图，本研究的估计次数M全部取250。具体程序执行过程，通过修改Wang等(2004a)提出的空间估计算法实现。变异函数是序列高斯协同模拟中最重要的参数之一，研究采用VARIOWIN 3.0(Pannatier, 1996)分析样地森林碳密度，采用Spherical标准化模型进行拟合，从中选择最优变异函数模型。

1.3.2 基于块的序列高斯协同模拟 为解决固定样地大小(25.8 m × 25.8 m)与MODIS空间分辨率(250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m)不一致的问题，本研究采用Wang等(2004a)提出的尺度转换方法——基于块的序列高斯协同模拟(sequential Gaussian block co-simulation, SGBCS)，首先将25.8 m × 25.8 m的地面样地分别上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3种尺度，实现地面样地数据和MODIS像元空间分辨率的匹配； 然后将上推后的样地数据与MODIS数据结合，利用序列高斯协同模拟算法依次开展250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3种尺度的森林碳密度估算； 最后将最优算法用于湖南省森林碳密度估算，探讨其在大区域森林碳密度估算领域的适用性。地面样本的划分同Landsat5数据，2/3数据用于建模，1/3数据用于验证。为深入探讨基于块的序列高斯协同模拟算法的适用性，本研究将其用于湖南省森林碳密度的估算研究，其中3/4(2 169块样地)数据作为建模样本，余下的1/4数据(723块样地)作为验证样本。

研究包括4种形式的森林碳密度模拟： 1) 基于固定样地和攸县Landsat5影像的SGCS(30 m × 30 m)； 2) 基于固定样地和攸县Landsat5影像的SGBCS(250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m)； 3) 基于上推样地和攸县MODIS影像的SGCS(250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m)； 4) 基于上推样地和湖南省MODIS影像的SGCS(500 m × 500 m)。

(3)

式中：cov(ui,uj)表示像元ui和uj的协方差，从该累积分布中，随机抽取一个值作为该块的模拟值；σ2(sck)表示像元ui的协同克里格方差。

基于块的协同克里格模拟中，每个块的条件累积分布取决于块内像元的预测值和克里格方差以及二者之间的协方差。基于块的协同克里格模拟算法克服了已有尺度上推算法的缺陷，可以较好地运用连续分布的遥感影像开展森林碳密度的空间估计研究，算法的效果取决于用户估计的精度、方差和运算时间。

1.3.3 精度评价 基于固定样地的森林碳密度估算在由单木水平、样地水平推算到区域水平的过程中，存在很大的不确定性，忽略这些不确定性将导致区域森林碳密度的高估或低估。本研究采用判定系数(R2)、均方根误差(root mean square error, RMSE)和相对误差(relative error, RE)3个指标，分别对Landsat5和MODIS像元水平的森林碳密度模拟结果进行不确定性分析。

(4)

(5)

(6)

R2可反映估测值与对应实测值之间的趋势线拟合程度； RMSE虽能很好地反映估测模型的可靠性，但其是一个绝对数，大小与均值相关，无法直观反映结果的准确程度；RE具有相对性，不仅考虑了样本估计值与实测值之间误差的大小，同时也兼顾了样本本身的大小，其值越小模型的估测精度越高。

2 结果与分析

2.1 空间变异函数

研究针对攸县25.8 m × 25.8 m、250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 4种尺度的样地数据进行半方差函数拟合(图2)。变程表示空间变异和自相关函数的最大距离，观察图2可知模型的变程分别为15.62、14.91、17.04和17.76 km，在此距离之内，观测值是空间相关的； 在此范围之外，观测值本质上相互独立。

图2 森林碳密度数据的标准化变异函数γ(h)Fig.2 Spatial autocorrelation γ(h) of forest carbon using spherical model for standardized datah表示距离h is distance(km).

其标准化Spherical模型为：

(7)

(8)

(9)

(10)

研究区样地间距在东西向和南北向分别为8 km和4 km，森林碳密度空间变异的动态范围(变程)大于样地间距，样本具有足够的空间信息量。

2.2 样地数据分析

分析研究区4种尺度的样地碳密度统计结果(表2)可知，攸县25.8 m × 25.8 m样地碳密度在0 ～ 39.27 t·hm-2之间变化，标准差为8.39，且变异系数大于1.00，属于强变异，说明攸县森林碳密度数值分布合理，相互之间存在较大差异。值得注意的是，随着样地空间尺度增大，样地碳密度平均值、标准差、变异系数呈逐渐减小的趋势，说明研究区的森林样地空间分布不均匀。

表 2 攸县固定样地碳密度统计结果①

① SD: Standard deviation; CV:Coefficient of variation.

2.3 相关性分析

采用SPSS 20.0软件分别计算样地碳密度(包括25.8 m × 25.8 m固定样地和3种上推样地)与遥感变量间的Pearson相关性，结果发现基于Landsat5和MOD09A1的遥感变量与样地碳密度的相关性较高，分别在-0.455～0.497和-0.626～0.763之间变化，表3列出了相关系数排在前6位的遥感变量。对于Landsat5数据，当显著水平为0.01时，与森林碳密度达到显著相关的变量有1/TM3、TM1var、Elevation、1/TM2等15个因子，其中1/TM1、1/TM2、1/TM3和TM2mean4个变量相关系数在0.450以上； 对于MOD09A1数据，当显著水平为0.01时，与森林碳密度相关系数达到显著的变量有1/Band1、1/Band3、1/Band7、NDVI等28个因子，其中1/Band1、1/Band3 和1/Band4 3个变量相关系数在0.700以上。

表3 样地森林碳密度与遥感变量的Pearson相关性①

① **： 在0.01水平(双侧)上显著相关，*： 在0.05水平(双侧)上显著相关； TM2mean表示TM2灰度共生矩阵的平均值。Correlation is significant at the 0.01 level ** and the 0.05 level*; TM2meanis the average value of gray-level co-occurrence matrix； Band67=Band6/Band7.

MOD13Q1和MOD15A2数据衍生的遥感变量与森林碳密度的相关系数均较低，分别介于-0.230～0.312和-0.390～0.395之间，当显著水平为0.01时，二者与森林碳密度相关系数达到显著的变量分别为NDVImean和LAImean、LAIhom、LAIdis、FPARmean、LAIcon，其中相关系数最高的因子为NDVImean和LAImean。因此，研究依次选择相关系数最高的4个因子： 1/TM3(Landsat5数据)、NDVImean(MOD13Q1数据)、1/Band1(MOD09A1数据)和LAImean(MOD15A2数据)，参与攸县森林碳密度的空间模拟，相关系数分别为0.497、0.312、0.763和0.395，显著水平均达0.000。

2.4 基于Landsat5的森林碳密度模拟

基于固定样地和Landsat5影像，利用序列高斯协同模拟算法导出攸县30 m × 30 m像元水平的森林碳密度分布图，同时采用验证样本对其进行精度评价。结果表明，模型预测值和实测值之间具有良好的线性拟合关系(图3a)，决定系数R2达0.82，估测精度为82.02%，且残差图(图3b)散点分布均匀。张茂震等(2014)基于浙江省仙居县固定样地和Landsat5影像，采用序列高斯协同模拟算法开展了森林碳估算研究，结果表明估算效果较好，进一步佐证了该算法的适用性。

图3 基于Landsat5的攸县森林碳密度模拟值精度验证Fig.3 Accuracy assessments for simulated values of forest carbon density in You county based on Landsat5a.预测值与实测值比较Comparison of predicted and observed values; b.残差分析Residual analysis.

为直观反映模拟结果与地面样地值的吻合程度，研究将森林碳密度模拟结果与地面样地值叠加(图4)，结果发现攸县森林碳密度整体上与地面样地实测值的分布趋势保持一致： 东部由于经营管理完善，分布大面积质量较好的针叶成熟林，森林碳密度高； 中部和西南部主要用于农业种植和城镇建设，森林碳密度相对较低。采用ArcGIS 10.0软件进行统计得出，攸县2009年平均森林碳密度达15.06 t·hm-2，具有较高的实用性和可信度。

研究采用的固定样地大小为25.8 m × 25.8 m，但用于大尺度反演的遥感数据MODIS像元大小分别为250 m × 250 m、500 m × 500 m、1 000 m ×1 000 m，本文在前期数据处理的基础上，通过引入Wang等(2004a)设计的基于块的序列高斯协同模拟，结合Landsat5数据和地面样地数据，将地面调查样地大小从25.8 m × 25.8 m分别上推到250 m × 250 m、500 m × 500 m 和1 000 m × 1 000 m(图5)。基于尺度上推的序列高斯协同模拟程序，通过修改Dr Journel的协同克里格程序和Dr Almeida的协同模拟程序实现(Pannatier, 1996; Deutschetal., 1992)。

为进一步验证Landsat5数据上推结果的有效性，研究采用ArcGIS10.0软件空间分析模块的“波段集统计”工具，分别计算3种尺度上推结果与30 m × 30 m模拟结果重采样后森林碳密度的相关性，结果显示相关系数分别为0.996、0.987和0.964，说明上推效果较好，可用于后续基于MODIS的碳密度估算。

图4 攸县基于Landsat5数据和序列高斯协同模拟算法的森林碳密度空间预测指标与样地实测值的比较Fig.4 Spatial distribution of predicted indexes for forest carbon compared with the plot values at spatial resolution of 30 m×30 m by sequential Gaussian co-simulation using Landsat5 image in You countya.预测均值Predicted values； b.预测均方差Predicted variance； c.预测值大于均值概率Probability for predicted values larger than sample means.

图5 攸县森林碳密度基于块的序列高斯协同模拟预测结果Fig.5 Spatial distributions of simulated forest carbon density using sequential Gaussian block co-simulation algorithm and Landsat5 imagea. 250 m × 250 m; b.500 m × 500 m; c.1 000 m × 1 000 m

2.5 基于MODIS的森林碳密度模拟

利用固定样地点分别提取250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3种水平的森林碳密度，同时结合MOD13Q1、MOD09A1和MOD15A2影像，开展相应空间分辨率的森林碳密度估算。表4列出了3种MODIS数据森林碳密度模拟值的精度检验结果。

分析表4可知，MOD09A1数据的模拟效果最好，R2达0.78，其次是MOD13Q1数据(R2=0.57)，MOD15A2数据效果最差(R2仅0.53)； 均方根误差(RMSE)和估测精度与决定系数(R2)表现一致，其中模拟误差分别为4.65、2.56和4.26 t·hm-2，估测精度分别为58.58%、75.64%和54.84%。3种检验指标共同表明多光谱数据MOD09A1的估算结果明显优于植被指数产品MOD13Q1和MOD15A2，在区域森林碳密度估算方面显示出明显的优势。

同Landsat5数据，研究基于MOD09A1影像的森林碳密度模拟结果包括3部分： 碳密度模拟平均值分布(图6a)、碳密度模拟平均方差分布(图6b)、碳密度模拟值大于均值的概率分布(图6c)，空间分辨率为500 m × 500 m。

表4 3种MODIS数据森林碳密度模拟值的精度检验结果

图6 攸县基于MOD09A1影像和序列高斯协同模拟算法的森林碳密度空间预测指标分布Fig.6 Spatial distribution of predicted indexes for forest carbon density using the sequential Gaussian co-simulation algorithm based on MOD09A1 image in You Countya.预测均值Predicted values; b.预测均方差Predicted variances; c.预测值大于均值概率Probability for predicted values larger than sample means.

分析可知，基于MOD09A1数据模拟的森林碳密度空间分布(图6a)与研究区内基于Landsat5模拟的森林碳密度具有较好的一致性，即样地模拟值较高的局部区域，相应MODIS数据的模拟值也较高。这说明利用MODIS数据结合基于块的序列高斯协同模拟算法，可以真实再现研究区森林碳密度的空间分布格局，模拟结果的方差分布(图6b)和模拟结果大于平均数的概率分布(图6c)，与森林碳密度模拟结果的平均值具有较好的空间一致性。

2.6 湖南省森林碳密度模拟

为进一步探讨基于块的序列高斯协同模拟算法的适用性，研究将其用于湖南省森林碳密度的估算研究。首先借助ArcGIS10.0软件的sample工具，抽取研究区3/4数据(2 169块样地)开展基于Landsat5影像和固定样地的森林碳密度估算研究； 然后采用预留的1/4数据(723块样地)进行精度验证(图7)。结果表明，模型预测值和实际观测值之间具有良好的线性拟合关系，决定系数(R2)达0.84，估测精度为87.02%，且残差散点分布均匀，说明模型的拟合效果较好，可用于后续基于MOD09A1的森林碳密度估算。

在前期数据处理的基础上，采用基于块的序列高斯协同模拟算法，将湖南省固定样地从25.8 m × 25.8 m上推到500 m × 500 m，利用固定样地点提取相应水平的森林碳密度值(图8a)。同时结合MOD09A1影像，采用序列高斯协同模拟算法开展湖南省森林碳密度的估算研究； 样本数据的划分同Landsat5。

采用决定系数(R2)、均方根误差(RMSE)和估测精度(EA)3个指标检验湖南省森林碳密度的估算精度，具体检验结果见表5。分析可知，基于MOD09A1数据的湖南省森林碳密度模拟效果较好，决定系数(R2)达0.75，估测精度为74.06%，进一步佐证了基于块的序列高斯协同模拟算法的有效性，可用于区域森林碳密度的估算。

图7 湖南省森林碳密度样地模拟值精度验证Fig.7 Accuracy assessments for simulated values of forest carbon density in Hunan Provincea.预测值与实测值比较Comparison of predicted and observed values; b.残差分析Residual analysis.

图8 湖南省森林碳密度空间分布Fig.8 Spatial distribution map of forest carbon density in Hunan Provincea.上推样点空间分布Spatial distribution for up-scaling plots using spatial block co-simulation algorithm;b.模拟值空间分布 Spatial distribution for simulated values using sequential Gaussian co-simulation algorithm.

表5 基于MOD09A1数据的湖南省森林碳密度精度检验结果

图8b显示了基于MOD09A1影像的湖南省森林碳密度分布，分析可知，湖南省森林碳密度与地面样地上推值保持一致的分布趋势： 高值区主要分布在湘东、湘南的郴州、湘西的怀化、吉首、张家界等地，湘中地区的森林碳密度较小。因为研究区东、南、西三面环山，中部分割为丘陵性河谷盆地，北部是沉积平原和吞噬湖泊，特殊的地形条件使得湘东、西、南三面山区分布着大面积的森林植被，由于经营管理完善，森林质量较好，碳密度大； 湘北和湘中地区主要用于农田种植和水体覆盖，森林碳密度相对较小。沅江以及资水附近的大面积连续森林碳密度高值分布，充分反映了森林植被随地表水分布的规律； 而中南部和东北部洞庭湖附近的大面积连续低值分布，则反映出人类活动对森林植被覆盖的影响，森林碳密度整体上与研究区植被分布规律一致。该结论进一步说明利用基于块的序列高斯协同模拟算法可以再现湖南省森林碳密度的空间分布格局。

3 讨论

3.1 相关性分析

作为反映变量之间相关性密切程度的统计指标，森林碳密度和遥感因子之间的相关系数与遥感影像的空间分辨率和传感器有关。本研究中4种遥感数据(Landsat5、MOD13Q1、MOD09A1和MOD15A2)衍生的变量均与森林碳密度具有较高的相关性，且相关系数均在0.312以上。Landsat5衍生变量与森林碳密度的相关系数在-0.455～0.497之间波动，且1/TM3与森林碳密度的相关系数最高(达0.497)，这与Wang等(2014)的研究结果一致(相关系数0.312)； 比较而言，MOD09A1衍生变量与森林碳密度的相关系数在-0.626～0.763之间波动，其中1/Band1与森林碳密度的相关系数最高。从某种程度上说，这一结论与Landsat5数据的研究结果相似，主要原因是Landsa5的 TM3波段与MOD09A1的Band1波段具有相似的波谱区间，分别在0.63～0.69 μm和0.62～0.67 μm之间波动。此外，MOD13Q1数据的衍生变量NDVImean和MOD15A2数据的LAImean分别与森林碳密度的相关性最高，因此二者分别用于250 m × 250 m和1 000 m × 1 000 m 2种空间分辨率的森林碳密度制图。

3.2 不同尺度森林碳密度模拟

为克服地面样地大小和遥感影像空间分辨率不一致的缺陷，研究采用Wang等(2004a)设计的基于块的序列高斯协同模拟算法，将样地碳密度从25.8 m × 25.8 m分别上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3种尺度，结果发现该算法可以很好地刻画森林碳密度的空间分布格局。同250 m × 250 m的MOD13Q1和1 000 m × 1 000 m的MOD15A2产品相比，多光谱数据Landsat5和MOD09A1在森林碳密度模拟方面显示出明显的优势，因为来自可见光、近红外和中红外波段的光谱变量可以提取指示森林碳密度信息的光谱变量(Chenetal., 2002)。Landsat5数据的模拟结果与样地实测值之间具有很好的一致性，且由其和固定样地数据上推的250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m 3种尺度的森林碳密度精度可靠，可作为相应尺度森林碳密度模拟的真值数据。虽然由MOD09A1数据模拟的森林碳密度在区域和国家尺度是可以接受的，但不能用于局部地区森林碳密度的高精度模拟，主要是因为MOD09A1数据在刻画空间细节方面不如Landsat5数据精细； 另外，受Landsat5影像质量和序列高斯协同模拟算法数据量的限制，本文只是通过随机抽样获取部分样地上推到500 m × 500 m，然后将上推后的样地与MOD09A1数据结合，开展湖南省森林碳密度的空间模拟，利用Landsat5和固定样地数据开展湖南省森林碳密度的高精度模拟将是下一步的研究重点。

3.3 不同信息源森林碳密度模拟

研究采用基于块的序列高斯协同模拟算法，将地面样地(包括250 m × 250 m的地面样地和3种上推后的地面样地)分别与不同空间分辨率的遥感影像结合，开展多种空间尺度的森林碳密度模拟，该算法克服了地面样地大小与遥感影像空间分辨率不一致的缺陷。虽然MODIS数据提供的模拟结果可用于大尺度的森林碳密度制图，但会损失部分森林碳密度的估计精度和空间变异性。比较而言，Landsat5数据在模拟森林碳密度的空间分布和变异性方面更为精确，比MOD09A1数据呈现的森林冠层结构更为复杂，这与Yang等(2007)的研究结论一致。另外，Landsat5的TM3波段和MOD09A1的Band1波段在提供森林碳密度的预测变量方面显示出较大的潜力。

3.4 森林碳密度模拟精度的影响因素分析

正确估算森林碳密度是研究全球气候变化的关键途径。本研究应用序列高斯协同模拟算法对攸县和湖南省森林碳密度的空间分布进行了估算，同利用光谱值和森林参数建立的回归模型相比，序列高斯协同模拟算法能够更多地考虑森林参数同光谱值之间的非线性依赖关系(Trotteretal., 1997)。SGCS是基于像元的估计方法，且模拟结果的最小值接近0，更能真实反映非林地的实际情况(沈希等, 2011)。森林碳密度的模拟精度受多种因素的影响，如模拟次数L、影像特征波段的选取和采样密度等：1) 相关研究表明，模拟次数L对模拟结果精度有重要影响，模拟次数越多，其平均值就越接近该位置的数学期望(郭含茹等, 2016)。本文当L达到250后，各像元位置森林碳密度的估计值趋于稳定，因此将其作为最终的模拟次数。2) 由于影响SGCS森林碳密度最优波段的选择往往随研究区和数据源的不同而变化，因此最优波段需要通过具体试验才能获取。本文采用的Landsat5和3种MODIS数据，虽然研究区相同，但是最优波段因数据源的不同而变化。3) 众所周知，随着采样密度的增加，抽样误差随之减小，但到达一定程度后，逐渐趋于稳定。本文主要研究森林资源清查样地与不同传感器数据相结合的模拟，然而关于不同采样密度下的森林碳密度估计考虑不多。郭含茹等(2016)通过分析4种采样密度对浙江省仙居县森林碳密度模拟精度的影响，指出利用SGCS、SGBCS和森林资源清查样地开展森林碳密度估计时，可适当降低对采样密度的要求。

4 结论

随着遥感技术的发展，基于固定样地数据的森林碳密度遥感制图得到了广泛应用，然而地面样地大小和遥感影像空间分辨率不一致的缺陷阻碍了制图精度的提高。鉴于此，本研究以地面固定样地和Landsat5影像为数据源，采用基于块的序列高斯协同模拟算法，将25.8 m × 25.8 m的地面样地分别上推到250 m × 250 m、500 m × 500 m和1 000 m × 1 000 m； 然后将上推样地和MODIS影像结合，探讨MODIS数据在区域森林碳密度遥感估算领域的潜力。结果表明： 1) 采用基于块的序列高斯协同模拟算法，可以实现由地面样地到不同空间分辨率MODIS像元之间的转换； 2) Landsat5和MODIS数据与森林碳密度的敏感因子具有高度的相似性，排在前3位的分别为TM3、TM2、TM1和Band1、Band4、Band3； 3) 同植被指数产品MOD13Q1和MO15A2相比，多光谱数据Landsat5和MOD09A1在模拟森林碳密度方面显示出巨大潜力，因为来自可见光、近红外和中红外波段的光谱变量可以提取指示森林碳密度信息的遥感因子，模拟精度分别82.02%和75.64%，但受空间分辨率的限制，MOD09A1数据在刻画空间细节方面不如Landsat5精细； 4) 为进一步探讨基于块的序列高斯协同模拟方法的适用性，研究将其推广到基于MOD09A1数据的湖南省森林碳密度模拟，结果表明估算效果较好，估测精度为74.06%，可用于区域森林碳密度的模拟。

崔丽华, 刘善军, 张艳博, 等. 2009. IDL程序去除MODIS遥感影像镶嵌中的黑色条带. 河北理工大学学报： 自然科学版，31(3): 79-82.

(Cui L H, Liu S J, Zhang Y B,etal. 2009. Using IDL to remove black belt from the mosaic of MODIS remote sensing images. Journal of Hebei Polytechnic University： Natural Science, 31(3): 79-82. [in Chinese])

郭含茹, 张茂震, 徐丽华, 等. 2016. 不同采样密度下县域森林碳储量仿真估计. 生态学报, 36(14): 4373-4385.

(Guo H R, Zhang M Z, Xu L H,etal. 2016. Simulation of regional forest carbon storage under different sampling desities. Acta Ecologica Sinica, 36(14): 4373-4385. [in Chinese])

李海奎. 2010. 中国森林植被生物量和碳储量评估. 北京： 中国林业出版社.

(Li H K. 2010. Estimation and evaluation of forest biomass carbon storage in China. Beijing: China Forestry Publishing House. [in Chinese])

刘 畅, 李凤日, 甄 贞. 2014. 空间误差模型在黑龙江省森林碳储量空间分布的应用. 应用生态学报, 25(10): 2779-2786.

(Liu C, Li F R, Zhen Z. 2014. Prediction of spatial distribution of forest carbon storage in Heilongjiang Province using spatial error model. Chinese Journal of Applied Ecology, 25(10): 2779-2786. [in Chinese])

戚玉娇, 李凤日. 2015. 基于KNN 方法的大兴安岭地区森林地上碳储量遥感估算. 林业科学, 51(5): 46-55.

(Qi Y J, Li F R. 2015.Remote sensing estimation of aboveground forest carbon storage in Daxing’an Mountains based on KNN method. Scientia Silvae Sinicae, 51(5): 46-55. [in Chinese])

冉有华, 李 新. 2009. 基于块克里金的土壤水分点观测向像元尺度的尺度上推研究. 冰川冻土, 31(2): 275-283.

(Ran Y H, Li X. 2009. Up Scaling of point soil moisture measurements to oixel averages based on block Kriging. Journal of Glaciology and Geocryology, 31(2): 275-283.[in Chinese])

沈 希, 张茂震, 祁祥斌. 2011. 基于回归与随机模拟的区域森林碳分布估计方法比较. 林业科学, 47(6): 1-8.

(Shen X, Zhang M Z, Qi X B. 2011. Comparison of regional forest carbon estimation methods based on regression and stochastic simulation. Scientia Silvae Sinicae, 47(6): 1-8.[in Chinese])

吴 丹, 邵全琴, 刘纪远, 等. 2011.1985—2030年江西泰和县森林植被碳储量的时空动态. 应用生态学报, 2011, 22(1): 41-46.

(Wu D, Shao Q Q, Liu J Y,etal. 2011. Spatiotemporal dynamics of forest carbon storage in Taihe County of Jiangxi Province in 1985—2030. Chinese Journal of Applied Ecology, 22(1): 41-46. [in Chinese])

严恩萍, 林 辉, 王广兴, 等. 2015. 基于MODIS混合像元分解的湖南省森林碳密度反. 应用生态学报, 26(11): 3433-3442.

(Yan E P, Lin H, Wang G X,etal. 2015.Estimation of human forest carbon density based on spectral mixture analysis of MODIS data. Chinese Journal of Applied Ecology, 26(11): 3433-3442. [in Chinese])

张茂震, 王广兴, 葛宏立, 等. 2014. 基于空间仿真的仙居县森林碳分布估算. 林业科学, 50(11): 13-22.

(Zhang M Z, Wang G X, Ge H L,etal. 2014. Estimation of forest carbon distribution for Xianju county based on spatial simulation. Scientia Silvae Sinicae, 50(11): 13-22. [in Chinese])

张茂震, 王广兴, 周国模, 等. 2009. 基于森林资源清查、卫星影像数据与随机协同模拟尺度转换方法的森林制图. 生态学报, 29(6): 2919-2928.

(Zhang M Z, Wang G X, Zhou G M,etal. 2009. Mapping of forest carbon by combining forest inventory data and satellite images with co-simulation based up-scaling method. Acta Ecologica Sinica, 29(6): 2919-2928. [in Chinese])

Chen J M, Pavlic G, Brown L,etal. 2002. Derivation and validation of Canada-wide coarse-resolution leaf area index maps using high-resolution satellite imagery and ground measurements. Remote sensing of environment, 80(1): 165-184.

Deutsch C V, Journel A G. 1992. GSLIB: geostatistical software library and user’s guide. 2nded. Oxford University Press.

Lu D S, Chen Q, Wang G X,etal. 2012. Aboveground forest biomass estimation with Landsat and lidar data and uncertainty analysis of the estimates. International Journal of Forestry Research, Article ID 436537, 16 pages, doi:10.1155/2012/436537.

Pannatier Y. 1996. Variowin: software for spatial statistics analysis in 2D. New York, Springer.

Tian X, Su Z, Chen E,etal. 2012. Estimation of forest above-ground biomass using multi-parameter remote sensing data over a cold and arid area. International Journal of Applied Earth Observation and Geoinformation, 14(1): 160-168.

Trotter C, Dymond J, Goulding C. 1997. Estimation of timer volume in a coniferous planation forest using Landsat TM. International Journal of Remote Sensing, 18(10): 2209-2223.

Wang G X, Gertner G Z, Anderson A B. 2004a.Spatial-variability-based algorithms for scaling-up spatial data and uncertainties. IEEE Transactions on Geoscience and Remote Sensing, 42(9): 2004-2015.

Wang G X, Gertner G Z, Anderson A B. 2004b. Mapping vegetation cover change using geostatistical methods andbi-temporal Landsat TM images. IEEE Transactions on Geoscience and Remote Sensing, 42(3): 632-643.

Wang G X, Tonny O, Zhang M Z,etal. 2009. Mapping and spatial uncertainty analysis of forest vegetation carbon by combining national forest inventory data and satellite images. Forest Ecology and Management, 258(7): 1275-1283.

Wang G X, Zhang M Z. 2014. Upscaling with conditional co-simulation for mapping above-ground forest carbon. Scale Issues in Remote Sensing,doi:10.1002/9781118801628.ch06.

Xiang H B, Liu J S, Cao C X,etal. 2013. Algorithms for moderate resolution imaging spectroradiometer cloud-free image compositing. Journal of Applied Remote Sensing, 7(1): 073486.doi:10.1117/1.JRS.7.073486.

Yan E P, Lin H, Wang G X,etal. 2015. Improvement of forest carbon estimation by integration of regression models and spectral unmixing of Landsat data. IEEE Geoscience and Remote Sensing Letters, 12(9): 2003-2007.

Yang P, Ryosuke S, Wu W B,etal. 2007. Evaluation of MODIS land cover and LAI products in cropland of North China Plain using in situ measurements and Landsat images. IEEE Transactions on Geoscience and Remote Sensing, 45(10): 3087-3097.

Zhao P, Lu D, Wang G,etal. 2016. Examining spectral reflectance saturation in Landsat imagery and corresponding solutions to improve forest aboveground biomass estimation. Remote Sensing, doi:10.3390/rs8060469.

(责任编辑 石红青)

Estimation of Forest Carbon Density Based on Geostatistics and Multi-Resource Remote Sensing Data

Yan Enping1,2Zhao Yunlin1,2Lin Hui1,2Mo Dengkui1,2Wang Guangxing1,2,3

(1.KeyLaboratoryofForestryRemoteSensingBigData&EcologicalSecurityforHunanProvinceCentralSouthUniversityofForestry&TechnologyChangsha410004; 2.CollegeofForestry,CentralSouthUniversityofForestry&TechnologyChangsha410004; 3.DepartmentofGeography,SouthernIllinoisUniversityCarbondaleIl629012USA)

【Objective】 Due to the inconsistency of spatial resolutions between sample plots and image pixels, in this study the estimation of forest carbon density was conducted by using geostatistics method and multi-resource remote sensing data, which aims to provide reference for the application of MODIS data in the regional estimation of forest carbon density.【Method】 Firstly, the spatial block co-simulation algorithm was employed to scale up the sample plots of forest carbon density in You county of Hunan Province from the spatial resolution of 25.8 m×25.8 m to the spatial resolutions of 250 m × 250 m, 500 m × 500 m and 1 000 m × 1 000 m respectively. Then, MODIS images with three spatial resolutions corresponding to those mentioned above, were applied to map forest carbon density for this county using sequential Gaussian co-simulation algorithm. Finally, the best model was applied in the estimation of forest carbon density for Hunan Province.【Result】 There were highly similarities for sensitive factors of forest carbon density between Landsat5 and MOD09A1 data, according to theresult of Pearson product moment correlations, the top three sensitive factors were 1/TM3, 1/TM2, 1/TM1 for Landsat5 and 1/Band1, 1/Band4, 1/Band3 for MOD09A1 respectively; compared to the vegetation product of MOD13Q1 and MOD15A2, multi-spectral data such as Landsat5 and MOD09A1 showed great potential in the simulation of forest carbon density with the accuracy of 82.02% and 75.64%, respectively; there is good application for the sequential Gaussian co-simulation algorithms based on the image of MOD09A1, which can be used in the spatial simulation of forest carbon density for Hunan Province with the accuracy of 74.07%.【Conclusion】 The spatial block co-simulation algorithm could be used to realize the conversion of spatial resolutions from sample plots to the pixels of MODIS images. It was also found that the MODIS derived maps were more smoothed than those from Landsat5 due to the limitation of spatial resolutions, especially in the terms of capturing the spatial variability of forest carbon density. The adopted method was well suited for regional estimation of forest carbon density based on the combination of forest inventory sample plot data and remotely sensed images, especially for the areas that the plot sizes and images pixels were inconsistent.

forestry remote sensing； forest resource inventory； multi resource remote sensing； sequential Gaussian block co-simulation； forest carbon density

10.11707/j.1001-7488.20170708

2016-03-21；

2016-09-15。

国家“十二五”高技术发展研究计划项目(2012AA102001)； 国家自然科学基金面上项目(31470643)； 林学重点学科开放基金拟资助项目(2016YB08)。

S718.5; S757.2

A

1001-7488(2017)07-0072-13

*林辉为通讯作者。