利用經(jīng)驗(yàn)-過程混合建模方法優(yōu)化華山松過程模型的參數(shù)*

薛海連 田相林 曹田健

(1.西北農(nóng)林科技大學(xué)理學(xué)院 楊凌 712100； 2.赫爾辛基大學(xué)林學(xué)系 赫爾辛基 FI-00014；3.西北農(nóng)林科技大學(xué)生態(tài)仿真優(yōu)化實(shí)驗(yàn)室 楊凌 712100)

森林生長模型是森林生態(tài)系統(tǒng)動(dòng)態(tài)模擬與森林經(jīng)營規(guī)劃的重要工具，依據(jù)模型構(gòu)建方式，可分為經(jīng)驗(yàn)?zāi)Ｐ?、過程模型和混合模型(Kimmins, 2004)。傳統(tǒng)的經(jīng)驗(yàn)?zāi)Ｐ徒Ｒ蕾嚧罅繕拥卣{(diào)查數(shù)據(jù)，采用統(tǒng)計(jì)回歸技術(shù)預(yù)估模型參數(shù)，準(zhǔn)確可靠地預(yù)測林分生長和收獲，但外推至建模數(shù)據(jù)之外的環(huán)境條件時(shí)應(yīng)用受到限制。相比之下，過程模型的參數(shù)通過生理學(xué)觀測獲得，模型形式基于理論性植物與環(huán)境的關(guān)系，可外推計(jì)算環(huán)境變化對(duì)樹木生長的影響(Weiskitteletal., 2011)，但其結(jié)構(gòu)復(fù)雜、參數(shù)繁多且參數(shù)值難以實(shí)測確定?；旌夏Ｐ筒捎冒虢?jīng)驗(yàn)半機(jī)理的建模方式，一方面在機(jī)理建模的基礎(chǔ)上結(jié)合部分經(jīng)驗(yàn)?zāi)K，另一方面對(duì)部分機(jī)理模塊進(jìn)行簡化，如基于生理學(xué)過程的3-PG(physiological principles in predicting growth)模型(Landsbergetal., 1997)和CABALA(CArbon BALAnce)模型(Battagliaetal., 2004)均是通過集成經(jīng)驗(yàn)的異速生長模塊來簡化樹木生長的機(jī)制模擬。雖然混合模型可簡化模型結(jié)構(gòu)且能夠降低建模數(shù)據(jù)需求，但仍包含動(dòng)輒數(shù)十或上百個(gè)參數(shù)。無論是過程模型還是混合模型，高維變量條件下的參數(shù)預(yù)估方法設(shè)計(jì)與選取一直是模型建模研究的焦點(diǎn)問題。

過程模型參數(shù)預(yù)估通常涉及樹木生理、結(jié)構(gòu)、功能性狀等信息，這些信息往往具有一定的不確定性，有些參數(shù)甚至無法測量，只能間接推斷。目前，過程模型參數(shù)預(yù)估方法主要包括同時(shí)擬合法(Siev?nenetal., 1993; Vanclayetal., 1997)、Monte Carlo法(Hornbergeretal., 1985)、貝葉斯法(van Oijenetal., 2005; Siviaetal., 2006)和試錯(cuò)法(M?kel?, 1997； 張廷龍等， 2011)等，其中，同時(shí)擬合法、Monte Carlo法和貝葉斯法均需先給定參數(shù)的參考值或可能的變化區(qū)間等先驗(yàn)信息，而后搜索計(jì)算似然值(M?kel?etal., 2000； M?kel?， 1988； van Oijenetal., 2005； Ewen， 2013)。如Minunno等(2019)研究CROBAS(a growth model based on crown and base dynamics)模型參數(shù)的不確定性時(shí)，基于芬蘭國家森林資源連續(xù)清查數(shù)據(jù)，采用貝葉斯法校正歐洲赤松(Pinussylvestris)、挪威云杉(Piceaabies)和垂枝樺(Betulapendula)的CROBAS模型； 然而，為收斂至參數(shù)的完整聯(lián)合后驗(yàn)分布，貝葉斯法在處理不可解析的高維模型時(shí)消耗大量時(shí)間，并需進(jìn)行收斂性的反復(fù)測試(van Oijenetal., 2005; Hartigetal., 2012)。試錯(cuò)法主要基于實(shí)測數(shù)據(jù)、擬合數(shù)據(jù)、文獻(xiàn)數(shù)據(jù)和專家經(jīng)驗(yàn)等信息系統(tǒng)測試模型中的參數(shù)，是多數(shù)過程模型進(jìn)行參數(shù)本地化的主要方法，如3-PG(Gonzalez-Beneckeetal., 2014； Xieetal., 2017)、BIOME-BGC(何麗鴻等， 2016； 曾慧卿等， 2008)和Tree-BGC(Koroletal., 1995)模型等。M?kel?(1997)基于碳平衡原理開發(fā)的歐洲赤松CROBAS原模型，部分參數(shù)也采用試錯(cuò)法確定，并與林分生長經(jīng)驗(yàn)?zāi)Ｐ皖A(yù)測結(jié)果進(jìn)行對(duì)照。試錯(cuò)法相對(duì)簡單易行，但變量過多時(shí)難以同時(shí)降低多個(gè)變量預(yù)測值與實(shí)測值的離差，且預(yù)測值累積誤差會(huì)隨著預(yù)測周期增加而逐漸增大。

