基于分層貝葉斯方法的高山松單木生物量模型

黃金君，舒清態(tài)，王柯人，席 磊，孫 楊，羅 浩

(西南林業(yè)大學 林學院，云南 昆明 650224)

準確估算森林生物量對于研究固碳調控和評估森林生態(tài)系統(tǒng)的健康至關重要[1]。單木生物量模型是進行林分和區(qū)域尺度生物量估測的基礎，使用不同的擬合方法對于提升單木模型參數(shù)的估計已成為當前森林生物量研究重要任務之一。傳統(tǒng)的單木生物量模型一般采用較成熟的異速生長方程，通過將易測的胸徑或樹高等指標結合以提供一種高效簡單的估算森林生物量方法[2-5]。然而使用傳統(tǒng)方法擬合異速生物量模型存在過多局限，其只將模型參數(shù)視作固定變量，未能綜合考慮存在的多種生物或非生物因素；且若樣本數(shù)據(jù)存在地域差異，精度會明顯降低[6-8]。因此選用綜合參數(shù)與非參數(shù)模型優(yōu)勢的擬合方法，如混合模型法與分層貝葉斯方法有助于解決這一問題?；旌夏Ｐ屯ㄟ^引入隨機效應因子反映了整體與個體之間的相關性與差異性，有效地彌補了傳統(tǒng)統(tǒng)計學方法存在的缺陷[9]。Y.Zhangetal[10]在擬合不同組分的單木生物量模型時，得出使用含有隨機效應的非線性混合模型法可以有效提高模型估測精度；陳哲夫等[11]以樣地為隨機效應構建馬尾松(Pinusmassoniana)次生林單木生長模型，結果為含有隨機效應混合模型法的擬合效果與預測精度均優(yōu)于傳統(tǒng)模型。貝葉斯學派利用概率分布描述模型中的未知變量，并將樣本信息與參量視作隨機變量[12]。在當前精準提升森林質量的背景下，貝葉斯方法正逐漸成為林業(yè)領域中學者們研究的熱點，包括樹木死亡率模型[13]、杉木(Cunninghamialanceolata)人工林最大密度線[14]、單木樹高-胸徑模型[15-16]等。這些研究表明，貝葉斯統(tǒng)計法可以提高模型參數(shù)的穩(wěn)定性，使得預測效果更可靠。分層貝葉斯方法是貝葉斯統(tǒng)計學中的一種，通過分離不同層面中多個參數(shù)的復雜關系，從而充分發(fā)揮具有層次結構模型的優(yōu)勢[17]。Q.M.Ketteringsetal[18]通過引入具有隨機效應參數(shù)和固定效應參數(shù)的混合模型法與分層貝葉斯方法，可以有效解決區(qū)域差異對森林生物量的影響；黃興召等[19]證實設置區(qū)域為隨機效應的混合模型法與分層貝葉斯方法，可以顯著提升估算林分生物量的精度。本研究以香格里拉市高山松(Pinusdensata)林為對象，以胸徑和樹高作為解釋變量，以傳統(tǒng)較成熟的異速生長方程為基礎模型，分別采用分層貝葉斯方法、非線性混合模型法、貝葉斯方法和非線性最小二乘法進行擬合，構建高精度單木地上生物量估測模型，闡述區(qū)域差異對估算單木地上生物量的影響，為低緯度高海拔地區(qū)典型森林碳匯經(jīng)營提供技術支持。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)地處“三江并流”核心區(qū)，即云南省迪慶州香格里拉市(26°52′-28°50′N，99°23′-100°18′E)。平均海拔3 459 m，多年平均氣溫4.7～16.5 ℃，年平均降水量649.4 mm。香格里拉市的森林覆蓋率已達76%，根據(jù)《云南植被》劃分標準，植被分布類型包括落葉闊葉林、暖性針葉林、溫性針葉林、硬葉常綠闊葉林、高山湖泊水生植被以及亞熱性常綠闊葉林等十多種。本研究區(qū)的高山松分布廣泛，是高寒山區(qū)最為典型的高山松和先鋒喬木樹種,具有耐貧瘠、更新快和環(huán)境適應性好等特征。

