多源數(shù)據(jù)對(duì)林分動(dòng)態(tài)預(yù)測(cè)的影響及不確定性分析*

田相林 廖梓延 孫帥超 薛海連 王 彬 曹田健

1.西北農(nóng)林科技大學(xué)林學(xué)院 生態(tài)仿真優(yōu)化實(shí)驗(yàn)室 楊凌 712100； 2.赫爾辛基大學(xué)林學(xué)系 赫爾辛基 FI-00014； 3.中國科學(xué)院成都生物研究所 成都 610041； 4.中國科學(xué)院大學(xué) 北京 100039； 5.福建農(nóng)林大學(xué) 福州 350002； 6.西北農(nóng)林科技大學(xué)理學(xué)院 楊凌 712100； 7.青海大學(xué)農(nóng)林科學(xué)院 西寧 810016)

在森林立地、密度效應(yīng)、林分生長(zhǎng)以及經(jīng)營措施等方面的理論研究中，森林建模往往要求林分觀測(cè)數(shù)據(jù)是一種控制試驗(yàn)的產(chǎn)物，并假定試驗(yàn)是被設(shè)計(jì)在嚴(yán)格可控的環(huán)境中以檢驗(yàn)特定變量對(duì)林分生長(zhǎng)動(dòng)態(tài)的影響。然而，受限于控制森林生態(tài)系統(tǒng)環(huán)境的能力、調(diào)查所能承擔(dān)的時(shí)間和空間范圍以及觀測(cè)到特定變量的能力等，森林動(dòng)態(tài)預(yù)測(cè)研究通常依賴于有限數(shù)量并伴隨復(fù)雜誤差結(jié)構(gòu)的數(shù)據(jù)，使得經(jīng)典統(tǒng)計(jì)模型的假設(shè)情境往往難以被滿足(Clark， 2007)。在森林生態(tài)系統(tǒng)研究中，林分生長(zhǎng)收獲分析需要綜合多種數(shù)據(jù)類型，如臨時(shí)樣地、固定樣地、解析木和生長(zhǎng)錐木芯等，這些信息由于采用不同抽樣設(shè)計(jì)或觀測(cè)方法等原因常常具有極其復(fù)雜的誤差結(jié)構(gòu)特征，且其復(fù)雜程度會(huì)隨不同時(shí)期或類型的新數(shù)據(jù)加入而不斷累積。同時(shí)，在林分動(dòng)態(tài)預(yù)測(cè)與生長(zhǎng)收獲模型應(yīng)用研究中，來自非控制條件下部分觀測(cè)過程的信息往往需要作為統(tǒng)計(jì)推斷的基礎(chǔ)(Dietze， 2017)，諸多不可控的環(huán)境或機(jī)制因素都會(huì)導(dǎo)致森林調(diào)查數(shù)據(jù)與林分動(dòng)態(tài)預(yù)測(cè)之間形成復(fù)雜的信息交互體系，最終使得模擬與預(yù)測(cè)的可靠性無法被量化評(píng)價(jià)。

隨著生態(tài)系統(tǒng)探討問題的復(fù)雜性不斷升高，現(xiàn)代模擬和計(jì)算技術(shù)被廣泛應(yīng)用于相關(guān)領(lǐng)域研究中，其中，貝葉斯框架的潛在假設(shè)使得貝葉斯方法不再受模型數(shù)據(jù)量和控制條件等因素制約，且具備綜合多源數(shù)據(jù)信息的能力(Clarketal.， 2007)，其應(yīng)用在20世紀(jì)90年代后飛速發(fā)展(Gelfandetal.， 1990； Berger， 2000)，并逐漸被林業(yè)研究者所關(guān)注(Greenetal.， 1996； Bullocketal.， 2007； Clarketal.， 2007； Metcalfetal.， 2009)。貝葉斯估計(jì)值同時(shí)依賴于先驗(yàn)和后驗(yàn)分布，復(fù)雜問題的后驗(yàn)信息往往由迭代數(shù)據(jù)隨機(jī)取樣來確定，即產(chǎn)生參數(shù)的完整后驗(yàn)分布，這意味著許多統(tǒng)計(jì)值可從生成的樣本中獲得。盡管經(jīng)典的最大似然方法和貝葉斯方法會(huì)產(chǎn)生相似的參數(shù)估計(jì)(Lietal.， 2011)，但貝葉斯方法不需要像最大似然方法那樣的附加條件和特定的分布假設(shè)，貝葉斯框架下林分模型的構(gòu)建，可極大限度減小統(tǒng)計(jì)學(xué)假設(shè)對(duì)數(shù)據(jù)的限制，并能對(duì)模型本身和林分生長(zhǎng)的不確定性進(jìn)行合理量化分析(Hartigetal.， 2012)?；赥ylor展開式或Monte Carlo抽樣方法推斷或預(yù)測(cè)的不確定性可以解析為不同的來源，包括數(shù)據(jù)的不確定性、模型結(jié)構(gòu)的不確定性和模型參數(shù)的不確定性，進(jìn)而對(duì)不確定的傳遞過程進(jìn)行量化(Lebaueretal.， 2013； Dietze， 2017)，其中多維相關(guān)性基礎(chǔ)上的參數(shù)不確定性是生態(tài)系統(tǒng)模型開發(fā)與應(yīng)用中難以回避的內(nèi)容(Van Oijen, 2017)。

在森林動(dòng)態(tài)預(yù)測(cè)研究中，貝葉斯方法也逐步得以應(yīng)用，如地上生物量模型(Zapata-Cuartasetal.， 2012； Zhangetal.， 2013)、直徑分布模型(Bullocketal.， 2007)和直徑生長(zhǎng)模型(Clarketal.， 2007)等。貝葉斯方法的核心優(yōu)勢(shì)在于其對(duì)數(shù)據(jù)與模型關(guān)系的描述，如估計(jì)多層次/體系模型(Gelmanetal.， 2007； Dietzeetal.， 2008)、利用有限數(shù)據(jù)信息校正復(fù)雜模型(Van Oijenetal.， 2005)、從不同模型中合并預(yù)測(cè)值(Radtkeetal.， 2006)以及克服數(shù)據(jù)稀疏問題(Metcalfetal.， 2009)等，而目前森林建模研究中，貝葉斯框架下數(shù)據(jù)與模型間信息融合與復(fù)雜誤差結(jié)構(gòu)的討論往往集中于過程機(jī)制模型而非經(jīng)驗(yàn)?zāi)Ｐ?Hartigetal.， 2012； Dietze, 2017； Van Oijen， 2017)。森林生長(zhǎng)收獲經(jīng)驗(yàn)?zāi)Ｐ偷闹髁魅匀皇腔谙到y(tǒng)性設(shè)計(jì)的長(zhǎng)期觀測(cè)(Pretzsch， 2010)，并對(duì)數(shù)據(jù)的誤差體系有著詳盡考量(Kangasetal.， 2006)； 然而這種大量建模數(shù)據(jù)的需求，在國內(nèi)多數(shù)林區(qū)難以滿足。由于缺少系統(tǒng)性的長(zhǎng)期觀測(cè)，且數(shù)據(jù)誤差隨不同抽樣策略或數(shù)據(jù)源而劇烈波動(dòng)，使得經(jīng)典回歸分析理論中的許多假設(shè)并不成立，因此將誤差結(jié)構(gòu)納入模型體系尤為重要(唐守正等， 1995； 1998； 符利勇等， 2014)。