模型預(yù)測準(zhǔn)確性決定模型的應(yīng)用價(jià)值，3-PG和CROBAS等過程模型之所以應(yīng)用較廣，也正是因其預(yù)測準(zhǔn)確性較高(Weiskitteletal., 2011)。森林生物統(tǒng)計(jì)模型的可靠性還常通過不同模型的結(jié)果對(duì)比來評(píng)價(jià)。相對(duì)于Monte Carlo法和貝葉斯法，參數(shù)優(yōu)化方法的優(yōu)點(diǎn)是在處理復(fù)雜模型參數(shù)化時(shí)的收斂速度更快。不同于經(jīng)典的非線性規(guī)劃，種群算法(population-based algorithms)無需設(shè)定決策變量的經(jīng)驗(yàn)初始值，對(duì)目標(biāo)函數(shù)無可導(dǎo)可微的要求，特別是種群算法中的差分演化(differential evolution, DE)算法(Stornetal., 1997)，結(jié)構(gòu)簡單、控制參數(shù)少、具有較強(qiáng)的全局收斂性和魯棒性，更適合于解決復(fù)雜環(huán)境中的優(yōu)化問題(盛萬興等， 2017)。目前，差分演化算法已被應(yīng)用于包括森林經(jīng)營規(guī)劃(Pukkala, 2009; Arias-Rodiletal., 2015)在內(nèi)的諸多領(lǐng)域(Lpez Cruzetal., 2003; Onwuboluetal., 2006)，但其是否適用于過程模型建模仍有待證明。

本研究針對(duì)過程模型普遍存在的結(jié)構(gòu)變量復(fù)雜、生理學(xué)數(shù)據(jù)獲取困難的問題，以華山松過程模型CROBAS-PA(PA=Pinusarmandii)的參數(shù)優(yōu)化為例，提出一種結(jié)合經(jīng)驗(yàn)?zāi)Ｐ团c過程模型的混合建模方法，探索在建模數(shù)據(jù)不足情況下預(yù)估復(fù)雜過程模型參數(shù)的有效途徑。

1 研究區(qū)概況

華山松林分生長數(shù)據(jù)源于秦嶺南坡中段的火地塘林區(qū)。該林區(qū)地處陜西省寧陜縣境內(nèi)(108°21′—108°39′E，33°18′—33°28′N), 海拔1 420～2 474 m，年均氣溫6～13 ℃，年降水量700～1 200 mm，屬北亞熱帶和暖溫帶過渡地帶。土壤主要為火成巖、變質(zhì)巖、石灰?guī)r母質(zhì)上發(fā)育的弱酸性山地棕壤，厚度50～70 cm?；鸬靥亮謪^(qū)華山松林面積占林地面積的22%(李悅黎等， 1993)，主要分布在海拔1 600～2 200 m之間。

2 研究方法

2.1 數(shù)據(jù)來源

