随着人类活动及现代化工业的迅速发展, 特别是化石燃料的燃烧、森林的滥伐、草原开垦等, 对地球生态系统产生着巨大的影响^[1]。作为全球陆地生态系统中的最大有机碳库^[2], 森林在调节全球碳平衡、减缓大气中温室气体浓度上升和维持全球气候稳定等方面具有不可替代的作用。作为反映森林生态系统生产力的重要指标, 森林生物量在计算过程中含有众多不确定性来源, 生物量的估算无论是在区域尺度上还是国家尺度上都存在极大的不确定性^[3]。如何计算这些不确定性以提高森林生物量估算的准确度已成为森林生物量研究的一大挑战。

不确定性是指不精确性、模糊性、不明确性等概念总称^[4]。在森林资源监测体系中, 森林资源连续清查提供了可靠的样地、样木定期观测数据。基于这些数据, 通过单木生物量模型可以获得样地森林生物量, 并可以推算到各尺度。但是从数据采集、模型建立到尺度上推过程中都包含一定的不确定性^[5]。估算结果精度是否满足要求, 必须量化和分析这些不确定性。生物量估算通常包括三方面的不确定性：测量不确定性、模型不确定性、抽样不确定性。Shettles等^[6]研究指出3种不确定性中, 模型不确定性所占比重较大, 约占总不确定性的70%。鉴于此, 本文仅研究不确定性较大的模型不确定性, 而测量和抽样不确定性, 本研究暂不考虑。模型不确定性来源主要有4个方面：输入变量本身的不确定性^[7-8]、模型函数形式的错误设定^[9]、模型残差变异^[10]以及模型参数误差^[11]等。输入变量本身的不确定性主要指变量(胸径、树高等)的测量误差, 测量误差主要受仪器精度、测量条件、以及人为等因素影响。与其余误差相比, 测量误差对单木生物量估算的影响相对较小^[12]。模型形式的错误设定主要是因为缺少合适的模型验证数据^[9]或者缺乏建模技术^[13]使得生物量模型的形式设定错误。随着生物量模型研究的不断深入, 生物量模型的形式也正逐渐确定。如Araújo等^[14]利用127棵实测样地数据对常用的14个生物量模型进行了验证分析, 结果表明方程式FW=aD^b/(1－M)和FW=aD^bH^c(式中FW为生物量鲜重, D为胸径, M为含水率, H为树高)对于数据的拟合效果最好。国内王轶夫等^[15], 蔡会德等^[16]也对不同模型的精度进行了对比。与前两种不确定性来源相比模型残差变异以及模型参数误差研究相对较少, 在国内更是少见。模型残差变异主要与模型的拟合精度有关, 通常模型的拟合精度以决定系数(R²)来表示, R²越接近于1, 表示模型的精度越好。模型的残差变异可以以残差的标准差来度量, 但是由于不同研究的研究区及数据不同对于研究结果也有一定的差异^[17-18]。与模型的残差变异相比模型的参数误差引起的不确定性相对较小, 而且模型参数误差与建模样本的数量有关^[17]。目前国内对于模型不确定性对森林生物量估算的影响研究较少, 傅煜等^[19-20]对森林生物量计算过程中的抽样误差和模型误差进行了研究。

本文以浙江省临安市森林资源连续清查数据为数据源, 以临安市为研究区, 通过建立了特定形式的一元、二元生物量模型来估算研究区杉木(Cunninghamia lanceolata)生物量, 并分析了两种模型形式下模型残差变异、模型参数误差带来的不确定性及模型总不确定性, 以期为提高森林生物量估计精度提供依据。

1 研究区概况与数据来源 1.1 研究区概况

临安市(118°51′—119°52′E, 29°56′—30°23′N)地处浙江省西北部, 隶属杭州市, 东西长约100 km, 南北宽约70.8 km, 总面积3126 km², 研究区地理位置如图 1所示。

