树木死亡是森林树木生长过程中极其重要的组成部分之一, 影响着整体林分动态的变化；它是一个复杂的自然生态过程, 不仅受自身条件的影响, 同时也受外界环境条件的影响^[1—2]。而树木死亡可以描述为规律性和非规律性死亡两类；规律性死亡主要是由于树木在发育过程中对生存资源所产生的竞争而导致的树木死亡, 一般发生在局部范围内；非规律性死亡则主要是由于人为或自然灾害干扰(干旱、火灾、病虫害等)所导致的树木死亡, 一般发生在较大范围内, 如大面积的火灾、病虫害爆发会改变森林的演替进程^[3—4]。

近年来, 关于树木死亡模型的研究主要集中在单木死亡模型和林分死亡模型两类上。单木死亡模型是通过计算每棵树的死亡概率, 然后经过阈值分析, 判定树木死亡与否^[5]。大多数学者采用Logistic回归方法构建该模型, 分析树木存活状态与其调查因子之间的关系^[6—9]。虽然, 单木死亡模型的预测效果不错, 但其存在的问题在于概率阈值的选择尚无公认的标准^[10]。而林分死亡模型则是针对林分死亡情况进行相关建模与分析^[11]。从林分死亡株数的数据结构来看, 它具有计数数据的特点；而且数据中存在不少的零值, 可能是由于部分调查样地上并没有发生树木死亡所致^[12]。其中, 关于计数模型在林木研究上的应用主要集中在森林火灾方面^[13—16], 而在林分死亡建模研究方面的应用较少。如张雄清等^[17]利用20块落叶松重复调查的固定样地数据, 通过添加林分因子, 然后应用Poisson模型、NB模型、ZIP模型和HP模型方法构建林分死亡模型, 发现林木死亡株数受林分年龄、林分密度和林分断面积的影响, 且HNB模型的拟合效果最好。韩培武^[18]利用49块落叶松重复调查的固定样地数据, 以ZINB模型为基础模型, 通过添加林分因子和气候因子构建林分死亡模型；发现林木死亡株数与林分断面积、胸径、温度和降水有关, 且加入气候变量后的ZINB模型的拟合效果更好。刘军等^[19]利用1793块杉木样地数据, 通过添加立地、林分和气候因子, 以ZIP模型和HP模型为基础模型对杉木林的林分死亡进行研究, 发现林木死亡株数与海拔、林分年龄和密度显著相关, 且ZINB模型为最佳的林分死亡模型。

此外, 构建林分死亡模型时, 考虑样地的随机效应问题不容忽视^[20], 可以通过在基础计数模型上构建混合效应模型, 进一步提高模型的模拟精度^[21]。然而, 目前在林分死亡建模研究中考虑样地随机效应的有关报道寥寥无几。国外就Zhang等^[22]基于样地的随机效应, 利用Poisson模型、ZIP模型和HP模型方法研究林分死亡；发现林木死亡株数与直径和优势高有关, 且栅栏混合效应模型为最优的林分死亡模型。而国内只有李春明^[23]和周潇^[24]在研究林分死亡时考虑了样地的随机效应；前者的研究指出林木死亡株数受株数、林分断面积和温度的影响, 且含随机效应的ZINB模型为最佳林分死亡模型；后者的研究指出林木死亡株数受林分断面积、平均直径和林分优势高的影响, 且含随机效应的HP模型为最优林分死亡模型。所以, 在以往林分死亡建模研究中, 是否考虑样地的随机效应, 研究者们得出的结论也不尽相同。

综上所述, 在以往林分死亡模型的研究中, 关于林分死亡模型的建模方法, 主要还是集中在不考虑随机效应的基础计数模型方法上, 而加入随机效应的计数模型方法明显能够提高模型的拟合精度, 但目前在林分死亡模型的研究应用上鲜少。其次, 有关林分死亡模型的研究中, 混合效应模型是以泊松分布为主的计数模型(Poisson模型、ZIP模型和HP模型)或是以负二项分布为主的计数模型(NB模型、ZINB模型和HNB模型)基础上加入随机效应进行建模, 然后与相对应的无随机效应的基础计数模型进行比较。以此为切入点, 考虑利用常见的6种基础计数模型(Poisson模型、NB模型、ZIP模型、ZINB模型、HP模型和HNB模型), 在此基础上添加随机效应, 构建相应的混合效应模型, 然后与相对应的无随机效应的基础计数模型进行比较, 以便筛选出最优的模型, 不失为一种全面的建模方法。

