基于树种分类的阔叶混交林单木枝下高预测模型

Prediction Model of Individual Height to Crown Base in Mixed Broad-leaved Forest Based on Tree Species Classification

摘要

帽儿山阔叶混交林以其复杂的林分结构、丰富的物种多样性以及树种间错综复杂的相互作用对枝下高的生长产生影响。精确构建枝下高模型对于指导林分经营管理、提升生长预测准确性有重要参考价值。选择28块在帽儿山地区有阔叶混交林代表性的样地，并将样地内广泛分布的18个树种按照树种胸径变异系数、树种平均胸径、树种平均高径比、树种占比、软阔硬阔树种分类以及树种耐荫性等指标进行分组，通过该方法解决树种复杂性以及单一树种样本量不足而导致的模型构建的难题。采用维科夫（Wykoff）模型作为基础模型。考虑林分、竞争以及物种多样性等因子表达林分混交度、竞争等情况，建立广义枝下高模型。考虑样地、树种组间枝下高的差异，构建枝下高混合效应模型。采用“留一法”进行检验预测，分析抽样对结果的影响。结果表明帽儿山常见的树种可以分为4个树种组，除胸径、树高外，优势木平均高、种内竞争大于对象木断面积之和以及Shannon指数对枝下高有显著影响。样地和树种组混合效应的基础和广义枝下高模型具有较高的拟合精度，R²为0.638和0.627，RMSE分别为1.880和1.909。每个样地、每个树种组随机抽取1株样本校正混合模型时，广义枝下高混合模型相较于固定效应模型，其MAE和MAPE分别降低了7.51%和13.51%，表现出更好的预测效果。分析帽儿山地区树木种类多样以及树木各生态功能对枝下高生长的影响，为预测帽儿山阔叶混交林不同树种的枝下高提供参考。

Abstract

The mixed broad-leaved forest at Maoer Mountain is characterized by its complex stand structure， rich species diversity， and intricate interactions between tree species， all of which influence height to crown base. Therefore， accurately constructing a height to crown base model can provide important reference value for guiding forest management and improving the accuracy of growth prediction. In this study， 28 representative sample plots of mixed broad-leaved forests in Maoer Mountain area were selected， and 18 tree species widely distributed in the sample plots were divided into groups according to the variability coefficient of tree species diameter， average tree species diameter， average height-diameter ratio， proportion of tree species， classification of soft broad-leaved and hard broad-leaved species， and shade tolerance of tree species， thereby solving the problem of tree species complexity and the insufficiency of sample size of a single tree species in model construction. In this study， Wykoff model was adopted as the basic model. At the same time， factors such as stand， competition and species diversity were considered to express the mixed degree and competition of stand， and a generalized model of height to crown base was established. At the same time， the difference of height to crown base among plots and tree species groups was considered to build a model of subbranch height mixed effect. The influence of sampling on the result was analyzed by using the method of ‘remaining one’. The results showed that the common tree species in Maoer Mountain could be divided into four tree species groups. Besides DBH and tree height， the average height of dominant trees， the sum of the section area of intra-species larger than that of the target trees and the Shannon index had significant effects on height to crown base. Considering the mixed effects of plots and tree species groups， the basic and generalized subbranch height models had high fitting accuracy R² of 0.638 and 0.627， and RMSE of 1.880 and 1.909， respectively. In addition， when one sample was randomly selected from each tree group in each plot to correct the mixed model， compared with the fixed effect model， the MAE and MAPE of the generalized height to crown base model decreased by 7.51% and 13.51%， respectively， showing better prediction effect without wasting manpower， material and financial resources. This study analyzed the effects of various tree species and ecological functions on the growth of height to crown base in Maoer Mountain， and provided some reference for predicting the height to crown base of different tree species in broad-leaved mixed forest in Maoer Mountain.

Graphical abstract

关键词

阔叶混交林 / 枝下高 / 树种分组 / 混合效应模型 / 模型校正

Key words

Mixed broad-leaved forest / height to crown base / grouping of tree species / mixed effect model / calibration of model

0 引言

我国现有森林资源的面积为2.20亿hm²，其中天然林面积约为1.22亿hm²，蓄积122.96亿m³，分别占我国森林面积总量和森林总蓄积量的58%和83%。由此可以看出天然林在我国森林资源中的重要性^［1］。天然林有丰富的物种多样性和稳固的生态系统功能^［2］，然而，其森林类型的复杂性、由多树种组成以及树种间错综复杂的相互作用关系，都为构建精准的生长模型增加了难度^［3］。现有许多研究在解决天然林建模方面已有成效。Temesgen等^［4］利用科或属作为分类标准将样地的林木分类并建立树高-胸径模型，但被证明该方法仅适合生物多样性的应用。对于多层次、多树种的天然混交林生长关系的描述相比纯林要困难和复杂，刘耀凤等^［5］为了描述树种间直径的生长差异，以树种作为哑变量建立单木直径生长模型，成功证实了该方法可以提高模型的预测能力，但是当树种较多时，需要引入大量的哑变量，这可能导致模型变得复杂且难以解释。相比之下，Zhao等^［6］使用层次聚类分析，根据树种直径结构、树木生长、死亡和耐荫性等特征作为分类指标对样地内树种进行分组，并成功应用于美国温带天然林，该方法具有较高的参考价值^［7］。

