基于Cox 比例風(fēng)險(xiǎn)函數(shù)及混合效應(yīng)的落葉松云冷杉混交林林木枯損模型研究

李春明

(中國(guó)林業(yè)科學(xué)研究院資源信息研究所 北京 100091)

林木的枯損是森林生長(zhǎng)過程中十分重要的組成部分，但又是一個(gè)知之甚少的森林動(dòng)態(tài)過程[1-3]。對(duì)林木枯損進(jìn)行準(zhǔn)確的預(yù)測(cè)是森林生長(zhǎng)收獲預(yù)估系統(tǒng)中十分重要的工作[4]。最早預(yù)測(cè)林木的枯損是通過建立枯損與協(xié)變量的回歸方程來預(yù)測(cè)枯損概率，之后一些學(xué)者也采用此方法估計(jì)林木的枯損概率[5]。由于林木的枯損是典型的二分量數(shù)據(jù)，只包括是否枯損兩種結(jié)局，因此此類方法的模擬效果并不理想。Logistic 回歸分析方法對(duì)于二分量數(shù)據(jù)具有先天的優(yōu)勢(shì)，Hamilton[6]和Monserud[7]引入了此方法來模擬林木的枯損，目前已被廣泛應(yīng)用[8-10]。但是上述方法在模擬林木枯損時(shí)，缺乏把樹木枯損的時(shí)間特性、假設(shè)檢驗(yàn)、刪失數(shù)據(jù)的觀測(cè)和協(xié)變量的影響同時(shí)考慮[11]。即不能夠給出樹木具體的生存或死亡時(shí)間，也不能預(yù)測(cè)未來林木的發(fā)展趨勢(shì)[12]。林木枯損數(shù)據(jù)的這些特性使得利用傳統(tǒng)的統(tǒng)計(jì)方法處理具有困難，但是忽略它們將降低模型估計(jì)的準(zhǔn)確性[13]。

生存分析方法既考慮觀測(cè)對(duì)象的事件結(jié)局“生存”或“死亡”，又能夠考慮觀測(cè)對(duì)象的生存時(shí)間長(zhǎng)短，可為構(gòu)建新的林木枯損模型提供支撐。是唯一考慮刪失數(shù)據(jù)，并且包括時(shí)間協(xié)變量，另外能夠處理非正態(tài)分布數(shù)據(jù)的一類方法[11]。Waters[14]第一次建議利用生存分析方法來描述林木的枯損問題，但僅限于同齡林分。隨著統(tǒng)計(jì)學(xué)和計(jì)算機(jī)技術(shù)的快速發(fā)展，生存分析方法逐漸地被用到分析和構(gòu)建林木的枯損模型中。例如，Woodall 等[11]采用生存分析方法中的生命表法（lifetime table），F(xiàn)an 等[15]利用Kaplan-Meier 方法，von Gadow 等[16]利用服從Weibull 時(shí)間分布的參數(shù)方法，Uzoh 等[17]利用Cox 比例風(fēng)險(xiǎn)函數(shù)模型來分析上層林木的枯損。這些學(xué)者把林木的枯損和生存時(shí)間結(jié)合起來，取得了很好的模擬和預(yù)測(cè)結(jié)果。

在構(gòu)建枯損模型時(shí)，很多數(shù)據(jù)都是基于永久性固定樣地的多次定期觀測(cè)，每個(gè)林分包括多個(gè)樣地，每個(gè)樣地包括多株枯損樹木，這些數(shù)據(jù)存在著多層嵌套結(jié)構(gòu)。另外林木在生長(zhǎng)過程中，受競(jìng)爭(zhēng)及干擾等諸多因素的影響，導(dǎo)致了林木枯損具有很高的時(shí)空變化性[18]。在利用生存分析方法構(gòu)建枯損模型時(shí)，忽略這些問題將會(huì)造成大的誤差。生存分析方法結(jié)合混合效應(yīng)模型方法是一種新的思路，可解決上述諸多問題，并且能夠提高模型的模擬精度。目前把兩者結(jié)合起來構(gòu)建林木枯損的模型還未見諸于文獻(xiàn)。