1.2 數(shù)據(jù)采集

所用數(shù)據(jù)來自公益性林業(yè)科研專項經(jīng)費項目(201404309)的立木地上生物量實測數(shù)據(jù)，以香格里拉市為研究區(qū)共設置面積為30 m×30 m的樣地56個。在樣地調查的基礎上展開樣木調查，分別在56個不同立地條件標準樣地中按照徑階分布(2～3株)選取區(qū)域Ⅰ(99.759 57°～99.762 90°E，27.625 47°～27.628 42°N)60株樣木和區(qū)域Ⅱ(100.060 64°～100.066 34°E，27.723 27°～27.727 18°N)55株樣木，共計115株健康且具有代表性的高山松(表1)。高山松樣木地上不同組分的生物量(如樹干、樹枝、樹冠、樹皮以及樹葉生物量等)采用分層切割法取樣測定，并帶回實驗室分別在105 ℃和72 ℃下恒溫烘干，對單木各部分鮮重和干重進行換算求和，最終得到高山松單木地上總生物量。

1.3 基礎模型構造

選擇異速生長方程作為高山松單木生物量的基礎模型，其通過破壞少量有限樣木，建立地上生物量與較易獲取參數(shù)指標(冠幅、林分狀況、年齡、樹高和胸徑等)之間的異速生長關系。計算公式如下。

M=aDbHc

(1)

式中：M代表立木地上生物量，D代表胸徑，H代表樹高，a、b、c為模型代估參數(shù)。

1.4 模型擬合方法

1.4.1 非線性最小二乘法 最小二進位乘法算法是在我國經(jīng)典統(tǒng)計學中一種最常見和普遍采用的一種逼近算法，其基本原理是通過不斷地替換和迭代線性函數(shù)來逼近非線性函數(shù)這一步的過程，可以求出誤差平方和最小，從而求出參數(shù)的最優(yōu)解。公式為：

(2)

表1 高山松實測樣木各變量統(tǒng)計值

1.4.2 非線性混合模型法 非線性混合模型法通過在非線性回歸模型中引入反映整體中不同部分個體變化的隨機效應因子，其自變量固定效應和隨機效應與因變量之間的聯(lián)系為非線性[20-21]，因此被視作是線性混合模型法與非線性回歸模型法的折中。表達式為：

(3)

1.4.3 貝葉斯方法 相較于經(jīng)典統(tǒng)計學的方法(最小二乘法)，貝葉斯方法增加了對先驗信息的考慮，且把樣本與未知參數(shù)都視為隨機變量。該方法結合未知參數(shù)θ的先驗分布與傳統(tǒng)似然函數(shù)，通過MCMC算法和GIBBS抽樣生成馬爾科夫鏈，對馬爾科夫鏈進行不斷迭代最終得到參數(shù)的后驗估計。貝葉斯公式的密度函數(shù)表達式為

(4)

1.4.4 分層貝葉斯方法 分層貝葉斯方法的核心思想是在原來的貝葉斯方法中加入分層先驗信息，構造分層先驗可有助于消除先驗分布對結果的過度影響。具體思路是：首先將式(1)中的未知參數(shù)a、b和c定義為可以使用概率分布的形式表達出來的先驗分布θ～T1(θ/λ)，其中λ代表超參數(shù)。當λ無法確定時可以再構造一個超先驗分布T2(λ)，由先驗θ～T1(θ/λ)與超先驗T2(λ)共同組成的新先驗構成了分層貝葉斯的先驗分布T(θ)，其表達式為

(5)

