油松幼樹樹高生長(zhǎng)預(yù)測(cè)的不確定性貝葉斯分析*

王 彬 田相林 曹田健

(西北農(nóng)林科技大學(xué)林學(xué)院 生態(tài)仿真優(yōu)化實(shí)驗(yàn)室 楊凌 712100)

森林的更新潛力源于幼樹的數(shù)量和生長(zhǎng)狀況，幼樹通過更新、枯損、生長(zhǎng)和進(jìn)界等一系列過程進(jìn)入喬木層改變林分結(jié)構(gòu)，參與森林生態(tài)系統(tǒng)的動(dòng)態(tài)演替。構(gòu)建幼樹樹高生長(zhǎng)模型是了解幼樹生長(zhǎng)動(dòng)態(tài)的有效途徑之一(Boisvenueetal., 2004)，但模型預(yù)測(cè)過程是在多種不確定性前提條件下進(jìn)行的，如幼樹當(dāng)前的個(gè)體高度、立地條件、地理位置、競(jìng)爭(zhēng)等(Cobbetal.， 1993； Oliveretal.， 1996); 同時(shí)，在幼樹樹高建模中，數(shù)據(jù)的測(cè)量、參數(shù)的選取以及模型誤差等也均存在不確定性，影響模型的精準(zhǔn)性和預(yù)測(cè)能力。因此，弄清楚模型預(yù)測(cè)中不確定性來源的類型及其對(duì)輸出的影響，對(duì)解釋模型預(yù)測(cè)結(jié)果具有重要意義。

不確定性，即不精確性、模糊性、不明確性等(Shi， 2009)，模型的不確定性來源一般可歸納為模型形式及結(jié)構(gòu)的不確定性、輸入變量的不確定性、模型誤差的不確定性和模型參數(shù)的不確定性4方面(Lessardetal.， 2001； Kitaharaetal.， 2009； Stahletal.， 2010)。不確定性的來源和量化分析是模型開發(fā)和應(yīng)用的重要環(huán)節(jié)。粗糙的森林調(diào)查數(shù)據(jù)通常會(huì)使模型中參數(shù)的解釋效應(yīng)偏低，參數(shù)再將這種不確定性傳遞給模型輸出，會(huì)進(jìn)一步擴(kuò)大預(yù)測(cè)的不確定性，降低模型預(yù)測(cè)精度。不確定性量化評(píng)估是衡量模型統(tǒng)計(jì)推斷意義的重要依據(jù)，缺少不確定性的估計(jì)，就無法真正了解模型預(yù)測(cè)和觀測(cè)數(shù)據(jù)之間的內(nèi)在聯(lián)系和規(guī)律 (LeBaueretal.， 2013)，繼而限制模型的應(yīng)用能力(Clarketal.， 2001)。分析模型在不同條件下預(yù)測(cè)結(jié)果的差異能夠幫助確定模型的合理應(yīng)用范圍，明確模型的改進(jìn)目標(biāo)(Williamsetal.， 2009)。

由于數(shù)據(jù)獲取的難易程度和林分生長(zhǎng)收獲模擬的需要，多數(shù)樹高生長(zhǎng)模型的研究對(duì)象為大樹(DBH≥5 cm)(Curtis， 1967； Monserudetal.， 1977； Dolph， 1988； Lappietal.， 1988； Huangetal.， 1994)，但這類模型很難準(zhǔn)確預(yù)測(cè)幼樹樹高生長(zhǎng)(Boisvenueetal., 2004)。曾有學(xué)者以線性回歸和非線性回歸為主要模型形式對(duì)幼樹樹高生長(zhǎng)進(jìn)行預(yù)測(cè)(Golseretal.， 1997； Hasenaueretal.， 2002； Nighaetal., 1999； 國(guó)紅等， 2009)。為提高模型適用性，Wykoff等(1982)以當(dāng)前樹高代替年齡，同時(shí)考慮立地因子(坡度、坡向)、競(jìng)爭(zhēng)因子(樹冠競(jìng)爭(zhēng)因子、大于對(duì)象木的斷面積之和)對(duì)幼樹樹高的影響，并將生境或生長(zhǎng)區(qū)域作為啞元變量添加到模型中，結(jié)果發(fā)現(xiàn)增加立地因子和競(jìng)爭(zhēng)因子可提高模型預(yù)測(cè)精度。Boisvenue等(2004)基于坡度、坡向、樹冠競(jìng)爭(zhēng)因子等變量，比較幾種線性和非線性模型對(duì)英屬哥倫比亞東南區(qū)域混交林分中幼樹5年的樹高生長(zhǎng)預(yù)測(cè)結(jié)果，發(fā)現(xiàn)非線性模型的預(yù)測(cè)精度相對(duì)較高。然而，目前對(duì)幼樹樹高生長(zhǎng)模型的研究仍停留在模型變量和模型形式表達(dá)方面，缺乏模型預(yù)測(cè)中不確定性的研究，對(duì)于幼樹樹高生長(zhǎng)模型輸出與數(shù)據(jù)之間的關(guān)系和差異尚待探索。

