異速參數(shù)先驗(yàn)信息對生物量方程擬合的影響

黃興召 豆玉萍 黃慶豐 陳東升

(1. 安徽農(nóng)業(yè)大學(xué)林學(xué)與園林學(xué)院，合肥 230036；2. 中國林業(yè)科學(xué)研究院林業(yè)研究所，北京 100091)

立木地上生物量是森林生物量的重要組成部分，如何準(zhǔn)確地估算立木地上生物量是研究森林生態(tài)系統(tǒng)物質(zhì)循環(huán)的基礎(chǔ)[1?2]。異速生物量模型作為預(yù)估生物量的常用方法，已被廣泛應(yīng)用于各類研究中[3?6]。異速生物量方程參數(shù)a和b隨著研究區(qū)域和樹種的變化而變化，通過收獲法，利用異速生物量方程預(yù)估生物量時，只有單獨(dú)建模才能提高預(yù)估精度[5?8]。貝葉斯方法作為其中的一種擬合方法，通過收集已建立方程的參數(shù)信息，將該信息轉(zhuǎn)化為參數(shù)的概率分布，以參數(shù)的均值和協(xié)方差矩陣為先驗(yàn)信息進(jìn)行擬合，越來越受到人們的關(guān)注[9?12]。但是，利用先驗(yàn)信息對貝葉斯方法進(jìn)行擬合時，研究者收集參數(shù)a和b的先驗(yàn)信息具有隨機(jī)性和異質(zhì)性，由此對模型擬合和預(yù)估效果產(chǎn)生的影響尚未得到驗(yàn)證，也沒有闡述如何選取參數(shù)a和b的先驗(yàn)信息來提升方程的擬合和預(yù)估效果。

收集已發(fā)表文獻(xiàn)異速生物量方程參數(shù)a和b的信息，是分析先驗(yàn)信息異質(zhì)性對立木生物量預(yù)估影響的基礎(chǔ)。Ter-Mikaelian和Korzukhin[13]建立了北美65個樹種的生物量預(yù)估方程、Zianis和Mencuccini[3]建立了279個異速生物量預(yù)估方程、Henry等[14]收集了850個立木及各組分的生物量預(yù)估方程、Zapata-Cuartas等[8]收集了134個的立木地上異速生物量方程、Rojas-García等[15]收集了487個立木及各組分的生物量預(yù)估方程。這些方程都是通過常規(guī)方法擬合得到參數(shù)a和b的值，為貝葉斯方法提供了豐富的先驗(yàn)信息。分析發(fā)現(xiàn)，立木地上部分生物量預(yù)估方程參數(shù)b與緯度呈顯著相關(guān)，且在不同屬之間存在顯著差異[8,13]。研究者在利用貝葉斯方法預(yù)估立木地上部分生物量的過程中，討論了樣本大小對方程擬合效果和預(yù)估效果的影響，缺乏反映樹種特征的先驗(yàn)信息的擬合比較[8,10,16]。樹木生長與生物量的積累受環(huán)境[17,18]和生物學(xué)特性[19?20]的影響，在方程擬合時先驗(yàn)信息也會受到影響[21]。因此，檢驗(yàn)參數(shù)a和b先驗(yàn)信息的氣候區(qū)和生物學(xué)特性對立木生物量預(yù)估的影響顯得尤為重要。

通過分析已發(fā)表文獻(xiàn)關(guān)于立木地上部分異速生物量方程的參數(shù)a和b信息，將參數(shù)的先驗(yàn)信息按主要?dú)夂騾^(qū)和樹種進(jìn)行分類，以日本落葉松（Larix kaempferi）立木地上部分生物量為擬合數(shù)據(jù)，利用貝葉斯方法擬合異速生物量方程，比較先驗(yàn)信息的氣候差異和樹種差異對立木生物量擬合和預(yù)估的影響，探討先驗(yàn)信息對立木生物量預(yù)估的作用規(guī)律，為利用貝葉斯方法擬合異速生物量方程時選擇先驗(yàn)信息提供了理論和實(shí)踐的依據(jù)。

