基于混合效应模型的塞罕坝华北落叶松人工林单木去皮胸径生长预测

洪李斌，卿蕴贤，田佳赫，康洁敏，卢伟

(河北农业大学 林学院，河北 保定 071000)

胸径作为林业中最基本的测树因子之一，是林分生长状况和生长潜力判定的重要依据，也是估算森林蓄积量、生物量、碳储量等指标的重要因子，为预估森林碳汇潜力，提高森林经营管理水平，实现我国碳达峰碳中和的林业发展目标提供重要参考[1-7]。塞罕坝机械林场拥有全球最大的人工林基地，是世界生态文明建设史上的典范，2018年塞罕坝林场实现了第一笔造林碳汇交易，标志着塞罕坝林场在碳汇项目上迈出实质性的步伐，有力诠释了“绿水青山就是金山银山”的理念，实现了森林生态效益和经济效益双赢。华北落叶松(Larixprincipis-rupprechtii)是塞罕坝地区主要造林树种，占林场有林地总面积的70.6%[8]，对塞罕坝地区水源涵养、防风固沙和维护该地区生态系统健康具有至关重要的作用[9]。因此，构建塞罕坝地区华北落叶松胸径预测模型对估算本地区森林碳汇储量，预测森林碳汇潜力具有重要意义。

针对单木胸径预测模型已有大量研究，于士涛基于塞罕坝林内160株华北落叶松解析木数据构建了华北落叶松估测胸径传统经验模型[10]；张冬燕在塞罕坝地区利用固定样地调查数据基于哑变量方法构建包含不同龄组的混合效应胸径生长模型[11]；周泽宇基于塞罕坝地区的标准地调查数据构建了华北落叶松人工林直径分布线性混合效应模型[12]。树皮生长量与木材生长量的速度并不相同，利用胸径尺测量数据构建胸径预测模型会引入树皮因子的误差，Finch研究表明如果不考虑树皮生长量，那么树木的带皮胸径生长量和材积的估算就会产生0～6%的误差[13]。传统经验模型多假设立地条件不变，而树木之间的竞争也会使同一年龄的林木胸径产生一定的差异，竞争能力与其对周围生存资源的占有及利用存在联系[14-15]。以往研究还表明树木生长和周围树木的相对空间位置和相互之间的距离显著相关，即林木的胸径生长受到邻木大小和密度制约[16-17]。混合效应模型是一种同时考虑时间变化和区域分布效应的有效模型，现已在单木以及林分等研究中取得了广泛的应用[18-23]。因此开展对华北落叶松单木去皮胸径预测模型的研究可以为林业调查提供技术支持，对华北落叶松的经营管理十分必要。

本研究以河北省塞罕坝机械林场的华北落叶松人工林为研究对象，于2018年7月在37块20 m×20 m的固定样地中钻取木芯数据，与传统胸径尺测量胸径相比，这种方法能够避免树皮因子的影响，可以得到确切的年龄与去皮胸径数据。通过引入样地竞争因子林分密度和单木竞争因子相对直径来构建塞罕坝地区华北落叶松单木去皮胸径预测的混合效应模型，以期能为塞罕坝地区华北落叶松的生长和经营提供参考依据。

1 数据与方法

1.1 研究区概况和数据来源

河北省塞罕坝机械林场(E 116°53′～117°39′，N 41°92′～42°36′)地处阴山山脉、大兴安岭余脉向西南延伸与燕山余脉的交合点，地形上分南部坝下和北部坝上2部分。该区域属寒温性大陆季风气候，海拔1 010～1 939.9 m。年平均气温-1.3 ℃，年均降水量460 mm，蒸发量1 339.2 mm，年均无霜期64 d，积雪时间长达7个月。土壤主为灰色森林土，该区域的主要乔木树种有华北落叶松、樟子松、白桦、云杉等[24]。

在塞罕坝机械林场华北落叶松纯林中设置37块固定样地，样地大小为20 m×20 m。于2018年7月在每个样地中选择长势良好的华北落叶松，使用口径5.5 mm的生长锥在距离地面1.3 m高度处，从2个不同的方向钻取2个样芯，并原则上保证每个样芯通过髓心以及保留树皮。将取好树芯粘到木质凹槽中，经晾干、打磨后进行交叉定年。测量并记录每个年龄的去皮胸径。建模数据概况见表1。