近年來，森林動(dòng)態(tài)預(yù)測(cè)中的不確定性逐漸得到研究者重視。如在森林生物量估算中，傅煜等(2015)、秦立厚等(2017)通過模型殘差變異和一階泰勒級(jí)數(shù)對(duì)模型的不確定性進(jìn)行量化，發(fā)現(xiàn)杉木(Cunninghamialanceolata)生物量模型不確定性以參數(shù)的不確定性對(duì)預(yù)測(cè)結(jié)果影響最大，其次是模型殘差變異的不確定性。趙菡等(2017)采用蒙特卡洛方法，通過對(duì)不同立地條件下馬尾松(Pinusmassoniana)的幾種生物量模型進(jìn)行評(píng)估，得出帶有樹高因子的生物量模型誤差相對(duì)較小。也有學(xué)者運(yùn)用全局敏感性分析(GSA)方法進(jìn)行模型參數(shù)的不確定性分析，主要有Morris篩選法、Sobol方法、傅里葉幅度敏感性檢驗(yàn)法(FAST)、GLUE方法等(Bevenetal.， 1992)，其中Sobol方法基于方差分解原理，可用于非線性、非單調(diào)的數(shù)學(xué)模型，結(jié)果穩(wěn)健可靠，且能給出參數(shù)的各階敏感性系數(shù)，因而得到廣泛應(yīng)用。Sobol方法的基礎(chǔ)是蒙特卡洛方法，通過蒙特卡洛模擬隨機(jī)抽樣計(jì)算敏感性系數(shù)，但該方法采用的準(zhǔn)隨機(jī)數(shù)字序列由系統(tǒng)軟件直接產(chǎn)生，對(duì)不確定性估計(jì)有所偏差； 而貝葉斯框架與全局敏感性分析結(jié)合比一般的參數(shù)搜索方法更具優(yōu)勢(shì)，因?yàn)槠淠軌蛱峁┡c試驗(yàn)設(shè)計(jì)相關(guān)的參數(shù)空間，預(yù)測(cè)結(jié)果更準(zhǔn)確(Erikssonetal., 2019)。貝葉斯方法能夠采用馬爾可夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)抽樣方法估計(jì)出每個(gè)參數(shù)的后驗(yàn)分布空間，然后從參數(shù)的后驗(yàn)分布空間直接抽樣量化模型參數(shù)的不確定性(Ellison， 2004； 李向陽(yáng)， 2005； Van Oijenetal.， 2005； Zhangetal.， 2009； 2014； 張雄清等， 2014； Simoetal.， 2017)，參數(shù)的不確定性再傳遞給模型輸出，最終實(shí)現(xiàn)模型預(yù)測(cè)不確定性的量化。這種方法將抽樣與數(shù)據(jù)緊密結(jié)合，約束了模型的參數(shù)空間，對(duì)參數(shù)及參數(shù)組合作用與模型輸出不確定性的關(guān)系進(jìn)行深入分析，從而確定不同參數(shù)對(duì)模型輸出不確定性的影響(Saltellietal.， 2018)。Tang等(2009)在一個(gè)現(xiàn)有的陸地生態(tài)系統(tǒng)模型(TEM)研究中得出，運(yùn)用貝葉斯推斷技術(shù)能夠約束和降低TEM預(yù)測(cè)的不確定性。Laloy 等(2009)采用貝葉斯方法結(jié)合全局敏感性分析技術(shù)模擬農(nóng)田生態(tài)系統(tǒng)中間作管理對(duì)徑流和侵蝕的影響，提取了對(duì)徑流和侵蝕有重要影響的參數(shù)，如模型中土壤參數(shù)可傳遞較大的不確定性。Eriksson等(2019)以參數(shù)的后驗(yàn)概率分布來完成敏感性分析，獲得了對(duì)模型輸出不確定性有重大貢獻(xiàn)的參數(shù)，并指出該方法既可用作試驗(yàn)設(shè)計(jì)，亦可用于模型構(gòu)建。綜上可知，貝葉斯方法結(jié)合全局敏感性分析技術(shù)在處理模型參數(shù)的不確定性傳遞方面具有特定優(yōu)勢(shì)。

本研究以秦嶺松櫟林地帶性樹種油松幼樹為研究對(duì)象，結(jié)合貝葉斯統(tǒng)計(jì)框架與全局敏感性分析技術(shù)，對(duì)油松幼樹5年的樹高生長(zhǎng)進(jìn)行預(yù)測(cè)，通過MCMC抽樣方法獲得模型參數(shù)的后驗(yàn)分布空間，量化模型輸入的不確定性及傳導(dǎo)過程，從模型誤差的不確定性、輸入變量(自變量測(cè)量誤差)的不確定性、模型參數(shù)的不確定性等入手，分析模型中不確定性的來源、參數(shù)對(duì)模型輸出不確定性的貢獻(xiàn)和影響，以期為提高幼樹樹高生長(zhǎng)建模的可靠性提供理論依據(jù)。

1 研究區(qū)概況

研究區(qū)位于秦嶺南坡中段和西段，研究樣地設(shè)在西北農(nóng)林科技大學(xué)火地塘教學(xué)林場(chǎng)、陜西省寧東林業(yè)局旬陽(yáng)壩林場(chǎng)和寶雞市辛家山林場(chǎng)?；鸬靥亮謭?chǎng)(108°21′—108°39′E，33°18′—33°28′N)地處秦嶺南坡中段安康市寧陜縣，處于北亞熱帶北緣，屬我國(guó)北亞熱帶和暖溫帶的過渡地帶，年平均氣溫8～10 ℃，海拔1 420～2 474 m，坡度20°～50°，年降雨量900～1 200 mm，年蒸發(fā)量800～950 mm，年日照時(shí)數(shù)1 668.4 h，無霜期216天，生長(zhǎng)期約6個(gè)月。土壤以山地棕壤、暗棕壤和山地草甸土為主。森林植被屬溫帶針闊混交林和寒帶針葉林?，F(xiàn)有森林是原生植被在20世紀(jì)60、70年代主伐后恢復(fù)起來的以華山松(Pinusarmandii)林、松櫟林和樺木林為主的次生林，其中松櫟林的主要樹種組成為銳齒櫟(Quercusalienavar.acuteserrata)、 白樺(Betulaplatyphylla)、油松(Pinustabulaeformis)、 華山松、青扦(Piceawilsonii)等。旬陽(yáng)壩林場(chǎng)(103°58′—109°48′E，32°29′—33°13′N)與火地塘林場(chǎng)相鄰，地貌以中山為主，地勢(shì)南高北低，平均海拔1 300 m，年均氣溫10 ℃，年均降雨量1 133 mm，年均蒸發(fā)量1 221.9 mm，年日照時(shí)數(shù) 1 638.3 h，無霜期199天。土壤為礦礫質(zhì)壤黏土，呈微酸性。森林植被屬暖溫帶落葉闊葉林和針闊混交林向北亞熱帶常綠落葉混交林的過渡帶。辛家山林場(chǎng)(106°26′—106°38′E，34°10′—34°20′N)地處秦嶺南坡西段的陜西省鳳縣，坡度 25°～50°，北高南低，海拔在1 400～2 700 m 之間，年平均氣溫11 ℃，年降雨量600～900 mm，土壤以棕壤、黃棕壤和山地草甸土為主。林區(qū)屬暖溫帶半濕潤(rùn)山地氣候區(qū)，以天然次生林群落為主，資源豐富，主要樹種有銳齒櫟、華山松、油松、漆樹(Toxicodendronvernicifluum)、山楊(Populusdavidiana)、云杉(Piceaasperata)等。

2 研究方法

2.1 樣地設(shè)置

數(shù)據(jù)采集自2013—2018年在秦嶺南坡火地塘林場(chǎng)、旬陽(yáng)壩林場(chǎng)和辛家山林場(chǎng)調(diào)查的松櫟混交林固定樣地(20 m×20 m)和臨時(shí)樣地(角規(guī)樣地)共152塊。對(duì)固定樣地內(nèi)所有達(dá)到起測(cè)胸徑(DBH≥5 cm)的喬木樹種進(jìn)行每木檢尺，調(diào)查樹高、胸徑、密度、枝下高、冠幅等信息，記錄樣地經(jīng)緯度坐標(biāo)、海拔、坡度、坡向、坡位和林分郁閉度。在固定樣地四周和中心布置2 m×2 m小樣方，在臨時(shí)樣地中心布置半徑4 m圓形小樣方，調(diào)查小樣方中所有樹高≥0.3 m、胸徑≤5 cm的幼樹。由于油松幼樹在樣地中更新較少，最終收集132株油松幼樹樣木數(shù)據(jù)進(jìn)行建模，采用計(jì)數(shù)輪生枝方法獲得油松樹高連年生長(zhǎng)量。樣地位置及油松幼樹分布情況見圖1，調(diào)查樣地基本信息見表1。