枝下高指树冠第一轮活枝距离地面的高度^［8］，是用来计算冠长和描述树冠大小的重要指标^［9］。此外，枝下高还可以反映树木生长活力、生产力和树木在林分中的竞争水平，对于树木生长的量化以及推断树木长期的竞争状态有重要的意义^［10-11］。在林分郁闭度过高或林况复杂的实地测量中，枝下高的测量会消耗大量的人力、物力和财力^［12］。为解决这一问题，建立枝下高模型用来预测树木枝下高尤为重要。国内外学者对枝下高的研究始终保持着动态发展的态势，不断取得新进展。Logistic是模拟枝下高最常用的模型形式，相关学者^［13-15］使用了Logistic及其多种变形拟合了美国14个树种的枝下高、云冷杉针阔混交林枝下高和红松枝下高模型，并获得较好的拟合精度和预测效果。这些枝下高研究中，林木和林分因子是影响枝下高的主要因子，林木因子包括胸径、树高和冠幅^［14］等，林分因子^［15］包括优势木胸径、优势木树高等。同时，竞争因子与枝下高显著相关，常用林分断面积^{［10，13］}和大于对象木断面积之和^［16］等因子表示，并且能显著提高模型的拟合效果和预测精度^［16］。目前，对天然林枝下高的研究较少，Fu等^［16］将坡度和坡向作为立地因子加入到冠幅模型中，拟合效果得到了显著的提高。De等^［18］考虑种内和种间的断面积大于对象木之和、种内和种间的断面积代表竞争因子提高树高-胸径模型拟合精度并获得了较好的预测效果。曹晓梅等^［20］将代表物种多样性因子的Shannon指数加入到帽儿山阔叶混交林树高-胸径模型并获得了较好的拟合效果和预测精度。

简单的模型使用普通最小二乘法进行估计，该方法需要遵循随机和独立观测的假设，然而建立枝下高模型，所需要的数据通常是从若干个样地中选择单个树木收集的，这些样地也嵌套在同一块林分内^［16］。针对具有嵌套结构且缺乏独立性的数据，众多学者已经提出了建立混合效应模型的方法，该方法不仅有效地解决了嵌套结构数据带来的挑战还提高了模型的预测精度^［17］。现在混合效应模型广泛应用于树高曲线^［20］、削度方程^［21］和林分生长模型^［21］等，Fu等^［16］等建立了区组和样地两水平的非线性混合效应模型与传统枝下高模型进行对比，并成功证明含随机效应的混合效应模型预测能力更加显著。易达等^［9］建立了长白落叶松枝下高非线性混合效应模型和分位数组合模型，并得出混合效应模型可以更加准确地预测枝下高。段光爽等^［10］将区组和嵌套在区组里的样地分别作为随机效应，构建两水平枝下高模型，并证明其预测精度相对于传统模型得到了明显提升。

帽儿山地区的林分类型主要是天然次生林，林分结构不稳定，树种组成复杂，在该地区建立生长模型有助于了解帽儿山地区林木生长规律与生长潜力，预测森林资源未来的发展趋势，因此在帽儿山地区建立枝下高生长模型具有重要意义。近年来，围绕帽儿山地区建立的生长模型大多为单木生长模型、断面积生长模型以及树高-胸径生长模型，但针对枝下高建模的文章却鲜有报道。因此，将样地中常见的18个树种按照树种生物学特性以及样地实测数据进行分组，再考虑林分、竞争以及物种多样性等因子对枝下高的影响，建立适用于帽儿山地区的枝下高混合效应模型，并选择最合适的二次抽样样本数，为帽儿山地区森林资源管理和可持续经营提供了一定的科学依据和技术支持。

1 研究数据与方法

1.1 研究区概况

本研究数据来自黑龙江省南部地区的尚志市帽儿山镇，地理坐标为127°18′0″—127°41′6″E，45°2′20″—45°18′16″N，该地区虽属温带季风气候，但大陆性明显，平均海拔为419 m，年平均气温2.4 ℃，年平均降水量为700 mm，土壤以暗棕壤为主，有较高的有机质含量和植物生长所需营养元素。森林资源丰富，植被种类繁多，植被属于长白山植物区系范畴，是由于人工干预形成的东北天然次生林区，该林区主要的林分类型有硬阔混交林、软阔混交林和白桦林等，林区内的主要树种有色木槭（Acer pictum）、水曲柳（Fraxinus mandshurica）、黄檗（Phellodendron amurense）、胡桃楸（Juglans mandshurica）和白桦（Betula platyphylla）等。

1.2 试验数据

本研究数据采用帽儿山实验林场不同林分设置的28块天然混交林样地复测数据，共9 796棵样木，样地面积在0.075~0.21 hm²。在2012、2017、2022年分别对样地内树木进行每木检尺，分别实测样木（胸径≥5 cm）的树种、状态、树高（H）、胸径（DBH）、枝下高（HCB）以及样木的冠幅（CW）。在林木结构复杂的天然林中，调查林分空间信息需要消耗巨大的人力、物力和财力，与距离有关的竞争指数和与距离无关的竞争指数相比，已被证明前者对林分生长的改善效果较小^［30］。因此，本研究采用的是与距离无关的竞争变量：大于对象木断面积之和（BAL）、林分每公顷断面积（BAS）。天然林中树种种内、种间的相互作用显著影响树木获得资源与能量的能力^［18］，计算了2个种内竞争变量（BAL_in为大于对象木的相同树种的每公顷断面积之和、BAS_in为与对象木为相同树种的每公顷断面积）和2个种间竞争变量（BAL_out为大于对象木的不同树种的每公顷断面积之和、BAS_out为与对象木不同树种的每公顷断面积），考虑到天然林中树种的多样性以及林分结构的复杂性，本研究计算了Shannon指数（Sn）、Simpson指数（Sp）、Pielou均匀度（Pi）指数用来描述物种多样性，用基尼（Gini）系数来描述结构多样性。此外，本实验数据共有18个树种，计算了所有各树种每公顷断面积占比（BAP）、各树种平均高径比（\overline{\mathrm{H}\mathrm{D}\mathrm{R}}）、胸径变异系数（CV）。计算单木、林分、竞争和样地多样性因子见表1，树种信息见表2。

