基于貝葉斯模型平均法和逐步回歸法構建杉木單木胸徑生長模型*

魯樂樂 王 震 張雄清 張建國

(1.中國林業科學研究院林業研究所 國家林業和草原局林木培育重點實驗室 北京100091；2.南京林業大學南方現代林業協同創新中心 南京210037)

單木胸徑是估算立木材積、評價森林生產力的重要基礎因子(Lhotkaetal.， 2011)。影響單木胸徑生長變化的因子很多，如氣候、立地、林齡、林分密度、林分結構等(Clarketal.， 2014)，總體上可以分為內部因子(主要指林分因子)和外部因子(主要指氣候因子)兩大類(Zhangetal.， 2017)。

氣候因子是單木胸徑生長變化的重要驅動因子，已引起學者廣泛關注(Allenetal.， 2015； Hanewinkeletal.， 2012)。Henderson 等(2009)建立長葉松(Pinuspalustris)胸徑生長與氣候因子之間的關系模型，結果表明在所有氣候變量中，溫度是最弱的驅動因子。黎敬業等(2019)分析馬尾松(Pinusmassoniana)徑向生長與氣候因子的關系，結果發現高海拔地區馬尾松徑向生長對氣候因子年際波動敏感性較強，與生長季前冬季光溫條件和生長季內7月降雨呈正相關，生長-氣候關系在不同樣點間表現出較強一致性。王延芳等(2020)基于 Vaganov-Shashkin 模型模擬青海云杉(Piceacrassifolia)徑向生長，結果表明降雨對祁連山中部低海拔地區青海云杉徑向生長具有決定性作用。

對于林分因子，包括林分密度(Zhangetal.， 2015)、立地(Aragaoetal.， 2009)、年齡(Fosteretal.， 2014)等對單木胸徑生長變化的驅動性也不小。通常來說密度越大，單木胸徑生長量越小，立地條件越好，單木胸徑生長量越大(DeRoseetal.， 2009； Zhangetal.， 2015)。王冬至等(2017)分析華北落葉松(Larixprincipis-rupprechtii)-白樺(Betulaplatyphylla)針闊混交林胸徑生長，結果表明年齡和競爭指數與胸徑生長量呈負相關關系，胸徑和林分優勢高與胸徑生長量呈正相關關系。張海平(2017)構建基于氣象因子的天然白樺林單木胸徑生長模型，結果發現溫度和降雨有利于白樺生長。

綜上可知，氣候因子和林分因子對單木胸徑生長均有一定驅動性。在進行模型分析選擇變量時，一般大部分研究都是主觀指定包含某些變量的某個特定模型為最優模型，由此導致所建立模型的預測能力減弱，這是頻率統計模型的一個弱點。如逐步回歸法(stepwise regression, SR)根據變量顯著性(P)，選擇顯著的變量保留在模型中； 而實際上事先并不知道包含哪些變量的模型為最優模型，即模型本身存在不確定性(Rafteryetal.， 1991； Draper， 1995)。當基于研究數據建立回歸模型時，若忽略模型本身的不確性僅以某一特定模型進行分析，一方面會低估模型的不確定性，高估預測結果，限制模型的適用范圍； 另一方面會比正常情況下更趨向于拒絕無效假設，從而產生誤導性結果(Viallefontetal.， 2001； 張志杰等， 2007)。

貝葉斯模型平均法(Bayesian model averaging，BMA)是近些年發展較快的一種多因子和模型不確定性分析方法(Picardetal.， 2012)。與傳統分析方法不同，BMA并非僅指定一個最優模型，而是考慮模型空間內所有可能的模型以及各模型的后驗概率，將模型的后驗概率加權平均得到貝葉斯平均模型，并根據模型中變量估計值不等于 0 的后驗概率進行變量重要性排序，判斷各模型出現的可能性(羅劍鋒， 2003)。近年來，BMA 在林業領域也有些研究，如Bullock等(2007)成功地將 BMA 應用于林木胸徑分布中、Picard等(2012)利用 BMA分析了熱帶森林地面生物量、Zhang 等(2014)利用 BMA 構建了杉木(Cunninghamialanceolata)林分斷面積生長模型、Lu等(2019)利用BMA構建了枯損率模型等； 但是，在錯綜復雜的影響因子中如何更好地選擇因子，并構建不確定性胸徑生長模型還未見報道。