临安市地形复杂、海拔高低悬殊, 因此气候基本呈垂直变化分布, 光、温、水地域差异明显。该市地处亚热带季风气候区, 温暖湿润, 四季分明, 光照充足, 雨量充沛, 全年平均气温16.4℃, 全年日照时数1847.3 h, 全年降雨量1628.6 mm。适宜各种林木生长, 植物资源丰富。

全市林业用地面积26.1047万hm², 活立木总蓄积量830.0101万m³, 森林覆盖率76.55%, 其中有林地面积24.0389万hm², 有林地蓄积量8230.5676万m³。森林资源分布格局为西高东低。植被区划属亚热带东部常绿阔叶林区, 森林类型主要有常绿落叶阔叶混交林、针叶阔叶混交林、常绿阔叶林、马尾松林、杉木林、毛竹林以及灌丛等类型。

1.2 样地清查数据

研究区地面调查数据为2009年临安市森林资源连续清查样地数据, 全市共有固定样地109个, 样地为面积0.08 hm²的正方形。所有样地共有样木8179棵, 其中杉木2535棵, 其统计特征见表 1。由于样地数据中缺乏树高因子, 因此文中树高通过陈丽聪等^[21]的树高曲线来求得。

1.3 生物量建模数据

建立杉木地上部分生物量模型所用数据为52棵浙江省杉木地上部分生物量实测数据。数据收集采用收获法对样木生物量进行测定。全部样木都实测胸径、地径, 将样木伐倒后, 测量其树干长度(树高), 分树干、树枝、树叶称鲜重, 并分别抽取样品带回实验室, 在85℃恒温下烘干至恒重, 称量样本干重, 算出含水率, 然后换算出样木各组分的干重, 汇总得到立木生物量。建模数据统计情况见表 2。

2 研究方法 2.1 模型残差变异不确定性计算

为了计算模型残差变异引起的不确定性, 本研究首先需要用52株杉木生物量实测数据拟合杉木单木生物量模型。再通过生物量预测值与实测值的差值来计算由残差变异引起的不确定性。设单木生物量实测值为B, 生物量预测模型为, 则B可用下式表示：

(1)

上式中, 为生物量模型中所含的参数预测值, x为生物量模型中的变量(胸径、树高), ε为生物量预测值与实测值的差即残差。对于f(x, α)常用的模型形式有一元模型和二元模型。本文所用的f(x, α)形式为：f(x, α)=aD^b和f(x, α)=aD^bH^c。得到各模型参数后, 残差ε可表示为：

(2)

对于模型残差变异引起的不确定性, 可以通过残差的标准差来衡量。本研究发现残差的标准差随单木生物量的增加而增加, 呈一定的线性关系。因此, 可以通过拟合残差的标准差与单木生物量的线性关系来计算残差变异不确定性, 即：

(3)

式中, σ_ε为残差的标准差, θ为拟合参数。在计算σ_ε时, 本文采用Hosmer和Lemeshow ^[22]所使用的六步法(six-step procedure)来计算。具体步骤可概括为：(1)将生物量实测数据按升序进行排列, 计算单株样木的生物量预测值, 并计算残差ε；(2)对测量数据进行分组, 每n棵树为一组, 若最后一组不足n棵则计入上一组。由于本文中建模数据为52棵, 样本数量较少, 因此分组时每5棵一组, 共分10组。(3)计算每组数据生物量预测值的平均值及残差的标准差σ_ε, 并根据与σ_ε的关系进行线性拟合。对于每组数据生物量均值及残差的标准差可以通过下式计算：

(4)

(5)

(6)

式中, 为每组单木生物量预测值的均值, 为每组生物量残差均值, 最后残差变异不确定性可以用以下模型表示：

(7)

θ为要拟合的参数。

2.2 模型参数不确定性计算

对于模型参数引起的不确定性, 可以采用一阶泰勒级数线性化处理对参数的不确定性进行计算, Box^[23]运用泰勒级数一阶展开项计算了非线性模型的参数误差, 取得较好的效果。若函数f(x)在x₀处存在n阶导数, 当x趋近于x₀时, f(x)可表示为：