1.3 研究方法

1.3.1 树种分组

由表2可知样地内树种信息，帽儿山阔叶混交林具有较高的树种多样性，在建立生长模型时要考虑树种多样性对枝下高的影响。因此，探索一种科学有效的树种分组方法，以精准预测枝下高，成为本研究的重要目标。参考Zhao等^［7］的分类变量依据数据实际情况选取树种胸径变异系数、树种平均胸径、树种平均高径比、树种占比、阔叶树种（软、硬）分类以及树种的耐荫性等指标，使用层次聚类分析中的凡尔德（Ward）算法对树种进行分类。Ward算法的定义是计算2个类之间所有变量的方差分析的平方和，其倾向于把有少量观测值的类聚合到一起，并且倾向于产生与观测值个数大致相等的类。

1.3.2 枝下高基础模型的选择和广义模型的建立

选择5个常见的枝下高基础模型作为备选，见表3，利用所有实测数据对5个模型进行拟合，通过比较拟合精度，选择一个最优的模型用于非线性混合效应枝下高模型的构建。

除了树高和胸径对枝下高的影响，还应考虑林分、竞争和物种多样性等因子对枝下高的影响建立广义枝下高（H_CB）模型（M6）。其模型形式为

H_{\mathrm{C}\mathrm{B}}=f(\mathrm{T}\mathrm{r}\mathrm{e}\mathrm{e},\mathrm{S}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{d},\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n},\mathrm{D}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{y})

式中：f\left(\bullet \right)代表最优模型形式；Tree表示单木因子如胸径（DBH）、树高（H）和高径比（HDR）等；Stand表示林分因子如林分年龄（Sg）、优势木平均高（HT）和优势木平均胸径（Ddom）等；Competition表示竞争因子如大于对象木断面积之和（BAL）、林分断面积（BAS）等；Diversity表示物种多样性因子如香农指数（Shannon index）、辛普森指数（Simpson index）和皮卢指数（Pielou index）等。

通过再参数化的方法对以上因子作为自变量进行筛选并引入到枝下高基础模型中，考虑模型变量之间的相关性以及各变量之间的共线性（使用多重共线性检验选择方差膨胀因子（VIF）小于5）。通过对模型的调整后决定系数（R²）、均方根误差（RMSE）和赤池信息准则（AIC）等拟合指标的比较，确定最优的变量组合，用于构建广义枝下高模型。

1.3.3 枝下高混合效应模型

确定好基础枝下高模型和广义枝下高模型之后，基于对树种的分组，考虑各变量在样地和树种组中的差异对模型的影响，最后在基础枝下高模型和广义枝下高模型参数上引入样地、树种组随机效应，对所有随机参数的组合进行拟合得到两水平非线性混合效应枝下高模型，模型的一般表达式为

\begin{array}{c}{H}_{\mathrm{C}\mathrm{B}ijk}=f({\mathit{\phi }}_{ijk},v_{ijk})+{\epsilon }_{ijk},i=\mathrm{1},\dots ,M_{i,}k=\mathrm{1},\dots ,n_{ij}\mathrm{。}\\ {\mathit{\phi }}_{ijk}={\mathit{X}}_{ijk}\mathit{\beta }+{\mathit{Z}}_{i,jk}{\mathit{b}}_i+{\mathit{Z}}_{ijk}{\mathit{b}}_{ij}\mathrm{。}\\ b_i~N\left(\mathrm{0},{\mathit{D}}_{\mathrm{1}}\right),b_{ij}~\left(\mathrm{0},{\mathit{D}}_{\mathrm{2}}\right),{\epsilon }_{ij}~N\left(\mathrm{0},R_{ij}\right)\mathrm{。}\end{array}（1）

式中：H_{\mathrm{C}\mathrm{B}ijk}表示地i块样地第j个树种组内的第k个枝下高观测值；f(·)为枝下高基础模型；v_{ijk}表示其他预测变量的线性组合；{\mathit{\phi }}_{ijk}为参数向量；{\epsilon }_{ijk}表示服从正态分布的组内误差；{\mathit{X}}_{ijk}、{\mathit{Z}}_{i,jk}\mathrm{和}{\mathit{Z}}_{ijk}分别是\mathit{\beta }\mathrm{、}{\mathit{b}}_i\mathrm{和}{\mathit{b}}_{ij}对应的设计矩阵；\mathit{\beta }为固定参数向量；{\mathit{b}}_i\mathrm{和}{\mathit{b}}_{ij}为样地和树种组水平的随机参数向量，假定服从期望为0的方差-协方差矩阵 D₁和 D₂的正态分布；M为样地的数量；M_i为第i个样地内的树种组数；n_ij 为第i个样地中第j个树种组中的枝下高观测值数量；\mathit{\epsilon }、{\mathit{b}}_i\mathrm{和}{\mathit{b}}_{ij}相互独立本研究使用方差-协方差矩阵{\mathit{R}}_{ij}，其结构为

{\mathit{R}}_{ij}={\sigma }^{\mathrm{2}}{\mathit{G}}_{ij}^{\mathrm{0.5}}{\mathit{\Gamma }}_{ij}{\mathit{G}}_{ij}^{\mathrm{0.5}}。（2）

式中：{\sigma }^{\mathrm{2}}代表误差扩散的比例因子（即模型的误差方差值）；{\mathit{G}}_{ij}为样地数据方差异质性的对角矩阵；{\mathit{\Gamma }}_{ij}为样地数据组内误差相关性结构，目前本研究的样地实测数据并没有发现自相关结构，因此{\mathit{\Gamma }}_{ij}假定为单位矩阵。