1 研究方法

1.1 異速生物量方程

異速生物量方程作為預(yù)估立木生物量的常用模型，其方程表達(dá)式為：

其中，Mi為第i株樣本的地上部分生物量；Di為第i株樣本的胸徑；a和b為方程參數(shù)；ei為誤差項(xiàng)。

貝葉斯方法將公式（1）中的參數(shù)a、b定義為θ，需要知道θ的概率密度函數(shù)。接著，根據(jù)先驗(yàn)信息確定參數(shù)a和b的先驗(yàn)分布 π(θ)，再結(jié)合樣本信息可得出后驗(yàn)分布 π(θ|D)。 π(θ|D)的計(jì)算公式為：

1.2 先驗(yàn)信息

利用貝葉斯方法擬合異速生物量方程時，以參數(shù)a、b的均值和協(xié)方差矩陣為先驗(yàn)信息。為研究不同先驗(yàn)信息對異速生物量方程擬合效果的影響，本文通過分析[8,13]立木地上部分異速生物量方程參數(shù)a和b的數(shù)據(jù)集，按照氣候和屬進(jìn)行分類（表1）。以主要?dú)夂騾^(qū)為單位，將參數(shù)a和b分為寒帶、溫帶和熱帶三種類型，分別計(jì)算它們的均值和協(xié)方差矩陣。同時，以屬為單位，將參數(shù)a和b分為櫟屬（Quercus）、樺木屬（Betula）、楊屬（Populus）、槭屬（Acer）、桉屬（Eucalyptus）和松屬（Pinus）等六個屬，計(jì)算其均值和協(xié)方差矩陣。以氣候區(qū)和屬為單位進(jìn)行分類，是進(jìn)一步研究先驗(yàn)信息異質(zhì)性對異速生物量方程擬合效果影響的關(guān)鍵。

表1 異速生物量方程先驗(yàn)信息分類Table1Priorinformationclassificationofallometric biomassequation

2 方程擬合

2.1 擬合數(shù)據(jù)