(8)

那么, 对于生物量模型B= f(x, α)泰勒级数一阶展开式可以表示为：

(9)

式中, 为单木生物量模型的预测值, x为生物量模型中的变量, 为模型中参数估计值, f(x, α)中α可表示参数的真实值, 为参数的偏导数。式(9)中反映的是由参数引起的模型估测误差σ_p, σ_p可以近似表示为：

(10)

式中, , Beck和Arnol^[24]定义为灵敏度系数矩阵。var(α)为生物量方程中估计参数α的协方差矩阵。

2.3 不确定性合成

如果间接测量的量是直接测量所得到的各个测量值的函数。则间接测量的不确定度可表示为各个直接测量不确定度的合成不确定度^[25]。若间接测量的函数关系为：

(11)

则Z的不确定度可以表示为：

(12)

式中m_z表示总不确定度, m_xi表示第i个变量的不确定度。根据此式, 单木生物量计算总不确定性σ_tree可以表示为：

(13)

3 结果与分析 3.1 一元生物量模型不确定性计算

通过52棵杉木生物量实测数据对一元生物量模型进行参数求解, 最终拟合的一元生物量模型为：B=0.074D^2.371, 模型的决定系数R²=0.975。基于一元生物量模型估算的临安市杉木生物量密度为6.94 Mg/hm²。图 2左图为一元生物量模型残差标准差与单木生物量预测值的散点图, 从图中可以看出两者基本成线性关系, 残差标准差随着单木生物量的增加而增加。对于一元生物量模型残差变异引起的不确定性, 最后拟合模型为：, 模型的决定系数R²=0.937, 有较好的拟合效果。这表明对于一元生物量模型, 由模型残差变异引起的不确定性约为11.1%。对于模型参数引起的误差, 通过计算参数的方差-协方差矩阵并结合式(10)可以求得。经过计算一元生物量模型参数的方差协方差矩阵为, 将其带入式(10)中可得一元生物量模型参数误差引起的不确定性约占总生物量的14.4%。通过公式(13)可知, 一元生物量模型估算合成不确定性为18.12%。

对于二元生物量模型, 最后拟合的生物量模型为：B=0.086D^1.979H^0.419, 模型的决定系数R²=0.984。基于二元生物量模型估算的临安市杉木生物量密度为7.71 Mg/hm²。图 2右图为二元生物量模型残差标准差与单木生物量预测值的散点图。从图中可以看出, 二元生物量模型残差标准差与单木生物量预测值的关系与一元生物量模型有相同的趋势。残差变异不确定性模型为: , 模型决定系数R²=0.727, 表明对于二元生物量模型, 由模型残差变异引起的不确定性约为7.0%, 不确定性要小于一元生物量模型。对于二元生物量模型参数引起的误差, 二元生物量模型参数的方差协方差矩阵为, 将其带入式(10)中可得二元生物量模型参数误差引起的不确定性约占总生物量的8.53%。综合考虑模型残差及参数不确定性, 二元生物量模型不确定性约为11.03%。

3.2 样木数量对模型估计不确定性的影响

为了体现建模样本数量对模型引起的生物量估算结果的影响, 本文随机抽取了建模样本N=30, N=40, N=52时模型参数不确定性。表 3为建模样本分别为30, 40, 52时模型参数引起的不确定性。由表可知当建模样本由30增加到40再增加到52时, 一元模型参数不确定性相应为20.26%, 16.19%和14.4%。二元模型参数不确定性由13.09%减少到9.4%再减少到8.53%。说明建模样本对于模型参数引起的不确定性影响较为显著。随着样本数量的增加模型参数不确定性逐渐降低。

表 4为不同建模样本下模型残差变异引起的生物量不确定性。在计算不同建模样本数对模型残差变异不确定性的影响时, 不同建模样本间保持相同的分组数, 通过改变每组的样本进行讨论。表 4显示, 当建模样本分别为30, 42, 48时, 一元生物量模型的残差变异引起的不确定性分别为15.2%, 12.3%, 11.7%, 二元生物量模型的不确定性为13.3%, 9.4%, 8.7%。可见, 随着建模样本的增加模型残差变异不确定性也相应变小。

