基于廣義線性混合效應模型的蒙古櫟林單木枯損建模及影響因子分析

李春明

(1. 中國林業科學研究院資源信息研究所 北京 100091；2. 國家林業和草原局森林經營與生長模擬重點實驗室 北京 100091)

作為生態過程中不可或缺的組成部分，林木枯損被認為是森林生態系統中極其關鍵，卻又知之甚少的森林動態過程[1]。單木枯損模型能夠對每株樹木的生長狀況進行模擬，諸多學者利用多種數據和方法構建了大量的單木枯損模型。傳統的經驗方程，如指數方程和Weibull方程等，主要是通過建立枯損與協變量的回歸方程來預測未來樹木的枯損概率[2-3]。由于每株樹木只包括枯損和存活兩種結局，是典型的二分類數據，上述傳統的經驗方程在估計枯損時效果并不理想。而Logistic回歸方法對于二分類數據具有先天的優勢，與傳統的經驗方程相比，模擬林木的枯損有較好的適用性，已經被廣泛應用[4-5]。

在構建單木枯損模型時，以往基于Logistic回歸方法的枯損模型只反映枯損林木總體的平均變化規律，忽略了樣地之間及林木個體之間的差異，也沒有考慮數據間的時空序列相關性，勢必會對模型的準確性造成影響[6]。混合效應模型由固定效應和隨機效應兩部分組成，既可以反映總體的平均變化趨勢，又可以提供數據方差、協方差等多種信息來反映個體之間的差異。另外，在處理連續觀測數據存在的序列相關性等問題時具有其他模型無法比擬的優勢[7]。因此一些學者開始嘗試利用廣義線性混合效應模型來構建單木枯損模型[7-9]。研究結果表明：與傳統的Logistic回歸方法相比，基于廣義線性混合效應模型方法的枯損模型明顯提高了模擬和預測精度。

影響單木枯損的因素很多，主要包括單木因子、林分因子、立地因子、氣象因子以及這些因子的交互作用。單木因子主要包括林木本身的生長和遺傳特性。林分因子方面，主要包括林分密度、競爭指數、林分生產力及林分結構等，多數情況下被同時考慮或考慮其中幾個因子[10-13]。立地因子，主要考慮坡度、坡向、海拔及其數學變化形式對枯損的影響[8]。氣象因子方面，主要考慮溫度、降水量、霜凍及干旱等對枯損的影響[1,14-16]。在實際情況中，樹木的枯損可能會由于多種原因同時發生，例如林分因子、立地因子及氣象因子的相互作用等[17]。

蒙古櫟（Quercus mongolicaFisch. ex Ledeb）是我國東北林區中主要的天然樹種[18]。研究蒙古櫟林的枯損和存活狀態，構建蒙古櫟林的單木枯損模型對于構建林分生長和收獲模型具有重要意義。馬武等[19]以蒙古櫟天然林為研究對象，基于吉林省汪清林業局195塊固定樣地數據，采用二分類的Logistic回歸方法建立了蒙古櫟天然林的單木枯損模型，但沒有考慮枯損數據存在的樣地和個體之間的差異和多次測量的時間序列相關性。本研究以吉林省東南部地區的295塊蒙古櫟林固定樣地多次觀測數據為研究對象，把數據分成兩部分，一部分為模擬數據（約80%），共計236塊樣地38 096條單木記錄，一部分為驗證數據（約20%），共計59塊樣地9 674條單木記錄。基于Logistic回歸及廣義線性混合效應模型方法，考慮單木和林分因子，立地因子和氣象因子的影響，構建單木枯損模型，并與不考慮隨機效應的傳統方法進行比較，然后采用驗證數據進行精度驗證。

1 數據來源

