基于分位数回归的杉木人工林地位级划分方法研究

张博，陈科屹，周来，Sajjad Saeed，张雅馨，孙玉军＊

(1.北京林业大学森林资源和环境管理国家林业和草原局重点实验室，北京 100083；2.中国林业科学研究院林业科技信息研究所，

北京 100091；3.中国林业科学研究院森林生态环境与保护研究所，北京 100091)

森林立地质量是林地更新、树种选择、地力维持、生产力评估和经营管理等林业工作和研究的基础[1-2]。作为评价林分生长类型和林地生产力的重要依据，立地质量评价对研究森林生长收获规律、预估林地生产力和制订相应营林措施具有重要的指导意义[3-5]。目前，利用林分生长量的数据来定量评价立地质量的方法包括：地位级法、地位指数法和立地形法等[1,6]。3种方法分别根据林分平均树高与林分平均年龄、优势木平均高与标准年龄以及优势木平均高与基准胸径的关系来反映立地条件对树高生长的影响。在森林经理经营调查中，由于林分平均高和林分平均年龄是必备调查因子，地位级法被国内外学者广泛关注[7-9]。

传统的地位级表编制方法基于生长模型、误差调整和简单的图形解释虽然能客观反映总体的立地水平及其差异，但采用均值回归模型和标准差调整法编制地位级表对建模数据本身不具有描述性，林分地位级的估测结果可能出现不同程度的高估或低估现象，或者现实林地的生长状况可能根本达不到相应的林分条件平均高[10]。因此，本研究提出一种基于分位数回归模型（Quantile Regression Model）构建地位级表的方法。分位数回归是一种估计因变量或特定分位数函数的完全条件分布的方法，由Koenker等人[11]提出，是统计学和计量经济学常用的回归分析方法之一[12]。一般的均值回归估计方法仅针对协变量的条件均值或中心效应[13]，而分位数回归方法可灵活地描述相应条件分位数自变量和因变量之间关系并得到的回归曲线。目前，分位数回归模型已被应用于多个林业研究中，如：森林资源量的估算[14]，林分密度[15]，直径分布预测[16]，直径增长[17]，削度[18]和森林的地上生物量[19]等。

根据分位数模型可得出各分位数对应的条件估计值和不易受到极端值影响的特点，该模型可用于描述树高生长过程与立地质量的关系[20-21]，但目前该模型在立地质量评价研究中鲜有报道。基于此，本研究以福建省三明市将乐国有林场的杉木（Cunninghamia lanceolata(Lamb.)Hook.）纯林为例，构建基于分位数回归模型的地位级表，并与标准差调整法进行比较，对将乐县国有林场立地质量进行评价与分析，为进一步提高地位级分级策略的效率和立地质量评价的准确性提供理论依据和参考。

1 数据与研究区概况

数据来源于福建省三明市将乐国有林场2012 至2017年期间调查的418 块杉木纯林小班调查样地数据，主要分布于将乐县南口乡、完全乡、黄潭镇、白莲镇、水南镇、余坊乡、光明乡和古镛镇。样地均为面积0.06 hm2的方形样地。主要调查内容包括各样地地理坐标、坡度、坡向、坡位、海拔等地形因子，土壤类型、土壤厚度、腐殖质层厚度等土壤因子，每木检尺测定样地内每株树木的胸径、树高、冠幅等，通过计算得出各样地的林分平均年龄、林分平均胸径、林分平均树高等林分因子。样地各龄级的信息统计见表1。

表1 各龄级样地基本信息Table 1 Summary of basic information statistics of sample plots for each age-class

2 研究方法

2.1 基础模型选择

基于已有的研究成果选择了7个常用的生长模型（表2）作为基础模型[22-25]。

为了比较模型的拟合优度，用赤池信息准则（AIC）、贝叶斯信息准则（BIC）、对数似然值（logLik）、均方根误差（RMSE）、决定系数（R2）和平均绝对误差（MAE），选择拟合模型，所有计算均使用R（Version 3.6.1）软件进行。

