基于零膨胀模型和栅栏模型的赣南杉木林林分枯损模型

刘 军，潘 萍，，欧阳勋志，*，臧 颢，郭 杨，游景晖

（1.鄱阳湖流域森林生态系统保护与修复国家林业和草原局重点实验室，江西 南昌 330045；2.江西农业大学 林学院，江西 南昌 330045）

【研究意义】枯损是森林动态变化过程中的一个重要表现，也是森林生长过程中不可或缺的组成部分，其受内在竞争、人为干扰以及气候变化等因素的综合影响，是森林演替、更新的重要内容[1-2]。分析林分的枯损情况，探讨树木死亡的影响因素，有助于提升森林质量及生态效益，为造林、育林和护林等生产活动提供科学依据[3]。【前人研究进展】由于枯损模型存在研究尺度间的差异性，因此可划分为单木水平和林分水平两类。当前，国内外诸多学者对于单木水平枯损模型已进行了较多研究[4-5]。虽然单木水平枯损模型的预测结果比较精确，但是需要花费大量的时间、精力和财力对单木、林分和立地因子等进行较为全面的调查[6]，而通过分析林分枯损株数与林分因子、立地因子和气候因子等环境因子之间的关系构建林分水平枯损模型则可以较为有效地克服这一缺点[7]。考虑到林分枯损株数为典型的计数数据，因此不少学者采用泊松模型、负二项模型、零膨胀模型和栅栏模型等计数模型对林分枯损进行了相关研究，如张雄清等[7]基于吉林省汪清林业局金沟岭林场的20块面积在0.077 5～0.25 hm2的长白落叶松（Larix olgensis）样地，得出栅栏负二项模型的拟合效果优于其他模型，其林分枯损与林龄和相对植距呈负相关，与林分断面积呈正相关。然而，目前大多研究都建立在较小范围尺度上，而基于较大区域尺度对林分枯损的研究还相对较少，且研究得出的较优林分水平枯损模型以及影响林分枯损的因素也不尽相同，这可能与树种的生物学特征及研究的尺度等因素有关。【本研究切入点】杉木（Cunninghamia lanceolata）是我国亚热带地区特有的优良用材树种。具有较喜光、不耐严寒与干旱等特点，在满足我国木材需求、维持生态安全和支持天然林保护等方面发挥着至关重要的作用[8]。近年来，国内外不少学者对杉木的枯损影响因素展开了相关研究，如郎奎建等[9]基于江西省16 块0.05 hm2的杉木人工林构建枯损模型，得出林分枯损株数与林分密度和立地因子密切相关；Zhang 等[10]以江西、福建、广西和四川等省份建立的60 块0.06 hm2样地为研究对象，构建逻辑回归模型，认为竞争、林龄和林分结构是导致杉木枯损的主要因素。然而，由于全球气候多变以及经营措施的不同，在不同区域导致杉木林林分枯损的主要影响因素可能存在差异。【拟解决的关键问题】为此，本研究以江西赣南（赣州市）为研究区，以地类为纯林的杉木样地为研究对象，选取零膨胀模型和栅栏模型，构建并筛选出最优林分水平枯损模型，并分析其影响因素，以期为杉木林的有效经营管理提供参考依据。

1 材料与方法

1.1 研究区概况

赣州市（24°29′～27°09′ N，113°54′～116°38′ E）位于江西省南部，地处赣江上游，总面积3.9×104km2，占全省总面积的23.6%。地貌以山地和丘陵为主，土壤以地带性红壤、黄壤为主。气候属亚热带季风气候区，年均温为19.1～20.8 ℃，年均降水量为1 580 mm，无霜期288 d，年均日照时数1 636 h，森林覆盖率达76.2%，主要森林类型有针叶林、常绿阔叶林等[11]。其中杉木林面积70.40 万hm2，占有林地总面积的26.2%。

1.2 数据来源与处理