經(jīng)驗(yàn)?zāi)Ｐ徒?shù)據(jù)來源于457塊角規(guī)繞測臨時(shí)樣地(角規(guī)系數(shù)為1)(表1)，包括火地塘教學(xué)試驗(yàn)林場1990(244塊)、2005(143塊)和2013(70塊)森林二類調(diào)查樣地。2014年采集該林區(qū)31塊20 m × 20 m固定樣地的每木檢尺數(shù)據(jù)(表1)，按地位級(jí)(Ⅰ級(jí)13塊，Ⅱ級(jí)14塊，Ⅲ級(jí)4塊)(吳恒等， 2015)分類，用于參數(shù)優(yōu)化后的過程模型檢驗(yàn)。

表1 華山松樣地概況Tab.1 Stand attributes of P. armandii plots

2.2 經(jīng)驗(yàn)?zāi)Ｐ图捌鋮?shù)預(yù)估

本研究建立的華山松經(jīng)驗(yàn)?zāi)Ｐ蚎UASSI 1.0系統(tǒng)由式(1)～(6)構(gòu)成。因缺乏優(yōu)勢木樹高數(shù)據(jù)，故選取基準(zhǔn)年齡時(shí)的林分平均高，即地位級(jí)指數(shù)SCI(site class index)作為立地質(zhì)量評(píng)價(jià)指標(biāo)。華山松林木材積與樹高連年生長量在林齡30年前已明顯下降(陳存根等， 1994)，本研究以30年為基準(zhǔn)年齡。式(1)為立地級(jí)導(dǎo)向曲線，采用Hossfeld方程，已知林分初始年齡和樹高后，可由該式推算SCI并折算不同地位級(jí)指數(shù)下的樹高生長過程(McDilletal., 1992)；式(2)為林分密度指數(shù)(Stand density index, SDI)模型(Reineke, 1933)，基于林分中立木度最高的部分樣地通過參數(shù)擬合得到；式(3)和(4)是林分胸高斷面積(G, m2·hm-2)模型和平均胸徑(D，cm)模型，屬于Schumacher模型形式(Schumacher, 1939)；林分蓄積量(V, m3·hm-2)模型(式6)為華山松材積方程(陳存根等， 1996)：

(1)

(2)

(3)

(4)

(5)

(6)

式中：H為樹高(m)；T為林齡(a)；a1、a2、a3、b1、b2、b3、b4、c1、c2、c3、c4和β均為系數(shù)；N為林分密度(trees·hm-2)。

本模型構(gòu)建數(shù)據(jù)來自表1描述的樣地，參數(shù)估計(jì)值與檢驗(yàn)見表2。

為保證檢驗(yàn)樣本獨(dú)立性，按調(diào)查年份將數(shù)據(jù)分為3組(表1)，采用交叉驗(yàn)證方法測試模型表現(xiàn)，最終以合并數(shù)據(jù)所獲得的模型參數(shù)(表2)模擬林分生長。本研究篩選樣地均為人工起源的華山松相對(duì)純林，林分內(nèi)部個(gè)體年齡差異基本在一個(gè)齡級(jí)范圍內(nèi)，按同齡林進(jìn)行建模。式(1)為樹高導(dǎo)向曲線，未考慮除年齡以外引起樹高數(shù)據(jù)變異的因素，導(dǎo)致較低的決定系數(shù)。而胸徑預(yù)測模型(式4)決定系數(shù)較低，是因?yàn)閮H考慮了年齡、立地和林分密度3個(gè)因素，多數(shù)樣地?zé)o法獲取撫育、間伐、擇伐的歷史信息，導(dǎo)致數(shù)據(jù)中較高比例變異未被解釋。

表2 華山松經(jīng)驗(yàn)?zāi)Ｐ蚎UASSI 1.0參數(shù)估計(jì)值Tab.2 Parameter estimates of the empirical model QUASSI 1.0 for P. armandii

2.3 過程模型

華山松過程模型CROBAS-PA的原模型為M?kel?(1997)開發(fā)的CROBAS，該模型在模擬單木各器官碳分配過程時(shí)提出3個(gè)假設(shè)(式7—9)： 1) 樹冠表面積(Ac, m2)與葉生物量(Wf, kg)異速相關(guān)； 2) 葉和細(xì)根生物量(Wr, kg)之比為常數(shù)； 3) 葉生物量和冠基處邊材面積之比為常數(shù)(αq, m2·kg-1)(Shinozakietal., 1964)。

