競爭和氣候及其交互作用對杉木人工林胸徑生長的影響*

臧 顥 劉洪生 黃錦程 張祖棟 歐陽勛志 寧金魁

1.江西農業大學林學院 南昌 330045； 2.崇義縣林業局 贛州 341300)

量化樹木生長和氣候的關系有助于評價氣候變化對森林生長的影響(韓士杰等， 2016)。一般認為，氣候因子如溫度和降水等會影響樹木的光合作用、呼吸作用和蒸騰作用，對樹木生長具有較強的驅動作用(Toledoetal.， 2011； 余黎等， 2014)，然而樹木生長和氣候的關系也會受到競爭的調節作用(Ettingeretal.， 2013； Fernndez-de-Uaetal.， 2015； Zhangetal.， 2015； Fordetal.， 2017)。競爭作為評價相鄰樹木之間擁擠程度的指標，是種群發展的主要因子，對林木個體生存、生長、繁殖及林分的水平分布和垂直分布有著重要影響(Das， 2012； 刁軍等， 2013)。競爭反映的是樹木之間相互爭奪資源而產生的壓力(Pretzschetal.， 2010)，而氣候因子中如降水等正是樹木爭奪的資源，因此，競爭可能會改變樹木生長對氣候的響應(Clarketal.， 2011； Xiangetal.， 2016)，即競爭和氣候對樹木生長存在交互作用(Clarketal.， 2014； Teietal.， 2014； Fordetal.， 2017)。Magruder等(2013)對美國密歇根州的紅松(Pinusresinosa)進行間伐試驗時發現，氣候對森林生長的影響會受間伐強度的制約。Ford等(2017)在研究美國華盛頓州的山地森林胸徑生長時發現，競爭會改變氣候對胸徑生長的影響。目前，盡管在單木和林分生長模型中同時考慮競爭和氣候因素的研究眾多(Bravo-Oviedoetal.， 2008； Crookstonetal.， 2010； Nunesetal.， 2011)，但探討競爭和氣候對樹木生長是否存在相互作用的報道仍較少。胸徑是反映樹木生長狀況的最重要指標(Sunetal.， 2019； 陳國棟等， 2020)，探討競爭和氣候交互作用對胸徑生長的影響有助于研究氣候變化背景下單木材積和生物量的估計及立地生產力的評價等。

杉木(Cunninghamialanceolata)是我國南方最重要的速生用材樹種之一，據第九次全國森林資源清查結果(國家林業和草原局， 2019)，杉木人工林總面積現已達990.20萬hm2、蓄積達7.55億m3，占全國人工林總面積的17.33%、總蓄積的22.30%。我國是一個木材資源嚴重缺乏的國家，木材對外依存度高(許傳德等， 2015)，在氣候變化和木材短缺的壓力下，發展人工用材林已成為保障國家生態安全和解決木材供需矛盾的戰略選擇。鑒于杉木在全國用材林中的重要地位，探索競爭和氣候對杉木人工林生長的交互作用，對在氣候變化背景下科學經營人工林、實現林業“雙增”目標具有重要意義。

鑒于此，本研究基于杉木人工林固定樣地數據，采用潛在生長量修正法構建包含競爭指標和氣候因子的胸徑生長模型，并嘗試分析競爭和氣候及其交互作用對杉木人工林胸徑生長的影響，以期為氣候變化背景下模擬撫育間伐、擇伐后保留木的生長變化奠定基礎，為森林適應性經營中科學合理地對杉木人工林進行間伐、擇伐提供依據。

1 研究區概況、數據與研究方法

1.1 研究區概況

研究區位于江西省贛州市的南康區(25°28′—26°14′N，114°29′—114°55′E)、崇義縣(25°24′—25°55′N，113°55′—114°38′E)和上猶縣(25°42′—26°10′N，114°01′—114°40′E)，羅霄山脈與大庚嶺相遇的南端，為諸廣山余脈的一部分。屬亞熱帶季風氣候，年均氣溫20.4 ℃，年均降水量2 265.5 mm，以低山、丘陵地貌為主，土壤主要為黃壤、紅壤等。研究區種植有大量人工針葉林，主要樹種有杉木、馬尾松(Pinusmassoniana)、木荷(Schimasuperba)、南酸棗(Choerospondiasaxillaris)、絲栗栲(Castanopsisfargesii)、米櫧(Castanopsiscarlesii)、米錐(Castanopsischinensis)等。