本研究選擇吉林省295塊蒙古櫟樣地數據，樣地大小為0.06 hm2。樣地中蒙古櫟為優勢樹種，其他主要樹種還包括云杉(Picea jazoensisNakai)、白樺(Betula platyphyllaSuk.)、紅松(Pinus koraiensisSieb. et Zucc.)、水曲柳(Fraxinus mandshuricaRupr.)和冷杉(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*對數似然值（-2LogL）這3個值經常被用來比較不同模型之間的擬合效果，這3個值越小，說明模型的模擬效果越好[22]。似然比卡方檢驗（LRT）被用來比較模型之間的差異顯著程度[23]。

AUC值是ROC曲線下面的面積，能夠描述Logistic回歸模型的辨別能力，被廣泛用來評價單木枯損模型的特異性和敏感性。AUC值越大，說明模型的特異性和敏感性越好，具體計算方法見Fielding and Bell[24]。

表 1 樣地各調查因子統計 Table 1 The statistics of plot factors

表 2 主要氣象變量統計 Table 2 The statistical table of main climate variables

2.3 閾值確定

在利用枯損模型來判斷林木是否枯損前需要科學合理的確定閾值。確定閾值的主要方法包括傳統閾值、隨機閾值、錯誤分類率（MCR）最低閾值和割點閾值等[9,25]。傳統閾值通常選擇0.5，該值僅適用于枯損木與存活木相近的情況，實際上枯損是小概率事件；隨機閾值選擇靈敏度與特異度曲線交叉時的閾值為最佳閾值；割點閾值選取林分各個樣地枯損率的平均值作為最佳閾值[26]。本研究選擇了隨機閾值、錯誤分類率最低閾值和割點閾值等3種閾值確定方法，來確定最佳閾值。

3 結果與分析

通過篩選，最后單木初始胸徑、大于對象木斷面積、林分公頃斷面積、林分公頃株數和林分平均年齡，氣象因子中的年平均溫度（MAT）、年平均夏季(5—9月)降水量（MSP）、上一年8月至當年7月的降雪量（PAS）等因子，被保留在基礎模型中。

在Logistic模型的基礎上，考慮樣地水平、單木水平及兩水平（樣地內嵌套單木水平）的隨機效應（分隨機截距效應和隨機系數效應兩種情況）。在考慮單木水平時只有在截距上收斂，但是其AIC、BIC和-2LogL值分別是18 559、18 639和18 539，遠遠高于固定效應模型。在考慮兩水平效應時，無論什么情況，模型均不能夠收斂。在考慮樣地水平隨機系數效應時，3個參數以上的隨機效應均不收斂，考慮2個參數時，只有兩種情況收斂，第一種是同時在截距和胸徑兩個參數上考慮隨機效應，第二種是同時在截距和大于對象木斷面積時，模型能夠收斂。具體模擬和計算的AIC、BIC、-2LogL、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、-2LogL值均比沒考慮樣地的隨機效應值要小，說明在構建單木枯損模型時考慮隨機效應的模擬效果要優于不考慮隨機效應的模型。利用LRT對不同的模型進行方差分析，考慮樣地隨機效應后，與固定效應模型相比，差異均達到極顯著水平（p＜0.001）；當同時考慮截距和隨機系數效應后比單獨考慮隨機截距效應，模型的模擬效果要好，差異也達到極顯著（p＜0.001）。在胸徑和截距上考慮隨機效應的AIC、BIC、-2LogL等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)）。

圖 1 蒙古櫟混交林枯損模型閾值點與分類率關系Fig. 1 Relationship of classification rate and threshold for the mortality model of Mongolian oak

表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值，并比較各模型的靈敏度。其中，M3的模型形式，由于考慮了樣地的隨機效應，則首先要求出驗證數據中樣地的隨機效應參數值，然后再計算枯損概率。隨機效應參數計算方法可參考Vonesh and Chinchilli[27]。

表 4 蒙古櫟混交林單木枯損模型閾值預測列聯表分析結果 Table 4 Confusion matrix result of tree mortality model of Mongolian oak

最后計算出M1模型的AUC值為0.763，M3模型的AUC值為0.766。計算結果表明，在截距和胸徑上同時考慮隨機效應的M3模型的AUC值高于M1模型的AUC值，說明M3模型的靈敏性和特異性要好于M1模型，這與模擬數據中的結論一致。

4 討論

（1）在本研究中，單木枯損模型不僅存在著樣地效應，而且還存在著單木水平上的效應。在模型中已經包含了單木的初始胸徑，再考慮單木水平的隨機效應，通常情況下（經驗）精度會有所提高，但在本研究中精度并沒有提高，相反還遠遠降低。造成的主要原因是在模型中已經包含了單木的初始胸徑，另外涉及到模型的復雜性問題，不一定能提高精度。因此再考慮單木水平的隨機效應，在實際中對于預測并不現實也沒有意義。另外，在預測單木枯損概率時，如果有先驗信息被利用，則就可以估計隨機效應參數值，這樣就能夠提高模型的估計和預測精度[28]。

