黃興召 豆玉萍 黃慶豐 陳東升
(1. 安徽農業大學林學與園林學院,合肥 230036;2. 中國林業科學研究院林業研究所,北京 100091)
立木地上生物量是森林生物量的重要組成部分,如何準確地估算立木地上生物量是研究森林生態系統物質循環的基礎[1?2]。異速生物量模型作為預估生物量的常用方法,已被廣泛應用于各類研究中[3?6]。異速生物量方程參數a和b隨著研究區域和樹種的變化而變化,通過收獲法,利用異速生物量方程預估生物量時,只有單獨建模才能提高預估精度[5?8]。貝葉斯方法作為其中的一種擬合方法,通過收集已建立方程的參數信息,將該信息轉化為參數的概率分布,以參數的均值和協方差矩陣為先驗信息進行擬合,越來越受到人們的關注[9?12]。但是,利用先驗信息對貝葉斯方法進行擬合時,研究者收集參數a和b的先驗信息具有隨機性和異質性,由此對模型擬合和預估效果產生的影響尚未得到驗證,也沒有闡述如何選取參數a和b的先驗信息來提升方程的擬合和預估效果。
收集已發表文獻異速生物量方程參數a和b的信息,是分析先驗信息異質性對立木生物量預估影響的基礎。Ter-Mikaelian和Korzukhin[13]建立了北美65個樹種的生物量預估方程、Zianis和Mencuccini[3]建立了279個異速生物量預估方程、Henry等[14]收集了850個立木及各組分的生物量預估方程、Zapata-Cuartas等[8]收集了134個的立木地上異速生物量方程、Rojas-García等[15]收集了487個立木及各組分的生物量預估方程。這些方程都是通過常規方法擬合得到參數a和b的值,為貝葉斯方法提供了豐富的先驗信息。分析發現,立木地上部分生物量預估方程參數b與緯度呈顯著相關,且在不同屬之間存在顯著差異[8,13]。研究者在利用貝葉斯方法預估立木地上部分生物量的過程中,討論了樣本大小對方程擬合效果和預估效果的影響,缺乏反映樹種特征的先驗信息的擬合比較[8,10,16]。樹木生長與生物量的積累受環境[17,18]和生物學特性[19?20]的影響,在方程擬合時先驗信息也會受到影響[21]。因此,檢驗參數a和b先驗信息的氣候區和生物學特性對立木生物量預估的影響顯得尤為重要。
通過分析已發表文獻關于立木地上部分異速生物量方程的參數a和b信息,將參數的先驗信息按主要氣候區和樹種進行分類,以日本落葉松(Larix kaempferi)立木地上部分生物量為擬合數據,利用貝葉斯方法擬合異速生物量方程,比較先驗信息的氣候差異和樹種差異對立木生物量擬合和預估的影響,探討先驗信息對立木生物量預估的作用規律,為利用貝葉斯方法擬合異速生物量方程時選擇先驗信息提供了理論和實踐的依據。
異速生物量方程作為預估立木生物量的常用模型,其方程表達式為:

其中,Mi為第i株樣本的地上部分生物量;Di為第i株樣本的胸徑;a和b為方程參數;ei為誤差項。
貝葉斯方法將公式(1)中的參數a、b定義為θ,需要知道θ的概率密度函數。接著,根據先驗信息確定參數a和b的先驗分布 π(θ),再結合樣本信息可得出后驗分布 π(θ|D)。 π(θ|D)的計算公式為:

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

表1 異速生物量方程先驗信息分類Table1Priorinformationclassificationofallometric biomassequation
以遼寧省清源縣大孤家林場(42°16?~42°23? N,124°48?~124°55? E)日本落葉松(Larixkaempeferi)立木地上部分生物量數據為擬合數據[22],該區屬于典型的溫帶季風氣候。基于樣地調查,在標準地附近選取標準木,利用收獲法測定樣木的樹干、樹皮、樹枝和樹葉各組分的生物量,求出立木地上部分生物量。本數據由幼齡林、中齡林、近熟林以及成熟林4個林齡段共45株立木的胸徑和地上部分生物量組成(表2)。本文在建模過程中,利用無重復抽樣,從45株總樣木中抽取30株樣木,進行1 000次無重復擬合,每次擬合未被抽取的15株作為驗證數據來檢驗預估效果。

