基于層次貝葉斯法的高山松生物量估測模型

黃金君， 舒清態(tài)， 席 磊， 孫 楊， 劉玥伶

( 1.廣西壯族自治區(qū)中國科學院廣西植物研究所，廣西桂林541006；2.西南林業(yè)大學林學院，云南昆明650224)

森林地上生物量(AGB)作為森林生態(tài)系統(tǒng)能量交匯以及生態(tài)健康功能評價的重要指標，對全球陸地與大氣間的碳物質(zhì)平衡影響巨大。根據(jù)《2019年中國國土綠化狀況公報》，中國森林覆蓋率已達到22.96%(2019年6月17日)，相比于全國第8次森林資源清查的覆蓋率提高了1.33個百分點。加強對喬木林不同器官生物量的研究，獲得更加精準的森林地上生物量數(shù)據(jù)，對于深入研究森林生態(tài)系統(tǒng)生產(chǎn)力以及揭示陸地生態(tài)系統(tǒng)結(jié)構(gòu)與生態(tài)環(huán)境的關(guān)系具有重要的科學與現(xiàn)實價值[1-2]。當前國內(nèi)外估測森林地上生物量的方法主要基于森林資源清查法與基于遙感估測法，基于森林資源清查法是指利用高密度樣地調(diào)查獲得的數(shù)據(jù)建立不同樹種異速生物量模型，進而估算不同尺度的森林地上生物量[3]。異速生物量模型通過將易測的胸徑或樹高等指標相結(jié)合以提供一種高效簡單的估算生物量方法，目前已被廣泛應(yīng)用[4-5]。然而使用傳統(tǒng)方法擬合異速生物量模型存在過多局限，其只將模型參數(shù)視作固定變量，未能綜合考慮存在的多種生物或非生物因素[6-7]。因此選用層次貝葉斯方法有助于解決這一問題，其利用概率分布描述模型中的未知變量，并將樣本信息與參量視作隨機變量[8]。在當前精準提升森林質(zhì)量的背景下，貝葉斯方法正逐漸成為林業(yè)領(lǐng)域中學者們研究的熱點。Green等[9]于1992年首次在林業(yè)應(yīng)用中引入層次貝葉斯方法，并于1996年建立了基于種植園杉木的貝葉斯模型[10]。Clark等[11]于2007年使用貝葉斯方法建立了樹木直徑生長模型，Metcalf等[12]于2009年建立了基于非參數(shù)貝葉斯法的樹木死亡率模型，張雄清等[13]于2014年建立了基于貝葉斯法的杉木人工林樹高增長模型，王冬至等[14]于2019年建立了基于貝葉斯法的針闊混交林樹高與胸徑混合效應(yīng)模型，姚丹丹等[15]于2020年建立了基于貝葉斯法的蒙古櫟林單木樹高-胸徑模型。已有研究結(jié)果表明，利用貝葉斯方法可以有效提高模型估測精度與可靠程度。層次貝葉斯法基于貝葉斯理論發(fā)展，通過分離不同層面中多個參數(shù)的復雜關(guān)系，從而充分發(fā)揮具有層次結(jié)構(gòu)模型的優(yōu)勢。本研究以高山松天然林為研究對象，以區(qū)域為分層單位，分別采用層次貝葉斯法與非層次貝葉斯法構(gòu)建高精度單木及不同組分生物量估測模型，探討貝葉斯方法在單木及不同組分生物量估測模型中的運用。

1 材料與方法

1.1 研究區(qū)概況

