异速参数先验信息对生物量方程拟合的影响

2021-09-23黄兴召豆玉萍黄庆丰陈东升

自然保护地 2021年1期

黄兴召豆玉萍黄庆丰陈东升

(1. 安徽农业大学林学与园林学院，合肥 230036；2. 中国林业科学研究院林业研究所，北京 100091)

立木地上生物量是森林生物量的重要组成部分，如何准确地估算立木地上生物量是研究森林生态系统物质循环的基础[1−2]。异速生物量模型作为预估生物量的常用方法，已被广泛应用于各类研究中[3−6]。异速生物量方程参数a和b随着研究区域和树种的变化而变化，通过收获法，利用异速生物量方程预估生物量时，只有单独建模才能提高预估精度[5−8]。贝叶斯方法作为其中的一种拟合方法，通过收集已建立方程的参数信息，将该信息转化为参数的概率分布，以参数的均值和协方差矩阵为先验信息进行拟合，越来越受到人们的关注[9−12]。但是，利用先验信息对贝叶斯方法进行拟合时，研究者收集参数a和b的先验信息具有随机性和异质性，由此对模型拟合和预估效果产生的影响尚未得到验证，也没有阐述如何选取参数a和b的先验信息来提升方程的拟合和预估效果。

收集已发表文献异速生物量方程参数a和b的信息，是分析先验信息异质性对立木生物量预估影响的基础。Ter-Mikaelian和Korzukhin[13]建立了北美65个树种的生物量预估方程、Zianis和Mencuccini[3]建立了279个异速生物量预估方程、Henry等[14]收集了850个立木及各组分的生物量预估方程、Zapata-Cuartas等[8]收集了134个的立木地上异速生物量方程、Rojas-García等[15]收集了487个立木及各组分的生物量预估方程。这些方程都是通过常规方法拟合得到参数a和b的值，为贝叶斯方法提供了丰富的先验信息。分析发现，立木地上部分生物量预估方程参数b与纬度呈显著相关，且在不同属之间存在显著差异[8,13]。研究者在利用贝叶斯方法预估立木地上部分生物量的过程中，讨论了样本大小对方程拟合效果和预估效果的影响，缺乏反映树种特征的先验信息的拟合比较[8,10,16]。树木生长与生物量的积累受环境[17,18]和生物学特性[19−20]的影响，在方程拟合时先验信息也会受到影响[21]。因此，检验参数a和b先验信息的气候区和生物学特性对立木生物量预估的影响显得尤为重要。

通过分析已发表文献关于立木地上部分异速生物量方程的参数a和b信息，将参数的先验信息按主要气候区和树种进行分类，以日本落叶松（Larix kaempferi）立木地上部分生物量为拟合数据，利用贝叶斯方法拟合异速生物量方程，比较先验信息的气候差异和树种差异对立木生物量拟合和预估的影响，探讨先验信息对立木生物量预估的作用规律，为利用贝叶斯方法拟合异速生物量方程时选择先验信息提供了理论和实践的依据。

1 研究方法

1.1 异速生物量方程

异速生物量方程作为预估立木生物量的常用模型，其方程表达式为：

其中，Mi为第i株样本的地上部分生物量；Di为第i株样本的胸径；a和b为方程参数；ei为误差项。

贝叶斯方法将公式（1）中的参数a、b定义为θ，需要知道θ的概率密度函数。接着，根据先验信息确定参数a和b的先验分布 π(θ)，再结合样本信息可得出后验分布 π(θ|D)。 π(θ|D)的计算公式为：

1.2 先验信息

利用贝叶斯方法拟合异速生物量方程时，以参数a、b的均值和协方差矩阵为先验信息。为研究不同先验信息对异速生物量方程拟合效果的影响，本文通过分析[8,13]立木地上部分异速生物量方程参数a和b的数据集，按照气候和属进行分类（表1）。以主要气候区为单位，将参数a和b分为寒带、温带和热带三种类型，分别计算它们的均值和协方差矩阵。同时，以属为单位，将参数a和b分为栎属（Quercus）、桦木属（Betula）、杨属（Populus）、槭属（Acer）、桉属（Eucalyptus）和松属（Pinus）等六个属，计算其均值和协方差矩阵。以气候区和属为单位进行分类，是进一步研究先验信息异质性对异速生物量方程拟合效果影响的关键。

