基于混合效应的飞播马尾松林单木冠幅预测模型

杨 阳 ，彭浩贤，潘 萍，欧阳勋志*，臧 颢，纪仁展，余 枭

（1.鄱阳湖流域森林生态系统保护与修复国家林业和草原局重点实验室，江西南昌 330045；2.江西农业大学林学院，江西南昌 330045；3.江西利森林业技术服务有限公司，江西南昌 330038）

【研究意义】林冠作为林分垂直分布的顶层，对林分内部养分循环、光能分配及降雨截留起着重要作用[1]，其生长状况决定了林木生产力及生活力的大小，冠幅是林冠结构中重要的特征因子[2]，冠幅的调查有助于调控林分生态效益、评估森林质量，同时对林分的最佳密度控制以及生物多样性的管理等具有重要意义[3]，而由于调查冠幅的成本较高且费时费力，因此，采用冠幅预测模型获取冠幅信息得到了较广泛的应用。【前人研究进展】构建单木冠幅模型的方法众多，多数学者利用传统回归方法构建了冠幅和林木因子（树高、枝下高、树高和胸径的比值等）[4-5]或林分因子（平均胸径、优势木平均高、株数密度等）[6-7]的冠幅-胸径模型。然而，建立回归模型时随机误差可能会产生异方差[8]，应用混合效应模型确定组内方差-异方差结构则能够消除模型异方差[9-10]。近年来，国内外不少学者采用混合模型建立冠幅模型，如Shar⁃ma 等[11-12]以挪威云杉（Picea excelsa）、欧洲山毛榉（Fagus sylvatica）以及欧洲赤松（Pinus sylvestris）为研究对象，在幂函数的基础上引入不同林分变量，将样地作为随机效应因子构建了3个树种的混合效应模型，结果显示混合模型的预测精度较基础模型有很大的提高；Xu 等[13]对比了杉木（Cunninghamia lanceolata）冠幅的线性混合效应模型与普通线性冠幅模型的拟合与预估效果，结果表明混合模型比普通模型有更好的拟合和估计效果；Fu 等[14]利用嵌套2 水平非线性混合效应模型构建的华北落叶松（Larix principisrupprechtii）冠幅预测模型也得到了较好的预测效果。总体来看，添加相关因子及随机效应均能有效提升冠幅模型的预测精度。【本研究切入点】马尾松（Pinus massoniana）是我国南方地区分布最广的针叶乡土树种之一，具有耐干旱、耐贫瘠、适应性强等特点。江西省兴国县曾是我国水土流失严重区域之一，曾被称为“江南沙漠”，自20世纪70年代到90 年代初进行了大面积的飞播马尾松造林，而这种基于曾经强烈水土流失背景下恢复起来的马尾松飞播林与其人工植苗及天然马尾松林相比，其生长、林分结构等都存在一定程度的差异，林分密度差异也较大，肖欣等[15]、褚欣等[16]对不同密度飞播马尾松林的林下植被特征、土壤质量及持水性能等研究得出不同密度组间存在不同程度的差异，但其冠幅生长与林木密度等因子间的关系目前尚不清楚。【拟解决的关键问题】本研究基于混合效应模型，构建飞播马尾松林单木冠幅预测模型，探讨林分密度等解释变量和随机效应对冠幅模型的影响，以为飞播马尾松林的结构调整及科学经营提供参考依据。

1 材料与方法

1.1 研究区概况

兴国县（26°03′～26°42′N，115°01～115°51′E）地处江西省中南部，赣州市北部，属亚热带季风湿润气候带，年均温18.9 ℃，年均降雨量1 539 mm，无霜期280～300 d。母岩主要为花岗岩、第四纪红色黏土、砂岩等，土壤类型主要为红壤。主要森林类型有马尾松林、杉木林、竹林和常绿阔叶林等，其中飞播马尾松占全县有林地面积的29.5%[17]，但由于受到人为干扰及病虫害的影响，早期的飞播林保留甚少，现存的飞播马尾松林以20～30年的为主。

1.2 标准地设置与调查

2014 年7 月，在对飞播马尾松分布区进行踏查的基础上，选择人为干扰程度低、下坡位20～30 年林分密度组I（900～1 500 株/hm2）、密度组II（1 500～2 100 株/hm2）、密度组III（2 100～2 700 株/hm2）、密度组IV（2 700～3 300株/hm2）分别设置标准地，每个密度组4个重复，共16个标准地，标准地面积为30 m×30 m。对胸径大于2 cm 的乔木进行每木调查，调查内容包括胸径（DBH）、树高（H）、枝下高（HCB）、冠幅（CW）、株数密度（M）、标准地优势木平均高（TH）以及立地等因子，在标准地中选择5 株树高长势较好的林木作为优势木，取其算术平均值作为TH。16 块标准地共计马尾松2 354 株，分别在各密度组随机选取3 块共12块标准地的1 732株马尾松作为建模数据，剩下的4块标准地的622株马尾松作为检验数据，各密度组标准地基本信息见表1。