杉木是我國速生高產用材樹種之一，也是我國南方亞熱帶地區最重要的鄉土栽培樹種。據第九次全國森林資源清查結果，杉木人工林總造林面積達 990.2萬hm2，占全國十大主要造林樹種面積的27.23%，蓄積占32.57%，均居各大造林樹種之首。本研究擬探索杉木人工林單木胸徑生長量變化的驅動因子，比較不同驅動因子的重要性，構建不確定性單木胸徑生長模型，以期為杉木經營管理者科學經營管理杉木人工林提供參考。

1 研究區概況與數據調查

研究區位于福建省邵武市衛閩林場(117°43′E，27°05′N)，為杉木中心產區，地貌主要為低山、高丘，海拔250～700 m，坡度25°～35°。屬亞熱帶季風氣候，年均氣溫17.7 ℃，1月平均氣溫6.8 ℃，7月平均氣溫28 ℃，極端低溫-7.9 ℃，年均降雨量1 768 mm，年均相對濕度82%。

研究區杉木密度試驗林為1982年使用1年生苗木造林，采用隨機區組試驗設計，設置5種密度： A，2 m×3 m(1 667 株·hm-2)； B， 2 m×1.5 m(3 333 株·hm-2)； C，2 m×1 m(5 000 株·hm-2)； D，1 m×1.5 m(6 667 株·hm-2)；E， 1 m×1 m(10 000 株·hm-2)。每種密度重復3次，共設立樣地15塊，每塊樣地大小為20 m×30 m。造林后1984—1990年逐年調查，之后到2010年隔年調查，每次調查均在秋冬季進行。測量樹高時，每塊樣地隨機抽取50株樹，計算最高5株樹的平均高作為林分優勢木平均高。

2 研究方法

2.1 變量選擇

內部因子主要指單木大小因子、競爭因子和林分因子。單木大小因子包括胸徑(diameter at breast height, DBH，cm)、年齡(age,A,a)； 競爭因子包括胸高斷面積(basal area, BA，m2·hm-2)、大于對象木的斷面積和(sum of basal areas of trees larger than the subject tree, BAL，m2·hm-2)； 林分因子包括林分密度(number,N，trees·hm-2)、優勢木平均高(dominant height,Hd，m)、林分平方平均胸徑(quadratic mean diameter,Dg，cm)。內部因子統計值見表1。

外部因子主要指氣候因子，通過ClimateAP軟件獲得(Wangetal.， 2012)，包括年均氣溫(mean annual temperature，MAT，℃)、年均降雨量(mean annual precipitation，AP，mm)、最冷月平均溫度(mean coldest month temperature，MCMT，℃)、年均干旱指數(annual heat-moisture index，AHM)、低于0 ℃天數(degree-days below 0 ℃, DD0, d)、夏季平均最高溫度(summer mean maximum temperature，SMMT，℃)、冬季平均最低溫度(winter mean minimum temperature，WMMT，℃)、最熱月平均溫度(mean warmest month temperature，MWMT，℃)和春季平均氣溫(spring mean temperature，SMT，℃)。氣候因子統計值見表2。

表2 氣候因子統計值①Tab.2 Summary statistics of climate factors

2.2 模型構建

根據樹木生長規律，胸徑生長量(annual increment of DBH, AID)大于或等于0; 但是當AID等于0時，其對數變換lnAID無法計算。為避免這些誤差，Calama等(2005)和李春明(2012)在對數變換之前給胸徑生長增加常數1，以胸徑生長量加1的對數值ln(AID +1)為模型因變量。本研究構建的胸徑生長量模型形式如下：

