基于立地与密度效应的湖南栎类天然林平均木胸径生长模型

何 静，朱光玉，张学余，王忠诚，刘洪娜，屠维亚

（1.中南林业科技大学 林学院，湖南 长沙 410004；2.青羊湖国有林场，湖南 宁乡 410627；3.金洞管理区林业局，湖南 永州 426191）

天然林是我国森林资源的重要组成部分，是人类社会生存和发展的重要物质基础[1]。天然混交林与人工林相比，往往具有组成树种多、异龄等特点，这些特点决定了在构建天然混交林生长模型时，必须解决天然混交林年龄、立地和林分竞争的表达问题[2]。为了充分了解林分结构以及林分的发展并为林分的经营措施提供依据，构建林分生长与收获模型就显得十分重要。单木生长模型又是林木生长与收获模型的重要组成部分[3-4]。

单木生长模型是分析树木生长过程的有效手段，它既可以对目前的生长状况进行评价，也可以对未来的生长趋势进行预估，并为森林的经营提供依据[5-6]。林分胸径作为林分基本测树因子，可以用围尺简单测定，具有较高的准确性，因此可以用胸径的大小评价林分生长状况。目前，对树木生长模型的研究也有许多,一些学者不考虑其他因素的影响，直接通过树木理论生长方程对单木胸径、材积、年龄进行拟合，构建单木生长模型[7-8]。也有学者在研究中认为了解树木生长趋势时的一个主要问题是难以分别量化树木大小和年龄的影响，所以他们提出建立一种新的二维对数正态生长模型，此模型通过明确划分年龄和大小的影响，控制协变量来评估树木生长随时间的趋势[9]。由于树木的生长对年平均风速、湿度和降雨量等气象因子有一定的敏感性，坡度、坡位、海拔等立地因子以及林分密度，施肥等都对林分胸径的生长产生影响[10-14]。传统的理论生长方程缺少反映密度、立地类型等影响因子的变量，无法预测不同条件下林分的生长变化。近些年，由于混合效应在林业中有十分广泛的应用，在考虑影响因素对林分生长的影响时，许多学者选择了构建混合效应生长模型[15-18]，将因子引入到参数上，构建基于密度水平的林分优势木断面积生长模型、基于气候条件的单木胸径生长模型等等。与传统的回归方法相比，基于混合效应构建的模型拟合效果明显优于传统方法，即混合效应能明显提高模型的拟合效果和预估精度。

栎类天然林是我国天然林的主要林分类型，作为一种阔叶用材林,栎类在湖南林业生产上发挥了重要的作用，但目前湖南栎类天然林存在着林分质量不高、林分平均胸径较小等问题[19]。建立生长模型是森林经营的基础，密度和立地类型对林分的生长有较大的影响，林分密度越大，林分胸径生长越受限制；不同林分的生长对立地条件也有所要求。目前，已有许多学者对栎类的生长模型进行了研究，但结合密度和立地因子建立胸径生长模型的研究较少[20]。本研究以湖南栎类天然林为研究对象，基于湖南栎类样地数据采用混合效应构建含林分密度和立地类型的胸径生长模型，为不同密度和不同立地条件下栎类天然林的生长和可持续经营提供经验模型和参考。

1 研究区概况

研究区位于湖南省，地处长江中游，是云贵高原向汉江平原过渡的地区，三面环山，与多个省份交界，全省年日照时长1 500 h 左右，多年平均降水量在1 450 mm，富含各种资源。土壤主要以红壤和黄壤为主。湖南省植物种类十分丰富，有水杉Metasequoia glyptostroboides、银杉Cathaya argyrophyllaChun et Kuang、单性木兰Kmeria septentrionalis、珙桐Davidia involucrata、伯乐树Bretschneidera sinensis等国家重点保护植物。同时，湖南是我国栎类植物的近分布中心，栎类资源丰富，面积大，分布广，种类多[21]。

2 研究方法

2.1 数据来源

本研究数据来源于湖南省五盖山林场、八大公山自然保护区、平江县芦头林场等5 个地点设置的51 块栎类天然林固定样地，对样地内达到起测径阶（D1.3＞5 cm）的林木进行每木检尺。样地调查内容包括基本树木测量因子如树高、胸径等和土层厚度、坡度、坡位、海拔等立地因子。对每一块样地的数据进行预处理，求得每个样地的株数密度、郁闭度以及样地内林木的平均胸径、平均高，找出每个样地平均标准木，通过生长锥钻取标准木，得到林分的平均年龄。以湖南栎类天然林为研究对象，对林分平均木胸径生长的影响因子进行分析，并进行建模分析。

2.2 数据整理

1）林分平均胸径的计算[22]

式中：Dg为林分平均胸径（cm），n为林木株数，di为第i株林木的胸径（cm）。

