区域森林碳储量及其分布已成为近年来森林生态学及其相关领域研究的热点(王效科等,2001; 杨洪晓等,2005; 黄东,2008; 黄从德等,2008)。对区域森林CO₂源汇功能进行评价需要森林碳储量分布信息,将来碳工业及碳市场的发展也要求提供区域、国家乃至全球尺度的森林碳分布(Wang et al.，2009)。近年来,相关研究主要集中在森林生物量和碳储量方面(聂道平等,1997; 王秀云等,2008; 刘爽,2009),森林碳分布研究仍然很少(吴仲民等,1998)。在REDD(reducing emissions from deforestation and degradation)背景下,由于区域森林碳分布存在不均匀性,且森林管理对碳储量和密度具有重要影响(Pussinen et al.,2002;Fahey et al.,2010),及时准确地获取区域森林碳储量及其分布信息具有重要意义(魏安世等,2006)。

区域森林碳分布的估计是一个空间估计问题。传统的统计学提供了总体特征值的估计方法,地统计学通过空间自相关理论提供了描述这些特征量空间分布的方法,在理论上解决了空间估计基础层面的问题(Wang,1993)。克里格方法和条件模拟算法是地统计学中进行空间插值的主要手段。协同克里格算法可以将各种不同类型、不同可靠程度的资料结合在一起进行线性回归,它是一种求最优、线性、无偏内插估计量的方法(Wang et al.，2009;Fuchs et al.，2009;Muukkonen et al.，2007)。森林碳估计值是一个区域化变量,而克里格方法对其具有较强的平滑效应,地统计学空间随机模拟法能克服克里格法的缺陷(Peter et al.，1982),它将数据作为一个整体来复原其整体的空间结构,追求的是模拟的真实性,尽可能地接近真实的空间分布,不像克里格方法,追求的是特定点位某个属性的局部最优估值(史舟等,2007)。已有研究表明,空间随机模拟方法在模拟森林碳分布方面具有较好效果(Wang,2009;张茂震等,2009;刘平等,2007)。

除空间随机模拟方法外,回归估计也是一种区域森林碳储量估计方法,基于遥感影像的像元进行回归估计,也可以实现区域森林碳空间分布估计。本研究选择以低山丘陵地貌为主的临安市为研究对象,分别采用一元二次非线性回归(回归估计)和空间随机模拟(仿真模拟)方法对研究区森林地上部分碳储量和密度(本文中提到的森林碳储量和碳密度,除特别指出外,均指森林地上部分碳储量和密度)进行空间估计,并对照地面样地调查结果对2种模拟结果在空间上的差异性进行分析,进而对2种方法作出评价。

1 研究区概况

临安市(118°51'—119°52'E,29°56'—30°23'N)属于中亚热带季风气候区,全年平均降水量1 628.6 mm,多年平均气温15.8 ℃,历年平均日照1 939 h,全市地势西北高东北低,海拔相差1 770 m,西北、西南部山区平均海拔1 000 m以上,而东部临近杭州部分是海拔50 m以下的河谷平原。西、北和南三面环山,东西长达100 km,气候在垂直方向上差异悬殊。由于水系的作用,形成了昌化、於潜、临安三大河谷平原,串珠状河谷地呈东西向和东北—西南向延伸。全市行政辖区总面积3 126.8 km²,林业用地面积26.104 7万hm²,活立木总蓄积量830.010 1万m³,森林覆盖率76.55%,其中有林地面积24.038 9万hm²,有林地蓄积量8 230.567 6万m³。森林资源分布格局为西高东低。植被区划属亚热带东部常绿阔叶林区,为浙皖山地青冈栎(Cyclobalanopsis glauca)和苦槠(Castanopsis sclerophylla)栽培区。

2 研究方法 2.1 数据来源 2.1.1 Landsat TM遥感影像数据

