基于零膨脹模型和柵欄模型的贛南杉木林林分枯損模型

劉 軍，潘 萍，，歐陽勛志，*，臧 顥，郭 楊，游景暉

（1.鄱陽湖流域森林生態(tài)系統(tǒng)保護與修復國家林業(yè)和草原局重點實驗室，江西 南昌 330045；2.江西農業(yè)大學 林學院，江西 南昌 330045）

【研究意義】枯損是森林動態(tài)變化過程中的一個重要表現(xiàn)，也是森林生長過程中不可或缺的組成部分，其受內在競爭、人為干擾以及氣候變化等因素的綜合影響，是森林演替、更新的重要內容[1-2]。分析林分的枯損情況，探討樹木死亡的影響因素，有助于提升森林質量及生態(tài)效益，為造林、育林和護林等生產活動提供科學依據[3]。【前人研究進展】由于枯損模型存在研究尺度間的差異性，因此可劃分為單木水平和林分水平兩類。當前，國內外諸多學者對于單木水平枯損模型已進行了較多研究[4-5]。雖然單木水平枯損模型的預測結果比較精確，但是需要花費大量的時間、精力和財力對單木、林分和立地因子等進行較為全面的調查[6]，而通過分析林分枯損株數(shù)與林分因子、立地因子和氣候因子等環(huán)境因子之間的關系構建林分水平枯損模型則可以較為有效地克服這一缺點[7]。考慮到林分枯損株數(shù)為典型的計數(shù)數(shù)據，因此不少學者采用泊松模型、負二項模型、零膨脹模型和柵欄模型等計數(shù)模型對林分枯損進行了相關研究，如張雄清等[7]基于吉林省汪清林業(yè)局金溝嶺林場的20塊面積在0.077 5～0.25 hm2的長白落葉松（Larix olgensis）樣地，得出柵欄負二項模型的擬合效果優(yōu)于其他模型，其林分枯損與林齡和相對植距呈負相關，與林分斷面積呈正相關。然而，目前大多研究都建立在較小范圍尺度上，而基于較大區(qū)域尺度對林分枯損的研究還相對較少，且研究得出的較優(yōu)林分水平枯損模型以及影響林分枯損的因素也不盡相同，這可能與樹種的生物學特征及研究的尺度等因素有關?！颈狙芯壳腥朦c】杉木（Cunninghamia lanceolata）是我國亞熱帶地區(qū)特有的優(yōu)良用材樹種。具有較喜光、不耐嚴寒與干旱等特點，在滿足我國木材需求、維持生態(tài)安全和支持天然林保護等方面發(fā)揮著至關重要的作用[8]。近年來，國內外不少學者對杉木的枯損影響因素展開了相關研究，如郎奎建等[9]基于江西省16 塊0.05 hm2的杉木人工林構建枯損模型，得出林分枯損株數(shù)與林分密度和立地因子密切相關；Zhang 等[10]以江西、福建、廣西和四川等省份建立的60 塊0.06 hm2樣地為研究對象，構建邏輯回歸模型，認為競爭、林齡和林分結構是導致杉木枯損的主要因素。然而，由于全球氣候多變以及經營措施的不同，在不同區(qū)域導致杉木林林分枯損的主要影響因素可能存在差異?！緮M解決的關鍵問題】為此，本研究以江西贛南（贛州市）為研究區(qū)，以地類為純林的杉木樣地為研究對象，選取零膨脹模型和柵欄模型，構建并篩選出最優(yōu)林分水平枯損模型，并分析其影響因素，以期為杉木林的有效經營管理提供參考依據。

1 材料與方法

1.1 研究區(qū)概況

贛州市（24°29′～27°09′ N，113°54′～116°38′ E）位于江西省南部，地處贛江上游，總面積3.9×104km2，占全省總面積的23.6%。地貌以山地和丘陵為主，土壤以地帶性紅壤、黃壤為主。氣候屬亞熱帶季風氣候區(qū)，年均溫為19.1～20.8 ℃，年均降水量為1 580 mm，無霜期288 d，年均日照時數(shù)1 636 h，森林覆蓋率達76.2%，主要森林類型有針葉林、常綠闊葉林等[11]。其中杉木林面積70.40 萬hm2，占有林地總面積的26.2%。

1.2 數(shù)據來源與處理

