基于空间结构参数的大兴安岭天然落叶松单木直径生长模型

吕沅杭，伊利启，王儒林，刘兆刚,2，董灵波,2＊

(1. 东北林业大学林学院，黑龙江 哈尔滨 150040；2. 东北林业大学林学院 森林生态系统可持续经营教育部重点实验室，黑龙江 哈尔滨 150040)

单木模型通常以树木个体的生长信息为基础，从林木生长的竞争机制出发，模拟林分内不同大小树木的生长过程[1-2]。与林分和径阶水平模型相比，单木模型不仅能够提供较为详细的林分结构特征及动态变化信息，而且还是研究林木生长量、收获量以及生物量的重要基础，因此单木生长模型是定量研究林木自身生长过程及其对各种干扰过程响应机制的有效手段。根据模型中是否包含距离信息，单木生长可分为与距离有关和与距离无关模型[3]，其中由于前者需要详细的距离信息，这极大程度上限制该类模型在林业生产中的应用，因此与距离无关的单木生长模型一直是国内外林业研究的重要内容。

传统研究认为林木直径生长主要受树木自身大小、竞争因子和立地条件的影响[1,4-5]。但近些年来，国内外学者逐渐关注到物种多样性和林分空间结构特征对树木生长的影响。在物种多样性方面，Vannoppen 等[6]发现林分物种多样性对北欧地区山毛榉（Fagus sylvaticaLinn.）和夏栎（Quercus roburLinn.）直径生长量的贡献量分别为2% 和21%；Wang 等[7]研究表明中国东北地区近熟林和过熟林内分别有66.7% 和64.3% 的样地内生物量和木本植物多样性间存在显著相关性；Liang 等[8]基于全球44 个国家约77 万个固定样地数据的研究结果表明：生物多样性与生产力间存在显著的正相关；谭凌照等[9]发现：在蛟河阔叶红松（Pinus koraiensisSieb. et Zucc.）林中物种多样性均匀度指数与生产力显著负相关。在林分空间结构特征方面，汤孟平等[10]研究表明：较高聚集程度、年龄隔离度以及较低的最邻近竹株数和竞争指数均有利于增加毛竹（Phyllostachys edulis(Carr.) H.de Lehaie）林单位面积生物量；吕延杰等[11]研究表明：当红皮云杉（Picea koraiensisNakai）和冷杉（Abies fabri(Mast.)Craib）处于随机分布时胸径生长量最大，且其随着混交度的增加而增加。总的来看，物种多样性与林分空间结构特征对单木直径生长的作用已得到广泛认同，但现阶段还鲜有将其整合到单木生长模型中的研究。

现阶段单木生长模型数据源主要包括2 类，即固定样地复测数据[1-2,4-5]和单木解析数据[12]。从统计学角度来看，2 种数据均存在明显的嵌套效应（即区域→样地→单木）、空间自相关性以及时间自相关性等问题，这些均严重违背传统统计学模型中对数据正态性、独立性和方差齐性的要求[13]，致使所建模型不能很好满足林业生产需要。与传统回归模型不同，混合模型只要求误差项为正态分布，而对模型误差项之间是否相互独立和是否等方差不做要求，所以混合模型能够处理更多类型的数据。目前，国内外学者采用混合模型对单木生长模型也进行了一些探索，如姜立春等[12]采用解析木数据建立了兴安落叶松（Larix gmelinii(Rupr.)Kuzen）单木直径和树高模型，但并未考虑林分变量、物种多样性和结构参数对生长的影响；彭娓等[4]基于固定样地复测数据建立了长白落叶松（Larix olgensisHenry）人工林的单木直径生长模型；王蒙等[5]则进一步考虑了抚育间伐对落叶松人工林单木直径生长的影响；祖笑锋等[14]采用解析木数据建立了美国黄松（Pinus ponderosaDougl. ex Laws.）林分优势高生长模型。这些研究均表明，混合模型不仅能够大幅度提升模型的拟合效果，同时还能有效反应区域环境、林分、经营等变量对林木生长过程的潜在影响。

