基于两水平非线性混合效应模型的长白落叶松削度方程构建

聂璐毅,董利虎,李凤日,苗 铮,谢龙飞

(东北林业大学林学院，黑龙江 哈尔滨 150040)

削度方程是以胸高直径、全树高及树干上各部位高度等变量为自变量，树干上各部位高度处的直径为因变量的数学函数。根据研究方法的不同，可划分为简单一次性拟合削度方程、分段削度方程和可变指数削度方程[1]。可变指数削度方程相对于简单削度方程和分段削度方程，其拟合精度更高、偏差较低[2]，虽然不可直接通过积分获得一致性材积方程，但是可以利用数值积分技术来达到高精度材积预测的效果。

树木的削度受树种、气候因素、年龄、树冠变量、林分密度等多种因素影响[3-5]。由于树冠与干形的生长发育关系密切，树冠大小通常被认为是描述干形变化的辅助变量[6]。一些研究提出树木个体的树冠结构会对树干各部位造成不同机械应力，导致树木的营养分配发生变化，从而影响树干的大小和形状[7]。Weiskittel等[8]提出树冠特征变量的加入不仅可改善模型拟合预测，并可确保模型更符合现实规律。Li等[9]在对北美针叶树种的干形进行研究时，发现引入活冠高能够提高削度方程对云杉(Piceaasperata)、冷杉(Abiesbalsamea)和白松(Pinusstrobus)的预测精度；Jiang等[10]在研究北美鹅掌楸(Liriodendrontulipifera)的干形变化时，将冠长率引入分段削度方程，发现冠长率能够提高模型对树干上部直径的预测精度。

在构建削度方程时，通常选择具有代表性的样地中多株树木，然后进行树干解析获取每棵树木不同位置处的直径。传统回归分析方法无法反映样地和样木的差异对树木生长的潜在影响，往往导致模型预测精度不高，而混合效应模型能有效地解决该问题[11]。国内外学者在对削度方程进行研究时，大多采用混合效应模型来建立削度方程。Shahzad等[12]在对白桦(Betulaplatyphylla)的干形进行研究时，使用混合效应模型解决数据的层次问题，从样地和样木两个水平对削度方程的建模过程进行了分析。Liu等[13]在分析气候因子对兴安落叶松干形影响的研究中，采用混合效应模型，引入指数函数和一阶连续自回归结构来定义组内误差结构，建立了气候敏感的削度方程。混合效应模型的个体预测需要通过二次抽样来对模型进行校正，以有效提高模型预测精度[14]。Cao等[15]通过胸高以上50%处的直径对模型进行校正，提出混合效应模型计算随机效应的校准方法显著提升了模型精度；Sharma等[16]在对黑云杉(Piceamariana)和北美短叶松(Pinusbanksiana)干形的研究中，提出以树干总高度的34%～38%测量树干上部直径进行模型校正最好。目前国内结合混合效应模型中的二次抽样对削度方程进行探讨的研究较少。

长白落叶松作为东北林区的主要树种之一，分布区域广泛，经济效益强，是建筑用材和工业用材的首选树种，因此对落叶松干形精准预测的研究具有重要意义。本研究以长白落叶松为对象，研究树冠变量与干形之间的关系，利用再参数化方法将树冠特征变量引入削度方程；结合混合效应模型并考虑样地、样木效应对树木干形的影响，建立长白落叶松削度方程，对比基础模型和混合效应模型的拟合检验精度；研究不同的二次抽样方案对混合效应模型预估精度的影响，为长白落叶松人工林干形的精准预测提供理论依据。

1 材料与方法

1.1 数据来源

本研究的数据分别来源于黑龙江省孟家岗林场、东京城林业局和林口林业局(129°13′～130°39′E，44°6′～46°25′N)。采集地区的气候为中温带大陆性季风气候，春秋短暂，冬夏分明，年平均降水量为550 mm左右，年平均气温2.7 ℃，极端最高气温35.6 ℃，山地土壤以典型暗棕壤为主。在此区域设置长白落叶松人工林固定样地共31个，对样地内树木进行每木检尺，按照等断面积标准木法在每块标准地附近选择5～7棵落叶松作为解析木。测定解析木的胸径、树高、冠幅、活枝高以及相邻木的树种、距离、树高和胸径等数据。伐倒后，测量从树冠梢头到树冠基部的距离，即冠长，并将样木基底至树梢的长度作为全树高，在相对高(0、0.02、0.04、0.06、0.08、0.10、0.15、0.20、0.25、0.30、0.40、0.50、0.60、0.70、0.75、0.80、0.90)处测量17个位置的带皮直径。长白落叶松人工林各标准地主要林分因子以及177棵解析木调查因子统计信息见表1。

