基于Cox 比例風險函數及混合效應的落葉松云冷杉混交林林木枯損模型研究

李春明

(中國林業科學研究院資源信息研究所 北京 100091)

林木的枯損是森林生長過程中十分重要的組成部分，但又是一個知之甚少的森林動態過程[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）分析模型的適應性，在適應的模型基礎上考慮樣地的隨機效應。并對傳統方法與考慮樣地隨機效應的模擬結果進行比較分析。最后按對枯損有影響的協變量的不同情況，預測不同時間內林木的生存率。

1 材料和方法

1.1 研究區概況

本研究選擇了位于吉林省汪清林業局金溝嶺林場境內的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。

1.2 研究方法

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

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

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

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

1.2.2 模型評價方法 在選擇Cox 半參數方法進行林木枯損及生存模擬時，首先要判斷實驗數據是否滿足Cox 比例風險函數模型的假設。擬合優度檢驗方法是十分普遍的方法[27]，其主要思想是首先計算待檢驗變量的Schoenfeld 的殘差，然后對非刪失數據的生存時間排秩，最后利用殘差圖來檢驗Schoenfeld 殘差與時間秩的相關性。如果兩者存在關聯性則認為該數據不滿足Cox 比例風險假設。具體的Schoenfeld 殘差公式見式（2）：

其中， xik為在ti時刻發生枯損的林木的第 k個協變量的實際取值， xmk為 ti時刻尚在風險集里的林木 m的第 k個協變量取值， pm為風險集中林木 m在 ti時刻發生枯損的概率。

在采用Schoenfeld 殘差圖判斷是否可以用Cox 比例風險函數后，本研究采用Wald 指標來評價模型本身的模擬效果。在比較不同模型的擬合效果優劣時，采用AIC 信息準則（Akaike information criterion，AIC）、BIC 信息準則（Bayesian Information Criterion，BIC）和-2*對數似然值（-2*Loglikelihood，-2LogL）這3 個值。這3 個值越小，說明模型的模擬效果越好[28-30]。LRT 卡方檢驗被用來比較模型之間的差異顯著程度[31]。

2 結果和分析

把所有備選的林分因子和立地因子作為協變量加入到Cox 比例風險函數模型中去，利用SAS9.4軟件進行模擬。考慮各因子之間存在的多重共線性（VIF＜5）和顯著影響（ α ＜ 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）。

表2 Cox 比例風險函數模型模擬結果Table 2 The result of Cox proportional hazard model

圖1 各協變量與觀測時間的殘差分布Fig.1 The residual figure between estimate and measure result of covariate variable

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

3 討論

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

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

圖2 10、15 和20 cm 初始胸徑樹木生存率（S (t)）趨勢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)）趨勢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 的樹木，在整個測量區間，隨著時間的推移，初始林分密度大的林分，林木的枯損概率逐漸增大并且生存率較初始密度小的林分快速降低。這充分說明，同一株林木，初值密度越大，對資源的競爭越激烈，與初始林分密度小的林分相比，林木的生存率越低。

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

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

Cox 比例風險函數在考慮了樣地的隨機截距效應后，模型的模擬精度顯著提高，并且差異達顯著程度。但是初始林分密度和大于對象木斷面積的影響差異并不顯著，在去掉這兩個協變量重新模擬后發現模型的模擬效果并沒有降低。因此，初始胸徑是最重要的影響林木本身枯損的指標[33-34]。

4 結論

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