氣候敏感的青岡櫟單木胸徑生長模型*

劉 帥 李建軍 卿東升 朱凱文 馬振燕

(中南林業科技大學 長沙 410004)

氣候變化是21世紀全球面臨的最嚴峻挑戰之一(Parmesanetal.，2003)。據政府間氣候變化專門委員會(IPCC，2001)預測，2100年全球平均溫度將比2000年升高1.4～5.8 ℃，約是20世紀升溫幅度的2～10倍。森林對溫度、降水和二氧化碳濃度改變通常較為敏感，未來氣候變化將不可避免地影響森林生長和收獲、樹種組成和分布、森林生產力及其各項功能(Battlesetal.，2008)，且這些影響會隨著森林類型、生長地點以及氣候場景不同而異(Oboiteetal.，2019)。因此，在氣候變化背景下，深入了解氣候對森林生長的影響有助于林業工作者應對未來氣候變化所帶來的挑戰(Temperlietal.，2012)。

基于森林生長和收獲模型預估未來林木生長動態，一般將競爭、立地條件和林分結構作為協變量(Bigingetal.，1995；Peng，2000；Portéetal.，2002；Bravo-Oviedoetal.，2008)，很少考慮氣候因素，即假定森林生長過程中氣候條件不發生改變，因此氣候與林木生長之間的關系尚不明確；但一旦未來氣候發生改變，將給森林經營帶來很大的不確定性(Subedietal.，2011；2013)。近年來，氣候與林木生長的關系開始得到越來越多關注，一些研究已將氣候變量納入到林分或單木生長建模中(Wangetal.，2007；Cortinietal.，2011；Soniaetal.，2012；Nighetal.，2012；Subedietal.，2011；2013；Sharmaetal.，2015；Aubry-Kientzetal.，2017)。Subedi等(2011；2013)研究發現，林木直徑生長與氣候變量之間存在復雜的非線性關系，年平均溫度和降水顯著影響加拿大東部北方森林中短葉松(Pinusbanksiana)和黑云杉(Piceamariana)的胸徑生長。Pokharel(2009)研究指出，年平均溫度是影響加拿大安大略省4種常見樹種斷面積生長的主要氣候因子。其他類似的研究還有很多，如Huang等(2011)基于廣義可加模型研究氣候變化對美國北部火炬松(Pinustaeda)的影響、余黎等(2014)探討氣候因子對長白落葉松(Larixolgensis)-云杉(Piceajezoensisvar.microsperma)冷杉(Abiesnephrolepis)林6種主要樹種單木胸徑生長的影響及不同樹種對氣候變化的響應等。可見，溫度和降水是影響林木直徑生長的關鍵氣候因子(Subedietal.，2011；2013；Sharmaetal.，2015；Oboiteetal.，2019)。

林木生長數據通常以單木、林分或區域進行嵌套或分組，一般具有縱向(多層)結構，且由于其獲取方式依賴于每隔一段時間對同一目標進行重復測量，數據往往存在較強的相關性，因此采用傳統回歸分析方法如最小二乘法可能導致模型參數有偏估計。而混合效應模型通過引入隨機效應描述組內和組間變異，為重復或嵌套數據提供了解決方案(Bates，1990；Cudecketal.，2002)。在森林生長建模中，混合效應模型相比傳統回歸模型具有更好的擬合效果和預測精度(Fangetal.，2001；Jordanetal.，2005；Yangetal.，2011；Subedietal.，2011；2013；Sharmaetal.，2015；Saudetal.，2019)。

第八次全國森林資源清查結果顯示，我國共有櫟類林約1 672萬hm2，占全國森林總面積的10.15%(Wangetal.，2019)。青岡櫟(Cyclobalanopsisglauca)是構成我國亞熱帶森林的主要成分之一，因具有良好的材質(如木質堅硬、耐久耐腐等)已成為我國南方地區最重要的用材樹種，且在水土保持、維持生物多樣性等方面發揮著重要的生態功能(梁貴等，2015)。目前，已有一些關于青岡櫟生長的研究，如劉洵等(2019)評價了湖南省天然櫟類林的立地質量，編制了天然櫟類林胸徑指數表；龍時勝等(2018)采用序類聚類方法建立了不同生長階段的青岡櫟生長方程。此外，也有關于青岡櫟生物量的研究(丁增發，2014)。但這些研究均未涉及氣候變化，關于青岡櫟生長過程中是否受氣候影響以及受氣候影響的程度如何依然不明確。鑒于此，本研究首先構建氣候敏感的青岡櫟混合效應模型，識別顯著影響青岡櫟胸徑生長的氣候變量，在此基礎上，量化和預測未來不同氣候場景下青岡櫟胸徑生長動態，以期為未來氣候變化條件下的青岡櫟林經營決策提供理論依據。