表1 长白落叶松人工林各样地林分因子和样木因子统计量

1.2 基础削度方程

削度方程按照形式可以分为简单削度方程、分段削度方程以及可变指数削度方程。本研究对各类型的削度方程进行拟合，通过比较各模型的拟合统计量(调整后决定系数、均方根误差)选择最优基础模型。最终选择Kozak[17]可变指数削度方程作为基础模型，模型在拟合过程中删除一个不显著参数项，具体形式如下：

(1)

式中：d为解析木在距地面高度h位置处的带皮直径，D为胸径，H为全树高，q=h/H，X=(1-q1/3)/[1-(1.3/H)1/3],ε为模型的误差项，a1～a8为模型待估计参数。

1.3 削度方程再参数化

为了评价树冠特征对削度方程的影响，采用模型再参数化[18]的方法构建包含树冠特征变量(冠长、冠幅、冠长率)的长白落叶松削度方程。将树冠特征变量加入模型中，根据变量在模型中参数的显著性，选出拟合优度和检验结果最优的形式作为最佳的再参数化模型，将与削度方程关系最密切的变量加入削度方程中。

1.4 混合效应模型

非线性混合效应模型(NLME)[19]，包含单水平和多水平形式。本研究针对已建立的基础模型，考虑样木效应、样地效应两个水平的随机效应来建立混合效应模型，其模型形式如下：

(2)

式中：yij为第i块样地第j棵树在距地面高度h处树干的直径观测值的向量；f为包含参数向量Φij和协变量向量vij的非线性函数；m为第1水平的分组数，mi为第2水平的分组数；Xij为设计矩阵；β为固定参数向量；Zi,j、Zij分别为第1水平和第2水平的随机效应设计矩阵；bi和bij分别为第1水平和第2水平的随机参数向量；D1、D2分别为第1水平和第2水平的随机参数方差-协方差矩阵；εij为随机误差矩阵；Rij为组内误差方差-协方差矩阵。

建立混合效应模型时，需要确定组内的方差-协方差结构(Rij)，其基本表达式如下：

Rij=σ2G0.5ijΓijG0.5ij。

(3)

式中：σ2为误差离散的比例因子，由估计模型的残差值得出；Gij为解释树内异方差的对角矩阵；Γij表示解释树内自相关结构的矩阵。由于干形数据来自同一株树的重复测量，这类数据中往往存在异方差和自相关问题[20-21]。因而在对模型进行拟合的过程中，必须考虑这两方面问题。结合本研究，最好的解释模型方差异质性的方差结构可以表示为Var(εij)=σ2exp(δXij)，σ2为误差离散的比例因子，δ为待估参数，Xij为模型变量，即Xij=(1-qij1/3)/[1-(1.3/Hij)1/3]；一阶连续自回归相关结构[CAR(1)]能够很好地模拟样木树内的空间自相关，表示为Cov(εk,εl)=ρhkl,Cov(εk,εl)为1棵树内两处不同位置的直径估计值的残差εk、εl的协方差；ρ为CAR(1)的估计参数；hkl=|hk-hl|，其中k≠l，表示在同一棵样木内不同的观测高的差值的绝对值。本研究混合模型的参数估计均在R软件的nlme包下采用限制极大似然法进行。

1.5 削度方程异方差消除

1.6 模型验证

(4)

1.7 模型拟合和检验指标

(5)

(6)

(7)

(8)

1.8 预测精度的估测

在实际生产中，树干上部的数据往往不易获得，这也导致了很难使用树干上部直径进行随机参数的计算，本研究结合数据获取的样本对不同相对高度组合的抽样方案进行比较，抽样方案为：方案Ⅰ，不考虑树干上部直径的测量难度，即不限定相对高的条件下，抽取任意位置的直径及其组合，当模型精度稳定时，确定此时二次抽样的数量和位置组合； 方案Ⅱ，将测量位置限制在人工可测量的高度下(相对高0.1以下)，从每棵树中抽取1～6个样本，找出每个抽样数量对应的最优组合，分析比较不同抽样数量之间预测精度的变化。

2 结果与分析

2.1 包含树冠变量基础模型的建立

以Kozak可变指数削度方程为基础进行再参数化，通过最小二乘法估计模型的参数，并对参数与冠长、冠幅、冠长率之间的关系进行相关性分析以及比较变量加入模型后相应参数的显著性、模型的稳定性和拟合效果，研究发现模型与冠长率有密切关系。将冠长率(CR)引入模型后，使用指数函数对其进行加权来消除异方差，模型的具体形式如下：

Xa7H(1-q1/3)+a8x+ε。

(9)

式中：CR为冠长率。

该模型参数均在95%置信水平上显著，且模型的拟合效果相较于原始模型也有一定的提升，因此引入冠长率的Kozak模型(9)，被选为本研究的基础模型。