表1 异速生物量方程先验信息分类Table1Priorinformationclassificationofallometric biomassequation

2 方程拟合

2.1 拟合数据

以辽宁省清源县大孤家林场（42°16ʹ～42°23ʹ N，124°48ʹ～124°55ʹ E）日本落叶松（Larixkaempeferi）立木地上部分生物量数据为拟合数据[22]，该区属于典型的温带季风气候。基于样地调查，在标准地附近选取标准木，利用收获法测定样木的树干、树皮、树枝和树叶各组分的生物量，求出立木地上部分生物量。本数据由幼龄林、中龄林、近熟林以及成熟林4个林龄段共45株立木的胸径和地上部分生物量组成（表2）。本文在建模过程中，利用无重复抽样，从45株总样木中抽取30株样木，进行1 000次无重复拟合，每次拟合未被抽取的15株作为验证数据来检验预估效果。

表2 日本落叶松样地信息Table2PlotsinformationofLarixkaempferiplantations

2.2 预估效果与预估精度检验

使用不同的先验信息，利用贝叶斯方法拟合异速生物量方程时，采用决定系数（R2）、平均偏差（Mean Bias , MB）和平均均方根差（Mean Root Mean Square Error, MRMSE）来检验方程的拟合结果和预估效果。R2、MB和MRMSE的表达式分别为：

其中，yi为实际观测值；y¯为样本平均值；yˆi为模型预测值；n为样本数。数据统计分析和绘图使用R（Development Core Team, 2009）和Excel软件。以不同区域和不同树种参数a、b的均值以及协方差矩阵为先验信息，应用R的MCMCglmm包[23]对异速生物量方程分别进行拟合。

3 结果与分析

3.1 先验信息气候区差异对异速生物量方程拟合的影响

为检验先验信息区域差异对异速生物量方程拟合的影响，分别以寒带、温带、热带以及总数据集参数a、b的均值及协方差矩阵为先验信息，利用无重复抽样，从45株总样木中抽取30株样木，进行1 000次无重复拟合（图1）。寒带和热带先验信息拟合后参数a、b值的变化范围相似；温带和总数据集先验信息拟合后参数a、b值的变化范围相似。以热带和寒带为先验信息进行1 000次无重复拟合，参数a的平均值均为−2.794，95%的置信区间分别为（−3.185, −2.386）和（−3.197,−2.368）；参数b的平均值均分别为2.640和2.638，95%的置信区间分别为（2.484, 2.788）和（2.476, 2.792）。以温带和总数据集为先验信息进行1 000次无重复拟合，参数a的平均值分别为−2.836和−2.826，95%的置信区间分别为（−3.240,−2.418）和（−3.228, −2.407）；参数b的平均值均分别为2.655和2.651，95%的置信区间分别为（2.496, 2.809）和（2.491, 2.805）（表3）。这说明使用不同区域的先验信息拟合异速生物量方程会影响参数a、b值的稳定性。

图1 不同气候区先验信息预估参数a和b的拟合结果Fig.1Simulationresultsfortheestimationofthe parametersaandbusingdifferentprior informationofregion

决定系数R2反映方程的拟合效果（表3），寒带、温带、热带以及总数据集四种先验信息1 000次无重复拟合，方程的平均决定系数均为0.934。平均偏差MB和平均均方根差MRMSE反映了方程的预估效果，四种先验信息1 000次无重复拟合，平均偏差分别为−1.392、−0.597、−0.929、−0.715，平均均方根差分别为0.347、0.344、0.347、0.345。由此可知，区域先验信息主要影响方程的预估效果。

表3 不同气候区先验信息拟合异速生物量方程的结果对比Table3Resultcomparisonallometricbiomassequationsindifferentpriorinformationofregions

