林分徑級分布預(yù)測模型參數(shù)的貝葉斯法估計(jì)1)
——以金溝嶺林場樣地?cái)?shù)據(jù)為例

于汝川 張青 岳朝方 亢新剛

(北京林業(yè)大學(xué)，北京，100083) (德國巴登—符騰堡州林業(yè)科學(xué)研究院) (北京林業(yè)大學(xué))

?

林分徑級分布預(yù)測模型參數(shù)的貝葉斯法估計(jì)1)
——以金溝嶺林場樣地?cái)?shù)據(jù)為例

于汝川 張青 岳朝方 亢新剛

(北京林業(yè)大學(xué)，北京，100083) (德國巴登—符騰堡州林業(yè)科學(xué)研究院) (北京林業(yè)大學(xué))

以吉林省金溝嶺林場云冷杉針闊混交異齡林26塊檢查法樣地的5次觀測數(shù)據(jù)為基礎(chǔ)，建立轉(zhuǎn)移矩陣模型對一定周期的林分徑級分布進(jìn)行預(yù)測。分別利用經(jīng)典統(tǒng)計(jì)學(xué)方法和貝葉斯方法對轉(zhuǎn)移矩陣模型的參數(shù)進(jìn)行估計(jì)，建立了固定參數(shù)的矩陣模型和貝葉斯矩陣模型，并對兩種模型的預(yù)測效果進(jìn)行對比。結(jié)果顯示，固定參數(shù)的轉(zhuǎn)移矩陣模型對林分徑級分布的預(yù)測值比實(shí)際值偏高，貝葉斯模型的預(yù)測結(jié)果更接近林分的實(shí)際徑級分布，證明了貝葉斯參數(shù)估計(jì)方法比固定參數(shù)平均的方法所建模型具有更高的預(yù)測精度。

轉(zhuǎn)移矩陣模型；林分徑級分布；貝葉斯統(tǒng)計(jì)；MCMC方法；Gibbs抽樣算法

貝葉斯法是基于貝葉斯定理而發(fā)展起來用于系統(tǒng)地闡述和解決統(tǒng)計(jì)問題的方法[1]。與經(jīng)典統(tǒng)計(jì)學(xué)相比，貝葉斯統(tǒng)計(jì)具有綜合利用先驗(yàn)信息和樣本信息、把樣本和參數(shù)都看作隨機(jī)變量、對參數(shù)或模型的構(gòu)造要求少、不要求樣本服從正態(tài)分布等優(yōu)點(diǎn)。許多學(xué)者對貝葉斯統(tǒng)計(jì)在林業(yè)上的應(yīng)用進(jìn)行了探索和研究，早在20世紀(jì)80年代，賈乃光[2]就曾指出，貝葉斯統(tǒng)計(jì)學(xué)在林業(yè)中的應(yīng)用前景是完全可以肯定的。張雄清等[1，3]利用貝葉斯方法估計(jì)杉木人工林林分?jǐn)嗝娣e生長模型及樹高生長模型的參數(shù)，取得了較好的建模效果，姚丹丹等[4]利用貝葉斯方法研究長白山落葉松林分優(yōu)勢高生長模型，證明了貝葉斯方法具有一定的優(yōu)越性。在林業(yè)研究中，每隔一定年限對樣地進(jìn)行的重復(fù)調(diào)查數(shù)據(jù)使得林業(yè)研究往往具有充足的歷史資料，這些資料所包含的信息是進(jìn)行統(tǒng)計(jì)預(yù)測不可缺少的一個(gè)因素。這些信息在貝葉斯統(tǒng)計(jì)中可作為先驗(yàn)信息得到充分的利用，每一次調(diào)查數(shù)據(jù)的分析結(jié)果都是下一次調(diào)查分析最合理的先驗(yàn)信息。因此，貝葉斯模型在林業(yè)數(shù)據(jù)的建模研究中具有一定的優(yōu)勢。

徑級分布模型是研究林分結(jié)構(gòu)較常用一種模型。將林分按照樹木的胸徑大小統(tǒng)計(jì)各徑級的立木數(shù)量即可得到徑級分布模型。該模型有效簡化了復(fù)雜的林分結(jié)構(gòu)，為許多林業(yè)經(jīng)營問題的研究提供了便利和基礎(chǔ)。一些關(guān)于天然異齡林混交林經(jīng)營的模型(如均衡曲線模型)都是以林分的徑級分布模型為基礎(chǔ)發(fā)展而來的[5]。

