基于混合效應(yīng)的湖北杉木冠長率模型

袁 慧，杜超群，趙小濤，肖志明，趙 虎

（1湖北省林業(yè)科學(xué)研究院，武漢430075；2陽新縣自然資源和規(guī)劃局，湖北陽新435200；3陽新縣自然資源和規(guī)劃局龍港自然資源所，湖北陽新435200）

0 引言

樹冠是樹木自身進(jìn)行同化作用、異化作用、水循環(huán)等生理過程的重要場所，直接影響著樹木的生長變化[1]、地上部分器官的垂直分布、同化器官和非同化器官的空間分配，又影響樹木截獲光照的程度，進(jìn)而影響樹木的能量交換和有機(jī)物的積累[2]。樹木的活枝條對凈光合作用具有積極貢獻(xiàn)[3]。

冠長率(crown ratio，CR)是活樹冠長度與樹高的比值，往往用來預(yù)測樹木或森林的生長和收獲量，是一種無量綱的測度[4]，是研究樹冠結(jié)構(gòu)非常重要的指標(biāo)之一。冠長率可以反映樹木活力、林木質(zhì)量、林分密度和競爭的影響[5-6]，也是部分經(jīng)濟(jì)樹種一個很重要的特征量[7]，同時對樹木局部修枝也有指導(dǎo)作用[8]。有研究結(jié)果表明，冠長率影響樹木的生長。歐建德等[9]對南方紅豆杉進(jìn)行修枝處理，發(fā)現(xiàn)當(dāng)修枝后冠長率降低14%時，胸徑、樹高、材積和胸高形數(shù)均顯著增加。袁從軍等[10]研究黔中喀斯特次生林中多個樹種冠型特征，得出冠長可以解釋胸徑和樹高的0.64～0.88。

冠長率或樹冠高度往往應(yīng)用線性或非線性異速生長方程，把樹木及立地特征因子作為模型回歸變量進(jìn)行擬合[11]。Ward采用多水平線性模型構(gòu)建紅橡樹冠長率[12]，分析了不同林齡冠長率差異。Jari利用非線性模型研究芬蘭地區(qū)間伐對樟子松冠長率的影響[11]，結(jié)果顯示間伐強(qiáng)度和間伐后持續(xù)的時間對冠長率影響較大。與線性回歸模型不同，由于模型固有結(jié)構(gòu)，大多數(shù)非線性冠長率模型的預(yù)估值為0～1，能夠得到合理的預(yù)估結(jié)果[13]。非線性模型應(yīng)用較廣泛，其中Logistics方程被應(yīng)用于很多研究中。Hasenauer等[14]研究了奧地利多個樹種的冠長率模型，分析了5000多個固定樣地中42000多株樹木的數(shù)據(jù)，利用Logistics方程進(jìn)行非線性擬合，把樹木大小特征量、林分密度和立地因子組合作為模型變量，得出樹木大小因子和競爭因子對林分主要樹種冠長率的影響程度最大，達(dá)到17%～41%的解釋程度，而立地因子影響程度較小，不足10%。符利勇等[15]研究東北地區(qū)蒙古櫟冠長率，構(gòu)建了多水平非線性混合效應(yīng)模型，結(jié)果顯示，立地和樹木大小因子對模型的預(yù)估作用最大，其中平均胸徑、優(yōu)勢木樹高和優(yōu)勢木胸徑變量是重要的模型預(yù)測變量。

湖北具有大面積的杉木人工林，分布于鄂東南、鄂西和鄂西北等地，各地的環(huán)境因子差異較大，對杉木人工林的生長影響也較大。本研究通過在湖北省多個杉木栽植地區(qū)設(shè)立臨時標(biāo)準(zhǔn)樣地，利用解析木測量數(shù)據(jù)擬合杉木冠長率，分析單木因子、林分因子和環(huán)境因子對冠長率的影響，預(yù)測冠長率變化，旨在為杉木人工林合理經(jīng)營提供參考。

1 材料與方法

1.1 研究區(qū)概況