以云南省香格里拉市小中甸鎮(zhèn)和洛吉鄉(xiāng)的高山松天然林為研究對象。小中甸鎮(zhèn)地理坐標為99°36′～99°59′E，27°20′～27°43′N，地處香格里拉市南部高寒壩區(qū)，多年平均氣溫5.8 ℃，森林覆蓋度69.42%，主要喬木樹種包括高山松(PinusdensataMast)、云杉(PiceaasperataMast)、冷杉[Abiesfabri(Mast.)Craib]、落葉松[LarixgmelinⅡ (Rupr.) Kuzen]等。洛吉鄉(xiāng)地理坐標為99°55′～100°19′E，27°38′～28°06′N，地處香格里拉市東部，海拔差異懸殊大，南部河谷區(qū)最低點海拔為1 503 m，北部高寒山區(qū)最高點海拔為4 495 m，南部河谷區(qū)年平均氣溫為13～15 ℃，北部高寒山區(qū)年平均氣溫為5.5 ℃，森林覆蓋度為78%，主要喬木樹種有高山松(PinusdensataMast)、云南松(PinusyunnanensisFranch)、冷杉[Abiesfabri(Mast.)Craib]、云杉(PiceaasperataMast)、落葉松[LarixgmelinⅡ (Rupr.) Kuzen]、核桃(JuglansregiaL.)等。將在小中甸鎮(zhèn)和洛吉鄉(xiāng)采集到的數(shù)據(jù)分為區(qū)域Ⅰ和區(qū)域Ⅱ，基本數(shù)據(jù)信息如表1所示。2個區(qū)域高山松的胸徑和樹高數(shù)據(jù)包含區(qū)域間與樣地間的嵌套關(guān)系，即分層數(shù)據(jù)[16]。數(shù)據(jù)分為2層，第1層是地域空間的胸徑和樹高數(shù)據(jù)，第2層為每個地域空間內(nèi)樣本的胸徑和樹高數(shù)據(jù)。通過對研究區(qū)樣本數(shù)據(jù)進行方差分析可知，單木樣本數(shù)據(jù)之間存在顯著的地域性差異。

表1 高山松基本數(shù)據(jù)信息

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

分別在2個區(qū)域內(nèi)同一森林類型的林分中設(shè)置56塊面積為30 m×30 m的高山松天然純林標準樣地，記錄樣地的中心點經(jīng)緯度坐標、海拔、坡度和坡向等空間位置信息，郁閉度、林齡和齡級等測樹學因子。最后于每塊樣地中隨機選擇1～3株具有代表性的高山松單木，伐倒后測定其生物學特性。伐倒樣木共計115株，包括區(qū)域Ⅰ的60株樣木和區(qū)域Ⅱ的55株樣木，每株高山松的數(shù)據(jù)被分成單木生物量、樹干生物量、木材生物量、樹皮生物量、樹枝生物量、樹冠生物量和樹葉生物量7個組分來分別進行測定。樹干生物量主要采用材積密度法測量，枝、葉生物量采用分級標準枝法測量。最后分別取樣稱量，烘干處理后計算單木生物量及各器官的生物量(表2)。

表2 單木及不同組分的高山松生物量數(shù)據(jù)

1.3 異速生物量模型的建立

異速生長方程反映了樹木各維量間按比例變化協(xié)調(diào)增長的關(guān)系，具有廣泛的適應(yīng)性、較高的靈活性和良好的擬合性等特性，因此該方程在森林生物量估算中被廣泛應(yīng)用[17-18]。模型表達式為：

G=aDbHc

(1)

式中：G代表高山松單木及不同組分的生物量，D代表胸徑，H代表樹高，a、b、c代表高山松單木及不同組分生物量模型的異速生長參數(shù)。

1.4 研究方法

1.4.1 非層次貝葉斯法 隨著馬爾科夫鏈蒙特卡洛(MCMC)方法的逐漸成熟，貝葉斯方法在多個領(lǐng)域的研究與應(yīng)用越來越廣泛[19-20]。與經(jīng)典統(tǒng)計學方法相比，貝葉斯方法包含1個最為重要的要素(參數(shù)的先驗信息)，該方法通過在已有知識的基礎(chǔ)上推測出新的信息，即得到所求模型參數(shù)的后驗分布。經(jīng)典統(tǒng)計學方法將未知參數(shù)視作固定變量，只考慮樣本與整體信息，而貝葉斯方法綜合考慮樣本信息、整體信息及先驗信息，并使用概率分布描述模型的未知參數(shù)。在未進行樣本數(shù)據(jù)抽樣時，研究人員根據(jù)已有的文獻資料或經(jīng)典統(tǒng)計學方法得到模型參數(shù)?的信息(如均值與方差)，且設(shè)定了1個概率分布R(?)來描述該信息的隨機性，此分布則成為參數(shù)?的先驗分布。將本研究的高山松樣本數(shù)據(jù)設(shè)為M=(M1，…，Mn)，似然函數(shù)為L(M|?)，后驗分布R(?|M)則是建立在先驗分布R(?)與似然函數(shù)的基礎(chǔ)上，樣本數(shù)據(jù)的邊緣概率密度分布為?L(M|?)R(?)d?，最后由貝葉斯公式對先驗分布、似然函數(shù)以及后驗分布進行組合，計算公式如下：

