基于貝葉斯模型平均法和逐步回歸法構(gòu)建杉木單木胸徑生長(zhǎng)模型*

魯樂樂 王 震 張雄清 張建國(guó)

(1.中國(guó)林業(yè)科學(xué)研究院林業(yè)研究所 國(guó)家林業(yè)和草原局林木培育重點(diǎn)實(shí)驗(yàn)室 北京100091；2.南京林業(yè)大學(xué)南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心 南京210037)

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

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

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

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

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

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

1 研究區(qū)概況與數(shù)據(jù)調(diào)查

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

研究區(qū)杉木密度試驗(yàn)林為1982年使用1年生苗木造林，采用隨機(jī)區(qū)組試驗(yàn)設(shè)計(jì)，設(shè)置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)。每種密度重復(fù)3次，共設(shè)立樣地15塊，每塊樣地大小為20 m×30 m。造林后1984—1990年逐年調(diào)查，之后到2010年隔年調(diào)查，每次調(diào)查均在秋冬季進(jìn)行。測(cè)量樹高時(shí)，每塊樣地隨機(jī)抽取50株樹，計(jì)算最高5株樹的平均高作為林分優(yōu)勢(shì)木平均高。

2 研究方法

2.1 變量選擇

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

外部因子主要指氣候因子，通過ClimateAP軟件獲得(Wangetal.， 2012)，包括年均氣溫(mean annual temperature，MAT，℃)、年均降雨量(mean annual precipitation，AP，mm)、最冷月平均溫度(mean coldest month temperature，MCMT，℃)、年均干旱指數(shù)(annual heat-moisture index，AHM)、低于0 ℃天數(shù)(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，℃)。氣候因子統(tǒng)計(jì)值見表2。

表2 氣候因子統(tǒng)計(jì)值①Tab.2 Summary statistics of climate factors

2.2 模型構(gòu)建

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

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年間隔的年均胸徑生長(zhǎng)量(為排除不同間隔期帶來的影響，僅研究每2年間隔的年均胸徑生長(zhǎng)量)；x0～x17為參數(shù)。

數(shù)據(jù)的對(duì)數(shù)變換會(huì)導(dǎo)致胸徑生長(zhǎng)量估計(jì)偏差，Baskerville(1972)給出了對(duì)這種偏差進(jìn)行修正估計(jì)的近似值：

(2)

2.3 BMA原理

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

(3)

式中：P(y|fk,Z)為給定調(diào)查數(shù)據(jù)Z和模型fk條件下y的后驗(yàn)分布；P(fk|Z)為給定調(diào)查數(shù)據(jù)Z條件下fk為最優(yōu)模型的概率。

由方程(3)可知，y的后驗(yàn)分布實(shí)際上是以后驗(yàn)概率P(fk|Z)為權(quán)重對(duì)所有模型后驗(yàn)分布進(jìn)行加權(quán)的一個(gè)平均值。在方程(3)中，假定fk為最優(yōu)模型，則概率分布為：

(4)

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

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

模型空間中所有模型數(shù)量很大，假設(shè)自變量為5個(gè)，那么即使不考慮變量間的交互作用，可構(gòu)建的模型也可達(dá)32個(gè)，然而在實(shí)際研究中只有滿足要求的模型才能被選進(jìn)來。Madigan等(1994)建議采用Occam窗法排除不當(dāng)?shù)哪Ｐ鸵詼p少模型數(shù)量： 當(dāng)一個(gè)模型后驗(yàn)概率小于最優(yōu)模型后驗(yàn)概率的5%時(shí)，則從模型空間中移除該模型。本研究采用Occam窗法排除不當(dāng)?shù)哪Ｐ?，排序得到后?yàn)概率大的模型。

根據(jù)貝葉斯理論，備擇假設(shè)H1為：Xk為杉木單木胸徑生長(zhǎng)量變化驅(qū)動(dòng)因子的后驗(yàn)概率有多大，即βk不等于0的后驗(yàn)概率有多大。BMA以所有包含變量Xk的杉木單木胸徑生長(zhǎng)模型后驗(yàn)概率的和作為βk不等于0的后驗(yàn)概率的估計(jì)：

(5)

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