以遼寧省清源縣大孤家林場（42°16?～42°23? N，124°48?～124°55? E）日本落葉松（Larixkaempeferi）立木地上部分生物量數(shù)據(jù)為擬合數(shù)據(jù)[22]，該區(qū)屬于典型的溫帶季風(fēng)氣候?；跇拥卣{(diào)查，在標(biāo)準(zhǔn)地附近選取標(biāo)準(zhǔn)木，利用收獲法測定樣木的樹干、樹皮、樹枝和樹葉各組分的生物量，求出立木地上部分生物量。本數(shù)據(jù)由幼齡林、中齡林、近熟林以及成熟林4個林齡段共45株立木的胸徑和地上部分生物量組成（表2）。本文在建模過程中，利用無重復(fù)抽樣，從45株總樣木中抽取30株樣木，進(jìn)行1 000次無重復(fù)擬合，每次擬合未被抽取的15株作為驗(yàn)證數(shù)據(jù)來檢驗(yàn)預(yù)估效果。

表2 日本落葉松樣地信息Table2PlotsinformationofLarixkaempferiplantations

2.2 預(yù)估效果與預(yù)估精度檢驗(yàn)

使用不同的先驗(yàn)信息，利用貝葉斯方法擬合異速生物量方程時，采用決定系數(shù)（R2）、平均偏差（Mean Bias , MB）和平均均方根差（Mean Root Mean Square Error, MRMSE）來檢驗(yàn)方程的擬合結(jié)果和預(yù)估效果。R2、MB和MRMSE的表達(dá)式分別為：

其中，yi為實(shí)際觀測值；yˉ為樣本平均值；y?i為模型預(yù)測值；n為樣本數(shù)。數(shù)據(jù)統(tǒng)計(jì)分析和繪圖使用R（Development Core Team, 2009）和Excel軟件。以不同區(qū)域和不同樹種參數(shù)a、b的均值以及協(xié)方差矩陣為先驗(yàn)信息，應(yīng)用R的MCMCglmm包[23]對異速生物量方程分別進(jìn)行擬合。

3 結(jié)果與分析

3.1 先驗(yàn)信息氣候區(qū)差異對異速生物量方程擬合的影響

為檢驗(yàn)先驗(yàn)信息區(qū)域差異對異速生物量方程擬合的影響，分別以寒帶、溫帶、熱帶以及總數(shù)據(jù)集參數(shù)a、b的均值及協(xié)方差矩陣為先驗(yàn)信息，利用無重復(fù)抽樣，從45株總樣木中抽取30株樣木，進(jìn)行1 000次無重復(fù)擬合（圖1）。寒帶和熱帶先驗(yàn)信息擬合后參數(shù)a、b值的變化范圍相似；溫帶和總數(shù)據(jù)集先驗(yàn)信息擬合后參數(shù)a、b值的變化范圍相似。以熱帶和寒帶為先驗(yàn)信息進(jìn)行1 000次無重復(fù)擬合，參數(shù)a的平均值均為?2.794，95%的置信區(qū)間分別為（?3.185, ?2.386）和（?3.197,?2.368）；參數(shù)b的平均值均分別為2.640和2.638，95%的置信區(qū)間分別為（2.484, 2.788）和（2.476, 2.792）。以溫帶和總數(shù)據(jù)集為先驗(yàn)信息進(jìn)行1 000次無重復(fù)擬合，參數(shù)a的平均值分別為?2.836和?2.826，95%的置信區(qū)間分別為（?3.240,?2.418）和（?3.228, ?2.407）；參數(shù)b的平均值均分別為2.655和2.651，95%的置信區(qū)間分別為（2.496, 2.809）和（2.491, 2.805）（表3）。這說明使用不同區(qū)域的先驗(yàn)信息擬合異速生物量方程會影響參數(shù)a、b值的穩(wěn)定性。

圖1 不同氣候區(qū)先驗(yàn)信息預(yù)估參數(shù)a和b的擬合結(jié)果Fig.1Simulationresultsfortheestimationofthe parametersaandbusingdifferentprior informationofregion

決定系數(shù)R2反映方程的擬合效果（表3），寒帶、溫帶、熱帶以及總數(shù)據(jù)集四種先驗(yàn)信息1 000次無重復(fù)擬合，方程的平均決定系數(shù)均為0.934。平均偏差MB和平均均方根差MRMSE反映了方程的預(yù)估效果，四種先驗(yàn)信息1 000次無重復(fù)擬合，平均偏差分別為?1.392、?0.597、?0.929、?0.715，平均均方根差分別為0.347、0.344、0.347、0.345。由此可知，區(qū)域先驗(yàn)信息主要影響方程的預(yù)估效果。

表3 不同氣候區(qū)先驗(yàn)信息擬合異速生物量方程的結(jié)果對比Table3Resultcomparisonallometricbiomassequationsindifferentpriorinformationofregions

為了進(jìn)一步比較寒帶、溫帶、熱帶以及總數(shù)據(jù)集四種先驗(yàn)信息對方程預(yù)估效果的影響，對BE和RMSE分別進(jìn)行方差分析（圖2）；總數(shù)據(jù)集先驗(yàn)信息的MB和MRMSE與溫帶的MB和MRMSE無顯著差異；寒帶與熱帶的MB和MRMSE也無顯著差異；其他先驗(yàn)信息的MB和MRMSE之間存在顯著差異（P<0.01）?？倲?shù)據(jù)集和溫帶先驗(yàn)信息的擬合效果優(yōu)于寒帶和熱帶，驗(yàn)證數(shù)據(jù)集所在區(qū)域?yàn)闇貛?，因此，利用貝葉斯方法擬合異速生物量方程時，應(yīng)充分考慮先驗(yàn)信息的氣候區(qū)特征。