表2 日本落葉松樣地信息Table2PlotsinformationofLarixkaempferiplantations
使用不同的先驗信息,利用貝葉斯方法擬合異速生物量方程時,采用決定系數(R2)、平均偏差(Mean Bias , MB)和平均均方根差(Mean Root Mean Square Error, MRMSE)來檢驗方程的擬合結果和預估效果。R2、MB和MRMSE的表達式分別為:

其中,yi為實際觀測值;yˉ為樣本平均值;y?i為模型預測值;n為樣本數。數據統計分析和繪圖使用R(Development Core Team, 2009)和Excel軟件。以不同區域和不同樹種參數a、b的均值以及協方差矩陣為先驗信息,應用R的MCMCglmm包[23]對異速生物量方程分別進行擬合。
為檢驗先驗信息區域差異對異速生物量方程擬合的影響,分別以寒帶、溫帶、熱帶以及總數據集參數a、b的均值及協方差矩陣為先驗信息,利用無重復抽樣,從45株總樣木中抽取30株樣木,進行1 000次無重復擬合(圖1)。寒帶和熱帶先驗信息擬合后參數a、b值的變化范圍相似;溫帶和總數據集先驗信息擬合后參數a、b值的變化范圍相似。以熱帶和寒帶為先驗信息進行1 000次無重復擬合,參數a的平均值均為?2.794,95%的置信區間分別為(?3.185, ?2.386)和(?3.197,?2.368);參數b的平均值均分別為2.640和2.638,95%的置信區間分別為(2.484, 2.788)和(2.476, 2.792)。以溫帶和總數據集為先驗信息進行1 000次無重復擬合,參數a的平均值分別為?2.836和?2.826,95%的置信區間分別為(?3.240,?2.418)和(?3.228, ?2.407);參數b的平均值均分別為2.655和2.651,95%的置信區間分別為(2.496, 2.809)和(2.491, 2.805)(表3)。這說明使用不同區域的先驗信息擬合異速生物量方程會影響參數a、b值的穩定性。

圖1 不同氣候區先驗信息預估參數a和b的擬合結果Fig.1Simulationresultsfortheestimationofthe parametersaandbusingdifferentprior informationofregion
決定系數R2反映方程的擬合效果(表3),寒帶、溫帶、熱帶以及總數據集四種先驗信息1 000次無重復擬合,方程的平均決定系數均為0.934。平均偏差MB和平均均方根差MRMSE反映了方程的預估效果,四種先驗信息1 000次無重復擬合,平均偏差分別為?1.392、?0.597、?0.929、?0.715,平均均方根差分別為0.347、0.344、0.347、0.345。由此可知,區域先驗信息主要影響方程的預估效果。

表3 不同氣候區先驗信息擬合異速生物量方程的結果對比Table3Resultcomparisonallometricbiomassequationsindifferentpriorinformationofregions
為了進一步比較寒帶、溫帶、熱帶以及總數據集四種先驗信息對方程預估效果的影響,對BE和RMSE分別進行方差分析(圖2);總數據集先驗信息的MB和MRMSE與溫帶的MB和MRMSE無顯著差異;寒帶與熱帶的MB和MRMSE也無顯著差異;其他先驗信息的MB和MRMSE之間存在顯著差異(P<0.01)。總數據集和溫帶先驗信息的擬合效果優于寒帶和熱帶,驗證數據集所在區域為溫帶,因此,利用貝葉斯方法擬合異速生物量方程時,應充分考慮先驗信息的氣候區特征。