(7)

Wr=αrWf;

(8)

Aq=αqWf,(q=s,b,t)。

(9)

式中：Aq為各部分的邊材面積;s、b、t分別代表樹干、樹枝、運(yùn)輸根；ξ′和z為參數(shù)；αr為細(xì)根與葉生物量比值。

模型鏈接了光合作用、呼吸作用、自然整枝、枯損等模塊。華山松過程模型CROBAS-PA以初始葉生物量Wf、初始第一活枝高度Hs、林分初始密度N0為輸入變量，計(jì)算模型的狀態(tài)變量(亦稱模型輸出變量)，包括林木各器官在20年間的邊材面積和生物量、冠長、第一活枝高度等。

本研究通過改變細(xì)根與葉生物量比值αr和細(xì)根衰老速率sr(表3，Sr=srWf，Sr為單樹細(xì)根每年的衰老量，單位為kg·a-1)反映立地條件變化，從而影響林木各器官的碳分配(M?kel?etal., 2016； M?kel?, 1997; Persson, 1980)，即以第Ⅱ地位級(jí)參數(shù)為基準(zhǔn)(αr=0.2，sr=1)，計(jì)算CROBAS-PA第Ⅰ地位級(jí)(αr=0.1，sr=0.5)和第Ⅲ地位級(jí)(αr=0.25，sr=1.5)的參數(shù)。該方法的依據(jù)是細(xì)根與葉生物比值αr隨立地質(zhì)量提高而變小(Vanninenetal., 1999)，即立地較差的林分會(huì)有更高比例的碳分配至細(xì)根，以提高養(yǎng)分水分吸收量。

2.4 參數(shù)優(yōu)化

參數(shù)優(yōu)化模型(式10—11)的目標(biāo)函數(shù)設(shè)為過程模型CROBAS-PA與經(jīng)驗(yàn)?zāi)Ｐ蚎UASSI 1.0在樹高和生物量預(yù)測上的離差。

續(xù)表3 Continued

min{x1,x2,…,x10|Z0}f=∑dRRMSEd；

(10)

s.t.X=(x1,x2,…,x10),xj∈[Lj,Uj],

j=1,2,…,10 。

(11)

式中：Z0指林分初始狀態(tài)；X=(x1,x2,…,x10)為決策變量x1、x2、…、x10所組成的向量，向量X中的10個(gè)元素是表3中的最后10個(gè)參數(shù)；Lj與Uj分別指xj的下界與上界；d為樹高或樹木各部分生物量；RRMSEd指各變量的相對(duì)均方根誤差(relative root mean square error, RRMSE)。

分別利用經(jīng)驗(yàn)?zāi)Ｐ蚎UASSI 1.0與過程模型CROBAS-PA模擬特定初始年齡、立地、密度條件下90塊華山松樣地生長過程，通過優(yōu)化過程模型參數(shù)，最小化二者間的差異。樣地具體信息為： 初始年齡20、30、40年，每個(gè)齡組中分別令地位級(jí)指數(shù)SCI=8、9、10、11、12、13 m，林分密度指數(shù)SDI=500、600、700、800、900。將初始年齡、SCI和SDI輸入全林模型后，計(jì)算90塊樣地的林分初始狀態(tài)及后續(xù)20年間的生長過程。20年間的葉、枝和根生物量變化可由CROBAS-PA計(jì)算，并與經(jīng)驗(yàn)生物量公式(陳存根等， 1996)對(duì)比。

參數(shù)優(yōu)化模型的決策變量選取過程模型CROBAS中基于生理學(xué)和物候?qū)W設(shè)定的10個(gè)參數(shù)，即樹葉衰老率、樹冠樹葉分形維數(shù)、比葉面積、最大光合速率、消光系數(shù)、葉表面積密度、自然整枝參數(shù)、樹枝邊材率、樹干邊材衰老率和樹枝邊材衰老率。約束條件為這10個(gè)決策變量的可行域(表3)。