本研究以1986 年在吉林省汪清林業(yè)局金溝嶺林場(chǎng)設(shè)置的20 塊落葉松云冷杉樣地?cái)?shù)據(jù)為例，利用生存分析方法中的半?yún)?shù)Cox 比例風(fēng)險(xiǎn)函數(shù)模型來進(jìn)行林木的枯損及生存模擬，并采用逐步回歸方法把對(duì)枯損具有重要影響的林分因子和立地因子作為協(xié)變量加入到林木枯損模型中去，選擇Schoenfeld 殘差方法（Schoenfeld residuals）分析模型的適應(yīng)性，在適應(yīng)的模型基礎(chǔ)上考慮樣地的隨機(jī)效應(yīng)。并對(duì)傳統(tǒng)方法與考慮樣地隨機(jī)效應(yīng)的模擬結(jié)果進(jìn)行比較分析。最后按對(duì)枯損有影響的協(xié)變量的不同情況，預(yù)測(cè)不同時(shí)間內(nèi)林木的生存率。

1 材料和方法

1.1 研究區(qū)概況

本研究選擇了位于吉林省汪清林業(yè)局金溝嶺林場(chǎng)境內(nèi)的1986 年設(shè)置的20 塊落葉松云冷杉混交林為研究對(duì)象，優(yōu)勢(shì)樹種主要包括長(zhǎng)白落葉松（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 年進(jìn)行了調(diào)查。調(diào)查的樣地因子包括坡度、海拔、坡向、土壤類型、土層厚度等內(nèi)容。調(diào)查因子包括每木檢尺（胸徑＞5 cm）、林木枯損情況、枯損時(shí)間、枯損原因、林下更新及林下植被等。其中2 塊樣地在1987 年和1993 年進(jìn)行了2 次間伐，其余18 塊只在1987 年進(jìn)行了間伐，因此只選擇了18 塊樣地，共3 477 株樹木進(jìn)行分析。外業(yè)調(diào)查結(jié)束后，對(duì)林分平均直徑、公頃株數(shù)、公頃斷面積、大于對(duì)象木斷面積（BAL）、公頃蓄積量和直徑比等指標(biāo)進(jìn)行了計(jì)算。其中大于對(duì)象木斷面積的定義為具體一個(gè)樣地內(nèi)，任意一株樹木，胸徑大于它的所有樹木的斷面積之和。從1987—2012 年整個(gè)觀測(cè)期間，共有1 352 株樹木發(fā)生枯損。1987 年伐后調(diào)查的樣地因子概況見表1。

1.2 研究方法

1.2.1 生存分析方法簡(jiǎn)介 生存分析方法的基本概念可參考文獻(xiàn)Allison[19]和Lawless[20]，常用概率密度函數(shù)（ f(t)） 、生存函數(shù)（ S(t)）以及風(fēng)險(xiǎn)函數(shù)（ h(t)）等3 個(gè)函數(shù)來描述生存過程。這3 種函數(shù)在數(shù)學(xué)上是等價(jià)的。如果給定其中一種函數(shù)，另兩種函數(shù)即可推導(dǎo)得出。在構(gòu)建林木的枯損模型時(shí)，生存時(shí)間雖為定量指標(biāo)，但往往不服從任何規(guī)則分布，要在這種情況下確定林木的枯損時(shí)間和各風(fēng)險(xiǎn)因素之間的定量關(guān)系，就不能夠采用參數(shù)方法來達(dá)到目的。Cox 比例風(fēng)險(xiǎn)函數(shù)模型可直接建立終點(diǎn)事件的發(fā)生風(fēng)險(xiǎn)與這些影響因素之間的函數(shù)關(guān)系[21]。本研究采用Cox 比例風(fēng)險(xiǎn)函數(shù)模型和樣地水平的隨機(jī)效應(yīng)，對(duì)應(yīng)的風(fēng)險(xiǎn)函數(shù)可以表示為式（1）：

h0(t)是僅與時(shí)間有關(guān)的風(fēng)險(xiǎn)函數(shù)，其分布沒有明確的假設(shè)，是模型的非參數(shù)部分， exp (βxT)是在h0(t)效應(yīng)基礎(chǔ)上產(chǎn)生作用的，估計(jì)值不受時(shí)間影響，是模型的參數(shù)部分。 bi為 樣地（ i= 1,2,···,18）間截距的隨機(jī)效應(yīng)值，服從均值為0，方差為 σ2的正態(tài)分布。

