基于混合效应的兴安落叶松单木断面积生长模型研究

覃鑫浩

(国家林业和草原局 调查规划设计院，北京 100714)

准确描述林分或林木的生长规律，从而更好地为森林的经营与利用提供参考依据。描述生长规律主要包括前人积累的定性经验以及林分生长与收获模型[1]。根据模拟目的和对象的不同将森林生长和收获模型划分为林分-径阶-单木3个水平。单木生长模型的模拟对象是单株树木，将林木大小，竞争以及立地条件表述为单木生长量的数学表达式，单木模型能够描述复杂的林分关系，不仅能适用于人工纯林，其对异龄混交林和天然林也具有较强的适应性[2-4]。此外，与径阶模型和全林分模型相比，单木模型针对不同经营措施模拟得到的林分生长动态更准确，为森林质量的精准提升提供理论依据[5]。国内学者对不同树种的单木生长模型做了大量研究，邵国凡以红松人工林为研究对象，构建了红松(Pinuskoraiensis)的单木生长模型，为单木模型在我国的普及和应用奠定基础[6]；孟宪宇等[7]采用生长量修正法构建了华北落叶松(Larixprincipis-rupprechtii)人工林与距离无关的单木生长模型；吕理兴[8]以栓皮栎(Quercusvariabilis)天然林为研究对象，选用简单竞争指数并建立以个体树木生长为基础的、与距离有关的单木生长模型。

传统建模方法为最小二乘法，其假设数据独立、等方差、正态分布等特点，然而林分生长和收获模拟数据为连续多次重复观测数据，数据之间普遍存在异方差和自相关，因而不满足普通回归分析中的独立性假设[9-10]。混合效应模型能够降低分组数据的异方差性和自相关性，从而提高预测精度并解释随机误差的来源[11-12]。雷相东等[13]以天然落叶松云冷杉林为研究对象，研究发现考虑层次结构的混合效应模型的拟合效果优于传统模型，误差、均方根误差及其相对值均显著地减少；杜志等[14]考虑样地效应的马尾松(Pinusmassoniana)天然次生林单木胸高断面积生长模型优于传统模型；龚召松等[15]考虑样地效应的楠木(Phoebezhennan)次生林林分断面积与蓄积相容性混合效应生长收获预估模型在精度、平均误差和均方根误差都优于传统模型。

兴安落叶松(Larixgmelinii)是我国最为典型的东部高寒地区的优势树种，喜光性强，对水分要求较高，在内蒙古大兴安岭地区分布广泛，是该地区的地带性植被[16]，同时他也是我国经济建设主要用材树种[17-18]。以兴安落叶松为主要建群树种的森林对内蒙古大兴安岭林区的水源涵养、空气净化等功能发挥着重要作用。科学地经营兴安落叶松林才能更好地发挥大兴安岭林区的生态功能和价值，在经营目标转变和森林功能多元化影响下，及时准确地掌握森林资源现状和动态，开展以森林提升为目标的兴安落叶松林密度调控技术，提高内蒙古大兴安岭林区森林资源的现代化经营管理水平。

1 研究区概况

内蒙古大兴安岭林区为大兴安岭的西半部分，南北长696 km，东西宽384 km，地理坐标119°36′20″-125°20′50″E，46°08′40″-53°20′00″E，海拔425～1 760 m[19]。地处欧亚大陆中高纬度地区，气候类型为寒温带的大陆性季风气候，春季干旱且多风，夏季炎热多雨，年平均气温-2～4 ℃，年均降水量450 mm。土壤类型主要包括漂灰土、黑土、暗棕壤、灰灰黑土和黑钙土[20]。主要乔木树种是兴安落叶松、白桦(Betulaplatyphylla)、山杨(Populusdavidiana)和栎类等，植被类型主体属于欧亚针叶林植物区，以中、低山森林植被占优势，共有维管束植物1 339种，其中种子植物1 297种，蕨类植物42种。

2 数据与方法

2.1 数据

采用第六次至第九次内蒙古大兴安岭重点国有林区森林资源连续清查数据，以兴安落叶松天然林为研究对象，以单水平线性混合效应模型的方法，构建单木胸高断面积生长的混合效应模型。样地面积0.067 hm2。外业调查因子包括林分调查因子(样地的地理坐标(经纬度)、坡向、坡位、坡地、地貌、起源、龄组、优势树种、平均年龄、平均胸径、林种和起源等)和样木调查因子(树种名称、林木的胸径)。样地分布如图1所示。

图1 兴安落叶松天然林样地分布Fig.1 Distribution map of Larix gmelinii forest sample plots