圖2 不同氣候區先驗信息擬合的平均偏差和平均均方根差比較Fig.2Simulationresultsforthemeanbiasandmean rootmeansquareerrorusingdifferent priorinformationofregion
以櫟屬、樺木屬、楊屬、槭屬、桉屬、松屬和總數據集為先驗信息,在樣木數量為30的情況下,也進行1 000次無重復擬合(圖3)。桉屬、櫟屬、松屬以及總數據集的參數a的變化范圍相似,其平均值分別為?2.859、?2.866、?2.855、?2.826;樺木屬和楊屬的參數a的變化范圍相似,其平均值分別為?2.728和?2.712;槭屬的參數變化范圍尤為突出,其平均值為?2.492。參數b的變化范圍也表現出相同的規律,桉屬、櫟屬、松屬以及總數據集參數b的平均值分別為2.663、2.666、2.662、2.651;樺木屬和楊屬參數b的平均值分別為2.616和2.605;槭屬參數b的平均值為2.528。

圖3 不同屬先驗信息預估參數a和b的擬合結果Fig.3Simulationresultsfortheestimationofthe parametersaandbusingdifferent priorinformationofgenus
以櫟屬、樺木屬、楊屬、槭屬、桉屬、松屬和總數據集為先驗信息,進行異速生物量方程擬合(表4)表4。槭屬的平均決定系數為0.935,其他屬的平均決定系數均為0.934。櫟屬、樺木屬、楊屬、槭屬、桉屬、松屬和總數據集的平均偏差分別為?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。由此可知,樹種先驗信息主要影響方程的預估效果。

表4 不同屬先驗信息擬合異速生物量方程的結果對比Table4Resultcomparisonallometricbiomassequationsindifferentpriorinformationofgenus
為了進一步比較六個屬以及總數據集七種先驗信息對方程預估效果的影響,對MB和MRMSE分別進行方差分析(圖4)。櫟屬和桉屬、松屬和桉屬之間的MB無顯著差異,其他屬之間的MB均存在顯著差異(P<0.01),即櫟屬和桉屬為先驗信息的MB>松屬>總數>樺木屬>楊屬>槭屬。櫟屬和桉屬、松屬和桉屬、總數據集和桉屬之間的MRMSE無顯著差異,其他屬之間的MRMSE存在顯著差異(P<0.01),即松屬和櫟屬為先驗信息的MRMSE>總數據集和桉屬>樺木屬>楊屬>槭屬。因此,樹種的生物學特性也是影響生物量預估的重要因子。