從森林生長(zhǎng)收獲模型的尺度上來說，全林模型只需要相對(duì)較少的信息來模擬林分生長(zhǎng)。相比徑階和單木模型，全林模型所模擬的基本單元是林分和立地參數(shù)，包括年齡、林分?jǐn)嗝娣e、立地指數(shù)和蓄積量等，雖然無法像徑階或單木模型一樣分析林分結(jié)構(gòu)的動(dòng)態(tài)變化，但全林模型可以保證較為準(zhǔn)確的未來林分蓄積生長(zhǎng)狀況信息，這在森林經(jīng)營實(shí)踐中非常重要(Weiskitteletal.， 2011)?？勺兠芏热帜Ｐ蛯⒘址置芏茸鳛槔^林分年齡、立地質(zhì)量的第3個(gè)變量，可以描述特定樹種和立地條件下任意林分的生長(zhǎng)，在密度控制和森林經(jīng)營規(guī)劃中有著廣泛應(yīng)用(張少昂， 1986； 李希菲等， 1988； 唐守正， 1991； 1993； 1999； 洪玲霞， 1993； 洪玲霞等， 2012； Weiskitteletal.， 2011)。

本研究以數(shù)據(jù)信息要求較低、應(yīng)用廣泛的可變密度全林模型為例，采用貝葉斯數(shù)據(jù)-模型框架，構(gòu)建森林資源連續(xù)調(diào)查中數(shù)據(jù)信息融合與模型學(xué)習(xí)的關(guān)系，探究隨著新一期數(shù)據(jù)不斷加入，模型參數(shù)概率分布和模型預(yù)測(cè)結(jié)果的變化規(guī)律； 在考慮誤差分布和異方差等特征的基礎(chǔ)上，比較角規(guī)臨時(shí)樣地、矩形固定樣地、樹干解析3種數(shù)據(jù)類型對(duì)模型構(gòu)建和林分動(dòng)態(tài)模擬的影響，以提高模型預(yù)測(cè)的準(zhǔn)確性和可靠性； 同時(shí)，基于不確定性的量化評(píng)估，為森林調(diào)查中的數(shù)據(jù)收集策略提供建議。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)

模型開發(fā)數(shù)據(jù)源自陜西省寧陜縣火地塘林區(qū)(108°21′—108°39′E，33°18′—33°28′N)和旬陽壩林區(qū)(103°58′—109°48′E，32°29′—33°13′N)。數(shù)據(jù)采集海拔范圍為1 200～1 800 m，坡度范圍為15°～42°。土壤為森林棕壤，平均厚度50 cm。該區(qū)域?qū)賮啛釒駶櫳降貧夂?，年均氣?～12 ℃，年均降雨量1 100 mm。這2個(gè)林區(qū)原始天然林于20世紀(jì)60—70年代進(jìn)行過采伐，目前多為天然次生林或人工林。研究區(qū)樣地位置分布見圖1。

研究共使用215塊角規(guī)繞測(cè)臨時(shí)樣地、22塊20 m×20 m矩形固定樣地和11株中等木解析木，所有樣地均為油松(Pinustabulaeformis)純林(油松斷面積占比高于90%)或相對(duì)純林(油松斷面積占比65%～90%)。角規(guī)斷面積系數(shù)為1 m2·hm-2。角規(guī)繞測(cè)臨時(shí)樣地調(diào)查數(shù)據(jù)包括1990年森林經(jīng)理調(diào)查(李悅黎等， 1993)獲得的58塊油松相對(duì)純林樣地和2005年森林經(jīng)理調(diào)查獲得的68塊油松相對(duì)純林樣地。為減少抽樣設(shè)計(jì)帶來的誤差，采用完全相同的方法于2012年補(bǔ)充調(diào)查45塊油松樣地，最終獲得3期角規(guī)樣地調(diào)查數(shù)據(jù)(表1)。除了角規(guī)繞測(cè)臨時(shí)樣地以外的 22塊矩形固定樣地和11株油松中等木樹干解析數(shù)據(jù)，被用于比較不同數(shù)據(jù)類型對(duì)模型表現(xiàn)的影響，全部數(shù)據(jù)對(duì)比見表2。

本研究多源數(shù)據(jù)包括森林生長(zhǎng)收獲建模中2種常見的數(shù)據(jù)狀況： 一為相同抽樣設(shè)計(jì)下不同時(shí)期的調(diào)查數(shù)據(jù)(表1)，數(shù)據(jù)誤差結(jié)構(gòu)不同多為設(shè)備、人員和質(zhì)量控制差異造成； 二為不同抽樣設(shè)計(jì)本身導(dǎo)致誤差的系統(tǒng)性差異，即本研究中點(diǎn)抽樣的臨時(shí)樣地、每木檢尺的矩形固定樣地、解析木(表2)在抽樣和林分調(diào)查因子測(cè)算原理上存在系統(tǒng)性區(qū)別。

表1 用于模型循環(huán)更新測(cè)試的3期角規(guī)樣地調(diào)查數(shù)據(jù)

表2 模型構(gòu)建中比較的3種數(shù)據(jù)類型

1.2 可變密度全林模型

可變密度全林模型僅需對(duì)基礎(chǔ)林分調(diào)查信息進(jìn)行參數(shù)化，可用于模擬給定林分的平均高H(m)、斷面積G(m2·hm-2)、平均胸徑D(cm)等。模型體系如下：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

根據(jù)油松生長(zhǎng)曲線和數(shù)據(jù)中的樹高變動(dòng)系數(shù)，設(shè)置A0=30 年(吳恒等， 2015)。通過調(diào)整SDI和SCI，可模擬不同林分密度和立地質(zhì)量的林分生長(zhǎng)過程。整個(gè)模型體系不依賴于固定樣地的連續(xù)復(fù)測(cè)，林分密度變化由斷面積與平均胸徑的關(guān)系獲得。