表1 不同密度组标准地主要调查因子Tab.1 Main investigation factors of different density groups

1.3 研究方法

1.3.1 基础模型的筛选 选取10个常用的冠幅-胸径模型作为备选基础模型（表2），分别进行拟合与检验，采用决定系数（R2）、均方根误差（RMSE）和平均绝对误差（MAE）选取拟合效果最好的模型作为基础模型。

1.3.2 广义模型的构建 采用再参数化的方法向基础模型中依次添加与冠幅相关性高的林木因子和林分因子构建广义模型。结合前人研究[6，18]，考虑的林木因子有H、HCB、树高与胸径的比值（HDR）、标准地中大于对象木胸径所有林木胸径之和（LDSD）及断面积之和（BAL）；林分因子包括标准地平均胸径（MD）、TH、M、密度组（lvM），考虑到密度组为定性因子，参考相关研究[19]，将密度组I、II、III、IV 赋值为1、2、3、4进行运算。

1.3.3 混合效应模型的构建 基础模型和广义模型确定后，进一步考虑随机效应对马尾松冠幅-胸径的影响，构建混合效应模型。混合效应模型的一般形式如下[18]：

表2 冠幅-胸径备选模型Tab.2 Crown-diameter alternative models

式（1）中，Yi是由冠幅实测值组成的因变量向量，Xi是由胸径及添加变量组成的自变量向量；f为关于向量参数φij和Xi的非线性函数；β是固定效应参数；ui是服从期望为0，方差-协方差矩阵为ψi且正态分布的独立随机效应向量；φij是与β、ui呈线性关系的形式参数；εij为随机误差项，假定服从期望为0，方差为Rij的正态分布，且假定ui与εij之间相互独立，Aij、Bij分别为β和ui的设计矩阵。

1.3.4 模型参数及方差协方差结构的确定 参考Pinheiro 等[20]的研究，根据所有固定效应参数组合添加随机参数对模型进行拟合，利用赤池信息量准则（AIC）、贝叶斯信息准则（BIC）和对数似然函数值（LL）在收敛的模型中确定最优随机效应模型，AIC、BIC值越小，LL值越大，模型拟合效果越好，同时使用似然比检验（LRT）避免模型过度参数化。

式（2）中，d为模型的参数个数；n与ln（n）分别为样本数量及其自然对数；L1、L2分别为比较的2 个模型的最大似然函数。

在模型构建过程中，数据间可能存在明显异方差和自相关，常用Ri描述方差异方差结构，表达式如下[21]：

式（3）中，σ2为模型的误差方差值；Γi为组内误差自相关结构；Gi为描述方差异致性的对角矩阵。

考虑到本研究数据为非连续测量数据，无需考虑时间序列上的相关性，数据间的异方差大多是通过对残差方差增加权重消除，故选用幂函数和指数函数来矫正模型中的异方差，其结构表达式如下[22]：

式（4）中，γ为待估参数。

1.3.5 模型检验与评价 基础模型与广义模型的检验只需要计算固定参数估计值，将检验样本带入模型中比较评价指标，而混合模型的检验包括固定效应和随机效应检验，其中固定效应部分相当于传统回归方法，随机效应部分的检验，则需要通过二次抽样计算随机参数检验随机效应，本研究使用最优线性无偏估计（EBLUP）来计算随机参数，计算公式如下[23]：

本研究采用决定系数（R2）、均方根误差（RMSE）、平均绝对误差（MAE）3个指标评价各模型的拟合与检验效果，计算公式参见文献[24-25]。

1.3.6 数据分析 用1stopt 软件对基础模型进行固定效应参数初步估计，基础模型的确定、混合模型的拟合及图表的绘制均是用R软件处理，其中混合效应模型通过nlme包计算。

2 结果与分析

2.1 基础模型

基于建模数据和检验数据，采用10 个基础模型进行拟合与检验，结果见表3。表3 可以看出，幂函数、二次型函数与Logistic2 模型的拟合效果较其他模型略优，其中幂函数预测效果最好，故选择幂函数作为基础模型：

式（7）中，CWij和DBHij分别为第i块标准地内第j株树的冠幅（m）和胸径（cm）；β0、β1为模型参数；εij为模型误差。

表3 基础模型建模及检验结果Tab.3 Results of the fitting and test of the basic models

2.2 广义模型