3.3 不同分组方式对模型残差变异不确定性的影响

由于本研究建模样本较少, 在分组计算模型残差变异时每组仅有5株, 样本数量较少。为了体现不同分组数据对模型残差变异的影响, 本文分别计算了每组8棵、10棵、15棵时残差变异不确定性。需要指出的是为了保证每种分组方式具有一定的组数, 分组时部分数据与上一组是重复的, 不同分组情况下重复的株数分别为3棵, 4棵, 5棵。不同分组情况下, 模型残差变异引起的不确定性如表 5所示。

从表 5可知, 当分组株数分别为8、10、15时, 一元生物量模型残差变异引起的不确定性分别为11.9%、11.9%、12.1%, 要高于每组5棵时的11.1%, 但相差不大。二元生物量模型残差变异引起的不确定性分别为9.2%、8.7%、9.4%, 高于每组5棵时的不确定性。不同的分组虽然对模型参数引起的不确定性有一定影响, 但影响较小。

4 结论与讨论

本文分别计算了生物量模型中所包含的2种不确定性, 从以上结果可得出如下结论：(1)在2种不确定性来源中模型参数不确定性对估计结果影响最大, 模型残差变异引起的不确定性较小。(2)模型残差变异不确定性、参数不确定性与建模样本有关, 随着建模样本的增加, 两种不确定性相应减少。在进行生物量建模时, 可以通过增加建模样本来提高单木生物量模型的估计精度, 减小模型不确定性对森林生物量不确定性的影响。

目前区域森林植被生物量估计研究是目前研究的热点问题, 该领域还面临着两大挑战^[26], 其中一个就是由于数据不准确、方法不得当, 使得森林生物量的估计存在众多不确定性。由于森林管理和决策的质量高低受森林生物量估计中所包含的不确定性影响, 生物量的估算迫切需要提高估测技术、数据质量来提高生物量估算的精度。本文分别计算模型各不确定性来源对森林生物量估算结果的影响, 其意义正在于此。与其他研究相比, 本研究还存在以下差异：(1)本研究模型残差变异的不确定性要小于国外现有研究。如Chen等^[17], Chave等^[18]曾对热带雨林地区生物量模型残差变异引起的误差坐了估算, 两者的结果分别为37.8%和31.3%。与本文的研究结果差距较大, 这可能由两方面原因造成：一是研究区与研究对象的不同, Chen、Chave中所用到的生物量模型为热带所有树种的生物量模型, 由于各树种之间木材密度或其他树种特征差异较大, 导致了较大的生物量残差。而本文生物量模型建模对象仅为单一树种——杉木。二是建模样本的差异。如Chen等^[17]所用的建模样本为4004棵, 胸径范围为5—212cm^[27], 而本研究中最大胸径为37.1cm远远小于Chen等建模样本。(2).本研究模型参数误差引起的不确定性较大。可能的原因为模型参数误差与建模样本数量有关。Chen等^[17]利用4004棵样木数据计算了模型形式为a(ρD²H)^b模型参数带来的不确定性, 研究表明当样木数据集足够大时由模型参数引起的不确定性很小仅为0.74%, 但样木数量减少到40棵时模型参数误差增加到5.69%。另外, 测算模型参数误差导致的不确定性有赖于生物量模型参数的协方差阵。本文杉木的建模数据52株, 由此推算整个研究区上万株林木的模型参数不确定性, 建模数据量可能偏小, 导致模型参数不确定性的度量结果偏大^[19]；而且本文生物量建模精度较高, 使得模型残差变异不确定性相应减小, 从而导致了模型参数误差引起的不确定性要大于模型残差变异引起的不确定性。

本研究生物量建模样本较少, 对于分组样本数对残差变异不确定性的影响本文最多讨论到每组15株的情况, 更高样本数对残差变异不确定性是否有影响尚需深入研究。