常曉敏,余新曉,*,賈國棟,孫立博,劉自強,楊照國,閆騰飛

1 北京林業(yè)大學,水土保持與荒漠化防治教育部重點實驗室, 北京 100083 2 豐寧縣大灘林場, 石家莊 068357 3 張北縣林業(yè)局, 張家口 076450

近年來全球氣候變化背景下大范圍樹木死亡逐漸引起人們的廣泛關注[1-2]。地區(qū)水平上大面積的樹木死亡不僅影響碳平衡,還影響森林生態(tài)系統(tǒng)服務功能[3]。我國的三北地區(qū)氣候也正經(jīng)歷著暖干化變化[4]。該地區(qū)三北防護林經(jīng)過40多年的建設,為治理風沙、改善京津冀環(huán)境質(zhì)量發(fā)揮了重要作用[5]。但是近十年來,楊樹防護林出現(xiàn)了大面積的枯梢和死亡現(xiàn)象?？萆伊址置娣e接近80%,枯死和瀕死的約占總面積的1/3[6],防護功能明顯下降。水分是該區(qū)域植物生長的主要限制因子[7],同時粗放經(jīng)營也是導致防護林生長發(fā)育不良甚至衰退的重要原因[4]。研究楊樹林的死亡原因和影響因素,對全球氣候變化背景下三北防護林的林分經(jīng)營管理是十分必要的。

國內(nèi)外很多學者從造林方式[8]、林分結(jié)構(gòu)[9-11]、立地條件[12-13]、氣候環(huán)境條件[14-16]等多方面對樹木死亡的原因開展了大量研究。Qiu等[9]建立了林木相對斷面積、胸徑、密度、土壤理化性質(zhì)和氣候因子與樹木死亡關系模型,認為林分密度過大,或者林木相對斷面積增長量達不到一定值時,樹木死亡率增加。陳清等[17]研究了林齡和林分斷面積對樹木死亡的影響。Fan等[18]研究了胸徑、密度、坡度、坡向、海拔等因子與樹木死亡的關系,發(fā)現(xiàn)胸徑是影響樹木死亡最重要的因子。De Toledo等[11]總結(jié)了影響樹木死亡的土壤和地質(zhì)因子。以上研究表明影響樹木死亡的原因有很多,包括胸徑、密度、林木相對斷面積、土壤理化性質(zhì)等。然而不同徑級的林木,導致其死亡的原因是否相同尚不清楚。

本文以張北壩上地區(qū)的小葉楊(PopulussimoniiCarr)為研究對象,以林木胸徑、林木相對斷面積、林齡、林分密度、不同土層土壤含水量和土壤容重等11個因子作為輸入變量,運用隨機森林建立不同徑級樹木死亡模型,對影響樹木死亡的變量重要性進行評估,揭示不同徑級楊樹死亡的主要原因,這對于今后壩上地區(qū)防護林的恢復重建和經(jīng)營管護具有重要的意義。

1 材料與方法

1.1 研究區(qū)概況

試驗地點位于河北省壩上地區(qū)的張北縣,屬溫帶大陸性氣候,多年平均氣溫3.2℃,極端最高氣溫30.8℃,極端最低氣溫-30.3℃,年降水量不足400mm。每年約有50—70d風力達到六級以上,以偏西風和西北風為主,土壤風蝕嚴重。地帶性土壤為栗鈣土和沙壤土。該地區(qū)楊樹防護林自20世紀60年代開始營造,小葉楊為主要樹種。近年來,楊樹出現(xiàn)大面積的枯梢枯死現(xiàn)象。喬木樹種還有榆樹、落葉松等,主要灌木有沙棘、檸條等。

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

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

張北的防護林以楊樹純林為主。2016年7—8月和2017年7—8月,在對張北全縣18個鄉(xiāng)鎮(zhèn)進行踏查的基礎上,目測立地類型、胸徑、樹高、林齡、密度等主要調(diào)查因子后,以此為輪廓設置了100塊20m×20m的標準地進行調(diào)查,基本上涵蓋了張北防護林不同立地類型、林分密度、林齡、胸徑分布范圍等。當調(diào)查對象為帶狀防護林時,為減少邊緣效應,設置矩形標準地,并根據(jù)實際情況改變標準地的大小,使得標準地內(nèi)的林木特征和立地條件盡量一致。