1.2.1 樣地數(shù)據與處理 樣地數(shù)據來源于贛州市2009 年森林資源二類調查，樣地面積為0.08 hm2（28.28 m×28.28 m），地理坐標以其西南角為準，對所有胸徑≥5 cm 的喬木進行每木檢尺，并調查了坡向、坡度等立地因子及林齡、平均胸徑等林分因子。參考《江西省森林資源二類調查操作細則（2009）》，基于樣地數(shù)據庫，剔除部分數(shù)據異常的樣地，篩選出地類為純林的杉木樣地1 973 塊，其樣地分布見圖1。按4∶1的比例隨機劃分為模擬樣本（1 579 塊樣地）和驗證樣本（394 塊樣地），分別用于模型的構建和驗證。通過篩選每木檢尺庫中檢尺類型為枯立木（已枯死但未倒伏的林木）的林木得出各樣地的枯損株數(shù)，其林分枯損株數(shù)頻數(shù)統(tǒng)計見圖2，由圖2可知，本研究數(shù)據結構中含有大量的零值，即處于零膨脹狀態(tài)，數(shù)據離散程度較大。

圖1 研究區(qū)位置及樣地分布Fig.1 Location of the study area and sample plots distribution

圖2 林分枯損株數(shù)頻數(shù)統(tǒng)計Fig.2 Frequency statistics of the number of stand-level mortality stems

1.2.2 氣候數(shù)據 利用ClimateAP v2.30軟件[12]，通過輸入樣地的經緯度和海拔，獲取各樣地2000—2009年的各類氣候因子數(shù)據，并對其求取年平均值。

1.3 研究方法

1.3.1 枯損影響因子選擇 根據樣地調查內容及參考相關文獻[6,13]，從林分、立地和氣候等方面選擇可能影響林分枯損的17個環(huán)境因子，其中林分因子包括林齡、平均樹高、斷面積和株數(shù)密度等4 個因子；立地因子包括海拔、坡向、坡位、坡度、土層厚度和腐殖質層厚度等6 個因子，且坡向、坡位等定性因子參考邵方麗等[14]的方法將其賦值量化，即用1、2、3、4、5 分別代表陰坡、半陰坡、半陽坡、陽坡和無坡向，用1、2、3、4、5、6、7分別代表山脊、上坡、中坡、下坡、山谷、平地和全坡；氣候因子包括年平均溫度、年平均降雨量、最暖月平均溫度、最冷月平均溫度、無霜期天數(shù)、前一年8 月至當年7 月的降雪量和夏季（6—8月）平均降雨量等7個因子。

1.3.2 林分水平枯損模型選擇 本研究結合數(shù)據情況，選擇零膨脹模型（Zero-inflated）和柵欄模型（Hurdle）構建林分水平枯損模型。

（1）零膨脹模型。零膨脹模型可以有效解決數(shù)據過度離散和零值過多等問題[15-16]，其數(shù)據結構中的所有零值既來源于Logit部分，也來源于截尾的泊松（Poisson）部分或負二項（NB）部分。模型主要由兩部分組成：第一部分為邏輯斯蒂模型（Logit），第二部分為Poisson 模型或NB 模型。前者大多用于解釋林分發(fā)生枯損的可能性（概率），即對應模型的零部分；后者主要用于預測林分枯損株數(shù)，即對應模型的離散部分?；谀Ｐ偷诙糠执嬖诘男问讲町悾煞譃榱闩蛎洸此赡Ｐ停╖ero-inflated-Poisson，ZIP）和零膨脹負二項模型（Zero-inflated-NB，ZINB）。其計算原理及方法詳見文獻[17]。

（2）柵欄模型。柵欄模型又稱兩部分模型，其優(yōu)點及對應形式與零膨脹模型相似，但其數(shù)據結構中的所有零值只來源于Logit部分。模型主要由2個階段構成：第一個階段判斷林分是否發(fā)生枯損，只有當林分發(fā)生枯損時，才會進入到第二個階段；第二個階段主要利用Poisson 分布和NB 分布預測林分枯損株數(shù)。同理，可分為柵欄泊松模型（Hurdle-Poisson，HP）和柵欄負二項模型（Hurdle-NB，HNB）。其計算原理及方法詳見文獻[18]。

1.3.3 模型評價 采用赤池信息準則（AIC）、貝葉斯信息準則（BIC）和-2 倍對數(shù)似然函數(shù)值（-2logL）比較各模型的擬合效果，其值越小，表明模型的擬合效果越好。利用驗證樣本檢驗模型精度時，選擇平均絕對誤差（MAE）和均方根誤差（RMSE），其值越小，表明模型的預測效果越好。

1.3.4 數(shù)據分析方法 數(shù)據的預處理在Excel 2016軟件中進行。利用R 4.1.2軟件完成模型的構建和驗證，首先通過R 4.1.2 的回歸診斷程序包car 中的VIF 函數(shù)進行多重共線性診斷；其次采用caret 包中的createDataPartition 函數(shù)將數(shù)據拆分為訓練樣本及測試樣本；然后利用pscl包中的zeroinfl和hurdle函數(shù)構建林分水平枯損模型；最后使用AIC函數(shù)、BIC函數(shù)和predict函數(shù)等進行模型評價。運用Origin 2019b軟件進行相關圖表的繪制。

