帽儿山天然林进界模型及影响因素分析

Analysis of the Ingrowth Model and Influencing Factors for the Natural Forests of Maoer Mountain

摘要

林分进界是林分生长动态变化过程中的一个重要环节，对于维持森林资源中的生物多样性和群落结构稳定至关重要。为此，在帽儿山实验林场设立的61块样地数据，从林分因子及生物多样性因子等方面进行考虑，通过肯德尔tau-b（Kendall-Tau-b）相关系数分析及考虑变量多重共线性选择最适变量，采用泊松（Poisson）模型、负二项模型（negative binomial model，NB）、零膨胀模型及障碍（Hurdle）模型构建进界模型。采用层次分割法分析变量贡献率，以找出影响进界模型的关键因素。研究结果表明，林分拥挤度（K）、林分算数平均胸径（d）、辛普森（Simpson）指数和林分平均高（mean height，MH）是影响每公顷进界株数（number of advance regeneration per hectare，Nn）的重要因子。通过赤池信息准则（akaike information criterion，AIC）、贝叶斯信息准则（bayesian information criterion，BIC）及对数似然值（logarithm likelihood，Loglike）指标对比，发现零膨胀负二项式（zero-inflated negative binomial，ZINB）模型与障碍负二项（hurdle negative binomial，HNB）模型显著优于其他模型。通过沃恩（Vuong）检验，发现负二项模型及其复合模型（ZINB、HNB）在拟合帽儿山天然林进界数量方面优于Poisson模型及其复合模型（零膨胀泊松模型（zero-inflated poisson model，ZIP）、障碍泊松模型（hurdle poisson model，HP）），且ZINB模型略优于HNB模型，故ZINB模型为拟合帽儿山天然林林分进界数量的最优模型，十折交叉检验也验证了这一结论。同时，通过层次分割法分析发现Simpson指数和林分平均高（MH）分别对最优进界模型（ZINB）的计数部分和零部分贡献度最高。所构建的帽儿山天然林进界模型具有一定的统计可靠性，可用于该地区进界生长预测，为当地天然林更新管理提供科学依据。

Abstract

Forest stand ingrowth is a critical component of the dynamic growth process of forest stands， essential for maintaining biodiversity and community structure stability in forest resources. Based on data from 61 plots established at the Maoer Mountain Experimental Forest Farm， this study considered factors such as stand characteristics and biodiversity. Through Kendall-Tau-b correlation coefficient analysis and selection of the most suitable variables considering multicollinearity among variables， models for ingrowth were constructed using Poisson， negative binomial （NB）， zero-inflated， and Hurdle models. The contribution rate of variables was analyzed using hierarchical partitioning to identify key factors influencing the ingrowth model. The results showed that stand density （K）， arithmetic mean diameter at breast height （d）， Simpson's index， and mean stand height （MH） were significant factors affecting the number of ingrowth trees per hectare （Nn）. Comparing models using AIC， BIC， and Loglike criteria， it was found that ZINB and HNB significantly outperformed other models. The Vuong test further revealed that the negative binomial models and their composite models （ZINB， HNB） performed better than Poisson models and their composites （ZIP， HP） in fitting the ingrowth quantity of natural forests in Maoer Mountain， with the ZINB model slightly outperforming the HNB model. Therefore， the ZINB model was the optimal model for fitting the ingrowth quantity of natural forest stands in Maoer Mountain， a conclusion also supported by ten-fold cross-validation. Additionally， hierarchical partitioning analysis indicated that the Simpson's index and mean stand height （MH） contributed most to the count and zero parts， respectively， of the optimal ingrowth model （ZINB）. The natural forest ingrowth model constructed by this research has a certain statistical reliability and can be used for ingrowth prediction in the Maoer Mountain area， providing a scientific basis for local natural forest regeneration management.

Graphical abstract

关键词

进界 / Kendall-Tau-b / 零膨胀模型 / Hurdle模型 / 层次分割

Key words

Ingrowth / Kendall-Tau-b / zero-inflated model / Hurdle model / hierarchical partitioning

0 引言

天然林在陆地碳汇平衡、减缓气候变化等方面具有十分重要的作用^［1］。而森林的自然更新和进界过程是维持天然林资源长期稳定的关键生态机制，对于保护生物多样性及保持群落结构的稳定性至关重要^［2-4］。进界指的是在特定时间间隔内树木达到一定生长标准（例如胸径达到5 cm）的现象。进界模型、生长模型和枯损模型构成了模拟森林群落动态变化的核心框架^［5］。由于数据较难获取及模型建立困难等，在研究林分生长模型时，对进界模型研究较少。实际上，考虑进界对于准确预测未来木材蓄积量非常重要；如果忽视这一环节，则可能导致对森林未来生长趋势及其潜在收获量的估计出现偏差^［6］。

