树高和胸径是衡量林分结构和质量的最重要的2个指标^[1]。其中，树高数据对估计树木材积、生物量，描述林分生长变化不可或缺^[2-3]。但树高的测量耗时耗力、误差较大^[4]，而胸径的测量相对容易、精确、花费少。故实际操作中，往往先测定少数树木的树高和胸径，建立H-D 的函数关系，再以胸径估计剩余树木的树高^[5-8]。由于立地、林龄和管理措施等因素影响，预测树高时需考虑H-D关系在不同林分间的差异^[9]。目前实现不同林分中树高预测的方法主要包括：1) 采用足够数量的树高实测值构建的H-D模型^[10]；2) 针对具体样地分别建立适用于小面积相似林分的H-D模型^[11]；3) 建立含有表征林分特征的林分变量或树木因子的H-D模型^{[1, 4, 8]}；4) 建立能够反映不同林分(或样地)间差异的混合模型^{[5-7, 12]}。前2种方法在森林规模较小时可以实现，但是在大范围林分中所需工作量太大，较难实现。后2种方法中，广义模型具有较好的生物学解释和较高的精度^{[11, 13-14]}，混合模型在处理纵向数据和重复测量数据方面优势明显，且精度较高^[15]。目前林业调查多使用后2种方法，或结合2种方法建立广义混合模型预测树高^{[6, 14]}。

大量关于H-D模型的研究表明：不同树种的H-D关系模型有明显差异^{[1-2, 5-8]}。多数是固定效应模型中的广义模型精度高于基础模型，混合模型精度高于固定效应模型^{[5, 14]}。但基础混合模型和广义混合模型间差异却随树种不同而变化。例如，对于欧洲云杉Picea abies、欧洲赤松Pinus sylvestris和白桦Betula pubescens是广义混合模型精度高于基础混合模型^[16]，但对于颤杨Populus tremuloides却是基础混合模型预测精度高于广义混合模型^[7]。还有研究发现对于红松Pinus koraiensis^[5]、杉木Cunninghamia lanceolata^[6]和落叶松Larix gmelinii^[17]，2种混合模型间差异不大。沙地樟子松Pinus sylvestris var. mongolica是我国北方地区主要造林树种之一，在防沙治沙和改善环境方面具有重要生态价值。迄今为止，该树种已经推广到我国北方地区的19省、自治区、直辖市^[18]。随着造林面积的扩大，如何合理经营沙地樟子松人工林？不断扩大的沙地樟子松人工林资源的清查工作如何快速高效地进行？解决这2个棘手的问题，必须有精确的沙地樟子松H-D关系模型，而该方面的研究鲜见报道。鉴于此，本文筛选合适的H-D基础模型，检验林分变量中优势木平均高 (H_T)、胸高断面积 (A_B) 和平方平均胸径 (D_QM) 对H-D模型的影响，建立H-D基础混合模型和广义混合模型，并确定计算随机参数时的抽样方案，以期为沙地樟子松林的生长预测和林分管理、资源清查提供科学依据。

1.   材料与方法

1.1   研究区概况

研究区位于辽宁省阜新市章古台镇 (42°39′~42°43′N, 122°22′~122°33′E)。该地区位于科尔沁沙地东南缘，属于半湿润气候，近40年(1973—2015)年平均气温6.99 ℃；年降水量平均483.63 mm，最大713.2 mm，最小285.9 mm，约67.0%集中在6—8月；年均蒸发量1 762 mm。最暖月份平均气温–12.91 ℃，最冷月份平均气温23.58 ℃，大于5 ℃有效积温2 570 ℃，全年无霜期150~160 d。土壤主要为生草固定风沙土和流动风沙土。土壤养分含量比较低，有机质质量分数为0.6~1.6 g/kg，全N质量分数为0.17~0.38 g/kg，P₂O₅质量分数为0.12~0.27 g/kg，pH为6.8~7.0^[18]。植物区系属于蒙古区系西辽河小区，主要代表性植物有色木槭Acer mono、山里红Crataegus pinnatifida var. major、榆树Ulmus pumila、大果榆U. macrocarpa、山杏Armeniaca sibirica、胡枝子Lespedeza bicolor、中华隐子草Cleistogenes chinensis等^[19]。