1.2 數據資料

1.2.1 固定樣地數據 固定樣地數據來自南康區、崇義縣和上猶縣39塊杉木人工林樣地。樣地布設時間為2011—2015年，每塊樣地進行過1～2次復測，調查間隔期1年。主要調查樣地的地理坐標、地形、坡度、海拔、土壤類型和種植時間，并對胸徑5 cm以上的樹木進行每木檢尺，記錄樹種，測量胸徑、樹木相對位置等因子，其中， 20塊樣地沒有調查樹高。依據胸徑測量數據，計算林木2次調查間的胸徑生長量，共計6 181個胸徑生長量觀測值。樣地調查因子的詳細統計量見表1，樣地分布示意見圖1。

表1 數據統計①

圖1 樣地分布示意

1.2.2 氣候數據 氣候數據依據地理坐標和海拔從Wang等(2017)編寫的ClimateAP軟件中提取，空間分辨率為4 km×4 km。考慮到固定樣地的布設和復測時間多在6月，為模擬氣候對胸徑生長可能存在的影響，將每年7月至次年6月視為調查間隔期，提取調查間隔期內的月值數據，匯總計算各樣地調查間隔期的平均溫度、最高溫度、最低溫度、降水量和大于5 ℃的積溫。詳細統計量見表1。

1.3 研究方法

潛在生長量修正法是構建單木生長模型的一種常用方法(孟憲宇， 2006)，其基本思路是： 模擬林木的潛在生長量(無競爭壓力時的生長量)，利用反映林木所受競爭壓力的修正函數對潛在生長量進行調整和修正，從而得到林木實際生長量，用數學形式可以表示為：

id=idpot×M(CI)。

(1)

式中： id為胸徑生長量； idpot為胸徑潛在生長量；M(CI)為反映林木所受競爭壓力的修正函數，且滿足0≤M(CI)≤1。

1.3.1 潛在生長量模擬 模擬潛在生長量的步驟如下：

1) 選擇基礎方程，以氣候和地形因子反映立地質量。考慮到林木潛在生長量受立地質量影響(劉洋等， 2012； 王冬至等， 2015； 雷相東等， 2018)，而立地質量可用環境因子評價，如氣候(張海平等， 2017； Jiangetal.， 2015； Shenetal.， 2015)和地形因子(杜紀山， 1999)，故本研究嘗試以氣候和地形因子反映立地質量。綜合數據分析和以往研究成果(Pretzschetal.， 2010； Weiskitteletal.， 2011； Burkhartetal.， 2019)，選用柯列爾方程(斯瓦洛夫， 1983； 孟憲宇， 2006)作為潛在生長量的基礎方程，采用再參數化方法，將柯列爾方程中的參數表示為關于5個氣候因子(調查間隔期的平均溫度、最高溫度、最低溫度、降水量和大于5 ℃的積溫)和2個地形因子(海拔和坡度)的函數，通過最小二乘法求解參數，并剔除不顯著因子(顯著性水平取0.05)，從而得到含環境因子的生長方程。柯列爾方程的數學形式如下：

y=a0xa1e-a2x。

(2)

式中：y和x分別為因變量和自變量；a0、a1和a2為模型參數。

2) 利用分位數回歸模擬含環境因子的生長方程。分位數回歸模型可擬合因變量任意分位數的預測值，并避免選取樣點的主觀因素，常應用于最大化和最小化問題(高慧淋等， 2016； Haoetal.， 2007)。潛在生長量，即理論生長最大值，可理解為某一特定樹種在特定立地條件下無競爭壓力的疏開木生長量(孟憲宇， 2006； Pretzschetal.， 2010)。本研究嘗試利用分位數回歸模擬含環境因子的生長方程，并將預測結果作為潛在生長量，分位點τ取0.90、0.95和0.99。模型形式如下：