由于长期以来过度的人为采伐和严重的火灾干扰，当地顶级群落兴安落叶松林已逐渐退化为以兴安落叶松和白桦（Betula platyphyllaSuk.）为主要优势树种的天然次生林，其结构、功能和生产力均与原始顶级群落存在显著差异。因此，通过建立该地区次生林中兴安落叶松单木直径生长模型，对揭示林木生长与林分因子、单木因子、物种多样性以及空间结构特征间的关系具有重要作用。为此，本研究以兴安落叶松天然次生林为研究对象，利用混合模型技术在传统单木生长模型基础上逐渐引入林分因子、单木因子、物种多样性以及林分空间结构参数，进而建立该地区兴安落叶松单木直径生长模型，为该地区次生林的有效恢复提供理论依据和技术支撑。

1 材料与方法

1.1 研究区概况

盘古林场（52°16′38″～52°47′04″ N、123°20′02″～124°21′40″ E）属于大兴安岭地区塔河林业局管辖，总经营面积15.2 万hm2。该区属于大兴安岭中部伊勒呼里山北麓，年均气温-3℃，极端最低温和最高温分别为-53℃和36℃，平均生长季约100 d；地带性土壤为棕色针叶林土，其中土壤有机质、全氮和全磷平均含量分别为330.48 g·kg-1、7.01 mg·kg-1和81.24 g·kg-1[15]；地带性植被为兴安落叶松，但现阶段主要为以兴安落叶松和白桦为主要树种的天然次生林，其伴生树种主要包括樟子松（Pinus sylvestrisvar.mongolicaLitv.）、红皮云杉、山杨（Populus davidianaDode）等。

1.2 数据收集

在全面踏查基础上，于2011 年7—8 月份在盘古林场选择不同年龄、不同密度、不同立地条件的林分设置固定样地18 块。各样地面积20 m × 30 m，采用相邻格网调查法，实测样地内直径 ≥ 5 cm 树木的树种、胸径、树高、枝下高、冠幅、坐标以及状态等指标。每块样地按等断面积径级标准木法，将林木大小划分为5 个等级，将各等级林木的平均胸径和树高作为选择解析木的标准，并进一步参考该样地树种组成在样地附近选伐样木。将解析木树干按1 m 区分段进行区分，分别在胸高处（1.3 m）和各区分处截取3 cm 厚圆盘，并在非工作面上标出南、北方向。在室内将圆盘的工作面刨光，并通过髓心分别画出东西、南北2 条直线，进而采用200 dpi 分辨率扫描成jpg 图像。将图像导入WinDENDRO 软件中，根据图像中早晚材色泽差异识别出各年轮具体位置，且要求各方向中识别出的年轮数必须严格相等，以各方向各年轮宽度的平均值作为该年龄下的平均直径。从18 块样地中共收集兴安落叶松解析样木50 株（其余40 株均为其他树种）。林分和单木统计数据如表1 所示。

表1 大兴安岭天然落叶松林分和单木调查因子统计Table 1 Descriptive statistics of stand- and tree-level of natural Lari xgmelinii stand in Daxing′anling Mountains

1.3 研究方法

1.3.1 基础模型 由于树木直径生长通常具有典型的非线性特征，因此选用一些常用的非线性曲线作为备选模型， 包括Mitscherlich 方程、Logistic方程、Weibull 方程、Gompertz 方程、Schumacher方程、Korf 方程、Hosfled 方程以及Richards 方程，各方程形式详见孟宪宇[3]，此处不再赘述。根据模型参数估计值、拟合效果以及残差分布图等特征，确定适于描述天然落叶松单木直径生长过程的基础模型。考虑到树木所处环境会对林木直径生长产生影响，因此除将常规林分因子和单木因子加入到模型中外，本研究还尝试引入林分物种多样性、林分空间结构参数指标，对基础模型的参数进行再参数化处理。林分因子主要包括林分平均胸径（Dg）、平均树高（Hg）、株数密度（N）、郁闭度（Dc）等，单木因子包括胸径（DBH）、树高（HT）、冠幅（CW）、冠长（CL）等，物种多样性指标包括Margalef 丰富度指数（Ma）、Simpson多样性指数（Sp）、Shannon-Wiener 指数（Sw）、Pielow 均匀度指数（Pw），林分空间结构指标包括混交度（M）、大小比（U）、角尺度（W）和Heygi 竞争指数（H）。此外，还进一步考虑上述各变量间的组合及其变换形式（如对数、平方等）。各物种多样性指数和结构参数具体定义详见谭凌照等[9]和吕延杰等[11]，此处不再赘述。