ln(AID+1)=x0+x1lnDBH+x2DBH+

x3BA+x4Dg+x5Hd+x6lnA+x7N+

x8BAL+x9MCMT+x10AP+x11AHM+

x12SMMT+x13WMWT+x14SMT+x15MAT+

x16DD0+x17MWMT。

(1)

其中： AID為 1984—2010 年每2年間隔的年均胸徑生長量(為排除不同間隔期帶來的影響，僅研究每2年間隔的年均胸徑生長量)；x0～x17為參數。

數據的對數變換會導致胸徑生長量估計偏差，Baskerville(1972)給出了對這種偏差進行修正估計的近似值：

(2)

2.3 BMA原理

假設y為研究感興趣的量(胸徑生長量)，Z為調查所得數據(氣候因子和內部因子)，f={f1,…,fk}(k=1,…,n)代表所有可能模型組成的模型空間，而包含哪些變量的模型為最優模型事先并不知道，即模型本身存在不確定性。通常情況下，可以考慮的模型數量很大，如有p個自變量，那么在僅考慮主效應時其模型數量為2p。BMA選取可能模型的子集，并用這些模型的后驗概率進行所有推理和預測(Kassetal., 1995)。根據貝葉斯模型平均法，y的后驗分布為：

(3)

式中：P(y|fk,Z)為給定調查數據Z和模型fk條件下y的后驗分布；P(fk|Z)為給定調查數據Z條件下fk為最優模型的概率。

由方程(3)可知，y的后驗分布實際上是以后驗概率P(fk|Z)為權重對所有模型后驗分布進行加權的一個平均值。在方程(3)中，假定fk為最優模型，則概率分布為：

(4)

式中：β=(β0,β1,…,βn)為模型fk回歸系數的向量。

對于變量的先驗分布，本研究利用無信息先驗分布(Boxetal.， 1973)。

模型空間中所有模型數量很大，假設自變量為5個，那么即使不考慮變量間的交互作用，可構建的模型也可達32個，然而在實際研究中只有滿足要求的模型才能被選進來。Madigan等(1994)建議采用Occam窗法排除不當的模型以減少模型數量： 當一個模型后驗概率小于最優模型后驗概率的5%時，則從模型空間中移除該模型。本研究采用Occam窗法排除不當的模型，排序得到后驗概率大的模型。

根據貝葉斯理論，備擇假設H1為：Xk為杉木單木胸徑生長量變化驅動因子的后驗概率有多大，即βk不等于0的后驗概率有多大。BMA以所有包含變量Xk的杉木單木胸徑生長模型后驗概率的和作為βk不等于0的后驗概率的估計：

(5)

式中：Ii為0/1指示變量，當βk在模型fi中，Ii=1，反之Ii=0；Ω為模型空間。

Kass等(1995)認為,BMA變量后驗概率推斷的規則如下：P(βk≠0|Z)<0.5，表示沒有證據表明Xk是杉木單木胸徑生長量變化的驅動因子； 0.5 ≤P(βk≠0|Z)<0.75，表示有弱的證據表明Xk是杉木單木胸徑生長量變化的驅動因子； 0.75 ≤P(βk≠0|Z)<0.95，表示有強的證據表明Xk是杉木單木胸徑生長量變化的驅動因子；P(βk≠0|Z)≥0.95，表示有很強的證據表明Xk是杉木單木胸徑生長量變化的驅動因子。

若Xk是杉木單木胸徑生長量變化的驅動因子，則根據BMA模型似然函數推導得出杉木單木胸徑生長模型中變量因子Xk的點估計的后驗概率均值為：

(6)

之后對該變量后驗概率均值進行排序，選擇杉木單木胸徑生長量變化的主要驅動因子，并構建相應的不確定性胸徑生長模型。

2.4 模型評價

采用決定系數(R2)和平均絕對差(mean absolute difference, MAD)作為模型預測的評價指標。R2越大、MAD越小，模型精度越高：

(7)

(8)

