基于分位数组合的杉木树高-胸径模型

余昆隆，谭 伟，杨 靖，王贵林，蒲秀青，姜仕昆

（贵州大学 a.林学院；b.林业信息工程研究中心；c.电气工程学院，贵州 贵阳 550025）

树高和胸径常被作为估测林分生产力和树木材积的重要因子之一[1-3]。胸径测量简单、准确、方便，而树高测量复杂、耗时、费力。利用胸径与树高之间的关系来估计树高已经成为林业生产实践中一种快速且准确的方法[4-5]。关于树高—胸径模型，在林业方面主要是通过非线性和混合效应模型计算参数[6]。非线性模型的参数计算主要是用最小二乘法（OLS），这种方法是基于均值回归，只能获得树高的平均值，但在计算具有多尺度的树高数据时会降低模型的精度和可靠性。而且，OLS 是基于符合独立正态分布假设的参数，在实际调查时，模型中参数的联合分布很遵循独立正态分布[7-8]。混合效应模型不仅可以通过固定效应参数反映群体的平均效应，还可以通过随机效应参数反映群体之间的个体差异，与最小二乘法相比，它能够更灵活地分析数据[9-10]。然而，单木树高的数据通常是由样地中的每木检尺，通常呈现层次嵌套结构[11]，这导致了模型出现许多问题，如异方差、自相关等。虽然加权回归和自相关函数可以改进模型，但模型的参数计算更加困难，甚至有时候效果还不佳。

Koenker 和Bassett[12]在1978 年在研究计量经济学模型时，提出了一种回归方法，即分位数回归（Quantile regression）。其特点是可以更全面地描述总体的分布，而且对残差的分布没有具体的要求，这使得参数估计更加灵活。同时，还可以评估不同分位点对变量的响应[13-15]。近年来，分位数回归在林业领域得到了广泛的应用，赵梦草[16]基于分位数回归，建立了关于树高-胸径、树高-活枝高之间的数学模型。孙拥康等[17]基于分位数回归的方法构建马尾松青冈栎混交林树高-胸径模型，并对比不同分位数模型与传统非线性回归模型的预测精度。王君杰等[18]验证了分位数回归比哑变量模型对于树高有更好的预测效果。他们对分位数回归的应用主要是拟合最佳分位数，然后将最佳分位数回归的参数估计代入模型进行预测。尽管有学者采用多分位数组合对模型进行预测，如马岩岩、Cao 等[19-20]基于树干削度方程的分位数回归方法，通过不同的分位数和分位数组合的方法，取得了很好的预测效果。但是都只是抽取一个数据对模型预测效果的影响，存在一定的局限性。后来 Özçelik 等[21]提出了抽取多个数据对模型进行分位数组合来预测树高，但是只局限于一种抽样方法，并没有深入研究不同的抽样方法对预测结果的影响。

杉木Cunninghamia lanceolata是贵州省的主要树种，具有材质好、成林快、耐腐蚀性强等优点[22]。而树高作为林分生长中的重要因子之一，因此，使用准确可靠的树高—胸径模型对杉木的经营管理具有重要的意义。目前，还没有报道使用分位数组合和不同的抽取方法来预测杉木的树高。本研究主要采用分位数组合来预测树高，将其与非线性回归、分位数回归进行比较，分析各模型之间的性能，并对模型在实际应用时进行评价，以便为今后的实际应用提供指导。

1 材料与方法

1.1 研究区概况

贵州省清镇市国有林场（26°21′00″～26°59′09″N, 106°07′06″～106°33′00″E）位于长江水系乌江上游南侧的清镇市境内，属亚热带季风湿润气候。年平均气温16.5℃，其中最高气温37.5℃，最低气温-7.8℃，海拔为1 143～1 390 m。土壤为黄壤、黄棕壤等，研究区的主要乔木树种有杉木、马尾松Pinus massoniana、枫香Liquidambar formosana等。

1.2 数据采集与预处理

研究数据来自2021 年在清镇市国有林场调查的杉木人工林49 块圆形标准地，单块面积为0.067 hm2。进行每木检尺，起测径阶为5 cm，记录胸径和树高，剔除异常值后，共获得3 795 株杉木数据。在每个样地中计算其胸高断面积和林分密度，并选择6 株优势木的平均值作为优势木平均高。随机选择39 块标准地，计3 054(80%)株杉木作为建模数据，剩余的10 块标准地，计741(20%)株杉木作为检验数据。数据统计量见表1～2。

表1 样地调查因子建模数据Table 1 Fitting data of survey factors in the sample plots

表2 样地调查因子检验数据Table 2 Validation data of survey factors in the sample plots

1.3 研究方法

1.3.1 基础模型选择

基于调查数据，采用最小二乘法拟合了7 个广泛使用且高精度的树高-胸径模型（表3，H为树高，D为胸径，a1、a2和a3为模型参数），并对拟合方程进行评价和比较，根据模型的拟合精度，选择出最优的基础模型。

表3 树高－胸径基础模型选择Table 3 Basic model selection of H-D relation

1.3.2 广义模型选择

