基于气候因子的杉木人工林削度方程构建

张中伟，周根苗，易 烜，齐战涛，吕 勇

（1. 中南林业科技大学 林学院，湖南 长沙 410004；2. 湖南省农林工业勘察设计研究总院，湖南 长沙 410007；3. 湖南省青羊湖国有林场，湖南 宁乡 410600；4. 湘潭市自然资源和规划局，湖南 湘潭 411100）

削度是描述树干直径沿其树干向上随干径位置的升高而逐渐减小变化程度的指标，削度方程是削度的数学表达式[1]，是以树干上不同部位的树干直径为因变量，全树高、树干直径及不同部位时距地面高等变量为自变量的数学函数，根据研究方法的不同，削度方程又可分为简单削度方程、分段削度方程和可变指数削度方程[2]。同时，削度方程能为森林经营者提供树干任意部位直径的高度、全树干材积以及从地面到任意高度的商品材材积[3]，因此构建树干削度方程已成为编制材种出材率表的首选方法和基础工作。

树干的削度受树种、年龄、林分密度、气候因素等多种因素的影响[4-6]。目前，大量学者研究了林分密度因子[7]和林地因子[8]对树干削度的影响，而对气候因子影响树干削度的研究较少。但树木的生长对温度、降水量等气候因子非常敏感[9]，因此，近年来有学者开始逐渐研究气候因子对树干削度的影响。Gordon 等[10]在研究气候因子对加拿大樟子松Pinus sylvestnisvar.mongolica树干削度的影响时，发现年平均降水量和无霜期均能显著影响树干削度的变化，且对树干上部削度的影响更大；Liu 等[11]在研究兴安落叶松Larix olgensis的树干削度时，考虑气候因子对树干干形的影响，将气候因子添加到削度方程中，建立了含气候因子的树干削度方程，使模型的精度有显著的提升。然而，在以往的杉木Cunninghamia lanceolata树干削度方程中，缺少代表气候因子的变量，故无法反映树干削度随气候因子变化而变化的过程。因此，构建含气候因子的杉木树干削度方程具有非常重要的意义。

随着数学建模方法的发展，在构建削度方程模型时，混合效应模型方法成为了更多学者的选择，因为传统的回归分析方法即最小二乘法是假定数据间相互独立且非异质性[12-13]，因此无法体现出不同因子之间的差异对树干生长的影响，从而使模型的精度不够高，但混合效应模型方法却能很好地解决上述问题[14-16]。如张森森[8]、聂璐毅等[17]使用了混合效应模型方法分别构建了江西省杉木和长白落叶松的树干削度方程，研究结果表明，使用混合效应模型方法所构建的模型精度均优于传统方法所构建的模型，且误差均有明显的降低。

杉木是我国主要的用材树种，也是我国南方人工造林的主要树种，其具有分布区域广、经济效益高等特点，广泛应用于建筑、家具等领域。本研究以湖南省永州市杉木人工林为研究对象，使用混合效应模型方法构建了含气候因子的杉木树干削度方程，研究气候因子对杉木树干削度的影响，可为森林经营者在气候影响的条件下采取的杉木经营措施提供理论依据。

1 材料与方法

1.1 数据来源

本研究的数据来源于湖南省永州市东安县、道县和零陵区的杉木人工林固定样地及临时样地（108°47′～114°15′E，24°38′～30°08′N）。采集地区的气候为大陆性亚热带季风湿润气候，雨量充沛，水热充足，年均温度为18 ℃，年均降水量为1 500 mm，土壤肥沃，主要树种有：杉木、马尾松Pinus massoniana、木荷Schima superba、青冈Cyclobalanopsis glauca等。

本次一共调查了75 块样地的解析木数据，样地大小均为20 m×20 m，其中东安县有20 块样地，道县有13 块样地，零陵区有42 块样地。每块样地选择一株平均木伐倒，按照树干解析的标准对其进行树干解析，测量其树高（H）、胸径（D）、树干上各部位的直径（d）和该干径位置距地面的高（h），同时记录对应样木的经纬度坐标和海拔（HB），共计550 组数据，其因子基本情况详见表1。