（2）一般來講在一個林分內，胸徑大的樹木在林分中處于競爭有利的地位，而胸徑小的樹木處于不利地位，就容易枯損[29]。但是在本研究中，M1和M2模型中單木初始胸徑與枯損概率呈正比，M3和M4模型中單木初始胸徑與枯損概率呈反比。除了數據本身和統計方面的原因外，也說明單純一個胸徑因子不足以表述單個林木在林分中的地位和狀態，例如一個林分內都是胸徑小的林木，或者是都是胸徑大的林木，則造成個別林木枯損的原因就不是林木本身的大小，而是其他原因，例如病蟲害會造成很多大樹死亡。大于對象木斷面積反映林木的競爭狀態，在同一林分內胸徑小的樹木，大于對象木斷面積值就大，處于競爭的劣勢。林分公頃斷面積把單木大小和密度結合起來，是一個很好的表達林分擁擠程度的指標，能夠反映系統性競爭[30]。與林齡類似，林分公頃斷面積與單木枯損呈反比，幼、中齡林株數多但斷面積小競爭激烈，近熟林以上就已經過了競爭激烈的階段，林木枯損逐漸變少。林分公頃株數是一個重要的影響因子，在給定年齡和地位指數條件下，林分密度越大，樹木越容易枯損[4]，這給森林經營者在制定森林采伐方案及時間上提供了很好的科學依據。本研究通過逐步回歸方法，林齡保留在模型中，意味著隨著林齡的增加枯損率降低，間接的說明枯損極易發生在密度大的幼齡林里，隨著林分的生長和穩定，枯損率逐漸降低。

（3）隨著全球氣候的變化及極端天氣的經常出現，氣象因子與枯損也具有直接的關系。與林分因子相比，本研究中的3個氣象因子對枯損的影響要略微小一些，這可能是由于林木本身的特性及林分密度折中了一部分氣象因子的影響。與很多學者的研究結論一致，年平均溫度與枯損率呈正比[31-32]。年平均夏季(5—9月)降水量（MSP）與枯損率呈反比，說明隨著溫度的升高，夏季降水量的降低，容易引起干旱，進而會增加樹木枯損的概率。一般情況下，冬季降雪量是一個降水量指標，其值的增加會降低枯損的發生，但本研究中恰好相反，PAS與枯損呈正比，可能的原因是吉林東南部隨著冬季降雪量的增加會增加局部凍害發生的概率，增加了樹木枯損的概率。

（4）本研究雖然考慮了立地因子的影響，但是在逐步回歸過程中，都因為不顯著而被移除出模型。不顯著的原因可能是數據主要來源于吉林省東南部地區，局部的海拔和坡度差別小。本研究由于沒有優勢木平均高數據，無法計算地位指數而被舍棄。但是并不能說明地位指數對枯損沒有影響。閾值分析是研究枯損時十分重要的工作，能夠有效提高模型的判別精度。最佳閾值點并不一定是枯損正確率最高的閾值點，不僅要考慮枯損判斷正確率，還要考慮存活樹木的判斷正確率，否則會有很多實際存活的樹木被錯誤地判斷為死亡樹木，反之亦然。本研究在確定閾值時，確定靈敏度與特異度的交點為最佳閾值點，這與Chen等[33]的結論一致。本研究中平均枯損率為0.06，與最佳閾值點明顯接近，在實際中為了方便也可直接作為最佳閾值點。

5 結論

研究表明，單木初始胸徑對枯損的影響比較復雜，要結合具體林齡等因素綜合考慮。大于對象木斷面積值越小，林分公頃株數越小，林分公頃斷面積和林齡越大，單木枯損率越低。年平均溫度（MAT）越高，年平均夏季(5—9月)降水量（MSP）越低，上一年8月至當年7月的降雪量（PAS）越大，單木枯損率越高。當模型考慮了樣地的隨機效應后，AIC、BIC、-2LogL值均比沒考慮樣地的隨機效應值要小；當同時考慮截距和隨機系數效應后比單獨考慮隨機截距效應模型的模擬效果要好，差異達到極顯著（ p＜0.001）。在胸徑和截距上考慮隨機效應的AIC、BIC、-2LogL等3個指標值最小，模型的模擬效果最好，而AUC值最大，達到0.852。驗證數據也支持上述結論。最佳閾值的確定能夠同時提高樹木枯損和存活的預測準確性。