区域尺度森林地上生物量估计与不确定性度量

傅 煜

(上海市园林科学规划研究院；上海城市困难立地绿化工程技术研究中心 上海 200232)

1 导言

作为《京都议定书》和《巴黎协定》的签署国之一，中国有义务向气候变化框架公约(UNFCCC)秘书处定期公布本国森林生物量、碳储量和碳汇能力估计，以及相应的不确定性度量方法、结果及降低措施[1]。而地上生物量(AGB)是五大森林碳库之一，是估算碳储量的计算基础[2]。

区域尺度生物量估计通常基于固定样地连续观测数据，通过异速生长方程预测样地内单木生物量，随后求和并推算至更大尺度的生物量[3]。该过程中存在多种误差来源，并造成了生物量估计的不确定性[4，5]。在实际生产中，这些不确定性被普遍忽略或笼统归类为一项随机误差[6]，而忽略这些不确定性会导致估计精度的降低[7]。因此，根据不同的误差来源，准确地分解和量化不确定性对提高精度和准确性，为科学、可靠地预测生物量和碳储量具有重要意义[8]。

估计不确定性的主要方法有三种：①误差传播法，通过使用泰勒级数一阶项作为不同来源误差的近似值，并将各误差源导致的不确定性加以求和来计算[9]。但各误差源间存在相关性，机械地相加的方式缺乏严谨性[10]。②模型分析法，基于系统抽样理论和泰勒级数原理，可分开量化抽样误差和模型误差引起的不确定性[11，12]。然而，由于估算模型误差导致的不确定性基于模型参数及其协方差矩阵，即该方法评估的不确定性较为依赖建模数据的质量，估计结果的稳健程度随着建模数据集的变化而降低。③蒙特卡洛法，通过反复模拟生物量估计过程，以获得模型参数、生物量估计和不确定性的概率分布[3]，可对生物量及其不确定性进行稳健估计[13]，但该方法较少用于分解不同误差源引起的不确定性[14]。

由于模型分析法可分解不同来源误差导致的不确定性但缺乏稳定性，而蒙特卡洛法稳定性高但缺乏分解误差源的相关研究，本文的研究目标是将上述两种方法结合起来，其基本原理是将模型分析法嵌套进蒙特卡洛法的每次循环中，以实现两种方法的优势最大化发挥。本研究的目的是：①采用通用的生物量模型估算区域尺度的AGB；②采用本文方法分解抽样误差和模型误差引起的不确定性。

2 材料与方法

2.1 数据来源

本研究调查数据基于江西省(118°29′N, 113°34′E)固定样地连续观测数据，以江西省主要树种杉木(Cunninghamialanceolata)为研究对象，数据收集于2009～2011年。研究区包括1 128个包含杉木的样地，样地面积为0.067 hm2。样地内每木测量胸径(起测直径为5 cm)。由于树高测量较为困难且精度低，本研究采用分级树高-胸径曲线模型[15]估计研究区内所有树高。调查数据统计特征如表1所示。

表1 研究区调查数据统计特征

建模数据采用江西、福建省150株杉木立木的地上生物量实测数据，采集时间为2009年6～9月。 样木数按5、7、11、15、19、25、31、37 cm以上8个径阶均匀分配，各径阶样木数按树高等级由低到高尽量均匀分布，并在大范围内具有广泛的代表性。全部样木现地实测胸径、地径和冠幅，将样木伐倒后，测量树高和冠长，分干、皮、枝、叶称鲜重。分别抽取样品带回实验室，在85 ℃下烘干至恒重，然后根据样品干重推算样木各器官干重并汇总得到地上部分干重。建模数据统计特征如表2所示。

表2 研究区建模数据统计特征

2.2 研究方法

2.2.1 生物量模型与残差

为选择拟和效果最佳的估计模型，本文尝试了三种通用的异速生长方程[16]，分别为β0dbhβ1、β0dbhβ1hβ2和β0·(dbh2·h)β1，最终选择了决定系数最高的公式(1):

gi=β0dbhiβ1hiβ2+εi

(1)

其中，i为建模数据中单木号，单木胸径dbh、树高h为自变量，地上生物量g为因变量。模型拟合采用非线性最小二乘法，β0，β1，β2为模型参数，εi为模型参差：

εi～N(0，exp(γ1+γ2·ln(dbh)))

(3)

2.2.2 Monte Carlo模拟步骤