1.3.2 混合模型 采用混合模型建立兴安落叶松单木直径生长模型，一是为通过引入随机因素而提高模拟精度，二是为了有效处理解析木数据中存在的异方差和时间自相关问题。根据Pinheiro 等定义[13]，非线性混合效应模型的基本形式为[13]：

式中：yi j为第i棵树的第j次观测值；f()为含有参数向量 φij和协变量vi j的 函数；Aij为已知固定效应的设计矩阵； β为未知的固定参数向量；Bij为已知随机效应的设计矩阵；bi为带有方差协方差结构矩阵的未知随机效应参数向量，其服从期望值为0、方差为D的正态分布；D为随机效应的协方差矩阵；ei j为 服从正态分布（期望为0、方差为 σ2Ri）的误差项；Ri为第i株树木的方差协方差矩阵。因各样地解析木数量相对较少，因此本研究仅考虑单木水平的随机因子。

对已确定的基础模型，构建混合模型时还需确定随机参数、方差协方差结构、时间序列结构以及异方差校正函数共4 类参数[12-13]。对随机参数，通常认为所有参数都有成为随机效应参数的可能，因此可对不同随机参数的组合模型进行拟合，通过多项统计指标来评价模型拟合效果。对混合模型中的方差协方差结构，本研究检验林业上3 种常用结构，即复合对称矩阵CS、对角矩阵UN(1) 和广义正定矩阵UN；为进一步提高模型拟合效果，采用复合对称结构CS、一阶自回归结构AR(1)、一阶自回归与滑动平均结构ARMA(1,1)以及高斯结构corGaussian 来描述解析木数据中的时间序列相关性。采用指数函数varExp、幂函数varPower 和常数加幂函数varConstPower 来校正解析木数据中的异方差现象。各模型拟合优度通过赤池信息准则（AIC）、贝叶斯准则（BIC）和负对数似然值（-LL）来确定。此外，对不同参数个数的模型，还需进行似然比检验（LRT）以最终确定最优模型。

1.3.3 模型检验 混合模型中包含固定参数β和随机参数bi两部分，其中固定效应参数与普通模型相同，而随机参数则需要根据模型拟合结果和检验数据重新计算，其具体公式为[13]：

各基础模型和混合效应模型的精度评价和检验均采用五折交叉检验法进行，即建模过程需采用全部样本数据，而检验过程则将数据随机分为5 份，每份单独进行评价，以其平均值作为最终结果。各模型采用调整系数（Ra2）、平均绝对误差（Bias）和均方根误差（RMSE）3 项指标进行评价和精度检验[4-5]，其中R2越大越好，Bias和RMSE均越接近0 越好。

1.3.4 数据处理 分别采用S-Plus 软件中的NLS 模块拟合各基础模型，采用NLME 模块拟合各混合效应模型；其余数据统计和图表绘制均采用EXCEL 完成。

表2 基础模型参数估计及拟合优度Table 2 The parameter estimation and goodness-of-fit for the base models

2 结果与分析

2.1 基础模型筛选

在8 种 候 选 模 型 中， 仅Mitscherlich、Schumacher、Logistic 和Gompertz 方程能够收敛，且各模型参数均能通过α= 0.05 水平的显著性检验（表2）。虽然各项指标均显示Gompertz 方程的拟合效果最好，但方差分析结果表明Mitscherlich方程和Gompertz 方程间的差异不显著（P= 0.09），且在后续再参数化过程中也较难收敛，因此在综合评估各模型拟合优度、模型参数个数以及残差分布图后，选择Mitscherlich 方程为大兴安岭地区兴安落叶松单木直径生长的最优基础模型。