idpot=(a00τ+∑a0jτEj)D(a10τ+∑a1jτEj)e-(a20τ+∑a2jτEj)D。

(3)

式中：D為胸徑；Ej為第j個環境因子；a00τ、a0jτ、a10τ、a1jτ、a20τ和a2jτ為第τ分位點模型參數。

第τ分位點模型參數求解可通過下式計算得到(Koenkeretal.， 1978)：

(4)

1.3.2 修正函數構建 為探討競爭對樹木生長的影響，本研究選擇6個備用競爭指標用于構建修正函數，以反映胸徑生長對競爭的響應： 1) 林分密度指標，包括公頃株數N、公頃斷面積BA和林分密度指數SDI； 2) 單木競爭指標，包括相對胸徑Rd、大于對象木的斷面積之和BAL、簡單競爭指數Hegyi(孟憲宇， 2006； Weiskitteletal.， 2011； Zhouetal.， 2019)。部分競爭指標計算公式如下：

(5)

(6)

(7)

(8)

式中：D0為標準平均直徑，本研究取10 cm(孟憲宇， 2006)；Dg為林分平均胸徑；Di和gi分別為對象木的胸徑和斷面積； Area為樣地面積；gk為樣地中大于對象木的第k株樣木斷面積；Dj為第j株競爭木胸徑； Distij為第j株競爭木與對象木的距離；m為競爭木株數。

采用3種方法確定競爭木(湯孟平等， 2007； Kuehneetal.， 2019)： 1) 以對象木為中心，將與對象木最近的1～10 株相鄰木作為競爭木； 2) 以對象木為中心，以1～10 m為半徑畫圓，將各半徑內所有非對象木視為競爭木； 3) 選用Voronoi圖確定競爭木。考慮到邊緣效應的影響，運用8鄰體平移法計算Hegyi(申瀚文等， 2012； Ripley， 1979)。

修正函數M(CI)采用指數函數ef(CI)形式。為求得修正函數中的參數，同時考慮到個別為0 cm的胸徑生長量無法取對數，參考以往研究成果(李春明， 2012； Pretzschetal.， 2010)，將式(1)整理為ln(id/idpot+1)=f(CI)形式，以避免出現負值，f(CI)為關于競爭指標的線性函數。為反映競爭和氣候可能存在的交互作用，在競爭指標構建的修正函數中添加競爭和氣候交互作用項，并與不添加競爭和氣候交互作用項的函數進行對比。不考慮交互作用和考慮交互作用的修正函數參數如下：

(9)

(10)

式中： CI為競爭指標；CLi為第i個氣候因子；b0、b1、ci為模型參數。

考慮到嵌套數據結構可能帶來的變量獨立性缺失問題(李春明， 2012； 彭娓等， 2018)，基于選定的修正函數形式，采用混合效應模型進行重建。混合效應模型一般形式為：

Yi=f(β,ui,Xi)+εi。

(11)

式中：Yi和Xi分別為因變量向量和自變量向量；εi為誤差向量；β和ui分別為固定效應參數向量和隨機效應參數向量，且εi～N(0,Ri)、ui～N(0,Ψ)，Ri和Ψ分別為樣地內的方差協方差矩陣和隨機效應參數的方差協方差矩陣。本研究中，隨機效應參數的方差協方差結構采用廣義正定矩陣形式。

對混合效應模型進行檢驗時，未包括在建模樣本中的樣地須計算其隨機效應參數。通過泰勒公式將已有非線性混合效應模型線性化，可以推導出隨機效應參數的經驗貝葉斯估計值(Davidianetal.， 1995)，公式如下：

(12)

1.3.3 模型評價 采用交叉驗證法評價模型，即以全部39塊樣地所有杉木的單木胸徑生長數據為基礎，每次將其中38塊樣地數據作為建模樣本，另外1塊樣地數據作為檢驗樣本進行模型評價，重復39次。使用平均絕對誤差(mean absolute error，MAE)、平均相對誤差絕對值(relative mean absolute error，RMAE)和平均預估誤差(mean predicted error，MPE)評價模型估計精度，計算公式(曾偉生等， 1999； 2011； Ouetal.， 2019)如下：