进界作为林分动态变化过程中的重要组成部分，一直备受国内外学者的高度关注。目前，构建进界模型主要采用静态模型和动态模型2种方法。静态模型在矩阵模型中应用较多，例如，Buongiorno等^［7］和Grimes等^［8］通过建立矩阵模型对林分进界进行预测。然而，随着研究的深入，相关学者发现林分进界数据存在大量的零数据，传统的线性或非线性模型难以准确预测。如果仅考虑发生进界的林分，可能会导致结果偏高。相比之下，动态模型考虑了林分的动态变化，可以较好地解决上述问题。常见用于拟合动态进界模型的方法为两阶段法^［9］。例如马武等^［10］以蒙古栎天然林为研究对象，通过两阶段法建立了蒙古栎天然林的进界生长模型，结果表明该模型具有较良好的预测精度。但是，这些方法都是基于假设进界数量服从正态分布的前提构建的，而实际上林分进界数量是典型的离散随机数据，不服从正态分布。因此，目前国内外学者主要通过计数模型来更好地模拟林分进界过程。如Zell等^［11］、何江等^［12］采用常用的计数模型如泊松（Poisson）模型、负二项模型（negative binomial model，NB）、零膨胀负二项模型（zero-inflated negative binomial，ZINB）、零膨胀泊松模型（zero-inflated poisson model，ZIP）、障碍（Hurdle）-负二项模型（hurdle negative binomial，HNB）、Hurdle-泊松模型（hurdle poisson model，HP）等拟合林分进界过程，精准预测了林分进界数量。本研究也采用上述模型对帽儿山地区天然林进界数量进行预测。

前人对于林分天然更新的影响因素进行了深入研究，主要集中于林分因子、立地因子及气候因子等方面。如郑嵘等^［13］依据湿地松人工林数据，建立进界平均胸径模型，结果表明坡向、坡位、林龄、样地林木数量是重要的影响因子。陈科屹等^［14］的研究也证实了林分因子和立地因子对进界木的数量和分布情况产生影响。此外，Russell等^［15］的研究发现，食草动物对森林植被的啃食压力及林分断面积是影响美国北部森林天然更新的重要因素。不同地区的气候条件各异，这些气候因子对林木进界过程具有显著影响。如在高山森林生态系统林木进界会出现温度敏感性沿纬度梯度增加的现象^［16］。在温带地区，影响林木进界的主要驱动因素是林木的耐阴性和林分结构^［17］。李春明等^［18］和Zell等^［11］的研究也证明气候因子会对进界产生影响。

目前，国内外学者在构建进界模型时，很少采用系统的方法进行变量选择。这主要是因为进界数据通常呈现零膨胀特性，使用随机森林或皮尔逊相关性系数等方法选择变量可能导致模型拟合效果不佳或无法收敛。然而，肯德尔-tau（Kendall-Tau）相关系数在分析零膨胀数据的变量相关性方面表现出色，有助于选择出关键因子并构建模型。例如，Guo等^［19］在对高速公路分流区的交通事故率进行建模之前，就采用了Kendall-Tau-b相关系数来评估变量之间的相关性，以减少多重共线性对模型估计的影响。Quero等^［20］在对亚得里亚海中部粪便污染研究过程中，也通过Kendall-Tau相关系数评估了大肠杆菌浓度和肠球菌浓度的关系。故本研究采用Kendall-Tau相关系数筛选进界模型的变量，使模型构建更具合理性。

帽儿山天然林作为重要的生态屏障，研究其天然更新过程及分析影响因素对有效保护和合理利用该地区林分至关重要。因此，本研究以帽儿山实验林场设立的61块样地，共22 848株样木为研究对象，从林分因子及生物多样性因子等方面考虑，通过Kendall-Tau-b相关系数分析选择变量，采用Poisson模型、负二项模型（NB）、零膨胀模型及Hurdle模型构建进界模型，对比选择出预测精度最高的模型，并通过层次分割法分析各因子对模型的贡献度，找出影响进界数量的关键因素，旨在为当地天然林更新管理提供科学依据。

1 研究区域概况与数据来源

1.1 研究区域概况

研究区位于黑龙江省尚志市的帽儿山实验林场，地理坐标为127°29'—127°44'E，45°14'—45°29'N。该区域地势由南向北逐渐升高，平均海拔300 m，最高山峰为帽儿山，海拔805 m。属于温带湿润地区，属大陆性季风气候，年平均气温2.8 ℃，年降水量为723 mm左右，森林资源丰富，主要树种有椴树（Tilia tuan）、枫桦（Betula costata）、胡桃楸（Juglans mandshurica）、山杨（Populus davidiana）、水曲柳（Fraxinus mandschurica）、榆树（Ulmus pumila）和柞树（Quercus mongolica）等。

1.2 数据来源

本研究数据来源于帽儿山61块天然林复测数据，共有25种树种、22 848株样木，包括进界木495株。其中28块样地设立于2007年，之后每隔5 a进行复测，直至2022年。另外33块样地设立于2017年，于2022年进行复测。本研究所用样地面积为0.06~0.210 hm²，设立好样地后对样地内的每个样木进行每木检尺，记录胸径、树高、冠幅和枝下高等。然后对林分起源、林分年龄、郁闭度、海拔因子进行调查，并在5 a后对样地进行复测，记录样地的进界情况。整个调查区间内的每公顷进界株数（number of advance regeneration per hectare，Nn）的频数分布直方图如图1所示。

2 模型构建

2.1 变量的选择

通过查阅相关文献并结合影响林分进界的生物学过程，本研究选择的与进界相关的因子主要有以下两类见表1。