關(guān)于本林區(qū)不同數(shù)據(jù)集和指標(biāo)對(duì)立地評(píng)價(jià)的具體影響以往研究已進(jìn)行了詳盡探討(吳恒等， 2015； 郭小陽等， 2017)，本研究之所以采用給定年齡下的林分平均高期望值，即地位級(jí)指數(shù)用于林分立地質(zhì)量評(píng)估，是因?yàn)?990和2005年的調(diào)查數(shù)據(jù)中許多樣地缺少優(yōu)勢(shì)木記錄，無法構(gòu)建常規(guī)立地指數(shù)曲線。林分平均高導(dǎo)向曲線、林分總斷面積模型和平均胸徑模型均采用Schumacher模型(Schumacher， 1939)的基本結(jié)構(gòu)。樹高模型擬合同時(shí)測(cè)試臨時(shí)樣地、固定樣地、解析木3種數(shù)據(jù)類型，胸徑和斷面積模型僅比較臨時(shí)樣地和固定樣地2種數(shù)據(jù)類型。林分密度指數(shù)推算基于完滿立木度林分密度(N)與平均胸徑(D)的關(guān)系，即需要完滿立木度林分?jǐn)?shù)據(jù)，采用逐步剔除(李希菲等， 1988)對(duì)樣地進(jìn)行篩選： 1) 利用全部數(shù)據(jù)擬合Reineke(1933)方程lnN= α-βlnD； 2) 舍棄所有低于該方程回歸線的低密度林分，并用剩余數(shù)據(jù)繼續(xù)擬合該方程； 3) 舍棄低于新回歸線的所有樣地。剩余樣地假定其為完滿立木度并重新擬合該方程，最終獲得式(1)中的參數(shù)β。分別完成地位級(jí)指數(shù)和林分密度指數(shù)模型擬合后，各樣地SCI和SDI通過計(jì)算得到，進(jìn)而依次完成林分?jǐn)嗝娣e和胸徑模型擬合。方程組聯(lián)立求解參數(shù)依1.4節(jié)操作。

1.3 貝葉斯數(shù)據(jù)-模型框架

貝葉斯方法提供數(shù)據(jù)信息融合與模型學(xué)習(xí)框架，新的數(shù)據(jù)信息可不斷加入模型中，該過程可以簡(jiǎn)單地表述為：

p(θ|D)∝p(D|θ)×p(θ)。

(8)

式中：p(θ|D)和p(θ)分別為模型參數(shù)θ的后驗(yàn)和先驗(yàn)分布；p(D|θ)為似然函數(shù)，即給定參數(shù)θ時(shí)數(shù)據(jù)D被觀測(cè)到的概率。

對(duì)于過程機(jī)制模型，先驗(yàn)分布通?；谏韺W(xué)觀測(cè)數(shù)據(jù)或理論推測(cè)的范圍(Van Oijenetal.， 2005)。對(duì)于參數(shù)數(shù)量較少且沒有生理學(xué)含義的經(jīng)驗(yàn)?zāi)Ｐ?，在首次校?參數(shù)擬合)時(shí)，采取無信息先驗(yàn)方式，參數(shù)聯(lián)合后驗(yàn)分布取決于似然函數(shù)，而之后的數(shù)據(jù)融合與模型更新過程中，將前一次擬合得到的參數(shù)聯(lián)合后驗(yàn)分布作為新數(shù)據(jù)加入時(shí)的先驗(yàn)。本研究中參數(shù)聯(lián)合后驗(yàn)分布不可解析，采用多維核密度估計(jì)方法得到近似分布，作為重新校正過程中的先驗(yàn)。

由于多種數(shù)據(jù)類型被應(yīng)用于貝葉斯校正過程，不同數(shù)據(jù)類型觀測(cè)誤差不同，因此需要對(duì)似然函數(shù)進(jìn)行分解。如對(duì)于樹高導(dǎo)向曲線[式(2)]，似然函數(shù)被分解為：

p(θ|D)∝p(DTP|θ)×p(DPP|θ)×

p(DSA|θ)×p(θ)。

(9)

式中：DTP、DPP和DSA分別表示角規(guī)臨時(shí)樣地(temporary plot，TP)、矩形固定樣地(permanent plot，PP)和解析木(stem analysis，SA)中的觀測(cè)數(shù)據(jù)，此時(shí)的p(θ|D)表示綜合考慮多種數(shù)據(jù)類型(multi-type/multi-source data)的參數(shù)后驗(yàn)概率。

考慮到林分調(diào)查中觀測(cè)誤差的不穩(wěn)定性，同時(shí)異常值很難有明確的判別標(biāo)準(zhǔn)，本研究采用重尾正態(tài)分布概率密度函數(shù)描述誤差容易受異常值影響而導(dǎo)致有偏的估計(jì)(Siviaetal.， 2006)：

(10)

式中：σi為噪聲尺度因子，用于描述第i次觀測(cè)隨機(jī)誤差的不確定性；Ri為第i次觀測(cè)值與對(duì)應(yīng)模型預(yù)測(cè)值的離差除以σi。

林業(yè)調(diào)查數(shù)據(jù)的另外一個(gè)重要特征就是異方差，即觀測(cè)值偏離均值的程度會(huì)受其他因素影響而非恒定。這種異方差特征一方面違反大部分經(jīng)典統(tǒng)計(jì)學(xué)假設(shè)，另一方面也導(dǎo)致模型預(yù)測(cè)中錯(cuò)誤的置信域估計(jì)。角規(guī)臨時(shí)樣地、矩形樣地和樹干解析信息需要依據(jù)噪聲尺度因子σi進(jìn)行融合，對(duì)于每種數(shù)據(jù)類型的每個(gè)響應(yīng)變量，均需要分別擬合其誤差分布。本研究假定觀測(cè)誤差的異方差性質(zhì)為：

(11)

1.4 聯(lián)立方程組估計(jì)下的參數(shù)聯(lián)合后驗(yàn)概率密度

可變密度全林模型理論上應(yīng)按1.2節(jié)介紹以特定順序單獨(dú)擬合，然而，由于全林模型內(nèi)各預(yù)估模型經(jīng)常會(huì)同時(shí)作為一個(gè)整體被使用，故本研究采用聯(lián)立方程組形式構(gòu)建模型，這就需要對(duì)似然函數(shù)進(jìn)行調(diào)整。以1990年采集的角規(guī)臨時(shí)樣地?cái)?shù)據(jù)擬合式(5)斷面積預(yù)測(cè)模型為例，在方程式單獨(dú)擬合時(shí)，其參數(shù)后驗(yàn)概率密度為：

p(c1,c2,c3,c4,gG,hG)。

(12)