对不同单木进行Mitscherlich 式拟合，结果表明模型收敛率达到82.86%，可以有效反映所有建模样本的整体特征。Mitscherlich 模型的参数a和b与林分变量、单木变量、物种多样性、结构参数的相关性分析（图1）。可以看出，参数a与林分变量中的平均树高（r= 0.528 9）、单木变量中的胸径（r= 0.787 1）、树高（r= 0.796 2）和冠幅（r=0.526 4）均显著正相关，而参数b则仅与林分变量中的平均树高（r= -0.504 2）显著负相关。为进一步揭示物种多样性和林分结构参数对直径生长的影响，结合相关性分析、变量共线性以及参数可获得性等问题，综合确定参数a的候选变量包括林分平均树高、单木胸径以及混交度和角尺度，而参数b的候选变量包括林分平均胸径、Pielow 均匀度指数（与Sw指数的r= 0.770 4）以及Heygi 竞争指数，其中平均树高和单木胸径为参数a的固定选项，平均胸径为参数b的固定选项。

图1 Mitscherlich 模型参数a 和b 与林分变量、单木变量、物种多样性以及林分结构参数的相关性Fig. 1 Relations between the parameters of Mitscherlich and the indices of stand-, tree-level, diversity and spatial structures

因候选变量个数较多，因此本研究采用相乘思想构造函数对Mitscherlich 模型中参数a和b进行预测，具体形式如下：

式中：a0,a1, ···,aN为Mitscherlich 模型参数a的待估系数，b0,b1,···,bN为Mitscherlich 模型参数b的待估系数，N为候选变量个数，X1X2···XN分别为各候选变量。

当加入林分和单木变量后，模型M0 的拟合优度较基础模型B1 显著增加约18%（以Ra2值为例，下同；表3）。对混交度、角尺度、Pielow 均匀度以及竞争指数4 个变量而言，当仅考虑1 个变量时，Heygi 竞争指数对模型的提升效果最为明显，其精度较模型M0 显著增加约7.48%；当考虑2 个指标时，混交度和Heygi 竞争指数的组合使模型精度提升约10.95%；当考虑3 个指标时，混交度&角尺度&Heygi 组合的效果最佳（11.21%）；当考虑4 个指标时，模型拟合效果较M0 增加约11.22%。方差分析结果表明：模型M0 与B1、M4与M1、M6 与M4、M14 与M6 均存在极显著差异（P＜ 0.000 1），而模型M14 与M15 间的差异则不显著（P= 0.415 7）。因此，在综合考虑各模型参数估计值、残差分布以及拟合优度等指标基础上，综合确定M14 为兴安落叶松单木直径生长的最优广义模型，其具体形式为：

式中：y(t)为第t年时树木的直径；DBH为解析木的胸径（伐倒时）；其余变量如前所述。

2.2 混合模型构建

在式5 基础上，引入随机参数、方差协方差结构、异方差校正函数和时间序列结构等构造兴安落叶松单木直径生长的混合效应模型（表4）。对随机参数而言，仅有8 种参数组合能够收敛；当考虑1 个随机参数时，模型E2 的AIC值最小（5 628.26）；当考虑2 个随机参数时，模型E8 的AIC值最小

（5 630.02），但较E2 略有增加（△AIC= 1.76），且似然比检验结果表明模型E8 与E2 间的差异也不显著（P= 0.627 6）。因此，将a1作为兴安落叶松单木直径生长混合效应模型的最优随机参数组合。

表3 基于林分、单木、物种多样性和结构参数的单木直径生长模型拟合优度Table 3 Goodness-of-fit for diameter growth model with stand-level, tree-level, diversity and structure variables

表4 兴安落叶松单木直径生长混合效应模型筛选结果Table 4 Selection results of diameter growth model of natural Larix gmelinii stand when using mixed-effects regression techniques

鉴于模型E2 中仅包含1 个随机参数，本研究不再进行混合模型方差协方差结构的选择。3 种异方差校正函数检测结果表明，常数加幂函数（模型E12）对异方差的校正效果显著优于幂函数（模型E10；P＜ 0.01）和指数函数（模型E11；P＜ 0.01）。但以常数加幂函数（模型E12）为基础构造的包含时间序列结构的混合效应模型均不能收敛，因此本研究选择次优的幂函数（即E10）为兴安落叶松单木直径生长混合效应模型的最适异方差校正函数。在模型E10 基础上，进一步测试4 种时间序列结构对模型拟合效果的影响，结果表明仅一阶自回归与滑动平均结构ARMA(1,1)和高斯结构corGaussian 能够收敛，且corGaussian 结构显著优于ARMA(1,1)结构（P＜ 0.01）；模型E13 与E14 间的差异同样达到极显著水平（P＜ 0.01）。因此，将模型E14作为兴安落叶松单木直径生长的最优混合效应模型。