2.2 研究方法

2.2.1 变量选择及基础模型的构建 林木生长的主要影响因子有林木的大小、林木所处的立地条件以及林木之间的竞争等，根据前人研究经验，本研究将林木生长描述为三者的线性函数[2-3]。采用林木大小、竞争及立地条件来描述单木胸高断面积的生长。其中林木大小为林木的胸径以及胸径的变形形式；林木竞争是林分生长必不可少的因素，对林木的生长趋势有重要的影响，本研究选择与距离无关的竞争指数参与模型构建，与距离无关的竞争指数主要包含以下4个变量，即每公顷株数(N,trees·hm-2),每公顷断面积(BA,m2·hm-2)，相对密度指数RD(对象木胸径与林分平均胸径之比：DBH/QMD)以及比对象木胸径大的断面积之和(BAL,m2·hm-2)；立地条件是影响林分生长的重要因子，本研究主要考虑海拔(EL)，坡度(SL)以及坡向(ASP)对于单木生长的影响。其中，采用Stage提出的公式对ASP和SL进行组合[1,13,21]，即SLcos=tan(SL)×cos(ASP)，SLsin=tan(SL)×sin(ASP)。

因变量的形式会对模型的表现产生很大的影响，在模型模拟前，分别采用不同的因变量对模型进行预测，根据残差图和评价指标，当自变量为ln(DBH22-DBH2)时，模型的模拟效果最好，因此选择ln(DBH22-DBH2)作为模型的因变量。

利用多元逐步回归方法，建立传统的兴安落叶松基础模型，根据树木的生长规律，选择对兴安落叶松生长影响性强的变量进入模型，并对模拟的结果采用多重共线性检验，最终选择对兴安落叶松生长影响显著，即回归系数显著(P<0.05)且方差膨胀因子(VIF)<5的变量进入模型。

筛选出兴安落叶松样地626块，样木株数60 279株。其中501块样地数据，兴安落叶松47 524株，用于构建兴安落叶松的胸高断面积预估模型，采用剩余的125块样地，兴安落叶松12 755株对构建的模型进行检验。建模数据统计如表1所示。

表1 兴安落叶松建模数据变量统计Table 1 Statistical table of L.gmelinii modeling data variables

由表1可以看出，兴安落叶松建模数据期初最小胸径为5 cm，最大胸径为74.1 cm；期末最小胸径为5 cm，最大胸径为79 cm；胸径的增量从0.1～5.0 cm；胸高断面积之和最大为32.08，最小为0.16，BAL最大值为2.14 m2，最小为0，林分每公顷株数最大为3 776，最小为15；坡度正切值变化为0.0～0.7，样地海拔分布266～1 650 m。

2.2.2 混合效应模型的构建 采用单水平线性混合效应模型的方法构建内蒙古大兴安岭重点国有林区兴安落叶松的单木胸高断面积生长预估模型，模型的一般形式如下[22]

(1)

式中:yi为第i个样地的单木断面积生长量；Xi为固定效应即与林木大小、竞争和立地条件相关的设计矩阵；β为由固定效应参数构成的向量；Zi为随机效应设计矩阵；bi为模型随机参数构成的向量；ei为误差向量；n为样地个数；D为随机效应的方差-协方差矩阵；σ2为模型方差；Ri为误差方差-协方差结构。

混合效应模型既有固定效应参数又有随机效应参数，其中随机效应包括随机截距模型、随机截距和斜率模型。本研究将基础模型的参数均作为随机效应参数，需计算随机截距和斜率模型，根据计算结果对不同的随机效应参数进行组合，并对不同组合进行模拟，最后筛选出最优的参数组合。

在混合效应的参数估算过程中，为了确定合适的方差协方差结构，对复合对称(CS)、一阶自回归[AR(1)]与一阶自回归移动平均结构[ARMA(1,1)]进行检验。采用指数函数(Exponent)、幂函数(Power)和常数加幂函数(ConstPower)对方差异质性进行描述。

2.2.3 模型的选择与检验

1)模型的筛选依据AIC(AIC)、BIC(BIC)和-2LL(-2Loglik,-2Loglik)越小越好原则，同时要对筛选出的含有不同随机参数个数的混合模型进行似然比检验(LRT)避免过多参数化问题的产生。

AIC=-2Loglik+2d

(2)

BIC=-2Loglik+dln(n)

(3)

LRT=2[log(L2)-log(L1)]

(4)

2)计算模型的随机效应参数值。

(5)

式中:D为随机效应方差-协方差矩阵；ZK为设计矩阵；RK为误差方差-协方差矩阵；ek为误差项。