在進(jìn)行生存分析時(shí)，影響林木枯損的因子是作為協(xié)變量加入到模型中去。影響林木枯損的因子很多，本研究只考慮立地因子和調(diào)查初始林分因子對(duì)林木枯損的影響。林分因子主要包括樹種本身的特性、林分密度、林分結(jié)構(gòu)及林木間的競(jìng)爭(zhēng)等[22-24]。本研究主要選擇的林分因子包括單木初始胸徑、單木初始胸徑與林分平均胸徑的比值、大于對(duì)象木斷面積、初始林分公頃株數(shù)、初始林分公頃斷面積、初始林分公頃蓄積、初始林分平均胸徑等指標(biāo)因子。立地因子影響森林環(huán)境中光照、溫度、水分、土壤等的分布，進(jìn)而影響林木的枯損狀況[25]。本研究主要選擇坡度、坡向及海拔等因子對(duì)林木枯損的影響[26]。

表1 樣地1987 年基本情況Table 1 The characteristics of experimental stands established plot

1.2.2 模型評(píng)價(jià)方法 在選擇Cox 半?yún)?shù)方法進(jìn)行林木枯損及生存模擬時(shí)，首先要判斷實(shí)驗(yàn)數(shù)據(jù)是否滿足Cox 比例風(fēng)險(xiǎn)函數(shù)模型的假設(shè)。擬合優(yōu)度檢驗(yàn)方法是十分普遍的方法[27]，其主要思想是首先計(jì)算待檢驗(yàn)變量的Schoenfeld 的殘差，然后對(duì)非刪失數(shù)據(jù)的生存時(shí)間排秩，最后利用殘差圖來檢驗(yàn)Schoenfeld 殘差與時(shí)間秩的相關(guān)性。如果兩者存在關(guān)聯(lián)性則認(rèn)為該數(shù)據(jù)不滿足Cox 比例風(fēng)險(xiǎn)假設(shè)。具體的Schoenfeld 殘差公式見式（2）：

其中， xik為在ti時(shí)刻發(fā)生枯損的林木的第 k個(gè)協(xié)變量的實(shí)際取值， xmk為 ti時(shí)刻尚在風(fēng)險(xiǎn)集里的林木 m的第 k個(gè)協(xié)變量取值， pm為風(fēng)險(xiǎn)集中林木 m在 ti時(shí)刻發(fā)生枯損的概率。

在采用Schoenfeld 殘差圖判斷是否可以用Cox 比例風(fēng)險(xiǎn)函數(shù)后，本研究采用Wald 指標(biāo)來評(píng)價(jià)模型本身的模擬效果。在比較不同模型的擬合效果優(yōu)劣時(shí)，采用AIC 信息準(zhǔn)則（Akaike information criterion，AIC）、BIC 信息準(zhǔn)則（Bayesian Information Criterion，BIC）和-2*對(duì)數(shù)似然值（-2*Loglikelihood，-2LogL）這3 個(gè)值。這3 個(gè)值越小，說明模型的模擬效果越好[28-30]。LRT 卡方檢驗(yàn)被用來比較模型之間的差異顯著程度[31]。

2 結(jié)果和分析

