林业科学  2014, Vol. 50 Issue (11): 13-22   PDF    
DOI: 10.11707/j.1001-7488.20141102
0

文章信息

张茂震, 王广兴, 葛宏立, 徐丽华
Zhang Maozhen, Wang Guangxing, Ge Hongli, Xu Lihua
基于空间仿真的仙居县森林碳分布估算
Estimation of Forest Carbon Distribution for Xianju County Based on Spatial Simulation
林业科学, 2014, 50(11): 13-22
Scientia Silvae Sinicae, 2014, 50(11): 13-22.
DOI: 10.11707/j.1001-7488.20141102

文章历史

收稿日期:2013-08-25
修回日期:2014-9-31

作者相关文章

张茂震
王广兴
葛宏立
徐丽华

基于空间仿真的仙居县森林碳分布估算
张茂震1, 2, 王广兴3, 葛宏立1, 2, 徐丽华1, 2    
1. 浙江农林大学环境与资源学院 临安 311300;
2. 浙江农林大学 浙江省森林生态系统碳循环与固碳减排重点实验室 临安 311300;
3. Department of Geography and Environmental Resources, Southern Illinois University at Carbondale Carbondale 62901
摘要:以浙江省仙居县为研究区,基于2008年森林资源二类调查样地(清查样地)数据和Landsat TM影像,用序列高斯协同仿真方法模拟全县森林碳储量及其分布。在此基础上,用总体估计值一致性(OEC)、仿真变动系数均值(ACV)和相对均方根误差(RRMSE)指标分析仿真精度; 用设置于清查样地周围的临时样地(验证样地)数据与LandsatTM数据进行森林碳序列高斯块协同仿真,分析清查样地的空间代表性和森林碳分布空间仿真的尺度上推方法。结果表明: 仙居县2008年森林总碳储量仿真估计值为2 667 878 Mg,大部分分布在南部和北部山区,中部东西向条带状低海拔区域分布较少; 区域碳密度仿真估计值为0~65.66 Mg·hm-2,无论是全部样地还是减少一半样地,仿真结果总体均值均在抽样估计置信区间以内; 基于清查样地与基于加密的验证样地森林碳仿真结果表明30 m × 30 m水平样地位置碳密度相关系数达0.95,以清查样地为中心1 km × 1 km块的碳密度相关系数为0.85,说明1 km × 1 km样地仍具有较好的代表性,块仿真效果满意; 以(1-RRMSE)/n定义成本效益,则使用一半样地得到的成本效益优于使用全部样地的结果,用此指标有望找到满足给定精度的最经济的样地数量。
关键词森林碳制图    空间协同仿真    Landsat TM    精度评估    尺度上推    
Estimation of Forest Carbon Distribution for Xianju County Based on Spatial Simulation
Zhang Maozhen1, 2, Wang Guangxing3, Ge Hongli1, 2, Xu Lihua1, 2    
1. School of Environmental & Resource Sciences, Zhejiang A&F University Lin'an 311300;
2. Zhejiang Provincial Key Laboratory of Forest Ecosystem Carbon Cycling and Sequestration, and Emission Reduction Zhejiang A & F University Lin'an 311300;
3. Department of Geography and Environmental Resources, Southern Illinois University Carbondale Carbondale 62901
Abstract: In this study, forest carbon stock and its spatial distribution of Xianju County, Zhejiang, were simulated by using a sequential Gaussian co-simulation algorithm based on data, collected in a forest resource inventory from sample plots in 2008, and the Landsat Thematic MapperTM image. The obtained forest carbon map was assessed using three accuracy measures of overall estimate consistency (OEC), average coefficient of variation (ACV), and relative root mean square error (RRMSE). Moreover, temporary sample plots were selected surrounding the permanent sample plots and the data were collected. The permanent and temporary plot data sets were respectively combined with the TM image to scale up forest carbon stock from a spatial resolution of 30 m × 30 m to 1 km × 1 km using a sequential Gaussian block co-simulation algorithm. The up-scaled map from the temporary plot data were used to assess the accuracy of the corresponding map from the permanent plot data and to analyze their spatial representatives. The results showed that the 2008's forest carbon stock estimated for the county was 2 667 878 Mg. Most of the stock was found in the southern and northern mountainous areas and less amount of the stock in the central west to the east narrow areas that had lower elevation. The forest carbon density varied from 0 to 65.7 Mg·hm-2. Whether all the permanent plots or half of them were used, the population estimates of forest carbon fell into the confidence intervals of the plot data at a significant level of 5%. When the data from the permanent sample plots were compared with those from the temporary plots, the coefficient of correlation for carbon density was 0.95 and 0.85 at the spatial resolutions of 30 m × 30 m and 1 km × 1 km, respectively. This result implied that when forest carbon stock was scaled up from a spatial resolution of 30 m × 30 m to 1 km × 1 km, the obtained map using the permanent plots re-produced the spatial distribution of the forest carbon density very well. In addition, if a cost efficiency was defined as (1-RRMSE)/number of sample plots, using half of the sample plots showed a higher cost efficiency than using all the sample plots, implying this cost efficiency indicator can be used to determine an appropriate number of field plots in sampling design for generating spatial distribution of forest carbon stock.
Key words: forest carbon mapping    spatial co-simulation    Landsat TM    accuracy assessment    scaling up    