采用Storn等(1997)提出的差分演化算法(DE)，該算法是基于“貪婪競爭”機(jī)制進(jìn)行全局尋優(yōu)的隨機(jī)算法，通過父代變異、雜交的方式產(chǎn)生子代，而后由競爭機(jī)制選擇最優(yōu)解。DE的變異操作充分利用其他個(gè)體信息，增加了種群多樣性，同時(shí)減少了單純個(gè)體變異帶來的隨機(jī)性和盲目性。利用最優(yōu)個(gè)體加隨機(jī)個(gè)體差分法的變異方法，既提高了搜索速度，也降低了陷入局部最優(yōu)的危險(xiǎn)；但是若種群分布密度過高，也可能使其陷入局部最優(yōu)。差分演化算法的初始群體從10個(gè)決策變量的可行域中隨機(jī)抽取，10個(gè)決策變量的值和表3中的前26個(gè)參數(shù)一起輸入CROBAS-PA模型進(jìn)行優(yōu)化計(jì)算，從中選取使目標(biāo)函數(shù)最小的一組參數(shù)值。通過以下6個(gè)步驟執(zhí)行優(yōu)化搜索。

1) 初始化。產(chǎn)生初始種群，種群個(gè)體Xi=(xi1,xi2,…,xijD),其中xij∈[Lj,Uj],i=1,2,…,Np,j=1,2,…,D，D為自變量的維數(shù)。種群大小Np=45，縮放因子F=0.6，雜交概率CR=0.8。迭代次數(shù)Gm=50。

2) 評(píng)價(jià)初始群體。計(jì)算初始群體中各個(gè)體的適應(yīng)值f(Xi)，并記錄最優(yōu)值與最優(yōu)個(gè)體。

3) 差分變異操作。在種群中隨機(jī)選取不同的個(gè)體Xa、Xc、Xd、Xe、Xf,進(jìn)行變異操作取得新的種群個(gè)體Vi：

rand/1Vi=Xi+F(Xa-Xc)；

rand/2Vi=Xa+F1(Xc-Xd)+F2(Xe-Xf)；

best/1Vi=Xb+F(Xa-Xc)；

best/2Vi=Xb+F1(Xa-Xc)+F2(Xd-Xe)；

current to best/1Vi=Xi+α(Xc-Xi)+F(Xa-Xc)。

式中：α為D維向量，其元素是[0,1]之間的隨機(jī)數(shù)；F1、F2均為縮放因子(Liuetal., 2010)。

4) 雜交操作(二項(xiàng)式雜交)。若αj[0,1]

5) 選擇操作。若f(Ui)

6) 檢驗(yàn)是否符合條件。如果當(dāng)前的迭代次數(shù)達(dá)到預(yù)定設(shè)置的最大迭代次數(shù)Gm，則停止迭代，輸出最優(yōu)解，否則轉(zhuǎn)到步驟3。

本研究在優(yōu)化算法選用時(shí)，比較了差分演化算法(DE)的6種變異方式(圖1)。其中rand/2變異方式收斂速度最慢，best/1收斂速度最快，但變異方式rand/2與best/1之間的最優(yōu)值竟相差5.8%。通過設(shè)置不同的初始種群，并測試其收斂結(jié)果，最后根據(jù)算法的穩(wěn)定性和收斂速度選用current to best/1變異方式。

圖1 算法DE的迭代過程Fig.1 Process of iteration of algorithms DE

2.5 模型檢驗(yàn)及敏感性分析

為檢驗(yàn)CROBAS-PA模型的準(zhǔn)確性，以31塊固定樣地(表1)作為一組獨(dú)立數(shù)據(jù)對(duì)參數(shù)優(yōu)化后的過程模型進(jìn)行驗(yàn)證。此外，設(shè)計(jì)不同立地(18種)、不同初始密度(5種)的90塊樣地(n=18×5=90)，分別采用CROBAS-PA和QUASSI 1.0預(yù)測樹高和胸徑生長過程，而后對(duì)二者在20年間的樹高和胸徑預(yù)測差異進(jìn)行統(tǒng)計(jì)學(xué)驗(yàn)證。使用經(jīng)驗(yàn)?zāi)Ｐ蚎UASSI 1.0生成的90塊虛擬樣地可以避免過程模型參數(shù)優(yōu)化因抽樣非均衡帶來的誤差。具體虛擬樣地包括5個(gè)林分密度指數(shù)(SDI=500、600、700、800、900)及3個(gè)地位級(jí)共18個(gè)地位級(jí)指數(shù)(第Ⅰ地位級(jí)SCI=12.3、12.7、13.1、13.5、13.9、14.3 m； 第Ⅱ地位級(jí)SCI=9.9、10.3、10.7、11.1、11.5、11.9 m； 第Ⅲ地位級(jí)SCI=7.7、8.1、8.5、8.9、9.3、9.7 m)。所用的統(tǒng)計(jì)量有平均誤差(mean error, ME)、平均絕對(duì)誤差(mean absolute error, MAE)和平均相對(duì)誤差(mean relative error, MRE)，3個(gè)統(tǒng)計(jì)量的絕對(duì)值越小，說明CROBAS-PA預(yù)測的樹高和胸徑越接近QUASSI 1.0 預(yù)測的樹高和胸徑，預(yù)測精度也越高。