p(α,β,b1,b2,c1,c2,c3,c4,d1,d2,d3,d4,glnN,hlnN,gH,

(13)

1.5 參數(shù)估計(jì)方法、收斂診斷與不確定性評(píng)估

由于參數(shù)聯(lián)合后驗(yàn)概率分布不可解析，故本研究采用集成Metropolis-Hastings算法的Markov Chain Monte Carlo(MCMC)技術(shù)(Hastings， 1970； Metropolisetal.， 1953)模擬計(jì)算參數(shù)分布，并運(yùn)用Gelman-Rubin收斂診斷技術(shù)(Brooksetal.， 1998； Gelmanetal.， 1992)確保MCMC收斂至正確的后驗(yàn)分布，該方法的優(yōu)勢(shì)在于整個(gè)過程既不涉及共軛計(jì)算，也不涉及條件概率和聯(lián)合后驗(yàn)分布解析式的推導(dǎo)計(jì)算(Van Oijenetal.， 2005)。參數(shù)聯(lián)合后驗(yàn)概率分布由MCMC的105次迭代生成，由于MCMC在開始階段未收斂至后驗(yàn)分布，故將鏈的前10%作為“燒入”(burn in)舍棄掉(Gilksetal.， 1996)。平行運(yùn)行3條MCMC鏈，計(jì)算其MPSRF(multivariate potential scale reduction factor)，該值較低時(shí)，代表不同起始值的子鏈最終獨(dú)立收斂至相同的后驗(yàn)分布。本研究設(shè)定MPSRF接受閾值為1.05，是一種相對(duì)嚴(yán)格的收斂診斷標(biāo)準(zhǔn)(Brooksetal.， 1998)。

不確定性評(píng)估采用Monte Carlo抽樣方法進(jìn)行計(jì)算。試驗(yàn)中截取MCMC記錄的后1/2，每隔10條記錄抽取1組參數(shù)，每條子鏈得到5 000組參數(shù)值，3條子鏈共計(jì)15 000組參數(shù)值。盡管參數(shù)聯(lián)合后驗(yàn)概率分布無法被解析，但無論是參數(shù)的邊際分布、相關(guān)性矩陣、標(biāo)準(zhǔn)差等信息，還是預(yù)測(cè)的貝葉斯可信區(qū)間，諸多不確定性指標(biāo)均可直接通過該樣本計(jì)算獲得。后續(xù)計(jì)算采用多維核密度估計(jì)方法得到近似分布，作為重新校正過程中的先驗(yàn)，該過程通過R工具包“BayesianTools”實(shí)現(xiàn)(Hartigetal.， 2019)。

采用Sobol指數(shù)量化分解各參數(shù)不確定性對(duì)最終預(yù)測(cè)不確定性的貢獻(xiàn)(Caribonietal.， 2007； Saltellietal.， 2008)。該方法可將不確定性分解到任意參數(shù)組，通過考慮參數(shù)相關(guān)性后單個(gè)參數(shù)所貢獻(xiàn)的總效應(yīng)(STi)對(duì)不確定性傳遞過程進(jìn)行評(píng)價(jià)，相應(yīng)計(jì)算公式為：

(14)

式中：x-i表示給定除xi以外的所有模型參數(shù)；Y表示模型輸出變量。

基于Monte Carlo輸出結(jié)果可用于計(jì)算上述指標(biāo)(Ioossetal.， 2020)。

各數(shù)據(jù)集本身規(guī)模不同，本研究采用自舉法分析數(shù)據(jù)規(guī)模對(duì)模型預(yù)測(cè)不確定性的影響。選取臨時(shí)樣地這一數(shù)據(jù)集，因?yàn)閿?shù)據(jù)集涵蓋林分條件，更為多樣化。依次設(shè)置樣本數(shù)量N的規(guī)模為10、20、30、…、200塊樣地，進(jìn)行以下計(jì)算： 1) 有放回的重抽樣，樣本數(shù)量為N； 2) 基于抽取到的樣地?cái)?shù)據(jù)對(duì)模型參數(shù)進(jìn)行擬合，采用擬合結(jié)果進(jìn)行預(yù)測(cè)； 3) 重復(fù)1、2過程1 000次。

2 結(jié)果與分析

2.1 不同時(shí)期調(diào)查信息的動(dòng)態(tài)融合

區(qū)別于經(jīng)典回歸參數(shù)符合正態(tài)分布的假設(shè)，本研究貝葉斯方法獲得的參數(shù)分布是不可解析的。盡管參數(shù)后驗(yàn)分布并不嚴(yán)格對(duì)稱，但除個(gè)別參數(shù)呈明顯偏態(tài)分布外，整體而言均值與最大后驗(yàn)概率估計(jì)值(maximum a posterior estimates,MAP)相差極小(表3、4)。圖2以參數(shù)b1為例，展示了隨著新數(shù)據(jù)加入?yún)?shù)邊際分布不斷更新的過程。樹高生長(zhǎng)的上漸近線參數(shù)b1(樹高生長(zhǎng)極大值)隨著新數(shù)據(jù)加入不斷升高，不確定性逐步降低；同時(shí)，參數(shù)b1與b2存在較高相關(guān)性，且其聯(lián)合分布表現(xiàn)出明顯正向關(guān)聯(lián)(圖3)，隨著數(shù)據(jù)更新，聯(lián)合分布峰值逐漸升高，參數(shù)不確定性也隨之下降。

表3 參數(shù)循環(huán)更新過程中后驗(yàn)概率密度分布的變化①

表4 基于多源數(shù)據(jù)的參數(shù)后驗(yàn)概率分布

圖2 馬爾科夫蒙特卡洛迭代及參數(shù)b1后驗(yàn)概率密度的變化過程

圖3 循環(huán)更新中參數(shù)b1與b2的聯(lián)合后驗(yàn)概率分布

圖4 林分平均高導(dǎo)向曲線在20、40、60年時(shí)的樹高預(yù)測(cè)分布

隨著參數(shù)概率分布改變，模型預(yù)測(cè)結(jié)果也隨著新數(shù)據(jù)加入不斷變化。貝葉斯框架下，參數(shù)不確定性最終會(huì)對(duì)預(yù)測(cè)的不確定性產(chǎn)生影響。基于1990年數(shù)據(jù)擬合模型預(yù)測(cè)的不確定性高于之后2次校正得到的模型(圖4、5)。對(duì)于95%貝葉斯可信區(qū)間而言，基于1990年數(shù)據(jù)擬合的模型整個(gè)預(yù)測(cè)周期內(nèi)樹高導(dǎo)向曲線的平均不確定區(qū)間為2.3 m(圖5a)，經(jīng)2005年數(shù)據(jù)一次校正后不確定性下降至1.9 m(圖5b)，而經(jīng)2012年數(shù)據(jù)二次校正后不確定性下降至1.1 m(圖5c)。同時(shí)，首次參數(shù)化時(shí)的模型在20年樹高導(dǎo)向曲線估計(jì)中傾向于給出更高的幼齡林樹高預(yù)測(cè)(圖4a)和更低的成過熟林樹高預(yù)測(cè)(圖4c)。