圖1 調(diào)查樣地位置及幼樹分布情況Fig.1 The locations of sampling plots and the distributions of height of the selected young trees in Qinling Mountains數(shù)字高程圖采用90 m分辨率，來自國(guó)際農(nóng)業(yè)研究磋商組織 (https:∥cgiarcsi.community)。秦嶺范圍基于世界自然基金會(huì)出版的生態(tài)區(qū)(https:∥ www.worldwildlife.org/biomes/)繪制。a、b、c中柱狀圖分別表示3個(gè)林區(qū)油松幼樹分布情況，其中縱軸表示幼樹高度級(jí)，橫軸表示不同高度的幼樹數(shù)量。The 90 m-DEM (digital elevation models) was downloaded from the CGIAR-CSI GeoPortal (http:∥srtm.csi.cgiar.org). Qinling Mountains was drawn based on ecological area published by the World Wildlife Fund (WWF). Subplot of Fig. 1 for a,b and c is frequency of young trees of P. tabulaeformis for three forest regions. The vertical axis represents height class, and the horizontal axis represent the frequency of young trees.

表1 油松幼樹調(diào)查樣地基本信息Tab.1 Statistics of young P. tabulaeformis plots

2.2 研究方法

2.2.1 幼樹樹高生長(zhǎng)模型 參考現(xiàn)有幼樹樹高生長(zhǎng)模型(Wykoffetal., 1982； Boisvenueetal., 2004)，應(yīng)用優(yōu)選出的擬合效果最佳的對(duì)數(shù)線性模型(表2中模型4)構(gòu)建秦嶺地區(qū)油松幼樹樹高生長(zhǎng)模型，以量化和分析模型預(yù)測(cè)中的不確定性來源。模型總體結(jié)構(gòu)如下：

ln(HTG)=f(SIZE,SITE,COMP)=

SIZE+SITE+COMP+ε。

(1)

該模型結(jié)構(gòu)表達(dá)為： 油松幼樹樹高5年生長(zhǎng)量(HTG)是林木大小(SIZE)、林分競(jìng)爭(zhēng)因子(COMP)和立地條件(SITE)的函數(shù)。在建模前采用逐步回歸法從初始變量林木大小(樹高)、競(jìng)爭(zhēng)因子(樹冠競(jìng)爭(zhēng)因子、大于對(duì)象木的林分?jǐn)嗝娣e和、光截留)和立地因子(坡度、坡向、坡位、海拔、林分優(yōu)勢(shì)高)中選擇最優(yōu)變量，并在建模時(shí)直接采用最終變量。

各變量公式及解釋如下：

SITE表示立地因子，坡度、坡向及其相互作用主要通過改變光照和水分條件影響幼樹樹高生長(zhǎng)，坡度改變會(huì)使土壤厚度發(fā)生變化，從而改變土壤水分和養(yǎng)分條件影響幼樹生長(zhǎng)(Boisvenueetal., 2004)。其公式為：

SITE=β0+β1×SL×cos(ASP)-β2×

SL×SIN(ASP)-β3×SL。

(2)

SIZE表示林木大小，采用當(dāng)前樹高表達(dá)，代表樹種生理特征和垂直結(jié)構(gòu)。有研究表明，當(dāng)前樹高代表樹種生長(zhǎng)潛力，作為一種重要的因子與其生長(zhǎng)量呈正比(Boisvenueetal.， 2004)。其公式為：

SIZE=β4×ln(HT)。

(3)

COMP表示林分競(jìng)爭(zhēng)因子，以樹冠競(jìng)爭(zhēng)因子和光截留來表示，主要反映樹種在林內(nèi)競(jìng)爭(zhēng)和光照條件下的樹高生長(zhǎng)變化(張仰渠等， 1989)。其公式為：

COMP=β5×ln(CCF)+β6×LI。

(4)

以上公式中:HTG為油松幼樹樹高5年增長(zhǎng)量(m)；SL為坡度的正切值；cos(ASP)為坡向的余弦值；sin(ASP)為坡向的正弦值；ln為自然對(duì)數(shù)；β0～β6為模型估計(jì)參數(shù)； ε為模型誤差項(xiàng)。

LI為光截留，代表林分內(nèi)的光照條件。其公式為：

LI=1-e-K×LAI。

(5)

式中：K為消光系數(shù)，取值0.14(劉勝， 2006； 宋子煒等， 2009)。

LAI為累計(jì)的葉面積指數(shù)，通過計(jì)算生物量方程獲得(陳存根等， 1996)：

LAI =1.345 + 0.501W,

W=2.458 44ln(DBH)- 5.869 02。

式中：W為油松幼樹生物量(1 t·hm-2)； DBH 為胸徑(cm)。

表2 油松幼樹樹高生長(zhǎng)模型①Tab.2 Height models for young P. tabulaeformis

2) 幼樹樹高生長(zhǎng)貝葉斯模型 本研究構(gòu)建的油松幼樹樹高生長(zhǎng)貝葉斯模型如下：

ln(HTGi)=β0+β1×x1i-β2×x2i-β3×x3i+

β4×x4i+β5×x5i+β6×x6i+εi。

(6)

式中： HGTi為油松幼樹樹高5年增長(zhǎng)量；β0～β6為模型估計(jì)參數(shù)；x1i=SLi×cos(ASPi)，x2i=SLi×sin(ASPi)，x3i=SLi，x4i=ln(HTi)，x5i為樹冠競(jìng)爭(zhēng)因子在考慮測(cè)量誤差情況下的值,x5i=ln(CCFi)+ui，ui～N(0,τ2)，x6i為光截留在考慮測(cè)量誤差情況下的值，x6i=LIi+νi；εi為模型誤差，εi～N(0,σ2)，σ2為模型殘差的方差。

測(cè)量誤差一般假定均值為0、方差為σ2的正態(tài)分布，也可以取一定區(qū)間上的均勻分布，本研究考慮正態(tài)分布更接近于實(shí)際問題中的誤差(李永慈等， 2005)，采用服從N(0,σ2)的隨機(jī)變量作為模型的測(cè)量誤差效應(yīng)。對(duì)于考慮測(cè)量誤差后的自變量取值，采用服從正態(tài)分布的模糊信息先驗(yàn)分布(Clark， 2007)。