森林是陆地生态系统的主要植被类型,碳储量占陆地生态系统的2/3,在全球碳循环中具有重要作用和地位(Dixon et al.,1994;Feng et al.,1999;Bolin et al.,2000)。准确估算区域森林碳储量及其分布是区域森林CO2源汇功能评价的基础,也是解释全球碳收支计算仍存在不平衡问题的一个关键因素(Fang,2000a)。但是,由于方法和数据缺乏,不同的学者估测的森林固碳量相差较大,导致森林碳汇功能评价具有较大的不确定性(Tans et al.,1990)。因此,研究如何有效估算区域森林碳储量及其分布,对于区域碳源汇评估、碳汇交易以及我国森林生态系统在全球碳平衡中的作用和地位评价具有重要意义。

当前,遥感技术在森林碳储量及其分布信息获取中发挥越来越重要的作用。利用遥感数据和地面样地数据进行森林碳储量反演,可在各种尺度上获取区域森林碳储量及其分布信息。在过去的20多年里,不少学者基于各种遥感影像数据和包括森林资源清查数据在内的地面调查数据,研究区域森林碳储量及其分布的估算方法,通过遥感反演获得森林碳储量的分布。

目前比较常用的基于遥感手段获取森林碳储量及其分布信息的方法有3类:基于遥感影像特征与地面森林碳密度的回归估计(RA)(Feng et al.,1999;Tans et al.,1990)、非参数神经网络(ANN)模拟(Fang,2000b;琚存勇等,2006;张超等,2012)和空间仿真(SS)模拟(冯益明等,2004;Wang et al.,2009)。

RA是广泛运用的预测方法。以遥感影像特征值、地理信息等作为自变量,森林碳密度调查值作为因变量建立回归模型和估计。尽管该方法建模和用模型估计总体总量一般能达到较好的精度,但所得区域森林碳储量的分布格局与现实相差较大,在一定程度上平滑了森林碳储量的空间分布(沈希等,2011)。回归模型的另一个主要缺陷是其估计值受地面样地数据的限制,导致其预测能力比较差,甚至产生森林碳密度负值(张茂震等,2009a)。

ANN是一种模仿生物神经网络的结构和功能的非参数估计方法,与传统的统计模型相比具有速度快、精度较高的优越性,而且不需要有关模拟过程和目标方程结构假定的预备知识。人工神经网络不会产生负的估计值。但是,由于人工神经网络是高度非线性的大型系统,其复杂性之高使得其难以分析各项性能指标(杨晓帆等,1994),无理论指导隐含层数和单元数的确定,学习与记忆具有不确定性(宰松梅等,2011),训练时间长,结果难以重复。回归估计和神经网络法估计结果的精度评估主要以均方根误差(RMSE)和残差(RE)为指标,基于地面样地位置进行(范文义等,2011)。当地面样地数量有限时,这种方法甚至连区域内各样地位置的估计精度都不能如实反映。对于人工神经网络方法来说,虽然预测精度往往多数都能达到90%以上,甚至99%,但预测值与实际值的相关系数却远低于这个值(琚存勇等,2006;张超等,2012)。

SS是对地理空间上真实事物或过程的模拟。它采用与地统计学方法相结合的随机算法模拟森林碳分布,通过对局部森林碳分布特征量的分析,得到局部森林碳分布函数,再用蒙特卡洛方法实现对局部的估计。它不像回归方法那样只考虑保证总体平均数的估计精度,而是最大限度地再现森林分布。该方法不仅能提供森林碳储量的空间分布和保证局部估计精度,而且能提供其估计值的方差(张茂震等,2009b)。