考虑样地内部数据的异方差性，因此本研究采用常用的幂函数来消除异方差。其公式为

\mathrm{v}\mathrm{a}\mathrm{r}\left({\epsilon }_{ijk}\right)={\delta }^{\mathrm{2}}{\widehat{H}}_{\mathrm{C}\mathrm{B}ijk}^{\mathrm{2}\delta }。（3）

式中：\delta为待估计的幂函数参数；{\widehat{H}}_{\mathrm{C}\mathrm{B}ijk}^{\mathrm{2}\delta }为固定效应参数的预测值。

1.3.4 模型拟合和检验评价

利用全部数据分别对基础枝下高模型和广义枝下高模型（采用R软件4.1.4中的nls包）、基础混合效应枝下高模型和广义混合效应枝下高模型（采用R软件4.1.4中的nlme）进行拟合^［9］，并使用调整好的决定系数（R_{\mathrm{a}}^{\mathrm{2}}）、均方根误差（RMSE，式中记为R_MSE）以及赤池信息准则（AIC，式中记为A_IC）等3个拟合统计指标来评价拟合效果，模型的检验采用“留一法”获得枝下高预测值^［20］，通过平均绝对误差百分比（MAPE，式中记为M_APE）、平均绝对误差（MAE，式中记为M_AE）等指标来评估构建的基础枝下高模型、广义枝下高模型以及非线性混合效应枝下高模型的预测性能。其公式为

R^{\mathrm{2}}=\mathrm{1}-\left(\frac{n-\mathrm{1}}{n-\lambda }\right)\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}(y_i-{\widehat{y}}_i{)}^{\mathrm{2}}}{{\displaystyle \sum _{i=\mathrm{1}}^n}(y_i{-\overline{y})}^{\mathrm{2}}}。（4）

R_{\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt[]{\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}(y_i-{\widehat{y}}_i{)}^{\mathrm{2}}}{n-\lambda }}。（5）

A_{\mathrm{I}\mathrm{C}}=-\mathrm{2}LL+\mathrm{2}\lambda。（6）

M_{\mathrm{A}\mathrm{P}\mathrm{E}}=\mathrm{100}\times {\displaystyle \sum _{i=\mathrm{1}}^n}\left|\frac{\left(y_i-{\widehat{y}}_{i,-i}/\overline{y}\right)}{n}\right|。（7）

M_{\mathrm{A}\mathrm{E}}={\displaystyle \sum _{i=\mathrm{1}}^n}\left|\frac{y_i-{\widehat{y}}_{i,-i}}{n}\right|。（8）

式中；n为总样本数量；\lambda为模型参数个数；y_i为第i个枝下高的观测值；{\widehat{y}}_i是全部样地数据拟合回归方程计算的预测值；\overline{y}是观测值的平均值；{\widehat{y}}_{i,-i}为剔除第i样本观测值个后拟合模型得到的y_i的预测值。

1.3.5 随机效应模型的抽样设计

非线性混合效应模型随机效应参数的检验需要计算随机参数值，通常利用预测样地i中至少存在的一个已知观测数据来计算随机效应参数。其公式为

{\widehat{\mathit{u}}}_i=\widehat{\mathit{D}}{\widehat{\mathit{Z}}}_i^{\mathrm{T}}{\left({\widehat{\mathit{Z}}}_i\widehat{\mathit{D}}{\widehat{\mathit{Z}}}_i^{\mathrm{T}}+{\widehat{\mathit{R}}}_i\right)}^{-\mathrm{1}}{\widehat{e}}_i。（9）

式中：{\widehat{\mathit{u}}}_i=\left({\mathit{u}}_i^{\mathrm{T}},{\mathit{u}}_{i\mathrm{1}}^{\mathrm{T}},{\mathit{u}}_{i\mathrm{2}}^{\mathrm{T}},\dots ,{\mathit{u}}_{i{M}_i}^{\mathrm{T}}\right)为\left(q_{\mathrm{1}}+M_i{q}_{\mathrm{2}}\right)维随机效应参数向量；M_i为第i个样地中的树种组数；q_{\mathrm{1}}为样地水平随机效应的个数；q_{\mathrm{2}}为树种组随机效应的个数；\widehat{\mathit{D}}为\left(q_{\mathrm{1}}+M_i{q}_{\mathrm{2}}\right)\times \left(q_{\mathrm{1}}+M_i{q}_{\mathrm{2}}\right)维方差协方差矩阵，{\widehat{\mathit{R}}}_i为n_i×n_i维的组内方差协方差矩阵；{\widehat{e}}_i为实际观测值与固定效应预测值的差值；{\widehat{\mathit{Z}}}_i为n_i\times \left(q_{\mathrm{1}}+M_i{q}_{\mathrm{2}}\right)维随机效应设计矩阵^［21］。

对混合效应枝下高模型使用随机抽样产生校正参数值，基于经验最优线性无偏预测（empirical best linear unbiased prediction，EBLUP）^［20］的方法，利用部分样地数据来校正模型。使用常见的随机抽样方案对非线性混合效应枝下高模型进行校正并选择最优结果，抽样方案为随机抽取1~8棵树（该方案重复500次，由R软件4.1.4中的Matrix包完成）

2 结果与分析

2.1 树种分类