2）林分调查因子统计

2.3 研究方法

2.3.1 平均木胸径生长基础模型筛选

2.3.2 林分调查因子分级及显著性因子筛选

1）林分调查因子分级

为了更方便的筛选显著性影响因子、划分立地类型和模型的建立，将海拔、坡度、坡向、坡位、土壤厚度和土壤类型等立地因子和株数密度进行分级处理，分级标准见下表：

2）数量化方法Ⅰ筛选因子：

数量化方法Ⅰ是指在自变量中包含定性因子的“回归”模型。表达式如下：

运用Forstat 统计之林软件进行数量化方法Ⅰ分析，对可能影响栎类平均木胸径生长的因子进行筛选。

2.3.3 混合效应模型构建

非线性混合效应模型（Nonlinear mixed effects models，简称NLMEMs），是依据回归函数依赖于效应参数非线性关系而建立的数学模型，在处理纵向数据、多水平数据及重复数据上有着突出的优势[23-24]。效应参数包括固定效应参数和随机效应参数。在构建混合效应模型时，首先要进行混合效应因子的选取：将筛选出的影响显著的立地因子进行分级组合，构成立地类型组合（ST）；再与其他影响显著的因子分别作为固定效应或随机效应加入到基础模型中，构建混合效应模型。

2.3.4 模型评价与检验

模型评价采用确定系数（R2）、平均绝对误差（MAE）及均方根误差（RMSE），3 个评价指标计算公式如下。其中，R2用于评价模型对建模样本的拟合效果，MAE 和RMSE 用于评价模型对检验样本的预测效果。

式中：yi为第i个样本的实测值，为样本的平均值，为第i个样本的预估值，n为所有样本数。

3 结果与分析

3.1 平均木胸径基础模型筛选

基于51 块样地的调查数据，分别以林分平均木胸径和平均林龄为因变量、自变量，构建林分平均木胸径生长的基础模型。通过R 语言，对表1中的理论生长模型进行拟合。由表4的拟合结果可知，4 个理论模型的拟合确定系数（R2）相差不大，都在0.3左右。其中拟合效果最佳的是Logistic模型，确定系数（R2）最大，达到0.338 97，平均绝对误差（MAE）和均方根误差（RMSE）最小，分别为2.606 19、3.278 92。所以最终选择的最优基础模型为Logistic 模型。具体模型表达式为：

表1 林分调查因子统计Table 1 Statistics of forest stand investigation factors

表2 平均木的4 种胸径-年龄曲线模型†Table 2 Four types of DBH-age curve models of average wood

表3 林分调查因子等级划分Table 3 The division of forest stand investigation factors

表4 4 种候选模型的参数拟合与精度评价Table 4 The parameter fitting and accuracy evaluation of four basic models

式中：D为林分平均木胸径，AGE 为林分平均木年龄，a、b、c为模型参数。

3.2 混合效应因子的构建

以林分的平均胸径为因变量，以调查因子为自变量，通过数量化方法Ⅰ，对影响胸径生长的因子进行筛选，显著性分析结果见表5。根据方差分析表中的“P＞F”值，对立地因子进行筛选，从而确认主要影响因子,由表可知海拔（ALT）、坡度（SLO）、坡向（SA）和株数密度（ND）对平均胸径的生长影响显著。将其中影响显著的立地因子（海拔、坡度、坡向）通过分级与组合来划分立地类型组合（ST）。经过组合可将51 个样本划分为34 个立地类型组合（ST）。

表5 林分调查因子的显著性检验†Table 5 Significance test of forest stand investigation factors

3.3 混合效应模型构建

在构建栎类天然林林分平均木胸径—年龄基础模型时发现，模型的精度较低，为了提高模型的精度，使模型更准确的预测林分胸径的生长规律，添加株数密度和立地类型作为混合效应。

3.3.1 添加株数密度作为固定效应

为了确定林分的密度对林分平均木胸径生长的影响，将株数密度作为固定效应加在参数上，分别构建的模型为：

将上述7 个模型用R 语言进行数据分析，得到结果如下：

由表中的评价指标可以得知，将株数密度作为固定效应加入到基础模型中后，模型的精度有所提高，MAE 和RMSE 都有所减小。当株数密度作为固定效应放在参数a 上时，确定系数R2由0.338 97提高到0.569 61；当固定效应放到参数c 上时，确定系数R2由0.338 97 提高到0.585 50；当固定效应放到参数a和c上时，确定系数R2由0.338 97提高到0.585 51；说明株数密度对林分胸径的生长有显著影响。

表6 固定效应模型参数估计Table 6 Parameter estimation of fixed effect model

3.3.2 添加立地类型组合作为随机效应