基于仿真的方法,尽管在反映局部特征、再现森林碳分布方面有其独特优势(Wang et al.,2009;沈希等,2011),但其总体估计精度、地面样本空间变异特性还少有研究。除此之外,国家和大区域的森林碳分布图经常要求低的空间分辨率,即大的像元,如1 km×1 km,而森林清查样地一般小于30 m×30 m,这就要求进行森林碳估计的尺度上推。由于缺乏大像元森林碳野外测定值,尺度上推后的森林碳估计值精度常常是未知的。本研究以浙江省仙居县为研究区,用2008年森林清查样地数据与同年度L and sat TM遥感数据进行森林地上部分碳储量(AGC,以下简称森林碳)的序列高斯协同仿真,分析评价森林碳储量及其分布空间仿真估计和尺度转换的精度。通过加密的临时样地(验证样地)与遥感数据协同仿真,评价研究区森林碳空间变异特征。

1 研究区概况

仙居县(120°17′16″—120°55′51″E,28°28′14″— 28°59′48″N)地处浙江省东南部,位于浙江第三大水系——椒江水系的源头,东西长63.6 km,南北宽57.3 km。全县总面积2 013.18 km2

仙居县属浙东丘陵地带,山系盘桓,溪流切割,低山和丘陵占全县总面积的86.4%。地表分割强烈,河谷切割深邃,海拔742~1 384.4 m,原始地表已被完全破坏。该县属典型亚热带季风气候,温暖湿润,四季分明,年均气温17.2 ℃,年均降雨量1 377 mm,年均蒸发量1 260.8 mm,年均日照时数1 932.6 h。

全县林业用地面积16.46万 hm2。林业用地中,有林地15.336 9 万hm2(其中竹林面积占6.73%),灌木林地62.89万hm2,未成林造林地0.257 9万hm2,无立木林地0.188 9万 hm2。森林覆盖率77.9%,森林总蓄积量555.5万 m3 (2008年)。

由于历史原因,境内大部分地区原生植被已经被反复利用和改造,代之以次生植被为主。境内植被分为暖性针叶林、常绿阔叶林、落叶阔叶林、常绿落叶阔叶混交林、针阔混交林、竹林和山地灌草丛7个类型。

2 研究方法 2.1 森林样地清查数据

地面调查数据为仙居县2008年森林资源调查结果,即清查样地数据。全县共有按系统抽样方法布设的固定样地327块,本研究中实际有效样地302块,样地形状为正方形,样地面积0.08 hm2,样地最大碳储量6.953 Mg。在模拟、分析过程中,样地碳储量全部采用碳密度表示,即每公顷碳储量(Mg·hm-2),如图 1

图 1 森林资源清查样地数据和Landsat TM影像 Fig. 1 Locations of inventory plots data and the Landsat TM image

清查样地数据随机分成2个数据子集:样地子集A和样地子集B,数据量分别为90%和10%。其中,样地子集A(90%)用于仿真,样地子集B(10%)用于检验。为了与2008年实地验证样地数据相区分,用于检验的样地数据称之为检验数据集。

2.2 遥感数据

研究区的遥感数据为2008年3月获取的L and sat TM数据。

该L and sat TM影像分幅号为119-40,仙居县辖区范围完全位于一景影像之内。L and sat TM影像采用北京54坐标系,6度分带,高斯克里格投影进行几何校正。仙居县经度范围为120°—121°E,6度带投影无跨带处理。影像在几何校正的基础上采用DOS法进行辐射校正处理。几何校正后影像几何误差小于1个像元,如图 2

图 2 森林碳样本空间变异标准化模型 Fig. 2 Standardized sample variogram of forest carbon

在此基础上,提取样地位置所对应的各波段的DN值,进行不同组合比值运算,计算TM1/TM4,TM2/TM4,TM3/TM4,TM4/TM3,TM5/TM4,TM7/TM4,(TM4-TM3)/(TM4+TM3),(TM3+TM5)/TM7,(TM3+TM5+TM7)/TM4和(TM2+TM3+TM5)/TM7共10种波段比值或波段组合比值,最后计算其比值与样地森林碳储量之间的相关系数,选取相关系数最大的波段比值图(TM4/TM3)参与仿真模拟。

2.3 验证样地设置与调查

2010年8月,在清查样地周围4个方向(东、南、西、北)或8个方向(北、东北、南、东南、南、西南、西、西北)建立和调查用于分析清查样地估计效果的验证样地。调查样地地类、林种、起源、郁闭度、优势树种、主要灌木种类、草本盖度、灾害影响、平均树高和平均胸径;调查样木立木类型、检尺类型、树种、胸径、树高和健康状况。设4个验证样地时,4个方向上验证样地与清查样地之间的距离均为250 m。设8个验证样地时,8个方向上验证样地距清查样地的距离分别为250,353,250,353,250,353,250和353 m(东、西、南、北4个方向上距清查样地250 m,东南、东北、西南、西北4个方向上距清查样地353 m)。验证样地用于分析以清查样地为中心的1 km2范围内森林碳密度空间变异特征和清查样地的代表性。

2.4 样地数据预处理