轉(zhuǎn)移概率矩陣模型[6-8]是用來預(yù)測林分徑級分布的一種模型。該模型可對被徑級分布模型簡化后的林分內(nèi)部結(jié)構(gòu)進(jìn)行建模研究，能反映出林分徑級之間的發(fā)展變化關(guān)系。轉(zhuǎn)移矩陣模型中有許多參數(shù)(即轉(zhuǎn)移概率)需要進(jìn)行準(zhǔn)確的估計(jì)，許多學(xué)者給出了不同的估計(jì)方法，Solomon et al.[9]提出兩級預(yù)測模型，把矩陣模型的參數(shù)作為林分指標(biāo)的函數(shù)用線性回歸進(jìn)行估測。由于林分參數(shù)之間并非線性關(guān)系，Kl?dtke和Yue[10-11]利用非線性實(shí)時(shí)參數(shù)模型模擬矩陣參數(shù)與林分密度和采伐強(qiáng)度的動(dòng)態(tài)關(guān)系。以上方法都是基于經(jīng)典統(tǒng)計(jì)學(xué)方法進(jìn)行的參數(shù)估計(jì)，筆者采用了貝葉斯方法對轉(zhuǎn)移矩陣模型的參數(shù)進(jìn)行估計(jì)研究，并與一些經(jīng)典統(tǒng)計(jì)方法進(jìn)行對比，以檢驗(yàn)貝葉斯方法在林業(yè)模型研究中的優(yōu)點(diǎn)和不足。

1 研究區(qū)概況

研究地區(qū)位于吉林省金溝嶺林場(43°22′N，130°10′E)。地貌屬低山丘陵，海拔300～1 200 m，坡度5°～25°。該地區(qū)屬季風(fēng)型氣候，全年平均氣溫3.9 ℃，年降水量600～700 mm，年生長期為120 d。研究區(qū)林分主要為天然針闊混交林，植被屬長白山植物區(qū)系，立地條件較好，植物種類繁多。主要樹種有云杉(Piceakoraiensis)、冷杉(Abiesholophylla)和紅松(Pinuskoraiensis)，其他樹種有：色木(Acermono)、水曲柳(Fraxinusmandschurica)、胡桃楸(Juglansmandshurica)、黃波欏(Phellodendronamurense)、白樺(Betulaplatyphylla)等[12-13]。

本研究采用金溝嶺林場第一大區(qū)第5小區(qū)的檢查法數(shù)據(jù)，共有26塊樣地，每塊樣地面積0.04 km2，總面積1.04 km2。如表1所示，樣本包含了從1993—2013年間共5次的調(diào)查觀測數(shù)據(jù)，相鄰兩次調(diào)查的時(shí)間間隔為5 a。因各年的調(diào)查月份不同，考慮到研究區(qū)生長季的問題，本研究按照Huang等[14-15]提出的北方森林生長季法，將各調(diào)查年份的數(shù)據(jù)統(tǒng)一轉(zhuǎn)換到當(dāng)年5月份。

表1 立木數(shù)量及胸徑統(tǒng)計(jì)

2 研究方法

2.1 模型設(shè)定

設(shè)定林分的起始徑級為4 cm，級差為4 cm，徑級數(shù)設(shè)定k=15，建立徑級分布模型。圖1展示了該研究區(qū)每公頃樣地各年份的林分徑級分布?？梢钥闯?，該林分各年的徑級分布大致服從反“J”型的指數(shù)分布。

2.2 轉(zhuǎn)移概率矩陣模型

對林分徑級分布的預(yù)測可以采用轉(zhuǎn)移概率矩陣模型，該模型描述如下：

設(shè)ω表示預(yù)測周期，k為林分的徑級數(shù)，Nt表示林分t時(shí)刻的觀測徑級分布，A表示轉(zhuǎn)移概率矩陣，M表示枯損率，I表示進(jìn)界率，則間隔ω時(shí)間預(yù)測的徑階分布Nt+ω為

Nt+ω=ANt-MNt+INt。

(1)

