馬愛云 李鳳日 董利虎

(東北林業(yè)大學，哈爾濱，150040)

樹冠是樹木的重要組成部分，它決定著樹木的生活力和生產(chǎn)力。在樹木生長過程中，樹冠反映了樹木長期競爭水平，是制定營林決策的重要樹木變量[1-3]。冠幅是樹冠研究中最常用的重要指標，可用來估計樹冠的表面積與體積[4]、模擬樹冠輪廓[5]、計算林木的競爭指數(shù)[6]等。此外，冠幅在林木生長與收獲模型中常作為預測變量出現(xiàn)。因此，建立精準的冠幅預測模型具有重要意義。

目前常采用傳統(tǒng)最小二乘法、啞變量方法以及線性和非線性混合效應模型方法等來構建冠幅預測模型[7-16]。Krajicek et al.[10]認為冠幅與胸徑呈線性關系，在一定范圍內，冠幅隨胸徑的增長而增大[10]；符利勇等[11]研究表明，基于Logistic、Richards非線性且具有生物學意義的理論模型能夠更好地描述冠幅與胸徑關系[12]；雷相東等[7]以落葉松云冷杉林為研究對象，采用多元逐步回歸方法建立了9個樹種的單木冠幅預測模型；張樹森等[8]通過將胸徑劃分為3個等級，并將這3個作為啞變量的方式建立疏開木冠幅模型。針對不同地區(qū)不同樹種，冠幅模型已經(jīng)從普通最小二乘法建模發(fā)展到采用線性或非線性混合效應方法進行建模，且冠幅預測模型的變量也從常用的單一變量胸徑逐漸加入各種林木變量和林分變量。

雜種落葉松(Larixkaempferi(Lamb.) Carr×LarixolgensisHenry)是落葉松屬中一個非常重要的速生樹種，具有高產(chǎn)、優(yōu)質、抗病蟲害等特性，是我國東北地區(qū)的重要用材樹種之一。當前，黑龍江省營造了大量的雜種落葉松示范林[13-17]。在雜種落葉松人工林后續(xù)的森林生產(chǎn)經(jīng)營管理理論研究中，肖銳等[18]針對不同初植密度雜種落葉松幼齡林的林分動態(tài)進行了模擬，王濤等[19]研究了雜種落葉松人工林幼齡林的枯損，但針對雜種落葉松的冠幅預測模型研究還未見報道。為了探索雜種落葉松冠幅預測模型，本文以黑龍江省江山嬌實驗林場1998年營造的雜種落葉松林為研究對象，具體為日本與長白落葉松的雜交種，從7種常用的冠幅與胸徑關系基礎模型中挑選較好的模型來進一步研究雜種落葉松冠幅預測模型，分析了雜種落葉松冠幅與胸徑及其他林木因子、林分因子的關系。最終將樣地作為隨機效應因子，構建雜種落葉松單木冠幅混合效應預測模型，為雜種落葉松人工林的可持續(xù)經(jīng)營及動態(tài)模擬提供依據(jù)。

1 研究區(qū)概況

黑龍江省江山嬌實驗林場，地處黑龍江省東南部，在黑龍江省寧安市境內，地理坐標為東經(jīng)128°53′16″～129°12′43″，北緯43°44′54″～43°54′12″。林場東南部與吉林省相接，西部以鏡泊湖為界，整個施業(yè)區(qū)處于東京城林業(yè)局之中，位于張廣才嶺南端，東西走向，東高西低，北高南低，屬于低山丘陵地區(qū)，海拔356～890 m。氣候屬亞寒帶大陸性氣候，年平均氣溫3.5 ℃，年平均降水500 mm。區(qū)域內主要林分類型為人工針葉林、闊葉混交林、針闊混交林和柞樹林4類，主要樹種有紅松(Pinuskoraiensis)、云杉(Piceaasperata)、落葉松(Larixgmelini)、水曲柳(Fraxinusmandshurica)等10余種，森林覆蓋率達95%。

2003年黑龍江省林業(yè)科學研究院以江山嬌實驗林場1998年種植的不同初植密度(2 500、3 300、4 400、6 600株/hm2)的雜種落葉松林為試驗地，在各初值密度的雜種落葉松人工林設置12塊固定樣地，共計48塊固定標準地，樣地面積為0.03～0.06 hm2，各固定樣地集中，立地條件基本相同，除一次輕微的人為修枝外，林分未進行過任何撫育間伐。2015年對48個樣地進行每木檢尺，用胸徑尺測量距地面1.3 m處的樹干直徑，即胸徑(DBH)，精確到0.1 cm，樹高(H)和第一活枝高(HCB)用超聲波測高器測量，精確到0.1 m，以樣木樹干為中心，用鋼尺分別對其東、西、南、北4個方向的樹冠半徑進行測量，精確到0.1 m。根據(jù)東、西、南、北4個方向的樹冠半徑，計算樹冠直徑平均值即樹木平均冠幅(CW)。除此之外，還計算了冠長率(CR)、高徑比(HDR)、林分斷面積(BAS)、優(yōu)勢木平均高(H0)、優(yōu)勢木平均胸徑(Dg)、大于對象木的胸高斷面積之和(BAL)等因子。48塊雜種落葉松林的林分因子統(tǒng)計信息見表1。