本研究中,清查样地和验证样地的调查时间分别为2008和2010年。为了2期数据具有可比性,以2008年清查样地调查时间为基准,修正2010年验证样地数据,即将验证样地的森林碳储量推算到2008年。2010年验证样地调查时基本未发现验证样地内近2年有采伐,因此根据各树种生长率和2010年验证样地调查数据推算验证样地2008年森林碳储量。其生长率取自《第八次全国森林资源清查浙江省森林资源清查成果(2009年)》(国家林业局华东森林资源监测中心、浙江省林业厅,2010.3)。

2.5 空间仿真

为获取研究区森林碳分布,采用基于地统计理论的仿真方法进行森林碳储量空间仿真。这种仿真是一种条件模拟方法,它把研究区划分成N个等大小的单元,采用序列高斯仿真算法导出每个单元的估计值。该算法假设每个单元的估计值是一个随机函数在该位置的随机变量Z(u)的实现,其概率分布假定为正态分布,并由该点周围的样地数据确定。当被估计位置的遥感影像数据及其周围已有的估计值也被用来确定这个条件分布时,此法称为序列高斯协同仿真(sequential Gaussian co-simulation,SGCS)。这个条件分布由一个统计平均数和方差来确定,而这个统计平均数和方差可以基于相邻位置的地面数据、已有的模拟数据和遥感影像数据用同位协同简单克里格估计来获得。

$\begin{align} & {{z}^{\text{sck}}}\left(u \right)=\sum\limits_{a=1}^{n\left(u \right)}{\lambda _{a}^{\text{sck}}\left(u \right)\left[ z\left({{u}_{a}} \right)-{{m}_{z}} \right]}+\\ & \lambda _{y}^{\text{sck}}\left(u \right)\left[ y\left(u \right)-{{m}_{y}} \right]+{{m}_{z}};\\ \end{align}$ (1)
$\begin{align} & {{\sigma }^{2\left(\text{sck} \right)}}\left(u \right)={{c}_{zz}}\left(0 \right)-\sum\limits_{a=1}^{n\left(u \right)}{\lambda _{a}^{\text{sck}}\left(u \right){{c}_{zz}}}\left({{u}_{a}}-u \right)- \\ & \lambda _{y}^{\text{sck}}\left(u \right){{c}_{zy}}\left(0 \right).\\ \end{align}$ (2)

式中: u为一个被估计的位置;zsck(u)为位置u的森林碳估计值; z(uα)为位置u局部范围内第α个样地的森林碳实测值;α=1,2,…, n(u)n(u)为在给定的搜索范围内所获得的样地数;y(u)为在像元u处的遥感影像数据;mz,my分别为地面样地数据和遥感影像数据的均值;λαsck(u)和λysck(u)分别为样地数据和影像数据的权重;Czz(0)为森林碳地面样地数据的方差;Czy(0)为森林碳和遥感影像数据的协方差。当h=uα-u时,Czy(h)为森林碳与遥感影像数据的交叉协方差函数。

当输出的森林碳分布图的像元与森林清查样地大小一致时,SGCS过程为: 1)用随机抽样方法设置一个估计每个像元的顺序;2)在每个像元的位置u,用同位协同简单克里格方法计算条件累积分布函数的平均数和方差;3)从此条件累积分布中进行随机抽样获得一个值,这个值作为随机变量在此像元的位置u的实现。重复步骤1)到3),直至所有像元都有模拟数据,由此获得一张整个研究区森林碳的分布图。每次用随机抽样方法设置估计每个像元的顺序,将这一过程执行L次,就可以得到L张分布图,最后以像元为单位计算森林碳分布图的平均数图、方差图和大于某一特定值的概率图。由于仿真是随机模拟,因此结果的精度在很大程度上取决于仿真次数。森林碳储量仿真试验表明,仿真模拟达到200次以后,各像元位置碳储量估计特征值趋于稳定。因此,本研究中仿真次数L全部取200。

当输出的森林碳分布图的像元(如1 km×1 km)比森林清查样地(如30 m×30 m)大时,其仿真过程需要用序列高斯协同块仿真SGBCS来完成尺度上推。这种SGBCS与当输出的森林碳像元和森林清查样地大小一致时的SGCS相似。其不同的是,SGBCS将每个大的像元看成由许多小像元所组成,每个小像元用同位协同简单克里格方法计算一个估计值,然后计算大像元内所有小像元估计值的平均值和方差,由此确定大像元森林碳累积概率分布函数,并完成空间模拟。

以上样地数据的权重λαsck(u)与样地至估计位置的距离和森林碳储量的变异函数相关。变异函数是随机函数空间变异或空间自相关的度量。相交的变异函数可以度量2个随机函数相互的空间相关关系。设变量z为森林碳储量,则z(u)为定义在二维空间u处的随机函数。随机函数的变异函数${{{\hat{\gamma }}}_{z}}\left(h \right)$可用下式来计算:

$\begin{align} & {{{\hat{\gamma }}}_{z}}\left(h \right)=\frac{1}{2N\left(h \right)}\sum\limits_{a=1}^{n\left(h \right)}{\left[ z\left({{u}_{a}} \right)\right.} \\ & {{\left.z\left({{u}_{a}}+h \right)\right]}^{2}}.\\ \end{align}$ (3)

式中: N(h)为在给定距离和方向上样本数据对的数量;z(uα)z(uα+h)分别为uαuα+h空间位置上的森林碳储量;h为矢量或给定方向上的距离差。通常样本半方差用下式模拟:

${{\gamma }^{\text{s}ph}}\left(h \right)=\left\{ _{{{c}_{0}}+{{c}_{1}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 其他.}^{{{c}_{0}}+{{c}_{1}}\left[ 1.5\frac{h}{a}-0.5{{\left(\frac{h}{a} \right)}^{3}} \right]\left(0<h\le a \right)} \right.$ (4)

式中: c0为块金值;c1为结构参数;c=c0+c1为基台值;a为变程——当模型值等于最大基台值时的距离。结构和基台值说明了空间变异的结构和总方差。块金值是在h=0时函数值的跳跃程度,可能代表野外调查或遥感测量中的微观变异、误差等因素。变程参数提供了随机变量空间相关的最大范围。如果森林碳储量有一个常数均值和方差,则它的协方差函数Cz(h)存在,且与变异函数有如下关系:

${{C}_{z}}\left(h \right)={{C}_{z}}\left(0 \right)-{{\gamma }_{z}}\left(h \right).$ (5)

受地面样地大小的限制,一般SGCS的尺度为30 m×30 m,在此基础上生成较大区域的森林碳分布图,需要进行尺度转换。本研究采用SGBCS方法,直接用30 m×30 m分辨率的数据源进行块仿真,得到分辨率为1 km×1 km的森林碳分布图。

本研究进行4种仿真: 1)基于清查样地的SGCS(30 m×30 m)(以下简称仿真);2)基于验证样地的SGCS(30 m×30 m)(以下简称验证仿真);3)基于清查样地的SGBCS(1 km×1 km)(以下简称块仿真);4)基于验证样地的SGBCS(1 km×1 km)(以下简称验证块仿真)。空间仿真程序基于ClaytonV.Deutsch和Andre G.Journel的GSLIB(geostatistical software library and user’s guide)(Deutsch et al.,1998)实现。空间变异半方差函数参数拟合工具为VARIOWIN2.2。

2.6 精度评估

用L and sat TM遥感影像和清查样地进行空间仿真和尺度转换,得到分辨率为30 m×30 m和1 km×1 km的森林碳估计值平均数、方差和分布概率图。基于均值、方差等统计特征数,可以对森林碳空间仿真进行精度评估。本研究采用3个指标从不同侧面描述估计精度:仿真结果变动系数均值ACV、相对均方根误差RRMSE和置信区间CI。ACV反映总体抽样变动情况,RRSME反映样地位置估计值的相对误差大小,CI可衡量区域分量之和与总量估计值的符合程度。

2.6.1 仿真结果变动系数均值(ACV)

在森林碳仿真模拟中,一个像元的估计值是一个随机函数在该位置的随机变量Z(u)的实现,该随机变量有均值和方差。随机变量通过从某个已知概率分布中随机抽取有限次来实现,得到代表该像元状态的均值与反映该估计值离散程度的方差。方差可以用作该位置估计值离散或偏差程度的度量。但由于每个像元位置的碳密度不同,方差的绝对值不具可比性,因此,这可用标准差与均值之比。将方差图开方并与均值图进行比值运算,可得每个像元位置的仿真结果变动系数cvi

$c{{v}_{i}}=\sqrt{{{\operatorname{var}}_{i}}}/{{\text{m}}_{i}}.$ (6)

式中: vari为第i像元位置碳储量估计值L次抽样的方差;mi为第i个像元位置碳储量估计值L次抽样的平均值;cvi代表第i个像元位置碳储量估计值的变动系数。

取cvj的平均值代表总体的仿真结果平均变动系数ACV。

$ACV=\frac{1}{n}\sum\limits_{i=1}^{n}{c{{v}_{i}}}.$ (7)

式中: n为研究区像元个数。

2.6.2 相对均方根误差(RRMSE)

从平均值图中提取样地位置的模拟值sci,用sci与样地实测值gci比较,可得均方根误差(RMSE)