树种分组的具体结果以及各树种组的密度分布曲线已分别展示在图1和图2中。其中，A组被归类为喜光树种组，并严格遵循王鹤智^［22］对东北地区阔叶树种的分类标准，即软阔树种与硬阔树种的划分，A组内的所有树种均属于硬阔树种。该组的树木主要由胸径在20~35 cm范围内的大树构成，且胸径分布在大约28 cm处达到峰值。B组和C组则以喜光树种为主，且两者均由软阔树种构成。具体而言，B组树木的胸径分布相对均匀，大多数树木的胸径值接近该组的平均胸径。相比之下，D组则包含了较多的耐荫树种，且林分主要由软阔树种组成。值得注意的是，C组和D组的树木主要由胸径<10 cm、年龄较小或刚完成更新的幼树构成。为了验证这4个树种组之间是否存在显著差异，本研究采用了R软件4.1.4版本中的multcomp包进行方差分析。结果显示，各分组之间存在极显著差异（P<0.001），这进一步证实了树种分组方法的合理性和有效性。

2.2 枝下高基础和广义模型

用所有数据对表3中5个常用的枝下高基础模型进行拟合，结果见表4。5个模型的决定系数均较小。最终选择R^{\mathrm{2}}最大，RMSE和AIC最小的模型M4为本研究的枝下高基础模型，并且M4计算的枝下高永小于树高。

由于枝下高不仅受胸径和树高的影响，本研究考虑单木、林分、竞争和物种多样性等4个方面对模型拟合性能的影响，通过不断引入其他变量进行筛选，并且考虑模型的实际意义，最终选取优势木平均高（HT，式中记为H_T）、种内竞争每公顷大于对象木断面积之和（BAL_in，式中记为B_ALin）以及Shannon指数（式中记为S_hannon）等3个变量来构建广义枝下高模型M6，由图3可以看出，所选的3个变量与模型各变量之间的相关性，并且选入的变量均通过共线性检验，VIF均小于（D_BH=2.83、H=2.95、H_T=1.09、B_{\mathrm{A}\mathrm{L}\mathrm{i}\mathrm{n}}=1.02、S_hannon=1.05）。广义枝下高模型（M6）形式为

\begin{array}{c}{H}_{\mathrm{C}\mathrm{B}}=H(\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}(a_{\mathrm{0}}+a_{\mathrm{1}}{D}_{\mathrm{B}\mathrm{H}}+a_{\mathrm{2}}{H}_{\mathrm{T}}+\\ a_{\mathrm{3}}{B}_{\mathrm{A}\mathrm{L}\mathrm{i}\mathrm{n}}+a_{\mathrm{4}}{S}_{\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{o}\mathrm{n}}))\end{array}。（10）

式中：a_{\mathrm{0}}、a_{\mathrm{1}}、a_{\mathrm{2}}、a_{\mathrm{3}}和a_{\mathrm{4}}为模型参数；H_{\mathrm{T}}为优势木平均高；B_ALin为种内竞争大于对象木断面积之和；S_{\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{o}\mathrm{n}}为Shannon指数；D_BH为胸径。

广义模型M6的拟合结果见表5，对比枝下高基础模型，可以看出加入了协变量后，模型拟合效果得到提升，R²提高了3.6%，RMSE下降了1.8%。所选变量中胸径（DBH））和树高（H）与枝下高呈负相关，优势木平均高（HT）、种内竞争大于对象木断面积之和（BAL_in）以及Shannon指数与枝下高呈正相关。

2.3 样地、树种组两水平枝下高混合效应模型

在基础枝下高模型M4的不同参数位置上添加树种组随机效应和样地随机效应，构建样地、树种组两水平基础枝下高混合模型M7。通过对模型拟合指标的对比，最终发现样地随机效应作用在a_{\mathrm{0}}、a_{\mathrm{1}}，树种组随机效应作用在a_{\mathrm{1}}上效果最好。基础枝下高混合模型M7模型形式为

H_{\mathrm{C}\mathrm{B}}=H(\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}(a_{\mathrm{0}}+a_{\mathrm{0}i}+(a_{\mathrm{1}}+a_{\mathrm{1}i}+a_{\mathrm{1}ij})D_{\mathrm{B}\mathrm{H}})（11）

式中：a_{\mathrm{0}i}、a_{\mathrm{1}i}为样地随机效应参数；a_{\mathrm{1}ij}为树种组随机效应参数。

其模型的参数估计值和拟合优度统计见表5，模型的参数估计值均显著，相比与基础枝下高模型，混合模型的R²为0.627（提升了21.2%），RMSE和AIC分别为1.909和41 614.23（降低了16.7%和4.5%）。结果表明基础枝下高混合模型M7有更小的预测误差。

在广义枝下高模型M6的基础上，本研究采用相同的方法构建两水平广义枝下高混合模型M8。通过对不同的参数加入随机效应后模型的拟合进行对比，最终发现样地随机效应作用于a_{\mathrm{1}}、a_{\mathrm{2}}、a_{\mathrm{3}}，树种组随机效应作用于a_{\mathrm{2}}、a_{\mathrm{3}}效果最好。广义枝下高混合模型M8形式为

\begin{array}{l}{H}_{\mathrm{C}\mathrm{B}}=H(\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}(a_{\mathrm{0}}+a_{\mathrm{1}i}{D}_{\mathrm{B}\mathrm{H}}+(a_{\mathrm{2}i}+a_{\mathrm{2}ij})H_{\mathrm{T}}+\\ (a_{\mathrm{3}}+a_{\mathrm{3}i}+a_{\mathrm{3}ij})B_{\mathrm{A}\mathrm{L}\mathrm{i}\mathrm{n}}+a_{\mathrm{4}}{S}_{\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{o}\mathrm{n}}))\end{array}。（12）

式中：a_{\mathrm{1}i}、a_{\mathrm{2}i}和a_{\mathrm{3}i}为样地随机效应参数；a_{\mathrm{2}ij}、a_{\mathrm{3}ij}为树种组随机效应参数