3)采用平均绝对误差(Bias)、均方根误差(RMSE)和决定系数(R2)对模型的拟合效果进行评价。

2.2.4 统计处理 样地数据的筛选与处理在Excel软件中完成，基础模型和混合效应模型的构建以及模型的检验在R语言软件中完成。

3 结果与分析

3.1 兴安落叶松基础模型结果

基于内蒙古大兴安岭重点国有林区第六次至第九次4期森林资源连续清查数据，筛选出兴安落叶松样木株数47524株用于基础模型的构建，得到最终的预估模型

ln(DBH22-DBH2)=β1+β21/DBH+β3NT+β4BAL+β5SLcos+eij

(6)

式中:DBH表示期初胸径，DBH2表示期末胸径，β1、β2、β3、β4和β5表示固定效应参数，NT表示每公顷株数，BAL表示样地内比对象木大的所有胸高断面积之和，即SLcos=SL×cos(ASP)，eij表示模型误差。

模型模拟结果如表2所示。从表2可以看出，1/DBH、NT、BAL和SLcos对兴安落叶松单木断面积的生长有显著影响。兴安落叶松断面积生长模型的决定系数R2为0.398。根据构建的模型使用R(ggplot2)绘制兴安落叶松的残差图(图2)，从残差图发现模型存在异方差性。

表2 传统方法模拟兴安落叶松单木胸高断面积预估模型的结果Table 2 The results of L.gmelinii basal area growth model based on traditional regression

图2 基础模型的残差图和QQ图Fig.2 Residual plot and quantile-quantile (QQ) plot of the basic model

3.2 混合效应模型

在基础模型的基础上，引入样地效应，通过对常数项(int)、1/DBH、NT、BAL和SLcos 5个变量不同随机参数组合进行拟合，构建兴安落叶松断面积生长混合效应模型，拟合收敛的情况共有31种(表3)。通过比较拟合的AIC、BIC以及-2LL发现，相比于传统模型，基于样地效应的混合效应模型的拟合效果更好。其中，只含有1个随机效应参数的拟合方程，常数项(模拟2)的拟合效果最好；当拟合方程含有2个随机效应参数时，常数项和1/DBH的组合(模拟7)拟合效果最好；当包含3个随机效应参数时，常数项、1/DBH和BAL的组合(模拟18)拟合效果最好；当含有4个随机效应参数时，常数项、1/DBH、BAL和SLcos的组合(模拟29)拟合效果最好；当含有常数项、1/DBH、NT、BAL和SLcos共5个随机效应参数时，得到模拟32。当拟合方程包含常数项、1/DBH、NT、BAL和SLcos(模拟32)的拟合效果优于其他任意一个组合的拟合效果，模拟精度最高。

对于包含不同参数的5种最优模型，似然比检验结果(表3)发现，模拟32>模拟29>模拟18>模拟7>模拟2，随着随机效应参数个数的增加，模型的模拟效果越好。最终选择模拟32(式7)作为最佳参数组合的混合效应模型。

表3 胸高断面积生长量线性混合模型收敛结果Table 3 The lack-of-fit statistics of the linear mixed effects individual-tree basal area growth model

ln(DBH22-DBH2)=(β1+b1)+(β2+b2)1/DBH+(β3+b3)NT+(β4+b4)BAL+(β5+b5)SLcos+eij

(7)

式中:β1～β5是固定效应参数；b1～b4是随机效应参数。

由于数据存在异方差性，本研究用指数函数、幂函数和常数加幂函数作为异方差结构消除数据的异方差性。表4显示，混合效应模型的3个异方差结构均显著改善了AIC值、BIC值和Loglik值，提高了模型的预估能力。从表4计算的LRT结果发现，指数函数在消除数据的异方差性中表现最佳，作为该模型最终的异方差函数。此外，由于数据为多期重复测量数据，因此存在自相关性。采用3种常用的自相关结构对数据的相关性进行表达，结果表明除ARMA(1,1)未收敛外，其余2种自相关结构均改善了模型的表型，其中AR(1)表现最好，确定为模型最优自相关结构。

表4 考虑异方差函数和自相关结构后混合效应模型结果比较Table 4 The lack-of-fit statistics of the linear mixed effects individual-tree basal area growth model using different error variance-covariance structures

确定了参数效应和误差方差-协方差结构后，基于限制性极大似然法(REML)，得到最终的兴安落叶松单木胸高断面积生长混合效应模型