1.4.5 先驗信息 貝葉斯方法中先驗分布信息的選擇至關重要，豐富的信息先驗影響著模型參數(shù)的估計[22]。在上述高山松單木生物量異速生長方程中，需要為所求參數(shù)構造適當?shù)南闰灧植?。有些研究人員建議采用無信息先驗，無信息先驗通常出現(xiàn)在方差無窮大與均值為0的高斯分布中，對參數(shù)不會造成太大的影響。也可采用有信息先驗(informative prior)作為貝葉斯方法的先驗分布信息，其來源于統(tǒng)計學方法計算得到的數(shù)據(jù)或前人的研究成果等。本研究采用非線性最小二乘法的估計結果作為貝葉斯方法與分層貝葉斯方法的先驗信息。

1.4.6 模型精度評價 選用決定系數(shù)(R2)、均方根誤差(RMSE，公式中用RMSE表示)和估測精度(E)作為分層貝葉斯方法、非線性混合模型法、貝葉斯方法和非線性最小二乘法擬合異速生長方程的統(tǒng)計量指標：

(6)

(7)

(8)

2 結果與分析

2.1 貝葉斯方法和非線性最小二乘法的擬合結果分析

使用貝葉斯方法和非線性最小二乘法的擬合結果見表2，以非線性最小二乘法計算得到的均值(Mean)和標準差(SD)作為貝葉斯方法中參數(shù)的先驗信息，因此這2種方法的均值和標準差均較為接近。貝葉斯方法和非線性最小二乘法的決定系數(shù)(R2)、均方根誤差(RMSE)和估測精度(E)分別為0.981 8、0.984 7，82.87%、84.29%和44.68，40.96。根據(jù)R2與E越大而RMSE越小，模型擬合效果越好的準則，非線性最小二乘法的模型擬合效果較優(yōu)于貝葉斯方法。但是通過對比這2種方法中模型參數(shù)的置信區(qū)間發(fā)現(xiàn)，貝葉斯方法擬合的b和c參數(shù)的區(qū)間范圍均小于非線性最小二乘法的b和c參數(shù)的區(qū)間范圍，這也說明了貝葉斯方法可以提高模型估算的穩(wěn)定性，使得參數(shù)分布更加集中。

表2 貝葉斯方法和非線性最小二乘法的擬合結果對比

2.2 分層貝葉斯方法和非線性混合模型法的擬合結果分析

使用模型擬合的參數(shù)估計值與評價指標見表3，a1、b1和c1以及a2、b2和c2分別代表分層貝葉斯方法中Ⅰ區(qū)和Ⅱ區(qū)單木地上生物量異速生長方程的參數(shù)估計值；ωa、ωb和ωc代表非線性混合模型法中參數(shù)a、b和c的固定效應參數(shù)值，μa1和μa2代表作用在參數(shù)a1和a2的區(qū)域隨機效應值，其相同的隨機構造變量在同一參數(shù)的效應參數(shù)值總和為0，μa12代表矩陣為1×1的隨機效應方差估計值。分層貝葉斯方法與非線性混合模型法均考慮了區(qū)域隨機效應因子，因此這2種方法的決定系數(shù)(R2)、均方根誤差(RMSE)和估測精度(E)較為接近。但仔細對比這2種方法的差異發(fā)現(xiàn)，分層貝葉斯方法相對于非線性混合模型法R2與E依次提高了0.000 1與0.07%，RMSE降低了0.18 kg，因此模型估計效果最好的是分層貝葉斯方法。

表3 分層貝葉斯方法和非線性混合模型法的擬合結果對比

2.3 4種方法的擬合結果分析

4種方法擬合高山松單木生物量模型的結果見表2和表3，分層貝葉斯方法和非線性混合模型法的R2分別為0.985 6、0.985 5，E分別為84.76%、84.69%，RMSE分別為39.75、39.93；貝葉斯方法和非線性最小二乘法的R2分別為0.981 8、0.984 7，E分別為82.87%、84.29%，RMSE分別為44.68、40.96。由此得出，加入了區(qū)域隨機效應的分層貝葉斯方法和非線性混合模型法的擬合效果均優(yōu)于未加入?yún)^(qū)域隨機效應的貝葉斯方法和非線性最小二乘法。同時，這4種方法的擬合精度R2均達到了0.98及以上且P<0.001，表明均適合該模型的擬合，其估計結果見圖1～圖4。