1.2   数据收集

数据来源于2016年设置于章古台镇35块临时标准地，标准地规模20 m×20 m。所选标准地涵盖了沙地樟子松不同的林龄 (A)、立地条件和林分密度 (N)。标准地的林分年龄为13~62 年，立地覆盖了平坦沙地(11块)及各个坡向的坡上(7块)、坡中(11块)和坡下(6块)条件。测量标准地内沙地樟子松的树高 (H)、胸径 (D)；以每公顷100棵胸径最大沙地樟子松的平均树高为优势木高^{[5, 7]}，按标准地面积比例确定株数。从测量的数据中，选择健康、未受损害且D≥5 cm的沙地樟子松样木数据用于模型建立。本研究共获得1 070组H-D数据，分为2个数据集，一个为由28块标准地中的866棵(80%)沙地樟子松构成的建模数据，另一个为由其余7块标准地的204棵(20%)沙地樟子松组成的检验数据，表1为数据基本概况。

表  1  沙地樟子松标准地的林分与树木指标统计

Table  1.  Summary statistics of stand and tree indices of Pinus sylvestris var. mongolica sample plot

变量1) Variable	总数据 Total data				建模数据 Modeling data				检验数据 Validation data
	$\overline x \pm {\rm{SE}} $	最小值 Min.	最大值 Max.		$\overline x \pm {\rm{SE}} $	最小值 Min.	最大值 Max.		$\overline x \pm {\rm{SE}} $	最小值 Min.	最大值 Max.
H/m	8.5±0.1	2.1	16.1		8.3±0.1	2.1	16.1		9.2±0.2	2.2	14.3
D/cm	14.6±0.2	5.0	29.1		14.3±0.2	15.0	29.1		15.8±0.3	5.1	28.6
Drange/cm	10.9±0.4	5.7	16.5		10.9±0.4	5.7	16.5		11.0±0.9	7.3	14.2
Dmax/cm	22.0±0.9	10.0	29.1		22.0±1.0	10.0	29.1		21.9±2.1	11.9	28.6
Dmin/cm	11.0±0.8	1.9	18.0		11.1±0.9	1.9	18.0		10.9±1.9	3.3	18.0
A/a	36±2	13	62		36±3	13	62		36±6	13	56
N/(trees·hm–2)	842±80	300	2 500		862±96	300	2 500		764±112	450	1 175
HT/m	9.7±0.5	3.2	13.6		9.7±0.6	3.2	13.6		10.0±1.1	4.1	12.7
AB/(m2·hm–2)	15.26±1.04	2.92	33.34		15.16±1.19	4.28	33.34		15.65±2.31	2.92	21.97
DQM/cm	16.7±0.9	7.0	24.0		16.6±1.0	7.0	23.9		16.8±1.9	7.7	21.9
1) Drange、Dmax 和Dmin 分别为样地内胸径的分布范围值、最大值和最小值；A：树龄；N：林分密度；HT：优势木平均高；AB：胸高断面积；DQM：平方平均胸径 　1) Drange, Dmax and Dminare the distribution range value, maximum and minimum of DBH in the plot; A: Age of stand; N: Stand density; HT: Dominant height; AB: Stand basal area; DQM: Quadratic mean diameter

下载: 导出CSV 

| 显示表格

1.3   基础模型

在以往报道的基础上选择了30个非线性H-D候选基础模型^{[8, 20]}，使用最小二乘法分别拟合总数据中沙地樟子松的树高–胸径关系，并根据评价指标：决定系数(R²)、平均绝对误差 (Mean absolute error, MAE)和均方根误差(Root mean square error, RMSE)选出较好的H-D基础模型。R²、MAE和RMSE的计算公式为：

式中：H_ij、 ${\hat H_{ij}}$ 和 ${\bar H_{ij}}$ 分别为树高观测值、预测值和平均观测值；m为标准地数量；n为样本木总量。

1.4   广义模型