1）林分因子。林分内部的生长发育情况常常影响着进界，故在构建模型时考虑林分因子是必要的。因此，考虑的指标主要包括：每公顷进界株数（N_n）、优势高（D_H）、林分断面积（B_A）、林分拥挤度（K）^［21］、算数平均胸径（d）、郁闭度、林分平均胸径（M_D）、每公顷株树（N_HA）、优势木平均胸径（D_g）、林分年龄、林分平均高（M_H）、海拔（E_L）等。

2）生物多样性因子。天然林拥有丰富的植物种类和复杂的生态系统，生物多样性是对天然林建模研究中常用的变量，本研究中选择了一些常见生物多样性因子如下：Shannon-Weiner多样性指数、Simpson指数、Pielou均匀度指数等。

2.2 Kendall-Tau相关系数

Kendall-Tau^［22］常用来比较2个有序数列的相关性。主要有2种形式：Tau-a和Tau-b。这2种形式的主要区别在于其如何处理数据中的平局（即2个或更多观测值在某一变量上的值相同的情况）。本研究通过Kendall-tau-b相关系数分析选择合适变量构建进界模型。对于2个变量X和Y，Kendall-Tau的基本计算方法如下：如果对于任意2个观测值（X_i,Y_i）和（X_j,Y_j），当且或者X_i>X_j且Y_i>Y_j时，称这样的数对为一致对（Concordant）；反之，如果但Y_i>Y_j或者X_i>X_j且，则称这样的数对为不一致对（Discordant）。Kendall-tau相关系数主要通过一致对和不一致对的个数来计算。其公式为

{\tau }_{\mathrm{a}}=\frac{N_{\mathrm{c}}-N_{\mathrm{d}}}{\frac{\mathrm{1}}{\mathrm{2}}n\left(n-\mathrm{1}\right)}。（1）

{\tau }_{\mathrm{b}}=\frac{N_{\mathrm{c}}-N_{\mathrm{d}}}{\sqrt[]{(N_{\mathrm{0}}-N_{\mathrm{1}})(N_{\mathrm{0}}-N_{\mathrm{2}})}}。（2）

式中：N_{\mathrm{c}}表示一致对的数量；N_{\mathrm{d}}表示不一致对的数量；n是样本数量；N_{\mathrm{0}}=\frac{\mathrm{1}}{\mathrm{2}}n\left(n-\mathrm{1}\right)；N_{\mathrm{1}}和N_{\mathrm{2}}分别表示变量X和Y中的平局对数。

2.3 进界模型

2.3.1 Poisson模型

Poisson模型常用于分析计数数据，其概率质量函数（PMF，式中记为P）为^［23］

F\left(y_i\right)=P\left(Y_i=y_i\right)=\frac{e^{-{\lambda }_i}{\lambda }^{y_i}}{y_i!}。（3）

{\lambda }_i=\mathrm{E}\mathrm{x}\mathrm{p}\left(X_i\mathit{\beta }\right)。（4）

式中：F（y_i ）为y_i的函数；P（Y_i =y_i ）为Y_i =y_i 时的概率；y_i为随机变量；Exp（）是以自然对数为底的指数函数，Poisson分布的期望和方差相等，即E\left(Y_i\right)=\mathrm{V}\mathrm{a}\mathrm{r}\left(Y_i\right)={\lambda }_i；X_i为自变量（林分拥挤度（K）、林分平均高（MH）、林分算数平均胸径（d）、Simpson指数等）； β 为参数向量。

2.3.2 负二项模型（NB）

在Poisson分布中，离散参数θ是固定的；然而，在负二项分布^［24］中，离散参数θ是一个随机变量，遵循伽玛函数（Γ）分布。负二项模型通过引入一个离散参数θ来增强Poisson模型的能力，使其能更好地处理数据的异质性，因此在某些情况下比Poisson分布更具优势。与Poisson分布的方差（λ_i ）不同，负二项分布的方差为{\lambda }_i(\mathrm{1}+\theta {\lambda }_i)。当θ接近0时，负二项分布简化为Poisson分布，表明事件的发生具有较大的随机性；而当θ远离0时，则意味着事件之间存在聚集性和依赖性。负二项模型的概率质量函数（PMF）表达式为

\begin{array}{c}F(y_i)=P(Y=y_i)=\\ \frac{\mathrm{\Gamma }\left(y_i+{\theta }^{-\mathrm{1}}\right)}{\mathrm{\Gamma }\left(y_i+\mathrm{1}\right)\mathrm{\Gamma }\left({\theta }^{-\mathrm{1}}\right)}{(\frac{{\theta }^{-\mathrm{1}}}{{\theta }^{-\mathrm{1}}+{\lambda }_i})}^{{\theta }^{-\mathrm{1}}}{(\frac{{\lambda }_i}{{\theta }^{-\mathrm{1}}+{\lambda }_i})}^{y_i}\end{array}。（5）

{\lambda }_i=\mathrm{E}\mathrm{x}\mathrm{p}\left(X_i\beta \right)。（6）

式中：Γ为伽玛函数；θ为离散参数。

2.3.3 零膨胀（Zero-inflated）模型

零膨胀模型^［25］是一种专门用于处理含有大量零值的计数数据的统计模型。模型将事件发生数分为2种可能的情形：第1种情况是零事件的发生，假定服从伯努利分布；第2种情况是非零事件的发生，假定服从Poisson分布或负二项分布。设有一个服从零膨胀分布的离散随机变量y（每公顷进界株数），p_i为零部分的概率，其概率质量函数（PMF）表达式为^［26］