把所有備選的林分因子和立地因子作為協(xié)變量加入到Cox 比例風(fēng)險(xiǎn)函數(shù)模型中去，利用SAS9.4軟件進(jìn)行模擬?？紤]各因子之間存在的多重共線性（VIF＜5）和顯著影響（ α ＜ 0.05）情況下，結(jié)果最后只有單木初始胸徑、大于對(duì)象木斷面積和初始林分公頃株數(shù)等林分因子對(duì)林木的枯損具有重要的影響。而立地因子的3 個(gè)具體因子對(duì)林木的枯損均沒有顯著影響，具體模擬結(jié)果見表2 的模型1（M1）。利用Schoenfeld 殘差分布圖來判斷實(shí)驗(yàn)數(shù)據(jù)及協(xié)變量是否滿足Cox 比例風(fēng)險(xiǎn)函數(shù)模型的條件假設(shè)，結(jié)果見圖1。從圖1 的3 個(gè)殘差圖可以看出，單木初始胸徑與觀測(cè)時(shí)間不呈線性關(guān)系、大于對(duì)象木斷面積、初始林分公頃株數(shù)與觀測(cè)時(shí)間不呈線性關(guān)系，說明本次研究采用的數(shù)據(jù)完全符合Cox 比例風(fēng)險(xiǎn)函數(shù)的條件假設(shè)，可以選擇Cox 比例風(fēng)險(xiǎn)函數(shù)來構(gòu)建林木的枯損模型。Wald 指標(biāo)表明，在考慮了3 個(gè)協(xié)變量后，3 個(gè)模型的模擬效果差異均極度顯著（p＜0.000 1），能夠滿足模型的精度要求。進(jìn)一步在Cox 模型的基礎(chǔ)上，考慮樣地的隨機(jī)效應(yīng)，本研究只考慮了樣地的截距效應(yīng)，結(jié)果為模型的AIC、BIC 和-2LogL3 個(gè)值均比不考慮隨機(jī)效應(yīng)值小，LRT 卡方檢驗(yàn)表明，差異達(dá)顯著水平（見表2 的模型2）。由于考慮了樣地的隨機(jī)效應(yīng)后，大于對(duì)象木斷面積和初始林分公頃株數(shù)兩個(gè)協(xié)變量差異不顯著（p＞0.05），因此去掉了這兩個(gè)協(xié)變量，然后重新模擬，結(jié)果見表2 的模型3（M3）。最后結(jié)果表明，當(dāng)考慮樣地水平的隨機(jī)效應(yīng)后，去掉大于對(duì)象木斷面積和初始林分公頃株數(shù)兩個(gè)協(xié)變量后，模型精度并沒有降低，兩者之間沒有顯著差異（LRT=4.2， p ＞ 0.05）。

表2 Cox 比例風(fēng)險(xiǎn)函數(shù)模型模擬結(jié)果Table 2 The result of Cox proportional hazard model

圖1 各協(xié)變量與觀測(cè)時(shí)間的殘差分布Fig.1 The residual figure between estimate and measure result of covariate variable

從表2 的結(jié)果可以看出，單木的初始胸徑與林木枯損的風(fēng)險(xiǎn)率呈反比，與生存率呈正比，初始胸徑大的林木，枯損的風(fēng)險(xiǎn)低于胸徑小的林木，進(jìn)而生存率高于胸徑小的林木。這與利用logistic 方法模擬枯損的結(jié)論基本一致[32]。

3 討論

在實(shí)際林分中，胸徑小的樹木在林分中處于競(jìng)爭(zhēng)的劣勢(shì)，對(duì)水、肥、氣、熱等資源的競(jìng)爭(zhēng)處于不利地位。而胸徑大的樹木處于競(jìng)爭(zhēng)的優(yōu)勢(shì)，更加有利于爭(zhēng)奪林分中的營(yíng)養(yǎng)。因此，胸徑小的樹木更容易發(fā)生枯損。為了直觀的表述初始胸徑與林木枯損的關(guān)系，進(jìn)一步選擇初始胸徑為10、15 和20 cm的3 株樹木，在初始林分密度均為2 000 株·hm-2的林分內(nèi)，在不考慮其他因子影響的前提下，利用Cox 比例風(fēng)險(xiǎn)函數(shù)模型來估計(jì)其生存率，并繪制生存曲線，見圖2。從圖中可以看出，隨著時(shí)間的推移，胸徑大的樹木其生存率始終高于胸徑小的樹木，這和上述模擬的結(jié)果一致。在同一林分內(nèi)，一般初始胸徑小的林木其大于對(duì)象木斷面積值就大，因此與初始胸徑相反，大于對(duì)象木斷面積與林木枯損的風(fēng)險(xiǎn)率呈正比，與生存率呈反比。