1 數據與方法

1.1 數據來源

青岡櫟單木生長數據來源于中南林業科技大學蘆頭林場(28°31′—28°38′N，113°51′—113°58′E)。該林場地處羅霄山脈北段，屬大陸性季風氣候，年均氣溫18.5 ℃，平均降水量1 450.8 mm，平均日照時數1 930 h。全場總面積5 308 hm2，地貌類型以中山、低山為主，平均坡度25°左右，土壤以紅壤和山地黃壤為主。

在林場內干擾較少的青岡櫟次生林中設置8塊固定樣地，樣地大小均為20 m×20 m。在樣地中選擇30株33～54年生、長勢良好、無病蟲害的青岡櫟，于2016年7月將其伐倒。采用橫剖法(周生祥等，1989)解析樹干，測量并記錄每個齡級的胸徑、樹高、材積以及樹齡對應的年度等，共獲得320組單木生長數據。限于篇幅考慮，本研究僅利用胸徑、樹齡數據探討青岡櫟生長。為便于后續處理，以樹號進行分組，組內按齡級排序。

氣候數據來源于ClimateAP氣候模型，該模型以無標度數據作為基線，根據緯度、經度和海拔計算指定位置的月、季和年度氣候變量，時間跨度為1901—2100年。通過對比亞太地區1 805個氣象站的觀測數據，驗證了氣候數據的準確性(Wangetal.，2015；2017)。本研究選擇與林木生長相關的常用氣候變量(Subedietal.，2011；2013；Sharmaetal.，2015；Saudetal.，2019；Oboiteetal.，2019)，包括年平均溫度(mean annual temperature，MAT)、最熱月均溫(mean warmest month temperature，MWMT)、最冷月均溫(mean coldest month temperature，MCMT)、MWMT與MCMT溫差(temperature difference between MWMT and MCMT，TD)、年平均降水(mean annual precipitation，MAP)、干燥指數(annual heat-moisture index，AHM)。采用ClimateAP氣候模型提取樹齡對應年度的氣候數據。林木生長和氣候數據統計如表1所示。

表1 林木生長和氣候數據統計

1.2 研究方法

1.2.1 基礎模型和再參數化模型 本研究通過樹齡將青岡櫟林木生長和氣候數據匯總成一個數據集，采用最小二乘方法(least squares method)擬合Mitscherlich、Gompertz、Korf、Richards和Logistic 5個候選指數生長方程，各生長方程的評價結果見附表1。根據赤池信息準則AIC(Akaike’s information criterion)、貝葉斯信息準則BIC(Bayesian information criterion)、對數似然Loglik(Log-likelihood)、殘差平方和RSS(residual sum of squares)、平均相對誤差MRE(mean relative error)、均方根誤差RMSE(root mean squared error)的評價結果，Mitscherlich生長方程性能最好。因此，本研究選擇Mitscherlich生長方程作為基礎模型。Mitscherlich生長方程及相關評價指標定義如下：

y=β0(1-e-β1t)；

(1)

(2)

(3)

(4)

式中：y為胸徑；t為樹齡；β0和β1分別為待估計的漸近線(asymptote)參數和生長速率(rate)參數；obji為實際測量值；esti為相應的估計值；n為樣本數。

在基礎模型中，林木生長完全依賴于樹齡(唯一協變量)。為了量化氣候變化對林木胸徑生長的影響，本研究將基礎模型中的參數(β0和β1)表示為氣候變量的線性組合，即再參數化(Wangetal.，2007；Bravo-Oviedoetal.，2008；Sharmaetal.，2015；Zangetal.，2016)。同時，為了避免模型過擬合及多重共線性，首先采用向后逐步回歸方法(backward stepwise regression)對變量進行初步篩選，然后計算剩余變量的方差膨脹因子(variance inflation factors, VIF)。將VIF最大的變量從模型中剔除，再重新擬合；重復該過程直至所有剩余變量的VIF均小于5(一般VIF低于5，可認為自變量之間無共線性)為止。

1.2.2 非線性混合效應模型 構建混合效應模型是為了處理重復和嵌套的測量數據，其一般形式如下：