表5为该模型的参数估计值及拟合优度统计，模型的参数估计值均显著。相比于枝下高广义模型M6，广义枝下高混合效应模型M8的拟合指标均得到了提升（R²提升了19.7%，RMSE和AIC均下降，分别是16.1%和3.8%）。模型的MAE和MAPE分别为1.410 m和27.3%，与广义枝下高模型M6相比分别降低了14.8%和18.3%。相比与基础枝下高模型M4，广义枝下高混合模型M8各拟合精度略高（R²提升了0.3%，RMSE和AIC均下降，分别是0.2%和2.2%，MAE和MAPE分别下降了0.1%和0.03%）。结果表明，引入样地随机效应和树种组随机效应建立两水平广义枝下高混合M8模型相比于其他模型得到了更小的预测误差。

通过利用广义枝下高混合模型M8的固定效应确定各协变量的参数取值，通过对胸径（DBH）、优势木平均高（HT）、种内竞争大于对象木断面积之和（BAL_in）以及Shannon指数的梯度取值，分析所选各因子在相同树高下对枝下高的影响，具体做法是对目标变量取等差梯度值，其他非目标协变量均取均值。图4表明，相同树高下，胸径越大的树木有更低的枝下高，而枝下高随着优势木平均高（HT）、种内竞争大于对象木断面积之和（BAL_in）以及Shannon指数的梯度增加而增加。

2.4 枝下高混合效应模型抽样校正

图5为对基础混合效应模型M7和广义混合效应模型M8进行随机抽样校正的结果。随机抽样的所有样本大小拟合结果相比于固定效应模型都有更高的预测精度，随机抽样样本大小由1到8，由于部分样地中某种树种组的树木不足8棵，对这种情况只能继续抽树种组中包含的树木。图5和表6可以看出，随着随机抽样样本数量增加，混合效应模型的预测精度更高，当每个树种组选择1棵树时MAE、MAPE的下降速率分别为-7.51%和-13.51%。当随机抽样样本数量大于1时模型拟合精度提升较小，因此认为随机抽取每个样地内每个树种组中的1棵树校正模型较为合理。

3 讨论

3.1 基于层次聚类分析的树种分组

在帽儿山的天然阔叶混交林样地中，树木种类繁多，形成了一个具有较高生物多样性的生态系统。在构建混交林枝下高模型的初期阶段，首要的步骤是对该样地中的树种进行科学且合理的分组^［8］。本研究运用了Ward算法进行层次聚类分析，成功将帽儿山常见的18个树种划分为4个树种组。选取了与曹晓梅等^［20］相同的树种分类指标包括CV、BAP、\overline{\mathrm{D}\mathrm{B}\mathrm{H}}和树种耐荫性等，还考虑了\overline{\mathrm{H}\mathrm{D}\mathrm{R}}、软阔叶和硬阔叶树种的分类。胸径变异系数作为衡量林分内胸径分布离散度与复杂性的有效指标^［23］，其变动直观反映了林分结构的多样性；而各树种占比的差异化分布^［20］，则直接决定了林分整体结构的繁复程度。林分平均胸径^［24］与高径比^［25］被选定为描述林分生长状态的两大要素。鉴于树种种类繁多，其胸径与高径比呈现出明显的多样性特征，尤其是那些生长速度较快的树种，通常具有较大的高径比。这一现象为深入理解林分生长动态提供了有价值的参考信息。耐荫性^［2］作为树种的关键生态特征之一，在混交林中扮演着分配与利用光资源的关键角色，深刻影响着林分的生态位配置。通过将帽儿山阔叶混交林的阔叶树种划分为硬阔与软阔两大类^［22］，进一步细化了树种的功能属性，为深入研究提供了理论基础。最终，本研究利用R软件4.1.4中的multcomp包对4个树种组进行了方差分析，结果显示各分组间存在显著差异，进一步验证了本研究分类方法的科学性与合理性。

3.2 枝下高基础模型和广义模型构建

本研究的最优基础模型为维科夫（Wykoff）^［14］模型，通过再参数化的方法将优势木平均高（HT）、种内竞争大于对象木断面积之和（BAL_in）以及Shannon指数引入模型，且对比基础模型获得了更好的预测精度。枝下高（HCB）和胸径（DBH）之间是负相关，树木胸径越大的树木有着更强的竞争能力，有更低的枝下高^［16］。优势木平均高（HT）已被证明与树木生长和林分发育之间有改善作用，因此本研究将其纳入模型。图4结果显示，优势木平均高（HT）的梯度增加表示更好的立地质量促进树木生长，使树冠向上移动且树冠下部的较大枝条会发生自然整枝，枝下高整体上移^［9，28］。树木拥挤造成的竞争对枝下高有较强的影响^［16-17］。本研究选择了种内竞争大于对象木断面积之和（BAL_in）来描述种内树木之间的竞争^［18］，结果表明枝下高（HCB）随着种内竞争大于对象木断面积之和（BAL_in）增加而增加，这与大多数文章相同^{［9，16-17］}。与同种树木相比，种内竞争大于对象木断面积之和（BAL_in）越大，树木直径越小，该树木竞争能力较弱，树木的竞争压力变大，枝下高增加^{［18，27］}。与其他天然混交林单木模型相同^［19-20］，本研究纳入了物种多样性指数Shannon，并且发现Shannon指数的梯度增加与枝下高呈正相关，树种的增加会调整天然阔叶混交林的整体结构，促进树木对资源的高效利用并提高树种的竞争与适应能力，从而更利于树木生长^［19］。本研究所选的变量均易获取并考虑了经济成本，帽儿山天然阔叶混交林以其植被组成丰富、林分结构鲜明和生态功能强大的特点，为本研究提供了建立枝下高模型的理想场所。