基础模型是在总体平均水平上对H-D关系的描述，不能反映H-D关系在不同林分中存在的实际差异^{[1, 17]}。基础模型中加入林分变量构建广义模型可以提升模型的生物学意义和预测精度。本文在确定基础模型后，检验了优势木平均高 (H_T)、胸高断面积 (A_B) 和平方平均胸径 (D_QM) 对H-D模型的影响(加入林分变量1~3个)。

1.5   混合模型

确定基础模型和广义模型后，采用非线性混合模型描述不同标准地内的H-D关系，模型的一般表达式为：

式中：y_i是第i块标准地的树高测量值，i=1，2，3，…，m；β为固定效应参数；u_i为随机效应参数，假定u_i服从Ψ正态分布，即u_i~N(0, Ψ)，Ψ为随机效应方差–协方差矩阵；x_i为自变量测量值；f为包含模型参数的x_i与y_i之间的函数关系；ε_i为误差项，假定ε_i~N(0, R_i)；假定u_i和ε_i之间相互独立；R_i为样地内方差–协方差结构矩阵；σ²为误差扩散比例因子，即残差方差；G_i为描述方差异质性的对角矩阵；Г_i为自相关结构，本文调查数据未涉及误差自相关性，故用单位矩阵代替。

1.6   模型的拟合与检验

固定效应模型和混合模型分别使用R软件中nls函数和nlme程序包中极大似然估计方法实现。拟合优度指标为Akaike信息量准则 (Akaike’s information criterion, AIC)、贝叶斯信息量准则 (Bayesian information criterion, BIC) 和似然函数值 (log-likelihood, LL)。

模型的预测类型包括：固定效应模型的平均水平预测 (FPA)，混合模型的总体平均响应预测 (MPA) 与主体响应预测 (MPS)。其中MPA仅使用混合模型的固定效应部分，预测时令随机参数为0。对未包含在建模数据中的林木进行MPS时，需使用经验线性最优无偏估计方法计算随机效应参数值^[11]。

式中： ${\hat u_i}$ 为随机参数最优无偏估计值， $\hat {\boldsymbol \varPsi} $ 为随机效应方差–协方差矩阵， ${\hat {\boldsymbol R}_i}$ 为标准地内方差–协方差结构矩阵；Z_i为随机参数的偏导数设计矩阵，T表示转置矩阵；e_i为MPA时所得的残差。

通常情况下，用于计算随机参数的样本木数量越多，预测效果越好，但调查成本会提高。为了选出工作量较少且精度较高的样本木抽样方案，本文对比了以下抽样方案：1)每块标准地随机抽取1~8棵树木，测量其树高；2)参照平均胸径，每块标准地抽取1~8棵平均木，测量其树高。

评价模型精度时，基于建模数据计算了不同模型不同预测类型时的R²、MAE和RMSE。使用检验数据计算MAE和RMSE来反映模型预测精度。

2.   结果与分析

2.1   构建的基础模型和广义模型

拟合总数据时，30个基础模型的R²、MAE和RMSE差异很小。依据参数最少、效果最好的原则，选择式(6)的Sibbesen模型^[17]作为樟子松的H-D基础模型。在式(6)基础上，将H_T、A_B和D_QM以不同的组合形式加入式(6)中，式(7)~(9)为分别加入不同数量林分变量时的最优形式，即广义模型。

式中：β₁~β₆为模型参数；H_ij和D_ij分别为第i块标准地中第j棵树的树高和胸径；H_Ti、A_Bi和D_QMi分别为第i块标准地优势木平均高、胸高断面积和林分平方平均胸径。

结果(表2)表明：加入林分变量后，广义模型的拟合优度极显著高于基础模型 (P<0.01)，且模型的精度明显提高，模型的残差方差 (σ²)、AIC和BIC均减小，LL值明显增大。只加入1个林分变量时，加入H_T的式(7)的效果最好，此时模型的统计指标变化较大。加入2个林分变量时，含H_T和A_B的式(8)精度进一步提高，且与式(7)差异极显著(P<0.01)。含3个林分变量的式(9)精度低于式(8)。

表  2  基础和广义模型各参数估计结果及统计量

Table  2.  Parameter estimates and evaluation statistics of basic and generalized models