表1 建模数据统计表Table 1 Statistics table of modeling data

1.2 研究方法

1.2.1 变量选择及基础模型的构建 林木胸径估测模型通常是基于林龄、树高、冠幅等自变量与胸径关系进行构建。然而，在实际测量中，树木高度和冠幅的数据采集成本相对较高且精度不高。相比之下，林龄是一个较为容易获取的数据，尤其在人工林中。因此，本文首先选择了14个以年龄为唯一自变量的年龄-去皮胸径模型来估测去皮胸径，从中选择合适的模型进行混合效应模型的构建。具体模型形式见表2。

表2 14个基于林龄-去皮胸径的胸径预测模型Table 2 The 14 models based on the age-inside bark diameter at breast height for the prediction of DBH

1.2.2 混合效应模型的构建 利用方差膨胀系数(Variance inflation factor，VIF) 检验模型多重共线性，依据VIF值<10原则，挑选出3个变量，分别为林分密度(D)、林木直径与林分平均直径之比(RH)和立地指数(SI)，由于变量的增加使混合效应模型过复杂，故依据Pearson相关系数删去对方程贡献较小的立地指数(SI)。选取林分密度(D)以及林木直径与林分平均直径之比(RH)作为随机效应构建华北落叶松去皮胸径估测模型。

单水平线性混合效应模型：

DBH=a0+a0m0i+(a1+m1i)A+εik(εik)～N(0，Ri) (m0im1i)T～N(0，Dm)

(1)

DBH=a0+a0u0j+(a1+u1j)A+εjk(εjk)～N(0，Rj) (u0ju1j)T～N(0，Du)

(2)

嵌套两水平线性混合效应模型：

图书馆阅读推广平台后台设计包括系统管理、用户管理、阅读推广服务、资源链接、读者积分管理、日志管理及其他管理模块。具体的设计如图2所示。

DBH=a0+a0m0i+a0u0ij+(a1+m1i+u1ij)A+εijk(εijk)～N(0，Rij) (m0im1i)T～N(0，Dm) (u0iju1ij)T～N(0，Du)

(3)

式中：A为林木年龄，i为密度水平编号，j为相对直径编号，k为样木编号，m0i、m1i为密度水平随机效应参数，u0i、u1i为相对直径水平随机效应参数；εik、εjk、εijk为残差向量，Dm、Du随机参数的方差-协方差矩阵。在混合效应模型构建的过程中需要明确3个问题。

(1)模型参数的确定

(4)

(5)

(6)

(7)

(2)误差方差-协方差结构的确定。

林业中常用于描述方差-协方差结构表达式为：

(8)

本研究采用广义正定矩阵(General positive-definite matrix)来作为随机效应的方差-协方差结构。其结构如下：

(9)

以上计算过程均在R软件(版本4.0.3)中完成，华北落叶松去皮胸径估测基础模型由R中的“lm”和“nls”函数完成；混合效应模型其建模和预测有R包的“NLME”完成。

2 结果与分析

2.1 基础模型的拟合

在R软件中将建模数据代入到14个林龄-去皮胸径估测模型中，对各个模型进行拟合，拟合结果见表3。

表3 基础模型的建模结果Table 3 Modeling results of base models

由表3可知，14个基础模型中除了Model.3以外的都达到了0.8以上，表明各模型对去皮胸径的拟合效果较好。其中，Model.1、Model.4、Model.5、Model.7、Model.12、Model.13、Model.14决定系数均达到了0.85，说明相对于其他7个模型拟合效果更好。其中，Model.1的赤池信息准则(AIC)值低于其余模型，同时考虑到实际工作中适用性，本研究选择Model.1作为基础模型。

2.2 混合效应模型

考虑不同随机参数的组合，对基础模型进行拟合，根据表4可见，基于竞争效应的华北落叶松去皮胸径估测混合效应模型共有8种收敛的情况，MEM.2、MEM.4、MEM.6、MEM.7以及MEM.8模拟效果优于传统模型，其中MEM.8模拟精度最高。具体MEM.8参数估计见表5。

表4 去皮胸径估测线性混合模型收敛结果Table 4 Simulated results of estimation inside bark diameter at breast height (DBH) on linear mixed model

