油松幼樹樹高生長預測的不確定性貝葉斯分析*

王 彬 田相林 曹田健

(西北農林科技大學林學院 生態仿真優化實驗室 楊凌 712100)

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

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

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

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

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

1 研究區概況

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

2 研究方法

2.1 樣地設置

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

圖1 調查樣地位置及幼樹分布情況Fig.1 The locations of sampling plots and the distributions of height of the selected young trees in Qinling Mountains數字高程圖采用90 m分辨率，來自國際農業研究磋商組織 (https:∥cgiarcsi.community)。秦嶺范圍基于世界自然基金會出版的生態區(https:∥ www.worldwildlife.org/biomes/)繪制。a、b、c中柱狀圖分別表示3個林區油松幼樹分布情況，其中縱軸表示幼樹高度級，橫軸表示不同高度的幼樹數量。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 油松幼樹調查樣地基本信息Tab.1 Statistics of young P. tabulaeformis plots

2.2 研究方法

2.2.1 幼樹樹高生長模型 參考現有幼樹樹高生長模型(Wykoffetal., 1982； Boisvenueetal., 2004)，應用優選出的擬合效果最佳的對數線性模型(表2中模型4)構建秦嶺地區油松幼樹樹高生長模型，以量化和分析模型預測中的不確定性來源。模型總體結構如下：

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

SIZE+SITE+COMP+ε。

(1)

該模型結構表達為： 油松幼樹樹高5年生長量(HTG)是林木大小(SIZE)、林分競爭因子(COMP)和立地條件(SITE)的函數。在建模前采用逐步回歸法從初始變量林木大小(樹高)、競爭因子(樹冠競爭因子、大于對象木的林分斷面積和、光截留)和立地因子(坡度、坡向、坡位、海拔、林分優勢高)中選擇最優變量，并在建模時直接采用最終變量。

各變量公式及解釋如下：

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

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

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

(2)

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

SIZE=β4×ln(HT)。

(3)

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

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

(4)

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

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

LI=1-e-K×LAI。

(5)

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

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

LAI =1.345 + 0.501W,

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

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

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

2) 幼樹樹高生長貝葉斯模型 本研究構建的油松幼樹樹高生長貝葉斯模型如下：

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

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

(6)

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

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

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

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

(7)

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

(8)

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

一階敏感性效應為：

(9)

二階敏感性效應為：

(10)

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

(11)

(12)

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

2.2.5 模型檢驗 1) 交叉檢驗法 將調查數據隨機等分為兩部分，第一部分作為訓練集進行建模，第二部分為驗證數據，同時再以第二部分數據作為訓練集進行建模，第一部分為驗證數據檢驗模型預測精度，通過均方根誤差(RMSE)表示。計算公式如下：

(13)

圖2 油松幼樹5年樹高生長預測的95%不確定性區間Fig.2 95% credible intervals of 5-year height increment model due to parameter uncertainty and total uncertaintyx軸表示每株油松幼樹樣本，黑色圓點表示每株油松幼樹的觀測值，灰色陰影部分表示預測的貝葉斯95%可信區間。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) 等效性檢驗法 基于回歸的等效性檢驗是由Robinson等(2005)提出的。不同于傳統的假設檢驗(差異性檢驗表明2個群體之間存在顯著性差異)，等效性檢驗假設模型與數據是不同的，檢驗的顯著性代表在某個研究區范圍內模型與數據是相似的，因此拒絕原假設是等效性檢驗的理想目標。基于回歸的等效性檢驗考慮均值差異和個體差異，通過以下6個步驟完成： (1) 預測值減去預測均值； (2) 確定一個截距與斜率可搖擺的等效性限制區域 (如±25%)； (3) 檢驗實測值與調整后預測值的線性回歸； (4) 通過計算2個單側置信區間得到線性回歸截距的等效性，同時與設定的等效性區間進行比較； (5) 通過計算2個單側置信區間得到回歸斜率的等效性，同時與設定的等效性區間進行比較； (6)根據置信區間是否落在等效區間內決定是否拒絕實測與預測的非相似性原假設。

3 結果與分析

3.1 模型整體表現

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

圖3 油松幼樹5年樹高生長模型的交叉檢驗Fig.3 Cross-validation of 5-year height increment model of young P. tabulaeformis交叉檢驗中調查數據隨機等分為兩部分，檢驗1以第一部分為建模數據，第二部分為驗證數據，而檢驗2以第二部分為建模數據，第一部分為驗證數據。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 油松幼樹樹高生長模型等效性檢驗結果①Tab.3 Summary of equivalence-based regression results for height model of young P. tabulaeformis

3.2 參數不確定性傳遞分析

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

表4 油松幼樹樹高生長模型的參數后驗分布Tab.4 Posterior probability distribution of parameters in the height model of young P. tabulaeformis

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

3.3 模型預測的不確定性

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

圖5 不同的HT、CCF、LI、SL條件下油松幼樹樹高5年生長預測Fig.5 Height increment on 5-year height increment of young P. tabulaeformis based on HT, CCF, LI, and SL圖中灰色陰影區域為貝葉斯預測95%可信區間，黑色垂直線為預測均值±1倍標準差，黑色圓點為MAP預測值。圖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 模型不確定性量化與分析

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

4.2 模型應用效果與預測分析

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

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

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

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

5 結論

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

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

3) 貝葉斯MAP估計和不確定性預測結果表明，隨著樹冠競爭因子、光截留和坡度增大，油松幼樹樹高5年生長量相應減小，表現出負相關關系； 隨著樹高增大其5年生長量相應增加，表現出正相關關系。結合參數的不確定性進一步分析得出，參數的不確定性越高，其控制的變量對模型預測結果的影響越不顯著。