湖北省地處中國長江中游，東經(jīng)108°21′42″—116°07′50″、北緯 29°01′53″—33°6′47″。全省總面積18.59萬km2，其中東西跨度740 km、南北距離470 km。地勢三面高、中間低，向南敞開，北部有不完整盆地，整體處于中國地勢第二級階梯向第三級階梯過渡地帶，呈北西向南東走向。最高峰為神農(nóng)架神農(nóng)頂，海拔3105 m，平原區(qū)海拔約為20～100 m。年均氣溫15～17℃，無霜期230～300天，年均降水量800～1600 mm。除高山地區(qū)外，大部分為亞熱帶季風(fēng)性濕潤氣候。

1.2 研究方法

1.2.1 數(shù)據(jù)收集與整理 以2016年6月—2017年10月杉木人工林標(biāo)準(zhǔn)地調(diào)查數(shù)據(jù)為基礎(chǔ)，樣地主要分布在湖北黃石、咸寧、襄陽、宣恩等杉木集中栽植地區(qū)。選擇具有代表性、不同立地條件且生長良好的杉木林分設(shè)置臨時標(biāo)準(zhǔn)樣地，樣地面積600 m2(20 m×30 m)。記錄樣地海拔、坡度等立地信息，樣地內(nèi)每木檢尺，起測胸徑為5 cm，分別測量樹高、胸徑、冠幅和活枝下高等，使用超聲波測高器測量樹高及枝下高。每個樣地內(nèi)選擇3株樹干通直、長勢良好的平均木和1株優(yōu)勢木，解析木伐倒后使用皮尺測量樹高、冠幅及冠長，伐株按2 m區(qū)分段獲取解析木圓盤。共設(shè)立標(biāo)準(zhǔn)樣地181塊，伐取解析木715株，按照總樣本量的30%做檢驗樣本，其中建模數(shù)據(jù)500株、檢驗數(shù)據(jù)215株。參與建模與檢驗的數(shù)據(jù)基本統(tǒng)計信息見表1。

表1 樣地和解析木基本信息統(tǒng)計

1.2.2 冠長率模型的構(gòu)建 按照冠長率的定義，冠長率的范圍為0～1，0表示沒有樹冠或者死亡，1表示全樹冠。由于Logistics方程受到[0,1]范圍的限制，可以很好地解釋冠長率的變化趨勢，已有許多研究應(yīng)用Logistics方程來描述針葉和闊葉樹種的冠長率[15-20]，基礎(chǔ)方程如式(1)。

式中，CR為冠長率，β×X表示一個參數(shù)為β、獨立變量為X的函數(shù)，e為自然常數(shù)。變量X通常由代表林木大小、競爭變量和立地條件等變量因子構(gòu)成。

本研究冠長率相關(guān)變量中的林木大小因子主要包括樹高、胸徑和高徑比，樹木競爭因子變量包括每公頃株數(shù)、大于對象樹木斷面積，立地條件因子包括樣地海拔、坡度。

在基礎(chǔ)模型基礎(chǔ)之上，考慮樣地水平的隨機(jī)效應(yīng)對冠長率的影響，構(gòu)建湖北人工杉木冠長率非線性混合效應(yīng)模型。根據(jù)構(gòu)建混合效應(yīng)模型的一般步驟[21-22]，首先確定模型的固定效應(yīng)參數(shù)與隨機(jī)效應(yīng)參數(shù)，再確定用于解釋不同樣地內(nèi)方差-協(xié)方差矩陣結(jié)構(gòu)，最后確定方差和協(xié)方差的隨機(jī)效應(yīng)矩陣。非線性混合效應(yīng)模型形式如式(2)。

式中，i、j分別代表第一水平和觀測值，本研究僅考慮樣地水平；CRij為第i個樣地中第j株樹冠長率的觀測值；M為樣地數(shù)量；ni為在第i個樣地上的株數(shù)；f為含有參數(shù)向量φij和協(xié)變量向量υij的函數(shù)；β為p×1維固定效應(yīng)參數(shù)向量，bi為帶有方差-協(xié)防差矩陣D的q×1維隨機(jī)效應(yīng)向量；Aij和Bij為設(shè)計矩陣；eij為模型服從正態(tài)分布的誤差項；σ2為方差；Ri為分組樣地內(nèi)的方差-協(xié)方差矩陣[24]；D為隨機(jī)效應(yīng)的方差-協(xié)方差矩陣；Gi為描述方差異質(zhì)性的對角矩陣；Γ為樣地內(nèi)誤差的相關(guān)性結(jié)構(gòu)。