2 結果與分析

2.1 多重共線性診斷

為提高模型模擬精度，對選取的17個可能影響林分枯損的環(huán)境因子，需盡可能地排除彼此間相關性較大的因子。本研究的篩選原則為，當1對因子間存在共線性時，只保留其中一個因子，優(yōu)先保留有研究表明對林分枯損影響顯著的因子，篩選后的診斷結果見圖3。由圖3可知，除最暖月平均溫度、最冷月平均溫度和無霜期天數(shù)外，其余14 個環(huán)境因子間均不存在多重共線性（VIF＜10），因此將其用于林分水平枯損模型的構建。具體環(huán)境因子統(tǒng)計見表1。

表1 環(huán)境因子統(tǒng)計Tab.1 Statistics of environmental factors

圖3 自變量多重共線性診斷Fig.3 Multicollinearity diagnosis of independent variables

2.2 林分水平枯損模型的構建

經過模型參數(shù)的檢驗，剔除不顯著的變量，得出零膨脹模型和柵欄模型的模擬結果見表2。由表2可知，在模型的零部分，各參數(shù)估計值均在0.001 水平上顯著，表明杉木林林分是否發(fā)生枯損主要與海拔、林齡和株數(shù)密度等3個因子顯著相關，而與氣候因子均沒有顯著影響。在模型的離散部分，海拔和林齡在ZINB 模型和HNB 模型上在0.01 水平上顯著，而在ZIP 模型和HP 模型上在0.001 水平上顯著；株數(shù)密度均在0.001水平上顯著，表明株數(shù)密度對林分枯損的影響最大。海拔和株數(shù)密度的參數(shù)估計值均為正值，說明隨海拔梯度的上升和株數(shù)密度的增大，林分枯損株數(shù)增加；而林齡的參數(shù)估計值均為負值，說明隨林齡的增加，枯損株數(shù)減少。

表2 林分水平枯損模型模擬結果Tab.2 The simulation results of stand-level mortality model

零膨脹模型和柵欄模型的擬合效果見表3。由表3 可知，ZIP 模型和HP 模型的AIC 值（63 048.66）、BIC 值（63 091.58）和-2logL 值（63 032.66）相一致；而ZINB 模型和HNB 模型的值相近，但ZINB 模型的AIC值（4 841.73）、BIC值（4 890.01）和-2logL值（4 823.73）均略小于HNB模型的AIC值（4 841.76）、BIC值（4 890.04）和-2logL 值（4 823.76）。由此可見，ZINB 模型的擬合效果最好，其次為HNB 模型，ZIP 模型和HP模型較差。

表3 林分水平枯損模型擬合效果Tab.3 Fitting effect of stand-level mortality model

2.3 模型驗證

用驗證樣本進行模型驗證，即使用394個樣本進行模型擬合計算得出零膨脹模型和柵欄模型的林分枯損株數(shù)，并與驗證樣本實際發(fā)生的林分枯損株數(shù)進行比較，結果見圖4。由圖4 可以看出，ZIP 模型、ZINB模型、HP模型和HNB模型的差別不是十分明顯。

圖4 零膨脹模型和柵欄模型的預測值與實際值Fig.4 Predicted value of zero-inflated model and Hurdle model vs.observed value

進一步采用平均絕對誤差（MAE）和均方根誤差（RMSE）等2 個指標量化比較4 個模型的預測效果，結果見表4。由表4 可知，ZINB 模型的MAE 值（39.429 3）和RMSE 值（116.508 9）依次小于HNB 模型的MAE 值（39.440 4）和RMSE 值（116.531 5）、ZIP 模型的MAE 值（39.584 0）和RMSE 值（116.950 7）及HP 模型的MAE 值（39.584 0）和RMSE 值（116.950 8），說明ZINB 模型的預測效果優(yōu)于其他模型，這也和前面的模型擬合效果相一致。綜合分析，本研究得出的最佳模擬模型為零膨脹負二項模型，其模型形式及概率分布形式如下。

表4 林分水平枯損模型預測效果Tab.4 Prediction effect of stand-level mortality model

模型形式：

概率分布形式：

式（1）、（2）中：AL為海拔，AG為林齡，DOT為株數(shù)密度。

3 討論

3.1 林分枯損影響因子