yi=f(α,β)+εi。

(5)

其中，

β=λ+γ；

γ～N(0,D)；

εi～N(0,Ri)；

式中：yi為測量值向量；f(.)為協變量向量α的非線性函數；β為待估計的混合效應參數向量，通常由λ和γ2部分組成：λ為固定效應參數向量；γ為隨機效應參數向量，假設隨機效應γ服從正態分布；D為其方差-協方差矩陣；εi為誤差項，假設服從正態分布；Ri為誤差項方差-協方差矩陣，用于表示組內誤差；σ2為模型誤差方差值；Φi為描述組內異方差(heteroscedasticity)的對角矩陣；Γi為描述組內相關性的自相關結構。

以上計算過程均在R軟件(版本3.5.2)中完成。青岡櫟單木生長基礎模型和再參數化模型采用最小二乘法擬合，由R包的“stats”完成；青岡櫟胸徑生長混合效應模型以單株樣木為隨機效應，其建模和預測由R包的“NLME”完成。

2 結果與分析

2.1 再參數化模型

經逐步回歸和多重共線性篩選，模型保留氣候變量包括年均溫MAT、最熱月均溫MWMT和最冷月均溫MCMT。因此，原基礎模型中的參數β0和β1可以表示為上述氣候變量的線性組合：

β0=λ00+λ01MAT+λ02MWMT+λ03MCMT；

β1=λ10+λ11MAT+λ12MWMT+λ13MCMT。

(6)

式中：λ00～λ03和λ10～λ13為待估計參數。

然而，并非每個氣候變量都對青岡櫟胸徑生長具有顯著影響，因此需對變量對應的參數進行顯著性檢驗。每次剔除P最大的不顯著變量，再重新擬合模型，重復以上過程直至所有變量均顯著(P< 0.05)為止。經處理，含氣候變量且參數統計顯著的再參數化模型如下：

yij=β0(1-e-β1Tij)。

(7)

其中，

β0=λ00+λ01MCMTij；

β1=λ10+λ11MCMTij。

式中：yij、Tij和 MCMTij分別為第i株林木第j次測量時的胸徑、樹齡和氣候變量。

2.2 非線性混合效應模型

yij=β0(1-e-β1Tij)+εij。

(8)

其中，β0=(λ00+γ00)+(λ01+γ01)MCMTij；

β1=(λ10+γ10)+λ11MCMTij；

εij～N(0,Ri)。

如表2所示，漸近線參數項中的MCMT系數λ01為負，表明MCMT負相關于青岡櫟胸徑生長，MCMT升高可能導致青岡櫟胸徑生長量減少；而指數項中的MCMT系數λ11為正，將加速上述過程。

表2 非線性混合效應胸徑生長模型的參數估計

2.3 模型評價

在逐步完善青岡櫟生長建模過程中，本研究依次構建了基礎模型【式(1)】、再參數化模型【式(7)】和混合效應模型【式(8)】。為便于比較，表3列出了3個生長模型的性能評價指標。可以看出，再參數化模型由于引入氣候變量，除BIC外的各項評價指標均略優于基礎模型，表明再參數化模型可用于解釋氣候變化對林木生長的影響；但基礎模型和再參數化模型均采用傳統最小二乘法擬合，未考慮林木生長數據的縱向和分組特征，無法修正數據間的相關性和殘差異方差，因此二者的性能包括擬合優度和誤差水平均遠遜于混合效應模型。與再參數化模型相比，混合效應模型的擬合優度指標AIC從1 070.342 0降至582.217 8，BIC從1 087.745 0降至620.504 9，LogLik從-530.171 1升至-280.108 9；其誤差指標RSS從1 165.440 9大幅減至57.733 2，MRE由0.263 0降至0.139 1，RMSE從2.208 2顯著降至0.491 5。

表3 模型擬合統計比較

經典的回歸分析要求誤差滿足同方差性，否則可能影響模型擬合與模型檢驗。圖1a、b所示為基礎模型和再參數化模型的殘差分布，可見存在明顯異方差性，直觀體現在隨著胸徑增加，殘差分布逐漸呈現喇叭狀開口；但在混合效應模型中，異方差性被消除(圖1c)。由于混合效應模型本身控制異方差效果較好，故本研究沒有在誤差項中進一步指定方差函數和自相關結構。

圖1 3個生長模型的殘差