兴安落叶松单木直径生长基础模型、广义模型和混合模型的参数估计以及拟合结果如表5 所示。除混合模型参数a0和b2外（P= 0.076 4 和0.095 3），其余各参数均达到显著水平（P＜ 0.05）。3 个模型统计结果表明：混合模型的Ra2分别较基础模型和广义模型增加约54%和17%，而平均绝对误差和均方根误差分别下降约72%、70%和59%、57%；混合模型不仅显著缩小了残差的范围，而且异方差现象也得到了显著削弱（图2），表明模型E14 能够满足兴安落叶松单木直径预测的需求。

表5 兴安落叶松单木直径生长基础模型、广义模型和混合模型的参数估计和拟合结果Table 5 Parameter estimation and fitting results of base-, general- and mixed-model for natural Larix gmelinii stand

2.3 模型检验与模拟

五折交叉检验结果表明，兴安落叶松单木直径生长混合模型的精度同样显著高于基础模型和广义模型（表6），能够满足林业生产的需要。结合表1，本研究假设林分平均胸径Dg= 10 cm, 林分平均高Hg= 12 m, 单木胸径DBH= 15 cm，角尺度W= 0.5，混交度M= 0.5 以及林分平均竞争指数H= 0.3 为一组基础的模拟情景，即在每次模拟中仅改变1 个变量的值而保持其他不变。以混合效应模型E14的固定效应参数为基础，模拟结果表明单木直径生长量随林分平均胸径、平均树高以及竞争指数的增大而减小，但随着单木胸径、角尺度和混交度的增加而增加（图3）。整体来看，各变量对单木直径大小的贡献量（定义为各变量每增加1 个单位时对应的胸径增长或减少量）依次表现为：混交度（+4.71 cm） ＞ 角尺度（+3.62 cm） ＞ 竞争指数（-0.73 cm） ＞ 单木胸径（+0.45 cm） ＞ 林分平均胸径（-0.34 cm） ＞ 林分平均树高（-0.21 cm）。

图2 兴安落叶松单木直径生长基础模型、广义模型和混合模型的残差分布Fig. 2 Residual distribution of base-, general- and mixed-model for natural Larix gmelinii stand in Daxing′anling Mountains

表6 兴安落叶松单木直径生长基础模型、广义模型和混合模型的独立样本检验Table 6 Validation of base-, general- and mixed-model for natural Larix gmelinii stand in Daxing'anling Mountains

3 讨论

在天然异龄林中，光照是限制树木生长的重要因素之一。显然，大树较小树具有显著的竞争优势，因此初始胸径的大小可充分代表树木在林分中的竞争态势，即一定范围内树木胸径越大，其直径生长量也越大。根据模型E14 测算结果，树木初始胸径每增加1 cm 时，单木直径生长量可平均增加4.50 cm 左右（生长至150 a 时，下同）；测算结果表明，林分平均胸径每增加1 cm，单木直径生长量将会下降约3.40 cm，而林分平均树高每增加1 m，单木直径生长量则会下降2.11 cm。综上所述，可通过适当的采伐来降低林分整体的竞争水平，这对促进林木直径生长具有重要作用。

除了上述基本林分和单木因子外，本研究还进一步测试了4 种物种多样性指数和4 种林分空间结构参数对模型拟合精度的影响（图1）。对天然林而言，林分内不同树种的搭配是经过长期的自然选择形成的，因此各树种间存在较强的优势互补关系，这是进行物种多样性与林分生产力关系研究的理论基础。结果表明，Mitscherlich 模型参数a与所有物种多样性指数均不相关，而参数b则与Sw呈负相关关系、与Pw呈正相关关系，这与谭凌照等[9]研究结果相似。鉴于Sw与Pw间显著相关（r=0.77），因此本研究仅尝试将Pw引入到模型中，但结果表明模型（M15）拟合效果并不理想，其与模型M14 间的差异不显著（P= 0.415 7）。因此，本研究所建模型中未包含任何生物多样性指标，但这可能与研究区域物种多样性水平整体偏低有关。在物种多样性丰富的地区，该结论是否成立仍有待于进一步检验。此外，本研究仅考虑了林分物种多样性指标，而与林木生长密切相关的结构多样性并未充分考虑[16]，这也有待于进步研究。

