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

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

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

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。香格里拉市的森林覆蓋率已達(dá)76%，根據(jù)《云南植被》劃分標(biāo)準(zhǔn)，植被分布類型包括落葉闊葉林、暖性針葉林、溫性針葉林、硬葉常綠闊葉林、高山湖泊水生植被以及亞熱性常綠闊葉林等十多種。本研究區(qū)的高山松分布廣泛，是高寒山區(qū)最為典型的高山松和先鋒喬木樹(shù)種,具有耐貧瘠、更新快和環(huán)境適應(yīng)性好等特征。

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

所用數(shù)據(jù)來(lái)自公益性林業(yè)科研專項(xiàng)經(jīng)費(fèi)項(xiàng)目(201404309)的立木地上生物量實(shí)測(cè)數(shù)據(jù)，以香格里拉市為研究區(qū)共設(shè)置面積為30 m×30 m的樣地56個(gè)。在樣地調(diào)查的基礎(chǔ)上展開(kāi)樣木調(diào)查，分別在56個(gè)不同立地條件標(biāo)準(zhǔn)樣地中按照徑階分布(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株樣木，共計(jì)115株健康且具有代表性的高山松(表1)。高山松樣木地上不同組分的生物量(如樹(shù)干、樹(shù)枝、樹(shù)冠、樹(shù)皮以及樹(shù)葉生物量等)采用分層切割法取樣測(cè)定，并帶回實(shí)驗(yàn)室分別在105 ℃和72 ℃下恒溫烘干，對(duì)單木各部分鮮重和干重進(jìn)行換算求和，最終得到高山松單木地上總生物量。

1.3 基礎(chǔ)模型構(gòu)造

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

M=aDbHc

(1)

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

1.4 模型擬合方法

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

(2)

表1 高山松實(shí)測(cè)樣木各變量統(tǒng)計(jì)值

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

(3)

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

(4)

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

(5)

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

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

(6)

(7)

(8)

2 結(jié)果與分析

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

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

表2 貝葉斯方法和非線性最小二乘法的擬合結(jié)果對(duì)比

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

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

表3 分層貝葉斯方法和非線性混合模型法的擬合結(jié)果對(duì)比

2.3 4種方法的擬合結(jié)果分析

4種方法擬合高山松單木生物量模型的結(jié)果見(jiàn)表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ū)域隨機(jī)效應(yīng)的分層貝葉斯方法和非線性混合模型法的擬合效果均優(yōu)于未加入?yún)^(qū)域隨機(jī)效應(yīng)的貝葉斯方法和非線性最小二乘法。同時(shí)，這4種方法的擬合精度R2均達(dá)到了0.98及以上且P<0.001，表明均適合該模型的擬合，其估計(jì)結(jié)果見(jiàn)圖1～圖4。

3 結(jié)論與討論

3.1 結(jié)論

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

3.2 討論

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

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