圖2所示為部分青岡櫟的胸徑-樹齡生長曲線。每幅子圖內均有2條曲線，其中，藍色曲線由再參數化模型(僅含固定效應)擬合而成，代表所有青岡櫟在生長期內的平均胸徑；紅色曲線由混合效應模型擬合而成，代表單株青岡櫟的胸徑生長。為了對比，每幅子圖中還標記有該株青岡櫟胸徑的實測值(散點)。可見，藍色曲線大多偏離于實測值，而紅色曲線則與實測值變化相吻合，更能反映林木個體生長情況，也體現了混合效應模型的低誤差水平。2條曲線差異是由單木水平(以單木為分組依據)的隨機效應決定的，通過引入單木隨機效應，可使得混合效應模型無論是誤差水平還是擬合優度相比基礎模型和再參數化模型均有實質性提升。

圖2 由固定效應和混合效應擬合的胸徑-樹齡生長曲線

2.4 模型預測

為了預測氣候變化下的青岡櫟胸徑生長，本研究采用3種典型濃度路徑(representative concentration pathways，RCP)作為未來氣候場景(Vuurenetal.，2011)，分別為RCP2.6、RCP4.5和RCP8.5，依次代表未來溫室氣體低、中、高3種排放標準(2.6、4.5和8.5 W·m-2)。所有氣候場景下的氣候數據均通過ClimateAP 氣候模型預測，由于ClimateAP氣候模型以30年為時間間隔批量產生氣候數據，因此將2011—2100年分為3個等間隔時間段(2011—2040年、2041—2070年和2071—2100年)，分別預測不同氣候背景下青岡櫟1～30年的生長過程。表4所示為3個時間段不同氣候場景下的氣候數據統計情況。為便于比較，當前氣候條件也作為一種氣候場景。

表4 不同氣候場景下2011—2100年氣候變量(MCMT)的預測值

采用非線性混合效應模型【式(8)】預測未來青岡櫟胸徑生長，結果如圖3所示。圖3分階段展示了2011—2100年不同氣候場景下青岡櫟胸徑-樹齡(1～30年)生長曲線，其中胸徑為所有單木預測值的平均值。2011—2040年，青岡櫟胸徑生長差異非常小(圖3a)；2041—2070年，青岡櫟胸徑生長隨著不同氣候場景開始出現明顯分化(圖3b)；2071—2100年，青岡櫟胸徑生長差異非常大(圖3c)。在樹齡相同情況下，3種典型濃度路徑尤其是RCP8.5場景下青岡櫟胸徑顯著低于氣候條件不變時，其生長隨著RCP濃度遞增依次減少。這表明，隨著碳排放量增加，未來氣候條件很可能不利于青岡櫟胸徑生長。氣候對青岡櫟胸徑生長的影響2071—2100年相比2041—2070年更顯著，當2100年樹齡30年時，RCP2.6、4.5和8.5場景下青岡櫟胸徑相比氣候條件不變時分別下降6.3%、15.6%和53.1%；而在2070年樹齡30年時，RCP2.6、4.5和8.5場景下青岡櫟胸徑相比氣候條件不變時僅下降3.7%、6.2%和21%。可見，未來隨著時間推移，青岡櫟胸徑生長受氣候因素的負面影響將更加強烈。

圖3 基于MCMT平均值的不同氣候場景下青岡櫟胸徑-樹齡生長曲線

3 討論

樹木生長是一個復雜的非線性過程，通常采用指數生長方程描述(Budhathokietal.，2008；Subedietal.，2013；Trasobaresetal.，2016；Timilsinaetal.，2013)。本研究中，Mitscherlich生長方程被選為基礎模型用于初步描述青岡櫟胸徑生長，其漸近線和生長速率參數被再參數化為氣候變量的線性組合，并采用混合效應建模方法擬合胸徑生長模型參數。結果表明，加入氣候變量的青岡櫟生長模型能夠響應氣候變化對林木生長的影響(Wangetal.，2007；Sharmaetal.，2015；Zangetal.，2016；Saudetal.，2019)，而隨機效應引入則進一步提升了生長模型的擬合優度和預測精度(Doradoetal.，2006；Timilsinaetal.，2013；Yangetal.，2013；Zhuetal.，2019)。