2.2.3 模型參數(shù)敏感性分析及不確定性度量 1)Sobol全局敏感性分析技術(shù) 俄羅斯數(shù)學(xué)家I.M. Sobol于20世紀(jì)90年代提出了一種全局敏感性分析技術(shù)，稱之為“Sobol 指數(shù)法”，其核心思想是方差分解，將模型分解為單個(gè)參數(shù)以及參數(shù)之間相互組合的函數(shù)，通過計(jì)算單個(gè)輸入?yún)?shù)或參數(shù)集的方差對(duì)總輸出方差的影響分析參數(shù)的重要性及參數(shù)之間的交互效應(yīng)。計(jì)算采用Saltelli(2002)提出的方法，通過固定單一參數(shù)，使模型輸出方差能減少多少來度量模型中該參數(shù)對(duì)模型輸出結(jié)果的敏感性。

2) 基于貝葉斯方法的敏感性分析及不確定性度量 第一步，采用MCMC抽樣方法，估計(jì)每個(gè)參數(shù)的后驗(yàn)分布空間。MCMC方法是在貝葉斯理論框架下，建立一個(gè)平穩(wěn)分布為p(θ|y)后驗(yàn)分布的隨機(jī)樣本，采用MCMC中DREAMZS算法來實(shí)現(xiàn)(Vrugtetal.， 2009)。第二步，從貝葉斯參數(shù)后驗(yàn)分布空間中抽樣，將抽取的參數(shù)樣本分為2部分，所得參數(shù)矩陣分別為A和B，如下：

(7)

矩陣B的第j列元素用矩陣A的第j列元素替代，獲得新的矩陣Cj：

(8)

根據(jù)以上抽樣方法，假定模型總輸出方差為V(Y)，模型輸出為Y，在固定參數(shù)Xi的情況下，模型輸出的總方差為：V(Y*)i=V[E(Y*|Xi)]+E[V(Y*|Xi)]，其中V[E(Y*|Xi)]為模型輸出的條件方差，E[V(Y*|Xi)]為模型輸出方差的均值。

一階敏感性效應(yīng)為：

(9)

二階敏感性效應(yīng)為：

(10)

不確定性度量采用參數(shù)不確定性傳遞的模型輸出變異系數(shù)(CV，%)和貝葉斯95%可信區(qū)間寬度(B)來表示(Xiongetal.， 2009； Maetal.， 2018)，公式如下：

(11)

(12)

2.2.4 貝葉斯參數(shù)估計(jì)及收斂性診斷 考慮到參數(shù)聯(lián)合后驗(yàn)分布不可解析，故本研究采用馬爾可夫鏈蒙特卡洛抽樣方法估計(jì)模型參數(shù)的聯(lián)合后驗(yàn)分布。運(yùn)用DREAMZS算法進(jìn)行參數(shù)估計(jì)，DREAMZS算法為DREAM自適應(yīng)抽樣方法，其結(jié)合了差分演化算法和自適應(yīng)Metropolis算法的優(yōu)點(diǎn)，可同時(shí)運(yùn)行多條馬爾科夫鏈進(jìn)行全局搜索，并可自動(dòng)調(diào)整建議分布的范圍和方向，在復(fù)雜高維、非線性和多峰目標(biāo)問題中表現(xiàn)出極高效率(Vrugtetal.， 2009)。DREAMZS算法是在已有樣本中產(chǎn)生候選樣本，從而大大減少算法中并行馬爾科夫鏈的數(shù)量，加速收斂到后驗(yàn)分布的速率。參數(shù)收斂性診斷選用Gelman等(1992)提出的指標(biāo)R(尺度縮減因子)，如果R<1.2，則認(rèn)為待估參數(shù)收斂于平穩(wěn)的后驗(yàn)分布；但本研究中，執(zhí)行嚴(yán)格的收斂診斷標(biāo)準(zhǔn)： 接受閾值R<1.05(Brooksetal.， 1998)，最終執(zhí)行3×107次迭代生成，在收斂性診斷和預(yù)測(cè)時(shí)舍棄掉10%的“燃燒值”(Gilksetal.， 1996)。

2.2.5 模型檢驗(yàn) 1) 交叉檢驗(yàn)法 將調(diào)查數(shù)據(jù)隨機(jī)等分為兩部分，第一部分作為訓(xùn)練集進(jìn)行建模，第二部分為驗(yàn)證數(shù)據(jù)，同時(shí)再以第二部分?jǐn)?shù)據(jù)作為訓(xùn)練集進(jìn)行建模，第一部分為驗(yàn)證數(shù)據(jù)檢驗(yàn)?zāi)Ｐ皖A(yù)測(cè)精度，通過均方根誤差(RMSE)表示。計(jì)算公式如下：

(13)

圖2 油松幼樹5年樹高生長(zhǎng)預(yù)測(cè)的95%不確定性區(qū)間Fig.2 95% credible intervals of 5-year height increment model due to parameter uncertainty and total uncertaintyx軸表示每株油松幼樹樣本，黑色圓點(diǎn)表示每株油松幼樹的觀測(cè)值，灰色陰影部分表示預(yù)測(cè)的貝葉斯95%可信區(qū)間。x axis represents each young P. tabulaeformis tree sample. The black dots represent observations of each young P. tabulaeformis, and the grey shade region represents the 95% Bayesian credible interval.

2) 等效性檢驗(yàn)法 基于回歸的等效性檢驗(yàn)是由Robinson等(2005)提出的。不同于傳統(tǒng)的假設(shè)檢驗(yàn)(差異性檢驗(yàn)表明2個(gè)群體之間存在顯著性差異)，等效性檢驗(yàn)假設(shè)模型與數(shù)據(jù)是不同的，檢驗(yàn)的顯著性代表在某個(gè)研究區(qū)范圍內(nèi)模型與數(shù)據(jù)是相似的，因此拒絕原假設(shè)是等效性檢驗(yàn)的理想目標(biāo)?；诨貧w的等效性檢驗(yàn)考慮均值差異和個(gè)體差異，通過以下6個(gè)步驟完成： (1) 預(yù)測(cè)值減去預(yù)測(cè)均值； (2) 確定一個(gè)截距與斜率可搖擺的等效性限制區(qū)域 (如±25%)； (3) 檢驗(yàn)實(shí)測(cè)值與調(diào)整后預(yù)測(cè)值的線性回歸； (4) 通過計(jì)算2個(gè)單側(cè)置信區(qū)間得到線性回歸截距的等效性，同時(shí)與設(shè)定的等效性區(qū)間進(jìn)行比較； (5) 通過計(jì)算2個(gè)單側(cè)置信區(qū)間得到回歸斜率的等效性，同時(shí)與設(shè)定的等效性區(qū)間進(jìn)行比較； (6)根據(jù)置信區(qū)間是否落在等效區(qū)間內(nèi)決定是否拒絕實(shí)測(cè)與預(yù)測(cè)的非相似性原假設(shè)。