考虑不同立地条件下，林分胸径生长可能存在差异，所以通过数量化方法Ⅰ对立地因子进行筛选，通过筛选结果可知：海拔（ALT）、坡度（SLO）、坡向（SA）3 个立地因子对林分平均木胸径生长影响最大，将这3 种立地因子按照表2进行等级划分后排列组合可以形成200 种立地类型组合（ST），本研究数据中只用到了其中的34 种，将组合后的立地类型组合作为随机效应加入到模型中。

在基础模型（3.1）的基础上添加随机效应，模型表达式为：

在基础模型（3.2）的基础上添加随机效应，模型表达式为：

在基础模型（3.4）的基础上添加随机效应，模型表达式为：

在基础模型（3.6）的基础上添加随机效应，模型表达式为：

上式中Dj表示第j种立地类型组合下的林分平均木胸径。

根据表7～8 的结果可知：当加入立地类型组合（ST）后，模型精度相比于基础模型6 均有大幅提高，说明立地类型对林分平均木胸径的生长有较大的影响。其中模型14 在基础模型6 的基础上添加立地类型组合（ST）作为随机效应后，确定系数R2由0.338 97 提高到0.883 91，MAE 和RMSE 均有所下降，其中MAE 由2.606 19 下降到0.990 47，RMSE 由3.278 92 下降到1.374 10；模型15 在将株数密度作为固定效应放在参数a上而构建的模型7 的基础上添加ST 为随机效应，确定系数R2由0.338 97 提高到0.930 16，MAE 和均RMSE 分别下降到0.836 70、1.065 80；模型16 在将株数密度作为固定效应放在参数c上而构建的模型9 的基础上添加ST 为随机效应，确定系数R2由0.338 97 提高到0.585 50。

表7 固定效应模型评价Table 7 Evaluation of fixed effect model

表8 随机效应模型参数估计Table 8 Parameter estimation of random effect model

根据这几个评价指标我们可知：将株数密度作为固定效应放在参数a上，再添加立地类型组合作为随机效应，最终得到的混合效应模型15 为林分平均木胸径生长的最优模型。由此也说明株数密度和立地类型对林分平均木胸径的生长影响显著。

3.3 模型残差

为了更好的反应混合效应模型的拟合效果，绘制每个模型的预测值-残差分布图如下：

由图1可知，基础胸径模型（模型6）预测值与残差的离散程度大于混合效应模型的离散程度，混合效应模型的残差相比于基础模型的残差基本在X 轴两侧，且分布更加均匀、更加集中，这说明胸径混合效应模型的拟合效果要优于基础模型的拟合效果。其中，模型15 的残差分布最为集中，最为均匀，说明模型15 的拟合效果最佳。

图1 模型预测值与残值图Fig.1 Prediction value and residual value diagram of model

续图1Continuation of Fig.1

表9 随机效应模型评价Table 9 Evaluation of random effect model

4 结论与讨论

以湖南栎类天然林51 块样地数据为基础，以林分平均木胸径为研究对象，采用数量化方法Ⅰ筛选出株数密度和海拔、坡度、坡向这3 个立地因子对林分平均木胸径生长影响显著。胸径的生长既受林分密度的影响，又与立地条件有关。林分的密度与胸径呈现负相关，密度越大，林分平均木胸径越小，反之亦然，与张海东等[25]的研究一致。

利用R 语言从4 个理论生长方程中筛选最优基础模型，在拟合平均胸径最优基础模型时，根据评价指标（确定系数R2、平均绝对误差MAE及均方根误差RMSE），从4 个基础模型中选定Logistic 为林分平均木胸径生长的最优基础模型，确定系数R2为0.338 97；MAE 和RMSE 分别为2.606 19、3.278 92。

在最优基础模型的基础上将株数密度和海拔、坡度、坡向组合成的立地类型分别作为固定效应和随机效应放入基础模型中，构建混合效应模型。最终确定最优生长模型为Dj=(aj+a1*ND)/(1+bj*exp(-cj*AGE))+ε，确定系数R2为0.930 16，MAE 和RMSE 分别为0.836 70、1.065 80。相比于基础模型，拟合精度有显著提高。说明在建立林木生长模型时，混合效应模型的拟合效果明显高于传统模型[26-27]，在基础模型中添加了密度和立地因子可以更加准确的预估林木在不同密度和立地条件下的生长状况。通过一系列的分析，确定湖南栎类生长的最佳条件，为提高栎类平均胸径的生长提供依据和参考。

由于本研究只采取了湖南省内的栎类样地数据，所以模型的适用范围有限。另外本研究只考虑了立地类型和竞争对平均胸径生长的影响，林分平均木胸径可能还受其他因素的影响，如气候因子、林分空间结构、林分类型等。所以，还可以进一步探讨其他因素的影响，构建含气候因子、林分类型的林分平均木生长模型，进一步优化模型，更好的预估林分平均木胸径的生长。