$RMSE=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{{{n}_{\text{p}}}}{{{\left(s{{c}_{\text{i}}}-g{{c}_{\text{i}}} \right)}^{2}}.}}$ (8)

式中: np为样地数。

虽然一般都以RMSE描述模拟结果的准确程度,但它是一个绝对数,大小与均值相关,无法直观反映结果的精确程度。本研究采用相对均方根误差(RRMSE)作为整个区域估计精度的主要评估指标之一,它表达了标准差相对平均值的变动程度。

$\text{RRMSE=RMSE}/\frac{1}{n}\sum\limits_{i=1}^{{{n}_{\text{p}}}}{g{{c}_{\text{i}}}.}$ (9)

式中:$\frac{1}{n}\sum\limits_{i=1}^{{{n}_{\text{p}}}}{g{{c}_{\text{i}}}}$为n个样地位置实测值的均值。

2.6.3 总体估计值一致性(OEC)

仅用地面样地数据虽然不能得到区域森林碳分布,但可以估计总体总量的估计均值和置信区间。通过仿真模拟得到区域森林碳分布(均值图),累加均值图各像元位置的估计值可以得到森林碳总体总量的仿真估计值。如果这个估计值落入由地面样地估计总体总量的置信区间(confidence interval),则仿真模拟结果在总量上可信。当地面样地数量满足大样本要求时,可以用总体总量的仿真估计值与样地估计值的符合程度来衡量森林碳分布估计的精度。当用仿真估计得到的总体总量落入基于样地估计的总体总量置信区间时,总体估计值一致性(OEC)可用仿真估计值与样地估计值置信区间(CI)的中值来度量。令森林碳总体总量为TC,样地估计结果为TCp(即置信区间的中值),仿真估计结果为TCs,则有:

$\text{T}{{\text{C}}_{\text{p}}}=\frac{A}{a}\frac{1}{m}\sum\limits_{j=1}^{m}{g{{c}_{j}};}$ (10)
$\text{T}{{\text{C}}_{\text{S}}}=\sum\limits_{i=1}^{n}{s{{c}_{i}}.}$ (11)

式中: A为总体面积;a为样地面积;m为样地数。

由此有:

$\text{OEC}=\text{T}{{\text{C}}_{\text{S}}}/\text{T}{{\text{C}}_{\text{p}}}.$ (12)

OEC代表仿真估计值与地面抽样总体估计值的接近程度,二者相等时,OEC=1。

3 结果与分析 3.1 空间变异函数

根据研究区2008年森林清查样地数据分析,302个样地森林碳密度的标准差为15.25,空间变异模型如图 2

其半方差函数为:

$\left\{ \begin{align} & \text{0}\text{.65+0}\text{.35}\left[ 1.5\frac{h}{4860}-0.5{{\left(\frac{h}{4860} \right)}^{3}} \right] \\ & 1\left(h\text{a} \right).\\ \end{align} \right.\left(0 <h\le a \right)$ (13)

式中:h为距离,rSph(h)为标准半方差,块金值c0=0.65,结构参数c10.35,变程a=4 860 m。

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

3.2 空间仿真

用302个森林清查样地数据与2008年L and sat TM遥感影像数据仿真,得到研究区30 m×30 m空间分辨率的森林碳分布图(图 3)。

图 3 森林碳空间仿真结果 Fig. 3 Results of forest carbon obtained using the conditional co-simulation

估计值的均值分布图简称均值图。均值图(图 3a)提供了每个像元碳密度仿真值的平均值。由均值图得全县森林碳密度的平均值为13.709 Mg·hm-2,标准差为8.327 Mg·hm-2,像元的碳密度为0~65.6587 Mg·hm-2,全县森林碳储量为2 667 878 Mg。

估计值的方差分布图简称方差图。方差图(图 3b)展示了各像元位置碳密度仿真值的方差。由于仿真过程所基于的碳密度分布函数来自像元周围有限个已知像元,因此,每个位置的方差代表了该位置碳密度的不确定性的大小。总体上看,南部和北部森林较多,碳储量变化也大,而中部变化小。方差最大值为575.075,将其转换成标准差与均值图的最大值比较,其比值为0.365 2,该数据也在一定程度上反映仿真结果中样本不确定性的大小。

各个像元位置估计值大于某一阈值的概率分布图简称概率图。概率图(图 3c)是每个像元位置碳密度大于平均值的概率分布图。图中每个位置的估计值是该位置200次仿真中大于均值的次数与总次数之比。

标准差(方差开平方)图与均值图的比值所形成的图即变动系数图(图 3d)。它反映了各空间位置从概率分布函数中多次抽样所形成样本空间的特征值的变异情况,它区别于地表两点之间特征(碳密度)的差异程度。

以上平均值图、方差图和概率图显示,森林碳估计值存在显著的空间自相关和空间分布的一致性,即估计值大的地方,其方差和大于平均值的概率也较大,反之亦然。然而,变动系数图的趋势相反,即估计值大的地方,其变动系数小,说明森林碳密度高的区域仿真结果稳定,而人口密集的城区、村镇区仿真估计值的相对误差大。

基于森林资源清查数据估算,浙江省2004年平均单位面积碳储量为11.02 Mg·hm-2(张茂震等,2009a),仙居县空间仿真平均森林碳密度高于全省平均数24.4%,据研究,空间仿真估计总体总量低于抽样估计总量(张茂震等,2009b),仙居县森林实际碳密度可能高于全省平均数50%左右。

3.3 精度分析 3.3.1 RRMSE

RMSE分析用前面提到的数据集B进行,检验10%样地位置的仿真估计值与其对应样地实测值的吻合程度。结果如图 4。据(8)式和(9)式分别计算均方根误差RMSE为8.617,相对均方根误差RRMSE为0.776。

图 4 仿真结果验证 Fig. 4 Comparison of simulated values with the field measurements
3.3.2 OEC和ACV

按照(6)式计算每个像元位置的抽样变动系数,然后用(7)式计算得到全部变动系数的平均值为1.172。这里像元数n为2 230 642,样地数为302,模拟次L数取200。

将模拟结果的像元值与样地调查数据比较,OEC为0.910,说明仿真结果的总体估计值与抽样的值基本接近。ACV为1.17,表明森林碳密度在研究区总体上空间变化仍然比较强烈,特别是在中部县城及附近盆地、河谷地区碳密度均值很小,但方差相对均值可能很大,导致变动系数(标准差比均值)较大。

3.3.3 样地数与成本效益

现从302个样地中随机抽取1/2样地组成样本,采用相同方法分别与遥感数据进行仿真模拟,与用全部样地仿真模拟结果对比,结果如表 1

表 1 不同样地数仿真模拟结果的比较 Tab.1 Comparison of results using different numbers of sample plots

无论是用全部样地(302个)或一半样地(153个),仿真模拟结果的全域森林碳总体总量均落入用样地数据估计的置信区间,而且比较稳定,仿真模拟结果与对应的置信区间中值(表 1)分别差9.02%和9.68%,可能是研究区空间一致性较好或者空间仿真过程的搜索半径增大的结果。

减少地面样地数量,估计精度随之下降,但满足一定估计精度的成本效益呈上升趋势。当使用302个样地时,以(1-RRMSE)/n定义的成本效益为0.000 695,而使用153个样地得到的成本效益为0.000 719,后者高于前者。即使用153个样地得到的成本效益较使用302个样地为优,但究竟使用多少样地成本效益最优,目前仍不清楚。

3.4 局部变异验证

在30 m×30 m空间分辨率水平上,分别用清查样地和验证样地数据与L and sat TM遥感数据进行高斯协同仿真,得到基于清查样地的森林碳分布均值图和基于验证样地的森林碳分布均值图,并将尺度上推到1 km×1 km水平(图 5)。

图 5 基于不同样地的块仿真均值 Fig. 5 Maps of up-scaled forest carbon using block co-simulation based on permanent and temporal sample plot data

在30 m×30 m空间分辨率水平上,比较基于清查样地的仿真均值图和基于验证样地的仿真均值图,结果具有高度一致性。在2个结果中分别按验证样地位置提取森林碳密度,其相关系数为0.94,分别按验证样地所涉及的清查样地位置提取森林碳密度,其相关系数为0.95,表明清查样地与周围的验证样地在森林碳密度上具有很好的一致性。

在1 km×1 km空间分辨率水平上比较显示,基于清查样地的森林碳分布均值和基于验证样地的森林碳分布均值的相关系数为0.85,清查样地用于森林碳分布的协同仿真估计具有良好的代表性。

最后比较基于块协同仿真法和基于平均数法的统计特征数,结合样地特征分析空间仿真结果精度和清查样地的代表性。表 2列出了22组验证样地与对应的清查样地所代表的各个1 km×1 km范围内森林碳密度仿真的方差。

表 2 基于验证样地的块内方差分析 Tab.2 Analysis of variance within the simulated block based on test plots

表 2以方差显示了2种估计结果的特征差异。其中块协同仿真结果的方差均值大于平均数法,反映块协同仿真法能更好地体现局部特征,而且最大值和最小值也体现了这一特点。平均数法的方差范围明显收窄,平均值降低。

4 结论与讨论

结合森林清查样地数据和中高空间分辨率遥感数据,用空间仿真方法可以有效估计县域范围森林碳储量及其空间分布,基于仿真方差可进行森林碳制图的精度评估。

仙居县森林总碳储量仿真值为2 667 878 Mg。空间仿真结果在总量上能较好地落入抽样估计的置信区间,在分布上与抽样样地对照相关系数达到0.95以上。

在一定范围内,减少样地数量对仿真估计精度的影响小于其对抽样估计精度的影响。

基于清查样地与基于验证样地的尺度上推仿真结果虽然在分布趋势上一致,但在碳密度值域分布上有一定差别,前者最大值明显较后者大,其原因是2010年调查的验证样地数据总体上较2008年的清查样地数据值低。造成验证样地数据在总体上低于清查样地的原因主要是2008年到2010年之间的森林采伐和林种改造,2009年浙江省森林资源连续清查报告显示,2005—2009年全省年均生长量1.977×107 m3,年均消耗量1.170×107 m3,消耗量中超过78%为采伐量。当采伐在一定程度上偏离清查样地时,就有可能造成清查样地数据偏大。自然灾害引起地类破碎化也是其中原因之一,但从布设的验证样地情况看,其影响比较小。

仿真结果与实际碳密度的吻合程度仍有提高的空间。理论上搜索半径越小,其仿真结果越接近实际,但实际上小到一定程度可能搜索不到样地,以至于无法估计或估计误差很大。而搜索半径太大又会平滑掉区域内的信息,最终结果与回归方法类似,达不到仿真的目的。因此,搜索半径的优化应该是下一步研究的重点。

以(1-RRMSE)/n定义的成本效益(cost efficiency)指标在本研究中与样地数之间显示了一定的规律,对于地面调查抽样设计有重要意义,但结论仍有待于更大范围的研究。

参考文献(References)
[1] 范文义, 张海玉, 于 颖, 等. 2011.三种森林生物量估测模型的比较分析. 植物生态学报, 35 (4): 402-410.(1)
[2] 冯益明,唐守正,李增援. 2004. 应用序列指示条件模拟算法模拟森林类型空间分布. 生态学报, 24(5):946-952.(1)
[3] 琚存勇, 蔡体久. 2006. 用泛化改进的BP 神经网络估测森林蓄积量. 林业科学, 42(12):59-62.(2)
[4] 沈 希,张茂震,祁祥斌. 2011. 基于回归与随机模拟的区域森林碳密度与储量估计方法比较.林业科学, 47(6):1-8.(2)
[5] 杨晓帆,陈廷槐. 1994. 人工神经网络固有的优点与缺点. 计算机科学,21(2): 23-26.(1)
[6] 宰松梅,郭冬冬,韩启彪,等. 2011. 基于人工神经网络理论的土壤水分预测研究.中国农学通报, 27(8):280-283(1)
[7] 张 超,彭道黎. 2012. 基于PCA-RBF神经网络的森林碳储量遥感反演模型研究. 中国农业大学学报, 17(4):148-153.(2)
[8] 张茂震,王广兴,刘安兴. 2009a. 基于森林资源连续清查资料估算的浙江省森林生物量及生产力.林业科学,45(9): 13-17.(2)
[9] 张茂震, 王广兴, 周国模,等. 2009b. 基于森林资源清查、卫星影像数据与随机协同模拟尺度转换方法的森林碳制图. 生态学报, 29(6): 2919-2928.(2)
[10] Bolin B, Sukumar R, Ciais P, et al. 2000. Special report on land use, land-use change and forestry. New York: Cambridge University Press, 30-51.(1)
[11] Deutsch C V, Journel A G. 1998. GSLIB: geostatistical software library and user's guide. New York: Oxford University Press.(1)
[12] Dixon R K, Solomon A M, Brown S. 1994. Carbon pools and flux of global forest ecosystems. Science, 263( 5144): 185-190.(1)
[13] Fang J Y. 2000a. Forest productivity in China and its response to global climate change. Acta Phytoecologica Sinic, 24( 5): 513- 517.(1)
[14] Fang J. 2000b. Forest biomass carbon pool of middle and high latitudes in the north hemisphere is probably much smaller than present estimates.Acta Phytoecologica Sinica, 24(5): 635-638.(1)
[15] Feng Z W, Wang X K, Wu G. 1999. The productivity and biomass of forest ecosystems in China. Beijing: Science Press.(2)
[16] Tans P P, Fung I, Takahasi T. 1990. Observational constraints on the global atmospheric carbon budget. Science, 247:1431-1438.(2)
[17] Wang G,Oyana T,Zhang M, et al. 2009. Mapping and spatial uncertainty analysis of forest carbon by combining national forest inventory data and satellite images. Forest Ecology and Management, 258(7):1275-1283.(2)