3 結論與討論

3.1 結論

以胸徑和樹高為解釋變量的異速生長方程為單木生物量基礎模型，采用分層貝葉斯方法、非線性混合模型法、貝葉斯方法和非線性最小二乘法求解基礎模型的異速生長參數(shù)，分析不同方法擬合模型參數(shù)的表現(xiàn)和估計精度。從擬合精度看，利用4種方法擬合單木生物量模型參數(shù)的決定系數(shù)與估計精度分別達到了0.98和82%以上，估計效果均較優(yōu)。通過對比不同方法的差異發(fā)現(xiàn)，加入了區(qū)域隨機效應的分層貝葉斯方法和非線性混合模型法的擬合效果均優(yōu)于未加入?yún)^(qū)域隨機效應的貝葉斯方法和非線性最小二乘法，證實了以區(qū)域為隨機效應與分層為基礎的分層貝葉斯方法和非線性混合模型法可以提高單木生物量模型的估計精度，為大尺度樣本數(shù)據(jù)模型參數(shù)估測方法提供新思路?？偟膩砜?，基于分層貝葉斯方法的高山松單木生物量模型擬合效果最優(yōu)，因此應增加分層貝葉斯方法在林業(yè)領域中的研究，且可通過完善先驗信息的選擇和考慮多種隨機效應來構造多層貝葉斯模型。

3.2 討論

研究中利用分層貝葉斯方法和非線性混合模型法的模型擬合效果均優(yōu)于貝葉斯方法和非線性最小二乘法，表明區(qū)域是造成單木生物量差異的主要原因。李春明等[24]利用陜西省10個區(qū)域39個樣地的胸徑與樹高數(shù)據(jù)，得出考慮區(qū)域效應或樣地效應的混合模型均比未考慮區(qū)域效應或樣地效應的固定模型的模型精度高，進一步證實了構造區(qū)域隨機效應因子對估算立木生物量存在顯著影響。所構造的隨機效應參量值或不同隨機效應因子之間的差異越大，4種方法的估算結果差異將會越明顯。研究中非線性最小二乘法的模型擬合效果較優(yōu)于貝葉斯方法，但是貝葉斯方法中參數(shù)的置信區(qū)間范圍更加集中，說明基于貝葉斯方法的模型擬合效果更具穩(wěn)定性。這是由于當樣本數(shù)量較大時，貝葉斯方法與經(jīng)典統(tǒng)計學方法(最小二乘法或最大似然法)的模型擬合效果相似；當樣本數(shù)量較少時，貝葉斯方法的擬合效果則優(yōu)于經(jīng)典統(tǒng)計學方法[25]。姚丹丹等[26]證實使用較少樣本數(shù)量時,貝葉斯方法的估計精度和穩(wěn)定性優(yōu)于經(jīng)典統(tǒng)計學方法。

分層貝葉斯方法是一種基于貝葉斯理論發(fā)展起來的有效描述復雜數(shù)據(jù)集和評估參數(shù)不確定性的方法[27]，現(xiàn)可以基于MCMC方法解決復雜高維的計算難題。D.Chenetal[28]通過采集湖北、甘肅、河北、遼寧、黑龍江、內蒙古等不同地區(qū)的數(shù)據(jù)，得出具有先驗信息的分層貝葉斯方法更適合擬合具有區(qū)域變化的生物量模型。貝葉斯理論認為在構建模型過程中，將從多種渠道獲取數(shù)據(jù)的先驗信息與樣本數(shù)據(jù)結合進而提高模型估測的質量，先驗信息的優(yōu)勢在數(shù)據(jù)較難獲取的情況下會表現(xiàn)得更明顯。在研究樹木增長率與死亡率之間聯(lián)系的過程中，其死亡率的數(shù)據(jù)獲取需耗費大量的時間與成本，P.H.Wyckoffetal[29]使用貝葉斯方法克服樣本量小的難題進而對林分死亡率進行推算。