應用貝葉斯模型平均法建立東北云冷杉針闊混交林蓄積生長模型1)

賈勃 王新杰

(北京林業大學，北京，100083)

針闊混交林是一種典型的森林植被類型，在東北地區森林生態系統中具有不可替代的作用，準確預估針闊混交林蓄積量對森林經營具有現實意義。隨著模型的愈加成熟，林業研究人員根據自變量和因變量的關系采用了多種模型方法進行預估蓄積量，參數回歸方法大致分為線性方程和非線性方程，具有解釋性強的優點；非參數回歸方法(例如機器學習)簡單易操作，但存在計算量大，對樣本數要求高，易發生過擬合的缺點。洪奕豐等[1]、劉瓊閣等[2]基于偏最小二乘法估測蓄積量；閔志強等[3]采用多元線性回歸的方法對云南松的蓄積量進行估算；唐文靜等[4]采用隨機森林的方法估測森林蓄積量及其動態變化；周宇飛等[5]、曾霞輝[6]采用無人機影像數據估測森林蓄積量；王海賓等[7]和李世波等[8]基于高分一號遙感估測森林蓄積量；曾偉生等[9]、袁鈺娜等[10]利用激光雷達數據估測森林蓄積量。

然而不斷發展的這些蓄積量模型，仍無法很好地解決模型存在較大的不確定性的問題[11]。解決模型不確定性的方法主要是將模型誤差歸結為模型參數的不確定性，進而以觀測數據為基準進行模型參數優化。但這種方法也忽略了模型結構的不確定性，因此，針對模型結構的不確定性，貝葉斯模型平均法(BMA)被眾多研究者采納應用，被認為可以解決模型的不確定性問題，較單一模型也能夠提高模型模擬和預測的精度[12]。貝葉斯模型平均法是利用多模式集合進行概率預報的處理方法，它允許枚舉和評估大量可能的模型，以清楚地表示它們與數據和參數的匹配程度。與典型的逐步回歸方法相比，貝葉斯模型平均法降低了在變量選擇上產生偏差的可能性。大量研究表明，貝葉斯模型平均法的預測結果比單一模型的結果更為可靠。貝葉斯模型平均法在醫學[13]、水利[14]等方面得到了廣泛應用。王慧亮等[15]運用貝葉斯模型平均法進行了河流污染點的模擬，巴歡歡等[16]、劉家琳等[17]、赫淑杰等[18]應用貝葉斯模型平均法對洪水發生的概率進行了預測；Lu et al.[19]利用貝葉斯模型平均法構建了杉木人工林樹木死亡率與立地、競爭和氣候因素的關系；白海峰等[20]運用貝葉斯模型平均法構建森林火災預測模型，有效提高林火預測的精度；王震等[21]構建了福建杉木人工林林分蓄積量生長模型，模型表現優于逐步回歸方法。

本文應用貝葉斯模型平均法針對東北云冷杉針闊混交林進行林分蓄積生長模擬，分析林分蓄積量與林分因子、氣候因子的關系，以期為東北云冷杉針闊混交林的合理經營提供依據。

1 研究區概況

研究地點位于吉林省汪清林業局金溝嶺林場，東經130°5′～130°20′，北緯43°17′～43°25′，屬吉林省東部山區長白山系老爺嶺山脈雪嶺支脈(見圖1)。地貌屬于低山丘陵，海拔300～1 200 m，平均坡度10°～25°。該地區屬溫帶大陸性季風氣候，年均氣溫3.9 ℃。1月份氣溫最低，平均最低氣溫-32 ℃；7月份氣溫最高，平均最高氣溫32 ℃；全年無霜期達138 d；年降水量600～700 mm，降水多集中在7月。土壤類型主要為暗棕壤，土層平均厚度40 cm。主要樹種有水曲柳(Fraxinusmandshurica)、魚鱗云杉(Piceajezoensis)、楓樺(Betulacostata)、紫椴(Tiliaamurensis)、紅松(Pinuskoraiensis)和色木槭(Acermono)等。

圖1 調查樣地位置示意

2 研究方法

2.1 數據來源

數據來源于2003、2008、2012、2017年對金溝嶺林場固定樣地部分調查數據，樣地大小均為20 m×20 m。對樣地內胸徑大于5 cm的活立木進行每木檢尺，記錄樹種名稱、胸徑、林分密度、坡向、坡度、海拔、土壤類型、坐標位置等信息。沿著坡的方向記錄坡向，從坡上看坡下，表達方式為：以正北為0°，順時針讀度數，直到360°。經過選擇符合要求的樣地并剔除缺失數據的樣地后，最終確定198塊固定樣地作為研究對象。