將所有可能的參數(shù)組合均作為隨機(jī)效應(yīng)進(jìn)行擬合，利用赤池信息準(zhǔn)則(akaike information criteria，AIC)、貝葉斯信息準(zhǔn)則(bayesian information criteria，BIC)和對數(shù)似然值(log likelihood，LL)統(tǒng)計指標(biāo)對模型進(jìn)行評價[15,22,24]。AIC、BIC值越小，對數(shù)似然值越大，模型擬合效果越好。分別選取含有不同個數(shù)參數(shù)效應(yīng)的最優(yōu)模型進(jìn)行似然比檢驗(likelihood ratio test，LRT)，最終選擇參數(shù)較少且模型顯著的作為最優(yōu)模型。

利用廣義正定矩陣和對角矩陣檢驗隨機(jī)效應(yīng)的方差-協(xié)方差矩陣，選擇AIC、BIC小且對數(shù)似然值較大的矩陣形式。通過消除樣地內(nèi)測量值間的誤差相關(guān)性和異方差問題確定組內(nèi)方差-協(xié)方差結(jié)構(gòu)，若觀測數(shù)據(jù)是連續(xù)型的，則要考慮誤差相關(guān)性。本研究中冠長率數(shù)據(jù)是一次調(diào)查，只需考慮異方差問題，采用指數(shù)函數(shù)及冪函數(shù)對模型進(jìn)行矯正，并用AIC、BIC、對數(shù)似然值比較矯正效果，指數(shù)函數(shù)和冪函數(shù)形式如式(3)～(4)。

式中，μij為基于固定效應(yīng)參數(shù)的預(yù)測值，α、β為指數(shù)函數(shù)和冪函數(shù)的參數(shù)。

模型擬合及檢驗除使用AIC、BIC、log-likelihood統(tǒng)計指標(biāo)，還采用均方根誤差(RMSE)、平均誤差(Bias)和擬合精度(R2)指標(biāo)。本研究中基礎(chǔ)模型及非線性混合模型均采用R軟件(3.6.1)中nlme程序包實現(xiàn)。

2 結(jié)果與分析

2.1 冠長率與林木及林分因子間相關(guān)性

冠長率是樹冠長度與樹木高度間比值。按照冠長率范圍0～1，把(0,0.1]的冠長率劃分為等級1，以此類推共劃分為10個等級。本研究區(qū)杉木人工林平均冠長率為0.49，范圍在0.1～0.9之間。由圖1可以看出，冠長率總體呈正態(tài)分布，等級主要集中分布在4～7，占總體的85%。冠長率與林木及林分因子的相關(guān)性如表2所示，胸徑和樹高與冠長率呈顯著正相關(guān)，而每公頃株數(shù)、大于對象木斷面積、高徑比與冠長率呈顯著負(fù)相關(guān)，海拔、坡度和年齡因子對冠長率無顯著影響，說明林木大小及競爭因子對冠長率影響較大，而立地因子影響較小。樹木間競爭程度大抑制樹冠的生長，致冠長率降低。

表2 冠長率與林木、林分因子的Person相關(guān)性分析

圖1 冠長率等級分布

2.2 基礎(chǔ)模型建立

根據(jù)冠長率與樹木及林分因子的相關(guān)性分析，將相關(guān)性顯著的DBH（胸徑）、HT（樹高）、SD（每公頃株數(shù)）、BAL（大于對象樹木斷面積）、HDR（高徑比）因子引入到Logistics方程作為變量。因此，杉木冠長率模型的具體形式如式(5)。

式中，a、b、c、d、e、f均為模型參數(shù)。使用R軟件nlme程序包中nls函數(shù)，擬合得到a=1.055，b=-0.693，c=-0.056，d=0.001，e=-0.01，f=-0.121，模型決定系數(shù)R2為0.3231。

2.3 混合模型參數(shù)的確定