P\left(Y_i=y_i\right)=\left\{\begin{array}{c}{p}_i+\left(\mathrm{1}-p_i\right)f(\mathrm{0})\\ \left(\mathrm{1}-p_i\right)f\left(y_i\right)\end{array}\right.。（7）

式中：0＜p_i＜1，用于解释模型中的多零部分，零部分常用Logistic的Logit形式来拟合。其公式为

\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}\left(p_i\right)=\mathrm{L}\mathrm{o}\mathrm{g}\left(\frac{p_i}{\mathrm{1}-p_i}\right)=X_i\mathit{\delta }。（8）

式中：Log（）是以自然对数为底的对数函数； δ 为参数向量。

1）零膨胀Poisson模型（ZIP）

式（7）中若模型的计数部分采用Poisson分布拟合，就可以得到ZIP模型。ZIP模型的概率质量函数（PMF）公式为

P\left(Y_i=y_i\right)=\left\{\begin{array}{cc}{p}_i+\left(\mathrm{1}-p_i\right)\mathrm{E}\mathrm{x}\mathrm{p}\left(-{\lambda }_i\right),& y_i=\mathrm{0}\\ \frac{\left(\mathrm{1}-p_i\right)\mathrm{E}\mathrm{x}\mathrm{p}\left(-{\lambda }_i\right){\lambda }_i^{y_i}}{y_i!},& y_i>\mathrm{0}\end{array}\right.。（9）

对于ZIP分布，其期望\left(E\left(Y_i\right)\right)和方差\mathrm{V}\mathrm{a}\mathrm{r}\left(Y_i\right)分别为

E\left(Y_i\right)=\left(\mathrm{1}-p\right){\lambda }_i。（10）

\mathrm{V}\mathrm{a}\mathrm{r}\left(Y_i\right)=E\left(Y_i\right)\left(\mathrm{1}+{\lambda }_i-E\left(Y_i\right)\right)。（11）

2）零膨胀负二项模型（ZINB）

若模型的计数部分用负二项模型来拟合，得到ZINB模型。其概率质量函数（PMF）公式为

P\left(Y_i=y_i\right)=\left\{\begin{array}{cc}{p}_i+\left(\mathrm{1}-p_i\right){\left(\frac{{\theta }^{-\mathrm{1}}}{{\theta }^{-\mathrm{1}}+{\lambda }_i}\right)}^{\mathrm{1}/\theta },& y_i=\mathrm{0}\\ \left(\mathrm{1}-p_i\right)\frac{\mathrm{\Gamma }\left(y_i+{\theta }^{-\mathrm{1}}\right){{\lambda }_i}^{y_i}}{\mathrm{\Gamma }\left(y_i+\mathrm{1}\right)\mathrm{\Gamma }\left(\theta \right)}{\left(\frac{{\theta }^{-\mathrm{1}}}{{\theta }^{-\mathrm{1}}+{\lambda }_i}\right)}^{\mathrm{1}/\theta }\left(\frac{{\lambda }_i}{{\theta }^{-\mathrm{1}}+{\lambda }_i}\right)y_i,& y_i>\mathrm{0}\end{array}\right.。（12）

对于ZINB分布，其期望和方差分别为

E\left(Y_i\right)=\left(\mathrm{1}-p\right){\lambda }_i。（13）

\mathrm{V}\mathrm{a}\mathrm{r}\left(Y_i\right)=E\left(Y_i\right)\left(\mathrm{1}+{\lambda }_i\left(\mathrm{1}+\theta \right)-E\left(Y_i\right)\right)。（14）

2.3.4 障碍（Hurdle）模型

Hurdle模型由Mullahy^［27］提出，该模型由2部分组成：第1部分用于模拟零值的出现，通常采用logit模型；第2部分用于模拟大于零的计数值，通常采用泊松模型或负二项模型。Hurdle模型的概率质量函数（PMF）公式为^［5］

P\left(Y_i=y_i\right)=\left\{\begin{array}{ccc}p& ,& y_i=\mathrm{0}\\ (\mathrm{1}-p)\frac{f\left(y_i\right)}{\mathrm{1}-f(\mathrm{0})}& ,& y_i>\mathrm{0}\end{array}\right.。（15）

1）Hurdle-Poisson模型（HP）

在Hurdle模型中用Poisson模型拟合计数部分，可得到HP模型，其概率质量函数（PMF）公式为

P\left(Y_i=y_i\right)=\left\{\begin{array}{cc}p,& y_i=\mathrm{0}\\ \frac{\left(\mathrm{1}-p\right)\mathrm{E}\mathrm{x}\mathrm{p}\left(-{\lambda }_i\right){\lambda }_i^{y_i}}{\left(\mathrm{1}-\mathrm{E}\mathrm{x}\mathrm{p}\left(-{\lambda }_i\right)\right)y_i!},& y_i>\mathrm{0}\end{array}\right.。（16）

2）Hurdle-NB模型（HNB）

Hurdle模型的另一种形式是Hurdle-NB模型（HNB），其概率质量函数（PMF）公式为

P\left(Y_i=y_i\right)=\left\{\begin{array}{ccc}P& ,& y_i=\mathrm{0}\\ \left(\mathrm{1}-p\right)\frac{\mathrm{\Gamma }\left(y_i+{\theta }^{-\mathrm{1}}\right)}{\mathrm{\Gamma }\left({\theta }^{-\mathrm{1}}\right)\mathrm{\Gamma }\left(y_i+\mathrm{1}\right)}{\left[{\left(\mathrm{1}+\theta {\lambda }_i\right)}^{\theta -\mathrm{1}}-\mathrm{1}\right]}^{-\mathrm{1}}{\left(\frac{{\lambda }_i}{{\lambda }_i+{\theta }^{-\mathrm{1}}}\right)}^{y_i}& ,& y_i>\mathrm{0}\end{array}\right.。（17）

通过软件R 3.4.2.中pscl软件包采用极大似然估计法对上述6个模型的参数进行拟合。

2.4 模型评价与检验

对于Poisson模型、NB模型、ZINB模型、ZIP模型、HNB模型和HP模型的拟合情况采用赤池信息准则（akaike information criterion，AIC）、贝叶斯信息准则（bayesian information criterion，BIC）以及对数似然值（logarithm likelihood，Loglike）这3个指标来进行比较。AIC、BIC值越小，Loglike值越大代表拟合效果越好。沃恩（Vuong）检验在比较一般计数模型与复合模型的分析中表现出较高的效能^［28］。当AIC值相近时，采用Vuong检验来比较拟合模型。模型采用十字交叉检验来评估各模型的预测精度，使用均方根误差（RMSE）和平均绝对误差（MAE）2个指标对模型进行评价，RMSE和MAE越小，模型预测精度越高。其中，RMSE（式中记为R_MSE）和MAE（式中记为M_AE）的公式为

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

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

式中：y_i为实测值；{\widehat{y}}_i为模型预测值；n为样本数；p为模型参数个数。

2.5 基于层次分割法分析影响因素

通过各指标对模型评价与检验后，选择出最优模型。因不同变量对帽儿山天然林进界数量的相对贡献问题较为复杂，故采用层次分割法的分析获得每个自变量对每公顷进界株数（N_n）的贡献率，以量化不同 变量对进界数量的影响，更深入理解各变量间的相 互作用和相对重要性，以上分析通过使用R包中的glmm.hp^［29-30］计算得到。

3 结果

3.1 模型变量选择

通过图2可以看出林分因子中算数平均胸径（d）、林分拥挤度（K）、优势高（D_H）、林分平均胸径（M_D）、优势木平均胸径（D_g）和林分平均高（M_H）与每公顷进界株数（N_n）的相关性较高，生物多样性因子中则是Simpson指数与每公顷进界株数（N_n）的相关性最高。但由于林分算数平均胸径（d）与林分平均胸径（M_D）和优势木平均胸径（D_g）、林分平均高（M_H）与优势高（D_H）之间的相关性较高，通过Kendall-Tau-b相关系数分析及考虑变量之间的多重共线性，最终选择林分拥挤度（K）、林分算数平均胸径（d）、Simpson指数、林分平均高（M_H）作为备选变量。

3.2 模型参数估计

通过Kendall-Tau-b相关系数分析及考虑变量多重共线性后，将筛选出的对进界模型有显著影响的备选变量输入到计数模型中得到结果：表2、表3列出了6个基础模型（Poisson、NB、ZIP、ZINB、HP、HNB）的参数估计及评价统计量，除截距外各参数估计均在0.05水平上显著。此外，除了Poisson模型和NB模型只包含计数部分外，其余模型均包含两部分，即离散部分和零部分。在离散部分，每公顷进界株数（N_n）与林分拥挤度（K）和Simpson指数呈正相关，与林分算数平均胸径（d）呈负相关；在零部分与林分平均高（M_H）呈正相关，与林分拥挤度（K）呈负相关。

3.3 模型评价与检验

通过表2、表3模型评价结果可以看出：NB模型（AIC为695.350、BIC为707.900、Loglike为-342.670）的拟合效果优于Poisson模型（AIC为4 036.460、BIC为4 046.510、Loglike为-2 014.230）；ZINB模型（AIC为623.990、BIC为644.080、Loglike为-304.000）的拟合效果要优于ZIP模型（AIC为1 890.080、BIC为1 907.650、Loglike为-938.040）；HNB模型（AIC为624.400、BIC为644.480、Loglike为-304.200）的拟合效果要优于HP模型（AIC为1 890.080、BIC为1 907.660、Loglike为-304.200）。当模型之间的AIC值差异大于7时，说明模型间存在显著差异^［31］，同时发现ZINB模型及HNB模型拟合效果显著优于其他4个模型，且二者相差不多。故采用Voung检验进一步分析，通过表4可以看出，NB模型及其复合模型（ZINB、HNB）在Vuong检验中，拟合效果较Poisson模型及其复合模型（ZIP、HP）更好，同时ZINB模型与HNB模型的Vuong检验结果显示ZINB模型略优于HNB模型。此外，通过图3也可以看出，ZINB模型为拟合效果最接近实测值的模型。采用十折交叉法对进界模型进行检验，其中RMSE及MAE越小，则模型精度越高，通过检验结果（表5）可发现，ZINB模型为预测精度最高的模型。故ZINB模型为预测帽儿山地区天然林林分进界的最优模型。

3.4 进界数量影响因素分析

通过层次分割法获得不同自变量对最优模型（ZINB）的贡献率结果见表6，由表6表明，对帽儿山天然林更新数量相对重要的变量在计数部分由大到小依次为生物多样性因子Simpson指数（57.64％）、林分算数平均胸径（22.97%）、林分拥挤度（19.39%）；在零部分由大到小依次为林分平均高（56.37％）、林分拥挤度（43.63%）。基于模型参数估计结果，确定各协变量的参数取值，通过对林分拥挤度（K）、林分平均高（M_H）和林分算数平均胸径（d）的梯度取值，分析所选各因子对进界数量的影响及进界概率的影响，具体做法是对目标变量取等差梯度值，其他非目标协变量都取均值。图4表明，相同Simpson指数下，林分每公顷进界株数（N_n）随林分拥挤度（K）的增加而增加，随林分算数平均胸径（d）及林分平均高（M_H）的增加而减少。在相同的林分平均高（M_H）下，进界概率随林分拥挤度（K）的增加而增加。

4 讨论

林分进界数量影响着森林演替的进程，在森林动态变化过程中起着至关重要的作用。通过建立模型和分析影响森林更新的因素及规律，可以为提升森林资源的整体质量创造条件。本研究基于帽儿山天然林清查数据，通过Kencall-Tau-b相关性分析及考虑变量多重共线性，最终选择林分拥挤度（K）、林分算数平均胸径（d）、Simpson指数、林分平均高（M_H）作为变量构建进界模型。

4.1 进界株数广义线性模型构建

在本研究中，采用6种模型进行拟合对比，通过AIC、BIC、Loglike及Vuong检验发现，在对帽儿山地区林分进界拟合时，NB模型及其复合模型（ZINB、HNB）优于Poisson模型及其复合模型（ZIP、HP）。这是由于进界数据过于离散（图1），而Poisson模型中要求期望和方差相等，但进界数据过于离散不符合这一条件。而负二项模型的方差大于期望，能够更好地解释数据的异质性^［32］。ZINB模型和HNB模型优于NB模型的现象在Cunningham等^［33］的研究中有出现过，这与本研究中呈现的结果相同。这是由于数据存在较多的零值，与负二项模型相比，负二项复合模型（ZINB模型及HNB模型）运用了2个独立的回归方程来拟合整个林分进界的过程，在拟合零数据的部分表现更好。通过ZINB模型与HNB模型的Vuong检验发现，ZINB模型略优于HNB模型，这可能与二者零数据的来源有关，零膨胀数据假设零值来源于2个过程：一个过程产生结构性零值，另一个过程产生采样性零值。这样，模型可以更准确地捕捉到不同类型的零值，从而提高拟合效果；而Hurdle模型则假设所有零值都来自同一个过程，并且用一个二元模型来决定是否为零，然后用截断的计数分布来建模非零值。这种方法在处理结构性零值和采样性零值的区别时可能不够灵活。因此，处理拥有大量零数据的数据集时，零膨胀模型的表现更好^［34］。综上，ZINB模型是拟合进界数量最优的模型。雷渊才等^［5］以长白落叶松为研究对象建立进界模型的研究中，ZINB模型同样为预测精度最高的模型，与本文的研究结果相同。

4.2 进界株数影响因素分析

进界数量是评价林分天然更新能力的主要标准，探究林分进界的影响因素及规律，可为森林资源管理提供科学依据。生物多样性因子是影响林分进界的主要因素之一，通过层次分割结果可以发现，在模型的计数部分Simpson指数对林分每公顷进界株数（N_n）有显著影响，且贡献度最高。Simpson指数是反映林分中物种多样性的一个指标，其与林分每公顷进界株数（N_n）呈显著正相关，表明随Simpson指数增加林分物种多样性增加，在多样性高的林分中，不同树种具有不同的生态位和功能特征，能够更高效地利用有限的资源（如光、水和养分等），使更多的树种有了进界机会，这可能是导致林分每公顷进界株数（N_n）增加的原因。在肖晨等^［35］的研究中也发现了生物多样性的增加有利于林分的更新，这一研究结果和本文相同。通过表6，可以看到Simpson指数在模型中贡献度最高，通过促进资源高效利用、生态位分化、树木之间的正相互作用等多种机制，显著影响了林分的进界数量。相比之下，林分拥挤度（K）虽然可以反映林分内的竞争强度，但并不能全面解释进界数量的变化。因此，进界数量更多地由Simpson指数影响，而不是由林分拥挤度（K）决定的。

林分因子同样是影响林分进界的主要因素，对林分进界产生直接影响。通过模型拟合结果可以看出林分拥挤度（K）与林分每公顷进界株数（N_n）呈显著正相关。在惠刚盈等^［21］的研究中发现随K增大，林分逐渐稀疏，林下竞争减少，光照增加，这些条件可能使林分进界增加。与之相反，同为林分因子的林分算数平均胸径（d）则与每公顷进界株数（N_n）呈显著负相关，这可能是在本研究中的帽儿山天然林处于近熟林或成熟林，随林分算数平均胸径（d）的增大，林分年龄逐渐增大，林分更趋于稳定，导致进界减少。有研究认为，林分平均高（M_H）的增加会导致林分更新密度的减少^［36］，这与本研究的结果相类似。林分平均高（M_H）的增加导致林下阳光变少，竞争增加，从而导致了林分进界的减少。通过分析各影响因子对帽儿山天然林进界数量的影响，有助于更好地保护和利用当地森林资源，更加科学地管理帽儿山地区天然林更新。

5 结论

本研究以帽儿山地区天然林为研究对象，基于林分因子和生物多样性因子，通过Kendall-Tau相关性分析发现林分拥挤度（K）、林分算数平均胸径（d）、Simpson指数及林分平均高（M_H）是影响帽儿山地区天然林林分进界数量的重要因子。采用Poisson模型、负二项模型（NB）、零膨胀模型及Hurdle模型构建进界模型。对比之下，零膨胀负二项模型（ZINB）的模型精度最高，拟合效果最好，可以准确地预测帽儿山地区天然林的更新情况。此外，通过层次分割法分析所选因子贡献度，发现生物多样性因子Simpson指数和林分平均高（M_H）分别在计数部分和零部分对最优进界模型（ZINB模型）的贡献率最高。综上所述，通过所构建进界模型及分析林分变量和生物多样性变量对进界的影响，能够更好地了解在不同环境条件下林分的进界情况，为帽儿山地区天然林更新管理提供科学依据。

参考文献

[1] GRISCOM B W， BUSCH J， COOK-PATTON S C，et al.National mitigation potential from natural climate solutions in the tropics［J］.Philosophical Transactions of the Royal Society B，2020，375（1794）：20190126.

[2] 李武赫.大兴安岭东部落叶松-白桦混交林进界生长模型［D］.哈尔滨：东北林业大学，2021.

[3] LI W H.Ingrowth model of Larix gmelinii-Betula platyphylla mixed forest in eastern Daxing'an Mountains［D］.Harbin：Northeast Forestry University，2021.

[4] 韩有志，王政权.森林更新与空间异质性［J］.应用生态学报，2002（5）：615-619.

[5] HAN Y Z， WANG Z Q.Spatial heterogeneity and forest regeneration［J］.Chinese Journal of Applied Ecology，2002（5）：615-619.

[6] 阎玺羽，张广奇，臧丽鹏，等.基于2012—2021年英文文献计量学的森林幼苗更新研究进展分析［J］.山地农业生物学报，2024，43（4）：23-30.

[7] YAN X Y， ZHANG G Q， ZANG L P，et al.Based on the 2012—2021 English literature metrology analysis of the forest seedling research progress update［J］.Journal of Mountain Agriculture and Biology，2024，43（4）：23-30.

[8] 雷渊才，张雄清.长白落叶松林分进界模型的研究［J］.林业科学研究，2013，26（5）：554-561.

[9] LEI Y C，ZHANG X Q，Tree recruitment model of Larix olgensis ［J］.Forest Research，2013，26（5）：554-561.

[10] ANDREASSEN K.Development and yield in selection forest［J］.Meddelelser fra Skogforsk，1994，47（5）：1-37.

[11] BUONGIORNO J， MICHIE B R.A matrix model of uneven-aged forest management［J］.Forest Science，1980，26（4）：609-625.

[12] GRIMES R F， PEGG R E.Growth data for a spotted gum-ironbark forest in south-east Queensland［J］.Queensland，Department of Foestry，Technical Paper，1979.

[13] FORTIN M， DEBLOIS J.Modeling tree recruitment with zero-inflated models：The example of hardwood stands in southern Quebec，Canada［J］.Forest Science，2007，53（4）：529-539.

[14] 马武，雷相东，徐光，等.蒙古栎天然林生长模型的研究——Ⅳ.进界生长模型［J］.西北农林科技大学学报（自然科学版），2015，43（5）：58-64.

[15] MA W， LEI X D， XU G，et al.Growth models for natural Mongolian oak forests——Ⅳ.recruitment model［J］.Journal of Northwest A&F University（Natural Science Edition），2015，43（5）：58-64.

[16] ZELL J， ROHNER B， THÜRIG E，et al.Modeling ingrowth for empirical forest prediction systems［J］.Forest Ecology and Management，2019，433：771-779.

[17] 何江，覃林.气候敏感的杉木天然林林分进界模型［J］.林业科学，1-15［2024-11-29］.

[18] HE J， TAN L.Climate-sensitive tree recruitment model for natural Cunninghamia lanceolata forests［J］.Scientia Silvae Sinicae，1-15［2024-11-29］.

[19] 郑嵘，肖淑萍，杨小纯，等.峡江县人工湿地松中幼龄林进界胸径预估模型［J］.江西林业科技，2004（2）：15-16.

[20] ZHENG R， XIAO S P， YANG X C，et al.Estimate model of ingrowth diameter for the middle and young artifical forest of Pinus elliottii in Xiajiang County［J］.Jiangxi Forestry Science and Technology，2004（2）：15-16.

[21] 陈科屹，张会儒，张博，等.基于地理加权回归的天然次生林进界木空间分布模拟［J］.北京林业大学学报，2021，43（2）：1-9.

[22] CHEN K Q， ZHANG H R， ZHANG B，et al.Spatial distribution simulation of recruitment trees of natural secondary forest based on geographically weighted regression［J］.Journal of Beijing Forestry University，2021，43（2）：1-9.

[23] RUSSELL M B， WESTFALL J A， WOODALL C W.Modeling browse impacts on sapling and tree recruitment across forests in the northern United States［J］.Canadian Journal of Forest Research，2017，47（11）：1474-1481.

[24] SHI H， ZHOU Q， HE R，et al.Temperature sensitivity of tree recruitment at alpine treelines increases along latitudinal gradients［J］.Ecography，2023，2023（9）：e06700.

[25] KÄBER Y， MEYER P， STILLHARD J，et al.Tree recruitment is determined by stand structure and shade tolerance with uncertain role of climate and water relations［J］.Ecology and Evolution，2021，11（17）：12182-12203.

[26] 李春明，李利学.基于零膨胀模型及混合效应模型相结合的蒙古栎林林木进界模拟研究［J］.北京林业大学学报，2020，42（6）：59-67.

[27] LI C M， LI L X.Simulating study on tree recruitment of Quercus mongolica based on zero-inflated model and mixed effect model methods［J］.Journal of Beijing Forestry University，2020，42（6）：59-67.

[28] GUO Y， LI Z， LIU P，et al.Modeling correlation and heterogeneity in crash rates by collision types using full Bayesian random parameters multivariate Tobit model［J］.Accident Analysis and Prevention，2019，128：164-174.

[29] QUERO M G， GUICCIARDI S， PENNA P，et al.Increasing trends in faecal pollution revealed over a decade in the central Adriatic Sea （Italy）［J］.Water Research，2024，262：122083.

[30] 惠刚盈，张连金，胡艳波，等.林分拥挤度及其应用［J］.北京林业大学学报，2016，38（10）：1-6.

[31] HUI G Y， ZHANG L J， HU Y B，et al.Stand crowding degree and its application［J］.Journal of Beijing Forestry University，2016，38（10）：1-6.

[32] KENDALL M G.A new measure of rank correlation［J］.Biometrika，1938，30（1/2）：81-93.

[33] PERDECK A C.Poisson regression as a flexible alternative in the analysis of ring-recovery data［J］.European Union for Bird Ringing Newsletter，1998，2：37-39.

[34] EVANS M， HASTINGS N， PEACOCK B.Statistical distributions［M］.New York：John Wiley，2000.

[35] LAMBERT D.Zero-Inflated poisson regression，with an application to defects in manufacturing［J］.Technometrics，2012，34（1）：1-14.

[36] ESKELSON B N I， TEMESGEN H， BARRETT T M.Estimating cavity tree and snag abundance using negative binomial regression models and nearest neighbor imputation methods［J］.Canadian Journal of Forest Research，2009，39（9）：1749-1765.

[37] MULLAHY J.Specification and testing of some modified count data models［J］.Journal of Econometrics，1986，33（3）：341-365.

[38] 郭福涛，胡海清，金森，等.基于负二项和零膨胀负二项回归模型的大兴安岭地区雷击火与气象因素的关系［J］.植物生态学报，2010，34（5）：571-577.

[39] GUO F T， HU H Q， JIN S，et al.Relationship between forest lighting fire occurrence and weather factors in Daxing'an Mountains based on negative binomial model and zero-inflated negative binomial models［J］.Chinese Journal of Plant Ecology，2010，34（5）：571-577.

[40] LAI J， ZOU Y， ZHANG S，et al.glmm.hp：an R package for computing individual effect of predictors in generalized linear mixed models［J］.Journal of Plant Ecology，2022，15（6）：1302-1307.

[41] LAI J， ZHU W， CUI D，et al.Extension of the glmm.hp package to zero-inflated generalized linear mixed models and multiple regression［J］.Journal of Plant Ecology，2023，16（6）：rtad038.

[42] BARRY C S， WELSH A.Generalized additive modelling and zero inflated count data［J］.Ecological Modelling，2002，157（2/3）：179-188.

[43] YAACOB W， FAIROS W Y W， LAZIM M A，et al.A practical approach in modelling count data［C］.Proceedings of the Regional Conference on Statistical Sciences.Malaysia，2010：176-183.

[44] CUNNINGHAM B R， LINDENMAYER B D.Modeling count data of rare species：Some statistical issues［J］.Ecology，2005，86（5）：1135-1142.

[45] FENG C X.A comparison of zero-inflated and hurdle models for modeling zero-inflated count data［J］.Journal of Statistical Distributions and Applications，2021， 8（1）：8.

[46] 肖晨，田栋元，马榕，等.兴安落叶松天然林更新数量相容性预测模型［J］.应用生态学报，2023，34（9）：2345-2354.

[47] XIAO C， TIAN D Y， MA R，et al.Compatibility predictive model for regeneration quantities of Larix gmelinii natural forest in Daxing’anling Mountains，China［J］.Chinese Journal of Applied Ecology，2023，34（9）：2345-2354.

[48] 田栋元.大兴安岭地区主要林分类型更新空间格局及数量模型［D］.哈尔滨：东北林业大学，2022.

[49] TIAN D Y.Spatial pattern and quantitative model for the regeneration of major forest types in the Daxing'anling region［D］.Harbin：Northeast Forestry University，2022.

基金资助

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

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