為探究不同初植密度單木胸徑生長響應差異，本研究分不同密度以及所有樣本進行分析處理。在不同處理中，隨機抽取總數據的60%作為建模數據，剩余40%作為檢驗數據。同時，分別利用BMA和SR(雙向選擇)2種方法分析杉木單木胸徑生長量與內部和外部因子的關系。BMA分析通過R軟件BMA包實現(Rafteryetal.， 2005)，SR分析通過step函數實現。

3 結果與分析

3.1 BMA和SR確定的最佳模型比較

在E密度中，SR確定的模型與BMA確定的最優模型相似。在A和D密度中，SR確定的模型與BMA確定的第二優模型相似。在B、C密度和全樣本數據中，SR確定的模型不在BMA確定的后驗概率較大的前幾個模型中(表3)。B、C密度和全樣本數據3種處理在模型選擇上表現出明顯差異。

表3 BMA和SR確定的單木胸徑生長模型和BMA模型的后驗概率①Tab.3 Individual tree diameter growth model by BMA and SR methods and posterior probability of BMA model

除此之外，2種方法確定的模型R2和MAD在5種密度水平和所有密度水平的集合數據集中幾乎相同，模型預測效果沒有明顯差異(表4)。

表4 基于BMA和SR的單木胸徑生長模型評價Tab.4 Evaluation of individual tree diameter growth model based on BMA and SR methods

3.2 BMA和SR的參數估計值和P比較

A密度：BMA模型中，單木胸徑生長量隨DBH、MWMT、MCMT、AP和AHM增加而增加，隨N、Dg、lnDBH、lnA、SMMT和WMMT增加而減小；SR模型中，單木胸徑生長量隨BA、DBH、MWMT、MCMT、AP、AHM、DD0和SMT增加而增加，隨N、lnDBH、lnA、MAT和WMMT增加而減小(表5、圖1)。

B密度：BMA模型中，單木胸徑生長量與Hd、DBH、MWMT、MCMT和AP呈正相關關系，與N、Dg、lnDBH、lnA、BAL、SMMT和WMMT呈負相關關系；SR模型中，單木胸徑生長量與DD0和SMT呈正相關關系(表5、圖1)。

C密度：BMA模型中，單木胸徑生長量與BA、DBH、lnDBH、MWMT、MCMT、AP和DD0呈正相關關系，與N、Dg、lnA、BAL、SMMT和WMMT呈負相關關系；SR模型中，單木胸徑生長量與Hd、MAT呈負相關關系，與AHM和SMMT呈正相關關系(表5、圖1)。

D密度：BMA模型中，單木胸徑生長量與BA、Hd、DBH和MAT顯著正相關，與N、Dg、lnA、BAL、AHM、DD0、WMMT和SMT顯著負相關；SR模型中，單木胸徑生長量與AP和SMMT顯著負相關，與MWMT顯著正相關(表5、圖1)。

E密度：BMA模型中，單木胸徑生長量與BA、Dg、DBH、MWMT、MCMT、DD0和SMT呈正相關關系，與lnDBH、lnA、BAL、MAT、SMMT和WMMT呈負相關關系；SR模型中，單木胸徑生長量也與AHM呈負相關關系(表5、圖1)。

全樣本數據處理中，BMA選擇的變量與D密度相似，較D密度多選擇了MWMT、MCMT和AP，且顯著正相關，少選擇了AHM；逐步回歸模型中，單木胸徑生長量與lnDBH顯著負相關，與AHM顯著正相關(表5、圖1)。

在多數處理中，BMA和SR的年均胸徑生長量隨DBH、BA、MWMT、MCMT、AP增加而增加，隨N、lnA、BAL、WMMT增加而減小，且隨期初胸徑增大，年均胸徑生長量逐漸變小(圖2)。Dg在A、B、C、D密度和全樣本數據中與胸徑生長量表現為負相關關系，但是在E密度中表現為正相關關系。與DBH不同，lnDBH只在C密度中與胸徑生長量表現為正效應關系，在A、B、E密度中表現為負效應關系，在另外2種處理中與胸徑生長量的關系不顯著。MAT在中低密度(A、B、C密度)中對胸徑生長量的影響不大，在D密度中表現為正效應關系，在E密度和全樣本數據中表現為負效應關系。對于AHM、SMMT和SMT，BMA只在較少處理中認為與胸徑生長量顯著相關，但是SR在大多數處理中都選擇了這些變量(表5)。

