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

李春明

(1. 中國林業(yè)科學(xué)研究院資源信息研究所 北京 100091；2. 國家林業(yè)和草原局森林經(jīng)營與生長模擬重點(diǎn)實(shí)驗(yàn)室 北京 100091)

作為生態(tài)過程中不可或缺的組成部分，林木枯損被認(rèn)為是森林生態(tài)系統(tǒng)中極其關(guān)鍵，卻又知之甚少的森林動態(tài)過程[1]。單木枯損模型能夠?qū)γ恐陿淠镜纳L狀況進(jìn)行模擬，諸多學(xué)者利用多種數(shù)據(jù)和方法構(gòu)建了大量的單木枯損模型。傳統(tǒng)的經(jīng)驗(yàn)方程，如指數(shù)方程和Weibull方程等，主要是通過建立枯損與協(xié)變量的回歸方程來預(yù)測未來樹木的枯損概率[2-3]。由于每株樹木只包括枯損和存活兩種結(jié)局，是典型的二分類數(shù)據(jù)，上述傳統(tǒng)的經(jīng)驗(yàn)方程在估計枯損時效果并不理想。而Logistic回歸方法對于二分類數(shù)據(jù)具有先天的優(yōu)勢，與傳統(tǒng)的經(jīng)驗(yàn)方程相比，模擬林木的枯損有較好的適用性，已經(jīng)被廣泛應(yīng)用[4-5]。

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

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

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

1 數(shù)據(jù)來源

本研究選擇吉林省295塊蒙古櫟樣地數(shù)據(jù)，樣地大小為0.06 hm2。樣地中蒙古櫟為優(yōu)勢樹種，其他主要樹種還包括云杉(Picea jazoensisNakai)、白樺(Betula platyphyllaSuk.)、紅松(Pinus koraiensisSieb. et Zucc.)、水曲柳(Fraxinus mandshuricaRupr.)和冷杉(Abies nephrolepis(Trautv.) Maxim)等。調(diào)查的樣地因子主要包括起源、優(yōu)勢樹種、平均林齡、平均樹高、郁閉度、海拔、坡度、坡向及坡位等因子。單木因子包括對胸徑大于5 cm的樹木進(jìn)行每木檢尺，枯損木的標(biāo)記和記錄。外業(yè)調(diào)查結(jié)束后進(jìn)行林分公頃株數(shù)、公頃斷面積、大于對象木斷面積等指標(biāo)的計算。建模數(shù)據(jù)采用1999年和2009年觀測的林分因子數(shù)據(jù)。枯損樹木的記錄分別是1999年觀測（沒有枯損）—2004年觀測（枯損）和2009年觀測（沒有枯損）—2014年觀測（枯損）的數(shù)據(jù)。而驗(yàn)證數(shù)據(jù)采用的是2004年觀測數(shù)據(jù)。具體的單木因子、樣地因子及立地因子見表1。氣象因子利用ClimateAP v2.10軟件[20]，通過輸入各樣地的空間位置坐標(biāo)和海拔，然后獲取各樣地每年的各類氣象因子變量，最后對每5年的變量進(jìn)行平均計算。氣象因子分別對應(yīng)的是1999—2004年和2009—2014年平均數(shù)據(jù)，具體因子見表2。

2 研究方法

2.1 主要影響因子選擇

首先把單木、樣地、立地和氣象等因子加入到模型中，逐步剔除差異不顯著（p＞0.05）以及存在多重共線性問題（方差膨脹因子 ＞ 5）的影響因子。本研究主要的候選因子包括單木初始胸徑、單木初始胸徑的自然對數(shù)、單木初始胸徑的平方、單木初始胸徑與林分平均胸徑的比值、大于對象木斷面積、林分算術(shù)平均直徑、林分平方平均直徑、林分公頃斷面積、林分公頃株數(shù)、平均年齡、郁閉度、坡度（slope）、坡向（aspect）、海拔（elevation）、SIC值（SIC= tan(slope) × cos(aspect)）、CE值（CE= cos（aspect） × ln （elevation））、年平均溫度（MAT）、最暖月平均氣溫（MWMT）、最冷月平均氣溫（MCMT）、年平均降水量（MAP）、年平均夏季(5—9月)降水量（MSP）、無霜期天數(shù)（NFFD）和上一年8月至當(dāng)年7月的降雪量（PAS）等23個影響因子。