(2)

1.4.2 層次貝葉斯法 層次貝葉斯法作為現(xiàn)代貝葉斯方法中極為典型的代表，通過對不同層次上多個參數(shù)間的復雜關(guān)系進行分離，解決了普通貝葉斯法對于復雜參數(shù)分布估計的困難[21]。與普通貝葉斯法相比，層次貝葉斯法允許不同的參數(shù)出現(xiàn)在不同層面上，且采用了分層先驗信息，有助于消除先驗信息對模型估計結(jié)果的過度影響，分層先驗的表達公式如下：

R(?)=?R1(?|τ)R2(τ)dτ

(3)

首先對未知參數(shù)?設(shè)定1個已知先驗分布?～R1(?|τ)，當給定的先驗分布中超參數(shù)τ難以確定時，再給出1個超先驗R2(τ)，因此由先驗與超先驗共同構(gòu)成新的先驗，稱為分層先驗。由參數(shù)?和超參數(shù)τ組成的后驗分布函數(shù)為：

R(?,τ|M)∝R(M|?,τ)R1(?|τ)R2(τ)

(4)

在本研究中不僅采用分層先驗信息，同時還對數(shù)據(jù)進行了分層。假設(shè)樣本數(shù)據(jù)被分成K組，即把單木及不同組分生物量的數(shù)據(jù)分成了I區(qū)和Ⅱ區(qū)，樣本觀測值Mij(i=1，2，3，…，K；j=1，2，3，…，Ni)代表第i組中第j個樣本觀測值，I區(qū)和Ⅱ區(qū)之間樣本觀測值差異較大，每個組內(nèi)的樣本觀測值差異較小。

1.4.3 先驗信息 貝葉斯方法中豐富的先驗信息影響著模型參數(shù)的估計，需要為所求參數(shù)構(gòu)造適當?shù)南闰灧植肌Ｓ行┭芯咳藛T建議采用無信息先驗，無信息先驗通常出現(xiàn)在方差無窮大與均值為0的高斯分布中，對參數(shù)不會造成太大的影響。也可采用有信息先驗(Informative prior)作為貝葉斯方法的先驗分布信息，其來源于統(tǒng)計學方法計算得到的數(shù)據(jù)或前人的研究成果等[22]。本研究采用非線性最小二乘法的估計結(jié)果作為非層次貝葉斯法與層次貝葉斯法的先驗信息[23]。

1.5 軟件處理

使用OpenBUGS軟件[24]實現(xiàn)層次貝葉斯法與非層次貝葉斯法的模型參數(shù)估計。OpenBUGS軟件以MCMC方法為基礎(chǔ)，通過Gibbs和Metropolis算法從參數(shù)的后驗分布中抽樣，從而完成馬爾科夫鏈的迭代以及參數(shù)的估計，極大促進了貝葉斯方法的應(yīng)用。

1.6 模型精度評價方法

采用十折交叉驗證方法，將高山松樣本數(shù)據(jù)集分成10份，輪流將數(shù)據(jù)集中的9份作為訓練數(shù)據(jù)，剩下的1份作為測試數(shù)據(jù)。通過選用決定系數(shù)(R2)、均方根誤差(RMSE)和平均絕對誤差(MAD)3個指標，對層次貝葉斯法與非層次貝葉斯法擬合不同組分異速生物量模型的精度進行評價，每個指標的計算公式如下：