为了进一步比较寒带、温带、热带以及总数据集四种先验信息对方程预估效果的影响，对BE和RMSE分别进行方差分析（图2）；总数据集先验信息的MB和MRMSE与温带的MB和MRMSE无显著差异；寒带与热带的MB和MRMSE也无显著差异；其他先验信息的MB和MRMSE之间存在显著差异（P<0.01）。总数据集和温带先验信息的拟合效果优于寒带和热带，验证数据集所在区域为温带，因此，利用贝叶斯方法拟合异速生物量方程时，应充分考虑先验信息的气候区特征。

图2 不同气候区先验信息拟合的平均偏差和平均均方根差比较Fig.2Simulationresultsforthemeanbiasandmean rootmeansquareerrorusingdifferent priorinformationofregion

3.2 先验信息属差异对异速生物量方程拟合的影响

以栎属、桦木属、杨属、槭属、桉属、松属和总数据集为先验信息，在样木数量为30的情况下，也进行1 000次无重复拟合（图3）。桉属、栎属、松属以及总数据集的参数a的变化范围相似，其平均值分别为−2.859、−2.866、−2.855、−2.826；桦木属和杨属的参数a的变化范围相似，其平均值分别为−2.728和−2.712；槭属的参数变化范围尤为突出，其平均值为−2.492。参数b的变化范围也表现出相同的规律，桉属、栎属、松属以及总数据集参数b的平均值分别为2.663、2.666、2.662、2.651；桦木属和杨属参数b的平均值分别为2.616和2.605；槭属参数b的平均值为2.528。

图3 不同属先验信息预估参数a和b的拟合结果Fig.3Simulationresultsfortheestimationofthe parametersaandbusingdifferent priorinformationofgenus

以栎属、桦木属、杨属、槭属、桉属、松属和总数据集为先验信息，进行异速生物量方程拟合（表4）表4。槭属的平均决定系数为0.935，其他属的平均决定系数均为0.934。栎属、桦木属、杨属、槭属、桉属、松属和总数据集的平均偏差分别为−0.205、−1.287、−2.926、−3.726、−0.367、−0.494、−0.715；平均均方根差分别为0.343、0.352、0.354、0.391、0.344、0.343、0.345。由此可知，树种先验信息主要影响方程的预估效果。

表4 不同属先验信息拟合异速生物量方程的结果对比Table4Resultcomparisonallometricbiomassequationsindifferentpriorinformationofgenus

为了进一步比较六个属以及总数据集七种先验信息对方程预估效果的影响，对MB和MRMSE分别进行方差分析（图4）。栎属和桉属、松属和桉属之间的MB无显著差异，其他属之间的MB均存在显著差异（P<0.01），即栎属和桉属为先验信息的MB>松属>总数>桦木属>杨属>槭属。栎属和桉属、松属和桉属、总数据集和桉属之间的MRMSE无显著差异，其他属之间的MRMSE存在显著差异（P<0.01），即松属和栎属为先验信息的MRMSE>总数据集和桉属>桦木属>杨属>槭属。因此，树种的生物学特性也是影响生物量预估的重要因子。

图4 不同属先验信息拟合的平均偏差和平均均方根差比较Fig.4Simulationresultsforthemeanbiasandmean rootmeansquareerrorusingdifferent priorinformationofgenus

4 讨论

4.1 区域差异对异速生物量方程拟合的影响

利用贝叶斯方法拟合生物量方程，先验信息尤为重要[8−9,16]，通过收集立木地上部分异速生物量方程参数信息，分析表明参数a和b符合二元正态分布。拟合异速生物量方程过程中，贝叶斯原理使用先验信息，而传统方法只能使用立木数据。因此，贝叶斯方法可以实现少量的样木高精度预估生物量，黄兴召等[10]和Zell等[16]均证明，当样木的数量小于50株时，贝叶斯方法的预估效果明显优于传统方法。立木地上部分异速生物量方程参数a和b的值有明显的区域特征，参数b与纬度呈线性关系[10]。本文将参数a和b划分为寒带、温带和热带三个主要气候区，基于信息量度的相似概率分布[22,24−25]，以三个气候区和总数据集参数a、b的均值和协方差矩阵为先验信息，利用无重复抽样，从总样木中抽取30株立木生物量数据，进行1 000次无重复拟合，得出每次拟合结果中参数a、参数b、R2、MB和MRMSE的值。