表1 雜種落葉松林基本因子統(tǒng)計

2 研究方法

2.1 基礎模型構建

2.2 線性混合效應模型

在線性模型中，所有的參數(shù)都是固定效應參數(shù)，線性混合效應模型是線性模型的進一步擴展，既包含固定效應又包含隨機效應，同時誤差有更為靈活的結構，還可以改變隨機效應的方差—協(xié)方差結構來反映個體間的差異[11-12，24，28，31-32]。本研究數(shù)據(jù)來自于不同初植密度的48塊樣地，為考慮樣地對冠幅的隨機干擾，以樣地為隨機效應因子，構建樣地單水平線性混合模型，模型形式如下：

(1)

表2 冠幅基礎模型

圖1 冠幅—胸徑相關關系

混合效應模型的構建，在R軟件nlme包的lme模塊中實現(xiàn)，模型的參數(shù)估計采用的是限制極大似然法(REML)。在未知隨機效應矩陣的情況下，考慮所有參數(shù)都可能存在隨機效應，對所有不同隨機效應參數(shù)組合的模型都進行擬合，在收斂的模型中，通過比較赤池信息準則、貝葉斯信息準則和負2倍對數(shù)似然值來確定最優(yōu)隨機效應組合。對于隨機效應的方差—協(xié)方差矩陣，本文采用常見的廣義正定矩陣、復合對稱和對角矩陣3種方差—協(xié)方差結構分別擬合，選擇赤池信息準則、貝葉斯信息準則和負2倍的對數(shù)似然值最小的矩陣結構作為隨機效應方差—協(xié)方差矩陣結構[21]。

2.3 模型的評價與檢驗指標

本研究采用修正決定系數(shù)、均方根誤差、赤池信息準則、貝葉斯信息準則、負2倍的對數(shù)似然值等指標對模型的擬合優(yōu)度進行評價比較，采用平均絕對偏差、平均相對偏差絕對值兩個偏差統(tǒng)計量和模型擬合效率一個檢驗指標對模型進行檢驗[25-27]。赤池信息準則、貝葉斯信息準則和負2倍的對數(shù)似然值由R軟件直接給出，其余指標計算公式如下。

(2)

均方根誤差(RMSE)：

(3)

平均絕對偏差(MAE)：

(4)

平均相對偏差絕對值(MAPE)：

(5)

模擬效率(EF)：

(6)

本文利用全部樣本來構建冠幅混合效應模型，采用留一交叉驗證法[24]進行模型檢驗，即每次取出一個樣地數(shù)據(jù)作為檢驗樣本，其余47個樣地數(shù)據(jù)作為建模數(shù)據(jù)，如此共循環(huán)檢驗48次，計算48次檢驗結果的平均值來衡量模型的預測性能。對混合效應模型隨機參數(shù)的檢驗采用經(jīng)驗線性無偏最優(yōu)預測法(EBLUP)來計算隨機效應參數(shù)值[24，27-28]。具體計算公式如下：

(7)

2.4 不同樣本量對隨機效應的校正

混合效應模型在應用時，對隨機效應參數(shù)的估計是最為關鍵的，本研究在檢驗模型時，用了各樣地所有樣木冠幅實測值來計算隨機效應參數(shù)，但在實際的森林資源調查中，要調查樣地內所有的樣木冠幅，是不切實際的，因此需要分析小量的樣本對隨機效應的校正精度的影響，確定較為合適的調查樣本量。根據(jù)前人的研究表明，不同數(shù)量的樣木對混合效應模型的預測有影響，預測偏差一般隨樣本量的增多而減小[29]，通常每個樣地抽取樣本量為4時，足以校正隨機效應[30]，Xie et al.[24]的研究表明抽取接近樣地優(yōu)勢木數(shù)量的6個樣本量時較為合理。因此本研究中，基于上述的留一法，采用隨機抽樣，在沒有替換的情況下，對每個樣地分別抽取2、3、4、…、30棵樹，分別重復抽取100次，組成不同抽樣株數(shù)的子樣本，用于計算各樣本量下的隨機效應參數(shù)，并得到混合模型對應的預測偏差，從而計算不同樣本量下的檢驗統(tǒng)計量，進行對比分析，最終得到較合適的隨機效應校正樣本量。

以上所有的計算過程均在R軟件4.0中實現(xiàn)，主要使用nlme包[25]。

3 結果與分析

