长白落叶松人工林单木冠幅模型1)

邱思玉 孙玉军

(森林资源和环境管理国家林业局重点开放性实验室(北京林业大学)，北京，100083)

树木进行生长发育时很大程度要依靠树冠来吸收阳光和各种养分[1]，树冠的形状和大小会影响树木生长和生产力的分配，还会影响树干生长的质量和质地[2-3]。同时，从树冠的生长情况我们能够分析出林木之间的竞争状态，也能反映出林分的健康水平[4-6]。在用来形容冠层生长的各项林分因子中，冠幅显得尤其重要。预测冠幅生长有利于计算林冠郁闭度等指标，有利于评估森林生物在栖息地的适宜程度，有利于预测火灾风险等级和林下植被的更新情况[7]。同时，冠幅也可以用来预测树冠的表面积和体积等[8-9]。此外，冠幅也是林木可视化的一项重要参数[10]。可见冠幅有巨大的应用价值，但是测量每株单木的冠幅会耗费巨大的人力物力，因此，研究冠幅的相关预测模型对于研究树木生长具有重要意义[11]。

早期，关于冠幅的研建主要是利用传统回归分析方法构建冠幅与林分因子的线性关系，添加多种林分变量提高模型精度[12-16]，但是实际工作中这种关系通常不只是线性关系所能表达的，可能会由于重复观察数据和多水平观察数据造成异方差和自相关等问题，相比之下，非线性混合效应模型能够解决这类问题，且其组合形式多样，更贴近实际工作[17-18]。Sanchez et al.[19]将样地作为随机效应，构建了西班牙栓皮栎的冠幅非线性混合效应模型,符利勇等[17]基于立地指数和样地二水平，构建了杉木的冠幅混合模型,Sharma et al.[20]考虑样地的随机效应对捷克的挪威云杉和欧洲山毛榉构建了单木冠幅模型,Raptis et al.[21]以幂函数为基础模型构建了黑松天然纯林的非线性混合效应模型,吕乐等[22]将模型中所有参数看成混合效应参数，拟合不同随机参数混合的模型构建了黑龙江省东部地区天然椴树的单木冠幅预测模型。目前利用这种方法构建长白落叶松人工林单木冠幅模型的研究报道相对较少，本研究一定程度上弥补了该方面的空白，有助于预测树木生长、生产力和评价木材品质等级等，能够为东北地区长白落叶松的科学经营管理提供依据，在林业生产上有重要意义。

1 研究区概况

研究区位于黑龙江省伊春市朗乡林业局东折棱河林场(128°55′30″～129°15′21″E，46°31′58″～46°49′38″N)。该林场属低山丘陵区，阳坡短且陡，阴坡长且缓，平均坡度12°，平均海拔约600 m。土壤主要是暗棕壤，土层厚度平均60 cm，土壤肥沃。气候属北温带大陆湿润性季风气候，冬季寒冷干燥持续时间较长，夏季温暖湿润持续时间短，年平均气温1.9°，年平均降水量约620 mm，该林场森林覆盖率达92%，其中落叶松人工林以长白落叶松为优势树种[23]。

2 材料与方法

2.1 数据来源

数据源于2008—2019年黑龙江省伊春市朗乡林业局东折棱河林场，长白落叶松人工林固定样地数据35块(面积600 m2的方形样地25块，半径7.32 m的圆形样地10块)，临时样地数据14块(面积600 m2)。对样地内树木进行每木检尺，起测径阶(D)为5 cm，还调查了林分起源、林分年龄(A)、林木胸径(D)、树高(H)、第一活枝高(HCB)、东西南北4个方向的冠幅。其中，树高是用超声波测高仪测量，冠幅是用激光测距仪站在每个方向的树冠最大投影距离处测量到树干的垂直距离。此外，还计算了冠长(CL)、冠长率(CR)、高径比(HDR)等。共调查长白落叶松2 499株，将49块样地数据按3∶1的比例随机选取37块样地(1 807株)作为建模数据，12块样地(692株)作为检验数据，统计信息见表1。

本研究中的单木冠幅是指东西南北4个方向冠幅直径的平均值[24]，计算公式：

CW=(CWE+CWW+CWS+CWN)/2。