(13)

(14)

(15)

2 結果與分析

2.1 潛在生長量與環境因子的關系

采用再參數化方法擬合得到含環境因子的潛在生長量方程如下：

idpot=(a00+a01Tmin+a02PPT+

a03Elevation)Da1e-a2D

(16)

式中：Tmin、PPT和Elevation分別為調查間隔期的最低溫度、降水量和樣地海拔。

剔除不顯著的環境因子后(顯著性水平=0.05)，最小二乘法的參數估計結果見表2。考慮到共線性可能對參數估計造成的影響，采用方差膨脹因子分析3個環境因子的共線性程度，最低溫度、降水量和海拔的方差膨脹因子分別為4.35、4.59和4.84，即可認為最低溫度、降水量和海拔之間不存在強共線性。

基于式(16)，利用分位數回歸(τ取0.90、0.95和0.99)模擬胸徑潛在生長量，參數估計結果見表2。調查間隔期的最低溫度、降水量和海拔對杉木人工林胸徑潛在生長量存在顯著影響，且最低溫度和降水量與潛在生長量呈正相關關系，海拔則呈負相關關系。

2.2 競爭對胸徑生長的影響

不同競爭指標構建的修正函數對胸徑生長量的擬合效果見表3。可以發現，含單木競爭指標的模型精度優于含林分密度指標的模型，與距離無關的單木競爭指標的生長模型估計精度略優于與距離有關的單木競爭指標的生長模型，含BAL的模型評價指標最優，其中，MAE降低0.71%～20.07%，RMAE降低12.52%～33.63%，MPE降低3.21%～25.16%。

表2 式(16)的參數估計

表3 含不同競爭指標的模型評價①

圖2 2種方法計算的Hegyi競爭指數對胸徑生長量的影響

通過對比不同分位數擬合杉木人工林胸徑潛在生長量對模型精度的影響(圖2)發現，分位數取0.90時，模型評價指標表現最差，隨著分位數增加，MAE、RMAE和MPE逐漸降低，僅分位數取0.95時含BA和Rd模型的MPE大于分位數取0.90時的模型。對比3種分位數下所有競爭指標模型，總體來看，分位數取0.99時含BAL的模型評價指標最優，故將其作為基礎模型進一步分析競爭和氣候交互作用對杉木人工林胸徑潛在生長量的影響。

2.3 競爭和氣候交互作用對胸徑生長的影響

結合上述分析，以分位數取0.99時含BAL的模型為基礎模型，進一步在修正函數中添加BAL與5個氣候因子(調查間隔期的平均溫度、最高溫度、最低溫度、降水量和大于5 ℃的積溫)的交互作用項，各模型的評價指標見表4。與不含交互作用項的模型(表3)對比，除RMAE外，考慮交互作用項后模型的評價指標更優，其中，MAE降低0.60%～18.69%，MPE降低0.12%～9.72%。考慮BAL與平均溫度、最高溫度、最低溫度交互作用的模型RMAE和MPE最小，即同時考慮BAL與平均溫度、最高溫度、最低溫度交互作用的模型最優。

表4 含不同交互作用項的模型評價

此外，考慮到嵌套數據結構可能帶來的變量獨立性缺失問題，采用混合效應模型重新模擬最優模型的修正函數部分。由于部分樣地僅有2期調查數據，即只有1期生長量數據，故未考慮單木水平的隨機效應。最終構建的胸徑生長模型形式如下：

(17)

式中： idijk為第i塊樣地中第j株樣木第k次調查間隔期的胸徑實測生長量； BALijk和Dijk分別為第i塊樣地中第j株樣木第k次調查間隔的期初BAL和胸徑；Tmeanijk、Tmaxijk、Tminijk和PPTijk分別為第i塊樣地中第j株樣木第k次調查間隔期的平均溫度、最高溫度、最低溫度和降水量； Elevationi為第i塊樣地海拔；a00τ、a10τ、a20τ、a0jτ、a1jτ、a2jτ為第τ分位數回歸模型參數(τ取0.99)；b0、b1、b2、b3、b4為固定效應參數；βi為樣地水平的隨機效應參數；εijk為誤差項。