林分公頃株數(shù)是一個(gè)重要的影響因子，屬于林分層次，反映了林分的密度情況。本研究中初始林分公頃株數(shù)與林木枯損的風(fēng)險(xiǎn)率呈正比，與生存率呈反比。這說明，林分密度越大，林木枯損的風(fēng)險(xiǎn)越大，生存率越低，越容易枯損[8,26]。為了更直觀的表示林分公頃株數(shù)對(duì)枯損的影響，進(jìn)一步選擇胸徑為15 cm 的樹木在初始林分密度為1 000，1 500和2 000 株·hm-2的3 個(gè)林分情況下，在不考慮其他因子影響的前提下，利用Cox 比例風(fēng)險(xiǎn)函數(shù)模型來計(jì)算其生存率，并繪制生存曲線，見圖3。

圖2 10、15 和20 cm 初始胸徑樹木生存率（S (t)）趨勢(shì)Fig.2 Thesurvival curve of tree of initial diameter with 10, 15 and 20 cm

圖3 初始林分密度為1 000,1 500 和2 000 株·hm-2 的林木生存率（S (t)）趨勢(shì)Fig.3 The survival curve of tree with initial stand density of 1 000, 1 500 and 2 000 trees·hm-2

從圖3 可以看出，胸徑為15 cm 的樹木，在整個(gè)測(cè)量區(qū)間，隨著時(shí)間的推移，初始林分密度大的林分，林木的枯損概率逐漸增大并且生存率較初始密度小的林分快速降低。這充分說明，同一株林木，初值密度越大，對(duì)資源的競(jìng)爭(zhēng)越激烈，與初始林分密度小的林分相比，林木的生存率越低。

立地因子對(duì)林木的枯損沒有任何影響，造成這一原因可能是選擇的20 塊落葉松云冷杉樣地分布在一個(gè)林場(chǎng)內(nèi)，坡度、坡向和海拔本身的差別小，難以體現(xiàn)在對(duì)枯損的影響上。

從圖2 和圖3 可以看出，林木在觀測(cè)15 a 至16 a 之間生存率有一個(gè)明顯的下降。從數(shù)據(jù)本身來看，從1986—2012 年整個(gè)觀測(cè)期內(nèi)共枯損1 352株樹木，而在1999—2002 年觀測(cè)期間就枯損了405 株，差不多占整個(gè)枯損量的三分之一。造成這一期間大量枯損的原因可能是極端氣候引起的（干旱，霜凍），由于本研究沒有考慮氣象因子對(duì)枯損的影響，因此無法獲得確切的原因，今后需要進(jìn)一步考慮。

Cox 比例風(fēng)險(xiǎn)函數(shù)在考慮了樣地的隨機(jī)截距效應(yīng)后，模型的模擬精度顯著提高，并且差異達(dá)顯著程度。但是初始林分密度和大于對(duì)象木斷面積的影響差異并不顯著，在去掉這兩個(gè)協(xié)變量重新模擬后發(fā)現(xiàn)模型的模擬效果并沒有降低。因此，初始胸徑是最重要的影響林木本身枯損的指標(biāo)[33-34]。

4 結(jié)論

本研究利用Cox 比例風(fēng)險(xiǎn)函數(shù)描述落葉松云冷杉林木的枯損情況，并考慮了樣地的隨機(jī)效應(yīng)，取得了很好的效果。林木初始胸徑、大于對(duì)象木斷面積及初始林分公頃株數(shù)對(duì)林木的枯損具有重要的影響。在考慮樣地的隨機(jī)效應(yīng)后，模型的模擬效果顯著提高，說明樣地的枯損存在著明顯的高可變性和隨機(jī)性。在實(shí)際中，林木初始胸徑是不容易控制的指標(biāo)，但是初始林分密度可通過科學(xué)的采伐來調(diào)節(jié)。如果想降低林木的枯損率，可適當(dāng)降低林分的密度，進(jìn)而能夠促進(jìn)林木本身胸徑的生長(zhǎng)。本研究基于生存分析方法為林木的枯損研究提供了可借鑒的新思路。Cox 比例風(fēng)險(xiǎn)函數(shù)模型的使用，使得森林經(jīng)營(yíng)者在制定森林采伐方案及采伐時(shí)間上提供了很好的科學(xué)依據(jù)。