BMA認為，Hd在B、D密度和全樣本數據中與胸徑生長量呈正相關關系，與SR略有不同，SR認為在C密度中胸徑生長量隨Hd增加而減小。在A、B、C密度和全樣本數據中，BMA認為胸徑生長量隨AP增加而增加；而在D密度中，SR認為AP與胸徑生長量呈負效應關系。

此外，SR選擇的變量比BMA更多。2種方法選擇的內部因子幾乎相同，BMA只在A密度中多選擇了Dg，SR在A密度中多選擇了BA，在全樣本數據中多選擇了lnDBH，但是在每種處理中，SR都比BMA多選擇了氣候因子。由此可見，SR更容易選擇冗余變量到模型中(表5)。

4 討論

4.1 BMA與SR的比較

本研究基于BMA與SR構建單木胸徑生長模型，SR根據自變量對因變量的意義選擇變量，但是容易忽略模型本身的不確定性； BMA的優點已在幾種不同類型模型中進行了評估，如線性回歸、生存模型等(Raftery， 1996； Volinskyetal.， 1997； Murphyetal.， 2001)，均表明BMA可提高模型性能。本研究結果同樣表明，BMA在變量選擇上優于SR。

采用SR擬合的模型在3種處理下位于BMA選擇的后驗概率高的前幾個模型中，僅在E密度中具有與BMA選擇最佳模型相同的預測變量(表3)，這表明模型本身存在不確定性，SR忽略了這一點。

表5 SR和BMA構建單木胸徑生長模型的參數估計值①Tab.5 The parameter estimates of individual tree diameter growth model by BMA and SR methods

圖1 不同密度處理下BMA模型各自變量后驗概率Fig.1 The posterior probability of variables in individual tree diameter growth model determined by BMA method

圖2 年均胸徑生長量與期初胸徑的關系Fig.2 The relationship between annual increment of DBH (AID)and diameter at the beginning of growth period

單木胸徑生長模型中，BMA認為，Hd在B、D密度和全樣本數據中與單木胸徑生長量呈正相關關系，而SR認為，在C密度中單木胸徑生長量隨Hd增加而減小。一般情況下，胸徑生長量應隨優勢高增加而增加。BMA認為，在A、B、C密度和全樣本數據中單木胸徑生長量隨年均降雨量增加而增加，而在D密度中，SR認為年均降雨量與胸徑生長量呈負效應關系。有研究發現，降雨量增加可以促進樹木生長(余黎等， 2014； 張海平， 2017)。可見，BMA的結論比SR具有更大的可解釋性。

此外，SR選擇的變量比BMA更多，即選擇冗余變量較多。單木胸徑生長模型中SR選擇了更多的氣候因子，余黎等(2014)、歐強新等(2019)發現，單木胸徑年均生長量受氣候因子影響較小，而競爭因子和單木大小因子是影響單木胸徑年均生長量的主要因子，這與BMA的結果是一致的。Genell等(2010)利用模擬數據表明，BMA不選擇冗余變量的概率高于SR，且選擇真實預測因子的概率與SR相似。

需要說明的是，2種方法在模型預測性能上的差異較小。Zapata-Cuartas等(2012)研究表明，對于小樣本，貝葉斯方法的性能優于經典方法。Wang等(2019)將分層貝葉斯方法與非線性混合效應模型進行比較，并估計地上生物量，結果發現2種方法在大樣本時表現相似，而對于小樣本貝葉斯方法要好得多。

4.2 單木胸徑生長量的影響因子