临安市TM遥感影像由119 /39和120 /392景组成,其主要区域在119 /39景内。受遥感数据质量限制,本次取东部119 /39景所覆盖的范围作为研究区,该范围约占全市行政区面积的80%(图 1)。研究区2004年Landsat TM影像数据包括6个波段,即波段1,2,3,4,5和7,空间分辨率均为30 m。几何校正后影像几何误差小于1个像元。辐射校正采用DOS校正法。在此基础上,提取样地位置所对应的6个波段影像亮度值,并进行不同组合比值运算,分析影像亮度值与样地森林碳之间的相关性,选取相关系数最高的波段4,5和7,用波段比值TM4/(TM5+TM7) 参与2类模型的模拟。

2.1.2 地面抽样调查数据

本研究采用2004年临安市森林资源地面抽样调查数据。该调查以整个临安市行政辖区为总体,按系统抽样方法、以1km×3 km的间距设置正方形样地,样地面积0.08 hm²(28.28 m × 28.28 m),调查因子包括地貌特征、土壤类型、植被类型以及与森林生物量相关的树种、胸径等测树因子。实际调查有效样地968个,其中林分样地506个,样地最大蓄积量25.37m³。由于西部约占全市总面积1/6的区域遥感数据缺失,本研究仅在东部地区,该区域的地面样地总数为778个。

2.2 森林碳分布估计方法 2.2.1 一元二次非线性回归

非线性回归是一种以假定因变量与自变量之间存在某种非线性规律而模拟因变量与自变量之间关系的方法。一般来说,非线性回归方法比线性回归方法能更正确地反映客观现象之间的相互联系,而一元二次非线性回归方法既有非线性的特点、又有复杂度低的优点。

本研究假设森林碳密度与Landsat TM遥感影像波段比值TM4/(TM5+TM7) 呈二次关系,建立模型C=av²+bv+c,C为碳密度,v为TM遥感影像4,5和7波段比值: TM4/(TM5+TM7),a,b和c分别为待估参数。

用地面样地数据和与之相对应的TM影像数据进行拟合,得到用于像元水平森林碳密度估计(30m×30m)的一元二次非线性回归(回归估计)模型,C=6 050.136 2v²-1 460.697 3v+152.531 0,最后利用该模型和研究区TM遥感影像第4,5和7波段数据,估算整个研究区碳密度分布,其分布图的分辨率为30 m × 30 m。

2.2.2 序列高斯协同模拟

序列高斯协同模拟(仿真模拟)以研究区地统计变异函数为基础,结合地面样地和TM遥感数据,用随机模拟方法估计森林碳密度,即导出每个像元(块)的森林碳密度估计值。序列高斯协同模拟算法假设:每个像元(块)的碳密度估计值是一个随机函数在该位置的随机变量Z(u)的实现,其过程是从周围的样地数据以及已有估计值所确定的条件累积分布中的一次随机抽样。典型的序列高斯协同模拟过程为: 1) 对于全部像元(块),随机设置一个估计顺序; 2) 在每个像元(块)的位置u,用同位协同简单克里格方法计算条件累积分布函数的均值和方差; 3) 在此条件累积分布中随机抽取一个值,作为随机变量在此像元或块位置u的实现。重复步骤1) 到3),直至所有像元(块)都有碳密度估计值,最后获得整个研究区森林碳密度分布图。将这一过程重复执行L次,就可以得到L个分布图,最后计算森林碳密度分布图的均值图和方差图。

序列高斯协同模拟中最重要的参数之一是变异函数,本研究采用VARIOWIN(Pannatier,1996)分析样地森林碳密度,其Spherical标准化模型如下:

式中: γ _s(| h |)为标准半方差,h为距离。

如图 2所示,模型的变程参数为2 739.6 m,代表空间变异和自相关函数的最大距离,在此范围之内,观测值是空间相关的,在此范围之外,观测值在本质上相互独立。图中在1 500 m左右和2 200 m左右分别有一个点半方差较大,表明在这2点间隔距离上的样地地上部分森林碳密度差异波动较大,可能是随机因素所致,仍属正常。