选择与冠幅相关性较高的H、HCB、HDR、M、lvM和TH作为添加变量，依次向所筛选出的基础模型中添加以上变量构建马尾松冠幅-胸径广义模型，为了避免变量过多导致多重共线性，运用逐步回归法剔除回归系数不显著（P>0.05）的解释变量，筛选出拟合优度最高的模型。收敛的拟合结果见表4，从表4可以看出，考虑1 个变量时，增加HCB的MHCB拟合优度最高（R2=0.599 2，RMSE=0.669 2，MAE=0.534 0），与基础模型差异极显著（F=106.47，P<0.01）；考虑2 个变量时，增加HCB和lvM的MHCB+lvM拟合效果最好（R2=0.601 0，RMSE=0.667 6，MAE=0.531 8），F检验结果为MHCB+lvM与MHCB差异极显著（F=8.137 1，P<0.01）；当考虑的变量超过2个时，收敛的模型与MHCB+lvM均不显著（P>0.05），故选择MHCB+lvM作为最优广义模型，其模型表达式为：

2.3 混合效应模型

将不同组合形式的随机效应添加到基础模型与广义模型中，构建混合效应模型，拟合结果见表5。从表5中可以看出，考虑1个随机参数时，在基础模型中加入u0的混合模型拟合效果较好，AIC、BIC和LL值分别为3 466.007、3 487.835、-1 729.004，似然比检验表明其与基础模型差异极显著（LRT=169.390 0，P<0.001）；在广义模型中加入u0的混合模型拟合效果较好，AIC、BIC和LL值分别为3 347.13、3 379.87、-1 667.57，较广义模型差异极显著（LRT=180.63，P<0.001）；考虑2 个随机参数时，含u0、u1的基础混合模型拟合效果较好，较含u0的基础混合模型差异极显著（LRT=16.0030，P<0.001），而广义模型中添加2 个以上随机参数时，收敛的模型与考虑1 个随机参数的混合模型均无显著性差异（P>0.01），以上结果可以说明随机参数可以显著提高模型的拟合优度。综上所述，考虑混合效应后对应的最优基础混合模型和广义混合模型表达式如下：

表4 广义模型拟合结果Tab.4 Fitting results of generalized models

式（9）中，u0、u1为随机变量。

表5 混合模型拟合结果Tab.5 Fitting results of mixed models

通过最优基础混合模型和广义混合模型的残差分布图（图1）发现，混合模型在拟合过程中残差呈现喇叭状分布，即模型存在异方差，通过对比2种加权残差方差模型与初始模型的评价指标可知（表6），考虑异方差结构明显提高了模型的预估精度（P<0.001），且均以指数函数效果最好，故选其作为混合模型的残差方差模型。

表6 各残差方差模型评价指标Tab.6 Evaluation indexes for each residual variance models

图1 冠幅-胸径模型的残差分布Fig.1 Residual distribution of the crown-diameter models

2.4 模型参数估计及评价

对各模型进行参数估计与评价（表7）。从表7 可以看出，未考虑随机效应时，广义模型的拟合精度较基础模型有所提高，较基础模型的R2提高了5.57%，RMSE和MAE分别降低了4.37%和2.99%；考虑随机效应时，混合模型的拟合精度较相应基础模型均有显著提升，基础混合模型较基础模型的R2提高了11.33%，RMSE和MAE分别降低了8.38%和8.28%，广义混合模型较广义模型的R2提高了8.29%，RMSE和MAE分别降低了6.44%和7.33%。

2.5 模型检验

对于基础模型和广义模型的检验，将检验数据直接代入到模型中，而混合效应模型的检验由二次抽样数据估计模型的随机参数，再将各参数代入模型中检验，多数研究表明选取样本中的4 株平均木作为检验样本对模型预测精度的提升最明显[26]，且随着样本株数的增加RMSE与MAE的降低速度缓慢[27]，故本研究从每个标准地中随机抽取4 株平均木，计算模型的随机效应参数值，将各参数代入各模型进行模型预测，检验结果见表8。由表8 可知，各模型的检验指标优度由大到小依次为广义混合模型、基础混合模型、广义模型和基础模型，说明增加林分变量与随机效应均能提高模型的预测精度。

表7 模型参数、方差估计值及评价指标Tab.7 Parameters，variance estimates and evaluation indexes for models

模型拟合与检验结果均表现为广义混合模型最好，故选择广义混合模型作为最优混合模型。将固定效应参数值分别代入其中，得到最优马尾松冠幅-胸径混合模型的表达式为：

表8 模型检验结果Tab.8 Results of models test

2.6 协变量对预测模型的影响