3 結(jié)果與分析

3.1 模型整體表現(xiàn)

對(duì)油松幼樹樹高5年生長(zhǎng)量進(jìn)行模擬，考慮模型預(yù)測(cè)時(shí)模型預(yù)測(cè)誤差的不確定性、輸入變量(自變量測(cè)量誤差)的不確定性和模型參數(shù)的不確定性，用貝葉斯方法對(duì)3種不確定性進(jìn)行定量描述(圖2)。圖2中陰影區(qū)域灰度由深到淺分別代表由參數(shù)傳遞的模型輸出不確定性(parameters uncertainty)、由參數(shù)和輸入變量(自變量測(cè)量誤差)共同傳遞的模型輸出不確定性(parameters and input uncertainty)以及模型預(yù)測(cè)總體不確定性(total uncertainty)，不確定性的量化以貝葉斯95%可信區(qū)間寬度表示。3種不確定性中，模型預(yù)測(cè)誤差的不確性對(duì)模型輸出不確定性的影響最大，其次是模型參數(shù)的不確定性，輸入變量(自變量測(cè)量誤差)傳遞的不確定性最小。通過計(jì)算模型預(yù)測(cè)的CV對(duì)各部分不確定性進(jìn)行量化，得到由模型參數(shù)傳遞的不確定性占總體不確定性的43%，自變量(光截留和樹冠競(jìng)爭(zhēng)因子)測(cè)量誤差傳遞的不確定性占總體不確性的6%，模型預(yù)測(cè)誤差傳遞的不確定性占51%。從模型的不確定性區(qū)間和數(shù)據(jù)的分布程度看出，由模型參數(shù)傳遞的不確定性區(qū)間覆蓋了60%以上的實(shí)測(cè)數(shù)據(jù)，模型總體不確定性區(qū)間基本覆蓋了97%以上的實(shí)測(cè)數(shù)據(jù)(圖2)，可見，所構(gòu)建模型較準(zhǔn)確地表達(dá)了預(yù)測(cè)中的不確定性。為評(píng)價(jià)模型預(yù)測(cè)效果，采用交叉檢驗(yàn)法(圖3)和等效性檢驗(yàn)法(表3)檢驗(yàn)?zāi)Ｐ皖A(yù)測(cè)精度。圖3中貝葉斯模型的預(yù)測(cè)值和實(shí)測(cè)值擬合曲線接近于理想狀態(tài)(斜率為45°的直線，反映實(shí)測(cè)值與預(yù)測(cè)值相等的擬合效果)，經(jīng)計(jì)算得到貝葉斯模型預(yù)測(cè)的RMSE為0.14，意味著模型具有較好的不確定性區(qū)間覆蓋度和擬合度。同時(shí)，等效性檢驗(yàn)結(jié)果表明，在對(duì)截距的檢驗(yàn)中其等效性區(qū)間[0.44, 0.74]包含模型預(yù)測(cè)的置信區(qū)間[0.58, 0.72]，斜率的檢驗(yàn)也表現(xiàn)出相同結(jié)果，根據(jù)等效性檢驗(yàn)原理，截距和斜率檢驗(yàn)均拒絕原假設(shè)，其模型預(yù)測(cè)數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)之間無顯著性差異。

圖3 油松幼樹5年樹高生長(zhǎng)模型的交叉檢驗(yàn)Fig.3 Cross-validation of 5-year height increment model of young P. tabulaeformis交叉檢驗(yàn)中調(diào)查數(shù)據(jù)隨機(jī)等分為兩部分，檢驗(yàn)1以第一部分為建模數(shù)據(jù)，第二部分為驗(yàn)證數(shù)據(jù)，而檢驗(yàn)2以第二部分為建模數(shù)據(jù)，第一部分為驗(yàn)證數(shù)據(jù)。The cross validation randomly split data into two parts. Test 1 used first half for fitting and second half for validation, while test 2 used second half for fitting and first half for validation.

表3 油松幼樹樹高生長(zhǎng)模型等效性檢驗(yàn)結(jié)果①Tab.3 Summary of equivalence-based regression results for height model of young P. tabulaeformis

3.2 參數(shù)不確定性傳遞分析

本研究采用MCMC抽樣方法估計(jì)每個(gè)參數(shù)的貝葉斯95%可信區(qū)間(貝葉斯可信區(qū)間與經(jīng)典統(tǒng)計(jì)中的置信區(qū)間有所不同，貝葉斯可信區(qū)間是指參數(shù)屬于該區(qū)間的概率，而置信區(qū)間是指隨機(jī)區(qū)間能把參數(shù)包在該區(qū)間的概率)和最大后驗(yàn)概率(MAP)(表4)，從參數(shù)后驗(yàn)分布空間抽樣量化模型參數(shù)的不確定性和敏感性。結(jié)果表明，在油松幼樹5年生長(zhǎng)量預(yù)測(cè)中，β3(坡度)、β4(樹高)、β5(樹冠競(jìng)爭(zhēng)因子)和β6(光截留)對(duì)模型輸出結(jié)果較敏感，敏感性指數(shù)均在1%以上。4個(gè)參數(shù)Sobol敏感性指數(shù)的中位數(shù)和[25%，75%]的分位數(shù)依次為0.05、[0.04，0.06]，0.06、[0.03，0.09]，0.40、[0.39，0.45]和0.04、[0.02，0.05](圖4a)，其平均敏感性指數(shù)分別為5%、1%、41%和4%。4個(gè)參數(shù)中，油松幼樹樹高預(yù)測(cè)結(jié)果對(duì)β5最敏感，其次為β6和β3，β4敏感性最低。通過模型輸出變異系數(shù)和貝葉斯95%可信區(qū)間寬度來度量每個(gè)參數(shù)或參數(shù)組合傳遞給模型輸出不確定性的貢獻(xiàn)(圖4b、c)。結(jié)果表明，對(duì)油松幼樹樹高預(yù)測(cè)不確定性貢獻(xiàn)最大的是控制樹冠競(jìng)爭(zhēng)因子的參數(shù)，其模型輸出變異系數(shù)和95%可信區(qū)間寬度均值分別為49.03%和0.69，中位數(shù)、[25%，75%]的分位數(shù)分別為49.56%、 [48.03%，51.20%]和0.68、 [0.62，0.77]；其次為控制坡度的參數(shù)，其模型輸出變異系數(shù)和95%可信區(qū)間寬度均值分別為13.03%和0.32；控制光截留的參數(shù)傳遞的模型輸出不確定性指標(biāo)平均值分別為8.76%和0.20，其輸出較低的不確定性；控制樹高的參數(shù)對(duì)模型輸出不確定性的貢獻(xiàn)也較小，變異系數(shù)和95%可信區(qū)間寬度均值分別為1.24%和0.04。在參數(shù)相互作用中，β3和β5相互作用對(duì)模型輸出不確定性的貢獻(xiàn)CV為6.61%，其他參數(shù)兩兩之間對(duì)模型輸出不確定性的貢獻(xiàn)均很小，低于1%。在模型參數(shù)不確定性中，依據(jù)95%不確定性區(qū)間計(jì)算相對(duì)不確定性，得到對(duì)模型輸出不確定性貢獻(xiàn)最大的是控制樹冠競(jìng)爭(zhēng)因子的參數(shù)，占參數(shù)總體不確定性的64.87%； 其次是控制坡度和光截留的參數(shù)，分別占參數(shù)總體不確定性的15.88% 和10.02%； 控制樹高的參數(shù)對(duì)模型輸出不確定性的貢獻(xiàn)僅為1.78%，其他參數(shù)貢獻(xiàn)的不確定性均低于1%。