將模型所有參數(shù)及其組合形式（49種）作為隨機(jī)效應(yīng)進(jìn)行擬合，由于模型參數(shù)較多，部分組合形式擬合無法收斂。對含有不同隨機(jī)參數(shù)個數(shù)的隨機(jī)效應(yīng)模型分別選擇1個最優(yōu)的模型進(jìn)行比較，結(jié)果如表3所示。

表3 基于樣地效應(yīng)的混合模型擬合結(jié)果

當(dāng)模型沒有隨機(jī)參數(shù)時，AIC=-975.154，BIC=-867.5326。當(dāng)隨機(jī)參數(shù)增加1個最優(yōu)參數(shù)a時，AIC(-1005.499)和BIC(-961.0544)的值減小。當(dāng)增加2個隨機(jī)參數(shù)a和f時，AIC(-1006.777)和BIC(-968.920)的值繼續(xù)減小。當(dāng)考慮其他隨機(jī)參數(shù)組合時，增加3個隨機(jī)參數(shù)b、c、e和4個隨機(jī)參數(shù)a、b、c、e的AIC和BIC的值相較無參數(shù)均降低，但與1個參數(shù)a和2個參數(shù)a、f的結(jié)果比較卻增大。似然比檢驗表明，無隨機(jī)參數(shù)和1個參數(shù)a模型有顯著不同(LRT=10.32，P=0.0221)，參數(shù)a模型擬合和參數(shù)a和f模型擬合有顯著不同(LRT=9.46，P=0.0314)。說明通過增加隨機(jī)參數(shù)可以顯著提高模型擬合精度。綜上分析，選擇含有隨機(jī)參數(shù)a和f為最優(yōu)混合模型。

2.4 方差-協(xié)方差矩陣的確定

確定混合模型樣地間方差-協(xié)方差矩陣，對比廣義正定矩陣和對角矩陣對混合模型擬合效果的影響。結(jié)果（表4）顯示，廣義正定矩陣的擬合效果較好，AIC(-1006.778)和BIC(-961.053)值小于對角矩陣，似然比檢驗表明兩者差異顯著(LRT=5.28，P=0.0216)。因此選擇廣義正定矩陣作為隨機(jī)效應(yīng)的方差-協(xié)方差矩陣形式。

表4 基于隨機(jī)參數(shù)效應(yīng)不同方差協(xié)方差結(jié)構(gòu)混合模型擬合結(jié)果

2.5 模型參數(shù)估計和獨立樣本檢驗

利用冪函數(shù)和指數(shù)函數(shù)來描述混合模型產(chǎn)生的異方差現(xiàn)象，擬合效果見表5。指數(shù)函數(shù)的AIC(-1015.075)和BIC(-964.7797)值小于冪函數(shù)，說明指數(shù)函數(shù)擬合效果較好，似然比檢驗也表明指數(shù)函數(shù)明顯提高了模型擬合精度。

表5 基于不同方差函數(shù)混合效應(yīng)模型擬合結(jié)果比較

混合效應(yīng)模型采用指數(shù)函數(shù)擬合后，由模型的殘差結(jié)果可知，殘差分布范圍明顯減?。▓D2～3）?；A(chǔ)模型與混合模型的參數(shù)估計結(jié)果見表6?；A(chǔ)模型的Bias、RMSE、R2分別為-0.0862、0.1458、0.9190，混合模型Bias、RMSE、R2分別為-0.0095、0.1175、0.9475，可知混合模型擬合效果較好。

表6 基礎(chǔ)模型與混合模型參數(shù)估計結(jié)果

圖2 基礎(chǔ)模型預(yù)估值殘差

圖3 混合模型預(yù)估值殘差

混合模型應(yīng)用獨立樣本進(jìn)行檢驗，一般包括固定效應(yīng)部分檢驗和隨機(jī)效應(yīng)部分檢驗。隨機(jī)效應(yīng)部分的檢驗需要二次抽樣來計算隨機(jī)參數(shù)值，固定效應(yīng)部分檢驗相當(dāng)于傳統(tǒng)回歸分析檢驗，即采用固定參數(shù)估計值，本研究只涉及固定效應(yīng)部分檢驗。利用215株樹木樣本對基礎(chǔ)模型及混合模型進(jìn)行檢驗（表7），混合模型的檢驗精度大于基礎(chǔ)模型，且均方根誤差和平均誤差均小于基礎(chǔ)模型，因此混合模型擬合效果顯著優(yōu)于基礎(chǔ)模型。