圖2 不同氣候區(qū)先驗(yàn)信息擬合的平均偏差和平均均方根差比較Fig.2Simulationresultsforthemeanbiasandmean rootmeansquareerrorusingdifferent priorinformationofregion

3.2 先驗(yàn)信息屬差異對異速生物量方程擬合的影響

以櫟屬、樺木屬、楊屬、槭屬、桉屬、松屬和總數(shù)據(jù)集為先驗(yàn)信息，在樣木數(shù)量為30的情況下，也進(jìn)行1 000次無重復(fù)擬合（圖3）。桉屬、櫟屬、松屬以及總數(shù)據(jù)集的參數(shù)a的變化范圍相似，其平均值分別為?2.859、?2.866、?2.855、?2.826；樺木屬和楊屬的參數(shù)a的變化范圍相似，其平均值分別為?2.728和?2.712；槭屬的參數(shù)變化范圍尤為突出，其平均值為?2.492。參數(shù)b的變化范圍也表現(xiàn)出相同的規(guī)律，桉屬、櫟屬、松屬以及總數(shù)據(jù)集參數(shù)b的平均值分別為2.663、2.666、2.662、2.651；樺木屬和楊屬參數(shù)b的平均值分別為2.616和2.605；槭屬參數(shù)b的平均值為2.528。

圖3 不同屬先驗(yàn)信息預(yù)估參數(shù)a和b的擬合結(jié)果Fig.3Simulationresultsfortheestimationofthe parametersaandbusingdifferent priorinformationofgenus

以櫟屬、樺木屬、楊屬、槭屬、桉屬、松屬和總數(shù)據(jù)集為先驗(yàn)信息，進(jìn)行異速生物量方程擬合（表4）表4。槭屬的平均決定系數(shù)為0.935，其他屬的平均決定系數(shù)均為0.934。櫟屬、樺木屬、楊屬、槭屬、桉屬、松屬和總數(shù)據(jù)集的平均偏差分別為?0.205、?1.287、?2.926、?3.726、?0.367、?0.494、?0.715；平均均方根差分別為0.343、0.352、0.354、0.391、0.344、0.343、0.345。由此可知，樹種先驗(yàn)信息主要影響方程的預(yù)估效果。

表4 不同屬先驗(yàn)信息擬合異速生物量方程的結(jié)果對比Table4Resultcomparisonallometricbiomassequationsindifferentpriorinformationofgenus

為了進(jìn)一步比較六個屬以及總數(shù)據(jù)集七種先驗(yàn)信息對方程預(yù)估效果的影響，對MB和MRMSE分別進(jìn)行方差分析（圖4）。櫟屬和桉屬、松屬和桉屬之間的MB無顯著差異，其他屬之間的MB均存在顯著差異（P<0.01），即櫟屬和桉屬為先驗(yàn)信息的MB>松屬>總數(shù)>樺木屬>楊屬>槭屬。櫟屬和桉屬、松屬和桉屬、總數(shù)據(jù)集和桉屬之間的MRMSE無顯著差異，其他屬之間的MRMSE存在顯著差異（P<0.01），即松屬和櫟屬為先驗(yàn)信息的MRMSE>總數(shù)據(jù)集和桉屬>樺木屬>楊屬>槭屬。因此，樹種的生物學(xué)特性也是影響生物量預(yù)估的重要因子。

圖4 不同屬先驗(yàn)信息擬合的平均偏差和平均均方根差比較Fig.4Simulationresultsforthemeanbiasandmean rootmeansquareerrorusingdifferent priorinformationofgenus

4 討論

4.1 區(qū)域差異對異速生物量方程擬合的影響