(1)

其中：CW为总冠幅，CWE、CWW、CWS、CWN分别为东、西、南、北4个方向的冠幅。

表1 建模数据和检验数据统计信息

2.2 基础模型

从10个常见的冠幅胸径模型中(见表2)[24-27]用赤池信息量准则(AIC)、剩余均方根误差(RMSE)和调整决定系数(Radj2)3个指标对模型进行评价，选用精度最好的模型作为构建改进模型的基础模型，赤池信息量准则用R语言的AIC函数实现，其余两个指标的计算公式见：

(2)

(3)

计算10个模型的AIC、RMSE和Radj2，为避免模型参数过多，选出精度最好的二参数模型和三参数模型作为备选模型，用方差分析判断二者显著性是否差异，若显著选择三参数模型，不显著则选择二参数模型作为基础模型[28]。

表2 备选基础模型

2.3 林分变量

除胸径外，冠幅也受其他单木或林分变量的影响，如树高、冠长等林分测树因子，立地因子，优势木平均高等林分竞争因子，因此，在模型中增加其他变量有助于提高模型精度，本研究考虑的因子如下。林分测树因子：树高(H)、枝下高(HCB)、冠长(CL)、冠长率(CR)、高径比(HDR)、林分年龄(A)、林分密度(M)；立地因子：立地指数(SI)；林分竞争因子：样地算术平均胸径(AMD)、样地平方平均胸径(QMD)、优势木平均高(HDOM)、优势木平均胸径(DDOM)、相对植距(RSI，见式(4))、每公顷断面积(BA)[29-32]。

(4)

式中：RSI为相对植距；M为林分密度(株/hm2)；HDOM为优势木平均高。

2.4 随机效应模型

本研究根据样地的立地类型将数据分为7个区组，考虑区组和区组内嵌套的样地两个水平，将其作为随机效应，用R软件的nlme包进行拟合，模型表达式如下：

yijk=f(φijk,vijk)+εijk，

i=1、…、M；j=1、…、Mi；k=1、…、nij。

(5)

式中：yijk是第i区组第j样地第k株树的冠幅观测值；f是含有参数向量φijk和协变量向量vijk的非线性函数；M是区组数，Mi是第i区组的样地数，nij是第i区组第j样地的单木株数；εijk是表示组内方差和相关性的误差项，假定其是含有零值和Rij(方差协方差结构，由参数向量λ表示)的独立正态分布，表示为εijk～N(0，Rij(λ))。

φijk=Aijkβ+Bijkui+Mijkuij，ui～N(0,ψ1)，

uij～N(0,ψ2)。

(6)

式中：Aijk、Bijk和Mijk为设计矩阵；β是表示固定效应的p×1维矩阵；ui是表示第一水平随机效应的q1×1维矩阵，是含有零值和ψ1(方差协方差矩阵)的独立正态分布；uij是表示第二水平随机效应的q2×1维矩阵，是含有零值和ψ2(方差协方差矩阵)的独立正态分布。

2.5 模型求解和检验

本研究中模型的计算是基于实测样地数据，应用R语言的nls和nlme函数实现，利用建模数据对各模型参数进行拟合，利用检验数据进行检验。应用AIC、BIC、-2LL等指标对模型进行估计，用平均偏差(Bias，见式(7))、RMSE以及Radj2对检验数据进行检验，通过比较分析得出构建该地区长白落叶松人工林的非线性混合效应模型。

(7)

3 结果与分析

3.1 基础模型评价

将样地数据代入10个常见的冠幅胸径模型，利用R语言的nls函数拟合求出模型参数值。用AIC、RMSE以及Radj2对建模数据进行评价，用Bias、RMSE和Radj2对检验数据进行评价(结果见表3)。其中M4、M7、M10为三参数模型，其余为二参数模型。AIC和RMSE的值越接近0，Radj2的值越接近1，模型的精度越好，故二参数模型中M1模型精度最高，三参数模型中M4模型精度最好。

对这两个模型进行方差分析(F=13.454，p<0.01)，可知二者差异性显著，故选M4作为本研究建模的基础模型，形式为式(8)：

CW=a/(1+bexp(-cD))。

(8)