本文選擇林分因子(林分平均胸徑、林分斷面積、林分密度)、地形因子(海拔、坡度、坡向)和氣候因子作為構建模型的初始變量。

利用汪清地區一元材積表法[22]獲得單木材積，匯總單木材積可得樣地蓄積量，進一步推算出每公頃林分蓄積量(M)。樣地林分因子和地形因子統計量見表1。

表1 林分因子統計表

氣候變量是根據各樣地經緯度坐標及海拔數據利用Climate AP軟件[23]獲得，本文獲取的因子主要有年平均溫度(MAT)、年均降雨量(MAP)、最暖月平均溫度(MWMT)，最冷月平均溫度(MCMT)、平均氣溫差(TD)、水分指數(AHM)、<0 ℃積溫(DD_0)、>5 ℃積溫(DD5)、<18 ℃積溫(DD_18)、>18 ℃積溫(DD18)、無霜期(NFED)、降雪量(PAS)、參考蒸發(EreF)、氣候水分虧損(CMD)等，氣候因子統計見表2。

表2 氣候因子統計值

統計量>5℃積溫/℃<18℃積溫/℃>18℃積溫/℃無霜期/d年降雪量/mm參考蒸發量/mm氣候水分虧損/mm平均值1848.195363.47157.28171.7953.40701.78203.68標準差47.82215.4830.764.4519.2621.0571.29

2.2 模型構建

2.2.1 逐步回歸(SR)

逐步回歸是回歸分析方法中常用的一種方法，運用逐步回歸法進行建模，即選擇一些線性關系較強的因子與因變量進行建模，從而得到較為理想的估測模型。逐步回歸的一般表達式為：

Y=a0+a1x1+a2x2+a3x3+…+anxn+ε。

式中：Y為解釋變量，a0為常數項，a1、a2、…、an為回歸系數，x1、x2、…、xn為自變量，n為自變量個數。

逐步回歸模型根據模型赤池信息量準則(AIC)值的變化及變量顯著與否，剔除與因變量不顯著的因子，最終確定模型形式。

2.2.2 貝葉斯模型平均法(BMA)

貝葉斯模型平均法(BMA)是由Rafery et al.[23]提出的一種利用多模型組合進行概率預報的統計方法，通過估算潛在變量(X)的所有可能組合的模型構建加權平均模型。在貝葉斯模型平均法中，模型的數量相當大，如果有n個自變量，模型個數為2n。貝葉斯模型平均法通過選擇可能模型的子集，并使用這些模型的模型后驗概率執行所有推斷和預測。貝葉斯模型平均法指定先驗分布，以模型后驗概率(PMP)為權重對可能的單項模型進行加權平均，以解釋變量的后驗包含概率(PIP)的大小作為選擇解釋變量的客觀標準。

一般認為，P<0.5表示該變量對因變量沒有影響；0.50≤P<0.75表示該變量對因變量有弱的影響；0.75≤P<0.95表示該變量對因變量有較強的影響；P≥0.95表示該變量對因變量有極強的影響。

2.2.3 模型評價

本研究評價模型的擬合效果和預估精度采用決定系數(R2)和均方根誤差(RMSE)。計算公式如下：

本研究選取數據的70%作為模擬數據，剩余的30%作為預測數據。分別使用逐步回歸法和貝葉斯模型平均法構建林分蓄積量模型。所有的計算均在R軟件中完成，利用step函數進行逐步回歸，利用貝葉斯自適應抽樣包BAS包實現貝葉斯模型平均法模型的構建[24]。在貝葉斯模型平均法的計算過程中，參數的先驗分布采用Zellner-Siow柯西先驗分布[25]。

3 結果與分析

3.1 變量選擇

本文采用14種氣候因子及7種林分因子進入模型擬合過程。貝葉斯模型平均法和逐步回歸選擇的變量如表3所示。采用逐步回歸方法最終確定坡向、海拔、林分密度、林分斷面積、林分平均胸徑、最暖月平均溫度、<0 ℃積溫、<18 ℃積溫、參考蒸發量等9個因子為模型影響因子，貝葉斯模型平均法(BMA)前5個模型選擇的變量各不相同，以此建立與蓄積量的關系。

表3 模型選擇

3.2 模型參數估計

逐步回歸根據變量顯著與否得出模型最佳形式，在一定程度上忽略了模型本身的不確定性。BMA模型考慮了所有可能的變量組合，考慮了林分蓄積量模型的不確定性。分別給出逐步回歸和BMA法的參數估計值。