为此, 本研究基于美国德克萨斯州东部森林的四次清查周期(第七、八、九和十周期)的固定样地数据, 将其按4:1的比例进行随机抽样划分成两部分, 即训练集(80%)作为模型的建模数据, 验证集(20%)作为模型的验证数据。基于模型的建模数据, 以林木死亡株数作为模型的因变量, 立地因子、林分因子和气候因子作为模型的自变量, 利用Poisson模型、NB模型、ZIP模型、ZINB模型、HP模型和HNB模型6种基础计数模型和加入样地随机效应的混合效应模型方法(Poisson-mixed模型、NB-mixed模型、ZIP-mixed模型、ZINB-mixed模型、HP-mixed模型和HNB-mixed模型)构建林分水平死亡模型。使用方差膨胀因子检验法排除各因子间的多重共线性, 采用AIC、BIC和-2logL三种评价指标筛选出拟合效果最优的模型。然后, 使用验证数据进行模型的检验, 采用均方根误差和平均绝对误差2种指标评估模型的预测效果。本研究拟解决以下问题：(1)通过比较各模型的拟合精度, 德州东部森林的最优林分死亡模型是怎样的？(2)加入样地的随机效应后, 模型的拟合精度是否提高？(3) 立地因子、林分因子和气候因子对德州东部森林的树木死亡影响如何？

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

研究区域为美国德克萨斯州东部森林(图 1)。该区域气候属于中纬度亚热带气候, 夏季炎热, 冬季温和；年平均气温大约在17—21℃的范围内变化, 年平均降水量大约在680—1700mm之间, 无霜期约为220—365d。最高海拔约为200m, 地势较为平坦, 土壤有淋溶土、壤土、粘土等类型；森林树种多以针叶和阔叶树种为主。

1.2 数据来源