用矩陣表示為

(2)

式中：ak為由徑級k到k+1的進(jìn)階生長率；bk為徑級k中留余立木數(shù)量的比例；ck為由徑級k到k+2的進(jìn)階生長率；mk為徑級k的立木枯損率；ik為徑級k的進(jìn)界率；nk,t為t時(shí)刻徑級k的立木數(shù)量。該矩陣模型的參數(shù)可根據(jù)實(shí)際樣本數(shù)據(jù)進(jìn)行估計(jì)，參數(shù)一旦確定，其在時(shí)間上就保持不變。

圖1 不同年份的林分徑級分布

2.3 模型參數(shù)的估計(jì)方法

2.3.1 固定參數(shù)模型

固定參數(shù)模型的參數(shù)可直接利用林分在一個(gè)預(yù)測周期前后的兩次觀測數(shù)據(jù)直接求解和估計(jì)：

(3)

2.3.2 貝葉斯模型

貝葉斯統(tǒng)計(jì)思想認(rèn)為模型的所有參數(shù)及預(yù)測值均可看作服從某種分布的隨機(jī)變量。設(shè)θ為模型的參數(shù)向量，則模型預(yù)測值y的密度函數(shù)可記為f(y|θ)，其中f為模型參數(shù)與預(yù)測值的連接函數(shù)。給定參數(shù)的先驗(yàn)密度函數(shù)為g(θ)，若觀測到某個(gè)預(yù)測值樣本后，則參數(shù)θ的后驗(yàn)密度g(θ|樣本)∝g(θ)f(樣本|θ)，貝葉斯模型參數(shù)的估計(jì)可根據(jù)其后驗(yàn)密度函數(shù)得到。

對于轉(zhuǎn)移概率矩陣模型，根據(jù)貝葉斯模型參數(shù)估計(jì)方法，可建立貝葉斯模型如下：

在一個(gè)預(yù)測周期ω過后，林分各徑級的立木數(shù)量

(4)

其中，

μ1=n1,t(b1-m1+i1)；μ2=n2,t(b2-m2+i2)+n1,ta1；

μp=np,t(bp-mp+ip)+np-2,tcp-2+np-1,tap-1,

p=3、…、k-2；

μk-1=nk-3,tck-3+nk-2,tak-2,μk=nk-2,tck-2。

(5)

用向量和矩陣表示為：

(μ1,μ2,…,μk)=μ=ANt-MNt+INt。

(6)

在轉(zhuǎn)移概率矩陣A中，設(shè)定參數(shù)服從正態(tài)分布

(7)

在向量M和I中，設(shè)定參數(shù)服從正態(tài)分布

(8)

由于模型參數(shù)均為概率值，因此將其正態(tài)分布均值的先驗(yàn)密度設(shè)定為0到1上的均勻分布，即

(9)

上述所有正態(tài)分布的隨機(jī)方差均從反伽馬分布中進(jìn)行抽樣，即

(10)

2.3.3 貝葉斯模型求解

2.4 模型評價(jià)

采用均方誤差MSE、均方根誤差RMSE以及平均絕對百分誤差MAPE來衡量模型的預(yù)測效果，其計(jì)算公式如下：

(11)

3 結(jié)果與分析

由于原始樣本共包含5次觀測年份的數(shù)據(jù)，建模過程中利用前4次(1993年、1998年、2003年和2008年)觀測數(shù)據(jù)作為訓(xùn)練樣本來估計(jì)模型的參數(shù)，利用第5次(2013年)觀測數(shù)據(jù)作為測試樣本來測試和評價(jià)所建立模型的預(yù)測精度和效果。模型的預(yù)測周期ω設(shè)定為5 a。

3.1 固定參數(shù)模型的參數(shù)估計(jì)

利用R語言編程計(jì)算固定參數(shù)模型的參數(shù)，訓(xùn)練樣本中的相鄰兩次測量可估計(jì)出一組參數(shù)值，4次測量共可計(jì)算出3組參數(shù)估計(jì)值。為充分利用訓(xùn)練樣本中所包含的信息，取3組參數(shù)估計(jì)值的平均值作為最終的固定參數(shù)模型的參數(shù)估計(jì)，結(jié)果如表2所示。