本研究發(fā)現(xiàn)，海拔、林齡和株數(shù)密度是影響杉木林林分枯損的重要因子。林分因子中，林齡和株數(shù)密度對杉木林林分枯損具有顯著影響。隨林齡的增加，枯損株數(shù)減少，這與Caspersen[19]基于密歇根州森林資源清查數(shù)據得出的研究結果一致，而與王濤等[20]通過對黑龍江省江山嬌實驗林場的48 塊樣地進行研究，得出的林分年齡大體與雜種落葉松（Larix gmelinii）人工幼齡林林分枯損株數(shù)呈正相關的結果相反，這可能是因為本研究的杉木林多為人工林，隨造林年限的增長，撫育間伐等經營措施使得林分內競爭減小而促使枯損株數(shù)減少。反之，在造林初期時，由于株數(shù)密度較大，使得林木對于光照、水分等的競爭加劇而導致枯損株數(shù)增加，與Zhao等[21]的研究結果相似。

立地因子中，海拔對林分枯損具有重要影響。這可能是因為杉木在生長過程中，較喜光、喜溫暖濕潤的氣候，隨海拔梯度的上升，溫度逐漸降低，不適宜杉木的生長，這與張凌宇等[22]的研究結果一致。但也有學者得出與本研究不一致的結果，如Radcliffe等[23]通過對俄亥俄州東南部的82塊0.05 hm2固定樣地進行研究，得出北美紅櫟（Quercus rubra）的枯損與坡位、土層厚度有關；郭韋韋[24]的研究表明，林下腐殖質層厚度與色木槭（Acer pictum）幼樹枯損顯著相關，紅松（Pinus koraiensis）幼樹的枯損則主要受坡度的影響。然而在本研究中，除海拔外，其他立地因子在模型中均不顯著，這可能是因為杉木人工造林時多從木材生產考慮而選擇在立地條件較好的地塊，因此在本研究所使用的1 973 塊樣地中，有70%左右的樣地坡位為中下坡，80%以上的樣地坡度介于16°～35°，2/3 左右的樣地土層厚度大于80 cm，2/3 左右的樣地腐殖質層厚度大于10 cm，樣地在坡位、坡度、土層厚度和腐殖質層厚度之間差異不大而導致的，其難以體現(xiàn)在對枯損的影響上，這與李春明等[6]利用吉林省295塊蒙古櫟（Quercus mongolica）樣地數(shù)據得出的結果相類似。

此外，不少研究表明，氣候因子與林分枯損狀況有關。如Kim等[25]基于韓國國家森林資源清查數(shù)據，利用MSN曲線構建林分水平枯損模型，認為針葉林的枯損株數(shù)隨最冷月平均溫度的升高有增加的趨勢；Zhang等[10]得出杉木人工林的枯損株數(shù)隨年平均溫度的升高而減少。但本研究中，氣候因子對杉木林林分枯損的影響均不顯著，可能是由于本研究區(qū)杉木分布區(qū)的主要氣候因子間的差異較小，也有可能是由于海拔等地形因子折射出了一部分氣候因子的影響，如Wu等[26]認為地形因子之間的差異會間接影響森林環(huán)境中光照和降水等氣候因子的分布情況，進而影響林分的生長和枯損狀況。

3.2 模型模擬效果

本研究表明，杉木林林分水平枯損模型的擬合效果和預測效果均以ZINB 模型最好，其次為HNB 模型、ZIP 模型，HP 模型最差。在模型的擬合效果上，李春明等[6]和Li 等[27]也得到了與本研究相同的結果，而張雄清等[7]利用汪清林業(yè)局金溝嶺林場的長白落葉松林分數(shù)據，得出HNB 模型的擬合效果最好。在模型的預測效果上，與郭韋韋[24]對冷杉（Abies fabri）幼樹枯損的研究結果一致，但也有學者得出不同的結果，如Crecente-Campo 等[28]基于西班牙巴斯克地區(qū)的森林資源清查數(shù)據，結果表明ZIP 模型在預測枯立木株數(shù)上效果最好。這些研究所得的最佳林分枯損模型存在一定差異，這可能是受研究的對象、區(qū)域、尺度以及環(huán)境因子等因素的影響[29]。

4 結論

本研究以江西贛南的杉木林為研究對象，以林分枯損株數(shù)為因變量，以林分因子、立地因子和氣候因子等14個環(huán)境因子為自變量，構建杉木林林分水平枯損模型，探討影響林分枯損的主要因素。結果表明，杉木林林分水平枯損模型的擬合效果和預測效果均以ZINB 模型最好，其次為HNB 模型、ZIP 模型，HP 模型最差；海拔、林齡和株數(shù)密度是影響杉木林林分枯損的重要因子，而氣候因子對其沒有顯著影響。在今后的研究中，有待進一步考慮預估期間林木采伐、土壤理化性質和病蟲害、人為干擾等因素展開更為深入的分析，以預測林分枯損株數(shù)變化的動態(tài)過程。

致謝：亞洲開發(fā)銀行CCF（氣候變化基金）項目、江西贈款項目（0229-PRC）和江西省林業(yè)科技創(chuàng)新專項（創(chuàng)新專項[2021]33號）同時對本研究給予了資助，謹致謝意！