-
林木的枯损是森林生长过程中十分重要的组成部分,但又是一个知之甚少的森林动态过程[1-3]。对林木枯损进行准确的预测是森林生长收获预估系统中十分重要的工作[4]。最早预测林木的枯损是通过建立枯损与协变量的回归方程来预测枯损概率,之后一些学者也采用此方法估计林木的枯损概率[5]。由于林木的枯损是典型的二分量数据,只包括是否枯损两种结局,因此此类方法的模拟效果并不理想。Logistic回归分析方法对于二分量数据具有先天的优势,Hamilton[6]和Monserud[7]引入了此方法来模拟林木的枯损,目前已被广泛应用[8-10]。但是上述方法在模拟林木枯损时,缺乏把树木枯损的时间特性、假设检验、删失数据的观测和协变量的影响同时考虑[11]。即不能够给出树木具体的生存或死亡时间,也不能预测未来林木的发展趋势[12]。林木枯损数据的这些特性使得利用传统的统计方法处理具有困难,但是忽略它们将降低模型估计的准确性[13]。
生存分析方法既考虑观测对象的事件结局“生存”或“死亡”,又能够考虑观测对象的生存时间长短,可为构建新的林木枯损模型提供支撑。是唯一考虑删失数据,并且包括时间协变量,另外能够处理非正态分布数据的一类方法[11]。Waters[14]第一次建议利用生存分析方法来描述林木的枯损问题,但仅限于同龄林分。随着统计学和计算机技术的快速发展,生存分析方法逐渐地被用到分析和构建林木的枯损模型中。例如,Woodall等[11]采用生存分析方法中的生命表法(lifetime table),Fan等[15]利用Kaplan-Meier方法,von Gadow等[16]利用服从Weibull时间分布的参数方法,Uzoh等[17]利用Cox比例风险函数模型来分析上层林木的枯损。这些学者把林木的枯损和生存时间结合起来,取得了很好的模拟和预测结果。
在构建枯损模型时,很多数据都是基于永久性固定样地的多次定期观测,每个林分包括多个样地,每个样地包括多株枯损树木,这些数据存在着多层嵌套结构。另外林木在生长过程中,受竞争及干扰等诸多因素的影响,导致了林木枯损具有很高的时空变化性[18]。在利用生存分析方法构建枯损模型时,忽略这些问题将会造成大的误差。生存分析方法结合混合效应模型方法是一种新的思路,可解决上述诸多问题,并且能够提高模型的模拟精度。目前把两者结合起来构建林木枯损的模型还未见诸于文献。
本研究以1986年在吉林省汪清林业局金沟岭林场设置的20块落叶松云冷杉样地数据为例,利用生存分析方法中的半参数Cox比例风险函数模型来进行林木的枯损及生存模拟,并采用逐步回归方法把对枯损具有重要影响的林分因子和立地因子作为协变量加入到林木枯损模型中去,选择Schoenfeld残差方法(Schoenfeld residuals)分析模型的适应性,在适应的模型基础上考虑样地的随机效应。并对传统方法与考虑样地随机效应的模拟结果进行比较分析。最后按对枯损有影响的协变量的不同情况,预测不同时间内林木的生存率。
HTML
-
本研究选择了位于吉林省汪清林业局金沟岭林场境内的1986年设置的20块落叶松云冷杉混交林为研究对象,优势树种主要包括长白落叶松(Larix olgensis Henry)、云杉(Picea jazoensis Nakai)和冷杉(Abies nephrolepis (Trautv.) Maxim)等,其他伴生树种包括红松(Pinus koraiensis Siebold et Zuccarini)、白桦(Betula platyphylla Suk.)、枫桦(Betula costata Trautv.)和水曲柳(Fraxinus mandshurica Rupr.)等。所有样地从1987年开始,分别在1988、1990、1992、1994、1997、1999、2002、2004、2006、2008、2009和2012年进行了调查。调查的样地因子包括坡度、海拔、坡向、土壤类型、土层厚度等内容。调查因子包括每木检尺(胸径>5 cm)、林木枯损情况、枯损时间、枯损原因、林下更新及林下植被等。其中2块样地在1987年和1993年进行了2次间伐,其余18块只在1987年进行了间伐,因此只选择了18块样地,共3 477株树木进行分析。外业调查结束后,对林分平均直径、公顷株数、公顷断面积、大于对象木断面积(BAL)、公顷蓄积量和直径比等指标进行了计算。其中大于对象木断面积的定义为具体一个样地内,任意一株树木,胸径大于它的所有树木的断面积之和。从1987—2012年整个观测期间,共有1 352株树木发生枯损。1987年伐后调查的样地因子概况见表1。
样地号
Plot number面积
Area/hm2海拔
Elavation/m坡向
Aspect坡度
Slope平均胸径
The mean diameter/cm公顷株数
The stem number/(N·hm−2)公顷断面积
Basal area per hectare/(m2·hm−2)301 0.077 5 760 东北 10 14.93 955 16.83 302 0.077 5 760 东北 10 14.71 1 572 26.69 303 0.13 760 东北 10 14.26 1 098 17.64 304 0.097 5 760 东北 10 15.88 883 17.56 305 0.2 780 西 18 15.07 1 158 20.64 306 0.2 780 东北 7 13.66 2 008 29.42 307 0.2 780 西 18 15.02 1 153 20.50 308 0.2 780 东北 10 15.82 846 16.39 309 0.25 660 西北 6 15.61 1 189 22.86 310 0.25 670 西北 10 16.18 860 17.74 311 0.25 670 西北 6 17.20 833 19.39 312 0.25 680 西北 10 15.49 914 17.30 315 0.112 5 660 北 7 13.58 1 694 24.44 316 0.1 645 北 7 14.87 1 127 19.66 317 0.1 615 东北 7 13.51 1 286 18.65 318 0.112 5 610 东北 7 14.39 919 15.03 319 0.1 605 东北 9 15.17 820 14.84 320 0.1 600 东北 9 13.65 1 670 24.41 Table 1. The characteristics of experimental stands established plot
-
生存分析方法的基本概念可参考文献Allison[19]和Lawless[20],常用概率密度函数(
$f(t)$ )、生存函数($S(t)$ )以及风险函数($h(t)$ )等3个函数来描述生存过程。这3种函数在数学上是等价的。如果给定其中一种函数,另两种函数即可推导得出。在构建林木的枯损模型时,生存时间虽为定量指标,但往往不服从任何规则分布,要在这种情况下确定林木的枯损时间和各风险因素之间的定量关系,就不能够采用参数方法来达到目的。Cox比例风险函数模型可直接建立终点事件的发生风险与这些影响因素之间的函数关系[21]。本研究采用Cox比例风险函数模型和样地水平的随机效应,对应的风险函数可以表示为式(1):${h_0}(t)$ 是仅与时间有关的风险函数,其分布没有明确的假设,是模型的非参数部分,$\exp (\beta {x^T})$ 是在${h_0}(t)$ 效应基础上产生作用的,估计值不受时间影响,是模型的参数部分。${b_i}$ 为样地($i = 1,2, \cdot \cdot \cdot, 18$ )间截距的随机效应值,服从均值为0,方差为${\sigma ^2}$ 的正态分布。在进行生存分析时,影响林木枯损的因子是作为协变量加入到模型中去。影响林木枯损的因子很多,本研究只考虑立地因子和调查初始林分因子对林木枯损的影响。林分因子主要包括树种本身的特性、林分密度、林分结构及林木间的竞争等[22-24]。本研究主要选择的林分因子包括单木初始胸径、单木初始胸径与林分平均胸径的比值、大于对象木断面积、初始林分公顷株数、初始林分公顷断面积、初始林分公顷蓄积、初始林分平均胸径等指标因子。立地因子影响森林环境中光照、温度、水分、土壤等的分布,进而影响林木的枯损状况[25]。本研究主要选择坡度、坡向及海拔等因子对林木枯损的影响[26]。
-
在选择Cox半参数方法进行林木枯损及生存模拟时,首先要判断实验数据是否满足Cox比例风险函数模型的假设。拟合优度检验方法是十分普遍的方法[27],其主要思想是首先计算待检验变量的Schoenfeld的残差,然后对非删失数据的生存时间排秩,最后利用残差图来检验Schoenfeld残差与时间秩的相关性。如果两者存在关联性则认为该数据不满足Cox比例风险假设。具体的Schoenfeld残差公式见式(2):
其中,
${x_{ik}}$ 为在${t_i}$ 时刻发生枯损的林木的第$k$ 个协变量的实际取值,${x_{mk}}$ 为${t_i}$ 时刻尚在风险集里的林木$m$ 的第$k$ 个协变量取值,${p_m}$ 为风险集中林木$m$ 在${t_i}$ 时刻发生枯损的概率。在采用Schoenfeld残差图判断是否可以用Cox比例风险函数后,本研究采用Wald指标来评价模型本身的模拟效果。在比较不同模型的拟合效果优劣时,采用AIC信息准则(Akaike information criterion,AIC)、BIC信息准则(Bayesian Information Criterion,BIC)和-2*对数似然值(-2*Loglikelihood,-2LogL)这3个值。这3个值越小,说明模型的模拟效果越好[28-30]。LRT卡方检验被用来比较模型之间的差异显著程度[31]。
1.1. 研究区概况
1.2. 研究方法
1.2.1. 生存分析方法简介
1.2.2. 模型评价方法
-
把所有备选的林分因子和立地因子作为协变量加入到Cox比例风险函数模型中去,利用SAS9.4软件进行模拟。考虑各因子之间存在的多重共线性(VIF<5)和显著影响(
$\alpha < 0.05$ )情况下,结果最后只有单木初始胸径、大于对象木断面积和初始林分公顷株数等林分因子对林木的枯损具有重要的影响。而立地因子的3个具体因子对林木的枯损均没有显著影响,具体模拟结果见表2的模型1(M1)。利用Schoenfeld残差分布图来判断实验数据及协变量是否满足Cox比例风险函数模型的条件假设,结果见图1。从图1的3个残差图可以看出,单木初始胸径与观测时间不呈线性关系、大于对象木断面积、初始林分公顷株数与观测时间不呈线性关系,说明本次研究采用的数据完全符合Cox比例风险函数的条件假设,可以选择Cox比例风险函数来构建林木的枯损模型。Wald指标表明,在考虑了3个协变量后,3个模型的模拟效果差异均极度显著(p<0.000 1),能够满足模型的精度要求。进一步在Cox模型的基础上,考虑样地的随机效应,本研究只考虑了样地的截距效应,结果为模型的AIC、BIC和-2LogL3个值均比不考虑随机效应值小,LRT卡方检验表明,差异达显著水平(见表2的模型2)。由于考虑了样地的随机效应后,大于对象木断面积和初始林分公顷株数两个协变量差异不显著(p>0.05),因此去掉了这两个协变量,然后重新模拟,结果见表2的模型3(M3)。最后结果表明,当考虑样地水平的随机效应后,去掉大于对象木断面积和初始林分公顷株数两个协变量后,模型精度并没有降低,两者之间没有显著差异(LRT=4.2,$p > 0.05$ )。协变量
Covariatevariable模型1(M1)
Model 1模型2(M2)
Model 2模型3(M3)
Model 3估计值(标准差)
Estimate result(Standard error)P值
P value估计值(标准差)
Estimate result(Standard error)P值
P value估计值(标准差)
Estimate result(Standard error)P值
P value单木初始胸径
Initial DBH/d−0.052 7(0.007 6) <0.000 1 −0.077 8(0.012 4) <0.000 1 −0.089 5(0.006 5) <0.000 1 大于对象木断面积(BAL)
Basal area of the trees larger
than the subject tree/(m2·hm−2)0.104 7(0.021 1) <0.000 1 0.056 4(0.048 5) 0.244 7 初始林分公顷株数
Initial stem number/(N·hm−2)0.656 5(0.066 4) <0.000 1 0.535 8(0.311 2) 0.085 2 Wald指标
Wald value304.078 9 <0.000 1 528.049 0 <0.000 1 523.407 5 <0.000 1 AIC 21 164.6 20 944.6 20 948.8 BIC 21 180.3 20 960.3 20 964.5 −2LogL 21 158.6 20 938.6 20 942.8 LRT M2/M1LRT=220($p < 0.001$) M3/M2LRT=4.2($p > 0.05$) Table 2. The result of Cox proportional hazard model
从表2的结果可以看出,单木的初始胸径与林木枯损的风险率呈反比,与生存率呈正比,初始胸径大的林木,枯损的风险低于胸径小的林木,进而生存率高于胸径小的林木。这与利用logistic方法模拟枯损的结论基本一致[32]。