式(17)中潛在生長量參數估計結果見表2，修正因子參數估計結果見表5。

添加隨機效應參數后，模型MAE=0.071 1 cm，RMAE=20.37%，MPE=4.90%。相比無隨機效應參數模型(表4)，MAE提升35.71%，RMAE降低4.36%，MPE降低37.26%。

基于最終胸徑生長模型，各徑階殘差分布如圖3所示。總體來看，該模型能夠較好預測各徑階的胸徑生長量，但徑階小于18 cm時，對部分林木實際生長量存在低估情況； 而徑階大于20 cm時，胸徑估計結果的精度和穩定性更好。

表5 最終模型修正函數的參數估計①

基于最終模型固定效應參數，不同環境條件下胸徑潛在生長量和修正函數的變化如圖4、5所示。就潛在生長量而言，胸徑20 cm時，潛在生長量達到最大值，且因環境條件改變潛在生長量的最大值也隨之變化，其中，調查間隔期的最低溫度影響最大，降水量影響最小。就修正因子而言，隨著林木所受競爭壓力不斷增加(BAL增加)，修正函數呈遞減趨勢，且在BAL>8 m2·hm-2時，減小速度變緩，并逐漸趨于穩定。此外，BAL在0～30 m2·hm-2之間時，平均溫度降低會導致修正函數減小，BAL在0～23 m2·hm-2之間時，最高溫度增加會導致修正函數減小，BAL在0～15 m2·hm-2之間時，最低溫度增加會導致修正函數減小，即平均溫度降低、最高溫度增加和最低溫度增加均會加劇競爭對胸徑生長的影響，且該影響隨著BAL增加呈先加強后減弱的趨勢，BAL≥30 m2·hm-2后不再有明顯影響，因此可認為： 對林分中承受競爭壓力較小(BAL接近0 m2·hm-2)和競爭壓力較大(BAL≥30 m2·hm-2)的林木而言，競爭是對胸徑潛在生長量修正的主要因子，氣候影響不大，而對林分中BAL在0～30 m2·hm-2之間的林木而言，競爭對胸徑潛在生長量的修正效應會受氣候影響。

圖3 各徑階殘差分布

圖4 不同環境條件下胸徑潛在生長量的預測

圖5 不同環境條件下修正函數的平均響應

3 討論

潛在生長量修正法是構建單木生長模型的一種常用方法，其關鍵在于構建林木潛在生長量方程及對潛在生長量進行調整的修正函數。理論上，潛在生長量方程應基于始終處于無競爭和外界壓力情況下自由生長的疏開木調查數據，但這一點在實際中較難滿足。分位數回歸模型可擬合因變量任意分位數的預測值，其在國內外林業研究中已有應用，且在最大化問題如計算最大密度限上應用較多(Zhangetal.， 2005； Xueetal.， 2015； 高慧淋等， 2016)，因此，采用分位數回歸模型模擬胸徑生長量時，通過設定高分位點(如0.99)即可得到該分位點的胸徑生長量預測值，并可將其視為胸徑潛在生長量。從本研究結果看，最終胸徑生長量模型式(17)估計精度良好，因此可認為高分位點對含環境因子的胸徑潛在生長量描述較好。Pretzsch等(2010)利用分位數回歸模型模擬中歐4個樹種的潛在生長量，結果發現，該模型能較好反映出不同立地條件下的胸徑潛在生長量。