3.3 枝下高混合效应模型构建与抽样校正分析

本研究参考曹晓梅等^［20］的样地随机与树种组两水平随机效应树高-胸径模型，在基础枝下高模型M4^［15］和广义枝下高模型M6中引入样地随机效应和树种组随机效应，充分考虑了样地嵌套结构和不同树种特性对枝下高的影响。结果显示，广义枝下高混合模型M8拟合结果均优于其他3种模型。对混合效应模型进行随机抽样校正的文章较多^{［16-17，31］}，研究表明预测精度与样本量大小呈正相关，但是影响幅度逐渐平缓。本研究对2种混合效应模型的抽样校正作对比^［15］结果显示，两模型的随机抽样结果相比于固定效应都有更精准的预测效果^［20］，本研究的抽样方法与Fu等^［16］一致，且Crecente等^［30］指出，在含有样地和树种组随机效应的模型，可以通过抽取样地内每一树种组的1棵树对模型进行校正。表6显示，当随机抽取样地内每个树种组的1棵树木校正基础枝下高混合模型M7和广义枝下高混合模型M8时，两模型的MAE和MAPE增长率分别为-4.59%、-7.51%和-9.47%、-13.51%，相比于其他样本数量的增长率明显较大。因此，选择随机抽取样地内树种组的1棵树校正模型较为合理，且广义枝下高混合模型M8预测精度表现较好。

4 结论

本研究结合样地数据选择了树种分类标准，并使用Ward算法将18个树种分为4个树种组，作为随机效应加入到基础枝下高混合效应模型M7和广义枝下高混合效应模型M8中，不仅解决了某种树木数量较少无法建模的问题，还考虑到树种多样性对枝下高的影响。进一步引入了优势木平均高（HT）、种内竞争大于对象木断面积之和（BAL_in）以及Shannon指数等协变量，旨在提升广义枝下高模型M6的预测精度。通过将帽儿山阔叶混交林全部样地数据带入模型拟合，广义枝下高混合模型M8比基础枝下高混合模型M7表现出更好的拟合和预测能力。因此，本研究从每个样地的每个树种组中随机抽取1棵树对广义枝下高混合模型M8进行校正，并对枝下高进行预测。

参考文献

[1] 乔丹.我国天然林保育类型划分及保育策略研究［D］.北京：北京林业大学，2018.

[2] QIAO D.Study on classification and conservation strategies of natural forest conservation in China［D］.Beijing：Beijing Forestry University，2018.

[3] DEL RÍO M， PRETZSCH H， ALBERDI I，et al.Characterization of the structure，dynamics，and productivity of mixed-species stands：review and perspectives［J］.European Journal of Forest Research，2016，135（1）：23-49.

[4] 马志海.东北林区天然混交林单木生长模型的研究［D］.哈尔滨：东北林业大学，2006.

[5] MA Z H.Individual tree growth models of natural mixed forest in Northern China［D］.Harbin：Northeast Forestry University，2006.

[6] TEMESGEN H， ZHANG C H， ZHAO X H.Modelling tree height-diameter relationships in multi-species and multi-layered forests：A large observational study from Northeast China［J］.Forest Ecology and Management，2014，316：78-89.

[7] 刘耀凤，董利虎，郝元朔，等.竞争对天然针阔混交林树木直径生长的影响［J］.东北林业大学学报，2023，51（9）：1-7.

[8] LIU Y F， DONG L H， HAO Y S，et al.Effect of competition on tree growth in natural coniferous and broad-leaved mixed forest［J］.Journal of Northeast Forestry University，2023，51（9）：1-7.

[9] ZHAO D H， BORDERS B， WILSON M.Individual-tree diameter growth and mortality models for bottomland mixed-species hardwood stands in the lower Mississippi alluvial valley［J］.Forest Ecology and Management，2004，199（2/3）：307-322.

[10] ZHAO D H， BORDERS B， WILSON M，et al.Modeling neighborhood effects on the growth and survival of individual trees in a natural temperate species-rich forest［J］.Ecological Modelling，2006，196（1/2）：90-102.

[11] KUPREVICIUS A， AUTY D， ACHIM A，et al.Quantifying the influence of live crown ratio on the mechanical properties of clear wood［J］.Forestry，2013，86（3）：361-369.

[12] 易达，李凤日，马爱云，等.基于混合效应模型和分位数回归的长白落叶松枝下高模型构建［J］.应用生态学报，2023，34（4）：1035-1042.

[13] YI D， LI F R， MA A Y，et al.Model construction for height to crown base of Larix olgensis based on mixed-effects model and quantile regression［J］.Chinese Journal of Applied Ecology，2023，34（4）：1035-1042.

[14] 段光爽，李学东，冯岩，等.基于广义非线性混合效应的华北落叶松天然次生林枝下高模型［J］.南京林业大学学报（自然科学版），2018，42（2）：170-176.

[15] DUAN G S， LI X D， FENG Y，et al.Generalized nonlinear mixed-effects crown base height model of Larix principis-rupprechtii natural secondary forests［J］.Journal of Nanjing Forestry University （Natural Sciences Edition），2018，42（2）：170-176.

[16] RUSSELL M B， WEISKITTEL A R， KERSHAW J A.Comparing strategies for modeling individual-tree height and height-to-crown base increment in mixed-species Acadian forests of northeastern North America［J］.European Journal of Forest Research，2014，133（6）：1121-1135.