林分動(dòng)態(tài)預(yù)測(cè)過程中，數(shù)據(jù)與模型的關(guān)系是相輔相成的： 一方面，模型的開發(fā)和參數(shù)化需要基于數(shù)據(jù)； 另一方面，模型的表現(xiàn)也可直接反映數(shù)據(jù)特征信息。在數(shù)據(jù)信息融合與模型學(xué)習(xí)關(guān)系構(gòu)建時(shí)，模型預(yù)測(cè)的最大概率曲線在幼齡林和成過熟林模擬上不斷進(jìn)行細(xì)微調(diào)整(圖5)，模型對(duì)數(shù)據(jù)信息的不斷吸收，更顯著的作用是不斷更新模型預(yù)測(cè)的不確定性?；?990年數(shù)據(jù)進(jìn)行模型擬合后，林分平均高、斷面積和平均胸徑的預(yù)測(cè)不確定性區(qū)間均在25～30年最低，而在其他年齡階段的不確定性則相對(duì)更高(圖5a、5d、5g)。經(jīng)2005和2012年數(shù)據(jù)校正后，各齡階的不確定性預(yù)測(cè)逐漸平衡。以林分平均胸徑預(yù)測(cè)為例， 1990年數(shù)據(jù)對(duì)應(yīng)最低不確定性出現(xiàn)在林齡25年(圖5d)，經(jīng)2005年數(shù)據(jù)一次校正，最低不確定性出現(xiàn)在林齡34年(圖5e)，而經(jīng)2012年數(shù)據(jù)二次校正，最低不確定性出現(xiàn)在林齡48年(圖5f)。

圖5 林分平均高(導(dǎo)向曲線)、平均胸徑和斷面積預(yù)測(cè)的更新變化過程(SCI=11, SDI=600)

2.2 不同數(shù)據(jù)類型對(duì)林分模擬的影響

角規(guī)繞測(cè)臨時(shí)樣地、矩形固定樣地、解析木數(shù)據(jù)在抽樣原理、測(cè)量誤差等方面的區(qū)別，導(dǎo)致模型參數(shù)擬合也存在差異。綜合多種數(shù)據(jù)類型的模型校正基于多個(gè)似然函數(shù)相乘來描述數(shù)據(jù)與模型的關(guān)系，即同時(shí)考慮數(shù)據(jù)量和數(shù)據(jù)準(zhǔn)確性對(duì)似然的影響。參數(shù)估計(jì)時(shí)，對(duì)于不同數(shù)據(jù)集中觀測(cè)誤差大小、變異幅度及其方差異質(zhì)性的假設(shè)，會(huì)基于模型與數(shù)據(jù)的匹配程度進(jìn)行自動(dòng)調(diào)整，因此并不需要預(yù)設(shè)各數(shù)據(jù)集在總體似然函數(shù)中所占比重。單從邊際分布考慮，綜合多源數(shù)據(jù)的模型參數(shù)大多是對(duì)角規(guī)臨時(shí)樣地、固定樣地和解析木對(duì)應(yīng)模型參數(shù)的折中平均，且整體上與角規(guī)臨時(shí)樣地的擬合結(jié)果最為接近(圖6e、6h)。相比樣地調(diào)查，解析木數(shù)據(jù)傾向于給出更高的參數(shù)b1估計(jì)，即更大的樹高極大值(圖6g)。從聯(lián)合后驗(yàn)分布角度來看，樣地調(diào)查數(shù)據(jù)給出的參數(shù)估計(jì)比解析木參數(shù)有著更強(qiáng)烈的相關(guān)性(圖7)。在樹高導(dǎo)向曲線擬合中，林齡20年時(shí)幾種類型數(shù)據(jù)給出的預(yù)測(cè)結(jié)果是相似的(圖8a)，林齡40年(圖8b)和60年(圖8c)時(shí)基于解析木的預(yù)測(cè)則明顯高于其他數(shù)據(jù)組。

圖6 不同數(shù)據(jù)類型中馬爾科夫蒙特卡洛迭代及參數(shù)b1后驗(yàn)概率密度的變化過程

相比之下，基于多源數(shù)據(jù)構(gòu)造的全林模型，在林分平均高、平均胸徑和斷面積預(yù)測(cè)中均表現(xiàn)出最低不確定性。固定樣地?cái)?shù)量遠(yuǎn)低于臨時(shí)樣地，調(diào)查年份均為2012—2015年，且多為坡度較緩的林分，這使得其在林分?jǐn)嗝娣e生長(zhǎng)預(yù)測(cè)和不確定性估計(jì)上與臨時(shí)樣地明顯不同。基于固定樣地的預(yù)測(cè)結(jié)果表明，林分?jǐn)嗝娣e增長(zhǎng)在林齡35年后并不會(huì)大幅減緩(圖9e)，但這種預(yù)測(cè)也伴隨著極高不確定性(±30.1%)。各參數(shù)對(duì)模型預(yù)測(cè)不確定性的貢獻(xiàn)并不相同，斷面積預(yù)測(cè)模型參數(shù)c4與胸徑預(yù)測(cè)模型參數(shù)d4在式(5)和式(6)中傳遞的不確定性明顯低于其他相關(guān)參數(shù)(表5)。同時(shí)，由于參數(shù)聯(lián)合后驗(yàn)概率分布中參數(shù)相關(guān)性較高，各參數(shù)Sobol總效應(yīng)指數(shù)加和后的數(shù)值也會(huì)遠(yuǎn)高于1。

圖8 基于不同數(shù)據(jù)類型林分平均高導(dǎo)向曲線在20、40、60年時(shí)的樹高預(yù)測(cè)分布

圖9 基于不同數(shù)據(jù)類型的林分平均胸徑和斷面積預(yù)測(cè)的更新變化過程(SCI=11, SDI=600)

表5 基于方差的不確定性量化分解

模型預(yù)測(cè)的不確定性與建模數(shù)據(jù)本身規(guī)模直接相關(guān)。以林分?jǐn)嗝娣e為例(圖10)，當(dāng)林齡為20、40、60年時(shí)，若抽選10塊樣地建模，1 000次重復(fù)后模型預(yù)測(cè)值波動(dòng)標(biāo)準(zhǔn)差為1.77、1.40、1.90，而當(dāng)樣地?cái)?shù)量增加至200塊時(shí)，標(biāo)準(zhǔn)差則下降至0.44、0.27、0.30，這也直接解釋了圖9e中的較高不確定性實(shí)際是由固定樣地?cái)?shù)據(jù)量較小造成的。隨著數(shù)據(jù)量增加，不確定性的降低速度逐步減緩，從10個(gè)樣本增加至20個(gè)樣本時(shí)，標(biāo)準(zhǔn)差減小14%； 從190個(gè)樣本增加至200個(gè)樣本時(shí)，標(biāo)準(zhǔn)差減小2%。

圖10 隨數(shù)據(jù)規(guī)模增加林分?jǐn)嗝娣e模型預(yù)測(cè)不確定性的變化過程(SCI=11, SDI=600)

3 討論

3.1 數(shù)據(jù)類型與抽樣方法對(duì)林分模擬的影響

