凉水自然保护区阔叶红松林冠下4个主要树种更新幼树树高生长模型1)

李武赫 董利虎 李凤日

(东北林业大学，哈尔滨，150040)

天然更新是天然林演替的重要生态过程，影响着森林群落结构、功能和生物多样性。如果没有足够幼苗幼树的补充，天然林上层林木组成的森林结构就会不稳定，更不会形成连续覆盖、稳定健康的森林生态系统。林冠下不同乔木的更新数量、分布格局和生长决定着未来林分的树种组成和林分结构。因此，研究天然林中幼树的生长情况对天然林演替和可持续经营的研究有着极为重要的意义[1]。

目前，在天然更新的研究中，研究对象大多是天然更新株数。国内外学者的研究表明,天然更新株数与多个林分因子有关，如郁闭度、土壤条件、坡度、坡向、海拔、林分密度等[2-6]，这些研究为天然更新株数的预测提供了参考，但无法从中了解成功更新幼树的生长情况。为了进行更深一步的研究，需要以树木生长的一个指标来量化林冠下更新幼树的生长情况。树木的树高和胸径是林业常用的调查因子，而低龄幼树树高往往不达1.3 m，没有胸径。虽然有以地径替代胸径作为幼树调查因子的做法，但由于树木干形的特点，这种做法并不能很好地反映幼树的生长情况。此外，生物量也可反映幼树的生长情况，但对其的研究需要挖取并破坏幼树[7]，不适用于原始林。可见树高是反映幼树生长状况的最佳数量指标。有学者通过建立幼树树高混合效应模型来描述幼树的生长过程，指出林分优势高、林分胸高断面积和林分平均胸径对不同种幼树有不同程度的影响[8]，但幼树树高生长规律总体不明显，随机效应(即幼苗自身的优劣)对幼树树高生长的影响较大[9-10]。混合效应模型近代发展较快且在林业中应用广泛。其在传统生长模型中引入可以体现样地或样木自身差异性的任意多个参数为随机效应参数，既可以反映总体平均变化趋势，也可以反映个体之间的差异，在观测值充足的情况下可充分预测树高生长过程[11-12]。依此法建立的生长模型比传统模型的预测精度大幅度提高。

阔叶红松林是中国东北地区的地带性顶级群落，在中国东北广泛分布。因此，原始阔叶红松林冠下幼树的生长情况对东北地区天然次生林的经营具有重要的参考意义。红松、红皮云杉、冷杉、色木槭是黑龙江省凉水自然保护区原始阔叶红松林林冠下更新的常见树种。它们在幼年时耐阴性较强，易于在遮阴条件下存活，对今后林分结构有着重要的影响[13-15]。本研究以黑龙江省凉水自然保护区原始阔叶红松林为研究对象，分析其林冠下红松、红皮云杉、冷杉、色木槭等主要更新幼树的树高生长规律，建立4个树种DB，H(胸径)<5 cm幼树树高生长的基础模型和混合效应模型。目前，国内外研究单木生长模型主要针对胸径大于5 cm的树木，对林冠下未达到起测径阶(5 cm)幼树的树高生长模型研究较少。本研究可以揭示该地区主要树种幼树树高生长规律，完善该地区主要树种单木生长模型，为天然林的生长模拟和森林经营优化提供基础。

1 研究地区概况

黑龙江省凉水国家级自然保护区位于我国小兴安岭山脉的东南段，归黑龙江省伊春市带岭区管辖。地理坐标为东经128°47′8″～128°57′19″，北纬47°6′49″～47°16′10″。该地区地带性土壤为暗棕壤。保护区内森林类型丰富，天然林包括处于演替顶级阶段的原始阔叶红松林、兴安落叶松林、云冷杉林和因受干扰而处于不同演替阶段的次生白桦林、山杨林、白桦山杨林、硬阔叶林、杂木林等；同时还存在红松、落叶松、云杉、樟子松等树种的人工林。地带性植被是以红松为主的温带针阔混交林，优势种为红松(Pinuskoraiensis)，主要伴生树种有冷杉(Abiesnephrolepis)、红皮云杉(Piceakoraiensis)、鱼鳞云杉(Piceajezoensis)、水曲柳(Fraxinusmanshuric)、色木槭(Acermono)、青楷槭(Acertegmentosum)、紫椴(Tiliaamurensis)、裂叶榆(Ulmuslaciniata)、枫桦(Betulacostata)、白桦(Betulaplatyphylla)、赤杨(Alnussibirica)、暴马丁香(Syringareticulata)等。

2 研究方法

2.1 数据来源

数据来源于2016年黑龙江省伊春市带岭区凉水国家级自然保护区调查数据。在该地区不同林分条件(森林类型、立地条件、林分密度等)的原生阔叶红松林中设置14块20 m×60 m临时标准地，所选标准地林木皆生长正常且具有代表性。测定各标准地郁闭度、坡度、坡位、坡向、海拔等因子，记录标准地坐标。