本研究基于生理學(xué)意義選取12個(gè)參數(shù)(樹冠樹葉分形維數(shù)、消光系數(shù)、比葉面積、最大光合速率、樹葉衰老率、葉表面積密度、自然整枝參數(shù)、樹枝邊材率、樹干邊材衰老率、樹枝邊材衰老率、細(xì)根與葉生物量比值和細(xì)根衰老速率)進(jìn)行全局敏感度分析，該方法考慮了單個(gè)參數(shù)及多個(gè)參數(shù)間的相互作用對(duì)模型輸出變量的影響。

應(yīng)用Sobol法(Cukieretal., 1975)將過程模型全局靈敏度的目標(biāo)函數(shù)設(shè)為CROBAS-PA與QUASSI 1.0在預(yù)測樹高或胸徑變化上的Nash-Sutcliffe效率(NSE)(Nashetal., 1970)。NSE常用于復(fù)雜模型靈敏度分析的目標(biāo)函數(shù)(Yang, 2011; Nossentetal., 2011)。NSE的理論范圍為-∞到1，通常0

本研究用于靈敏度分析計(jì)算的樣地共36塊(n=3×4×3=36，由QUASSI 1.0模擬)，初始年齡分別設(shè)為20、30、40年，地位級(jí)指數(shù)SCI=6.5、8.5、10.5和12.5 m, 林分密度指數(shù)SDI=500、800、1 100。NSE具體計(jì)算公式如下(式12)：

(12)

3 結(jié)果與分析

3.1 過程模型參數(shù)本地化

3.2 過程模型與經(jīng)驗(yàn)?zāi)Ｐ捅容^

過程模型CROBAS-PA與經(jīng)驗(yàn)?zāi)Ｐ蚎UASSI 1.0的20～40年樹高生長預(yù)測對(duì)比(圖2)顯示，CROBAS-PA預(yù)測的樹高相較QUASSI 1.0低估0.1～1.0 m，且這種差異在地位級(jí)高的林分中更明顯。在初始密度為2 000和1 500株·hm-2的樣地中，CROBAS-PA預(yù)測的胸徑比QUASSI 1.0低估0～1.1 cm，而在初始密度1 000株·hm-2中，CROBAS-PA預(yù)測的胸徑比QUASSI 1.0高估0.6～1.2 cm(圖3)。

圖3 基于兩種模型的不同立地與林分初始密度條件下林分胸徑預(yù)測值Fig.3 Comparisons between two models for DBH developments under different site and initial density conditions

圖4 經(jīng)驗(yàn)生物量公式(陳存根等， 1996)與過程模型CROBAS-PA預(yù)測的葉生物量、枝生物量和根生物量Fig.4 Comparison between empirical biomass equations(Chen et al.，1996) and CROBAS-PA for foliage biomass, branch biomass, and root biomass predictions

圖5 第Ⅱ地位級(jí)不同初始密度下的林分密度、第一活枝高度和冠長率的比較Fig.5 Comparisons of stand density, crown base height and crown ratio between stands of different initial densities in site classⅡ

本研究以中等立地與密度條件的林分為例(圖4，第Ⅱ地位級(jí)、初始密度1 500株·hm-2)，比較CROBAS-PA和經(jīng)驗(yàn)生物量公式(陳存根等， 1996)預(yù)測的各器官生物量，結(jié)果表明，林齡在20～40年間，CROBAS-PA模型的葉生物量預(yù)測值偏低，枝生物量和根生物量的預(yù)測結(jié)果分別在林齡35年之前和34年之前偏低，但之后偏高。