ln(DBH22-DBH2)=(4.477 5+b1)+(-11.440 8+b2)1/DBH+(-0.000 2+b3)NT+(0.749 2+b4)BAL+(0.066 6+b5)SLcos+eij

(8)

图3为混合效应模型的残差图和QQ图。通过对比发现，模型残差图在引入随机效应之后得到了改善。

图3 混合效应模型的残差图与QQ图Fig.3 Residual plot and QQ plot of our final mixed effects model

3.3 模型检验

采用剩余的125块样地、12 755株兴安落叶松对构建的混合模型与传统模型进行检验，比较模型的均方根误差、平均绝对误差和决定系数。基础模型的检验只需要考虑固定效应，将自变量参数代入到模型中，与实测值进行比较即可；考虑到混合效应模型的固定效应和样地的随机效应，本研究将获得的固定效应参数与随机效应参数求和，带入模型得到拟合值并与实测值进行比较，判断模型的模拟效果。结果如表5和图4所示，传统最小二乘方法的3个统计值计算结果分别为平均绝对误差0.671 9，均方根误差0.820 7，模型的决定系数R2为0.353 7。考虑样地效应的混合效应模型的平均绝对误差为0.588 9，均方根误差为0.730 0，模型的决定系数R2为0.488 6。结果显示考虑样地效应的混合效应模型的平均绝对残差与均方根误差都比基础模型的低，模型的决定系数R2由0.353 7提高到0.488 6。模型检验结果证明考虑样地效应的混合效应模型的模拟效果优于一般线性模型。

图4 兴安落叶松胸高断面积基础模型与混合效应模型观测值与预测值Fig.4 The observed and predicted values of the basic model and mixed effects model for L.gmelinii

表5 传统模型与混合效应模型验证拟合统计值比较Table 5 The lack-of-fit statistics of model validation using the basic model and mixed effects model

4 结论与讨论

以内蒙古大兴安岭地区兴安落叶松天然林为研究对象，构建了兴安落叶松天然林中兴安落叶松的胸高断面积生长模型。林木的大小即胸径的倒数(1/DBH)、每公顷株数(NT)、比对象木大的胸高断面积之和(BAL)、坡度和坡向的组合(SLcos)4个参数明显影响了兴安落叶松天然林的生长效果，其中林分生长与胸径的倒数呈负相关，即林木胸径随个体胸径倒数的减小而增大，因此提升了林木获取资源的能力，促使林木生长量增大。杜志等[14]以湖南地区马尾松为研究对象，王建军等[1]以福建省三明市杉木为研究对象，Zhao等[23]以硬阔树种为研究对象，R.Cannell等[24]以云杉(Piceaasperata)和松树为研究对象构建各个树种的生长模型发现了相同的规律。林分密度(NT)表示单位面积林木株数，它能间接反映单位面积林木间的竞争，胸高断面积生长量与每公顷株数成反比，单位面积林木株数越多，获取光和水分的几率就越小，树木间的竞争压力越大，胸高断面积生长量越小。BAL是比对象木大的所有林木的胸高断面积之和，BAL值越大则表明林分中比对象木大的树木越多，对象木受到的竞争压力越大，因此与胸高断面积生长量成反比。坡度和坡向的组合(SLcos)对林分的生长表现显著。由于森林监测数据具有连续重读观测的特点，与传统方法会得到有偏估计，本研究通过对不同的随机效应参数进行组合，筛选最优的参数组合和误差方差-协方差结构提高了模型预测精度。与传统模型相比混合效应的决定系数从0.353 7提高到0.488 6，误差、均方根误差及其相对值均显著减少，建立的模型有一定的实用性。本研究以兴安落叶松天然林为研究对象，探究主要树种胸高断面积生长模型，发现林木大小、林地立地条件以及林木之间的竞争显著影响着林木的生长，这个与很多专家的研究结果一致，并且上述因变量在林业调查工作中易于获取。因此，本研究的兴安落叶松胸高断面积生长模型具有一定的生物学意义，统计结果可靠。可为以后兴安落叶松林的采伐、经营等提供基础。

由于建模的数量较大，数据跨距的南北范围较广，温度、降水、土壤类型等环境因子可能对树木生长存在潜在的影响。因此，在未来的研究中有必要考虑气候因子，以更精准、更全面地描述树木的生长。此外，构建模型的竞争指数时只考虑了林木大小的相关参数没有把树种组成进行讨论在以后研究中可以将树种组成加进来。另外，混合效应模型没有明确的区域划分，在以后研究中要考虑区域水平、样地水平以及2水平的混合效应模型，比较不同混合效应的模型，找出最优的混合效应模型。