表4 油松幼樹樹高生長(zhǎng)模型的參數(shù)后驗(yàn)分布Tab.4 Posterior probability distribution of parameters in the height model of young P. tabulaeformis

圖4 參數(shù)傳遞不確定性量化(變異系數(shù)和95%可信區(qū)間寬度)及參數(shù)敏感性Fig.4 Uncertainty quantification (CV and 95% credible interval) and sensitivity of 5-year height increment model due to parameter uncertainty

3.3 模型預(yù)測(cè)的不確定性

本研究通過改變其中1個(gè)變量并控制其他變量，在均值水平上采用貝葉斯參數(shù)后驗(yàn)分布空間和參數(shù)MAP即最大后驗(yàn)概率預(yù)測(cè)油松幼樹樹高5年生長(zhǎng)量(圖5a-d)。結(jié)果表明， 樹高與5年樹高生長(zhǎng)量呈正相關(guān)，模型預(yù)測(cè)的不確定性區(qū)間基本覆蓋模型預(yù)測(cè)的誤差范圍(均值±標(biāo)準(zhǔn)差)；同時(shí)，從預(yù)測(cè)值的變化幅度看出，樹高對(duì)模型預(yù)測(cè)結(jié)果影響較顯著。當(dāng)樹高為2 m時(shí)，樹高生長(zhǎng)量預(yù)測(cè)值和不確定性區(qū)間分別為0.42 m、[0.21 m, 0.86 m]，當(dāng)樹高達(dá)到3 m時(shí)，樹高生長(zhǎng)量預(yù)測(cè)值和不確定性區(qū)間分別為0.78 m、[0.39 m，1.6 m]，隨著樹高增大，其樹高5年生長(zhǎng)量預(yù)測(cè)值和不確定性區(qū)間均增加。光照條件(光截留)與樹高5年生長(zhǎng)量呈負(fù)相關(guān)，且影響較顯著。當(dāng)光截留從0.4增至0.8時(shí)，樹高生長(zhǎng)量預(yù)測(cè)值和不確定性區(qū)間相應(yīng)從0.71 m、 [0.35 m，1.47 m]降至0.54 m、 [0.27 m，1.10 m]。坡度與樹高5年生長(zhǎng)量呈負(fù)相關(guān)，隨著坡度從緩坡(10°)到陡坡(30°)，其樹高5年生長(zhǎng)量預(yù)測(cè)值和不確定性區(qū)間分別從0.79 m、 [0.38 m,1.64 m]降至0.66 m、[0.33 m，1.34 m]。樹冠競(jìng)爭(zhēng)因子與樹高5年生長(zhǎng)量也呈負(fù)相關(guān)關(guān)系，樹冠競(jìng)爭(zhēng)因子從50增至350，其樹高5年生長(zhǎng)量預(yù)測(cè)值和不確定性區(qū)間從0.70 m、 [0.34 m,1.45 m]降至0.57 m、 [0.28 m,1.17 m]。從變量模擬效果可看出，當(dāng)前樹高對(duì)模型預(yù)測(cè)結(jié)果變化影響更明顯(圖5a)，其次為光截留(圖5b)，而樹冠競(jìng)爭(zhēng)因子對(duì)模型預(yù)測(cè)結(jié)果變化影響最小(圖5d)。結(jié)合參數(shù)的不確定性分析得出，對(duì)模型輸出不確定性貢獻(xiàn)大的參數(shù)，其控制的變量對(duì)模型輸出結(jié)果的影響較?。环粗?，對(duì)模型輸出不確定性貢獻(xiàn)小的參數(shù)，其控制的變量對(duì)模型輸出結(jié)果的影響較大。由此可見，參數(shù)不確定性引起的模型輸出不確定性越低，其相應(yīng)的變量對(duì)模型輸出結(jié)果的影響越顯著。

圖5 不同的HT、CCF、LI、SL條件下油松幼樹樹高5年生長(zhǎng)預(yù)測(cè)Fig.5 Height increment on 5-year height increment of young P. tabulaeformis based on HT, CCF, LI, and SL圖中灰色陰影區(qū)域?yàn)樨惾~斯預(yù)測(cè)95%可信區(qū)間，黑色垂直線為預(yù)測(cè)均值±1倍標(biāo)準(zhǔn)差，黑色圓點(diǎn)為MAP預(yù)測(cè)值。圖5a中，ASP=200°，CCF=150，LI=0.6，SL=35°； 圖5b中，ASP=200°，CCF=150，HT=2 m，SL=35°； 圖5c中，ASP=200°，CCF=150，HT=2 m，LI=0.6； 圖5d中，ASP=200°，SL=35°，HT=2 m，LI=0.6。The grey shaded region represented 95% credible interval of Bayesian prediction, black lines represented mean ± one standard deviation, and black dot represented prediction by Bayesian MAP estimation. For Fig.5a, the aspect (ASP) is 200°, the crown competition factor(CCF)is 150, the light interception(LI)is 0.6, slope (SL) is 35°. For Fig.5b, the aspect (ASP) is 200°, the crown competition factor(CCF)is 150, the current height(HT)is 2 m, and the slope (SL) is 35°. For Fig.5c, the aspect (ASP) is 200°, the crown competition factor(CCF)is 150, the current height(HT)is 2 m, and the light interception(LI)is 0.6. For Fig.5d, the aspect (ASP) is 200°, the current height(HT)is 2 m, the light interception(LI)is 0.6, and the slope (SL) is 35°.

4 討論

4.1 模型不確定性量化與分析