經(jīng)驗(yàn)?zāi)Ｐ屯ㄟ^林分密度控制來描述競爭，過程模型可通過第一活枝高和冠長率的變化反映林木競爭(圖5)。林分枯損模塊是基于林分密度指數(shù)進(jìn)行構(gòu)建的(Reineke, 1933)即密度高于給定SDI值的林分會(huì)自然枯損。CROBAS-PA基于林冠結(jié)構(gòu)和光截留計(jì)算林冠內(nèi)的蔭蔽環(huán)境和自然整枝過程。密度較大的林分對(duì)空間和光照的競爭更激烈，第一活枝高度較高且冠長率較小(圖5)。

3.3 預(yù)測誤差與參數(shù)敏感性

固定樣地實(shí)測值與過程模型CROBAS-PA預(yù)測值的相關(guān)分析表明，CROBAS-PA模型的樹高和胸徑預(yù)測值與實(shí)測值基本呈線性相關(guān)(圖6)。即便31個(gè)用于對(duì)比的固定樣地林分狀況差別較大，如林齡從20到52年，林分密度從289株· hm-2到4 616株· hm-2，但樹高和胸徑的實(shí)測值和模型預(yù)測值的R2仍分別達(dá)到了0.714 8和0.643 3。

虛擬樣地20年間的生長統(tǒng)計(jì)檢驗(yàn)結(jié)果顯示，CROBAS-PA模型在樹高和胸徑預(yù)測上存在低估，但誤差較小(表4)。CROBAS-PA模型樹高預(yù)測值的平均誤差與平均絕對(duì)誤差均低于1.03 m，平均相對(duì)誤差小于5.59%，其中第Ⅰ地位級(jí)的預(yù)測誤差略高，第Ⅲ地位級(jí)的預(yù)測誤差則較低(表4)。CROBAS-PA胸徑預(yù)測值的平均誤差、平均絕對(duì)誤差和平均相對(duì)誤差分別小于0.36 cm、1.19 cm和2.59%，其中平均相對(duì)誤差最高(2.59%)，出現(xiàn)在第Ⅰ地位級(jí)； 最高平均誤差(0.36 cm)發(fā)生于第Ⅱ地位級(jí)； 而平均絕對(duì)誤差最高(1.19 cm)的為第Ⅲ地位級(jí)(表4)。

全局靈敏度Sobol分析表明(表5)，最大光合速率P、樹冠樹葉分形維數(shù)z、消光系數(shù)k和比葉面積An對(duì)樹高生長影響的總靈敏度系數(shù)排序最高，分別為0.444 3、0.361 0、0.312 9和0.301 0。胸徑的總靈敏度系數(shù)排序最前的4個(gè)參數(shù)與樹高一致。樹高和胸徑對(duì)各參數(shù)的一階靈敏度系數(shù)之和分別為0.646 6、0.588 5，均小于1，說明該模型是不可加的，即非線性的。

圖6 樹高、胸徑的實(shí)測值和預(yù)測值比較Fig.6 Comparison between observed and simulated values for tree height and DBH

表4 CROBAS-PA樹高、胸徑預(yù)測值的檢驗(yàn)結(jié)果Tab.4 Model validation of CROBAS-PA on tree height and DBH predictions

表5 樹高和胸徑的一階靈敏度系數(shù)和總靈敏度系數(shù)Tab.5 Sobol’s first-order indices and total-effect indices for height and DBH

4 討論

參數(shù)優(yōu)化后的過程模型CROBAS-PA可以實(shí)現(xiàn)華山松林木各器官生物量的預(yù)估，并預(yù)測不同立地與密度條件下林分胸徑和樹高的生長。靈敏度分析顯示胸徑和樹高生長對(duì)比葉面積An、最大光合速率P、消光系數(shù)k和樹冠樹葉分形維數(shù)z較為敏感。