为了更直观地反映冠长率对树木干形的影响，选取数据中3棵不同大小的样木(D=27.0 cm，H=21.5 m；D=14.5 cm，H=14.6 m；D=5.4 cm，H=7.25 m)，在不同冠长率(0.20，0.45，0.97)的条件下对干形进行模拟，见图1。

图1 不同冠长率的树干曲线模拟结果

从图1可以看出，对于不同大小的树木，冠长率对直径的预测均能产生一定的影响，冠长率的增大会导致树干削度的增大。对于同一棵树，从图1中可以看出树干根部削度是最大的，随着树干高度的增大削度逐渐减小，干形曲线趋向于一条向中心线微凹的曲线。

2.2 混合效应模型的建立

在基础模型(9)的基础上，建立考虑样地、样木效应的混合效应模型(10)见表2。经过对所有不同随机参数组合的模型进行拟合，对于模型(10)在随机参数组合(样木效应a2i、a6i、a8i和样地效应a2i，j、a6i，j)时取得最好的拟合效果(AIC为5 864.76，BIC为5 978.42，LogLik为-2 913.4)。表2给出了模型(9)和模型(10)的似然比检验结果，根据AIC、BIC和LogLik等指标，混合效应模型(10)优于基础模型(9)，且似然比检验的显著性(P<0.01)，说明考虑样地、样木两个水平的随机效应显著提高了模型的拟合精度。

表2 模型(9)和模型(10)的拟合统计量及LRT检验指标

考虑样地、样木效应的双水平混合效应模型并不能够消除模型的异方差和自相关。本研究采用指数函数来消除异方差，采用一阶连续自回归结构[CAR(1)]来描述树木内误差的相关性。从表3可以看出添加CAR(1)结构的模型(11)比无结构的双水平混合效应模型(10)拟合效果更优，AIC降低了278.21，BIC降低了272.24；在模型(11)的基础上添加指数方差函数的模型(12)，相较于模型(10)AIC降低了438.38，BIC降低了426.43，取得更优的拟合效果；且似然比检验的结果表明CAR(1)结构和指数方差函数能显著提高混合效应模型的拟合效果，因此模型(12)为本研究最佳的混合效应模型。

表3 不同方差-协方差结构下的模型拟合统计量及LRT检验指标

2.3 模型拟合结果

表4 模型(9)—(12)的参数估计及拟合指标统计表

基础模型(9)和最优混合效应模型(12)的残差见图2，残差图显示采用权函数加权的基础模型残差分布较为均匀，没有明显的异方差性；而最优混合效应模型(12)残差的分布更为聚集，表现出比基础模型(9)更好的拟合效果。

图2 基础模型(9)和最优混合效应模型(12)的预测值-标准化残差分布图

2.4 模型检验和预测

2.4.1 模型的检验

模型的检验采用全部数据(混合效应模型的二次抽样采用建模时的全部数据)，使用留一交叉验证法对模型进行检验。模型(9)—(12)模型检验的结果见表5，可以看出检验结果和拟合结果相差不大：混合效应模型(10)—(12)均优于基础模型(9)；最优的混合效应模型(12)的检验精度要高于基础模型(9)的精度(MAE降低了20.15%、MAPE降低了15.41%)。

表5 模型(9)—(12)的检验指标统计表

2.4.2 混合效应模型二次抽样设计对模型精度的影响

混合效应模型的预测采用二次抽样计算随机效应参数，基础模型(9)的预测精度及混合效应模型(12)的两种二次抽样方案(方案Ⅰ、方案Ⅱ)的预测精度的差异见表6。

表6 基础模型(9)和最优的混合效应模型(12)的全局最优抽样方案以及相对高0.1以下的抽样方案的检验指标统计

1)方案Ⅰ：在不限定相对高的条件下，随着抽样数量的增加，模型预测的精度不断增加，当抽样数量达5时，检验精度提升趋于稳定，此时需额外测量相对高0、0.10、0.50、0.60、0.75处的直径，包含了树干下部和上部的拐点，对于误差的解释最为充分，取得最好的检验精度(MAE为0.470 0,MAPE为4.62%)，相对于基础模型MAE降低了17.9%、MAPE降低了13.2%。

2)方案Ⅱ：随着树干下部直径抽样数量的增加预测的精度不断提升，当抽样数量达4时提升幅度趋于稳定，相对于基础模型取得了一定的提升(MAE降低了9.8%,MAPE降低了3.5%)，而抽样数量达6时，则出现精度的下降。同时可以观察到当抽样数量从1增加到6的过程中，统计指标差别不大(MAE的范围0.520 3～0.536 6，MAPE的范围5.14%～5.22%)。在考虑到测量时间和成本的情况下，方案Ⅱ在抽样数量为1时(额外测量树干基部直径)也能够获得较好的预估精度。