1.2.1 样地数据与处理 样地数据来源于赣州市2009 年森林资源二类调查，样地面积为0.08 hm2（28.28 m×28.28 m），地理坐标以其西南角为准，对所有胸径≥5 cm 的乔木进行每木检尺，并调查了坡向、坡度等立地因子及林龄、平均胸径等林分因子。参考《江西省森林资源二类调查操作细则（2009）》，基于样地数据库，剔除部分数据异常的样地，筛选出地类为纯林的杉木样地1 973 块，其样地分布见图1。按4∶1的比例随机划分为模拟样本（1 579 块样地）和验证样本（394 块样地），分别用于模型的构建和验证。通过筛选每木检尺库中检尺类型为枯立木（已枯死但未倒伏的林木）的林木得出各样地的枯损株数，其林分枯损株数频数统计见图2，由图2可知，本研究数据结构中含有大量的零值，即处于零膨胀状态，数据离散程度较大。

图1 研究区位置及样地分布Fig.1 Location of the study area and sample plots distribution

图2 林分枯损株数频数统计Fig.2 Frequency statistics of the number of stand-level mortality stems

1.2.2 气候数据 利用ClimateAP v2.30软件[12]，通过输入样地的经纬度和海拔，获取各样地2000—2009年的各类气候因子数据，并对其求取年平均值。

1.3 研究方法

1.3.1 枯损影响因子选择 根据样地调查内容及参考相关文献[6,13]，从林分、立地和气候等方面选择可能影响林分枯损的17个环境因子，其中林分因子包括林龄、平均树高、断面积和株数密度等4 个因子；立地因子包括海拔、坡向、坡位、坡度、土层厚度和腐殖质层厚度等6 个因子，且坡向、坡位等定性因子参考邵方丽等[14]的方法将其赋值量化，即用1、2、3、4、5 分别代表阴坡、半阴坡、半阳坡、阳坡和无坡向，用1、2、3、4、5、6、7分别代表山脊、上坡、中坡、下坡、山谷、平地和全坡；气候因子包括年平均温度、年平均降雨量、最暖月平均温度、最冷月平均温度、无霜期天数、前一年8 月至当年7 月的降雪量和夏季（6—8月）平均降雨量等7个因子。

1.3.2 林分水平枯损模型选择 本研究结合数据情况，选择零膨胀模型（Zero-inflated）和栅栏模型（Hurdle）构建林分水平枯损模型。

（1）零膨胀模型。零膨胀模型可以有效解决数据过度离散和零值过多等问题[15-16]，其数据结构中的所有零值既来源于Logit部分，也来源于截尾的泊松（Poisson）部分或负二项（NB）部分。模型主要由两部分组成：第一部分为逻辑斯蒂模型（Logit），第二部分为Poisson 模型或NB 模型。前者大多用于解释林分发生枯损的可能性（概率），即对应模型的零部分；后者主要用于预测林分枯损株数，即对应模型的离散部分。基于模型第二部分存在的形式差异，可分为零膨胀泊松模型（Zero-inflated-Poisson，ZIP）和零膨胀负二项模型（Zero-inflated-NB，ZINB）。其计算原理及方法详见文献[17]。

（2）栅栏模型。栅栏模型又称两部分模型，其优点及对应形式与零膨胀模型相似，但其数据结构中的所有零值只来源于Logit部分。模型主要由2个阶段构成：第一个阶段判断林分是否发生枯损，只有当林分发生枯损时，才会进入到第二个阶段；第二个阶段主要利用Poisson 分布和NB 分布预测林分枯损株数。同理，可分为栅栏泊松模型（Hurdle-Poisson，HP）和栅栏负二项模型（Hurdle-NB，HNB）。其计算原理及方法详见文献[18]。

1.3.3 模型评价 采用赤池信息准则（AIC）、贝叶斯信息准则（BIC）和-2 倍对数似然函数值（-2logL）比较各模型的拟合效果，其值越小，表明模型的拟合效果越好。利用验证样本检验模型精度时，选择平均绝对误差（MAE）和均方根误差（RMSE），其值越小，表明模型的预测效果越好。

1.3.4 数据分析方法 数据的预处理在Excel 2016软件中进行。利用R 4.1.2软件完成模型的构建和验证，首先通过R 4.1.2 的回归诊断程序包car 中的VIF 函数进行多重共线性诊断；其次采用caret 包中的createDataPartition 函数将数据拆分为训练样本及测试样本；然后利用pscl包中的zeroinfl和hurdle函数构建林分水平枯损模型；最后使用AIC函数、BIC函数和predict函数等进行模型评价。运用Origin 2019b软件进行相关图表的绘制。

2 结果与分析

2.1 多重共线性诊断