記錄標準地的地理位置、海拔、地形等因子。由于張北的防護林都在平坦區(qū)域,因此沒有記錄坡度、坡向、坡位等信息。對標準地內(nèi)的楊樹進行每木檢尺,記錄胸徑、樹高、林木枯梢情況、林木死亡情況等。楊樹防護林大多為同齡林,因此每塊標準地選取3株標準木,應用生長錐測定林木的年齡,取其平均值代表該林分的年齡。在每個標準地內(nèi)按照一定方向選取3個有代表性的位置挖掘土壤剖面,分別在每一剖面按照0—10、10—20、20—60cm采集土壤樣本,并把同一層次土壤樣本進行混合。利用烘干法測定土壤含水量,環(huán)刀法測定土壤容重、孔隙度等土壤物理性質(zhì)。由于土壤含水量是變化值,因此對樣地進行了多次調(diào)查并取其平均值。該地區(qū)降水頻率和強度都較小,為研究一般情況下土壤含水量對樹木死亡的影響,調(diào)查時間選擇在長期未降水的時期,避免降水對土壤含水量造成的影響。

1.2.2數(shù)據(jù)處理

通過實地調(diào)查發(fā)現(xiàn),張北的楊樹林以成熟林和過熟林居多,即大徑級林木較多。為了使不同徑級的林木樣本量相近,去除掉部分大徑級的林木和部分生長健康無枯死木的林分,只保留部分生長健康的林分,最終選取26塊樣地,基本涵蓋了張北地區(qū)楊樹林的徑級范圍,并考慮到不同的立地類型和林分密度等因子。

研究顯示與樹木死亡有顯著相關性的因子有林木徑級大小、林齡、密度、土壤理化性質(zhì)、立地等[9,18]。張北的楊樹防護林分布于平坦區(qū)域,海拔的分布范圍較窄。研究表明張北防護林的退化死亡情況與土壤養(yǎng)分沒有明顯的相關關系[6]。因此選擇表征林齡、林木徑級大小、林木疏密程度、土壤物理性質(zhì)的因子。表征林木徑級大小的因子有胸徑和林木相對斷面積,表征林木疏密程度的因子有林分密度和林分斷面積。土壤物理性質(zhì)包括0—20cm土壤含水量、20—60cm土壤含水量、土壤平均含水量、0—20cm土壤容重、20—60cm土壤容重、土壤平均容重。林木相對斷面積即林木斷面積與林分平均斷面積的比。

應用主成分分析(PCA)的方法對11個解釋變量降維,經(jīng)過降維后得到4個主成分,累計方差貢獻率為87%。PC1主要由土壤含水量和土壤容重決定,說明土壤物理性質(zhì)；PC2主要由林木相對斷面積決定；PC3主要由密度決定；PC4說明林齡。將林木分為5個徑級(<10cm、10—15cm、15—20cm、20—25cm、>25cm),利用邏輯斯蒂廣義線性混合效應模型分別不同徑級建立樹木死亡與4個主成分因子的回歸模型。因變量為二元分類變量,枯死記為1,無枯死記為0。由于空間自相關性的影響,距離相近的林木其生長與死亡可能互相影響,因此將林木所屬的樣地號作為模型的隨機效應。固定效應系數(shù)b表示模型中固定效應對因變量影響大小的參數(shù),將固定效應系數(shù)和模型顯著性相近的劃分為一類,以此作為劃分徑級大小的標準。

應用隨機森林研究解釋變量對樹木死亡的相對重要性,建立不同徑級樹木死亡預測模型。

所有的統(tǒng)計分析應用R軟件實現(xiàn)。

1.3 研究方法

1.3.1邏輯斯蒂廣義線性混合模型