公式 Formula	β1	β2	β3	β4	β5	β6	σ2	AIC	BIC	LL
(6)	14.242 3	–27.039 5	1.784 9				3.305 1	3 498.3	3 517.4	–1 745.2
(7)	–0.105 1	–5.428 3	1.442 1	1.073 6			0.724 7	2 184.8	2 208.6	–1 087.4
(8)	0.190 2	–4.472 6	1.217 7	1.112 6	0.055 7		0.709 8	2 167.7	2 196.3	–1 077.9
(9)	0.404 3	–4.646	1.222 3	1.177 7	0.045 9	–0.034 5	0.709 3	2 168.2	2 201.5	–1 077.1

下载: 导出CSV 

| 显示表格

2.2   构建的混合模型

拟合基础混合模型时，只含1个随机参数的模型拟合优度最高。拟合广义混合模型时，模型含有2个随机参数时的拟合优度略高于含1个随机参数的模型，但无显著差异 (P>0.05)。因此，4个混合模型中均只含1个随机参数，具体见式(10)~(13)。

式中：u_i为随机效应参数；其余同上。图1为式(10)的残差分布图[式(11)~(13)与此相似而未列出]，从图中未发现模型存在明显的方差异质性，因此不考虑对混合模型进行异方差处理。式(10)~(13)的AIC、BIC及LL的变化规律与4个固定模型相一致。式(10)的随机参数方差 (σ_i²) 在4个混合模型中最大，说明式(10)中参数β₁最不稳定，在不同标准地间变化最大。

图  1  混合模型残差分布图

Figure  1.  Distribution of residuals for mixed model

下载: 全尺寸图片 幻灯片

由表2和表3可知，混合模型的拟合精度均高于相应的固定效应模型。加入随机参数对H-D基础模型改良效果最明显，混合模型式(10)的AIC和BIC较基础模型式(6)分别降低了34.70%和34.38%，LL值升高了34.84%。说明混合模型可以提高H-D模型的预测精度，随机参数可以描述不同立地条件、样地间树高和胸径关系的差异。广义混合模型与对应的广义模型间都存在极显著差异 (P<0.01)，但精度提高的幅度较小(AIC、BIC和LL变化幅度低于1.50%)。

表  3  混合模型参数估计及统计量

Table  3.  Parameter estimates and evaluation statistics of mixed model

公式 Formula	β1	β2	β3	β4	β5	β6	σi2	σ2	AIC	BIC	LL
(10)	10.258 1	–7.630 1	1.565 7				9.617 6	0.682	2 284.3	2 308.2	–1 137.2
(11)	0.211 3	–6.414 4	1.491 4	1.046 9			7.98×10−3	0.681	2 163.3	2 191.9	–1 075.7
(12)	0.248 3	–5.571 8	1.378 5	1.030 0	0.043 0		5.81×10−3	0.680	2 158.2	2 191.6	–1 072.1
(13)	0.375 3	–5.693 8	1.380 3	1.085 9	0.034 9	–0.029 1	5.65×10−3	0.680	2 159.3	2 197.4	–1 071.6

下载: 导出CSV 

| 显示表格

2.3   模型评价与对比

2.3.1   模型拟合与预测精度评价

拟合结果(图2)表明，式(10)的MPS在所有模型中精度最高，但与广义模型的FPA、MPA和MPS差异均不明显。除了式(12)MPA时的MAE略小于对应FPA外，其余指标表明，4对模型精度均为：MPS> FPA> MPA。模型只含预测变量D时，FPA和MPS精度明显较其MPA提高，R²分别升高了39.51%和93.00%，MAE分别降低了27.68%和68.88%，RMSE分别降低了20.69%和64.38%。含林分变量时，各模型的拟合精度较高但差异较小。相较于对应MPA，各模型FPA和MPS时精度变化不明显，FPA时R²、MAE和RMSE改良范围不足1%，MPS时精度变化稍微增大但低于5%。

图  2  基于建模数据计算所得不同模型的不同预测类型指标统计

A：式(6)和(10)；B：式(7)和(11)；C：式(8)和(12)；D：式(9)和(13)

Figure  2.  Indices calculated based on model fitting data for different types of predictions in different models

A: Eq.(6) and (10); B: Eq.(7) and (11); C: Eq.(8) and (12); D: Eq.(9) and (13)

