賈勃 王新杰

(北京林業(yè)大學(xué)，北京，100083)

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

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

本文應(yīng)用貝葉斯模型平均法針對(duì)東北云冷杉針闊混交林進(jìn)行林分蓄積生長(zhǎng)模擬，分析林分蓄積量與林分因子、氣候因子的關(guān)系，以期為東北云冷杉針闊混交林的合理經(jīng)營(yíng)提供依據(jù)。

1 研究區(qū)概況

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

圖1 調(diào)查樣地位置示意

2 研究方法

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

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

本文選擇林分因子(林分平均胸徑、林分?jǐn)嗝娣e、林分密度)、地形因子(海拔、坡度、坡向)和氣候因子作為構(gòu)建模型的初始變量。

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

表1 林分因子統(tǒng)計(jì)表

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

表2 氣候因子統(tǒng)計(jì)值

統(tǒng)計(jì)量>5℃積溫/℃<18℃積溫/℃>18℃積溫/℃無霜期/d年降雪量/mm參考蒸發(fā)量/mm氣候水分虧損/mm平均值1848.195363.47157.28171.7953.40701.78203.68標(biāo)準(zhǔn)差47.82215.4830.764.4519.2621.0571.29

2.2 模型構(gòu)建

2.2.1 逐步回歸(SR)

逐步回歸是回歸分析方法中常用的一種方法，運(yùn)用逐步回歸法進(jìn)行建模，即選擇一些線性關(guān)系較強(qiáng)的因子與因變量進(jìn)行建模，從而得到較為理想的估測(cè)模型。逐步回歸的一般表達(dá)式為：

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

式中：Y為解釋變量，a0為常數(shù)項(xiàng)，a1、a2、…、an為回歸系數(shù)，x1、x2、…、xn為自變量，n為自變量個(gè)數(shù)。

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

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

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

一般認(rèn)為，P<0.5表示該變量對(duì)因變量沒有影響；0.50≤P<0.75表示該變量對(duì)因變量有弱的影響；0.75≤P<0.95表示該變量對(duì)因變量有較強(qiáng)的影響；P≥0.95表示該變量對(duì)因變量有極強(qiáng)的影響。

2.2.3 模型評(píng)價(jià)

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

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

3 結(jié)果與分析

3.1 變量選擇

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

表3 模型選擇

3.2 模型參數(shù)估計(jì)

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

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

表4 逐步回歸的參數(shù)估計(jì)

表5 基于貝葉斯模型的參數(shù)估計(jì)值及后驗(yàn)概率

通過兩個(gè)模型的對(duì)比評(píng)價(jià)可知，逐步回歸的決定系數(shù)(R2)為0.95，均方根誤差(RMSE)為17.53。而BMA的決定系數(shù)(R2)也為0.95，均方根誤差(RMSE)為37.51。

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

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

圖3可知，根據(jù)BMA方法預(yù)測(cè)的混交林林分蓄積量殘差波動(dòng)相對(duì)平穩(wěn)，逐步回歸預(yù)測(cè)的林分蓄積量殘差波動(dòng)較大。

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

4 討論

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

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

貝葉斯模型平均法的預(yù)測(cè)準(zhǔn)確率和逐步回歸基本一致，一方面是由于逐步回歸更容易選擇冗余變量，會(huì)簡(jiǎn)單地提升模型擬合效果，但貝葉斯平均法選擇的變量更加準(zhǔn)確，所選變量更加可靠；另一方面由于總樣本量太大，貝葉斯模型平均法的預(yù)測(cè)準(zhǔn)確率沒有充分體現(xiàn)，對(duì)于小樣本數(shù)據(jù)，貝葉斯方法的性能會(huì)顯著優(yōu)于傳統(tǒng)方法[28]。另外，逐步回歸所確定的模型與貝葉斯模型平均法所確定的前5種模型選擇變量不盡相同，因?yàn)橹鸩交貧w法相對(duì)于貝葉斯模型平均法更易選擇冗余變量[29]，且逐步回歸法忽略了模型的不確定性，造成結(jié)果出現(xiàn)差異。貝葉斯模型平均法通過考慮模型的不確定性，將不顯著因子剔除，模型結(jié)果更加準(zhǔn)確可靠。

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