設樹木死亡的概率為P(二項分類因變量Y=1),樹木無枯死的概率為1-P(二項分類因變量Y=0)。對P進行Logit變換,

(1)

概率P與解釋變量的關系為

(2)

1.3.2隨機森林算法

隨機森林算法是由Breiman于2001提出的一種基于分類回歸樹的機器學習算法。如果因變量有n個觀測值,有k個解釋變量與之相關,它利用Bootstrap重抽樣的方法從n個原數(shù)據(jù)中隨機抽取n個觀測值,選擇部分變量作為分類樹節(jié)點,從而構(gòu)建幾百甚至幾千個分類樹組合成隨機森林,將重復程度最高的樹作為最終結(jié)果。利用隨機森林算法可以對解釋變量的相對重要性進行評價,對解釋變量進行預測。平均準確率降低度是衡量將一個變量的取值變?yōu)殡S機數(shù)后隨機森林預測準確性降低程度的指標,該值越大,說明該變量的重要性越大。

2結(jié)果與分析2.1不同徑級樹木死亡與環(huán)境因子關系

根據(jù)環(huán)境因子對各徑級樹木死亡影響的相似性,可將樹木徑級劃分為<10cm、10—25cm、>25cm(表1)。各徑級林木死亡率與土壤含水量、土壤容重(PC1)和林齡(PC4)的關系較弱。徑級為10—25cm的林木,其死亡率與林木相對斷面積(PC2)具有極顯著的正相關關系(P<0.001),與密度(PC3)具有顯著的負相關關系(P<0.05),固定效應系數(shù)分別為1.011和-0.499。徑級<10cm和>25cm的林木,死亡率與解釋變量之間無顯著的相關關系。徑級為10—25cm時,固定效應解釋的變異最大,R2為0.297。不同徑級的回歸模型中,固定效應解釋的變異小于隨機效應解釋的變異。

表1 不同徑級樹木死亡對影響因子的響應

Table 1 Results of logistic generalized linear mixed-effects model(GLMM)relating gradients of four axes to tree mortality in five DBH size classes

徑級Classes of DBH /cm固定效應系數(shù)Fixed effects coefficients (b)PC1PC2PC3PC4固定效應R2Marginal R2隨機效應R2Random effect R2<10-0.0250.355-0.2240.3660.1050.5210—25-0.1271.011???-0.499?0.4520.2970.322>25-0.4861.997-1.9820.1710.1730.766

*P<0.05；***P<0.001;PC1主要由土壤含水量和土壤容重決定；PC2主要由林木相對斷面積決定；PC3主要由密度決定；PC4主要由林齡決定;DBH: 胸徑Diameter at breast height

2.2 隨機森林模型預測精度

運用隨機森林分別不同徑級建立樹木死亡預測模型。表2的混淆矩陣顯示了隨機森林模型的預測結(jié)果。不同徑級模型預測的總體誤判率為10.7%—31.25%,解釋了樹木死亡70%—90%的變異。其中,模型對于楊樹防護林無枯死的誤判率為4.8%—13.1%,對于枯死的誤判率為30.2%—44.8%。預測精度較為滿意。

表2 隨機森林模型混淆矩陣

2.3 不同徑級樹木死亡影響因子重要性評估

分別不同徑級,運用隨機森林評估11個影響因子對樹木死亡的相對重要性,并比較哪些變量的影響是相似的。由圖1可知,胸徑<10cm時,林木相對斷面積對死亡的影響遠遠大于其他因子,其次為胸徑和0—20cm土壤水分,密度和林分斷面積的影響相對較小。當林木相對斷面積<0.22或>0.45時,樹木死亡率較大(圖2)。

胸徑為10—25cm時,密度對樹木死亡的影響最大,其次20—60cm土壤水分、林木相對斷面積、林分斷面積對樹木死亡的影響相似。當密度>600株/hm2時,樹木死亡率隨著密度的增加明顯上升(圖2)。