3.2 貝葉斯模型的參數(shù)估計(jì)

貝葉斯模型的參數(shù)是根據(jù)某一觀測值出現(xiàn)后的后驗(yàn)概率密度進(jìn)行估計(jì)的，因此該模型的參數(shù)估計(jì)會(huì)隨著每次新觀測數(shù)據(jù)的出現(xiàn)而進(jìn)行更新調(diào)整，這樣調(diào)整后的后驗(yàn)密度信息中不但包含先驗(yàn)信息，而且包含了樣本中的全部信息，因此模型參數(shù)估計(jì)的精度會(huì)隨著樣本信息的不斷更新而逐漸提高。

首先設(shè)定模型參數(shù)的先驗(yàn)信息，即各參數(shù)服從正態(tài)分布，且分布的均值為0到1上均勻分布的隨機(jī)變量；然后結(jié)合1993年和1998年的樣本信息求得參數(shù)的后驗(yàn)密度，該后驗(yàn)密度中包含了先驗(yàn)信息和樣本信息；然后將得到的后驗(yàn)密度結(jié)果作為下一次模型更新的先驗(yàn)信息，再結(jié)合下一周期(1998和2003年)的樣本信息求得更新后的后驗(yàn)密度，以此類推重復(fù)上述過程，不斷加入新的樣本信息進(jìn)行模型參數(shù)的后驗(yàn)樣本更新，直到將訓(xùn)練樣本中4次測量的樣本信息全部包含到后驗(yàn)密度中，最后利用生成的參數(shù)后驗(yàn)樣本即可估計(jì)出最終的模型參數(shù)值。

該過程通過WinBUGS軟件的MCMC算法和Gibbs抽樣算法實(shí)現(xiàn)。在根據(jù)貝葉斯方法進(jìn)行參數(shù)的后驗(yàn)樣本估計(jì)時(shí)，為了保證迭代收斂以得到穩(wěn)定的參數(shù)后驗(yàn)概率值，本研究設(shè)定的迭代次數(shù)為1 000萬次，其中每隔1 000次迭代記錄一次樣本數(shù)據(jù)，共搜集了10 000個(gè)樣本。圖2展示了模型中部分參數(shù)(a1，a2，a3)的迭代軌跡，可以看出，大約經(jīng)過100次樣本記錄(即10萬次迭代)之后所有參數(shù)趨于收斂，此時(shí)前面不收斂的10萬次迭代稱為退火迭代。在選取后驗(yàn)樣本進(jìn)行參數(shù)估計(jì)時(shí)需要排除退火迭代的影響，因此本研究截取了迭代100萬次之后所記錄的共計(jì)9 000個(gè)后驗(yàn)樣本，得到的部分參數(shù)(a1，a2，a3，a4)的后驗(yàn)密度如圖3所示。利用后驗(yàn)樣本估計(jì)最終轉(zhuǎn)移矩陣模型的全部參數(shù)，結(jié)果見表3。

圖2 部分參數(shù)迭代軌跡

3.3 模型比較

利用上述兩種模型的參數(shù)結(jié)果所建立的轉(zhuǎn)移概率矩陣模型對該林分2013年的徑級分布進(jìn)行預(yù)測，將預(yù)測結(jié)果與測試樣本中2013年的真實(shí)結(jié)果進(jìn)行比較，以此判斷兩種模型預(yù)測的準(zhǔn)確性。

兩種模型預(yù)測的徑級分布與實(shí)際分布如圖4所示。可見，在小徑級(4～8 cm)和大徑級(48 cm以上)處，兩模型的預(yù)測值與真實(shí)值差異不大；對于中徑級(12～44 cm)立木數(shù)量的預(yù)測分布，貝葉斯模型的預(yù)測值與真實(shí)值較接近，而固定參數(shù)模型的預(yù)測值比貝葉斯模型的預(yù)測值和真實(shí)值都要偏大。因此，貝葉斯模型比固定參數(shù)模型的預(yù)測效果更好。圖5反映了兩種模型的預(yù)測值與林分分布真實(shí)值之間的差異，圖中的對角線表示預(yù)測值與實(shí)際觀測值的等值線?？梢钥闯觯潭▍?shù)模型的預(yù)測值絕大多數(shù)分布在等值線左上方，貝葉斯模型的預(yù)測值則在等值線附近輕微波動(dòng)。這進(jìn)一步證明了固定參數(shù)模型的預(yù)測值比真實(shí)值偏高的趨勢，而貝葉斯模型則較精確地預(yù)測了林分的真實(shí)分布。