Kass等(1995)認(rèn)為,BMA變量后驗(yàn)概率推斷的規(guī)則如下：P(βk≠0|Z)<0.5，表示沒有證據(jù)表明Xk是杉木單木胸徑生長(zhǎng)量變化的驅(qū)動(dòng)因子； 0.5 ≤P(βk≠0|Z)<0.75，表示有弱的證據(jù)表明Xk是杉木單木胸徑生長(zhǎng)量變化的驅(qū)動(dòng)因子； 0.75 ≤P(βk≠0|Z)<0.95，表示有強(qiáng)的證據(jù)表明Xk是杉木單木胸徑生長(zhǎng)量變化的驅(qū)動(dòng)因子；P(βk≠0|Z)≥0.95，表示有很強(qiáng)的證據(jù)表明Xk是杉木單木胸徑生長(zhǎng)量變化的驅(qū)動(dòng)因子。

若Xk是杉木單木胸徑生長(zhǎng)量變化的驅(qū)動(dòng)因子，則根據(jù)BMA模型似然函數(shù)推導(dǎo)得出杉木單木胸徑生長(zhǎng)模型中變量因子Xk的點(diǎn)估計(jì)的后驗(yàn)概率均值為：

(6)

之后對(duì)該變量后驗(yàn)概率均值進(jìn)行排序，選擇杉木單木胸徑生長(zhǎng)量變化的主要驅(qū)動(dòng)因子，并構(gòu)建相應(yīng)的不確定性胸徑生長(zhǎng)模型。

2.4 模型評(píng)價(jià)

采用決定系數(shù)(R2)和平均絕對(duì)差(mean absolute difference, MAD)作為模型預(yù)測(cè)的評(píng)價(jià)指標(biāo)。R2越大、MAD越小，模型精度越高：

(7)

(8)

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

3 結(jié)果與分析

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

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

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

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

表4 基于BMA和SR的單木胸徑生長(zhǎng)模型評(píng)價(jià)Tab.4 Evaluation of individual tree diameter growth model based on BMA and SR methods

3.2 BMA和SR的參數(shù)估計(jì)值和P比較

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

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

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

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

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

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

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

BMA認(rèn)為，Hd在B、D密度和全樣本數(shù)據(jù)中與胸徑生長(zhǎng)量呈正相關(guān)關(guān)系，與SR略有不同，SR認(rèn)為在C密度中胸徑生長(zhǎng)量隨Hd增加而減小。在A、B、C密度和全樣本數(shù)據(jù)中，BMA認(rèn)為胸徑生長(zhǎng)量隨AP增加而增加；而在D密度中，SR認(rèn)為AP與胸徑生長(zhǎng)量呈負(fù)效應(yīng)關(guān)系。

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

4 討論

4.1 BMA與SR的比較

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

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

表5 SR和BMA構(gòu)建單木胸徑生長(zhǎng)模型的參數(shù)估計(jì)值①Tab.5 The parameter estimates of individual tree diameter growth model by BMA and SR methods

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

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

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

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

需要說明的是，2種方法在模型預(yù)測(cè)性能上的差異較小。Zapata-Cuartas等(2012)研究表明，對(duì)于小樣本，貝葉斯方法的性能優(yōu)于經(jīng)典方法。Wang等(2019)將分層貝葉斯方法與非線性混合效應(yīng)模型進(jìn)行比較，并估計(jì)地上生物量，結(jié)果發(fā)現(xiàn)2種方法在大樣本時(shí)表現(xiàn)相似，而對(duì)于小樣本貝葉斯方法要好得多。

4.2 單木胸徑生長(zhǎng)量的影響因子

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

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

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

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

5 結(jié)論

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

內(nèi)部因子和氣候因子均會(huì)影響杉木單木胸徑生長(zhǎng)，但單木胸徑年均生長(zhǎng)量受氣候因子的影響較小，而競(jìng)爭(zhēng)因子和單木大小因子是影響單木胸徑年均生長(zhǎng)量的主要因子。單木胸徑生長(zhǎng)量隨林分密度、林分平方平均胸徑、大于對(duì)象木的斷面積和、年齡、冬季平均最低溫度增加而減小，隨期初胸徑、胸高斷面積、優(yōu)勢(shì)木平均高、最冷月平均溫度、最熱月平均溫度、年均降雨量增加而增加?？偟膩碚f，杉木單木胸徑生長(zhǎng)量隨競(jìng)爭(zhēng)增加而減小，隨溫度和降雨增加而增加。