2.2 地位级表的编制（标准差调整法）

以导向曲线为基础，按标准年龄时树高值和地位级距（C），采用标准差调整法，可形成地位曲线簇（即树高生长曲线簇）[1]。杉木在25 a 左右树高生长区域稳定，且杉木在20～30 a 时达到数量、经济成熟龄[26]。因此，本研究以林分树高生长量趋于稳定、杉木达到成熟龄确定基准年龄（A0）为25 a。

（1）拟合各龄级树高标准差方程

根据各龄级树高标准差（SH）与龄级平均年龄（Ai），利用SH=a+b×lg（Ai）式拟合龄级树高标准差方程。将各龄级代入，计算出各龄级树高标准差理论值（SA）。本研究方程为：

（2）导算地位级表

通常在基准年龄（A0）时，由于导向曲线的理论树高值可能不是地位级数值，因此需要根据基准年龄时的树高（H0）与标准差理论值（S0）的大小进行调整，公式如下：

式中，Hij为第i 龄级第j 地位级调整后的树高；Hik为第i 龄级的导向曲线树高；H0j为基准年龄时第j 地位级的树高；H0k为基准年龄时导向曲线树高；为基准年龄所在龄级树高标准差理论值；为第i 龄级树高标准差理论值。

以调整后的导向曲线为准，按地位级距C 逐龄级导算出各地位级曲线上的树高值，其余地位级的调整系数Kj 为：

2.3 分位数回归模型

分位数回归模型基于表2 中的线性和非线性模型来预测第τ 分位数的树高模型：

与最小二乘方法相比，分位数回归模型的参数通过最小化分位回归领域的损失函数（或称为检验函数）获得[27]。

在本研究中，分位数回归模型的构建通过使用R 语言中的“quantreg”包来完成[28]。

2.4 地位级分级和评价

计算每个样地林分平均高与各分级曲线树高预测值的差值平方和（或差值的绝对值），以确定每个样地的地位级，具有最小残差平方和的模型即为该样地所属的立地类型，从而确定该林分目前的地位级。分别采用传统方法与分位数回归模型方法统计出每个林分的地位级，对分级结果进行比较和差异分析。

3 结果与分析

3.1 导向曲线的拟合

根据数据资料整理，采用最小二乘法和非线性拟合技术拟合7个基础模型，其拟合结果如表3。依据AIC、BIC、RMSE 和MAE 最小，logLik 和R2值最大的原则，选出最优模型Mod.4（Logistic）为最优导向曲线方程：

表3 导向曲线模型拟合结果汇总Table3 Fitting statistics for Guide curve growth models

3.2 地位级表

（1）基准年龄及地位级距

根据基准年龄确定条件，本研究中杉木的标准年龄为25 a。杉木达到基准年龄时基准树高H0为14.24 m，树高变动范围为6.8～19.3 m。根据将乐地区杉木的编表资料，以及树高、胸径的绝对变动幅度和经营水平，确定地位级距C 为2 m，即地位级H0j分别为6 至20 的8个地位级。

（2）地位级表的编制

标准年龄A0（25 a）代入导向曲线方程得到树高理论值H0k（14.24 m）并计算调整系数Kj和各相应龄级树高值绘制地位级表（表4）。根据立地级表和导向曲线方程拟合8个地位级的生长趋势模型（图1），模型参数及统计量如表5。

表4 杉木人工林地位级及相应树高Table 4 Site class and corresponding tree height of Chinese fir

表5 地位级模型参数及统计量Table 5 Parameter estimation and fitting statistics for site class models

图1 传统方法的地位级分布和建模数据散点分布Fig.1 Distribution of site classes by traditional method and scatter plot of modeling data

3.3 分位数回归模型