本研究涉及的幾種數(shù)據(jù)類型各有其局限性，多源數(shù)據(jù)可避免某一數(shù)據(jù)類型缺陷過于嚴(yán)重傳遞至模型預(yù)測(cè)中。出于科研目的繼而提供精準(zhǔn)信息的固定樣地長(zhǎng)期觀測(cè)數(shù)據(jù)有著極高的經(jīng)濟(jì)、人力和時(shí)間采集成本，通?？色@取數(shù)據(jù)量較少，很難涵蓋林分條件的變異幅度。臨時(shí)樣地成本低、數(shù)據(jù)量大，但與室內(nèi)控制試驗(yàn)不同，野外調(diào)查臨時(shí)樣地?cái)?shù)據(jù)無法滿足抽樣均衡性，即不能保證各年齡階段樣地的立地質(zhì)量分布和變異是否相似，容易在某個(gè)齡階出現(xiàn)有偏預(yù)測(cè)。解析木數(shù)據(jù)對(duì)林分整體狀況的代表性很難保證。以樹高生長(zhǎng)為例，立地質(zhì)量對(duì)林分生長(zhǎng)動(dòng)態(tài)具有直接影響，這意味著多源數(shù)據(jù)中的變異有特定比例可以被立地差異所解釋。在林分平均高、平均胸徑、斷面積模型中，基于林齡和林分平均高得到的地位級(jí)指數(shù)被用作解釋性變量，以描述不同立地條件下的林分生長(zhǎng)收獲差異。為減小立地評(píng)價(jià)的系統(tǒng)性偏差，樹高導(dǎo)向曲線對(duì)抽樣均衡性有著極高要求，即各年齡階段樣地的立地質(zhì)量分布和變異應(yīng)是相似的，且單個(gè)樣地的立地質(zhì)量在較長(zhǎng)時(shí)期內(nèi)保持一致。這些要求在多數(shù)情況下往往無法被完全滿足。在樹高導(dǎo)向曲線擬合中，解析木與樣地調(diào)查數(shù)據(jù)在樹木生長(zhǎng)趨勢(shì)描述上有著明顯區(qū)別。理論上，解析木測(cè)量值比大多數(shù)樹高調(diào)查方法都更加準(zhǔn)確，且能夠提供完整時(shí)間序列上的樹高動(dòng)態(tài)信息，因此在樹高導(dǎo)向曲線中，本研究增加部分中等木和亞優(yōu)勢(shì)木的樹干解析信息，以輔助解決幼齡林階段數(shù)據(jù)不足的問題，同時(shí)避免各年齡階段樣地立地質(zhì)量分布差異導(dǎo)致的系統(tǒng)性偏差。本研究涉及林分為相對(duì)同齡而非嚴(yán)格同齡林，林木間存在一定程度年齡差異，作為樣本的中等木或亞優(yōu)勢(shì)木可能在生命周期中相當(dāng)長(zhǎng)的時(shí)間內(nèi)處于被壓迫狀態(tài)，而在樹高達(dá)到一定高度后從競(jìng)爭(zhēng)壓力中恢復(fù)過來(Hannetal.， 2002)。另外，解析木選取常常基于樹木的干形和活力特征，這些都使得其樹高生長(zhǎng)趨勢(shì)無法代表林分平均水平，且在成過熟林階段產(chǎn)生一定程度的高估(圖7g、7c)。當(dāng)多源數(shù)據(jù)中不同數(shù)據(jù)集跨越較長(zhǎng)時(shí)間尺度時(shí)，立地評(píng)價(jià)穩(wěn)定性也會(huì)受多種因素影響，包括樹木子代個(gè)體基因型的改變(Monserudetal.， 1990)、氣候條件的改變(Bontempsetal.， 2009)、特定林分經(jīng)營措施(Fox， 2000； Skovsgaardetal.， 2008)等。

盡管角規(guī)臨時(shí)樣地和矩形固定樣地均按照森林經(jīng)理調(diào)查時(shí)分層抽樣與隨機(jī)抽樣相結(jié)合的準(zhǔn)則設(shè)置，但在實(shí)際操作中，矩形固定樣地常常會(huì)被主觀建立在小班內(nèi)坡度相對(duì)平緩的區(qū)域。相比之下，由于操作的便捷性，角規(guī)臨時(shí)樣地可以實(shí)現(xiàn)在平緩的下坡位和陡峭的山脊或陡坡均衡抽樣，而山脊或陡坡的林分常常因土壤瘠薄等因素生產(chǎn)力較低(Schrothetal.， 2003)。角規(guī)調(diào)查的缺陷是林下灌木、幼樹和地形等形成的復(fù)雜結(jié)構(gòu)會(huì)造成繞測(cè)時(shí)視線遮擋(Van Laaretal.， 2007)。這些因素均會(huì)使得基于矩形固定樣地的斷面積生長(zhǎng)量預(yù)測(cè)高于角規(guī)臨時(shí)樣地推算的結(jié)果(圖9e、8f)。

3.2 不確定性的量化

復(fù)雜非線性模型回歸分析通常被認(rèn)為是優(yōu)化問題。相比僅僅給出參數(shù)的最優(yōu)估計(jì)值，貝葉斯方法通過MCMC抽樣技術(shù)獲得參數(shù)聯(lián)合后驗(yàn)分布，使得模型中的信息可以更好體現(xiàn)在模型中。但即使是最優(yōu)參數(shù)組合相同或相似，不同數(shù)據(jù)所給出的參數(shù)概率分布也可能存在很大區(qū)別，數(shù)據(jù)的異質(zhì)性、觀測(cè)誤差、樣本數(shù)量等會(huì)反映到參數(shù)聯(lián)合后驗(yàn)分布中，并最終影響模型預(yù)測(cè)的不確定性。對(duì)比圖5中2005和2012年數(shù)據(jù)組，或圖7、9中臨時(shí)樣地和固定樣地，最高概率曲線模擬的林分生長(zhǎng)趨勢(shì)非常相似，但數(shù)據(jù)中的信息(如不同齡階時(shí)預(yù)測(cè)的可信區(qū)間)卻并不相同?；?990年數(shù)據(jù)的預(yù)測(cè)在林分成過熟林階段有著很高不確定性，這是因?yàn)樵撾A段的觀測(cè)較少且離散較高，而2005和2012年數(shù)據(jù)加入補(bǔ)充了成過熟林階段的林分信息(圖5)。大部分模型結(jié)果在林齡15～20年時(shí)均有著較高不確定性，這是因幼齡林階段數(shù)據(jù)缺失造成的，尤其在林分?jǐn)嗝娣e預(yù)測(cè)中最為明顯(圖9d、9e、9f)。整體而言，在不增加林分異質(zhì)性信息的前提下，增大數(shù)據(jù)量是降低不確定性最直接的方法。相比林分平均高和平均胸徑，林分?jǐn)嗝娣e預(yù)測(cè)呈現(xiàn)出最高不確定性，且在數(shù)據(jù)量較低時(shí)差異尤為明顯(圖9e)。這說明對(duì)于林分?jǐn)嗝娣e而言，降低不確定性所需的數(shù)據(jù)量和抽樣調(diào)查的系統(tǒng)性都更高。影響林分動(dòng)態(tài)的其他因素，如人為干擾、林型特征等，需要被更多地考慮到斷面積模型構(gòu)建中，無論是2.1節(jié)的多期數(shù)據(jù)更新，還是2.2節(jié)的多種數(shù)據(jù)類型比較，綜合考慮全部信息的模型無論是參數(shù)還是預(yù)測(cè)均具有最高精度(最低不確定性)。