圖3 部分參數(shù)的后驗(yàn)密度

表3 貝葉斯模型的參數(shù)估計(jì)

為了更加科學(xué)地衡量兩種模型的預(yù)測效果，分別計(jì)算了兩種模型預(yù)測結(jié)果的均方誤差MSE、均方根誤差RMSE以及平均絕對百分誤差MAPE，結(jié)果見表4。

可見，貝葉斯模型的各項(xiàng)指標(biāo)均小于固定參數(shù)模型，說明貝葉斯模型的預(yù)測效果優(yōu)于固定參數(shù)的模型。

表4 模型預(yù)測效果指標(biāo)

圖4 兩種模型的預(yù)測分布與觀測分布

圖5 兩種模型預(yù)測值與觀測值的比較

4 結(jié)論與討論

采用固定參數(shù)平均的方法所建立的轉(zhuǎn)移矩陣模型在進(jìn)行實(shí)際林分結(jié)構(gòu)預(yù)測時(shí)容易產(chǎn)生比實(shí)際值偏高的預(yù)測結(jié)果，而采用貝葉斯方法所建立的模型則沒有這種趨勢。與經(jīng)典統(tǒng)計(jì)學(xué)方法相比，貝葉斯方法有效提高了模型的預(yù)測精度，是一種更好的模型參數(shù)估計(jì)方法。

從貝葉斯方法的參數(shù)估計(jì)過程可以看出，其利用先驗(yàn)信息和樣本信息計(jì)算后驗(yàn)分布的方法可有效利用定期更新的樣本數(shù)據(jù)不斷調(diào)整模型參數(shù)的后驗(yàn)分布，從而得到更符合實(shí)際的模型參數(shù)。該方法特別適合類似于森林資源連續(xù)調(diào)查的數(shù)據(jù)，每一次調(diào)查分析的結(jié)果都可以作為下一次調(diào)查分析的先驗(yàn)信息，使得調(diào)查數(shù)據(jù)得到了充分的利用。因此，貝葉斯方法很適合在一些林業(yè)數(shù)據(jù)的模型參數(shù)估計(jì)上進(jìn)行推廣應(yīng)用。

盡管如此，貝葉斯方法目前仍存在一些不足之處。其計(jì)算后驗(yàn)密度的方法需要依賴大量的隨機(jī)模擬實(shí)驗(yàn)，盡管MCMC法結(jié)合Gibbs抽樣算法可以有效減少計(jì)算后驗(yàn)密度的運(yùn)算量，但隨著模型參數(shù)數(shù)量的增加，隨機(jī)馬爾科夫鏈?zhǔn)諗康乃俣戎饾u變慢。因此，對于參數(shù)較多的模型采用貝葉斯方法進(jìn)行建模尚存在一定的局限性。這種局限性會(huì)隨著計(jì)算機(jī)計(jì)算能力的發(fā)展以及算法的優(yōu)化逐步被打破。

[1] 張雄清，張建國，段愛國.基于貝葉斯法估計(jì)杉木人工林樹高生長模型[J].林業(yè)科學(xué)，2014,50(3):69-75.

[2] 賈乃光.貝葉斯統(tǒng)計(jì)學(xué)及其應(yīng)用于林業(yè)的前景[J].北京林業(yè)大學(xué)學(xué)報(bào)，1988(4):91-96.

[3] 張雄清，張建國，段愛國.杉木人工林林分?jǐn)嗝娣e生長模型的貝葉斯法估計(jì)[J].林業(yè)科學(xué)研究，2015,28(4):538-542.

[4] 姚丹丹，雷相東，張則路.基于貝葉斯法的長白落葉松林分優(yōu)勢高生長模型研究[J].北京林業(yè)大學(xué)學(xué)報(bào)，2015,37(3):94-100.

[5] 張玉環(huán)，張青，亢新剛，等.異齡針闊混交擇伐林均衡曲線的確定方法：以金溝嶺林場樣地?cái)?shù)據(jù)為例[J].北京林業(yè)大學(xué)學(xué)報(bào)，2015,37(6):53-60.