本研究結(jié)合貝葉斯統(tǒng)計(jì)框架與Sobol全局敏感性分析技術(shù)，從模型參數(shù)后驗(yàn)分布空間中抽樣，量化每個(gè)參數(shù)或參數(shù)組合傳遞給模型輸出的不確定性，該方法比一般參數(shù)搜索方法更具優(yōu)勢(shì)，因?yàn)槠淇紤]到了模型參數(shù)空間的變化(Erikssonetal.， 2019)。結(jié)果表明，模型預(yù)測(cè)的不確定性區(qū)間基本覆蓋97%以上的觀測(cè)數(shù)據(jù)，所構(gòu)建模型很好表達(dá)了實(shí)測(cè)數(shù)據(jù)分布的隨機(jī)性。在所有不確定性成分中，模型參數(shù)的不確定性是眾多不確定性來源之一(McMahonetal.， 2009)，如果約束一個(gè)模型的參數(shù)為一固定值而非概率分布，就忽視了參數(shù)的不確定性，不能確保模型與實(shí)際數(shù)據(jù)相匹配，也很難評(píng)估模型預(yù)測(cè)的準(zhǔn)確性(Vanlieretal.， 2013)。采用貝葉斯統(tǒng)計(jì)框架結(jié)合全局敏感性分析技術(shù)比單獨(dú)全局敏感性參數(shù)抽樣方法具有更好的分析能力，因?yàn)槠洳捎门c模型數(shù)據(jù)緊密聯(lián)系的參數(shù)空間來分析模型的不確定性，能更多解釋變量的作用(Erikssonetal.， 2019)。了解不同參數(shù)對(duì)模型輸出不確定性的貢獻(xiàn)可以幫助建模工作者快速制定數(shù)據(jù)收集和模型改進(jìn)方案來降低模型預(yù)測(cè)的不確定性。此類研究在森林生物量模型中已有報(bào)道，如Chen等(2015)、Chave等(2004)對(duì)熱帶雨林地區(qū)生物量模型誤差的估算表明，模型誤差具有很大的不確定性；而秦立厚等(2017)在森林生物量估算模型不確定性分析中得出，模型參數(shù)的不確定性大于模型誤差的不確定性，原因可能是建模數(shù)據(jù)量偏小導(dǎo)致模型參數(shù)具有較高不確定性(傅煜等， 2014)。本研究表明，模型的不確定性主要來源于模型預(yù)測(cè)誤差的不確定性和模型參數(shù)的不確定性。在模型參數(shù)不確定性研究中，控制樹冠競(jìng)爭(zhēng)因子的參數(shù)對(duì)油松幼樹樹高生長(zhǎng)預(yù)測(cè)的不確定性貢獻(xiàn)很高，其次為坡度和光截留，而控制樹高的參數(shù)傳遞的不確定性較低。除了β3(坡度)和β5(樹冠競(jìng)爭(zhēng)因子)外，其他參數(shù)的相互作用對(duì)模型輸出不確定性的貢獻(xiàn)較小，這也反映了模型預(yù)測(cè)中參數(shù)的相互作用對(duì)模型輸出的不確定性具有一定意義。因此，在模型預(yù)測(cè)分析時(shí)參數(shù)相互作用的分析依然很重要，有時(shí)只考慮單獨(dú)參數(shù)對(duì)模型輸出不確定性的影響很難確定該參數(shù)對(duì)模型是重要還是無效 (LeBaueretal.， 2013)。

4.2 模型應(yīng)用效果與預(yù)測(cè)分析

目前，多數(shù)樹高生長(zhǎng)模型的研究對(duì)象以大樹(DBH≥5 cm)為主，專門針對(duì)幼樹樹高生長(zhǎng)模型的研究較少，應(yīng)用較成熟的幼樹樹高生長(zhǎng)模型代表為Stage(1973)和Wykoff等(1982)開發(fā)的Prognosis模型系統(tǒng)(Wykoff， 1986； Boisvenue， 2000； Boisvenueetal., 2004； Lacerteaetal.， 2006)。本研究在模型應(yīng)用前對(duì)3種類型(線性、對(duì)數(shù)線性和非線性)幼樹樹高生長(zhǎng)模型中的6個(gè)模型進(jìn)行比較，根據(jù)BIC和RMSE選出最優(yōu)的對(duì)數(shù)線性模型(表2中模型4)應(yīng)用于油松幼樹樹高生長(zhǎng)預(yù)測(cè)和不確定性分析；為了解模型適用性，對(duì)參數(shù)化后的模型進(jìn)行交叉檢驗(yàn)和等效性檢驗(yàn)。結(jié)果表明，模型預(yù)測(cè)數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)之間無顯著性差異，模型具有較高的預(yù)測(cè)精度，同時(shí)也表現(xiàn)出與實(shí)測(cè)數(shù)據(jù)隨機(jī)分布相匹配的優(yōu)良不確定性區(qū)間。因此，該模型適用于本研究區(qū)油松幼樹樹高的生長(zhǎng)預(yù)測(cè)和不確定性研究。

在建模時(shí)變量的使用主要考慮以下2方面： 1) 盡量保留原始模型中的變量，因?yàn)槟Ｐ徒?jīng)過學(xué)者們的不斷改進(jìn)和驗(yàn)證，其中的變量具有更可靠的解釋效應(yīng)； 2) 在此基礎(chǔ)上結(jié)合目標(biāo)樹種的生物學(xué)特性以及在研究區(qū)的生長(zhǎng)規(guī)律進(jìn)行適當(dāng)修正，即以具有廣泛生物生態(tài)學(xué)意義的關(guān)鍵性變量為主。目前，樹高生長(zhǎng)模型中常用的競(jìng)爭(zhēng)因子為樹冠競(jìng)爭(zhēng)因子、大于對(duì)象木的林分?jǐn)嗝娣e和， 立地因子為坡度、坡向、坡位、海拔、林分優(yōu)勢(shì)高(Wykoffetal., 1982； Boisvenueetal., 2004； 張西等， 2015； 衣旭彤等， 2017)。本研究考慮到油松樹種喜光的生物學(xué)特性，在競(jìng)爭(zhēng)因子中增加光照因子(光截留)，并通過逐步回歸法剔除了對(duì)模型沒有貢獻(xiàn)的變量，最終在競(jìng)爭(zhēng)因子中保留了樹冠競(jìng)爭(zhēng)因子和光截留，在立地因子中保留了坡度和坡向(Wykoffetal., 1982； Wykoff, 1986； Boisvenue， 2000; Boisvenueetal., 2004),與張仰渠等(1989)在秦嶺地區(qū)油松幼樹生長(zhǎng)的主要影響因子研究結(jié)果一致。因此，結(jié)合國(guó)內(nèi)外研究成果、油松樹種生物學(xué)特性和本研究變量逐步回歸綜合分析， 坡度和坡向是能夠很好表達(dá)油松幼樹樹高生長(zhǎng)的立地因子，而樹冠競(jìng)爭(zhēng)因子和光截留可以作為很好表達(dá)油松幼樹樹高生長(zhǎng)的競(jìng)爭(zhēng)因子。