本研究方法原理是将模型分析方法与蒙特卡洛模拟相结合，以获得这两种方法的优势，同时避免它们的弱势。该方法的核心是在Monte Carlo模拟的每次循环过程中实现模型分析方法，即能稳健估计生物量，又能有效分解抽样误差和模型误差的不确定性。

步骤 1:

步骤 2:

基于新的生物量“观测值” 和胸径、树高，重新拟合模型参数，公式如下：

g＇i=α0·dbhiα1·hiα2+εi

(4)

式中α0,α1,α2为新的拟合参数，拟合方法采用最小二乘法。

步骤 3:

采用公式(5)计算调查数据中所有单木地上生物量，公式如下：

(5)

对单木AGB求和获得样地水平Gj，公式如下：

(6)

式中ij表示第j号样地的第i棵树，且有j=1,…,n,i=1,…,nj,nj为第j号样地内某树种的样木棵数；x=(dbh,h)＇为自变量向量，α=(α0,α1,α2)＇为模型参数向量。因此单位面积生物量μk如公式(7)：

式中k为Monte Carlo第k次循环，Aj为第j号样地面积。

步骤 4:

不确定性通常以均方根误差(RMSE)为表达方式[7,10]，因此第k次循环中生物量不确定性为公式(8)：

由于抽样误差和模型误差通常被假设为相互独立[4]，因此，通过方程(9)和(10)可分别度量第k次循环中抽样误差和模型误差(RMSESkand RMSEMk)引起的不确定性[7,10]，公式如下：

步骤 5:

不断重复步骤1～4，直到生物量均值、抽样误差不确定性、模型误差不确定性、总不确定性和相对不确定性估计值趋于稳定，生物量均值和总不确定性计算公式如下：

本文生物量模型拟合、估计及各类不确定性度量采用R软件的实现[17]。

3 结果与分析

3.1 生物量模型与残差

表3为给方程(1)和(3)的拟合结果，方程(1)有较高的决定系数，说明本研究所采用的生物量估计模型具有良好的精度基础。参数γ1、γ2可在接下来的蒙特卡洛模拟程序中被用于生成虚拟残差。

表3 模型参数估计

3.2 生物量估计与不确定性度量

如图1所示，经过Monte Carlo模拟1 000次后，杉木平均AGB及其不确定性的变化趋势逐渐稳定。图中黑色曲线表示平均AGB，灰色指数表示相对不确定性，显然，两条曲线经过400次模拟后越来越光滑。

图1 地上生物量及其不确定性的模拟趋势

表4 地上生物量及其不确定性统计结果

*①采用传统的模型分析法；②采用本文方法

为对比本文方法的效果，本文也实践了Ståhl等[7,10]提出的传统模型分析法，表4为比较结果。如表所示，采用传统的模型分析法和本文方法，杉木生物量分别为17.36±1.09 t·hm-2和17.35±0.98 t·hm-2。两种方法下生物量均值和抽样误差的不确定性差异不大，但本文方法显著降低了模型误差不确定性，进一步降低了相对不确定性。

4 结论与讨论

本文将传统的模型分析方法引入Monte Carlo模拟，对区域尺度地上生物量及其不确定性进行了更加稳健的估计，并分解了抽样误差和模型误差导致的不确定性。该方法结合了Monte Carlo模拟法和模型分析法的双重优点，同时一定程度上降低了模型误差对生物量估计不确定性的影响。

本文方法的优势在于可稳健估计模型参数、生物量及其不确定性。因为采用传统的模型分析法，不确定性的估计精度较为依赖建模数据集的特性。Ståhl等[10]所述的方法中，用于建立生物量模型的数据集是固定的，因此只要拟合方法不变，模型参数和协方差矩阵也是不变的，那么最终得到的生物量和不确定性估计值也是不变的。然而，如果更换为另一组建模样本，估计结果将随之改变。因此，采用传统模型分析法获得的生物量和不确定性估计可能不够可靠和精确。本文将传统模型分析法引入Monte Carlo模拟，估计结果的变异性经反复多次的模拟后，模型参数及其协方差矩阵可变性逐渐降低。随着模拟次数的增加，生物量的估计和不确定性逐渐趋于稳定。

本文方法对区域尺度生物量均值的估算结果与传统模型分析法相差不大，但对生物量估算的不确定性有一定的影响，主要表现在显著降低了模型误差对不确定性的影响，从而降低总不确定性，原因可能是Monte Carlo模拟有效降低了模型参数和协方差矩阵的变异性。本文方法可用于估计区域尺度森林生物量及其不确定性，并分解不同误差来源导致的不确定性，为有针对性提高生物量、碳储量估计精度具有重要意义。