[17] TEMESGEN H， LEMAY V， MITCHELL S J.Tree crown ratio models for multi-species and multi-layered stands of southeastern British Columbia［J］.Forestry Chronicle，2005，81（1）：133-141.

[18] RITCHIE M W， HANN D W.Equations for predicting height to crown base for fourteen tree species in southwest Oregon［R］.Research paper-Forest Research Laboratory，Oregon State University（USA），1987.

[19] 李应涛，刘时良，孙海龙，等.云冷杉针阔混交林单木枝下高和冠幅模型构建［J］.森林与环境学报，2022，42（3）：289-296.

[20] LI Y T， LIU S L， SUN H L，et al.Modeling of height to crown base and crown width in spruce-fir conifer-broadleaf mixed forest［J］.Journal of Forest and Environment，2022，42（3）：289-296.

[21] 燕云飞.基于混合效应模型的人工红松枝下高模型研建［D］.哈尔滨：东北林业大学，2021.

[22] YAN Y F.Construction of the height to crown base mixed model for Korean pine［D］.Harbin：Northeast Forestry University，2021.

[23] FU L Y， ZHANG H R， SHARMA R P，et al.A generalized nonlinear mixed-effects height to crown base model for Mongolian oak in Northeast China［J］.Forest Ecology and Management，2017，384：34-43.

[24] SHARMA R P， VACEK Z， VACEK S.Individual tree crown width models for Norway spruce and European beech in Czech Republic［J］.Forest Ecology and Management，2016，366：208-220.

[25] DE PRADO D R， RIOFRÍIO J， ALDEA J，et al.Competition and climate influence in the basal area increment models for Mediterranean mixed forests［J］.Forest Ecology and Management，2022，506：119955.

[26] TIAN D Y， JIANG L C， SHAHZAD M K，et al.Climate-sensitive tree height-diameter models for mixed forests in Northeastern China［J］.Agricultural and Forest Meteorology，2022，326：109182.

[27] 曹晓梅，苗铮，郝元朔，等.基于树种分类的帽儿山阔叶混交林树高-胸径模型［J］.应用生态学报，2024，35（2）：307-320.

[28] CAO X M， MIAO Z， HAO Y S，et al.Height-diameter model of broad-leaved mixed forest based on species classification in Maoershan，North-east China［J］.Chinese Journal of Applied Ecology，2024，35（2）：307-320.

[29] 聂璐毅，董利虎，李凤日，等.基于两水平非线性混合效应模型的长白落叶松削度方程构建［J］.南京林业大学学报（自然科学版），2022，46（3）：194-202.

[30] NIE L Y， DONG L H， LI F R，et al.Construction of taper equation for Larix olgensis based on two-level nonlinear mixed effects model［J］.Journal of Nanjing Forestry University（Natural Sciences Edition），2022，46（3）：194-202.

[31] 王鹤智.东北林区林分生长动态模拟系统的研究［D］.哈尔滨：东北林业大学，2012.

[32] WANG H Z.Dynamic simulating system for stand growth of forests in Northeast China［D］.Harbin：Northeast Forestry University，2012.

[33] 李建，李晓宇，曹静，等.长白山次生针阔混交林群落结构特征及群落动态［J］.生态学报，2020，40（4）：1195-1206.

[34] LI J， LI X Y， CAO J，et al.Community structure and dynamics of secondary coniferous and broad-leaved mixed forest in Changbai Mountains［J］.Acta Ecologica Sinica，2020，40（4）：1195-1206.

[35] 娄明华，张会儒，雷相东，等.天然栎类阔叶混交林林分平均高与平均胸径关系模型［J］.北京林业大学学报，2020，42（9）：37-50.

[36] LOU M H， ZHANG H R， LEI X D，et al.Relationship model between stand mean height and mean DBH for natural Quercus spp.broadleaved mixed stands［J］.Journal of Beijing Forestry University，2020，42（9）：37-50.

[37] 严铭海，黄清麟，王金池，等.林木高径比研究综述［J］.世界林业研究，2023，36（1）：59-65.

[38] YAN M H， HUANG Q L， WANG J C，et al.Research on height to diameter ratio of forest trees：A review［J］.World Forestry Research，2023，36（1）：59-65.

[39] FU L Y， SUN H， SHARMA R P，et al.Nonlinear mixed-effects crown width models for individual trees of Chinese fir（Cunninghamia lanceolata） in south-central China［J］.Forest Ecology and Management，2013，302：210-220.

[40] SHARMA P， VACEK Z， VACEK S，et al.Modelling individual tree height to crown base of Norway spruce（Picea abies （L.） Karst.） and European beech（Fagus sylvatica L.）［J］.Plos One，2017，12（10）：e0186394.

[41] KUEHNE C， WEISKITTEL A R， WASKIEWICZ J.Comparing performance of contrasting distance-independent and distance-dependent competition metrics in predicting individual tree diameter increment and survival within structurally-heterogeneous，mixed-species forests of Northeastern United States［J］.Forest Ecology and Management，2019，433：205-216.

[42] XIE L F， WIDAGDO F R A， DONG L H，et al.Modeling height-diameter relationships for mixed-species plantations of Fraxinus mandshurica Rupr.and Larix olgensis Henry in Northeastern China［J］.Forests，2020，11（6）：610.

[43] CRECENTE-CAMPO F， CORRAL-RIVAS J J， VARGAS-LARRETA B， et al. Can random components explain differences in the height-diameter relationship in mixed uneven-aged stands？［J］. Annals of Forest Science， 2014， 71（1）： 51-70.

基金资助

“十四五”国家重点研发计划项目(2022YFD2201000)

国家自然科学基金面上项目(31971649)

黑龙江省博士后科学基金项目(LBH-Z23050)