地位等级通常分成5～7 级，为便于与传统地位级表编制方法比较，本研究分位数回归模型的地位级表也分为8个地位级。根据数据的分布和导向曲线方程（Logistic 模型）选择了8个分位数（0.01、0.05、0.15、0.30、0.70、0.85、0.95、0.99）（图2）。其中，分位数0.01 和0.99 接近于地位级的下限和上限，可作为地位级Ⅷ和Ⅰ，分位数0.05 和0.95是研究区地位级的最前5%和最后5%水平，分别记作地位级Ⅶ和Ⅱ，分位数0.15、0.30、0.70、0.85则依据数据的分布状况以及保持分位数回归曲线簇形状的相对均匀，分别记作地位级Ⅵ、Ⅴ、Ⅳ和Ⅲ，模型参数统计结果（表6）。

表6 分位数回归模型参数及统计量Table 6 Parameter estimation and fitting statistics for quantile regression models of site class

图2 分位数回归方法的地位级分布和建模数据散点分布Fig.2 Distribution of site classes by quantile regression method and scatter plot of modeling data

3.4 两种分级方法结果比较

从两种方法的地位级分级过程看，传统方法地位级分级曲线规则，能够反映出林分地位级的普遍规律，得到的地位级表示基准年龄时的林分平均高。而分位数回归的8个分位数点的选择是根据建模数据分布来确定的，地位级分级结果反映出的是不同特定条件下的变化规律。因此，拟合的分级曲线可能存在过于接近的现象，基准年龄较小时无法直接判断地位级分级状况，造成对地位级分级结果模糊或对普遍规律的理解不足的现象。

利用各模型的最小残差平方和来判断418 块样地的地位级状况，从两种方法的分级结果（表7）可以看出，传统地位级表中大多数样地地位级分布于10 至16 地位级，与分位数回归模型中多数样地分级分布于Ⅵ至Ⅲ地位级结果基本一致。对两种分级结果进行差异显著性检验，分位数回归模型的地位级分级效果与传统方法没有显著差异（F值0.130 3，P值0.719 3），表明分位数回归模型能够较为准确的描述研究区内林分的地位级的分布特征。两种方法分级的主要差异在于对最大和最小地位级的分级表现，分位数回归方法将更多的样地划分到最大和最小地位级中，反映了研究区各龄级林分平均高的完整分布状况。总体而言，传统方法的拟合曲线表现出的是平均状态和预设的变化规律，而分位数回归模型利用数据分布特征所反映出的是特定状态和条件对应的变化规律，其划分的地位级结果在很大程度上取决于建模样本的结构和数据质量。

表7 传统方法（左）和分位数回归模型（右）对样地地位级分级的结果Table 7 The result of site classes grouped by traditional method (left) and quantile regression model (right)

4 结论

本研究提出一种基于分位数回归模型的立地质量分级和评价方法，以福建省三明市将乐县国有林场小班调查样地数据构建基于分位数回归模型的地位级分级模型，并与传统地位级划分方法相比较。结果表明，基于不同分位数拟合曲线簇，可准确评估林地地位级水平；根据林分平均高与各分级曲线树高预测值的差值平方和（或差值绝对值）最小的原则，可迅速确定该林分的地位级和生长类型，实现快速、准确的立地质量评价，为精准评估森林生产力和进一步提升森林质量评价效率提供方法和依据。

我国森林资源清查数据中主要记录平均木的树高，采用林分平均树高构建地位级表可以使分位数回归方法与森林清查数据兼容。资源清查的小班调查数据虽然没有标准地调查数据准确，但其数据量大且覆盖广，能够从一定程度上反映立地质量对林分生产力的影响[7]。分位数回归模型可描述、分类、预测和验证小班数据林分生长与地位级之间的关联性，基于导向生长模型的分位数回归曲线簇能够更直观的反映出不同地位级下杉木树高的变化轨迹，从整体上描述了各小班林地生产力，并有效地简化了地位级划分的测算工作。因此，结合分位数回归方法和森林资源连续清查数据构建林分尺度和区域尺度的地位级分级体系，可以实现对数据的充分利用。但是，由于分位数回归模型受数据分布的影响，地位级的划分结果很大程度上取决于建模样本的结构和数据质量。今后的研究中需要考虑地位级划分和分位点选择的联系，以进一步提高分位数回归模型在生产力评价中的适用性。