圖4 不同屬先驗信息擬合的平均偏差和平均均方根差比較Fig.4Simulationresultsforthemeanbiasandmean rootmeansquareerrorusingdifferent priorinformationofgenus
利用貝葉斯方法擬合生物量方程,先驗信息尤為重要[8?9,16],通過收集立木地上部分異速生物量方程參數信息,分析表明參數a和b符合二元正態分布。擬合異速生物量方程過程中,貝葉斯原理使用先驗信息,而傳統方法只能使用立木數據。因此,貝葉斯方法可以實現少量的樣木高精度預估生物量,黃興召等[10]和Zell等[16]均證明,當樣木的數量小于50株時,貝葉斯方法的預估效果明顯優于傳統方法。立木地上部分異速生物量方程參數a和b的值有明顯的區域特征,參數b與緯度呈線性關系[10]。本文將參數a和b劃分為寒帶、溫帶和熱帶三個主要氣候區,基于信息量度的相似概率分布[22,24?25],以三個氣候區和總數據集參數a、b的均值和協方差矩陣為先驗信息,利用無重復抽樣,從總樣木中抽取30株立木生物量數據,進行1 000次無重復擬合,得出每次擬合結果中參數a、參數b、R2、MB和MRMSE的值。
將寒帶、溫帶和熱帶三個主要氣候區先驗信息的擬合結果與總數據集先驗信息的擬合結果進行比較,發現R2均為0.934無顯著差異,說明先驗信息的氣候差異對異速生物量方程的擬合效果無影響。MB和MRMSE差異分析表明,先驗信息的氣候差異對異速生物量方程的預估效果有顯著影響。本文中擬合數據的研究地點屬于典型的溫帶季風氣候,溫帶的先驗信息擬合效果(顯著)優于寒帶和熱帶,與總數據集無差異,驗證了先驗信息的氣候差異對立木生物量預估的影響,使用立木生長主要氣候區的先驗信息,有助于提升立木生物量的預估精度。
通過模型擬合,證實了先驗信息的氣候區差異性對立木生物量預估的影響。氣候是環境因子最直接的外在表現,其綜合了降雨、光照、地形地貌等環境因子,決定了森林生物量和生產力的大小[26?27]。利用異速生物量方程進行預估時,相同胸徑不同區域的立木生物量相差達到500倍[3,28],相同區域內的立木生物量相差也達到99倍[14,29]。Deans等[30]和Henry等[31]分別利用異速生物量方程估算熱帶森林生物量,發現同一氣候區內立木生物量預估方程也存在差異,如何通過先驗信息解釋預估方程的差異是下一步的研究重點。另外,氣候區之間的過渡帶森林的生物量和生產力均有增加的趨勢,進一步驗證先驗信息區域差異對過渡區生物量預估的影響,對闡述氣候對森林生物量和生產力的影響有重要的作用。
大量研究從針葉林和闊葉林[3,13]、天然林和人工林[13,15]、常綠林和落葉林[14,32]等森林類型來比較立木生物量和預估方程的差異。樹種之間生物量的差異多數從木材密度來分析對比不同演替過程中木材密度的差異對立木生物量的影響[14,33]。本文利用櫟屬、樺木屬、楊屬、槭屬、桉屬和松屬六個屬的先驗信息擬合異速生物量方程,從先驗信息的生物學特性分析方程的預估效果。結果表明,槭屬平均決定系數R2為0.935,其他的平均決定系數R2均為0.934,樹種先驗信息的差異對方程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,樹種先驗信息的差異對MB和MRMSE影響顯著(P<0.01),與先驗信息的氣候區差異的結果一致。
黃興召等[22]研究結果表明,a和b的值在針葉樹和闊葉樹之間有明顯差異,本研究結果顯示以松屬和櫟屬為先驗信息的擬合效果>總數據集和桉屬>樺木屬>楊屬>槭屬。櫟屬、樺木屬、楊屬、槭屬、桉屬和松屬分別屬于殼斗科(Fagaceae)、樺木科(Betulaceae)、楊柳科(Salicaceae)、槭樹科(Aceraceae)、桃金娘科(Myrtaceae)和松科(Pinaceae)。但是,落葉松屬與松屬的生物學關系>落葉松屬與桉屬>落葉松屬與槭屬>落葉松屬與楊屬>落葉松屬與櫟屬>落葉松屬與樺木屬。由此推斷,樹種先驗信息影響異速生物量的預估效果,與樹種之間的生物學關系不一致,其作用規律有待進一步研究。本文選擇先驗信息的氣候區差異和樹種差異對立木生物量預估的影響,驗證了氣候因子和生物學特性對先驗信息的影響和方程預估的作用,但先驗信息差異性對交錯區立木生物量預估的影響及與樹種間生物學特性的關系,還有待進一步研究。
先驗信息的區域差異和樹種差異對異速生物量方程的擬合效果R2影響不顯著,對方程的擬合效果MB和MRMSE影響顯著(P<0.01),先驗信息的差異性影響立木地上部分生物量的預估效果。
先驗信息區域差異對立木生物量預估的影響為總數據集和溫帶的MB和MRMSE優于寒帶和熱帶的MB和MRMSE。擬合的立木生物量研究地點位于遼寧,屬于典型的溫帶季風氣候,先驗信息的區域作用效果與立木所在區域一致。
先驗信息區域差異對立木生物量預估效果的影響程度為:松屬和櫟屬>總數據集和桉屬>樺木屬>楊屬>槭屬。但是,落葉松屬與松屬的生物學關系>落葉松屬與桉屬>落葉松屬與槭屬>落葉松屬與楊屬>落葉松屬與櫟屬>落葉松屬與樺木屬,先驗信息生物學特性的預估效果表現規律與樹種的生物學關系不一致。