为提高模型模拟精度，对选取的17个可能影响林分枯损的环境因子，需尽可能地排除彼此间相关性较大的因子。本研究的筛选原则为，当1对因子间存在共线性时，只保留其中一个因子，优先保留有研究表明对林分枯损影响显著的因子，筛选后的诊断结果见图3。由图3可知，除最暖月平均温度、最冷月平均温度和无霜期天数外，其余14 个环境因子间均不存在多重共线性（VIF＜10），因此将其用于林分水平枯损模型的构建。具体环境因子统计见表1。

表1 环境因子统计Tab.1 Statistics of environmental factors

图3 自变量多重共线性诊断Fig.3 Multicollinearity diagnosis of independent variables

2.2 林分水平枯损模型的构建

经过模型参数的检验，剔除不显著的变量，得出零膨胀模型和栅栏模型的模拟结果见表2。由表2可知，在模型的零部分，各参数估计值均在0.001 水平上显著，表明杉木林林分是否发生枯损主要与海拔、林龄和株数密度等3个因子显著相关，而与气候因子均没有显著影响。在模型的离散部分，海拔和林龄在ZINB 模型和HNB 模型上在0.01 水平上显著，而在ZIP 模型和HP 模型上在0.001 水平上显著；株数密度均在0.001水平上显著，表明株数密度对林分枯损的影响最大。海拔和株数密度的参数估计值均为正值，说明随海拔梯度的上升和株数密度的增大，林分枯损株数增加；而林龄的参数估计值均为负值，说明随林龄的增加，枯损株数减少。

表2 林分水平枯损模型模拟结果Tab.2 The simulation results of stand-level mortality model

零膨胀模型和栅栏模型的拟合效果见表3。由表3 可知，ZIP 模型和HP 模型的AIC 值（63 048.66）、BIC 值（63 091.58）和-2logL 值（63 032.66）相一致；而ZINB 模型和HNB 模型的值相近，但ZINB 模型的AIC值（4 841.73）、BIC值（4 890.01）和-2logL值（4 823.73）均略小于HNB模型的AIC值（4 841.76）、BIC值（4 890.04）和-2logL 值（4 823.76）。由此可见，ZINB 模型的拟合效果最好，其次为HNB 模型，ZIP 模型和HP模型较差。

表3 林分水平枯损模型拟合效果Tab.3 Fitting effect of stand-level mortality model

2.3 模型验证

用验证样本进行模型验证，即使用394个样本进行模型拟合计算得出零膨胀模型和栅栏模型的林分枯损株数，并与验证样本实际发生的林分枯损株数进行比较，结果见图4。由图4 可以看出，ZIP 模型、ZINB模型、HP模型和HNB模型的差别不是十分明显。

图4 零膨胀模型和栅栏模型的预测值与实际值Fig.4 Predicted value of zero-inflated model and Hurdle model vs.observed value

进一步采用平均绝对误差（MAE）和均方根误差（RMSE）等2 个指标量化比较4 个模型的预测效果，结果见表4。由表4 可知，ZINB 模型的MAE 值（39.429 3）和RMSE 值（116.508 9）依次小于HNB 模型的MAE 值（39.440 4）和RMSE 值（116.531 5）、ZIP 模型的MAE 值（39.584 0）和RMSE 值（116.950 7）及HP 模型的MAE 值（39.584 0）和RMSE 值（116.950 8），说明ZINB 模型的预测效果优于其他模型，这也和前面的模型拟合效果相一致。综合分析，本研究得出的最佳模拟模型为零膨胀负二项模型，其模型形式及概率分布形式如下。

表4 林分水平枯损模型预测效果Tab.4 Prediction effect of stand-level mortality model

模型形式：

概率分布形式：

式（1）、（2）中：AL为海拔，AG为林龄，DOT为株数密度。

3 讨论

3.1 林分枯损影响因子

本研究发现，海拔、林龄和株数密度是影响杉木林林分枯损的重要因子。林分因子中，林龄和株数密度对杉木林林分枯损具有显著影响。随林龄的增加，枯损株数减少，这与Caspersen[19]基于密歇根州森林资源清查数据得出的研究结果一致，而与王涛等[20]通过对黑龙江省江山娇实验林场的48 块样地进行研究，得出的林分年龄大体与杂种落叶松（Larix gmelinii）人工幼龄林林分枯损株数呈正相关的结果相反，这可能是因为本研究的杉木林多为人工林，随造林年限的增长，抚育间伐等经营措施使得林分内竞争减小而促使枯损株数减少。反之，在造林初期时，由于株数密度较大，使得林木对于光照、水分等的竞争加剧而导致枯损株数增加，与Zhao等[21]的研究结果相似。