胸徑>25cm時,20—60cm土壤水分和林分斷面積對樹木死亡的影響最大且相近,其次為林分密度、林木相對斷面積、0—20cm土壤容重的影響相近,胸徑對樹木死亡的影響最小。當20—60cm土壤質(zhì)量含水量>5%時,樹木死亡率明顯下降。

圖1 影響樹木死亡的自變量的重要性排序Fig.1 Variable importance measure 1：密度Density；2：林木相對斷面積Relative basal area；3：林分斷面積Stand basal area；4：胸徑Diameter at breast height；5：林齡Stand age；6: 0—20cm土壤容重，0—20 cm soil bulk density；7: 20—60cm土壤容重，20—60 cm soil bulk density；8：土壤平均容重Average soil bulk density；9: 0—20cm土壤水分，0—20 cm soil moisture content；10: 20—60cm土壤水分，20—60 cm soil moisture content；11：平均土壤水分Average soil moisture content。

圖2 不同徑級相對重要變量對樹木死亡的影響Fig.2 Partial effect of variables on tree mortality in different tree size classes

3 討論

3.1 不同徑級樹木死亡對環(huán)境因子的響應

本文的研究表明林木相對斷面積、密度和土壤水分對該區(qū)域樹木死亡的影響較大,可能與該區(qū)域林木可獲得的資源及受到的大風、凍害等惡劣的自然環(huán)境有關。林分中每株樹木獲得的資源一方面由其獲得資源的能力即競爭力決定,一方面由環(huán)境可提供的資源多少決定。胸徑和林木相對斷面積較大時,林木競爭資源的能力較大[19]。林分密度和林分斷面積較大時,環(huán)境可提供給每株樹木的資源較少[20]。本文的研究顯示不同徑級的林木,林木相對斷面積對死亡的影響大于胸徑對死亡的影響,說明林木相對斷面積可以更好的描述林木之間的競爭力,相對越大的樹木競爭力越強,這與Luo等[21]的研究結(jié)果相近。

不同徑級的林木,各影響因子對樹木死亡的影響權(quán)重不同。胸徑<10cm時,林木相對斷面積對樹木死亡的影響最大,其次為胸徑和0—20cm土壤水分,密度和林分斷面積的影響最小。說明該階段樹木大小是影響樹木死亡最主要的因子。Zhu等[12]的研究也發(fā)現(xiàn)樹木大小對小徑級樹木的死亡有顯著影響。本研究中當0.23<林木相對斷面積<0.43時,樹木死亡率最低(圖2)。當林木相對斷面積<0.23時,樹木死亡率較高。壩上地區(qū)屬半干旱地區(qū),水資源嚴重缺乏,當林木相對較小時,一方面其競爭資源的能力較弱,包括水資源、土壤資源、光照資源等。楊樹屬于速生物種,木材密度較低[22],若得不到充足的水分、營養(yǎng)物質(zhì)和光照,將難以維持碳平衡,根系和主干的碳儲存減少,樹干密度降低,抵御病蟲害的能力隨之降低,死亡率增加。另一方面,周圍高大的林木會造成陰影,不利于矮小的林木獲得光照。同時,有研究表明當樹木年均斷面積增長量達不到一定值時,死亡率會增加[9]。本研究中當林木相對斷面積>0.43時,林木死亡率增大。原因可能是壩上地區(qū)大風天氣頻繁,相對較高的林木更易受大風凍害的影響,大風造成落葉增加,葉損失會減少樹木光合作用產(chǎn)物,增加非結(jié)構(gòu)性碳消耗,減少樹木體內(nèi)的碳儲備,當非結(jié)構(gòu)性碳消耗到一定程度時會導致樹木因碳饑餓而死亡[23]。這與Coomes[24]和De Toledo[11]的研究結(jié)果相似,大徑級林木更易受大風、凍害等外界環(huán)境干擾,死亡率相對較高。