根据变异函数模型、TM影像数据和地面样地地上部分碳密度,在TM影像像元(30 m × 30 m)水平上进行序列高斯协同模拟,得到全市森林碳密度及其分布。模拟程序参照GISLab的FORTRAN版本编写。

2.3 样地碳储量计算

用单株生物量模型计算样地单株地上部分生物量,累加获得样地森林总生物量。由于目前文献报道的模型覆盖的树种有限,本研究将样木划分为杉木(Cunninghamia lanceolata)、松类(Pinus spp.)、硬阔和软阔4个树种(组),分别用与其对应的模型计算样地森林地上部分总生物量(王震等,2006; 方精云等,2007; 魏文俊等,2007; 焦秀梅等,2005; 孙世群等,2008)。并按碳/生物量的平均转换系数0.5将样地总生物量转换为碳储量,即得样地森林地上部分碳储量,根据样地面积将其转换为碳密度。

3 估计结果统计分析

基于仿真模拟和回归估计2种方法的森林地上部分碳分布如图 3所示,其块大小为30 m × 30 m。基于回归估计方法得研究区森林地上部分碳储量为2 365 404.37 t,最大碳密度为73.714 4 t·hm^-2,最小碳密度为0.715 6 t·hm^-2,平均碳密度为9.000 0 t·hm^-2,标准差为0.380 6。基于仿真模拟方法得研究区森林地上部分碳储量为3 291 659.83 t,最大碳密度为78.913 3 t·hm^-2,最小碳密度为0.083 3 t·hm^-2,平均碳密度为12.523 3 t·hm^-2,标准差为0.639 4。

根据2004年森林资源清查样地数据,按随机抽样方法估计研究区森林碳储量为2 708 897.90 t,样地碳密度平均值为10.306 5 t·hm^-2,其最大值为96.962 5 t·hm^-2,最小值为0.000 0 t·hm^-2; 在像元(30 m × 30 m)水平上分别对仿真模拟和回归估计结果进行特征数统计,并与基于地面样地按随机抽样方法估计(抽样估计)的结果对比,结果见表 1。

由表 1可知,仿真模拟结果和回归估计结果在总量上分布于基于样地实测数据估计结果的两边。回归估计较仿真模拟结果森林碳总量少28.14%。若以基于样地实测数据的估计结果为参照,则回归估计结果的相对差值为-12.7%,仿真模拟的相对差值为21.5%。以仿真模拟结果为参照,回归估计平均碳密度低3.523 3 t·hm^-2,最大碳密度低5.198 9 t·hm^-2,而最小碳密度高0.632 2 t·hm^-2,说明回归估计结果较仿真模拟结果碳密度分布范围变窄,均值下移,下限上移,整体上碳密度分布较平坦,如实反映无林地、非林地等无森林生物量地类的能力较弱。残差分析显示,序列高斯协同模拟残差呈以-3.94为对称轴的正态分布,残差平方和为130 650.274;回归估计呈以0为对称轴的正态分布,其残差的平方和为158 023.884(图 4)。与实测结果对比及残差统计分析结果均表明,仿真模拟结果更接近实际情况。由表 1可见,抽样估计标准差15.160 6与仿真模拟和回归估计标准差相差悬殊,其原因是样本单元数相差太大,一个是全部像元数,一个是地面样地数。

为了从总体上全面分析2种估计方法的特点,将2种方法估计结果进行差运算(像元水平)。图 5为序列高斯协同模拟与一元二次非线性回归碳密度差值图,由像元水平上对应地理位置的森林碳密度估计值相减得到。差值图的值域为-48.710 0 ~ 49.544 4,基本上以0为中心,说明总体上偏差不大,但正负方向的极值差距不小。在碳密度分布上,2种方法的值域不同,仿真模拟的估计值分布范围宽,能更好地表现多林地区和少林及无林地区的实际情况。