2.2 廣義線性混合效應(yīng)模型及評價方法

廣義線性混合效應(yīng)模型（GLMM）通過在線性預(yù)測部分引入隨機(jī)效應(yīng)參數(shù)而推廣的廣義線性模型（GLM）[21]。本研究基于Logistic回歸的廣義線性混合效應(yīng)模型在截距和參數(shù)系數(shù)上考慮隨機(jī)效應(yīng)，利用SAS 9.4 的PROC NLMIXED 模塊來估計廣義線性混合效應(yīng)模型的待估參數(shù)。

AIC信息準(zhǔn)則、BIC信息準(zhǔn)則和-2*對數(shù)似然值（-2LogL）這3個值經(jīng)常被用來比較不同模型之間的擬合效果，這3個值越小，說明模型的模擬效果越好[22]。似然比卡方檢驗(yàn)（LRT）被用來比較模型之間的差異顯著程度[23]。

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

表 1 樣地各調(diào)查因子統(tǒng)計 Table 1 The statistics of plot factors

表 2 主要?dú)庀笞兞拷y(tǒng)計 Table 2 The statistical table of main climate variables

2.3 閾值確定

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

3 結(jié)果與分析

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

在Logistic模型的基礎(chǔ)上，考慮樣地水平、單木水平及兩水平（樣地內(nèi)嵌套單木水平）的隨機(jī)效應(yīng)（分隨機(jī)截距效應(yīng)和隨機(jī)系數(shù)效應(yīng)兩種情況）。在考慮單木水平時只有在截距上收斂，但是其AIC、BIC和-2LogL值分別是18 559、18 639和18 539，遠(yuǎn)遠(yuǎn)高于固定效應(yīng)模型。在考慮兩水平效應(yīng)時，無論什么情況，模型均不能夠收斂。在考慮樣地水平隨機(jī)系數(shù)效應(yīng)時，3個參數(shù)以上的隨機(jī)效應(yīng)均不收斂，考慮2個參數(shù)時，只有兩種情況收斂，第一種是同時在截距和胸徑兩個參數(shù)上考慮隨機(jī)效應(yīng)，第二種是同時在截距和大于對象木斷面積時，模型能夠收斂。具體模擬和計算的AIC、BIC、-2LogL、LRT、AUC值和參數(shù)的估計值見表3。表3中，M1為傳統(tǒng)Logistic回歸模型，M2為在截距上考慮隨機(jī)效應(yīng)的混合效應(yīng)模型，M3為在截距和胸徑變量上同時考慮隨機(jī)效應(yīng)的混合效應(yīng)模型，M4為在截距和大于對象木斷面積變量上同時考慮隨機(jī)效應(yīng)的混合效應(yīng)模型。

表 3 基于廣義線性模型方法和廣義線性混合效應(yīng)模型方法的單木枯損模型模擬結(jié)果Table 3 The result of parameter estimate based on generalized linear model and generalized linear mixed model in tree mortality

表3的模擬結(jié)果表明，當(dāng)模型考慮了樣地的隨機(jī)效應(yīng)后，AIC、BIC、-2LogL值均比沒考慮樣地的隨機(jī)效應(yīng)值要小，說明在構(gòu)建單木枯損模型時考慮隨機(jī)效應(yīng)的模擬效果要優(yōu)于不考慮隨機(jī)效應(yīng)的模型。利用LRT對不同的模型進(jìn)行方差分析，考慮樣地隨機(jī)效應(yīng)后，與固定效應(yīng)模型相比，差異均達(dá)到極顯著水平（p＜0.001）；當(dāng)同時考慮截距和隨機(jī)系數(shù)效應(yīng)后比單獨(dú)考慮隨機(jī)截距效應(yīng)，模型的模擬效果要好，差異也達(dá)到極顯著（p＜0.001）。在胸徑和截距上考慮隨機(jī)效應(yīng)的AIC、BIC、-2LogL等3個指標(biāo)值最小，說明M3的模擬效果最好。在本研究中，M3的AUC值等于0.852 ＞ M4（0.850）＞M2（0.838）＞ M1（0.779），也說明M3模型是構(gòu)建蒙古櫟單木枯損模型的最優(yōu)模型。

3.1 單木和林分因子對枯損的影響