不確定性的解析式分解遵從Tylor展開式原理(Dietze， 2017)，貝葉斯方法以概率分布方式來描述參數(shù)和模型預(yù)測(cè)，由于本研究涉及的參數(shù)后驗(yàn)概率分布不可解析，不確定性評(píng)估以Monte Carlo抽樣方法進(jìn)行計(jì)算。需要強(qiáng)調(diào)的是，本研究涉及的數(shù)據(jù)和模型都較為簡(jiǎn)單，因而只考慮了參數(shù)不確定性對(duì)預(yù)測(cè)的影響，并沒有涉及數(shù)據(jù)、模型結(jié)構(gòu)、模型鏈接等部分的不確定性及其傳遞關(guān)系。在復(fù)雜模型的組合應(yīng)用中，不確定性的傳遞與解構(gòu)則是理解模型預(yù)測(cè)結(jié)果的關(guān)鍵性信息(Lebaueretal.， 2013)。

3.3 先驗(yàn)的設(shè)計(jì)

不同于過程機(jī)制模型，經(jīng)驗(yàn)?zāi)Ｐ蛥?shù)本身沒有實(shí)際生理學(xué)意義，很難給出準(zhǔn)確合理的先驗(yàn)設(shè)定。本研究沒有探討何種先驗(yàn)是“好的先驗(yàn)”，實(shí)際上討論的先驗(yàn)都是基于以往數(shù)據(jù)信息，即基于往期數(shù)據(jù)得到的后驗(yàn)作為模型更新時(shí)的先驗(yàn)。以1990年數(shù)據(jù)進(jìn)行首次擬合時(shí)，采用無信息先驗(yàn)，即等同于極大似然：

p(θ1990|D1990)∝p(D1990|θ1990)。

2005年數(shù)據(jù)的模型校正，以1990年模型后驗(yàn)作為先驗(yàn)，即：

p(θ2005|D2005)∝p(D2005|θ2005)×p(θ2005)∝

p(D2005|θ2005)×p(θ1990|D1990)∝

p(D2005|θ2005)×p(D1990|θ1990)。

2012年數(shù)據(jù)的模型校正，以2005年一次校正模型后驗(yàn)作為先驗(yàn)，即：

p(θ2012|D2012)∝p(D2012|θ2012)×p(θ2012)∝

p(D2012|θ2012)×p(θ2005|D2005)∝

p(D2012|θ2012)×p(D2005|θ2005)×p(D1990|θ1990)。

當(dāng)模型完成基于2012年數(shù)據(jù)的二次校正后，參數(shù)后驗(yàn)同時(shí)包含1990、2005和2012年的數(shù)據(jù)信息。在這一模型循環(huán)更新框架下，如果有明確證據(jù)表明1990年數(shù)據(jù)不該作為2005年數(shù)據(jù)的先驗(yàn)，那么就意味著這2期數(shù)據(jù)中有1期是錯(cuò)誤的且應(yīng)直接將其舍棄，但通常情況下，隨著數(shù)據(jù)積累，模型可靠性會(huì)逐步升高，新數(shù)據(jù)涵蓋了新的或更復(fù)雜的林分狀態(tài)信息。2005年數(shù)據(jù)校正和2012年數(shù)據(jù)校正模型變化幅度更小，則說明新增加數(shù)據(jù)提供的額外信息很少，不需要再投入過高成本對(duì)這些變化微小的變量繼續(xù)收集數(shù)據(jù)。本研究采用重尾正態(tài)分布降低異常值對(duì)模型擬合的影響，各期數(shù)據(jù)似然結(jié)構(gòu)相互獨(dú)立，以避免精確數(shù)據(jù)被粗糙數(shù)據(jù)所稀釋。貝葉斯方法被廣泛應(yīng)用于工程與信息科學(xué)的重要原因之一在于其高效的模型學(xué)習(xí)框架，這種高效性體現(xiàn)在當(dāng)獲取到新一期數(shù)據(jù)對(duì)模型參數(shù)進(jìn)行更新時(shí)，并不需要對(duì)新獲取數(shù)據(jù)的規(guī)模有任何特定要求(甚至1個(gè)觀測(cè)點(diǎn)即可)，模型更新也不需要重新提取整理往期數(shù)據(jù)，因?yàn)橥跀?shù)據(jù)信息已包含在先驗(yàn)中。

3.4 數(shù)據(jù)-模型關(guān)系

森林生態(tài)系統(tǒng)研究是一個(gè)漸進(jìn)的過程，模型在該過程中至關(guān)重要。將已知信息作為先驗(yàn)并與未來可能獲得的數(shù)據(jù)相結(jié)合，從而得到后驗(yàn)和預(yù)測(cè)分布(Westetal.， 1997)，隨著數(shù)據(jù)積累逐漸更新對(duì)林分動(dòng)態(tài)的理解，這一過程即數(shù)據(jù)-模型融合過程，其目的是使模型不斷吸收新的信息(Clark， 2007)。

模型的準(zhǔn)確性和適用范圍取決于建模數(shù)據(jù)本身。理想情況下，模型構(gòu)建和數(shù)據(jù)收集是不斷迭代的過程，模型設(shè)計(jì)決定數(shù)據(jù)需求，進(jìn)而決定外業(yè)調(diào)查內(nèi)容和方式； 然而森林生長(zhǎng)周期漫長(zhǎng)，獲取全部必要信息亦需漫長(zhǎng)觀測(cè)周期，同時(shí)，選用不同數(shù)據(jù)也會(huì)對(duì)模型擬合造成明顯差異(Raulieretal.， 2003)，數(shù)據(jù)的類型、尺度、質(zhì)量和規(guī)模等因素對(duì)模型構(gòu)建與基于模型衍生的推論都有著重要影響。森林資源連續(xù)調(diào)查中，每一次調(diào)查分析結(jié)果都是下一次調(diào)查分析的最合理先驗(yàn)(張雄清等， 2014)。森林生長(zhǎng)收獲模型開發(fā)常常需要收集多個(gè)數(shù)據(jù)集信息，不同來源的數(shù)據(jù)在抽樣設(shè)計(jì)、觀測(cè)誤差和觀測(cè)數(shù)量等方面都會(huì)有所區(qū)別，而這些信息很難反映在模型的一組最優(yōu)參數(shù)上。隨著信息數(shù)據(jù)加入，模型改變也并不局限于最優(yōu)參數(shù)組合的調(diào)整。綜合考慮各數(shù)據(jù)組的誤差分布、異方差等問題，需要模型開發(fā)者針對(duì)不同數(shù)據(jù)集采用獨(dú)立的似然結(jié)構(gòu)[式(9)]，不是簡(jiǎn)單地將數(shù)據(jù)進(jìn)行混合。而在數(shù)據(jù)不斷更新、模型學(xué)習(xí)并升級(jí)過程中，參數(shù)聯(lián)合后驗(yàn)分布不斷改變(圖2和3)，進(jìn)而導(dǎo)致預(yù)測(cè)的不確定性也在不斷發(fā)生變化(圖4和5)。在這種貝葉斯數(shù)據(jù)-模型框架中，數(shù)據(jù)與模型之間以概率分布形式融合為信息交互的循環(huán)體系。