由表4、表5可知，在逐步回歸方法中，林分蓄積量與坡位、林分斷面積、平均胸徑、最暖月平均溫度、<18 ℃積溫呈顯著正相關關系，與海拔、林分密度、<0 ℃積溫、參考蒸發呈顯著負相關關系。在BMA方法中，林分蓄積量與林分密度呈極顯著負相關關系(后驗概率P>0.95)，與林分斷面積呈極顯著正相關關系(后驗概率P>0.95),與海拔、<0 ℃積溫、參考蒸發呈弱負相關關系(后驗概率P<0.95)，與平均胸徑、<18 ℃積溫呈弱正相關關系(后驗概率P<0.95)。

表4 逐步回歸的參數估計

表5 基于貝葉斯模型的參數估計值及后驗概率

通過兩個模型的對比評價可知，逐步回歸的決定系數(R2)為0.95，均方根誤差(RMSE)為17.53。而BMA的決定系數(R2)也為0.95，均方根誤差(RMSE)為37.51。

圖2模型按照后驗概率從左側最好到右側最差排序，排名在頂部X軸上。模型中排除的變量以黑色顯示，而包含的變量則以顏色顯示，顏色與對數后驗概率相關。每列的顏色與該模型的后驗概率(X軸)的對數成正比。由此可知，林分蓄積量模型并非只是一個簡單的模型，而是由多種模型組合。

Intercept.截距；aspect.坡位；N.林分密度；Dg.林分平均胸徑；MWMT.最暖月平均溫度；TD.平均氣溫差；AHM.水分指數；DD5.>5 ℃積溫；DD18.>18 ℃積溫；PAS.降雪量；CMD.氣候水分虧損。

圖3可知，根據BMA方法預測的混交林林分蓄積量殘差波動相對平穩，逐步回歸預測的林分蓄積量殘差波動較大。

圖3 貝葉斯模型和逐步回歸法的林分蓄積量殘差

4 討論

本文應用逐步回歸法和貝葉斯模型平均法建立了東北針闊混交林蓄積生長模型，對比分析了逐步回歸和貝葉斯模型平均法對林分蓄積量的預估效果，進一步分析了貝葉斯模型平均法在模型組合方面的解釋的能力。

本文采用了氣候因子和林分因子共計21個變量進入模型。在逐步回歸中，最終選擇坡位、海拔、林分密度、林分斷面積、林分平均胸徑、最暖月平均溫度、<0 ℃積溫、<18 ℃積溫、參考蒸發9個因子進入模型。在BMA方法中，后驗概率高于0.5的只有海拔、林分密度、林分斷面積、平均胸徑、<0 ℃積溫、<18 ℃積溫、參考蒸發量等7個變量，說明這7個變量是影響林分蓄積量生長的主要因素。胸高斷面積和林分平均胸徑與云冷杉針闊混交林蓄積量呈正相關關系，與前人研究結果一致[26-27]。但是林分密度、海拔與蓄積量卻呈現負相關關系。這是因為隨著海拔增加氣溫降低，導致樹木生長緩慢，林分中小樹生長頻繁，林分密度不再主導林木生長；由于林木間的競爭關系，隨著林分密度的逐漸增大，林木發育逐漸受到抑制，蓄積呈現下降趨勢。氣候因子方面年均氣溫溫和年降水量對林木蓄積生長影響不顯著，因為東北獨特的氣候條件，巨大的晝夜溫差等間接影響林木生長，所以導致<0 ℃積溫、<18 ℃積溫等與溫度相關的因子同時可以影響林木生長發育。因此，在混交林的經營中，應合理調整林分密度及樹木徑級比例，從而促進樹木生長發育。

貝葉斯模型平均法的預測準確率和逐步回歸基本一致，一方面是由于逐步回歸更容易選擇冗余變量，會簡單地提升模型擬合效果，但貝葉斯平均法選擇的變量更加準確，所選變量更加可靠；另一方面由于總樣本量太大，貝葉斯模型平均法的預測準確率沒有充分體現，對于小樣本數據，貝葉斯方法的性能會顯著優于傳統方法[28]。另外，逐步回歸所確定的模型與貝葉斯模型平均法所確定的前5種模型選擇變量不盡相同，因為逐步回歸法相對于貝葉斯模型平均法更易選擇冗余變量[29]，且逐步回歸法忽略了模型的不確定性，造成結果出現差異。貝葉斯模型平均法通過考慮模型的不確定性，將不顯著因子剔除，模型結果更加準確可靠。

通過考慮模型的不確定性，貝葉斯模型平均法使模型更加準確。在不顯著降低預測精度的情況下，應該考慮如何利用貝葉斯模型平均法更全面地包括所有影響因素且不出現冗余現象是未來研究的方向，同時未來應該偏向小樣本數據且可拓展到森林生物量、碳儲量等方面進行研究。