內部因子中，年均胸徑生長量隨單木期初胸徑增加而增加，隨lnA增加而減小，這表明單木期初胸徑越大，單木生長優勢越大，且在幼齡時增長迅速，隨著年齡增加增長速度變緩，即隨著期初胸徑增大，年均胸徑生長量變小(圖2)。一些研究表明，期初胸徑達到一定大小后，即達到大徑階后，其增長逐漸緩慢(張海平， 2017)。但lnDBH只在C密度中對年均胸徑生長量有正效應影響，在3種處理(A、B、E密度)中與年均胸徑生長量呈負相關關系。

年均胸徑生長量隨林分密度、林分平方平均胸徑和大于對象木的斷面積和增加而減小，與以往研究結果相同(Brooksetal.， 1998； Huangetal.， 2010)。大于對象木的斷面積之和越小，林分內單木平均擁擠程度越小； 林分密度和林分平方平均胸徑越小，林分內單木對生長空間的平均占有和利用越大，故單木年均胸徑生長量越大(余黎等， 2014)。這都表明競爭越小，胸徑生長量越大。但是在中高密度(C、D、E密度)和全樣本數據中，胸高斷面積對年均胸徑生長量有正效應影響。在B、E密度和全樣本數據中，年均胸徑生長量隨優勢木平均高增加而增加，在其他處理中優勢木平均高對年均胸徑生長量沒有影響。一般來說，樹木生長與優勢木平均高呈正相關關系，王冬至等(2017)研究華北落葉松-白樺混交林發現，立地質量對胸徑生長量也具有一定影響，并與胸徑生長量呈正相關關系，但影響較小。

氣候因子中，最熱月平均溫度、最冷月平均溫度和年均降雨量均對年均胸徑生長量有促進作用。有研究發現，年均降雨量對冷杉(Abiesnephrolepis)胸徑生長量表現為正效應關系，但達到一定閾值后表現為負效應關系(余黎等， 2014)，即當水分成為限制因子時，降雨往往對生長表現出正效應，而當水分充足或過多時，逐漸表現出負效應。也有研究發現，生長季水分多少可直接或間接影響細胞的分裂速率(程瑞梅等， 2015)，而且還在一定程度上影響生長季中樹木體內營養物質的積累(姜倩倩等， 2012)。張海平(2017)研究發現，生長季最低溫對年均胸徑生長量表現為正效應關系。李文馨等(2015)發現，年均積溫和年均降雨對胸徑生長量有促進作用，該結果與本研究結果相反。

另外，在低密度(A密度)時，干旱指數越大年均胸徑生長量越大，在高密度(D密度)時，干旱指數越大年均生長量越小。余黎等(2014)發現，年均胸徑生長量隨年均氣溫與年均降雨量之比增大而下降，與本研究在D密度的結果是一致的，在A密度呈正相關的原因可能是低密度時降雨對每株樹都是充足的，所以干旱指數對年均胸徑生長量沒有抑制作用。

5 結論

本研究采用逐步回歸法和貝葉斯模型平均法構建杉木人工林單木胸徑生長模型，逐步回歸法獲得模型的后驗概率小于BMA獲得最佳模型(最高后驗概率)或逐步回歸模型不在BMA模型空間前幾個后驗概率高的模型中，由此表現出杉木單木胸徑生長模型本身的不確定性。此外，逐步回歸法更容易選擇冗余變量。

內部因子和氣候因子均會影響杉木單木胸徑生長，但單木胸徑年均生長量受氣候因子的影響較小，而競爭因子和單木大小因子是影響單木胸徑年均生長量的主要因子。單木胸徑生長量隨林分密度、林分平方平均胸徑、大于對象木的斷面積和、年齡、冬季平均最低溫度增加而減小，隨期初胸徑、胸高斷面積、優勢木平均高、最冷月平均溫度、最熱月平均溫度、年均降雨量增加而增加。總的來說，杉木單木胸徑生長量隨競爭增加而減小，隨溫度和降雨增加而增加。