胸徑為10—25cm時,密度對樹木死亡的影響最大,其次為20—60cm土壤水分、林木相對斷面積、林分斷面積對樹木死亡的影響相近。說明林分所能提供的資源多少成為該階段樹木生長的主要限制因子。由于根系的生長,這一階段林木主要利用20—60cm土壤水分。苗博等人[25]的研究中,在生長季節(jié)胸徑在該范圍的楊樹主要利用30—80cm土壤水分,與本文的研究結(jié)果相似。隨機森林的結(jié)果顯示,在這一階段密度越大,林木的死亡率越高。因為密度越大,來自周圍林木的資源競爭越強,林分所能提供給單株林木的資源越少,同時林木也會缺乏生長所需的空間。說明這一階段的楊樹都已成林,需要進行疏伐。然而邏輯斯蒂廣義線性模型的擬合結(jié)果顯示,這一階段密度與樹木死亡呈現(xiàn)負相關關系。可能是由于本研究是對樹木死亡和4個主成分之間進行了擬合回歸,PC3不能完全代表密度。同時也可能是由于該區(qū)域在冬春季節(jié)大風、凍害頻繁,當林分的密度較大時,有利于林分抵御不利的自然條件。因此該徑級范圍內(nèi),密度與樹木死亡的關系還需進一步的探究,從而得到樹木生長的最適密度。

胸徑>25cm時,20—60cm土壤水分和林分斷面積對樹木死亡的影響最大且相近。說明該階段林木主要利用20—60cm土壤水分,且林分所能提供的水資源多少成為導致樹木死亡的主要因子。在該階段,密度是導致樹木死亡的第三個主要因子,其影響小于林分斷面積,因為林分斷面積是一個隨著林分生長而變化的密度指標,不僅能表示林木的擁擠程度,還能反應林分對資源的需求。

在本研究中林齡對不同徑級樹木死亡的影響都較小,影響樹木生長的主要原因是水資源缺少,大風、凍害等惡劣的自然環(huán)境,以及人類不合理的經(jīng)營。這與部分學者認為的壩上地區(qū)樹木死亡是樹齡老化的觀點不一致[26],在以后需要對楊樹是否達到老齡化這個問題進行進一步的研究。

3.2 模型預測精度

本研究運用隨機森林算法建立的樹木死亡預測模型可以解釋樹木死亡70%—90%的變異,Botter等[27]、De Toledo等[11]、Luo等[21]建立的簡單線性模型或回歸樹模型可以解釋樹木死亡20%—80%的變異,說明本研究中建立的樹木死亡模型擬合效果較好。邏輯斯蒂回歸對自變量的多元共線性非常敏感,要求自變量之間相互獨立,因此本研究在運用模型之前將11個影響因子做了主成分分析,而隨機森林算法可以評估所有變量的重要性,不需要顧慮一般回歸問題面臨的多元共線性問題。

本研究只選取了表示林分結(jié)構(gòu)和土壤性質(zhì)的幾個因子作為解釋變量,在以后的研究中可以加入氣候因子。另外本文只研究了0—60cm土壤水分與樹木死亡的關系,以后可以增加更深層次土壤水分對樹木死亡的影響。

4 結(jié)論

本研究基于隨機森林,分析了壩上張北地區(qū)導致不同徑級樹木死亡的主要影響因子,并得到了關鍵因子的影響閾值。結(jié)果表明,林木相對斷面積、密度、林分斷面積、土壤水分是影響該區(qū)域樹木死亡的主要因子。林木胸徑<10cm時,當林木相對斷面積在0.23—0.43之間時,樹木死亡率較低；胸徑為10—25 cm時,密度>600株/hm2,樹木死亡率明顯上升；胸徑>25cm,當20—60cm土壤質(zhì)量含水量>5%時,樹木死亡率明顯下降。不同徑級樹木死亡的主要影響因子不同,但都表明水資源是影響該區(qū)域樹木生長的關鍵限制因子。不同階段的林木,有其不同的適宜密度?；谒Y源承載力的林分密度動態(tài)調(diào)控是該區(qū)域防護林經(jīng)營的關鍵。要加強森林的經(jīng)營管護,且要遵循自然規(guī)律,不同生長發(fā)育階段的林木經(jīng)營措施不同,才能獲得森林的可持續(xù)發(fā)展。