式中：CW为冠幅;D为胸径，a，b，c为参数。

表3 备选基础模型评价指标

3.2 模型改进

为提高模型精度，在模型中增加对冠幅影响较大的单木或林分变量，经Pearson相关性分析得出树高(H)、冠长率(CR)和高径比(HDR)与冠幅的相关性较大，故将三者加入以构建改进模型，共有6种参数分配形式，经分析得出下列参数拟合形式为最佳(见式(9))。

CW=(a+a1×H)/(1+(b+b1×CR)×

exp(-(c+c1×HDR)×D))。

(9)

式中：CW为冠幅，D为胸径，a=12.841 9，a1=-0.196 9，b=10.890 3，b1=-6.403 6，c=0.060 0，c1=0.033 9为参数。

3.3 随机效应模型参数

CWijk=(a+(a1+ua1i+ua1ij)×CR)/(1+(b+b1×HDR)×

exp(-(c+c1×H)×D))+εijk。

(10)

式中：a、a1、b、b1、c、c1为固定效应参数，ua1i是区组效应作用在a1的随机效应参数，ua1ij是区组和样地效应作用在a1的随机效应参数。

用Rstudio的nlme包对模型进行拟合，当随机效应作用在参数a1时，区组水平的随机效应方差为0.000 6，样地水平的随机效应方差为0.001 1，模型拟合的随机效应和固定效应参数以及相应指标如下，见表(4)和表(5)：

表4 随机效应拟合结果

表5 固定效应拟合结果

3.4 模型评价

绘制改进模型式(9)和嵌套二水平非线性混合效应模型式(10)的拟合效果图(图1)，可见，混合模型的拟合效果图较比改进模型来说更趋向于线性相关，散点分布更趋于一条直线，说明预测值和实测值数值更接近，因此拟合效果更准确。

为了分析模型的异方差性，对式(9)和式(10)分别绘制残差分布图(图2)，可见两个模型的残差都围绕零值均匀分布，区间在-2到2之间，不具有明显的喇叭形状，说明没有明显的异方差性。因为数据不存在重复测量，所以没有自相关性。同时，混合模型的残差散点较比改进模型来说分布更为集中，拟合效果更好。

图2 模型残差图

由图3分析两个模型的密度直方图，其中高耸的蓝色区域为混合模型冠幅预测值的密度分布，低缓的粉色区域为改进模型冠幅预测值的密度分布，由此可见，混合效应模型的冠幅预测值密度分布相对来说更为集中，拟合效果更好。

图3 密度直方图

对改进模型和混合模型进行精度比较(见表6)，可以看出式(10)的AIC、BIC、-2LL、RMSE和Bias均小于式(9)，Radj2则高，故式(10)的模型拟合精度更优，可知在考虑了区组和样地嵌套二水平后模型的拟合能力得到提高，因此非线性混合效应模型式(10)能够更好地应用于长白落叶松人工林的冠幅模拟。

表6 模型评价指标

4 结论与讨论

树冠在森林的生长发育过程中起着至关重要的作用，冠幅又是描述树冠的一项重要指标。本研究发现，在常用的10种冠幅胸径模型中，Logistic模型具有最优的拟合能力，这与其他学者的研究结果一致[17，22]。在常见的单木或林分变量中，树高、冠长率和高径比与冠幅的相关性最大，可以考虑到建模中。

选用区组和样地两个水平作为随机效应加入建模后，分析数据的异方差性和自相关性，利用实测值与预测值之间的拟合关系图分析两种模型的冠幅预测情况与实测值的关系，判断预测值的真实程度，可以得出混合模型的拟合效果更好。利用残差图分析模型是否具有异方差，同时，从残差图也能分析出模型间的对比情况，由图像可以得出模型不具有异方差性且混合模型的残差值分布更为均匀密集，效果更好。从密度直方图可以分析两种模型冠幅预测值的密度分布情况，判断模型预测的离散程度，可以看出混合模型的预测值分布更对称和集中。从拟合指标表可以分析对比出两种模型的精度大小情况，混合模型无论是建模数据还是检验数据，其指标精度值都更好，因此，本研究构建构建的冠幅模型可以较好的拟合长白落叶松人工林的冠幅生长情况。