FIA(Forest Inventory and Analysis)是美国林业局对全美森林资源开展的清查项目。德州东部FIA项目自1935年开始到1992年已完成6轮清查。2000年美国统一森林清查与分析的标准后, 完成了第七周期(2001—2003)的森林资源清查。自2003年开始, 美国对FIA的调查体系进行了调整, 将定期清查改为每年调查一次约20%的固定样地, 5年完成一轮清查；目前已完成第八周期(2004—2008)、第九周期(2009—2013)和第十周期(2014—2018)的森林资源清查^[25]。FIA数据可供研究人员公开使用^[26], 其全部数据可在美国林业局官方网站(http://fiatools.fs.fed.us)免费下载获取。

FIA极其重视数据的质量控制, 不仅对野外所调查的数据进行抽查和检验, 而且对录入计算机的野外调查数据也进行了实时检查；最后将准确无误的相关数据纳入有关的出版资料中, 所以其数据详实可靠^[27]。该项目的抽样则是使用三阶抽样方法；第一阶是利用遥感技术获取卫星图像进行林地和非林地的划分；第二阶则是进行地面样地调查, 调查的主要内容包含每木调查、生长、死亡、立地、林分情况等因子；第三阶则是进一步扩增了第二阶的样地调查内容, 如土壤、枯枝落叶、树冠、苔藓群落等调查因子^[28]。而本研究使用的样地调查数据主要源于第二阶调查；样地设计如图 2所示, 一个标准的FIA样地是由4组样圆构成；每组样圆又由3个半径不同的样圆构成, 分别为小样地(半径为2.07m)、子样地(半径为7.32m) 和大样地(半径为17.95m)；小样地是测量胸径小于12.7cm的树木, 子样地是测量胸径大于且等于12.7cm的树木, 大样地则是测量更大的树木(胸径>61cm或>76cm)^[29]；而1号子样地处于中心, 其他三个子样地(2、3、4)则位于距1号子样地(方向角分别是360°、120°和240°)的36.58m处。

本研究所使用的数据来源于德州东部森林四个清查周期(第七、八、九和十周期)的四次清查数据, 调查样地共计1798个, 每个样地系统的面积为0.0672hm², 调查样地的总面积为120.8256 hm²。设置筛选条件：(1)胸径≥12.7cm的树木；(2) 为了避免死树重复统计的情况, 选择各周期死树时以本周期为死树, 上一周期为活树的条件筛选。然后计算林分密度、树木基面积和林分平方平均胸径相关变量的值, 计算公式见(1)、(2)和(3)。最终选出的样地数为5249个；将其进行随机抽样, 则有4200个样地作为模拟样本(80%)用于模型的构建, 1049个样地作为验证样本(20%)用于模型的验证。林分死亡株数直方图见图 3。由图 3可以看出, 数据结构存在大量的零值, 呈现较为离散的状态。

(1)

(2)

(3)

式中, SD表示林分密度(株/hm²), N表示样地上所有树木的株数(株), S表示样地面积(hm²)；BA表示树木基面积(m²/hm²), B表示单位面积上所有树木的胸高断面积(m²)之和; SSMDBH表示林分平方平均胸径(cm), d_i²表示第i株树的胸径的平方。

2 研究方法 2.1 模型自变量的选取

林木死亡主要受立地、林分和气候因子的影响, 所以本研究将从立地、林分和气候三个方面选取模型的自变量。

2.1.1 立地因子

选取的立地因子为海拔(Elevation)、坡度(Slope)和坡向(Aspect)三个指标。这些因子影响着树木生长过程中需要的水热条件和土壤条件, 推动着森林资源的林分发展与演替^[30—32]。

2.1.2 林分因子

林木死亡主要受林分因子的影响, 而林分因子可以用林分年龄、林分密度、林分平方平均胸径和树木基面积等因子来量化^[33—35]。因此, 选择林分年龄、林分密度、林分平方平均胸径和树木基面积四个指标作为林分因子(表 1)。

2.1.3 气候因子

影响林木死亡的外在因素主要是气候, 气候因子通常由年平均温度、夏季平均温度、年平均降水量、夏季平均降水量、标准化降水蒸发指数(SPEI)和干旱长度等因子构成^[36—38]。因此, 选取年平均温度、夏季平均温度、年平均降水量、夏季平均降水量、SPEI和干旱长度六个指标作为气候因子(表 2)。而使用的气候数据来源于PRISM中提取的德州东部森林的1798个固定样地的气候数据, 并通过月平均温度和月平均降水量的数据计算得出所需的SPEI。

2.2 模型的选择

林分死亡株数的数据结构存在离散和零值较多的特点, 适用于计数模型进行拟合。而计数模型属于广义线性模型, 一般包括Poisson模型、NB模型、ZIP模型、ZINB模型、HP模型和HNB模型^[39—41]。

2.2.1 Poisson模型

泊松模型是用于分析计数类数据最基本的方法, 其概率质量函数(PMF)为：

(4)

式中, y_i表示随机变量, λ_i表示泊松分布的期望。

2.2.2 负二项(negative binomial, NB)模型

负二项模型和泊松模型的不同之处在于多了额外的离散参数来解释数据的异质性, 其概率质量函数(PMF)为：

(5)

式中, y_i表示随机变量, Γ表示伽马函数, θ表示离散参数。

2.2.3 零膨胀(zero-inflated)模型

当计数数据存在零膨胀现象时, 常用零膨胀模型进行处理。零膨胀模型主要由零部分和计数部分组成, 前者常用Logit模型进行拟合, 后者常用泊松模型或负二项模型进行拟合。因此, 零膨胀模型主要分为零膨胀泊松(zero-inflated Poisson, ZIP)模型和零膨胀负二项(zero-inflated negative binomial, ZINB)模型。零膨胀泊松(ZIP)模型和零膨胀负二项(ZINB)模型的概率质量函数(PMF)分别为：

(6)

(7)

式中, p表示零部分的概率。

2.2.4 Hurdle模型

零膨胀计数模型主要以零膨胀模型和Hurdle模型为主, 用于解决零值过多的数据问题。同零膨胀模型相似, Hurdle模型可以分为Hurdle-Poisson(HP)模型和Hurdle-NB(HNB)模型。Hurdle-Poisson模型和Hurdle-NB模型的概率质量函数(PMF)分别为：

(8)

(9)

2.2.5 连接函数

考虑在基础模型上加入样地水平的随机效应。其中, 泊松、零膨胀和Hurdle分布的连接函数, 以及负二项、零膨胀负二项和Hurdle负二项分布的连接函数分别为^[23]：

(10)

(11)

2.3 模型的评价

使用AIC、BIC和-2logL三个评价指标来比较模型之间的拟合效果。其中, AIC、BIC和-2logL的值越小, 说明模型的拟合效果越好。此外, 由于计数数据的特点不符合正态分布, 可以采用Pearson残差描述模型的拟合情况。利用验证数据检验模型时, 选择评价模型的指标为平均绝对误差(MAE)和均方根误差(RMSE)。MAE和RMSE的值越小, 说明模型的预测效果越好。

2.4 数据处理

首先, 利用Office 2016软件对数据进行处理。其次, 使用R 4.2.3软件中的MASS包、pscl包和glmmTMB包完成模型的构建；然后使用AIC()函数、BIC()函数以及modelr包等完成模型的评价。最后, 使用Origin 2022软件完成作图, 地图则由ArcGIS 10.8软件完成。FIA样地设计图则由Adobe Illustrator 2023软件完成。

3 结果与分析 3.1 多重共线性的检验

在选取模型的自变量时, 需要尽可能排除变量间的多重共线性, 因为一些变量间的多重共线性会影响模型模拟的精度。所以, 在构建模型前对选取的变量均进行了方差膨胀因子(VIF)检验。当VIF值< 10时, 说明各变量之间不存在明显的相关性。由图 4可知, 选取的13个变量之间均不存在多重共线性, 可用于模型的构建。

3.2 模型的构建

经过模型参数的检验, 剔除不显著的变量, 基础模型的模拟结果见表 3和表 4。由表 3可知, Poisson模型和NB模型的各参数均在0.05水平上显著。林木死亡株数与海拔、林分年龄、树木基面积、林分平方平均胸径、SPEI、干旱长度、年平均温度和夏季平均温度显著相关。其中, 林分年龄、树木基面积和夏季平均温度的参数估计值为正值, 呈现正相关, 说明随林分年龄的增加、树木基面积的增大和夏季平均温度的升高, 林木死亡株数增加。而海拔、林分平方平均胸径、SPEI、干旱长度和年平均温度的参数估计值为负值, 呈现负相关, 说明随海拔和年平均温度的升高、胸径的增大以及SPEI和干旱长度的增加, 林木死亡株数减少。

由表 4可知, 各模型的参数均在0.05水平上显著, 林木死亡株数与树木基面积、林分平方平均胸径、干旱长度、夏季平均温度、夏季平均降雨量、坡度、林分年龄和年平均温度显著相关。模型的计数部分结果显示, 树木基面积和夏季平均温度的参数估计值为正值, 呈现正相关, 说明随树木基面积的增大和夏季平均温度的升高, 林木死亡株数增加。而林分平方平均胸径、干旱长度和夏季平均降雨量的参数估计值为负值, 呈现负相关, 说明随胸径的增大、干旱长度和夏季平均降雨量的增加, 林木死亡株数减少。而模型的零部分结果显示, 林木死亡株数与坡度和林分年龄呈现正相关, 与年平均温度呈现负相关；说明随坡度和林分年龄的增加, 林木死亡株数增加, 随年平均温度的增加, 林木死亡株数减少。

在基础模型上考虑样地的随机效应, 各模型的模拟结果见表 5和表 6。由表 5可知, 泊松混合效应模型和负二项混合效应模型的各参数均在0.05水平上显著。林木死亡株数与海拔、林分年龄、树木基面积、林分平方平均胸径、SPEI、干旱长度、年平均温度和夏季平均温度显著相关, 与没加随机效应的泊松模型和负二项模型的模拟结果较为相近。由表 6可知, 零膨胀混合效应模型、零膨胀负二项混合效应模型、栅栏混合效应模型和栅栏负二项混合效应模型的各参数均在0.05水平上显著。林木死亡株数与树木基面积、林分平方平均胸径、SPEI、干旱长度、夏季平均温度、夏季平均降雨量、坡度、林分年龄、林分密度和年平均温度显著相关。其中, 与没有加入随机效应的ZIP模型、ZINB模型、HP模型和HNB模型的模拟结果有所差异, 模型的计数部分结果显示, 林木死亡株数与SPEI呈现显著负相关；而模型的零部分结果显示, 林木死亡株数与林分密度呈现显著负相关, 说明林木死亡株数随着林分密度的增大而减少。

图 5给出了模型的Pearson残差图。由图 5可知, 与没有加入随机效应的基础模型相比, 加入随机效应的各模型的残差波动范围均有所缩小, 说明加入随机效应的模型模拟效果更好。其中, Poisson模型的残差波动范围最大, 所以模型的拟合效果相对较差。而NB-mixed模型、ZINB-mixed模型和HNB-mixed模型的残差波动范围相对较小, 大部分落在-1到4之间, 所以这三个模型比较稳定, 拟合效果相对较好一点。

结合表 7, 我们发现与没有加入随机效应的基础模型相比, 加入随机效应的各模型的AIC值、BIC值和-2logL值均有所降低, 说明加入随机效应的模型模拟精度有所提高。其中, Poisson模型的拟合效果最差, 而NB模型、ZINB模型、HNB模型、NB-mixed模型、ZINB-mixed模型和HNB-mixed模型的AIC值、BIC值和-2logL值较为相近, 但ZINB-mixed模型的AIC值、BIC值和-2logL值略小于其他模型, 说明ZINB-mixed模型的拟合效果最好。因此, 本研究选择ZINB-mixed模型作为林分水平死亡模型的最优模拟模型。

3.3 模型的验证

利用验证数据进行模型拟合, 计算得出各模型的预估值, 并与验证样本的实际值进行比较。其中, 所用的验证样本删除异常值后为992个, 结果见图 6。由图 6可以看出, 在没有加入随机效应的6个基础计数模型(Poisson模型、NB模型、ZIP模型、ZINB模型、HP模型和HNB模型)中, 各模型的林木死亡株数的预估效果没有明显的差别；而在加入随机效应的6个混合效应模型中, Poisson-mixed模型、ZIP-mixed模型和HP-mixed模型的林木死亡株数的预估效果与NB-mixed模型、ZINB-mixed模型和HNB-mixed模型有明显的差别；其中, NB-mixed模型、ZINB-mixed模型和HNB-mixed模型的林木死亡株数的预估效果与无随机效应的6个基础模型较为相近。

进一步采用平均绝对误差(MAE)和均方根误差(RMSE) 两个指标来量化各模型的预测效果。由表 8可知, 在不考虑随机效应的6个基础模型中, Poisson模型和NB模型的MAE值和RMSE值相近且略大, 预测效果相对较差；ZIP模型、ZINB模型、HP模型和HNB模型的MAE值和RMSE值较为相近, 但ZINB模型的RMSE值(3.35)和MAE(1.53)值相对略小, 所以预测效果相对较好一点；而拟合效果较好的是ZINB模型(AIC值、BIC值和-2logL值均小于其他五个基础模型), 可见在不考虑随机效应的基础模型中, ZINB模型的拟合效果与预测效果是相一致的。在考虑随机效应的6个混合效应模型中, Poisson-mixed模型的MAE值和RMSE值相对较大, 预测效果相对较差；其次是NB-mixed模型；ZIP-mixed模型、ZINB-mixed模型、HP-mixed模型和HNB-mixed模型的MAE值和RMSE值较为相近, 但HP-mixed模型的RMSE值(3.58)和MAE(1.52)值相对略小, 所以预测效果相对较好一点；而拟合效果较好的是ZINB-mixed模型(AIC值、BIC值和-2logL值均小于其他五个混合效应模型), 可见在考虑随机效应的混合效应模型中, 预测效果最好的模型与拟合效果最好的模型有所差异。综合比较, 所有模型中ZINB模型的RMSE值和MAE值均为最小, 预测效果最好；而拟合效果最好的是ZINB-mixed模型, 可见预测效果最好的模型与拟合效果最好的模型并不一致。

4 讨论 4.1 影响树木死亡的因子

本研究发现, 德州东部森林的林木死亡株数与立地因子、林分因子和气候因子有关。立地因子方面, 林木死亡株数主要受海拔和坡度的影响。其中, 海拔与林木死亡株数呈现显著的负效应, 与刘军等^[19]在研究林木死亡时, 得出的海拔与林木死亡株数呈现正效应的结果相反；而周泽宇等^[42]在研究林木死亡株数时, 得出海拔与林木死亡株数的影响并不显著, 与本研究结果不一致；这可能是因为海拔影响着光照、水分及温度等的再分配^[43—44], 进而影响树木的生长与死亡；虽然随着海拔的升高, 光照、水分及温度可能会有所减少, 但是本研究选取的树是胸径大于12.7cm的成年树, 成年树可能在早期发育过程中经历生境过滤得以存活下来, 相比幼树, 成年树对其生存环境具有较强的适应性^[45]。坡度与林木死亡株数呈现显著的正效应, 这可能是因为随着坡度的增加, 地表截留的降水将会减少, 从而影响土壤的养分, 进而影响树木的生长^[46]。一些学者在研究树木死亡时, 认为树木死亡与海拔和坡度有关^[47—49]；这些与本研究结果类似, 说明海拔和坡度对树木死亡是有影响的。而李春明等^[23]的研究表明, 海拔和坡度对林木死亡的影响均不显著, 这与本研究结果相反, 可能是因为研究样地之间的海拔和坡度等因子的差异性不同所致。

林分因子方面, 林分年龄、树木基面积、林分平方平均胸径和林分密度是影响林木死亡株数的重要因子。其中, 林分年龄与林木死亡株数呈现显著的正效应, 这可能是因为随着成年树年龄的增加, 林分内竞争激烈, 导致树木无法得到足够的养分而发生死亡^[50]；这与Wang等^[51]的研究结果一致, 而与Caspersen等^[52]的研究结果不一致, 他们发现随着林分年龄的增加, 林木死亡株数减少；可能是因为林分年龄较大致使林分趋于稳定, 其林分年龄多数集中在5—250年, 而本研究的林分年龄则多数集中在5—100年。树木基面积与林木死亡株数呈现显著的正效应, 这可能是因为树木基面积在一定程度是反映了样地上树木的大小和密度^[53], 树木基面积越大, 代表林分内树木相对较多, 则对生长空间和养分的需求变大, 导致竞争而促使树木发生死亡；这与张雄清等^[17]的研究结果相同, 而与李春明等^[54]的研究结果相反, 他们发现随着树木基面积的增大, 树木死亡减少；可能是因为他们研究的树木中成熟木较多, 进入了平稳生长的阶段, 使得树木不易发生死亡。林分平方平均胸径与林木死亡株数呈现显著的负效应, 这可能是因为胸径的大小反映着竞争的强弱, 胸径较大时, 林木竞争资源的能力就强^[55], 处于优势地位, 导致越难发生死亡；这与闫明准^[56]和金凤伟等^[57]的研究结果类似, 他们发现树木死亡主要受直径的影响, 直径越小, 死亡越大。但是Li等^[58]的研究结果与本研究结果有所差异, 他们发现林木死亡株数与林分平均直径呈正相关, 产生差异的原因可能是因为研究角度的不同所致, 其是将胸径进行径阶等级划分去研究树木死亡, 而本研究是利用整体的胸径大小估测对树木死亡的影响。林分密度与林木死亡株数呈现显著的负效应, 这与Ma等^[59]的研究结果类似, 林分密度越大, 树木死亡率越低；而与一些学者^[60—62]的研究结果并不一致, 他们认为林分密度与树木死亡呈正相关；可能的原因是虽然林分密度较大时, 光照、水分和养分等可能会引起树木间激烈的竞争, 使得树木容易死亡, 但是本研究的对象是胸径大于等于12.7cm的树, 林分内的树木生长环境相对稳定, 不易发生林木死亡。

气候方面, SPEI、干旱长度、年平均温度、夏季平均温度和夏季平均降雨量是影响林木死亡株数的主要因子。其中, SPEI、干旱长度、年平均温度和夏季平均降雨量与林木死亡株数呈现显著的负效应, 这与臧颢^[63]的研究结果相似, 随温度升高和降雨增加, 林木死亡减少；可能是由于气温和降水量的增加, 改善了树木生长的环境, 降低了林分内竞争, 有利于树木的存活。但是闫明等^[64]的研究结果与本研究结果有所差异, 我们发现SPEI和干旱长度与树木存活存在负相关关系, 年降水量则成正相关关系, 可能是因为研究方法的不同和研究范围仅为德州东部的四个国家森林且包含幼树造成的。而干旱强度和夏季平均温度与林木死亡株数呈现显著的正效应, 这与Kim等^[2]的研究结果相似, 他们发现林木死亡株数与最冷月平均温度呈正相关, 可能是因为随着夏季温度的升高, 水分蒸发较快, 一方面影响着水分的供需, 导致树木缺水而亡；另一方面致使气孔关闭, 导致碳水化合物的供给不足, 进而使树木发生碳饥饿现象^[65]；但是刘军等^[19]的研究结果与本研究结果并不一致, 他们发现气候因子对林分死亡的影响并不显著, 可能是因为研究区域的气候因子间的差异性不同所致。

4.2 模型的效果

模型预测效果的评价指标选择的是较常用的平均绝对误差(MAE)和均方根误差(RMSE), 是一种准确检验模型的常见度量指标, 它是利用预测的林木死亡株数与实际的林木死亡株数之间的差异进行计算的^[18]。其中, 无论是在基础模型中还是混合效应模型中, 泊松模型和负二项模型的预测效果均较差, 可能是因为林木死亡数据中存在大量零值, 这两类模型并不太适用所致。而在混合效应模型中, 整体模型的预测效果并不太理想；可能的原因是：首先, 可能是加入随机效应后的基础模型低估了计数部分, 过度预测了零部分, 致使模型的预测效果不佳^[24]；其次, 可能是所用的验证样本量不足所致^[23]。

本研究的结果表明, 无随机效应的基础计数模型中, ZINB模型的拟合效果最好；这与刘军等^[19]的研究结果一致；而与张雄清等^[17]利用基础计数模型方法研究林分死亡模型时, 得出Hurdle-NB模型为最优模型的结果并不一致。其次, 与没有加入随机效应的基础模型相比, 加入随机效应的混合模型明显提高了整体模型的模拟精度；而含随机效应的ZINB模型的模拟效果最好, 可作为林分水平死亡模型的最优模型。李春明等^[23]基于混合效应模型和零膨胀模型方法研究林分死亡模型时, 发现加入随机效应的ZINB的模拟效果最好, 与本研究结果一致。但也一些学者得出和本研究不一致的结果, 如Zhang等^[22]通过利用泊松混合效应模型研究林分死亡模型时, 发现栅栏混合效应模型为最优林分死亡模型；周潇^[24]以华北落叶松天然林为研究对象, 考虑样地的随机效应构建林分死亡模型, 发现含随机效应的Hurdle模型为最优林分死亡模型。基于以上分析, 这些最优林分死亡模型的研究结果有所差异, 可能是研究对象、研究区域、研究方法以及林木死亡因子的选取等因素造成的。

5 结论

本研究基于美国德州东部森林连续清查数据和气候数据, 利用计数模型方法以及样地的随机效应, 构建基于立地因子、林分因子和气候因子的林分水平死亡模型, 探究影响树木死亡的重要因素, 并对所构建的模型进行精度比较, 筛选出最优的林分水平死亡模型。结果表明, 立地因子方面, 林木死亡株数与海拔呈显著的负效应, 与坡度呈显著的正效应；林分因子方面, 林木死亡株数与林分年龄和树木基面积呈现显著的正效应, 与林分平方平均胸径和林分密度呈显著的负效应；气候因子方面, 林木死亡株数与干旱长度、年平均温度和夏季平均降雨量均呈显著的负效应, 与干旱强度和夏季平均温度呈显著的正效应。其次, 在不考虑随机效应的基础计数模型中, ZINB模型的拟合效果最好；而加入随机效应后的基础计数模型, 明显提高了整体模型的拟合精度；在所有模型中, ZINB-mixed模型的模拟效果最好, 可作为林分水平死亡模型的最佳模型。本研究揭示了德州东部森林中林木死亡的影响因素；以及加入随机效应后的计数模型可以提高模型的拟合精度；筛选出了最优的林分水平死亡模型；可为有关树木死亡模型的研究提供一种可行性方法；同时, 构建的树木死亡模型有利于评估该研究区影响树木死亡的因子, 在一定程度上可以反映森林树木存活的情况, 对森林资源的监测与管理经营具有一定的参考价值。在今后的研究中, 可以考虑从分径阶的角度入手, 结合样地的随机效应, 构建林木死亡的径阶计数模型, 进一步深入研究树木死亡模型。