表1 各因子基本情况†Table 1 Basic information table of factors

气候数据是从Wang 等[18]编写的可获取亚洲区域气候数据的ClimateAP（2019）软件中[19-21]获得，其中包括年积温（DD5）、夏季平均最高温（Txmax）、最冷月平均温度（TMmin）等共11 类气候因子。

为保证模型构建的合理性，本研究先使用三倍标准差法对异常数据进行剔除，剔除后剩476组数据，再将这476 组数据按照2∶1 的比例分为建模数据317 组与检验数据159 组，用以杉木人工林树干削度方程的构建与检验。

1.2 研究方法

1.2.1 基础模型的筛选

在三类削度方程中，可变指数方程有着更好的拟合效果[8]，能更好地用来描述树干削度的变化[2]。参考前人的研究，本研究选择了4 个描述杉木树干削度较好的可变指数削度方程作为备选的基础模型[22-25]，各方程的表达式如下：

式中：h为不同部位时距地面高；d为不同部位的树干直径；D为树干胸径；H为总树高；b1、b2、b3、b4、b5、b6、b7、b8、b9为参数；T=h/H，p=1.3/H，

为了后续表达方便，将模型（1）称为M1，模型（2）称为M2，模型（3）称为M3，模型（4）称为M4。

1.2.2 自变量的筛选

本研究使用Forstat 中的数量化方法Ⅰ对气候因子进行显著筛选，以“P＞F”为标准，剔除P＞0.05 的因子，从而保留影响显著的因子。参考前人的研究，林分密度因子对树干削度的影响非常显著[7-8]，所以本文直接使用林分密度因子作为显著影响因子。

1.2.3 非线性混合效应模型

依据回归函数同时依赖于固定效应参数和随机效应参数的回归关系而建立的模型称为混合效应模型，其一般形式[26]如下：

式中：yi与xi分别为第i个样地的因变量向量和自变量向量；εi为误差项；β与µi分别为固定效应参数向量和随机参数向量[27]。

1.2.4 模型的评价指标

本研究选用调整决定系数Ra2、均方根误差RMSE、平均绝对误差MAE、赤池信息准则AIC与贝叶斯信息准则BIC 对构建的模型进行评价，表达式如下：

式中：yi为第i样本的实测值；是第i样本的预估值；为平均实测值；p为模型中参数的个数；n为样本数；l为模型极大似然函数值。

2 结果与分析

2.1 基础模型的筛选

利用R 软件的nls 板块将4 个备选基础模型进行拟合，以Ra2、MAE 和RMSE 等评价指标对各基础模型进行评价，结果详见表2。

表2 候选基础模型评价Table 2 Candidate basic models

通过拟合，对比Ra2、MAE、RMSE 这3 个指标，可以得出，在备选的4 个模型中，模型M3 的各项指标均为最优，其Ra2最大为0.959 5，RMSE 最小为0.894 7，MAE 最小为0.658 9，很多学者的研究也得出了一样的结论[28-29]，说明Kozak（2004）模型的普适性较好，在可变指数模型中很有代表性[30]。因此将模型M3 作为后续研究的基础模型。

2.2 因子显著性分析

从获取的11 个影响树干削度的气候因子中，根据取值范围的不同，参考齐战涛等[31]的划分标准，对气候因子进行等级划分，结果详见表3。

表3 气候因子等级划分Table 3 The division of climatic factors grades

使用Forstat 中的数量化方法Ⅰ对各已分好级的气候因子进行显著性分析，并以“P＞F”为标准，剔除P＞0.05 的因子。由于削度是描述树干直径沿其树干向上随干径位置的升高而逐渐减小变化程度的指标，为考虑距地面高（h）、树干胸径（D）不同部位的树干直径（d）的影响，在筛选时，选取不同部位时距地面高（h）、树干胸径（D）作为辅助变量进行分析。各因子分析结果详见表4。