对碳密度差值图(图 5)的单位地块(30 m × 30 m)进行统计分析,结果见表 2。2种方法的碳密度差值的平均数为3.524 4 t·hm^-2,为仿真模拟平均碳密度的28.14%,为回归估计的39.16%。对应单位地块差值的最大值为49.544 4 t·hm^-2,最小值为-48.710 0 t·hm^-2。差值大的区域均为森林密度大、生长好的区域,说明回归估计低估了森林碳密度。差值最小的区域,一般为无森林植被区,说明回归估计高估了这些区域的森林碳密度。

估计结果与地面样地实测值比较,是判断估计结果正确性的主要方法之一。将森林碳分布的回归估计和仿真模拟结果图分别与样地森林碳分布图叠加,得森林碳密度比较图(图 6,7)。图 6和7表明,2种方法模拟结果的碳密度分布与该范围内样地实测值具有较好的一致性。虽然2种估计方法均可以在一定程度上真实表现森林碳密度的空间分布格局,但相对而言,回归估计结果平均碳密度较小,为0.715 6 ~ 73.714 4 t·hm^-2,区间跨度为72.998 9 t·hm^-2。仿真模拟结果平均碳密度为0.083 3 ~ 78.913 3 t·hm^-2。仿真模拟结果平均碳密度及碳密度区间跨度均大于回归估计结果,各种无森林植被地类的碳密度接近零,更能体现不同地类森林碳密度差异。

将地面样地调查数据与同位置地块上的估计结果进行对比,可进一步分析不同方法估计碳分布的特征差异。从2个估计结果中分别提取与地面样地位置相对应地块(30 m × 30 m)的碳密度估计结果,以地面样地数据为真值,分别就2种方法的估计结果进行比较,主要以离差平方和来反映2种方法估计结果与相同位置上地面样地数据的差异,其统计量详见表 3。

表 3显示,仿真模拟结果在最大值与最小值之间的分布区间较回归估计结果有明显上移,其差值平方和以及平方和均值显示2种方法有明显差异,仿真模拟相对于回归估计的RMSE更小。

2种方法所得结果与实测值之间的相关性,虽然并不明显集中于y = a+bx函数附近,但整体趋势仍明显存在。实测值中碳密度为0的点过多,这在一定程度上影响了结果的准确性和相关性。一元二次非线性回归与样地观测值之间相关系数为33.93%,序列高斯协同模拟相关系数为54.17%,可见序列高斯协同模拟结果更接近地面样地实测值。

4 结论与讨论

利用森林资源清查数据和遥感影像数据进行区域森林碳分布估计,仿真模拟和回归估计都是可选方法。仿真模拟和回归估计在森林碳总量上均能得到较为合理的结果,但前者在碳密度分布上较后者有更宽的变化幅度,其最大值大于后者,最小值小于后者,森林碳总体总量及其分布均更接近实际。仿真模拟结果的最小值接近0,更能反映非林地及无覆盖地类的实际情况。回归估计结果的森林碳密度分布范围较窄,反映从森林到无植被地类的连续变化的能力较弱。

用仿真模拟和回归估计进行森林碳分布估计,方法上存在本质区别。回归估计方法利用样本拟合模型,再将模型用于估计总体中的每一个单元,估计过程中除了自变量以外没有新的信息加入。仿真模拟是一种仿真方法,它将样地实测值作为真值,根据已知点和周围的遥感数据计算待估点的条件概率分布,从该分布中随机抽取一个值作为该点的模拟实现。模拟过程与位置有关,充分利用了局部信息。可以说,2种方法的本质区别决定了本研究的结果。

2种方法都用到了样地实测数据与遥感影像数据。样地大小为28.28 m × 28.28 m,遥感图像分辨率为30 m × 30 m,匹配有一定误差,几何校正也存在一定偏差,这些在一定程度上影响了模拟结果的精确性和比较的准确性。