從表3可知，M1和M2中單木初始胸徑值與枯損概率呈正比，而M3和M4中單木初始胸徑值與枯損概率呈反比，說明單純一個胸徑因子不足以表述單個林木在林分中的地位和狀態(tài)。大于對象木斷面積與枯損概率呈正比，說明大于對象木斷面積值越小，處于競爭的優(yōu)勢，不容易發(fā)生枯損，這與林分中現(xiàn)實(shí)情況相符。林分公頃斷面積與單木枯損概率呈反比，反映了在給定同樣公頃株數(shù)的林分條件下，隨著林分公頃斷面積增加，林分內(nèi)的樹木胸徑更大，處于競爭的優(yōu)勢，則枯損概率降低。林分平均年齡與單木枯損概率呈反比，說明枯損容易在幼齡林中發(fā)生，而隨著林分的生長和成熟，林木胸徑變大，枯損概率降低，林木也不易發(fā)生枯損。林分公頃株數(shù)與單木枯損概率呈正比，說明在密度大的林分，林木之間對于光、熱、水肥等營養(yǎng)的競爭更加激烈，林木容易發(fā)生枯損。

3.2 氣象因子對枯損的影響

在表3的模擬結(jié)果中，年平均溫度（MAT）、年平均夏季(5—9月)降水量（MSP）、上一年8月至當(dāng)年7月的降雪量（PAS）對枯損有明顯的影響。MAT與單木枯損概率呈正相關(guān)，同樣PAS也與單木枯損概率呈正相關(guān)。而MSP與單木枯損概率呈負(fù)相關(guān)。氣象因子的影響說明，氣候變化引起的干旱或霜凍對單木枯損有重要的影響。

3.3 閾值分析

根據(jù)相應(yīng)的閾值標(biāo)準(zhǔn)，本研究枯損株數(shù)與總林木株數(shù)的比為0.06，則割點(diǎn)閾值選擇為A= 0.06。而傳統(tǒng)閾值0.5被舍去。圖1列出了隨機(jī)閾值和錯誤分類率最小的閾值。其中基礎(chǔ)模型（M1）的隨機(jī)閾值為B= 0.08，而錯誤分類率最低閾值為C=0.07（圖1(a)）?；旌夏Ｐ停∕3）的隨機(jī)閾值為B=0.07，而錯誤分類率最低閾值為C= 0.06（圖1(b)）。

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

表4表明，基礎(chǔ)模型（M1）中，當(dāng)閾值為0.06時枯損的準(zhǔn)確分類率（TPR）最高，為81.9%，當(dāng)閾值為0.08時，存活的準(zhǔn)確分類率（TNR）最高，為72.7%；混合模型（M3）中，當(dāng)閾值為0.06時枯損的準(zhǔn)確分類率（TPR）最高，為84.1%，當(dāng)閾值為0.07時，存活的準(zhǔn)確分類率（TNR）最高，為75.4%。在考慮混合效應(yīng)后，無論閾值為0.06、0.07或0.08，枯損和存活的準(zhǔn)確分類率均高于基礎(chǔ)模型，而枯損和存活的錯誤分類率均低于基礎(chǔ)模型。通過計算總正確分類率（ACR），基礎(chǔ)模型閾值為0.08時，總正確分類率最高，為72.6%；混合模型閾值為0.07時，總正確分類率最高，為75.7%。因此可以選擇靈敏度與特異度曲線交叉時的閾值（隨機(jī)閾值）作為最佳閾值。

3.4 模擬驗(yàn)證

根據(jù)表3的模擬結(jié)果，分別利用M1和M3模型形式對驗(yàn)證數(shù)據(jù)每株樹木進(jìn)行枯損概率的計算，最后計算各自的AUC值，并比較各模型的靈敏度。其中，M3的模型形式，由于考慮了樣地的隨機(jī)效應(yīng)，則首先要求出驗(yàn)證數(shù)據(jù)中樣地的隨機(jī)效應(yīng)參數(shù)值，然后再計算枯損概率。隨機(jī)效應(yīng)參數(shù)計算方法可參考Vonesh and Chinchilli[27]。

表 4 蒙古櫟混交林單木枯損模型閾值預(yù)測列聯(lián)表分析結(jié)果 Table 4 Confusion matrix result of tree mortality model of Mongolian oak

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

4 討論

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

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

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

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

5 結(jié)論

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