本研究分析是基于貝葉斯框架下全林模型以聯(lián)立方程組方式進(jìn)行參數(shù)估計(jì)的結(jié)果[將式(1)和式(3)代入式(5)和式(6)后進(jìn)行聯(lián)立]，由此導(dǎo)致模型誤差的傳遞，即地位級(jí)指數(shù)模型和林分密度指數(shù)模型存在誤差，之后又將誤差引入林分?jǐn)嗝娣e和平均胸徑模型中。貝葉斯框架下聯(lián)立方程組比單獨(dú)擬合有著更為復(fù)雜的似然結(jié)構(gòu)[式(12)、(13)]，相應(yīng)地，模型開發(fā)者可以根據(jù)自身對(duì)數(shù)據(jù)和模型的理解來調(diào)整誤差分布和異方差假設(shè)，在擬合生長(zhǎng)收獲模型參數(shù)的同時(shí)獲得描述誤差分布的參數(shù)估計(jì)值[式(10)、(11)]。聯(lián)立方程組模型中常見處理變量誤差和異方差的方式之一是進(jìn)行對(duì)數(shù)轉(zhuǎn)換(洪玲霞等， 2012)，但對(duì)變量進(jìn)行對(duì)數(shù)轉(zhuǎn)化本身并沒有生態(tài)學(xué)或生理學(xué)解釋，同時(shí)，基于模型輸出結(jié)果對(duì)變量進(jìn)行逆向變換后預(yù)測(cè)結(jié)果也常常是有偏的(Dietze， 2017)。Zellner(1962)的似乎不相關(guān)聯(lián)立估計(jì)在考慮不同方程誤差項(xiàng)相關(guān)性的同時(shí)，也兼顧了異方差性問題，但較難與本研究的貝葉斯框架兼容。貝葉斯框架下進(jìn)行似乎不相關(guān)回歸和聯(lián)立方程組回歸通常會(huì)圍繞基于協(xié)方差矩陣對(duì)聯(lián)立后參數(shù)條件概率與聯(lián)合后驗(yàn)的推導(dǎo)計(jì)算(Andoetal.， 2010)，本研究完全依賴于集成Metropolis-Hastings算法的MCMC技術(shù)來實(shí)現(xiàn)對(duì)聯(lián)合后驗(yàn)分布的抽樣和估計(jì)，從而最終實(shí)現(xiàn)將式(1)和式(3)代入式(5)和式(6)后所獲得聯(lián)立方程組的參數(shù)估計(jì)過程。考慮誤差傳遞的另一個(gè)重要方法是度量誤差效應(yīng)模型(唐守正等， 1996； 1998)，該模型在貝葉斯框架下亦可實(shí)現(xiàn)(Clark， 2007； Clarketal.， 2007)，不僅地位級(jí)指數(shù)和立地指數(shù)評(píng)估的誤差在胸徑和斷面積模型中傳遞，而且森林調(diào)查中的胸徑、樹高、斷面積觀測(cè)誤差也可以被考慮。本研究所涉及的參數(shù)不確定性及其對(duì)預(yù)測(cè)不確定性的影響，僅僅只是不確定性或誤差分析的一部分內(nèi)容，同時(shí)受限于數(shù)據(jù)規(guī)模和質(zhì)量，在測(cè)試數(shù)據(jù)代表性和描述模型開發(fā)中不確定性傳遞上仍有諸多不足。

本研究采用貝葉斯數(shù)據(jù)-模型框架，出發(fā)點(diǎn)是為了兼顧不同數(shù)據(jù)集所具有的不同抽樣和觀測(cè)誤差，進(jìn)而對(duì)不確定性進(jìn)行量化。貝葉斯框架將模型校正和預(yù)測(cè)中的所有要素以概率分布方式進(jìn)行處理(Clark， 2007； Dietze， 2017)，能夠方便實(shí)現(xiàn)以上目標(biāo)。但需要強(qiáng)調(diào)的是，貝葉斯并不是實(shí)現(xiàn)這些目標(biāo)的唯一方法，多源數(shù)據(jù)的觀測(cè)誤差處理可通過若干似然函數(shù)相乘或權(quán)重法設(shè)置目標(biāo)函數(shù)來實(shí)現(xiàn)(Marleretal.， 2010)，抽樣方法的系統(tǒng)性誤差也可通過多層級(jí)的混合效應(yīng)模型來評(píng)估(Bijleveldetal.， 1998； Wareetal.， 1996)，不確定性的量化則可依賴于自舉法來實(shí)現(xiàn)(Efron, 1979)。

4 結(jié)論

本研究以全林模型更新為例，描述貝葉斯框架下建模數(shù)據(jù)累積與生長(zhǎng)模型不確定性循環(huán)學(xué)習(xí)的信息融合過程。結(jié)果表明，隨著森林調(diào)查數(shù)據(jù)不斷積累，模型參數(shù)和預(yù)測(cè)不確定性逐漸下降。建模數(shù)據(jù)信息中存在的缺陷可被量化并以不確定性方式直接反映在模型參數(shù)和預(yù)測(cè)中，進(jìn)而為下一步模型改進(jìn)與數(shù)據(jù)采集指明方向。多源數(shù)據(jù)的使用弱化了單一數(shù)據(jù)源中非均衡抽樣和異常觀測(cè)誤差對(duì)模型的影響，可有效提高林分動(dòng)態(tài)預(yù)測(cè)質(zhì)量。這種數(shù)據(jù)-模型循環(huán)更新和兼容多源數(shù)據(jù)的貝葉斯框架，同樣也適用于森林建模研究的其他諸多方向，如森林生產(chǎn)力、森林立地質(zhì)量以及進(jìn)界、自疏伐等森林動(dòng)態(tài)與不確定性的研究和應(yīng)用。