由于树干直径是从高至低逐渐增大的，为了进一步分析抽样方案对树干各相对高处的预测精度的影响，本研究选择了基础模型(9)的基础方案和最优的混合效应模型(12)的方案Ⅰ-5(抽样数量为5)、方案Ⅱ-1(抽样数量为1)和方案Ⅱ-4(抽样数量为4)进行各相对高处的分段统计(图3)。

图3 基础模型(9)和最优的混合效应模型(12)的方案Ⅰ-5、方案Ⅱ-1和方案Ⅱ-4在各相对高处的评价指标柱状图

从图3可以看出，随着相对高的增大，模型预测的绝对误差百分比逐渐增大，在相对高为0.90时误差最大。对比不同的抽样方案，模型的预测能力在各分段的表现并不相同。

方案Ⅰ在抽样数量为5时，包含了树干上部和下部的拐点，对树干直径的预测，整体上取得了较好的预测效果，各相对高处的精度均优于基础模型。方案Ⅱ在抽样数量为1时，在树干下部(相对高为0和0.02处)，表现出比基础模型更优的预测精度，其他相对高位置的预测精度并没有出现明显的提升；方案Ⅱ在抽样数量为4时，模型精度的提升主要集中在树干下部(相对高0～0.2)，对于树干上部的预测精度提升并不明显。

3 讨 论

在Kozak模型的基础上，通过再参数化的方法将冠长率变量引入模型，构建的模型参数在5%水平上均是显著的，且再参数化的模型精度比基础模型高，这表明了树冠对干形有很大的影响。树冠不仅为树干的生长和发育提供物质和能量基础，同时树冠内枝和叶的排列决定了树干的形状。冠长率是树冠结构中的重要特征因子，不仅反映了树木的竞争能力，还意味着树木获取太阳辐射能力的强弱，冠长率越大说明树冠的发育越好，同时树干也表现得更为尖削[6]。正如Valentine等[31]在研究树冠对干形的影响时所提出的，发育良好的树冠会导致其内部树干的削度变化加剧。早期的研究表明，向削度方程中加入树冠特征变量不仅能够改善模型的拟合效果，还使得模型构建更具生物统计意义[10]。Ramazan等[32]在对黑松的研究中，认为冠长率在描述树冠与干形的关系时更具有代表性，并将冠长率成功应用于削度方程中。虽然Hann等[33]认为树冠特征变量只解释了很小一部分的树干变化规律，但是将树冠特征变量纳入削度方程还是很有必要的。本研究建立的模型在对树干进行模拟时，树干上部的直径随着冠长率的增大而减小，这与以往的研究结果一致[6-7]，造成这种现象的原因可能是因为较大的树冠导致树干上部营养分配减少[34]，从而导致上部直径的减小，同时模拟的树干曲线也随着冠长率的增大表现得更为尖削，更符合实际林分中干形的变化规律。

本研究采用非线性混合效应模型的方法，考虑了削度方程中普遍存在的自相关和异方差问题，建立了高精度的削度方程。在对混合效应模型预测能力进行验证时，结合了实际应用中数据获取的难度，比较了两种二次抽样方案对削度方程的预估精度的影响。第1种抽样方案(方案Ⅰ)是在不考虑测量成本的条件下，当抽样数量为5时(测量位置为相对高0、0.10、0.50、0.60、0.75处)，混合效应模型的预估精度的提升趋于稳定，此时校正直径的位置包含了树干上部和下部两个拐点，模型表现出较高的预测精度。第2种抽样方案(方案Ⅱ)为将相对高限制在0.1以下，这个高度范围一般是调查人员能直接测量的位置，抽样数量达4时模型预测效果提升幅度趋于稳定，而抽样数量达6时，出现精度下降的现象，造成这种情况的原因可能是使用较多的相对高0.1以下位置处的直径校准的随机效应参数，反而使树干相对高0.1以上处位置的预测偏差更大。在考虑二次抽样样本大小方面，方案Ⅱ随着样本量的变化保持轻微的减小，差异并不明显，如果从减小测量工作量的角度出发，选择抽样数量为1时，如额外测量树干的基底直径，可以提高模型的预测精度。此外，方案Ⅰ的精度普遍高于方案Ⅱ的精度，表明树干上部直径校准可以有效提高模型的精度，这与Manuel等[35]提出的观点是一致的。在以往的研究中[15-16,36-37]，通常采用上部直径对削度方程进行校准，但上部直径在人工测量的条件下往往难以获得，因此在考虑测量时间和成本的前提下，方案Ⅱ更具有实用性。