在现实林分中，树高除了胸径的影响外，还受到其他单木因子和林分变量的影响。为了减少林木水平或林分水平差异对树高的影响，从优势木平均高、胸高断面积和林分密度中引入能够提升模型的预测效果的变量作为协变量。在上述的最优基础模型中，加入2 个协变量构建树高—胸径广义模型。

1.3.3 分位数回归

分位数回归可以对数据的任意分位点进行响应估计，分位数一般用τ表示。其参数估计的目的是使损失函数（1）的目标值（Lmin）达到最小。

式中：τ分位点数值，y与分别为τ点的真实值和估计值。

τ介于0 与1 之间。本研究基于最优基础和广义模型，分别选取9 个分位点（τ=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9）对树高进行预测，当分位数τ=0.5 时的分位数回归定为中位数回归。

1.3.4 分位数组合应用方法

Bohora 等[23]于2014 年首次提出分位数组合方法，该方法主要是对检验数据进行抽样处理，根据抽样数据找到最近的两条分位数曲线进行插补，最后得到改进后的分位数曲线。本研究主要探究三分位数组合（τ=0.1,0.5,0.9）、五分位数组合（τ=0.1,0.3,0.5,0.7,0.9）、九分位数组合（τ=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9）下树高的预测效果。

当每个样地内抽样数量为1 时，找到最接近该数据的两条分位数回归曲线，即：(Hij)≤Hij≤最后通过插值改进后得到1 条树高曲线，其树高的估值计算为：

1.3.5 分位数组合抽样设计

本研究主要设计了4 种抽样方案，使用方案中的抽样数据对模型进行修正，并对未抽样的数据进行检验。具体如下：

1）每个样地中抽取胸径最大的树1～9 株；

2）每个样地中抽取胸径最小的树1～9 株；

3）每个样地中抽取最接近平均胸径的树1～9 株；

4）每个样地中随机抽取1～9 株。

1.4 模型的评价与检验

本研究选用3个常用的指标，即决定系数（R2）、均方根误差（RMSE）、平均绝对误差（MAE）对模型进行验证。R2指标值越接近1，则表明模型拟合效果越好；RMSE 和MAE 指标值越标值越接近0，则表示模型效果越好。3个指标的表达式如下:

式中：Hij与分别为第i个样地中第j株树高的实测值和预测值（m）；为树高实测值的均值（m）；m为样地的数量；n为样本的数量。

1.5 数据处理

所有数据分析均通过R 统计软件完成，其中，分位数回归由quantreg 软件包拟合，并使用Origin软件绘图。

2 结果与分析

2.1 基础模型选择

根据建模实测数据，利用R 软件拟合了7 个树高-胸径模型，根据R2、RMSE 和MAE 三个统计指标，筛选出能够最好体现树高与胸径关系的基础模型。基础模型的拟合结果如表4，这7 个基础模型的R2在0.553 3～0.590 4 之间，RMSE 在2.112 1～2.223 3 之 间，MAE 在1.710 5～1.791 5 之间，都能很好反映树高和胸径之间的关系。以Richards（M1）基础模型的拟合结果最好，其R2、RMSE 和MAE 分别为0.590 4、2.112 1 和1.710 5。所以本研究选择Richards 作为最优基础模型，其一般表达式为:

表4 基础模型拟合结果Table 4 Fitting results of basic models

式中：a1、a2和a3为模型参数。

2.2 广义模型拟合

故在Richards（M1）基础模型上，选取2 个变量进行不同组合的检验。经过计算，加入优势木平均高和胸高断面积的广义模型拟合效果最好，与基础模型（M1）相比，R2上升了0.115，RMSE和MAE 分别降低了0.344 3、0.310 5，进一步促进模型的拟合结果。综上，推导出最优广义模型一般表达式为：

式中：H0,i为第i个样地的优势木平均高（m）；BA,i为第i个样地的胸高断面积（m²·hm-2）；a1、a2、a3、a4和a5为模型参数。

2.3 分位数模型的参数估计

本研究具体采用9 个不同的分位点（τ=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9），使用R 统软件中的quantreg 软件包分别对基础和广义模型进行拟合。由表5 可以看出，基础胸径-树高模型在9 个分位点都能收敛且有一定的规律，其中参数a1和a2逐渐增大。可能是由于广义模型的参数较多，9 个分位点的参数估计值并没有表现出一定的规律。在基础和广义模型中MAE 和RMSE 先减少然后又增加，当τ=0.5 时，模型的拟合效果最好，可以作为最优分位数模型，而且τ离中位数（0.5）越远，模型的拟合效果越差。基础模型的中位数回归与非线性回归相比，MAE、RMSE分别降低了0.001 3、0.003 8；广义模型的中位数回归与非线性回归相比，MAE、RMSE 分别降低了0.022 7、0.003 7。

表5 非线性和分位数回归的参数估计Table 5 Parameter estimation for nonlinear and quantile regression

2.4 不同分位数组合预测树高