将寒带、温带和热带三个主要气候区先验信息的拟合结果与总数据集先验信息的拟合结果进行比较，发现R2均为0.934无显著差异，说明先验信息的气候差异对异速生物量方程的拟合效果无影响。MB和MRMSE差异分析表明，先验信息的气候差异对异速生物量方程的预估效果有显著影响。本文中拟合数据的研究地点属于典型的温带季风气候，温带的先验信息拟合效果（显著）优于寒带和热带，与总数据集无差异，验证了先验信息的气候差异对立木生物量预估的影响，使用立木生长主要气候区的先验信息，有助于提升立木生物量的预估精度。

通过模型拟合，证实了先验信息的气候区差异性对立木生物量预估的影响。气候是环境因子最直接的外在表现，其综合了降雨、光照、地形地貌等环境因子，决定了森林生物量和生产力的大小[26−27]。利用异速生物量方程进行预估时，相同胸径不同区域的立木生物量相差达到500倍[3,28]，相同区域内的立木生物量相差也达到99倍[14,29]。Deans等[30]和Henry等[31]分别利用异速生物量方程估算热带森林生物量，发现同一气候区内立木生物量预估方程也存在差异，如何通过先验信息解释预估方程的差异是下一步的研究重点。另外，气候区之间的过渡带森林的生物量和生产力均有增加的趋势，进一步验证先验信息区域差异对过渡区生物量预估的影响，对阐述气候对森林生物量和生产力的影响有重要的作用。

4.2 树种生物学特性对异速生物量方程拟合的影响

大量研究从针叶林和阔叶林[3,13]、天然林和人工林[13,15]、常绿林和落叶林[14,32]等森林类型来比较立木生物量和预估方程的差异。树种之间生物量的差异多数从木材密度来分析对比不同演替过程中木材密度的差异对立木生物量的影响[14,33]。本文利用栎属、桦木属、杨属、槭属、桉属和松属六个属的先验信息拟合异速生物量方程，从先验信息的生物学特性分析方程的预估效果。结果表明，槭属平均决定系数R2为0.935，其他的平均决定系数R2均为0.934，树种先验信息的差异对方程R2无显著影响。平均偏差MB分别为− 0.205、−1.287、−2.926、−3.726、−0.367、−0.494；平均均方根差MRMSE分别为0.343、0.352、0.354、0.391、0.344、0.343，树种先验信息的差异对MB和MRMSE影响显著（P<0.01），与先验信息的气候区差异的结果一致。

黄兴召等[22]研究结果表明，a和b的值在针叶树和阔叶树之间有明显差异，本研究结果显示以松属和栎属为先验信息的拟合效果>总数据集和桉属>桦木属>杨属>槭属。栎属、桦木属、杨属、槭属、桉属和松属分别属于壳斗科（Fagaceae）、桦木科（Betulaceae）、杨柳科（Salicaceae）、槭树科（Aceraceae）、桃金娘科（Myrtaceae）和松科（Pinaceae）。但是，落叶松属与松属的生物学关系>落叶松属与桉属>落叶松属与槭属>落叶松属与杨属>落叶松属与栎属>落叶松属与桦木属。由此推断，树种先验信息影响异速生物量的预估效果，与树种之间的生物学关系不一致，其作用规律有待进一步研究。本文选择先验信息的气候区差异和树种差异对立木生物量预估的影响，验证了气候因子和生物学特性对先验信息的影响和方程预估的作用，但先验信息差异性对交错区立木生物量预估的影响及与树种间生物学特性的关系，还有待进一步研究。