3.1 模型的構建

3.1.1 單木因子與林分因子

基于選定的最優(yōu)基礎模型M4，采用再參數(shù)化的方法引入一些常用的單木因子和林分因子。具體做法如下：1)分樣地擬合基礎模型，得到各樣地的模型參數(shù)；2)分析參數(shù)與單木因子和林分因子的相關性，確定β=f(x)的具體形式，其中β為模型中的參數(shù)，x為相關性強的單木因子和林分因子；3)將f(x)的具體形式代回到模型中。除此之外還需要求各變量在統(tǒng)計上顯著(p<0.05),且自變量之間存在較弱的共線性(方差膨脹系數(shù)小于10)。再參數(shù)化構建的冠幅最優(yōu)線性模型形式為：

CW=a0+a1DBH+a2CR+a3HDR+a4BAS。

(8)

3.1.2 混合效應模型

CW=(a0+b0)+(a1+b1)DBH+a2CR+(a3+b3)HDR+

a4BAS。

(9)

表3 模型(8)在不同隨機效應參數(shù)下的擬合比較

對模型(9)采用不同隨機效應方差—協(xié)方差矩陣結構進行模擬，選擇赤池信息準則、貝葉斯信息準則和負2倍對數(shù)似然值最小時的廣義正定矩陣結構為協(xié)方差矩陣結構。

3.2 模型的評價與檢驗

圖2 冠幅觀測值與模型(8)、模型(9)的預測值相關圖

表4 模型(8)與模型(9)的參數(shù)和隨機效應方差估計

3.3 不同樣本量對隨機效應的校正對比

用不同抽樣數(shù)量組成的子樣本計算混合效應模型(9)的隨機效應參數(shù)時，模型的預測效果見圖4。其MAE和MAPE值整體上隨著抽樣數(shù)量的增大而減小，與模型擬合效率隨著抽樣株數(shù)量增大而增大的結果相符，說明增加計算隨機效應參數(shù)的樣本量，模型的預測精度會提高，這與前人的研究一致[24，29-30]。當樣本量為5株及大于5株時，模型預測精度提高減緩，對不同樣本量下的模型預測的誤差進行顯著性檢驗也表明，5株及更多的抽樣株數(shù)之間沒有顯著差異，當樣本量達到12株左右時，模型預估精度趨于穩(wěn)定不變。

圖3 不同徑階下模型(8)與模型(9)的預測誤差分布圖

表5 線性模型與混合效應模型檢驗結果

4 結論與討論

圖4 不同抽樣數(shù)量計算隨機效應的預測精度變化

1)在變量DBH、CR和HDR上添加樣地隨機效應時，模型預測效果最佳，線性混合效應模型的擬合與檢驗結果均比最優(yōu)線性模型好。本文比較了不同抽樣株數(shù)對隨機效應的校正，表明模型的預測效果隨著樣本數(shù)的增加而提高，在樣本量為12及以上時，模型的預測指標趨于穩(wěn)定不變，但當抽取株數(shù)為5及以上時，便沒有顯著性差異。因此，在應用該模型進行預測時，根據(jù)抽樣調查成本和精度的要求，建議隨機抽取接近樣地優(yōu)勢木株數(shù)的5個樣本即可，這與Xie et al.[24]的研究一致。

2)模型中，變量DBH和CR的參數(shù)均為正數(shù)，這與以往研究相符[12，23，32]，表明冠幅與DBH、CR呈正相關；變量HDR的參數(shù)為負，說明冠幅與HDR成負相關，這與前人的研究一致[20，23，32]；變量BAS的參數(shù)為負數(shù)，因為本研究數(shù)據(jù)來自于不同初植密度的林分，且林分年齡相同、立地質量相似，所以在林分株數(shù)密度較大時，林分斷面積相對較大，此時的競爭也較大，從而導致冠幅較小，該變量以往很少被作為協(xié)變量加入冠幅預測模型，但在本研究中，變量BAS的加入具有一定的生物學意義和統(tǒng)計意義。

總之，本研究所構建的雜種落葉松單木冠幅線性混合效應模型有較好的預測精度，該模型對研究區(qū)域的雜種落葉松能進行較好的預測，且確定了每個樣地抽取5個樣本量即可校正隨機效應參數(shù)，從而減少了森林資源調查中的成本，具有實際應用價值，為雜種落葉松林的經(jīng)營管理提供了一定的依據(jù)。此外，該模型所選入的變量DBH和BAS在實際地面調查中較容易獲取，變量CR和HDR需要調查林木樹高，目前無人機雷達數(shù)據(jù)獲取樹高較方便，所以該模型有較好的實用性。但由于數(shù)據(jù)的局限性，對于其它區(qū)域的預測效果還需要進一步驗證，在今后的研究中可基于多區(qū)域多林齡的林分進行研究，擴大模型的應用范圍。