在不同的研究中，各学者对幼树的定义各不相同[16-18]，文中依据研究目的，不考虑树木的树龄，而是将胸径小于5 cm的树木全部作为幼树并研究其生长规律。对各标准地中具代表性的部分红松、云杉、冷杉、色木槭幼树(胸径<5 cm)进行每木检尺，测定幼树树高(h)、胸径(DB，H)、地径(DG，H)、冠幅(CW)，并测定幼树树龄(t)以及每一年的树高生长量(Δh)，具体的测定方法为，针叶幼树树龄和幼树年高生长根据其每年长一轮轮枝的规律判断和测量[19-20]，色木槭幼树则根据芽鳞脱落后留下的痕迹判断和测量[21-22]。对数据进行整理汇总，大树的胸径5.0～107.5 cm，幼树的胸径为0～4.91 cm，幼树的地径为0.11～7.37 cm，幼树的树高为5～480 cm。考虑到幼树生长过程中，已经长出的轮枝(针叶树)和芽鳞痕(阔叶树)高度不会发生变化，故文中将每一株幼树每一年的树高—树龄对应关系均用来建立模型(每一年对应的树高可利用树高生长量(Δh)推算)。各树种幼树的样木株数、用于建模的树龄—树高对应关系数(样本数)、树龄、树高以及树高生长量数据详见表1。

将4个树种的全部数据按4∶1的比例分成建模数据和检验数据，分别用于模型的建立和模型的独立检验，具体建模数据和检验数据的基本信息统计见表2。

表1 4个主要更新树种幼树的树龄、树高及树高生长量统计表

表2 4个树种幼树树高生长模型建模数据与检验数据树龄和树高统计量

2.2 树高生长基本模型的构建

作出各树种幼树的树高与年龄的散点图并加以分析，确定树高随年龄的变化规律。在此基础上，为了更深入的了解幼树树高的生长规律，需要建立树高生长模型来量化林下幼树的生长。由于幼树树高随着年龄的增大而增大且没有最大值的限制，目前常用于描述树高生长的理论方程(Richards方程等)不适用于幼树树高的研究[23]。文中参考前人对幼树树高生长的研究以及幼树树高生长的规律，选用以下5个方程作为4种树种幼树树高生长的备选方程(表3)。

表3 4个树种树高生长备选模型

2.3 树高生长混合模型的构建

利用传统方法拟合出的树高生长模型只能预测总体的平均水平，不能体现组间的差异，而混合效应模型则可以解决这个问题。由于混合效应模型可对组水平和个体水平的数据进行同时分析，在一个模型中可同时检验组变量和个体变量的效应[24]，混合效应模型在拟合树高生长模型上得到了广泛应用[25-28]。以选出的最优模型为基础，通过再参数化加入林分因子，并构建线性或非线性混合模型。拟采用样地单水平、样木单水平的线性或非线性混合模型，以单水平非线性混合模型为例，具体形式如下：

(1)

式中：yij为第i个区组中的第j次观察值；M为区组数量；ni为在第i个区组上观测的次数；f为含有参数向量(φij)和协变量向量(vij)的函数；β为(p×1)维固定效应向量；bi是带有方差—协方差矩阵D的(q×1)维随机效应向量；Aij和Bij为相应的设计矩阵；eij为服从正态分布的误差项；σ2为方差；Ri为区组i的方差—协方差矩阵；D为随机效应的方差—协方差矩阵。

对于已经确定的基础模型，构建混合效应模型首先需要确定参数效应，即分别确定仅含有固定效应的参数和既含有固定效应又含有随机效应的参数。文中认为所有参数都有成为随机效应参数的可能，尝试不同随机参数的组合分别拟合，比较模型拟合统计量，赤池信息准则(AI，C)越小、贝叶斯准则(BI，C)越小、对数似然值(LL)越大，则拟合效果越好，即为最优模型，不同随机参数个数的最优模型要进行似然比检验(LR，T)以选取最优模型。AI，C、BI，C、似然比检验的公式如下：

AI，C=-2LL+2d；

(2)

BI，C=-2LL+dlnn；

(3)

LR，T=-2ln(L2-L1)=2(lnL2-lnL1)。

(4)

式中：LL为最大似然函数的对数值；d为模型中估计参数个数；n为有效观察个数；L1、L2分别为2个模型的似然值，并且L2>L1，似然比检验近似卡方分布，自由度为参数个数差，若LR，T>χ2(df，α)，则L2相对应的模型更优。

2.4 模型检验与评价

对模型进行独立样本检验来评价模型的预测能力。混合效应模型的检验稍不同于普通模型的检验，其中固定效应部分参数与普通模型相同，随机效应部分参数则需根据检验数据重新计算，含有随机效应的参数需将固定效应参数与重新计算的随机效应参数相加得到。根据以下公式来计算随机参数(bk)：