林分空间结构参数是近年来国内外森林经营领域研究的热点。在此基础上，惠刚盈及其团队提出的结构化森林经营理论被认为是实现森林可持续经营的重要途径之一[17]，但现阶段还很少有将其整合到传统林木生长模型的研究中。汤孟平等[10]、吕延杰等[11]虽然揭示了林分空间结构参数与林木生长的关系，但并未给出实用可行的预测模型。本研究在他们的基础上，进一步研究发现Mitscherlich 模型参数a与混交度（r= -0.23）、大小比（r= -0.19）和角尺度（r= -0.26）间均存在较强的相关性，而参数b则仅与Heygi 竞争指数相关（r= -0.12）。经过综合评估，本研究试图将混交度和角尺度引入到参数a中，而将Heygi 竞争指数引入到参数b中，结果表明各模型Ra2值分别+3.77%、+0.30%和+7.48%。因此，最终将混交度、角尺度和Heygi竞争指数引入到模型中，模型模拟结果表明林木混交程度越高、水平分布越随机，其直径生长量越大；林木竞争越激烈，其直径生长量越小，这些均与前人研究结果一致[10-11]。模型E14 模拟结果表明，混交度、角尺度和Heygi 竞争指数的数值每增加10%，其对应的直径生长量平均分别为+4.71 cm、+3.62 cm 和-0.73 cm，但因不同的林分特征和生长阶段而异（图3）。

图3 不同变量对兴安落叶松林单木直径生长的影响Fig. 3 Effects of different variables on the diameter growth of natural Larix gmelinii Stand in Daxing'anling Mountains

与传统统计学模型相比，混合模型能够通过固定参数反应总体的变化趋势，而通过随机参数来反映树木个体生长的差异。此外，混合模型还能提供其他一些丰富的信息，如方差协方差结构能够反应不同林分条件下树木直径生长过程的差异，这对复杂的异龄混交林而言具有明显优势；模拟结果表明，兴安落叶松单木直径生长的混合模型仅包含1个随机参数，其方差异方差函数则能有效解决解析木数据中的时间自相关性，进而确保得到的参数估计值是无偏的；模拟结果表明兴安落叶松单木直径生长混合模型的最适异方差校正函数为幂函数。模型模拟和检验结果表明，相较于广义模型（M14），混合模型E14 的各项指标和残差分布图均得到显著改善，这与其他学者的结果一致[12-13]。但需要强调的是，混合模型的检验和应用过程较传统模型要复杂的多，必须借助少量样本、模型参数估计结果（表5）和式（2）实现不同林分条件下树木个体直径生长过程中随机参数的确定，因此开发具有针对性的应用软件是当务之急。此外，由于样本量不足的问题，本研究仅考虑了单木水平的随机因子，而没有同时考虑区域、样地和单木等多水平的随机效应，随着数据的积累，可在这方面继续深入研究。

4 结论

林分结构显著影响兴安落叶松单木的直径生长过程。当引入林分角尺度和竞争指数后，广义基础模型M0 的Ra2值较基础模型B1 增加约18%；当进一步引入混交度、角尺度和竞争指数3 个结构参数后，广义最优模型M14 的（Ra2）值进一步增加约11%；而混合模型（E14）则使模型精度进一步提升约17%，且模型异方差得到显著削弱。模拟结果表明：兴安落叶松单木直径生长随林分平均胸径、平均树高和平均竞争指数的增加而减小，但随单木胸径、林分角尺度（接近随机分布，即W=0.5）和林分混交度的增加而增加。各变量对单木直径大小的贡献量（各指标每增加10%）依次表现为：混交度（+4.71 cm） ＞ 角尺度（+3.62 cm） ＞竞争指数（-0.73 cm） ＞ 单木胸径（+0.45 cm） ＞林分平均胸径（-0.34 cm） ＞ 林分平均树高（-0.21 cm）。因此，在后续森林经营中，应通过采伐逐渐降低林分的竞争水平，注重引入阔叶树种增加林木的混交程度，并逐步调整林木的水平格局为随机分布模式。