运用广义混合模型模拟分析各协变量（HCB、lvM）对冠幅预测的影响，以检验数据为基础，分段增加其中一个协变量，其余协变量采用总数据的平均值，结果见图2。由图2（a）可知，当控制lvM时，冠幅随HCB的增大呈逐渐增大的趋势，且随着DBH的增大各段HCB的冠幅大小差距越大，直至DBH等于11 cm后变化趋势趋于平缓；由图2（b）可知，当控制HCB时，冠幅随lvM的增大表现出减小的趋势，且随着密度组的增大各密度组单株林木冠幅的变化幅度逐渐减小。

图2 控制协变量冠幅-胸径模拟图Fig.2 Crown width-DBH curves under control covariates

3 结论与讨论

本研究以江西省兴国县飞播马尾松林为对象，筛选出幂函数为飞播马尾松林单木冠幅-胸径的最优基础模型，并采用逐步回归法将与冠幅相关性较高的HCB和lvM作为解释变量添加到基础模型中构建广义模型，对比分析了标准地效应对基础模型与广义模型拟合精度的影响。研究发现，添加lvM、HCB等解释变量和标准地随机效应的广义混合模型能够较好的拟合马尾松林单木冠幅，其R2为0.652 8，RMSE为0.624 6，MAE为0.492 8。

冠幅预测模型的预测精度决定了模型的实用性，本研究通过比较不同模型的RMSE和MAE发现，添加解释变量与考虑标准地的随机效应都能有效提高模型的预测精度，这与构建的樟子松[28]和椴树[29]混合模型所得出的结果相同或相似，但也有学者认为基础混合模型与增加解释变量的广义混合模型预测效果相差不大[30]。这可能是由于密度组在标准地水平上对冠幅产生的随机效应不同而导致研究结果的差异。

混合效应模型中，基础模型的选取对后期随机效应的添加与模型预测有较大影响，而由于树种生长特性及立地条件的不同，符合研究对象的模型也不尽相同，如对马尾松人工林[31]、侧柏[32]、杉木[33]等树种的冠幅基础模型均有所不同。本研究对飞播马尾松冠幅预测模型的拟合结果得出，林木在生长过程中冠幅与胸径呈正相关关系，林木树冠横向生长的速率随林木个体的生长呈缓慢减缓的趋势，其生长趋势较符合幂函数的曲线导向，并用于混合模型的构建与模型的预测。构建混合效应模型时容易受到异方差结构的影响[34]，而指数函数对模型异方差的消除有显著作用[25,35]，同时，相关研究也表明，当不同密度林分的立地条件差异不大时，广义混合模型能够较好的反映林分密度和标准地差异对冠幅的影响，对冠幅预测的精度较高，能较好预测不同密度林分冠幅的生长[8,28]。

在构建单木冠幅-胸径模型时选择影响冠幅生长的因子不尽相同，如符利勇等[18]在构建嵌套2 水平非线性混合冠幅模型中选择了枝下高、树高和标准地优势木平均高作为添加变量；Yang 等[36]在比较加拿大白云杉（Picea glauca）不同复杂程度模型预测效果研究中，选择了冠长、竞争和高径比作为协变量比较不同变量组合对模型预测效果的影响。本研究选择HCB和lvM作为协变量对模型精度有较好的提升，并与冠幅存在极显著相关性，通过控制HCB与lvM分段分析协变量对冠幅的影响发现HCB与冠幅呈正相关关系，该结论与Fu 等[37]对杉木的研究结果相类似，但有些研究得出HCB与冠幅大小呈负相关关系[7,25]，即冠幅减小时枝下高增大，这可能是因为本研究对象的立地条件较差，养分供应有限，限制了林木生长，导致林木生长时冠幅与树高的增大加剧了林木自然整枝的进程。lvM能体现林木在林分中的竞争状况，随着竞争的增大冠幅呈衰退现象[38]，因此，lvM与冠幅呈负相关关系，随着lvM的增大，林木在林分中的竞争压力也随之增大，这与Sharma 等[39]的研究结果一致；Sun 等[40]对樟子松冠幅模型的研究中也表明林分密度会抑制冠幅的生长；此外，Sharma 等[11]和高慧淋等[41]运用哑变量方法将林分密度作为虚拟变量建立了挪威云杉的单木冠幅预测模型与樟子松的树冠轮廓模型，对基础模型精度的提升有不小帮助，并得出林分密度的增加会导致单木树冠在水平及垂直方向上受到更大的压制效果。由此可见，由于林木冠幅的生长受林分密度、立地条件及人工管理措施等因素的影响，针对不同的林分，如何选择影响冠幅生长的因子对构建单木冠幅模型的精度有较大的影响。此外，不少研究表明冠幅会随林龄的增加而增大[38,42-43]，而本研究对象处于同一龄级，构建的冠幅预测模型是否适用于其他龄级需要进一步研究。