立地因子中，海拔对林分枯损具有重要影响。这可能是因为杉木在生长过程中，较喜光、喜温暖湿润的气候，随海拔梯度的上升，温度逐渐降低，不适宜杉木的生长，这与张凌宇等[22]的研究结果一致。但也有学者得出与本研究不一致的结果，如Radcliffe等[23]通过对俄亥俄州东南部的82块0.05 hm2固定样地进行研究，得出北美红栎（Quercus rubra）的枯损与坡位、土层厚度有关；郭韦韦[24]的研究表明，林下腐殖质层厚度与色木槭（Acer pictum）幼树枯损显著相关，红松（Pinus koraiensis）幼树的枯损则主要受坡度的影响。然而在本研究中，除海拔外，其他立地因子在模型中均不显著，这可能是因为杉木人工造林时多从木材生产考虑而选择在立地条件较好的地块，因此在本研究所使用的1 973 块样地中，有70%左右的样地坡位为中下坡，80%以上的样地坡度介于16°～35°，2/3 左右的样地土层厚度大于80 cm，2/3 左右的样地腐殖质层厚度大于10 cm，样地在坡位、坡度、土层厚度和腐殖质层厚度之间差异不大而导致的，其难以体现在对枯损的影响上，这与李春明等[6]利用吉林省295块蒙古栎（Quercus mongolica）样地数据得出的结果相类似。

此外，不少研究表明，气候因子与林分枯损状况有关。如Kim等[25]基于韩国国家森林资源清查数据，利用MSN曲线构建林分水平枯损模型，认为针叶林的枯损株数随最冷月平均温度的升高有增加的趋势；Zhang等[10]得出杉木人工林的枯损株数随年平均温度的升高而减少。但本研究中，气候因子对杉木林林分枯损的影响均不显著，可能是由于本研究区杉木分布区的主要气候因子间的差异较小，也有可能是由于海拔等地形因子折射出了一部分气候因子的影响，如Wu等[26]认为地形因子之间的差异会间接影响森林环境中光照和降水等气候因子的分布情况，进而影响林分的生长和枯损状况。

3.2 模型模拟效果

本研究表明，杉木林林分水平枯损模型的拟合效果和预测效果均以ZINB 模型最好，其次为HNB 模型、ZIP 模型，HP 模型最差。在模型的拟合效果上，李春明等[6]和Li 等[27]也得到了与本研究相同的结果，而张雄清等[7]利用汪清林业局金沟岭林场的长白落叶松林分数据，得出HNB 模型的拟合效果最好。在模型的预测效果上，与郭韦韦[24]对冷杉（Abies fabri）幼树枯损的研究结果一致，但也有学者得出不同的结果，如Crecente-Campo 等[28]基于西班牙巴斯克地区的森林资源清查数据，结果表明ZIP 模型在预测枯立木株数上效果最好。这些研究所得的最佳林分枯损模型存在一定差异，这可能是受研究的对象、区域、尺度以及环境因子等因素的影响[29]。

4 结论

本研究以江西赣南的杉木林为研究对象，以林分枯损株数为因变量，以林分因子、立地因子和气候因子等14个环境因子为自变量，构建杉木林林分水平枯损模型，探讨影响林分枯损的主要因素。结果表明，杉木林林分水平枯损模型的拟合效果和预测效果均以ZINB 模型最好，其次为HNB 模型、ZIP 模型，HP 模型最差；海拔、林龄和株数密度是影响杉木林林分枯损的重要因子，而气候因子对其没有显著影响。在今后的研究中，有待进一步考虑预估期间林木采伐、土壤理化性质和病虫害、人为干扰等因素展开更为深入的分析，以预测林分枯损株数变化的动态过程。

致谢：亚洲开发银行CCF（气候变化基金）项目、江西赠款项目（0229-PRC）和江西省林业科技创新专项（创新专项[2021]33号）同时对本研究给予了资助，谨致谢意！