DOI:10.13275/jki.lykxyj.2020.06.013
李春明1,2*
(1. 中国林业科学研究院资源信息研究所 北京 100091；2. 国家林业和草原局森林经营与生长模拟重点实验室 北京 100091)
摘要：[目的]通过基于广义线性混合效应模型的方法，建立枯损模型，提高枯损模拟和预测精度。[方法]以1995年在吉林省布设的295块蒙古栎永久样地作为研究对象，基于logistic 回归及广义线性混合效应模型方法，考虑单木、林分、气象及立地等因子对枯损的影响，构建单木水平枯损模型。并与传统广义线性模型方法进行比较，采用验证数据进行精度验证。[结果]表明：与固定效应模型相比，考虑样地水平随机效应后，模型的模拟精度获得极大的提高，并且达到极显著程度；在截距和胸径变量上考虑随机效应的模型模拟精度最高；对单木枯损有重要影响的因子中，林分公顷断面积、林分平均年龄及年平均夏季(5—9月)降水量（MSP ）与单木枯损概率呈负相关，大于对象木断面积、林分公顷株数、气象因子中的年平均温度（MAT ）及上一年8月至当年7月的降雪量（PAS ）与单木枯损概率呈正相关，初始单木胸
径对枯损的影响比较复杂，要和其他因子结合起来分析。[结论]甄别影响单木枯损的不同来源及重要影响因素，有助于进一步制定科学的森林经营措施。
关键词：广义线性混合效应模型；单木；枯损；蒙古栎；logistic 模型
中图分类号：S757.2 文献标志码：A 文章编号：1001-1498(2020)06-0105-09
作为生态过程中不可或缺的组成部分，林木枯损被认为是森林生态系统中极其关键，却又知之甚少的森林动态过程[1]。单木枯损模型能够对每株树木的生长状况进行模拟，诸多学者利用多种数据和方法构建了大量的单木枯损模型。传统的经验方程，如指数方程和Weibull 方程等，主要是通过建立枯损与协变量的回归方程来预测未来树木的枯损概率[2-3]。由于每株树木只包括枯损和存活两种结局，是典型的二分类数据，上述传统的经验方程在估计枯损时效果并不理想。而Logistic 回归方法对于二分类数据具有先天的优势，与传统的经验方程相比，模拟林木的枯损有较好的适用性，已经被广泛应用[4-5]。
在构建单木枯损模型时，以往基于Logistic 回归方法的枯损模型只反映枯损林木总体的平均变化规律，忽略了样地之间及林木个体之间的差异，也没有考虑数据间的时空序列相关性，势必会对模型的准确性造成影响[6]。混合效应模型由固定效应和随机效应两部分组成，既可以反映总体的平均变化趋势，又可以提供数据方差、协方差等多种信息来反映个体之间的差异。另外，在处理连续观测数据存
在的序列相关性等问题时具有其他模型无法比拟的优势[7]。因此一些学者开始尝试利用广义线性混合效应模型来构建单木枯损模型[7-9]。研究结果表明：与传统的Logistic 回归方法相比，基于广义线性混合效应模型方法的枯损模型明显提高了模拟和预测精度。
影响单木枯损的因素很多，主要包括单木因
收稿日期：2020-04-08 修回日期：2020-06-08
基金项目：国家自然科学基金面上项目“基于混合效应模型的联立方程组及概率分布模型在模拟森林生长中的方法研究”（31570625）* 通讯作者: 李春明，E-mail ：****************
林 业 科 学 研 究2020，33（6）：105-113Forest Research www.lykxyj
子、林分因子、立地因子、气象因子以及这些因子的交互作用。单木因子主要包括林木本身的生长和遗传特性。林分因子方面，主要包括林分密度、竞争指数、林分生产力及林分结构等，多数情况下被同时考虑或考虑其中几个因子[10-13]。立地因子，主要考虑坡度、坡向、海拔及其数学变化形式对枯损的影响[8]。气象因子方面，主要考虑温度、降水量、霜冻及干旱等对枯损的影响[1, 14-16]。在实际情况中，树木的枯损可能会由于多种原因同时发生，例如林分因子、立地因子及气象因子的相互作用等[17]。
蒙古栎（Quercus mongolica Fisch. ex Ledeb）是我国东北林区中主要的天然树种[18]。研究蒙古栎林的枯损和存活状态，构建蒙古栎林的单木枯损模型对于构建林分生长和收获模型具有重要意义。马武等[19]以蒙古栎天然林为研究对象，基于吉林省汪清林业局195块固定样地数据，采用二分类的Logistic回归方法建立了蒙古栎天然林的单木枯损模型，但没有考虑枯损数据存在的样地和个体之间的差异和多次测量的时间序列相关性。本研究以吉林省东南部地区的295块蒙古栎林固定样地多次观测数据为研究对象，把数据分成两部分，一部分为模拟数据（约80%），共计236块样地38 096条单木记录，一部分为验证数据（约20%），共计59块样地9 674条单木记录。基于Logistic回归及广义线性混合效应模型方法，考虑单木和林分因子，立地因子和气象因子的影响，构建单木枯损模型，并与不考虑随机效应的传统方法进行比较，然后采用验证数据进行精度验证。
1 数据来源
本研究选择吉林省295块蒙古栎样地数据，样地大小为0.06 hm2。样地中蒙古栎为优势树种，其他主要树种还包括云杉(Picea jazoensis Nakai)、白桦(Betula platyphylla Suk.)、红松(Pinus koraiensis Sieb. et Zucc.)、水曲柳(Fraxinus mandshurica Rupr.)和冷杉(Abies nephrolepis (Trautv.) Maxim)等。调查的样地因子主要包括起源、优势树种、平均林龄、平均树高、郁闭度、海拔、坡度、坡向及坡位等因子。单木因子包括对胸径大于5 cm的树木进行每木检尺，枯损木的标记和记录。外业调查结束后进行林分公顷株数、公顷断面积、大于对象木断面积等指标的计算。建模数据采用1999年和2009年观
测的林分因子数据。枯损树木的记录分别是1999年观测（没有枯损）—2004年观测（枯损）和2009年观测（没有枯损）—2014年观测（枯损）的数据。而验证数据采用的是2004年观测数据。具体的单木因子、样地因子及立地因子见表1。气象因子利用ClimateAP v2.10软件[20]，通过输入各样地的空间位置坐标和海拔，然后获取各样地每年的各类气象因子变量，最后对每5年的变量进行平均计算。气象因子分别对应的是1999—2004年和2009—2014年平均数据，具体因子见表2。
2 研究方法
2.1 主要影响因子选择
p>0.05
首先把单木、样地、立地和气象等因子加入到模型中，逐步剔除差异不显著（）以及存在多重共线性问题（方差膨胀因子 > 5）的影响因子。本研究主要的候选因子包括单木初始胸径、单木初始胸径的自然对数、单木初始胸径的平方、单木初始胸径与林分平均胸径的比值、大于对象木断面积、林分算术平均直径、林分平方平均直径、林分公顷断面积、林分公顷株数、平均年龄、郁闭度、坡度（slope）、坡向（aspect）、海拔（elevation）、SIC值（SIC = tan(slope) × cos(aspect)）、CE值（CE = cos（aspect） × ln （elevation））、年平均温度（MAT）、最暖月平均气温（MWMT）、最冷月平均气温（MCMT）、年平均降水量（MAP）、年平均夏季(5—9月)降水量（MSP）、无霜期天
数（NFFD）和上一年8月至当年7月的降雪量（PAS）等23个影响因子。
2.2 广义线性混合效应模型及评价方法
广义线性混合效应模型（GLMM）通过在线性预测部分引入随机效应参数而推广的广义线性模型（GLM）[21]。本研究基于Logistic回归的广义线性混合效应模型在截距和参数系数上考虑随机效应，利用SAS 9.4 的PROC NLMIXED 模块来估计广义线性混合效应模型的待估参数。
AIC信息准则、BIC信息准则和−2*对数似然值（−2Log L）这3个值经常被用来比较不同模型之间的拟合效果，这3个值越小，说明模型的模拟效果越好[22]。似然比卡方检验（LRT）被用来比较模型之间的差异显著程度[23]。
AUC值是ROC曲线下面的面积，能够描述Logistic回归模型的辨别能力，被广泛用来评价单木枯损模型的特异性和敏感性。AUC值越大，说
106林 业 科 学 研 究第 33 卷
明模型的特异性和敏感性越好，具体计算方法见Fielding and Bell [24]。2.3 阈值确定
在利用枯损模型来判断林木是否枯损前需要科
学合理的确定阈值。确定阈值的主要方法包括传统阈值、随机阈值、错误分类率（MCR ）最低阈值和割点阈值等[9, 25]。传统阈值通常选择0.5，该值仅适用于枯损木与存活木相近的情况，实际上枯损
表 1 样地各调查因子统计
Table 1 The statistics of plot factors
表 2 主要气象变量统计
中国武术 泰拳
Table 2 The statistical table of main climate variables
第 6 期李春明： 基于广义线性混合效应模型的蒙古栎林单木枯损建模及影响因子分析
奥德修斯
107
是小概率事件；随机阈值选择灵敏度与特异度曲线交叉时的阈值为最佳阈值；割点阈值选取林分各个样地枯损率的平均值作为最佳阈值[26]。本研究选择了随机阈值、错误分类率最低阈值和割点阈值等3种阈值确定方法，来确定最佳阈值。
3 结果与分析
通过筛选，最后单木初始胸径、大于对象木断面积、林分公顷断面积、林分公顷株数和林分平均年龄，气象因子中的年平均温度（MAT）、年平均夏季(5—9月)降水量（MSP）、上一年8月至当年7月的降雪量（PAS）等因子，被保留在基础模型中。
在Logistic模型的基础上，考虑样地水平、单木水平及两水平（样地内嵌套单木水平）的随机效应（分随机截距效应和随机系数效应两种情况）。在考虑单木水平时只有在截距上收敛，但是其AIC、BIC和−2Log L值分别是18 559、18 639和18 539，远远高于固定效应模型。在考虑两水平效应时，无论什么情况，模型均不能够收敛。在考虑样地水平随机系数效应时，3个参数以上的随机效应均不收敛，考虑2个参数时，只有两种情况收敛，第一种是同时在截距和胸径两个参数上考虑随机效应，第二种是同时在截距和大于对象木断面积时，模型能够收敛。具体模拟和计算的AIC、BIC、−2Log L、LRT、AUC值和参数的估计值见表3。表3中，M1为传统Logistic回归模型，M2为在截距上考虑随机效应的混合效应模型，M3为在截距和胸径变量上同时考虑随机效应的混合效应模型，M4为在截距和大于对象木断面积变量上同时考虑随机效应的混合效应模型。
表 3 基于广义线性模型方法和广义线性混合效应模型方法的单木枯损模型模拟结果
Table 3 The result of parameter estimate based on generalized linear model and
generalized linear mixed model in tree mortality美洲虎攻击机
表3的模拟结果表明，当模型考虑了样地的随机效应后，AIC、BIC、−2Log L值均比没考虑样地的随机
效应值要小，说明在构建单木枯损模型时考虑随机效应的模拟效果要优于不考虑随机效应的模型。利用LRT对不同的模型进行方差分析，考虑样地随机效应后，与固定效应模型相比，差异均达
108林 业 科 学 研 究第 33 卷
p <0.001p <0.001到极显著水平（）；当同时考虑截距和随机系数效应后比单独考虑随机截距效应，模型的模拟效果要好，差异也达到极显著（）。在胸径和截距上考虑随机效应的AIC 、BIC 、−2Log L 等3个指标值最小，说明M3的模拟效果最好。在本研究中，M3的AUC 值等于0.852 > M4（0.850）>M2（0.838）> M1（0.779），也说明M3模型是构建蒙古栎单木枯损模型的最优模型。3.1 单木和林分因子对枯损的影响
从表3可知，M1和M2中单木初始胸径值与枯损概率呈正比，而M3和M4中单木初始胸径值与枯损概率呈反比，说明单纯一个胸径因子不足以表述单个林木在林分中的地位和状态。大于对象木断面积与枯损概率呈正比，说明大于对象木断面积值越小，处于竞争的优势，不容易发生枯损，这与林分中现实情况相符。林分公顷断面积与单木枯损概率呈反比，反映了在给定同样公顷株数的林分条件下，随着林分公顷断面积增加，林分内的树木胸径更大，处于竞争的优势，则枯损概率降低。林分平均年龄与单木枯损概率呈反比，说明枯损容易在幼龄林中发生，而随着林分的生长和成熟，林木胸径变大，枯损概率降低，林木也不易发生枯损。林分公顷株数与单木枯损概率呈正比，说明在密度大的林分，林木之间对于光、热、水肥等营养的竞争更加激烈，林木容易发生枯损。3.2 气象因子对枯损的影响
在表3的模拟结果中，年平均温度（MAT ）、年平均夏季(5—9月)降水量（MSP ）、上一年8月至当年7月的降雪量（PAS ）对枯损有明显的影响。MAT 与单木枯损概率呈正相关，同样PAS 也与单木枯损概率呈正相关。而MSP 与单木枯损概率呈负相关。气象因子的影响说明，气候变化引起的干旱或霜冻对单木枯损有重要的影响。3.3 阈值分析
根据相应的阈值标准，本研究枯损株数与总林木株数的比为0.06，则割点阈值选择为A = 0.06。而传统阈值0.5被舍去。图1列出了随机阈值和错误分类率最小的阈值。其中基础模型（M1）的随机阈值为B = 0.08，而错误分类率最低阈值为C =0.07（图1(a)）。混合模型（M3）的随机阈值为B =0.07，而错误分类率最低阈值为C = 0.06（图1(b)）。
0.30.2
0.1
假阴性率 False negative rate (FNR)假阳性率 False positive rate (FPR)错误分类率
Mistake classification rate (MCR)
假阴性率假阳性率错误分类率
Mistake classification rate (MCR)(B =0.08)
(B =0.07)
C =0.07)(C =0.06)
0.70.60.50.40.30.2
1962饥荒0.10阈值 Threshold
(a) (b) 随机效应 Random effect
图 1 蒙古栎混交林枯损模型阈值点与分类率关系
Fig. 1 Relationship of classification rate and threshold for the mortality model of Mongolian oak
2010年央视中秋晚会表4表明，基础模型（M1）中，当阈值为0.06时枯损的准确分类率（TPR ）最高，为81.9%，当阈值为0.08时，存活的准确分类率（TNR ）最高，为72.7%；混合模型（M3）中，当阈值为0.06时枯损的准确分类率（TPR ）最高，为84.1%，当阈值为0.07时，存活的准确分类率（TNR ）最高，为75.4%。在考虑混合效应后，无论阈值为0.06、0.07或0.08，枯损和存活的准确分类率均高于基础模型，而枯损和存活的错误分类率均低于基础模
型。通过计算总正确分类率（ACR ），基础模型阈值为0.08时，总正确分类率最高，为72.6%；混合模型阈值为0.07时，总正确分类率最高，为75.7%。因此可以选择灵敏度与特异度曲线交叉时的阈值（随机阈值）作为最佳阈值。3.4 模拟验证
根据表3的模拟结果，分别利用M1和M3模型形式对验证数据每株树木进行枯损概率的计算，最后计算各自的AUC 值，并比较各模型的灵敏
VBR
第 6 期李春明： 基于广义线性混合效应模型的蒙古栎林单木枯损建模及影响因子分析109
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议