下载: 全尺寸图片 幻灯片

2.3.2   不同混合模型抽样对比

图3为各模型FPA、MPA预测及不同抽样方案下MPS预测时MAE和RMSE的变化。进行FPA和MPA预测时，广义模型精度均高于基础模型，且FPA中式(8)的精度最高。进行MPS时，混合模型的MAE和RMSE变化趋势相同。其中，基础模型精度一直为：MPS>FPA>MPA。使用随机选择样本木计算模型随机参数时的模型预测效果好于使用平均木计算。每块标准地内随机测量3棵树的树高时，式(10)的MAE和RMSE的降低速率最大，此时模型MAE和RMSE较其MPA分别降低了57.97%和57.63%。广义混合模型[式(11)~(13)]的预测精度同样随样本数量增多而降低(1~8棵)，随机选择2个样本木和选择3个平均木时模型MAE和RMSE降低速率最大，但整体变化缓慢。

图  3  不同抽样方案下各模型平均绝对误差(MAE)和均方根误差(RMSE)随样本木数量变化情况

FPA中，A~D依次为式(6)~(9)；MPA及1~8抽样中，A~D依次为式(10)~(13)

Figure  3.  Relationship between mean absolute error (MAE), root mean square error (RMSE) and number of sample trees based on different sampling designs

FPA: A-D are eq.(6)-(9); MPA and different sampling design: A-D are eq. (10)-(13)

下载: 全尺寸图片 幻灯片

3.   讨论与结论

本文以辽宁章古台地区樟子松人工纯林调查数据为例，从30个常用模型中选出了Sibbesen模型为适合樟子松(D>5 cm)的H-D基础模型。以此为基础，建立不同的H-D广义模型和混合模型时发现，广义模型和混合模型都能提高模型精度，且广义模型和基础混合模型精度均较高。

研究表明，立地条件和林分密度是影响H-D关系最为明显的因素^{[1, 7]}。H_T可以代表林分的立地质量和发展阶段^[21]，王冬至等^[13]在研究华北落叶松Larix principis-rupprechtii和白桦H-D模型时的结果表明，包含H_T的模型拟合效果较好。本文的研究也表明，H_T对沙地樟子松H-D模型的影响最大，含H_T的式(7)模型精度明显高于式(6)的基础模型，且模型R²、MAE和RMSE的变化幅度高于30%。A_B可以表示林分密度和林分竞争，它可以反映其他因素(如林龄)对H-D关系的影响^[9]。当模型中含有H_T和A_B时，加入D_QM反而降低了模型精度。此外，H-D关系是一个动态关系，即林龄不同时期H-D关系模型的参数也不同，由于混合模型中随机参数综合反映了包含林龄和立地条件在内的多种差异^[7]，且广义模型中A_B和H_T均与林龄关系密切，故而本文未将树木按照林分阶段划分后分别建模。

模型中无随机参数时(MPA和FPA)，模型的预测精度主要受预测变量的影响^[22]，此时FPA与对应MPA随林分变量个数变化产生的趋势相同。还可能因为混合模型的MPA不能够充分代表其总体水平预测^[23]，而出现FPA优于MPA(尤其基础模型)。林分变量和随机参数都能反映不同林分间差异^[22]，4个混合模型间MPS差异不明显，说明林分变量对混合模型的作用不大^[17]。但并不意味广义模型无构建价值，因为广义模型具有较强的生物学意义且数据量较小时具有明显优势^[22]。

模型预测时，如何抽取样本木使得模型精度较高且工作量较小十分重要。使用经验线性最优无偏估计法计算随机参数时，式(10)在随机选择3棵样本木时，模型的MAE和RMSE降低速率最大，之后降低速率变慢。基础混合模型的MPS精度与广义模型各预测形式间差异不大。由于广义混合模型预测时，同时需要林分变量和先验信息，调查成本较高，但模型预测精度变化较小，不太实用。因此，本文建议使用基础混合模型，并通过随机测量每块样地内3棵树的树高，以此计算随机参数来预测树高，当无先验信息时，使用含有林分变量优势木平均高、胸高断面积的广义模型预测树高。