(5)

(6)

采用调整后的确定系数、均方根误差对基础模型的拟合结果进行比较并选出最优模型。模型检验的评价指标是平均误差、平均绝对误差、平均预测误差百分比。所有指标中，调整后的确定系数越大、均方根误差越小，说明模型的拟合优度越好；平均误差越靠近0；平均绝对误差和平均预测误差百分比越小说明所构建的预估模型精度越高。计算公式如下：

调整后的确定系数：

(7)

其中：

(8)

均方根误差

(9)

平均误差

(10)

平均绝对误差

(11)

平均预测误差百分比

(12)

3 结果与分析

3.1 幼树树高生长规律

4个树种幼树树高和树龄的关系见图1。由图1可以看出，幼树树高随着年龄的增大而增大，二者呈明显的正相关关系。色木槭的连年树高生长量较大，变异范围相对较小；红松、云杉、冷杉3个针叶树种连年树高生长量较小，变异范围相对较大。在幼树年龄较大时，冷杉的数据点稍显离散，但4个树种总体趋势都较为明显。此外，4个树种在幼树生长过程中，树高的年增长量差别不大。这说明Richards方程、Gompertz方程等目前应用较多的树高生长方程并不适合描述幼树的树高生长，而线性方程、幂函数等形式的方程相对而言更为适合。

3.2 最优基础模型

图1 4个树种幼树树高—树龄散点图

表4 4个树种幼树树高生长基本模型拟合统计量

3.3 单水平混合效应模型

以选出的最优模型(M5)为基础，试图通过再参数化的方法加入林分因子。结果表明，林分因子(每公顷株数、每公顷胸高断面积、郁闭度等)虽然与幼树树高生长存在较弱的相关性，但是加入模型后对模型的精度提升十分微弱甚至几乎没有影响。可知影响幼树生长的主要因素并不是林分因子，而是幼树自身的差异，即随机效应。为了反映幼树自身的差异，本研究分别构造样地单水平、样木单水平的混合效应模型。为了方便区分混合模型的固定效应参数和随机效应参数，将4个树种基础模型的固定效应参数a、b、c替换为β1、β2、β3，并以b1、…、bn表示混合模型的随机效应参数。建立混合模型后，采用赤池信息准则、贝叶斯信息准则、对数似然值和似然比检验等指标来分析不同水平模型的拟合效果，最后选择模型收敛并且精度较高的形式。由于此步骤包含样地单水平的混合模型，而本研究的建模、检验数据未按照样地来抽取，故此步骤利用所有数据进行建模，仅用于模型之间的比较和后续最优模型形式的选择。不同混合模型混合作用比较见表5。

表5 4个树种幼树基于不同随机效应的树高生长混合效应模型结果

续(表5)

每个树种可分别建立14个不同水平、不同随机效应参数组合的混合模型。其中，红松、冷杉、色木槭的14个混合模型均收敛；云杉的14个混合模型中有13个收敛。由于篇幅较大，表5仅列出4个树种在不同水平、不同参数个数下的最优模型。从拟合结果可以看出，当加入样地效应时，红松、云杉、冷杉、色木槭等4个树种的AI，C、BI，C值均有较小幅度的减少，LL值增大；而4个树种加入样木效应后的混合模型的AI，C、BI，C值都比无随机效应和仅加入样地效应的模型有较大幅度的减少，LL值大幅度增加。证明4个树种都是在随机效应水平为样木水平时的拟合效果最好。在此基础上分别选择4个树种的最优模型。

将红松样木水平混合模型中含一个随机效应参数的模型1-2-1和原模型(模型1-0)做似然比检验，发现模型之间有显著差异，且模型1-2-1的AI，C、BI，C值较小，LL值较大，拟合效果较好。将模型1-2-1与含两个随机参数的模型1-2-2做似然比检验，存在显著差异，模型1-2-2的AI，C、BI，C值较小，LL值较大，拟合效果相对较好。将模型1-2-2与含3个随机参数的模型1-2-3做似然比检验，存在显著差异，模型1-2-3的AI，C、BI，C值较小，LL值较大，拟合效果最好。由此确定，含3个随机效应参数b1、b2、b3的模型(模型1-2-3)为红松幼树树高生长的最优模型。同理选出云杉、冷杉、色木槭的最优混合效应模型，分别为模型2-2-3、3-2-3、4-2-3，都为含3个随机效应参数b1、b2、b3的模型。

加入样木效应的模型与加入样地效应的模型相比，AI，C值明显减小，这也证明了4个树种的随机效应主要源于样木之间的不同，即影响幼树树高生长的一个主要因素是幼树个体之间的自身差异。4个树种模型形式相同，如式13所示：