利用貝葉斯方法擬合生物量方程，先驗(yàn)信息尤為重要[8?9,16]，通過收集立木地上部分異速生物量方程參數(shù)信息，分析表明參數(shù)a和b符合二元正態(tài)分布。擬合異速生物量方程過程中，貝葉斯原理使用先驗(yàn)信息，而傳統(tǒng)方法只能使用立木數(shù)據(jù)。因此，貝葉斯方法可以實(shí)現(xiàn)少量的樣木高精度預(yù)估生物量，黃興召等[10]和Zell等[16]均證明，當(dāng)樣木的數(shù)量小于50株時，貝葉斯方法的預(yù)估效果明顯優(yōu)于傳統(tǒng)方法。立木地上部分異速生物量方程參數(shù)a和b的值有明顯的區(qū)域特征，參數(shù)b與緯度呈線性關(guān)系[10]。本文將參數(shù)a和b劃分為寒帶、溫帶和熱帶三個主要?dú)夂騾^(qū)，基于信息量度的相似概率分布[22,24?25]，以三個氣候區(qū)和總數(shù)據(jù)集參數(shù)a、b的均值和協(xié)方差矩陣為先驗(yàn)信息，利用無重復(fù)抽樣，從總樣木中抽取30株立木生物量數(shù)據(jù)，進(jìn)行1 000次無重復(fù)擬合，得出每次擬合結(jié)果中參數(shù)a、參數(shù)b、R2、MB和MRMSE的值。

將寒帶、溫帶和熱帶三個主要?dú)夂騾^(qū)先驗(yàn)信息的擬合結(jié)果與總數(shù)據(jù)集先驗(yàn)信息的擬合結(jié)果進(jìn)行比較，發(fā)現(xiàn)R2均為0.934無顯著差異，說明先驗(yàn)信息的氣候差異對異速生物量方程的擬合效果無影響。MB和MRMSE差異分析表明，先驗(yàn)信息的氣候差異對異速生物量方程的預(yù)估效果有顯著影響。本文中擬合數(shù)據(jù)的研究地點(diǎn)屬于典型的溫帶季風(fēng)氣候，溫帶的先驗(yàn)信息擬合效果（顯著）優(yōu)于寒帶和熱帶，與總數(shù)據(jù)集無差異，驗(yàn)證了先驗(yàn)信息的氣候差異對立木生物量預(yù)估的影響，使用立木生長主要?dú)夂騾^(qū)的先驗(yàn)信息，有助于提升立木生物量的預(yù)估精度。

通過模型擬合，證實(shí)了先驗(yàn)信息的氣候區(qū)差異性對立木生物量預(yù)估的影響。氣候是環(huán)境因子最直接的外在表現(xiàn)，其綜合了降雨、光照、地形地貌等環(huán)境因子，決定了森林生物量和生產(chǎn)力的大小[26?27]。利用異速生物量方程進(jìn)行預(yù)估時，相同胸徑不同區(qū)域的立木生物量相差達(dá)到500倍[3,28]，相同區(qū)域內(nèi)的立木生物量相差也達(dá)到99倍[14,29]。Deans等[30]和Henry等[31]分別利用異速生物量方程估算熱帶森林生物量，發(fā)現(xiàn)同一氣候區(qū)內(nèi)立木生物量預(yù)估方程也存在差異，如何通過先驗(yàn)信息解釋預(yù)估方程的差異是下一步的研究重點(diǎn)。另外，氣候區(qū)之間的過渡帶森林的生物量和生產(chǎn)力均有增加的趨勢，進(jìn)一步驗(yàn)證先驗(yàn)信息區(qū)域差異對過渡區(qū)生物量預(yù)估的影響，對闡述氣候?qū)ι稚锪亢蜕a(chǎn)力的影響有重要的作用。

4.2 樹種生物學(xué)特性對異速生物量方程擬合的影響