對(duì)參數(shù)化后的最大后驗(yàn)概率進(jìn)行預(yù)測(cè)分析，不同參數(shù)在解釋模型預(yù)測(cè)不確定性方面差異較大，保持其他變量在均值水平，通過改變某一變量進(jìn)行模型預(yù)測(cè)(Weiskitteletal.， 2011)，可以反映出該變量對(duì)模型預(yù)測(cè)的效果。本研究結(jié)果表明，油松幼樹樹高5年生長(zhǎng)量與當(dāng)前樹高正相關(guān)，主要是因?yàn)闃涓唧w現(xiàn)了樹種的光合能力和垂直結(jié)構(gòu)，對(duì)其生長(zhǎng)具有積極作用(Boisvenueetal., 2004)。油松幼樹樹高生長(zhǎng)與樹冠競(jìng)爭(zhēng)因子和光截留呈明顯負(fù)相關(guān)，隨著林內(nèi)競(jìng)爭(zhēng)和光照條件減弱，油松幼樹樹高生長(zhǎng)量顯著降低。Yu等(2013)在我國(guó)秦嶺松櫟混交林環(huán)境因子對(duì)幼苗更新的影響中得出油松種子萌發(fā)和幼苗初期需要更多的光照條件，應(yīng)證了本研究結(jié)果。張仰渠等 (1989)指出油松樹高生長(zhǎng)隨上層林冠郁閉度增加，光照減弱，光合強(qiáng)度大大降低，特別是隨著年齡增加對(duì)光照的需求越來越強(qiáng)，亦與本研究結(jié)果一致。在坡度對(duì)油松生長(zhǎng)的研究中，張仰渠等(1989)指出水熱條件較好的平緩坡有利于油松樹種生長(zhǎng)，而本研究也得出隨著坡度增加油松幼樹樹高5年生長(zhǎng)量減小，主要是由于坡度大小與土壤厚薄相關(guān)，坡度太大的地方土層較薄，土壤養(yǎng)分會(huì)發(fā)生變化，阻礙樹種生長(zhǎng)(Boisvenueetal., 2004)。從模型預(yù)測(cè)結(jié)果看出，樹高對(duì)模型預(yù)測(cè)結(jié)果影響最顯著，其次為光截留，而樹冠競(jìng)爭(zhēng)因子對(duì)模型預(yù)測(cè)結(jié)果影響最小。結(jié)合參數(shù)的不確定性不難發(fā)現(xiàn)，參數(shù)的不確定性越高，越難確定其控制的變量對(duì)模型預(yù)測(cè)結(jié)果的影響；相反，參數(shù)的不確定性越低，越容易確定其控制的變量對(duì)模型預(yù)測(cè)結(jié)果的影響。生態(tài)建模一般更趨于選擇對(duì)模型有統(tǒng)計(jì)學(xué)意義的參數(shù)，即傳遞給模型輸出不確定性小的參數(shù)； 但有學(xué)者提出，在模型中較弱的部分更是將來的研究重點(diǎn)(Erikssonetal.， 2019)。因此，在建模中對(duì)于具有一定生態(tài)學(xué)或生物學(xué)意義且傳遞給模型輸出不確定性較大的參數(shù)是未來需要深入研究的重要目標(biāo)。

本研究構(gòu)建的貝葉斯模型不確定性分析方法主要通過貝葉斯參數(shù)后驗(yàn)概率來表述模型傳遞的不確定性，量化模型中各種來源的不確定性； 而貝葉斯統(tǒng)計(jì)框架與全局敏感性分析技術(shù)結(jié)合有效約束了模型的參數(shù)空間范圍，明確了模型預(yù)測(cè)中參數(shù)獨(dú)立或相互作用對(duì)模型輸出不確定性的貢獻(xiàn)和作用，可為建模工作者進(jìn)行目標(biāo)參數(shù)校正和模型優(yōu)化提供技術(shù)支持。但是如何改進(jìn)和優(yōu)化對(duì)模型輸出所傳遞不確定性較大的參數(shù)有待進(jìn)一步研究，這將需要后期數(shù)據(jù)的不斷收集和補(bǔ)充、調(diào)查方法或模型形式的改進(jìn)等策略，對(duì)模型參數(shù)進(jìn)行不斷校正和優(yōu)化，獲得更理想的參數(shù)后驗(yàn)分布來降低模型預(yù)測(cè)的不確定性。

5 結(jié)論

1) 采用馬爾可夫鏈蒙特卡洛抽樣方法的貝葉斯技術(shù)估計(jì)模型參數(shù)的聯(lián)合后驗(yàn)分布和模型誤差，量化模型預(yù)測(cè)中各種來源傳遞的不確定性。結(jié)果表明，在模型不確定性分析中，模型預(yù)測(cè)誤差傳遞的不確定性最大，占總體不確定性的51%，其次為模型參數(shù)的不確定性，占總體不確定性的43%，輸入變量(自變量測(cè)量誤差)傳遞的不確定性最小，僅占總體不確定性的6%。

2) 貝葉斯框架與全局敏感性分析技術(shù)結(jié)合，進(jìn)一步量化每個(gè)來自后驗(yàn)分布空間的參數(shù)或參數(shù)組合傳遞給模型輸出的不確定性。結(jié)果表明，模型中的主要參數(shù)對(duì)變量具有重要解釋意義，不同參數(shù)對(duì)模型輸出不確定性的貢獻(xiàn)差異較大，其中控制樹冠競(jìng)爭(zhēng)因子的參數(shù)對(duì)模型輸出不確定性的貢獻(xiàn)最大，其次為控制坡度和光截留的參數(shù)，而控制樹高以及其他變量相關(guān)的參數(shù)對(duì)模型輸出不確定性的貢獻(xiàn)較小。

3) 貝葉斯MAP估計(jì)和不確定性預(yù)測(cè)結(jié)果表明，隨著樹冠競(jìng)爭(zhēng)因子、光截留和坡度增大，油松幼樹樹高5年生長(zhǎng)量相應(yīng)減小，表現(xiàn)出負(fù)相關(guān)關(guān)系； 隨著樹高增大其5年生長(zhǎng)量相應(yīng)增加，表現(xiàn)出正相關(guān)關(guān)系。結(jié)合參數(shù)的不確定性進(jìn)一步分析得出，參數(shù)的不確定性越高，其控制的變量對(duì)模型預(yù)測(cè)結(jié)果的影響越不顯著。