表7 基礎(chǔ)模型和混合模型檢驗

3 結(jié)論與討論

為保證樹冠長度數(shù)據(jù)準(zhǔn)確性，本研究所用樹冠長度數(shù)據(jù)由樣地調(diào)查中伐倒解析木實地測量得到。通過分析冠長率與林木及林分因子相關(guān)性，將相關(guān)性顯著因子引入Logistics方程作為變量，并建立非線性混合效應(yīng)模型。杉木冠長率主要集中在0.3～0.6范圍內(nèi)，占總樣本量的85%。對混合模型所有參數(shù)及其組合形式均作為隨機(jī)效應(yīng)進(jìn)行擬合，分析模型AIC、BIC值和似然比顯著性檢驗，與無隨機(jī)參數(shù)相比較，含有隨機(jī)參數(shù)a、f的混合模型擬合效果最優(yōu)。采用廣義正定矩陣來確定樣地間方差-協(xié)方差矩陣形式，利用指數(shù)函數(shù)來描述混合模型產(chǎn)生的異方差，最終得到混合模型參數(shù)?；旌闲?yīng)模型應(yīng)用獨立樣本進(jìn)行檢驗包括固定效應(yīng)和隨機(jī)效應(yīng)檢驗，本研究采用固定效應(yīng)參數(shù)估計值，即相當(dāng)于傳統(tǒng)回歸分析檢驗，沒有計算隨機(jī)參數(shù)值和樣地隨機(jī)抽樣檢驗。根據(jù)獨立樣本檢驗結(jié)果，混合模型相比基礎(chǔ)模型的決定系數(shù)明顯提高，均方根誤差和平均誤差均有所降低，殘差范圍也顯著減小，綜合表明非線性混合效應(yīng)模型擬合效果優(yōu)良。

從混合模型擬合參數(shù)可知，胸徑變量參數(shù)值為正，每公頃株數(shù)、大于對象樹木斷面積和高徑比的參數(shù)值為負(fù)，除樹高變量參數(shù)外，其余變量參數(shù)結(jié)果與冠長率相關(guān)性分析結(jié)果吻合。結(jié)果顯示，樹高的增加導(dǎo)致冠長率的降低，這與Hailemariam等[16]研究加拿大哥倫比亞地區(qū)的闊葉樹種冠長率結(jié)果一致；而郭孝玉等[20]研究長白落葉松冠長率的結(jié)果顯示，樹高與冠長無顯著相關(guān)性，這說明對于不同地區(qū)、不同樹種，同一因子對冠長率的影響結(jié)果也有差別。胸徑是一個重要的樹木生長變化的變量，通常用于解釋林分結(jié)構(gòu)、樹木活力和競爭能力[15]，本研究中胸徑變量與冠長率顯著正相關(guān)，混合模型參數(shù)也為正值，與盧軍等[19]研究次生林部分樹種冠長率結(jié)果一致；郭孝玉等[20]研究結(jié)果也得到類似的擬合結(jié)果，表明胸徑指標(biāo)與冠長率關(guān)系密切，是一個很重要的變量。競爭因子變量結(jié)果意味著林分密度增加、大于對象樹木的林木競爭能力的增強(qiáng)，林木冠長率變小，說明樹木周圍環(huán)境競爭壓力的增加會抑制樹木本身樹冠的生長，進(jìn)而影響樹干的發(fā)育。擬合的變量參數(shù)對模型具有一定解釋意義。

立地因子的差異對不同樹種冠長率影響也有所不同[15-16,19]。本研究中只涉及湖北省人工杉木，沒有做研究區(qū)其他樹種的對比分析，研究中立地因子（海拔、坡度）對冠長率無顯著影響，可能與本研究中立地因子的類型和數(shù)據(jù)量較少有關(guān)，還無法有效反映立地的影響程度。