大量研究從針葉林和闊葉林[3,13]、天然林和人工林[13,15]、常綠林和落葉林[14,32]等森林類型來比較立木生物量和預(yù)估方程的差異。樹種之間生物量的差異多數(shù)從木材密度來分析對比不同演替過程中木材密度的差異對立木生物量的影響[14,33]。本文利用櫟屬、樺木屬、楊屬、槭屬、桉屬和松屬六個屬的先驗(yàn)信息擬合異速生物量方程，從先驗(yàn)信息的生物學(xué)特性分析方程的預(yù)估效果。結(jié)果表明，槭屬平均決定系數(shù)R2為0.935，其他的平均決定系數(shù)R2均為0.934，樹種先驗(yàn)信息的差異對方程R2無顯著影響。平均偏差MB分別為? 0.205、?1.287、?2.926、?3.726、?0.367、?0.494；平均均方根差MRMSE分別為0.343、0.352、0.354、0.391、0.344、0.343，樹種先驗(yàn)信息的差異對MB和MRMSE影響顯著（P<0.01），與先驗(yàn)信息的氣候區(qū)差異的結(jié)果一致。

黃興召等[22]研究結(jié)果表明，a和b的值在針葉樹和闊葉樹之間有明顯差異，本研究結(jié)果顯示以松屬和櫟屬為先驗(yàn)信息的擬合效果>總數(shù)據(jù)集和桉屬>樺木屬>楊屬>槭屬。櫟屬、樺木屬、楊屬、槭屬、桉屬和松屬分別屬于殼斗科（Fagaceae）、樺木科（Betulaceae）、楊柳科（Salicaceae）、槭樹科（Aceraceae）、桃金娘科（Myrtaceae）和松科（Pinaceae）。但是，落葉松屬與松屬的生物學(xué)關(guān)系>落葉松屬與桉屬>落葉松屬與槭屬>落葉松屬與楊屬>落葉松屬與櫟屬>落葉松屬與樺木屬。由此推斷，樹種先驗(yàn)信息影響異速生物量的預(yù)估效果，與樹種之間的生物學(xué)關(guān)系不一致，其作用規(guī)律有待進(jìn)一步研究。本文選擇先驗(yàn)信息的氣候區(qū)差異和樹種差異對立木生物量預(yù)估的影響，驗(yàn)證了氣候因子和生物學(xué)特性對先驗(yàn)信息的影響和方程預(yù)估的作用，但先驗(yàn)信息差異性對交錯區(qū)立木生物量預(yù)估的影響及與樹種間生物學(xué)特性的關(guān)系，還有待進(jìn)一步研究。

5 結(jié)論

先驗(yàn)信息的區(qū)域差異和樹種差異對異速生物量方程的擬合效果R2影響不顯著，對方程的擬合效果MB和MRMSE影響顯著（P<0.01），先驗(yàn)信息的差異性影響立木地上部分生物量的預(yù)估效果。

先驗(yàn)信息區(qū)域差異對立木生物量預(yù)估的影響為總數(shù)據(jù)集和溫帶的MB和MRMSE優(yōu)于寒帶和熱帶的MB和MRMSE。擬合的立木生物量研究地點(diǎn)位于遼寧，屬于典型的溫帶季風(fēng)氣候，先驗(yàn)信息的區(qū)域作用效果與立木所在區(qū)域一致。

先驗(yàn)信息區(qū)域差異對立木生物量預(yù)估效果的影響程度為：松屬和櫟屬>總數(shù)據(jù)集和桉屬>樺木屬>楊屬>槭屬。但是，落葉松屬與松屬的生物學(xué)關(guān)系>落葉松屬與桉屬>落葉松屬與槭屬>落葉松屬與楊屬>落葉松屬與櫟屬>落葉松屬與樺木屬，先驗(yàn)信息生物學(xué)特性的預(yù)估效果表現(xiàn)規(guī)律與樹種的生物學(xué)關(guān)系不一致。