氣候被認為是林木生長的潛在驅動，Latta等(2010)指出，森林生長和生產力受環境條件特別是氣候條件的影響。一般而言，氣溫升高通過產生額外熱量促進林木生長(Monserudetal.，2008)，且溫暖的氣候條件還會延長樹木生長期，從而有利于林木徑向生長(Bronsonetal.，2009；Huangetal.，2010；2013)。但也有研究表明，氣溫升高引起的蒸散量增加導致土壤水分減少，可能會對樹木生長產生負面影響，如果降雨量跟不上氣溫上升速度，森林生產力反而可能下降(Soniaetal.，2012)。Subedi等(2013)研究發現，氣溫升高可能是導致北方針葉林中黑云杉胸徑生長減少的主要原因。Aubry-Kientz等(2017)指出，如果未來氣候變暖持續下去，內華達山脈中一些樹種的生長量可能會大幅下降。本研究亦發現，最冷月均溫與青岡櫟胸徑生長呈顯著負相關。

本研究選取3種典型濃度路徑(RCP2.6、4.5和8.5)作為未來氣候場景，預測2011—2100年青岡櫟胸徑生長。相比氣候條件不變時，高排放的RCP8.5對青岡櫟胸徑生長的不利影響更大，而低排放的RCP2.6對青岡櫟胸徑生長的負面影響相對較小，且這些影響在2071—2100年尤其顯著。結合模型預測結果，隨著溫室氣體排放增加，最冷月均溫對青岡櫟胸徑生長的負面效應還會加劇。類似地，Sharma等(2018)量化分析了氣候變化對加拿大安大略地區紅松(Pinusresinosa)林生產力的影響，結果發現，在RCP8.5和RCP2.6場景下紅松林林分高度均顯著下降，且隨不同地區和氣候差異而有所變化。Battles等(2008)認為，未來所有氣候場景都將減緩加利福尼亞內華達山脈的針葉混交林生長，并導致針葉樹種生產力顯著下降。

本研究中，最冷月均溫是唯一顯著影響青岡櫟胸徑生長的氣候變量，其他氣候變量如年平均溫度、最熱月均溫和年平均降水等在統計上均不顯著，因此暫時無法評估這些氣候變量對青岡櫟胸徑生長的影響；特別是降水，通常被認為是樹木生長的另一類關鍵氣候因素，當溫度適宜時，充沛的降雨可以促進樹木和森林生長(Wangetal.，2007；Sharmaetal.，2015；Zangetal.，2016)。本研究數據集中有2個與降水相關的氣候變量，即年平均降水和干燥指數，干燥指數結合了年平均溫度和年平均降水2個重要指標，非常適合描述一年中氣候干燥情況(Zhangetal.，2014)，但這2個與降水相關的氣候變量由于統計不顯著而被排除在生長模型之外。

林分結構、土壤和地形等有關變量有助于提高氣候變化影響林木生長的預測(Rubénetal.，2015；Mansoetal.，2015；Oboiteetal.，2019)。Lo等(2010)研究發現，只有氣候變量還不足以解釋林木胸徑生長變化，且胸徑生長與林分密度和相鄰木之間的競爭密切相關。然而，由于缺乏相應數據，本研究生長模型中未包含與林分結構、土壤和地形相關的協變量，如果將來收集到相關數據，可以考慮將這些變量納入到模型中。此外，由于本研究林木數據來源于同一樣地，由空間差異引起的氣候變化無法在生長模型中得到體現。同時，考慮到氣候變化強烈影響著森林生長、樹種組成和分布(Chiangetal.，2008；Schelleretal.，2008)，后續研究有必要擴大研究區域，收集更多與氣候相關的青岡櫟生長數據。

4 結論

全球氣候變暖正深刻影響著森林生態系統的各方面。本研究以我國南方典型闊葉樹種青岡櫟為研究對象，基于蘆頭林場青岡櫟解析木數據和氣候數據，構建青岡櫟單木混合效應模型，選取代表未來溫室氣體排放標準的典型濃度路徑RCP2.6、4.5和8.5作為氣候場景，預測青岡櫟在未來不同時期不同氣候場景下的生長趨勢。本研究構建的青岡櫟單木混合效應模型具有氣候敏感、統計可靠和預測有效等優點，3種典型濃度路徑尤其是RCP8.5場景下青岡櫟胸徑顯著低于氣候條件不變時。隨著時間推移和碳排放量增加，未來氣候條件很可能不利于青岡櫟生長。本研究對氣候變化影響青岡櫟生長作了有益探索，研究結果有助于林業工作者在經營實踐中應對未來氣候變化所帶來的挑戰。

附錄Appendix

附表1 5種常用生長方程的比較①

附表2 15個非線性混合效應生長模型的比較①