h=(β1+b1)t(β2+b2)exp[-(β3+b3)t]。

(13)

式中:β1、β2、β3为固定效应参数；b1、b2、b3为随机效应参数。

确定最优模型形式后，利用建模数据对数据进行建模。表6列出最优混合效应模型的参数结果。

表6 4个树种幼树树高生长混合效应模型参数估计值

作出4个树种基础模型、混合效应模型的预测值和实测值散点图(图2)。图2中同一行a、b两幅图分别对应同一树种的基础模型和混合效应模型的结果。从图2中可以看出，4个树种基础模型的数据点较为离散，而混合效应模型的数据点集中在截距为0、斜率为1的直线两侧，说明混合效应模型能较好的对各个树种的幼树年龄和树高的关系进行拟合。其数据点聚集集中，变异范围小，拟合效果更好。

a.基础模型；b.混合效应模型。

3.4 模型检验

采用独立检验数据分别对4个树种的基础模型和混合效应模型进行检验，计算平均误差(ME)、平均绝对误差(MA，E)和平均预测误差百分比(MP，S，E)，结果见表7。

从表7中可以看出，4个混合效应模型的平均误差都比基础模型的要接近0，其中红松、云杉、冷杉3个树种的平均误差为负值，说明混合效应模型对这3个树种的预测略微偏高；色木槭的平均误差接近0，预测误差很小。与基础模型相比，红松、云杉、冷杉、色木槭的预测精度都有较大程度的提高，混合效应模型的MA，E分别降低了9.63、9.28、13.37、6.85 cm；MP，S，E分别降低了26.49%、33.00%、49.51%、14.40%。结合表4、表6和表7可以看出，混合效应模型的拟合结果和检验结果都远好于基础模型。说明混合效应模型能够提高模型的拟合优度和预测精度。

表7 4个树种树高生长基础模型和混合效应模型检验结果

4 结论与讨论

森林的天然更新是天然林繁衍的重要手段。林冠下不同乔木的更新状况对未来林分的树种组成、森林群落结构等方面起着重要作用。因此，对其的研究就显的尤为重要。本研究结果表明,凉水地区阔叶红松林冠下幼树的树高生长主要受随机效应(其中主要为幼树个体之间的差异)影响，而受林分因子影响较小。在确定柯列尔方程为最优基础模型后建立主要树种树高生长的混合效应模型，较好地描述了凉水地区阔叶红松林冠下主要树种的幼树树高生长。

由红松、云杉、冷杉、色木槭4个树种的树高—树龄散点图得知，4个树种幼树树高随树龄增加而增大且总体趋势较为明显。此外，4个树种在幼树生长过程中，树高的年增长量差别不大。在对树高生长的研究中，较为常用的方程是Richards方程、单分子式、Logistic方程等存在一条渐近线的理论方程。此类方程拥有渐近线，可以解释当树木年龄无限大时树高几乎停止增加并趋近于最大值的现象。但由于树木树高在整个幼树阶段持续增加，这些理论方程不适合描述幼树生长过程。本研究针对幼树树高的生长特点，以不受最大值限制的5个备选模型拟合了幼树树高随树龄的变化过程。结果表明，指数方程的拟合效果相对较差，其余4个模型均可以较好地拟合幼树的树高生长，最优模型为柯列尔方程(M5)。

通过再参数化的方式分析林分因子对模型的影响，结果表明，林分因子对4个树种幼树的树高生长的影响不大，故本研究未在模型中添加林分因子。可知幼树树高随树龄变化的差异主要来自随机效应，而幼树自身的优劣是造成随机效应的最主要因素，与前人的结果一致。在模型中引入随机效应参数以反映幼树树高生长的随机效应，结果表明，样木水平的非线性混合效应模型最优。它相当于为每一株树设置一组参数，减小模型拟合和预测过程中样木个体间的差异，提高模型精度。模型的检验结果显示混合效应模型的精度较普通模型有很大提高，说明该混合效应模型可以很好的反映树木之间生长的差异，也再次证明了幼树的树高生长主要与幼树自身的优劣有关，而与林分因子关系相对较小。以往的研究中仅对胸径大于5 cm的大树建立相关生长模型，对于幼树树高生长模型的研究较少，而本研究补充了胸径小于5 cm部分幼树的生长过程，可完善该地区红松、云杉、冷杉、色木槭单木生长模型，并为天然林的生长模拟和森林经营优化提供科学的依据。

随机效应产生的原因是各树种幼树之间的个体差异，造成个体差异的最主要原因可能是各树种种子存储的营养不同，且在幼树阶段从外界吸收养分的能力差异较大，影响幼树的成活和幼树树高的生长。文中并未考虑造成个体差异的原因。对幼树树高研究的趋势可能是在树高生长模型中加入个体差异具体造成的影响，需要长期的观察与试验。