表5 MEM.8去皮胸径估测线性混合模型的参数估计Table 5 Parameter estimation of linear mixed model of MEM.8 inside bark diameter at breast height (DBH) estimation

2.3 传统模型和混合效应模型的比较

本研究通过对传统模型以及不同水平的混合效应模型进行比较分析，探寻能在塞罕坝地区更精确预测华北落叶松去皮胸径的模型形式,详见表6和图1。

表6 不同模型拟合统计量Table 6 Fitting statistics of different models

(a)

(b)图1 传统模型(a)与混合效应模型(b)的的残差图Figure 1 Residual plots of the traditional model (a) and the mixed effects model (b)

由表6可知，相对直径水平混合效应模型比密度水平混合效应模型更能体现竞争对林木生长的影响；与传统模型和单水平混合效应模型相比，嵌套两水平的混合模型拟合精度显著提高。根据图1可以看出，混合效应模型的残差离散程度较小，说明混合效应模型的拟合效果较好。

3 讨论

林木生长是一个复杂的过程，个体之的差异往往受到立地条件、光照、温度以及林木相互之间竞争等多方面的的影响，故仅通过年龄推断去皮胸径是不准确的。本研究在仅考虑年龄的情况下，其均方根误差(RMSE)达到2 cm左右，误差较大，这与周润凯等人的研究结果一致[26]。因此在构建去皮胸径预测模型时必须要考虑林木的生存环境，而刘泰瑞在关帝山的研究表明，华北落叶松受不同竞争压力后胸径生长量差异显著，由此得知林木胸径对竞争压力的响应往往超过其对立地条件、光照、温度等其他因素的响应[27-28]。考虑到所要构建模型能兼顾实用性、复杂性和精确性的要求，因此除年龄外，竞争是必须的建模因子。在引入竞争因子构建华北落叶松去皮胸径混合效应模型后，预测精度被大幅提高，尤其是嵌套两水平的混合效应模型，其均方根误差(RMSE)达到了0.29 cm，均方根误差降低了85%，本结论与王玉涛、杨惠滨的研究结果相一致[29-30]。混合效应模型的拟合精度优于传统经验模型，其原因主要是混合效应模型既可以反映总体的平均变化趋势，又可以提高数据方差、协方差等多种信息来反映个体之间的差异[31]。李春明、符利勇等人运用混合效应模型显著提高了单木胸径生长量的拟合精度[32-33]。

为进一步探讨竞争因子对华北落叶松去皮胸径-年龄关系的影响，分别建立密度水平、相对直径水平以及嵌套两水平的混合效应模型。结果表明，单水平混合模型中，相对直径水平拟合精度优于密度水平，而嵌套两水平的混合效应模型拟合效果达到了最优，明显提高了对去皮胸径的预测精度。说明在预测单木去皮胸径时，仅考虑林分尺度的竞争因子是不够的，更需要考虑单木之间竞争的差异。嵌套两水平混合效应模型拟合精度进一步提高，体现对林分尺度竞争因子分组后考虑单木竞争效应对去皮胸径生长影响更为准确，且嵌套两水平混合效应模型还可以分析不同随机效应对去皮胸径生长影响程度。本研究数据来源于塞罕坝机械林场，在实际应用中需要考虑区域性差异以及实际需要，选择不同层次进行模拟和预测。

本研究数据中缺乏大于42年的华北落叶松木芯数据，具有一定的局限性，后续应当对数据进行补全，完善模型。尽管在引入竞争因子后，模型达到了很高的精度，但桑杰通过建立气候影响下的青冈栎胸径生长预测混合效应模型，模拟不同气候条件下青冈栎的生长趋势，为华北落叶松胸径预测模型的研究提供了新的思路[34]。以后的研究需考虑将竞争指数作为参数直接用于构建胸径模型，再引入气候等因子进一步从机理上解释胸径的增长，探讨光照、温度等气候因子对预测模型的敏感性和贡献性。最后，通过研究华北落叶松去皮胸径随年龄变化的规律，可进一步研究华北落叶松带皮胸径与去皮胸径的关系，探寻树皮生长量随林龄增加的变化规律，达到改进华北落叶松胸径生长模型的目的。

4 结论