基于各分位数参数的估计，使用3 种不同的分位数组合和4 种抽样方案对树高进行预测，并计算预测误差RMSE 统计量。为了方便比较分位数组合与最优分位数回归（τ=0.5）的拟合效果，分位数组合的RMSE 用直方图表示，最优分位数回归的RMSE 用一条水平线表示（图2）。由图2可以看出，无论采取哪种抽样方式，3 种分位数组合的RMSE 基本上随着样本数量的增加而降低，而且它们之间的差异很小，且随着抽样样本的增加而模型精度提高。此外，不同的抽样方案也会对模型的预测结果产生显著的影响。其中，只抽大树和小树时，模型的预测精度没有明显地提升，少量的抽样数据的预测效果甚至不如中位数回归。根据RMSE 的下降趋势，最优抽样方案为：对基础模型的分位数组合抽取5 株平均木、对于广义模型的分位数组合抽取7 株平均木。虽然抽取更多的平均木可以进一步降低RMSE，但RMSE 的下降速率较慢，而且会增加更多的采样成本。

图1 基础模型（A）和广义模型（B）的5 个分位点（0.1、0.3、0.5、0.7 和0.9）和OLS 的回归曲线Fig.1 Regression curves for five quartiles (0.1, 0.3, 0.5, 0.7, and 0.9) and OLS for the base models (A) and generalized models (B)

图2 分位数组合模型的预测误差比较Fig.2 Comparison of prediction errors of the quantile groups

2.5 模型验证评价

为了更好地提高模型的预测精度，利用检验数据分别对各模型进行验证。从表6 可以看出，参数估计方法由高到低排序为：分位数组合＞分位数回归＞非线性回归。本研究基于三位数组合和上述的最优抽样方案设计的模型，与非线性回归和分位数回归模型相比，基础和广义模型的MAE 分别下降了0.409 0 和0.141 5，RMSE 分别下降0.379 4 和0.103 5，R2分别提高了0.131 0 和0.031 8。说明分位数组合模型能提升模型拟合能力。

表6 非线性回归、中位数回归和三分位数组合的拟合结果Table 6 Fitting results of nonlinear regression, median quantile regression and three quantiles groups

3 结论与讨论

3.1 讨 论

本研究以贵州省清镇市国有林场49 块样地中的3 795 株杉木为研究对象，利用7 个树高-胸径模型采用最小二乘法拟合，选择Richards 模型（R2=0.590 4，RMSE=2.112 1，MAE=1.710 5）作为本研究的最优基础模型，这与樊伟等[5]的研究结果一致。大量研究表明，树木生长受到立地条件、林分特征以及竞争的影响。本研究将优势木平均高和胸高断面积加入Richards 模型，可以显著提高模型的拟合和预测能力。拟合精度高于基础模型，表明林分变量对树高与胸径关系具有一定的影响，增加林分变量的树高-胸径模型具有更高的拟合和预测效果，这与臧颢、陈浩等[9-10]对树高模型进行研究时，得到了类似的结果。

分位数回归可以更加灵活地反映树高的分布，中位数回归相比于其他几个分位数点能够更准确地对树高进行预测，这与王君杰等[18]利用分位数对树高-胸径模型研究结果相同。本研究基于3 种分位数组合和4 种不同的抽样方案对树高进行预测。研究结果表明，不同的抽样方式对预测效果产生显著的差异，这与王君杰等[24]预测枝下高的结果一致。抽取平均木的抽样方法使分位数组合的预测效果最佳，这或许是由于平均木最具有代表性，能够全面地表达整个样地的信息。对基础和广义的三分位数组合，抽取5 和7 株平均木以上时，没有显著性差别。综合考虑下，本研究采用上述方法，相比于非线性回归,基础和广义模式的MAE 减少了0.409 0 和0.408 9，RMSE 下降了0.379 4 和0.290 4，R2提高了0.131 0 和0.086 2，相较于最优分位数回归，MAE 下降了0.409 5 和0.386 2，RMSE 下降了0.370 9 和0.286 7，R2提高了0.130 4 和0.086 2。结果表明，分位数组合通过少量的抽样数据能够产生较好的预测精度，具有很强的适用性。由于林分结构和生长动态的复杂性和多变性，如果只考虑某些因子如优势木平均高、胸高断面积等对树高-胸径模型的影响，而忽略了如气候、降水、立地条件等的影响，会导致模型精度不高、通用性不强等问题。而且本研究区域仅在一个林场，导致研究数据具有局限性，其他区域的杉木树高预测需要进一步验证。在未来的研究中，可以多因素和多区域对模型进行应用，为杉木人工林可持续发展及经营管理提供参考依据。

3.2 结 论

准确可靠的树高预测对于林业生产和实践具有重要的意义。本研究基于Richards 基础和广义模型，通过引入9 个分位数点和3 种分位数组合的4 种抽样方案对模型中的参数进行计算。在三分位数（τ=0.1,0.5,0.9）组合的基础和广义模型中，发现基于平均木的抽样方案的预测效果最好，在分别抽取5 株和7 株平均木时，模型预测精度提高幅度最大。在实际中使用分位数组合模型时，应注意同时既提高了模型预测的精度又要考虑到节约成本，建议使用广义模型三分位数组合、抽取7 株平均木的方法对每块样地的树高进行预测。