[6] LESLIE P H . On the use of matrices in certain population mathematics[J]. Biometrika，1945,33(5):183-212.

[7] USHER M B. A matrix approach to the management of renewable resources, with special reference to selection forests[J]. Journal of Applied Ecology，1966,3(2):347-348.

[8] BUONGIORNO J, PEYRON J L, HOULLIER F, et al. Growth and management of mixed-species, uneven-aged forests in the French Jura: implications for economic returns and tree diversity[J]. Forest Science，1995,41(3):397-429.

[9] SOLOMON D S, HOSMER R A, HAYSLETT H T . A two-stage matrix model for predicting growth of forest stands in the Northeast[J]. Canadian Journal of Forest Research，1986,16(3):521-528.

[10] YUE C H, KLDTKE J, LENK E. Ein verfahren zur bestimmung zielorientierter gleichgewichtskurven im plenterwald[J]. Berichte der Sektion Biometrie und Informatik im DVFF，1997,174:145-154.

[12] 趙俊卉，亢新剛，張慧東，等.長白山3個(gè)主要針葉樹種的標(biāo)準(zhǔn)樹高曲線[J].林業(yè)科學(xué)，2010,46(10):191-194.

[13] 劉洋，亢新剛，郭艷榮，等.長白山主要樹種直徑生長的多元回歸預(yù)測模型:以云杉為例[J].東北林業(yè)大學(xué)學(xué)報(bào)，2012,40(2):1-4.

[14] 胡云云，閔志強(qiáng)，高延，等.擇伐對天然云冷杉林林分生長和結(jié)構(gòu)的影響[J].林業(yè)科學(xué)，2011,47(2):15-24.

[15] HUANG S. Diameter and height growth models[D]. Edmonton: University of Alberta,1992.

[16] ALBERT J. Bayesian computation with R[M]. New York： Springer Science & Business Media,2009.

[17] 張偉.基于MCMC的貝葉斯模型平均組合預(yù)測方法及其在我國能源消費(fèi)預(yù)測中的應(yīng)用研究[D].重慶：重慶師范大學(xué),2011.

[18] SPIEGELHALTER D, THOMAS A, BEST N. WinBUGS user manual[M]. Cambridge: MRC Biostatistics Unit，2003.

[19] CHIB S, GREENBERG E. Understanding the Metropolis-Hastings algorithm[J]. American Statistician，1995,49(4):327-335.

[20] KABACOFF R,陳鋼，肖楠，等.R語言實(shí)戰(zhàn)[M].北京：人民郵電出版社,2016.

Bayesian Method in Estimating Model Parameters to Predict Diameter Class Distribution of Stand-Taking Sample Plot Data in Jingouling Forest Farm as Example//

Yu Ruchuan, Zhang Qing
(Beijing Forestry University, Beijing 100083, P. R. China);
Yue Chaofang
(Forest Research Institute of Baden-Wurttemberg, Freiburg 79100, F. R. Germany);
Kang Xingang
(Beijing Forestry University)//Journal of Northeast Forestry University，2017,45(6)：30-35.

Transition matrix model; Stand diameter class distribution; Bayesian statistics; MCMC method; Gibbs sampling algorithm

1)“948”國家林業(yè)局引進(jìn)項(xiàng)目(2013-4-66)。

于汝川，男，1991年9月生，北京林業(yè)大學(xué)理學(xué)院，碩士研究生。E-mail：yuruchuan@qq.com。

張青，北京林業(yè)大學(xué)理學(xué)院，教授。E-mail：zhangq@bjfu.edu.cn。

2017年2月27日。

S758.5

責(zé)任編輯：戴芳天。

With the observed data in five different years from 26 sample plots in Jingouling Forest Farm of Wangqing Forestry Bureau, Jilin Province, transition matrix models were established to predict the stand diameter class distribution in a given period. Parameters of the models were estimated by classical statistical method and Bayesian method, respectively. After the fixed parameter model and the Bayesian model were built, their prediction results could be compared to test the prediction accuracy of the two models. The predicted values of fixed parameter model tend to be higher than the true value. The prediction of Bayesian model is relatively more accurate than that of fixed parameter model.