CROBAS-PA模型的樹高、胸徑預(yù)測誤差與模型結(jié)構(gòu)簡化直接相關(guān)。Shcherbinina(2012)基于CROBAS預(yù)測短葉松(Pinusbanksiana)生長時(shí)，亦出現(xiàn)類似本研究中樹高被低估的結(jié)果，其原因是細(xì)根與葉生物量比值αr的變化對(duì)樹高生長的影響為間接過程(M?kel?, 1997; Shcherbinina, 2012)。另外，本研究顯示CROBAS-PA模型預(yù)測的樹高生長受林分初始狀態(tài)影響較小，但林分初始密度對(duì)CROBAS-PA模型預(yù)測的胸徑精度有一定影響，林分密度越大時(shí)胸徑預(yù)測偏差越大。由于本研究采用異速生長方程(陳存根等， 1996)約束華山松過程模型CROBAS-PA的生物量分配，簡化了林分密度對(duì)葉生物量的影響，而根據(jù)葉生物量計(jì)算胸徑生長，誤差由此而生(圖3)。

本研究的靈敏度分析證明比葉面積An、最大光合速率P、消光系數(shù)k對(duì)樹高和胸徑生長影響較高(表5)，與Ewen(2013)的研究結(jié)果一致。Weiskittel 等(2011)指出過程模型對(duì)光合作用、光截留方面的參數(shù)是非常敏感的。增加An、P、k意味著光合作用增強(qiáng)，進(jìn)而促進(jìn)了胸徑和樹高生長。本研究對(duì)葉表面積密度x的靈敏度分析結(jié)果表明，葉表面積密度x對(duì)胸徑和樹高生長影響較低，樹冠樹葉分形維數(shù)z的影響則高于葉表面積密度x(表5)，與歐洲赤松的CROBAS模型靈敏度分析結(jié)果類似(M?kel?, 1997)。因?yàn)闃涔跇淙~分形維數(shù)z是異速生長方程中的指數(shù)，其不僅影響葉生物量曲線，還影響該曲線的形狀。

本研究的CROBAS模型是與年齡無關(guān)(age-independent)的林分生長模型，因而可用于模擬復(fù)層或混交林生長動(dòng)態(tài)(Minunnoetal., 2019)。然而，這種與年齡無關(guān)的模型，通常無法描述樹木衰老引起的生理組織代謝速率和碳利用效率的改變。Vanninen等(1996)對(duì)歐洲赤松林木(18～212年生)生物量的研究顯示，細(xì)根與葉生物量在樹木整體生物量的相對(duì)占比隨年齡增加而減小，而與此同時(shí)，細(xì)根與葉生物量比值會(huì)隨年齡增大而增大，本研究中參數(shù)細(xì)根與葉生物量比值和細(xì)根衰老速率對(duì)胸徑與樹高預(yù)測的全局靈敏度值相對(duì)較低(表5)，因而未設(shè)定其與年齡的關(guān)系，但這種與年齡無關(guān)(age-independent)的模型簡化會(huì)導(dǎo)致對(duì)不同齡級(jí)林木細(xì)根生物量估計(jì)的誤差。因此，在構(gòu)建與年齡相關(guān)(age-dependent)的生長模型時(shí)，應(yīng)補(bǔ)充考慮這種變化。

5 結(jié)論

本研究通過參數(shù)優(yōu)化模型調(diào)試CROBAS過程模型參數(shù)，基于經(jīng)驗(yàn)?zāi)Ｐ秃瓦^程模型的混合建模方法構(gòu)建了華山松過程模型CROBAS-PA。研究表明，經(jīng)過生物量異速生長方程和全林經(jīng)驗(yàn)?zāi)Ｐ图s束的華山松過程模型CROBAS-PA，可準(zhǔn)確預(yù)測林分樹高和胸徑生長及林木各器官中的生物量分配。全局靈敏度分析顯示，在CROBAS-PA模型參數(shù)中，比葉面積An、最大光合速率P、消光系數(shù)k和樹冠樹葉分形維數(shù)z的敏感度較高。研究證明，在有限實(shí)測數(shù)據(jù)和文獻(xiàn)信息條件下，可在林分水平上通過經(jīng)驗(yàn)?zāi)Ｐ图s束保證過程模型預(yù)測的準(zhǔn)確性。這意味著對(duì)于數(shù)據(jù)缺失情境下的復(fù)雜過程模型的參數(shù)預(yù)估，這種經(jīng)驗(yàn)-過程混合建模方法有著特定的優(yōu)勢與應(yīng)用潛力。