因部分樣地沒有調查樹高，因此本研究在構建胸徑生長模型時，用環境因子描述立地質量對林木潛在生長量的影響，結果顯示，調查間隔期的最低溫度、降雨量和樣地海拔對杉木人工林胸徑潛在生長量具有顯著影響，且均呈正相關。很多學者也采用氣候和地形等環境因子反映立地質量，如杜紀山(1999)在構建落葉松(Larixgmelinii)胸徑生長模型時，使用海拔、坡率(坡度的正切值)、坡向及坡率與坡向之間的組合項反映立地質量，并發現坡度的效果最優。Ou等(2019)在模擬落葉松-云杉(Piceaasperata)-冷杉(Abiesfabri)混交林胸徑生長量時，考慮坡度、海拔、坡向3個地形因子及44個氣候因子，結果發現盡管地形因子能提升模型精度，但效果弱于氣候變量。

競爭是樹木生長的主要影響因素之一(彭娓等， 2018)，因此競爭指標常作為自變量用于胸徑生長模型。依據是否反映出單木的競爭狀況差異，可將競爭指標分為反映林分內單木平均競爭狀況的林分密度指標和描述各單木不同競爭狀況的單木競爭指標(孟憲宇， 2006)。本研究中，通過比較6個競爭指標對胸徑生長量的模擬后發現，描述林木競爭狀況差異的單木競爭指標明顯優于反映林木平均競爭情況的林分密度指標，其中大于對象木的斷面積之和擬合效果最好。盡管有研究顯示樹木胸徑生長還與其周圍樹木的距離具有很高相關性(劉洋等， 2012； 房曉娜等， 2015； Coomesetal.， 2007)，但也有學者發現與距離無關的競爭指標在擬合單木胸徑生長時效果更好(Kuehneetal.， 2019)，這可能與研究對象的差異和競爭指標的計算有關。本研究中，樣地邊緣林木的競爭木運用8鄰體平移法確定，這也可能導致計算得到的Hegyi指數與林木的真實競爭狀況存在差異，因而降低了對胸徑生長量的擬合效果。另外，由于不需要考慮林木間的相對位置，因而在生產實際中，與距離無關的競爭指標應用也更為方便。很多學者在構建胸徑生長模型時會選用與距離無關的競爭指標，且估計精度均較高(李春明， 2012； Zhaoetal.， 2013； Fordetal.， 2017)，也有一些學者研究發現，林分密度指標和單木競爭指標均能提高胸徑生長模型的估計精度，可同時加入模型中，如呂勇等(1999)在模擬杉木人工林單木胸徑生長時，選取1個林分密度指標(相對植距)和1個單木競爭指標(相對優勢度)； 馬武等(2015)采用逐步回歸法，以方差膨脹因子和調整決定系數為依據，篩選用于構建蒙古櫟(Quercusmongolica)天然林胸徑生長模型的自變量，備選的7個競爭因子(林分斷面積、林分密度指數、大于對象木的樹木斷面積之和、相對胸徑、大于對象木的所有林木胸徑平方和、對象木直徑與林分中最大林木直徑之比以及郁閉度)中有1個林分密度指標(林分斷面積)和1個單木競爭指標(大于對象木的樹木斷面積之和)被選入最終模型中。

4 結論

為分析胸徑生長受競爭和氣候條件的影響，本研究基于江西省贛州市南康區、崇義縣和上猶縣杉木人工林固定樣地數據，采用潛在生長量修正法構建胸徑生長模型，結果發現，海拔、調查間隔期的最低溫度和降水量對胸徑潛在生長量具有顯著影響，胸徑20 cm時，潛在生長量達到最大值。競爭作為林木個體生存、生長的重要影響因子之一，也對模型擬合好壞有較大影響，且單木競爭指標估計精度優于林分密度指標，與距離無關的單木競爭指標估計精度略優于與距離有關的單木競爭指標。通過分析競爭和氣候交互作用對潛在生長量的修正效果發現，競爭和氣候交互作用可在一定程度上提高估計精度，但不同競爭指標和氣候因子構成的交互作用效果并不一致，其中，同時考慮平均溫度、最高溫度和最低溫度與大于對象木的斷面積之和構成的交互作用的模型預估精度最好，且平均溫度降低、最高溫度和最低溫度增加會加劇競爭對胸徑生長量的影響。在此基礎上，下一步可考慮探討在氣候變化下如何設計間伐、擇伐才能最大化經營目的。