(5)

(6)

(7)

2 結(jié)果與分析

2.1 基于非層次貝葉斯法的參數(shù)估計結(jié)果分析

采用非層次貝葉斯法對高山松單木及不同組分異速生物量模型[公式(1)]進行擬合，發(fā)現(xiàn)模型R2均較高(樹葉生物量模型除外)且存在極顯著差異(P<0.01)，表明該方法適合該方程的擬合。為獲得準確的異速生長后驗參數(shù)估計值，本研究最開始設(shè)置1 000次抽樣以消除初始值對抽樣的影響，待模型參數(shù)擬合穩(wěn)定后迭代次數(shù)設(shè)為10 000次。迭代完成后便可得到非層次貝葉斯法中各個器官異速生長參數(shù)的均值、標準差和95%置信區(qū)間。如表3所示，異速生長參數(shù)a的標準差和95%置信區(qū)間范圍<異速生長參數(shù)b的標準差和95%置信區(qū)間范圍<異速生長參數(shù)c的標準差和95%置信區(qū)間范圍，說明異速生長參數(shù)a的估計值在模型抽樣迭代過程中比參數(shù)b和參數(shù)c更穩(wěn)定。

表3 基于非層次貝葉斯法的單木及不同組分生物量模型參數(shù)估計結(jié)果

2.2 基于層次貝葉斯法的參數(shù)估計結(jié)果分析

使用非線性最小二乘法的估計結(jié)果作為先驗信息，以區(qū)域進行分層，采用層次貝葉斯法對高山松單木及不同組分異速生物量模型[公式(1)]進行擬合，發(fā)現(xiàn)模型R2均較高(樹葉生物量模型除外)且存在極顯著差異(P<0.01)，表明該方法適合該方程的擬合。為獲得最優(yōu)的模型參數(shù)后驗估計值，本研究共設(shè)置11 000次迭代，得到層次貝葉斯法I區(qū)和Ⅱ區(qū)中各個器官異速生長參數(shù)的均值、標準差和95%置信區(qū)間。層次貝葉斯法的參數(shù)估計結(jié)果如表4所示，可見Ⅰ區(qū)和Ⅱ區(qū)中各個器官異速生物量模型均有不同的異速生長參數(shù)，且參數(shù)的均值差異較大，說明不同區(qū)域中各個器官的生物量存在差異。異速生長參數(shù)a的標準差和95%置信區(qū)間范圍整體上小于異速生長參數(shù)b和c，與非層次貝葉斯法的估計結(jié)果一致。

表4 基于層次貝葉斯法的單木及不同組分生物量模型參數(shù)估計結(jié)果

2.3 不同方法擬合高山松異速生物量模型的精度對比分析

用層次貝葉斯法與非層次貝葉斯法擬合高山松不同組分異速生物量模型的結(jié)果(表5)顯示，這2種方法的決定系數(shù)(R2)、均方根誤差(RMSE)和平均絕對誤差(MAD)呈現(xiàn)的規(guī)律一致：樹干生物量、木材生物量和單木生物量模型效果最優(yōu)，樹皮生物量、樹冠生物量和樹枝生物量模型效果較優(yōu)，樹葉生物量模型效果較差，層次貝葉斯法擬合模型的R2區(qū)間為[0.365 0,0.965 0]，非層次貝葉斯法擬合模型的R2區(qū)間為[0.437 0,0.964 7]。通過對比層次貝葉斯法與非層次貝葉斯法擬合模型的效果發(fā)現(xiàn)，層次貝葉斯法擬合模型的效果更好(樹枝與樹葉生物量模型除外)，用該法擬合的單木生物量、樹干生物量、木材生物量、樹皮生物量和樹冠生物量模型的R2分別較非層次貝葉斯法提高了0.012 0、0.000 3、0.000 1、0.006 5、0.005 6，RMSE分別降低了8.94 kg、0.28 kg、0.03 kg、0.30 kg、0.27 kg，MAD分別降低了3.31 kg、0.11 kg、0.42 kg、0.03 kg、0.30 kg。

