基于混合效應的湖北杉木冠長率模型

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

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

0 引言

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

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

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

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

1 材料與方法

1.1 研究區概況

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

1.2 研究方法

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

表1 樣地和解析木基本信息統計

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

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

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

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

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

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

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

式中，μij為基于固定效應參數的預測值，α、β為指數函數和冪函數的參數。

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

2 結果與分析

2.1 冠長率與林木及林分因子間相關性

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

表2 冠長率與林木、林分因子的Person相關性分析

圖1 冠長率等級分布

2.2 基礎模型建立

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

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

2.3 混合模型參數的確定

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

表3 基于樣地效應的混合模型擬合結果

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

2.4 方差-協方差矩陣的確定

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

表4 基于隨機參數效應不同方差協方差結構混合模型擬合結果

2.5 模型參數估計和獨立樣本檢驗

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

表5 基于不同方差函數混合效應模型擬合結果比較

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

表6 基礎模型與混合模型參數估計結果

圖2 基礎模型預估值殘差

圖3 混合模型預估值殘差

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

表7 基礎模型和混合模型檢驗

3 結論與討論

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

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

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