由表4 可得，生长积温（DD5）、夏季平均温度（Txave）、夏季平均最高温（Txmax）为影响树干上各部位的直径（d）变化的显著因子。由于Txave=（Txmax+Txmin）/2，因此夏季平均温度（Txave）与夏季平均最高温（Txmax）之间可能存在自相关、共线性问题，为避免这些问题对模型的构建造成影响，因此本次研究暂时不考虑夏季平均温度，仅选取生长积温、夏季平均最高温为显著影响树干上各部位的直径（d）变化的气候因子。其中生长积温的第一个等级为3 200 ≤DD5＜3 400；夏季平均最高温的第一个等级为25 ≤Txmax＜26。

2.3 含林分密度因子的非线性混合模型构建

针对选中的基础模型M3，使用非线性混合效应模型方法构建含林分密度因子的杉木树干削度方程，利用R 语言的nlme 包，将林分密度（DM）作为随机效应分别添加至参数b1、b2、b3、b4、b5、b6、b7、b8、b9上进行拟合，剔除不收敛的组合形式，根据AIC、BIC 对剩下的组合形式进行评选，最后筛选出最好的组合形式，结果详见表5。

表5 林分密度因子随机效应精度评价结果Table 5 Evaluation results of the random effects of stand density factors

由表5 可得，当林分密度（DM）作为随机效应添加在参数b8上时模型的AIC、BIC 值最低，即效果最佳。为后续表达方便，将构建的含林分密度因子的削度方程模型简称为M11。

因此，构建的含林分密度因子的混合效应模型表达式如下，简称M11：

式中：d、D、H、h、X、T分别为不同部位的树干直径、树干胸径、总树高、不同部位时距地面高、X、T值；a8i为DM的随机效应参数；参数b1、b2、b3、b4、b5、b6、b7、b8、b9为固定参数。

2.4 含气候因子和林分密度因子的非线性混合模型构建

为比较气候因子和林分密度因子对杉木树干削度的影响，针对选中的基础模型M3，使用非线性混合效应模型方法构建含双因子的杉木树干削度方程。

运用R 语言的nlme 包，分别以筛选所得的对树干削度有显著影响的2 个气候因子夏季平均最高温、生长积温和林分密度及其组合形式作为随机效应，将其分别添加至参数b1、b2、b3、b4、b5、b6、b7、b8、b9及其组合形式上进行拟合，剔除不收敛的组合形式，根据AIC、BIC 对剩下的组合形式进行评选，最后筛选出最好的组合形式，结果详见表6。

表6 气候因子和林分密度因子随机效应精度评价结果Table 6 Evaluation results of the random effects of climatic factors and stand density factors

从表6 可得，将Txmax添加到参数b5、DM添加到参数b7、DD5添加到参数b8上时其AIC、BIC值最低，即效果最好。因此，构建的含双因子的混合效应模型表达式如下，简称M12：

式中：d、D、H、h、X、T分别为不同部位的树干直径、树干胸径、总树高、不同部位时距地面高、X、T值；a5i、a7i、a8i分别为Txmax、DM与DD5的随机效应参数；参数b1、b2、b3、b4、b5、b6、b7、b8、b9为固定参数。

2.5 模型模拟分析

为分析模型M3、M11、M12 之间拟合表现的优劣，利用R 软件对模型M3、M11、M12 进行模拟，各拟合结果详见表7。

表7 模型的参数拟合与精度评价Table 7 Parameters fitting and precision evaluation of models

由表7 可得，模型M12 在AIC、BIC 上的拟合结果均低于模型M3 和M11，说明加入双因子的模型M12 的拟合效果要优于基础模型M3 和只加入林分密度因子的M11，拟合效果提升明显。从Ra2、RMSE、MAE 上的表现来看，模型M12无论是在建模样本还是检验样本上的拟合效果均优于模型M3 和M11。其建模样本调整后决定系数Ra2的排序为M12 ＞M11 ＞M3，其中M12 相较于M3 的增幅约为2.3%，相较于M11 的增幅约为1.1%；均方根误差RMSE 的排序为M12 ＜M11 ＜M3，其中M12 相较于M3 的降幅约为32.56%，相较于M11 的降幅约为20.31%；平均误差MAE 的排序为M12 ＜M11 ＜M3，其中M12相较于M3 的降幅约为34.56%，相较于M11 的降幅约为21.12%。

