應(yīng)用混合效應(yīng)法建立的雜種落葉松人工林單木冠幅預(yù)測模型1)

馬愛云 李鳳日 董利虎

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

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

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

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

1 研究區(qū)概況

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

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

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

2 研究方法

2.1 基礎(chǔ)模型構(gòu)建

2.2 線性混合效應(yīng)模型

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

(1)

表2 冠幅基礎(chǔ)模型

圖1 冠幅—胸徑相關(guān)關(guān)系

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

2.3 模型的評價(jià)與檢驗(yàn)指標(biāo)

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

(2)

均方根誤差(RMSE)：

(3)

平均絕對偏差(MAE)：

(4)

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

(5)

模擬效率(EF)：

(6)

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

(7)

2.4 不同樣本量對隨機(jī)效應(yīng)的校正

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

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

3 結(jié)果與分析

3.1 模型的構(gòu)建

3.1.1 單木因子與林分因子

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

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

(8)

3.1.2 混合效應(yīng)模型

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

a4BAS。

(9)

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

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

3.2 模型的評價(jià)與檢驗(yàn)

圖2 冠幅觀測值與模型(8)、模型(9)的預(yù)測值相關(guān)圖

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

3.3 不同樣本量對隨機(jī)效應(yīng)的校正對比

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

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

表5 線性模型與混合效應(yīng)模型檢驗(yàn)結(jié)果

4 結(jié)論與討論

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

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

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

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