表5 基于單木及不同組分異速生物量模型的不同方法精度結(jié)果對比

3 討論與結(jié)論

本研究主要探討了貝葉斯方法在異速生物量模型中的應(yīng)用，選取香格里拉市Ⅰ區(qū)和Ⅱ區(qū)共115株高山松單木進行解析，每株高山松的數(shù)據(jù)被分成單木生物量、樹干生物量、木材生物量、樹皮生物量、樹冠生物量、樹枝生物量和樹葉生物量7個組分來進行測定，重點分析了層次貝葉斯法與非層次貝葉斯法擬合不同組分異速生物量模型中存在的差異，證實不同區(qū)域?qū)ι锪慨a(chǎn)生了重要的影響。

貝葉斯方法的一個優(yōu)勢在于使用馬爾科夫鏈判斷模型參數(shù)的平穩(wěn)分布與收斂，而經(jīng)典統(tǒng)計學方法使用樣本的漸進方差來判定，無法考慮樣本數(shù)據(jù)中包含的不確定信息。在使用貝葉斯方法擬合不同組分異速生物量模型的過程中，除了樹枝與樹葉生物量模型，層次貝葉斯法的擬合效果整體優(yōu)于非層次貝葉斯法。該結(jié)果一方面說明了層次貝葉斯法有助于提高生物量模型的估測精度與準確度，另一方面說明，樹枝與樹葉生物量占總樹生物量的比例較小，所含有的不確定性因素較少[25]，因此使用非層次貝葉斯法效果更好。層次貝葉斯法適合分析包含多種不確定性因素的復雜數(shù)據(jù)，通過綜合考慮多層面的各種因素以及將模型中的參數(shù)看作隨機變量，來展現(xiàn)不同組間樣本數(shù)據(jù)的異質(zhì)性與同一組內(nèi)不同樣本數(shù)據(jù)的相異性，從而提高建模的靈活性和準確性。姚丹丹等[26]使用層次貝葉斯法、貝葉斯法和非線性最小二乘法擬合蒙古櫟林單木枯死模型，發(fā)現(xiàn)層次貝葉斯方法的模型結(jié)果最好。黃興召等[27]使用層次貝葉斯法擬合安徽省馬鬃嶺林場和福建省東安林場杉木的生物量轉(zhuǎn)換因子函數(shù)時，發(fā)現(xiàn)設(shè)置區(qū)域為隨機效應(yīng)有助于提高模型精度，層次貝葉斯法可以解決區(qū)域?qū)τ谀骋活愋蜕稚锪康挠绊憽?/p>

本研究采用層次貝葉斯法與非層次貝葉斯法擬合了高山松單木及不同組分的異速生物量模型，提供了1種在區(qū)域?qū)用嫔瞎烙嫺呱剿缮锪康挠行Х椒āＭㄟ^分析2種方法的精度發(fā)現(xiàn)，層次貝葉斯法擬合的單木生物量、樹干生物量、木材生物量、樹皮生物量和樹冠生物量模型的R2分別較非層次貝葉斯法提高了0.012 0、0.000 3、0.000 1、0.006 5、0.005 6，RMSE分別降低了8.94 kg、0.28 kg、0.03 kg、0.30 kg、0.27 kg，MAD分別降低了3.31 kg、0.11 kg、0.42 kg、0.03 kg、0.30 kg。通過對比2種方法對不同組分生物量模型擬合的效果發(fā)現(xiàn)，樹干生物量、木材生物量和單木生物量模型的效果最優(yōu)，樹皮生物量、樹冠生物量和樹枝生物量模型的效果較優(yōu)，樹葉生物量模型的效果較差。綜上所述，層次貝葉斯法可以有效地提高模型估測精度，且更適合分析在復雜環(huán)境中受多因素影響的數(shù)據(jù)，加強層次貝葉斯法在林業(yè)領(lǐng)域中的研究，將具有較大的參考意義。