为了检验基础模型M3 和混合模型M12 是否存在异方差效应，建立了模型M12 和M3 的残差分布（图1 ～2）。

图1 模型M12 残差分布Fig. 1 Residual distribution diagram of model M12

图2 模型M3 残差分布Fig. 2 Residual distribution diagram of model M3

由图1 ～2 可得， 无论是模型M3 还是M12，它们的残差都较集中且随机分布在横坐标轴两侧，说明模型M3 和M12 都没有明显的异方差效应，故不用考虑。

3 结论与讨论

3.1 讨 论

研究杉木人工林树干削度方程对于杉木人工林的生长与收获的预测和森林经营方案的编制具有重要意义。本研究使用混合效应模型方法构建了含双因子和单含林分密度因子的杉木树干削度方程，以对比分析气候因子对树干削度的影响。姜立春等[32]使用了同样的方法研究了样地效应和树木效应对落叶松树干削度的影响，虽然研究的对象不同，但模型的精度都有显著提升；李春明等[33]使用混合效应模型方法构建了红松、冷杉和云杉的树干削度方程，同时将其与传统拟合方法所构建的树干削度方程进行比较，最终得出无论是哪种树种，混合模型都要优于传统拟合方法所构建的模型，因此构建削度方程时使用混合效应模型方法非常合理。

最终添至模型中的气候因子为夏季平均最高温和生长积温，相较于基础模型，加入气候因子的混合模型的精度有明显提升，说明气候因子对杉木树干削度有显著影响。Liu 等[11]在研究气候因子对兴安落叶松树干削度的影响时，得出影响树干削度显著的气候因子为年平均温度和年平均降水量，造成这种不同的原因可能是不同树种对不同气候因子的敏感程度不同；Gordon 等[10]在研究加拿大松树干削度时，引入气候因子，发现气候对树干上部的削度影响最大，认为削度和气候之间的联系是通过树冠来实现的，而树冠又是影响削度变化的一个非常重要的因子[5,34]，因此气候因子对树干削度的影响非常显著。

为减少模型的参数数量，本研究只考虑了气候因子和林分密度因子，没有将样地因子考虑在内。在后续对树干削度的研究中，可以将样地因子也一并考虑，加入到削度方程中，或许可以使所构建模型的拟合精度有进一步的提升。

3.2 结 论

本研究以湖南省杉木人工林的树干削度为研究对象，基于永州市75 块杉木人工林样地的实测数据，选取4 个备选基础模型，通过拟合，确定Kozak（2004）模型为后续研究的基础模型（Ra2=0.959 5，RMSE=0.894 7，MAE=0.658 9）。为研究气候因子对树干削度的影响，通过数量化方法Ⅰ对气候因子进行显著性分析，结果显示，夏季平均最高温和生长积温为影响杉木树干削度显著的气候因子。使用非线性混合效应模型方法构建了含林分密度因子和含双因子的杉木人工林树干削度方程，确定了最优的随机效应参数构造形式。从Ra2、RMSE、MAE 上的表现来说，含林分密度因子的模型M11 和含双因子的模型M12 相较于基础模型M3 均有明显提升，其中模型M12的Ra2最大，为0.981 6，相较于M3 的增幅为2.3%，相较于M11 的增幅为1.1%；RMSE 最小为0.603 4，相较于M3 的降幅为32.56%，相较于M11 的降幅为20.31%；MAE 最小为0.431 2，相较于M3 的降幅为34.56%，相较于M11 的降幅为21.12%。上述研究结果表明，气候因子对杉木树干削度有显著影响，为今后在研究树干削度时考虑气候因子对树干削度的影响提供了理论依据，其既提高了模型的拟合精度，同时也有利于杉木人工林森林经营管理措施的制定和林业数表的编制。