婁明華 張會(huì)儒 雷相東 李春明 臧 顥

(1.中國(guó)林業(yè)科學(xué)研究院資源信息研究所 北京 100091; 2.寧波市農(nóng)業(yè)科學(xué)研究院 寧波 315040; 3.江西農(nóng)業(yè)大學(xué)林學(xué)院 南昌 330045)

?

基于空間自相關(guān)的天然蒙古櫟闊葉混交林林木胸徑-樹高模型*

婁明華1,2張會(huì)儒1雷相東1李春明1臧 顥3

(1.中國(guó)林業(yè)科學(xué)研究院資源信息研究所 北京 100091; 2.寧波市農(nóng)業(yè)科學(xué)研究院 寧波 315040; 3.江西農(nóng)業(yè)大學(xué)林學(xué)院 南昌 330045)

【目的】 考慮林木間的空間自相關(guān)，構(gòu)建基于空間自相關(guān)的林木胸徑-樹高模型，為可持續(xù)經(jīng)營(yíng)天然混交林提供理論依據(jù)?！痉椒ā?以天然蒙古櫟闊葉混交林為研究對(duì)象，選擇適宜的線性化林木胸徑-樹高模型為基礎(chǔ)模型(BM)，應(yīng)用3個(gè)同步自回歸(SAR)模型即空間滯后模型(SLM)、空間誤差模型(SEM)和空間Durbin模型(SDM)研究該混交林的林木胸徑-樹高模型。同時(shí)，將Delaunay三角網(wǎng)(DT)矩陣、逆距離一次冪(ID1)、逆距離二次冪(ID2)、逆距離五次冪(ID5)、球狀變異函數(shù)(SV)矩陣、高斯變異函數(shù)(GV)矩陣和指數(shù)變異函數(shù)(EV)矩陣共7個(gè)空間加權(quán)矩陣應(yīng)用于SAR模型中。利用普通最小二乘法(OLS)估計(jì)BM參數(shù)，極大似然法估計(jì)3個(gè)SAR模型參數(shù)，并對(duì)4個(gè)模型的回歸參數(shù)進(jìn)行T檢驗(yàn)，對(duì)3個(gè)SAR模型的自回歸參數(shù)進(jìn)行似然比檢驗(yàn)。選擇Moran’s I(MI)指數(shù)比較分析BM、SLM、SDM和SEM 4個(gè)模型的殘差空間自相關(guān)，選擇決定系數(shù)(R2)、均方根誤差(RMSE)和Akaike信息準(zhǔn)則(AIC)3個(gè)擬合指標(biāo)比較分析4個(gè)模型的擬合效果。【結(jié)果】 空間加權(quán)矩陣SV的BM殘差MI值大于1，因此以下結(jié)果分析中將不再考慮SV。其他6個(gè)空間加權(quán)矩陣的BM和SLM殘差MI值均顯著大于期望值I0，但SLM殘差MI值較相同空間加權(quán)矩陣的BM殘差MI值小。除了GV和ID1外，其他4個(gè)空間加權(quán)矩陣的SDM殘差MI值均與I0差異不顯著。除了ID1外，其他5個(gè)空間加權(quán)矩陣的SEM殘差MI值均與I0差異不顯著。3個(gè)SAR模型的3個(gè)擬合指標(biāo)均優(yōu)于BM。在相同的空間加權(quán)矩陣中，SDM和SEM的3個(gè)擬合指標(biāo)非常接近，但均優(yōu)于SLM。在SDM和SEM中，不同空間加權(quán)矩陣(除GV外)根據(jù)3個(gè)擬合指標(biāo)從優(yōu)至劣的排序?yàn)镮D2> DT > ID > ID5 > EV。無論采用哪個(gè)空間加權(quán)矩陣，3個(gè)SAR模型的回歸參數(shù)β1均與BM中的β1相似，且均顯著不為零。相比β1，SEM和BM中的β0相似，但SDM和SLM中的β0與BM中的β0不相似，并且隨著空間加權(quán)矩陣的變化而變化。應(yīng)用于SAR模型的所有空間加權(quán)矩陣中，利用ID1得出的自回歸參數(shù)ρ、γ和λ均明顯高于利用其他空間加權(quán)矩陣計(jì)算的值。GV只有在SEM中才能使自回歸參數(shù)λ顯著。除了GV外，利用其他5個(gè)空間加權(quán)矩陣得出的ρ、γ和λ均顯著?！窘Y(jié)論】 應(yīng)用于SAR模型的7個(gè)空間加權(quán)矩陣中，SV和ID1為不合理的空間加權(quán)矩陣。SLM只能降低模型殘差的空間自相關(guān)，改善模型擬合效果較SDM和SEM差。只要選擇合適的空間加權(quán)矩陣，SDM和SEM就可以消除模型殘差的空間自相關(guān)，提高模型擬合效果，其中ID2是最好的空間加權(quán)矩陣。利用ID2和SEM構(gòu)建以樹種為啞變量的胸徑-樹高模型，從而得出基于空間自相關(guān)的蒙古櫟、楊樺(山楊和白樺)、紅松的胸徑-樹高模型。關(guān)鍵詞： 空間自相關(guān)； 胸徑； 樹高； 自回歸模型； 空間滯后

Individual Diameter-Height Models for MixedQuercus

胸徑和樹高是森林調(diào)查及經(jīng)營(yíng)管理的2個(gè)重要因子(Leietal., 2009)。相對(duì)于胸徑，樹高的測(cè)量比較困難，存在成本高、耗時(shí)長(zhǎng)、觀測(cè)誤差大、視覺阻礙(如樹梢頂部容易遮擋，特別是天然混交林)等不利因素(Chhetrietal., 1996; Colbertetal., 2002; Meht?talo, 2004; Trincadoetal., 2007; Krisnawatietal., 2010; Sharmaetal., 2015)，因此，學(xué)者們經(jīng)常利用林木胸徑-樹高模型估計(jì)樹高，進(jìn)而估計(jì)林木材積、立地質(zhì)量、林分動(dòng)態(tài)以及其他有關(guān)森林生長(zhǎng)與收獲(Leietal., 2001; Peng, 2001; Calamaetal., 2004; Leietal., 2009)。

傳統(tǒng)林木胸徑-樹高模型均假定林木間是相互獨(dú)立的，通常采用普通最小二乘法(ordinary least squares, OLS)進(jìn)行估計(jì)。然而，隨著研究的逐步深入，學(xué)者們逐漸意識(shí)到林木間的關(guān)系并不是相互獨(dú)立的(Frodetal., 1981; Tomppo, 1986; Mateuetal., 1998)，這是因?yàn)榱帜鹃g的競(jìng)爭(zhēng)和微環(huán)境結(jié)合作用，導(dǎo)致林木間存在空間自相關(guān)(spatial autocorrelaton)(Magnussen, 1990; Foxetal., 2001)?？臻g自相關(guān)是指一些變量在同一分布區(qū)內(nèi)的觀測(cè)數(shù)據(jù)之間潛在的相互依賴性，一般認(rèn)為，相鄰林木間的空間自相關(guān)程度強(qiáng)于相離林木間的空間自相關(guān)程度(Legendre, 1993; Legendreetal., 1998),因此，忽略空間自相關(guān)將會(huì)使OLS回歸違背殘差獨(dú)立分布假設(shè)，導(dǎo)致犯第一類錯(cuò)誤(原假設(shè)為真，假設(shè)檢驗(yàn)拒絕了原假設(shè))的可能性變大(Legendre, 1993; Lennon, 2000; Legendreetal., 2002; Dormannetal., 2007; Kisslingetal., 2008)、模型參數(shù)標(biāo)準(zhǔn)差的有偏估計(jì)及回歸模型估計(jì)有效性降低(Kr?mer, 1980; Westetal., 1984; Gregoire, 1987; Foxetal., 2001; Zhangetal., 2004b)等。

鑒于林木間存在空間自相關(guān)，目前已有學(xué)者將空間自相關(guān)模型加入到林木胸徑-樹高模型中，如地理加權(quán)回歸(geographically weighted regression)模型(Zhangetal., 2004a; 2004b; 2005)、同步自回歸(simultaneous autoregressive, SAR)模型(Mengetal., 2009; Luetal., 2011)等，其中，SAR模型全面考慮了空間自相關(guān)過程發(fā)生在模型應(yīng)變量、模型解釋變量及模型誤差中的可能性，是常用的空間自相關(guān)模型(Kisslingetal., 2008; Mengetal., 2009; Luetal., 2011)。具體來說，SAR模型包括空間滯后模型(spatial lag model, SLM)、空間誤差模型(spatial error model, SEM)和空間Durbin模型(spatial Durbin model, SDM)，SLM假設(shè)空間自回歸過程只發(fā)生在應(yīng)變量中，SEM假設(shè)該過程只發(fā)生在模型誤差項(xiàng)中，SDM假設(shè)該過程同時(shí)發(fā)生在應(yīng)變量和解釋變量中(Cressie, 1993; Haining, 2003)。

蒙古櫟(Quercusmongolica)闊葉混交林是我國(guó)東北地區(qū)主要的天然次生林，蒙古櫟主干堅(jiān)硬耐腐、紋理美觀，是良好的造材原料，枝條發(fā)熱量高，可作為良好的薪炭材，葉可飼蠶和動(dòng)物，根系發(fā)達(dá)且適應(yīng)性強(qiáng)，具有涵養(yǎng)水源、 防風(fēng)固沙、 減少徑流、保持水土等作用，果實(shí)富含淀粉，經(jīng)加工提煉可轉(zhuǎn)化燃料乙醇及生物化工產(chǎn)品(田超等， 2011； 馬武， 2012)，對(duì)國(guó)民經(jīng)濟(jì)和生態(tài)環(huán)境均具有很高的價(jià)值(于順利等， 2001； 2004； 許中旗等， 2006； 陳新美等， 2010； 洪玲霞等， 2012)。為此，本文以天然蒙古櫟闊葉混交林為研究對(duì)象，考慮林木間的空間自相關(guān)，構(gòu)建基于空間自相關(guān)的林木胸徑-樹高模型，以期為可持續(xù)經(jīng)營(yíng)天然混交林提供理論依據(jù)。

1 材料與方法

1.1 研究區(qū)概況與樣地調(diào)查

研究區(qū)位于吉林省汪清林業(yè)局塔子溝林場(chǎng)，地理坐標(biāo)為129.971°—130.222°E，43.327°—43.492°N。屬長(zhǎng)白山系老爺嶺山脈雪嶺支脈，低山丘陵地貌，海拔300～1 200 m，坡度5°～25°。研究區(qū)屬季風(fēng)性氣候，年均氣溫3.9 ℃左右，年均降水量600～700 mm，且集中在夏季，占全年總降水量的80%。土壤主要以暗棕壤為主，平均厚度40 cm左右。植被屬長(zhǎng)白山植物區(qū)系，群落結(jié)構(gòu)復(fù)雜，植物種類較多，主要樹種有蒙古櫟、白樺(Betulaplatyphylla)、紅松(Pinuskoraiensis)、山楊(Populusdavidiana)、榆(Ulmuspumila)、椴(Tiliatuan)、黑樺(Betuladahurica)、長(zhǎng)白落葉松(Larixolgensis)、色木槭(Acermono)、紅皮云杉(Piceakoraiensis)、冷杉(Abiesfabri)、水曲柳(Fraxinusmandshurica)、黃檗(Phellodendronamurense)等。

2013年7—9月，在塔子溝林場(chǎng)選擇典型的天然蒙古櫟闊葉混交林，設(shè)置100 m×100 m固定樣地，樣地坡度8°，坡位為中坡，坡向?yàn)槲髌?，樣地中心海?05 m，樣地左下角坐標(biāo)為130.112°E，43.450°N。利用森林羅盤儀，將樣地劃分100個(gè)10 m×10 m的網(wǎng)格單元，在每個(gè)網(wǎng)格單元內(nèi)對(duì)胸徑5 cm以上的樹木進(jìn)行每木調(diào)查，記錄樹種，測(cè)量胸徑、樹高、枝下高、東西南北冠幅等因子，并利用徠卡激光測(cè)距儀測(cè)定林木坐標(biāo)(X，Y)。樣地林分特征詳見表1。

1.2 基礎(chǔ)林木胸徑-樹高模型

由表1可知，天然蒙古櫟闊葉混交林中有4個(gè)主要樹種，即蒙古櫟、白樺、紅松和山楊?？紤]到山楊和白樺具有相近的林木生長(zhǎng)特性，本文將山楊和白樺稱為楊樺，作為同一建模樹種，并建立全林分、蒙古櫟、楊樺和紅松4個(gè)胸徑-樹高模型。為不改變林木間的空間關(guān)系即林木間的空間自相關(guān)結(jié)構(gòu)，采用啞變量方法即以樹種為啞變量構(gòu)建分樹種的胸徑-樹高模型(簡(jiǎn)稱啞變量模型)。

當(dāng)前SAR模型均為線性的(Mengetal., 2009)，因此本文只考慮線性化的基礎(chǔ)林木胸徑-樹高模型。根據(jù)圖1所示的胸徑與樹高散點(diǎn)圖，選用3個(gè)備選線性化基礎(chǔ)模型，具體見表2。

1.3 同步自回歸模型

SAR模型包括空間滯后模型(SLM)、空間誤差模型(SEM)和空間Durbin模型(SDM)，其表達(dá)式(Cliffetal., 1981; Anselin, 1988; Haining, 2003; Fortinetal., 2005; Dormannetal., 2007)如下：

(1)

(2)

(3)

式中：Y表示應(yīng)變量矩陣；X表示解釋變量矩陣；W表示空間加權(quán)矩陣；β表示回歸系數(shù)；ε表示模型殘差；ρ、λ和γ表示空間自回歸系數(shù)；ξ表示空間相關(guān)殘差。

表1 固定樣地林分特征①Tab.1 Stand characteristics of permanent plots

①N表示林木株數(shù)； BA表示斷面積； SBA表示樹種斷面積比例；Dmean、Dmax和Dmin分別表示平均胸徑、最大胸徑和最小胸徑；Hmean、Hmax和Hmin分別表示平均樹高、最大樹高和最小樹高。Ndenote tree number; BA denote basal area; SBA denote species basal area ratio;Dmean,DminandDmaxdenote mean, minimum and maximum diameter at breast height;Hmean,HminandHmaxdenote mean, minimum and maximum tree height.

圖1 不同樹種的胸徑-樹高散點(diǎn)示意Fig.1 Scatter plot of diameter at breast height against tree height for different species

模型序號(hào)ModelNo．模型表達(dá)式Modelexpression參考文獻(xiàn)ReferencesAlnH=β0+β1lnD+εCurtis，1967BlnH=β0+β1D+εHuangetal．，1992；李希菲等，1994．ClnH=β0+β1D-1+εCurtis，1967

①H表示樹高；D表示胸徑；β0和β1表示回歸系數(shù)；ε表示模型殘差。Hdenotes tree height;Ddenotes diameter at breast height;β0andβ1denote regression coefficients;εdenotes model error term.

WY、WX和Wξ為3個(gè)空間滯后變量。當(dāng)ρ=0時(shí)，式(1)變?yōu)榫€性回歸模型； 當(dāng)λ=0時(shí)，式(2)變?yōu)榫€性回歸模型； 當(dāng)ρ和γ同時(shí)為零時(shí)，式(3)變?yōu)榫€性回歸模型。采用極大似然法(maximum likelihood)估計(jì)3個(gè)SAR模型的參數(shù)，通過R統(tǒng)計(jì)語言spdep包的lagsarlm和errorsarlm函數(shù)實(shí)現(xiàn)。

1.4 空間加權(quán)矩陣

空間加權(quán)矩陣(spatial weight matrix)是SAR模型的重要組成部分，本文選取常用的空間加權(quán)矩陣，即相接鄰近(contiguous neighbors)矩陣、逆距離冪(inverse distances raised to some power)矩陣和地統(tǒng)計(jì)(geostatistical)矩陣(Getisetal., 2004; Luetal., 2011)。

相接鄰近矩陣包括Rook相接矩陣、Queen相接矩陣、Bishop相接矩陣(Anselin, 2001)和Delaunay三角網(wǎng)(Delaunay triangulation, DT)矩陣(Smirnovetal., 2001)等。其中，Rook相接矩陣、Queen相接矩陣和Bishop相接矩陣一般用于規(guī)則的空間相接關(guān)系，DT矩陣可處理不規(guī)則的空間相接關(guān)系，因?yàn)榱帜鹃g的空間關(guān)系一般是不規(guī)則的，因此本文采用DT矩陣。DT矩陣中，林木間為相接鄰近關(guān)系時(shí)，其空間加權(quán)值為1；否則為0。同時(shí)，利用R統(tǒng)計(jì)語言spdep包的trib2nb函數(shù)計(jì)算空間加權(quán)值。

逆距離冪矩陣包括逆距離一次冪(inverse distances raised to one power, ID1)、逆距離二次冪(inverse distances raised to two power, ID2)和逆距離五次冪(inverse distances raised to five power, ID5)3種形式(Getisetal., 2004)，其表達(dá)式如下：

(4)

式中：dij表示林木i與其相鄰木j的距離，wij表示林木i與其相鄰木j的空間加權(quán)值，當(dāng)n=1時(shí)，則為ID1； 當(dāng)n=2時(shí)，則為ID2； 當(dāng)n=5時(shí)，則為ID5。

Getis等(2004)利用滿足本征假設(shè)(intrinsic stationarity)的變異函數(shù)構(gòu)建了3個(gè)地統(tǒng)計(jì)矩陣，即球狀變異函數(shù)(spherical variogram, SV)矩陣、高斯變異函數(shù)(Gaussian variogram, GV)矩陣和指數(shù)變異函數(shù)(exponential variogram, EV)矩陣，其表達(dá)式如下：

(5)

(6)

(7)

式中：r表示變程，wij和dij意義同上。

利用R統(tǒng)計(jì)語言nlme包的corExp、corGaus和corSpher函數(shù)分別計(jì)算SV、GV和EV矩陣的空間加權(quán)值。

本文采用DT、ID1、ID2、ID5、SV、GV和EV共7個(gè)空間加權(quán)矩陣，應(yīng)用于3個(gè)SAR模型中。

1.5 空間自相關(guān)檢驗(yàn)指標(biāo)

選擇常用的Moran’s I(MI)指數(shù)檢驗(yàn)?zāi)Ｐ蜌埐畹目臻g自相關(guān)，其表達(dá)式(Paradis, 2014)如下：

(8)

其原假設(shè)為I等于期望值I0=-1/(n-1)，表明不存在空間自相關(guān)。如果I顯著(α=0.01)大于I0，則模型殘差為正自相關(guān)，說明相鄰林木間的相似程度高； 反之如果I顯著(α=0.01)小于I0，則模型殘差為負(fù)自相關(guān)，說明相鄰林木間的不相似程度高。利用R統(tǒng)計(jì)語言spdep包的moran.test函數(shù)進(jìn)行檢驗(yàn)。1.6 模型擬合指標(biāo)

為便于比較模型之間的擬合程度，本文選用決定系數(shù)(coefficient of determination,R2)、均方根誤差(root mean square error, RMSE)和Akaike信息準(zhǔn)則(Akaike information criterion, AIC)3個(gè)擬合指標(biāo)比較分析模型的擬合效果。

2 模型比較分析

從3個(gè)備選模型中選擇全林分模型和啞變量模型建模效果均最優(yōu)的線性化基礎(chǔ)模型作為全林分和啞變量的基礎(chǔ)模型(base model, BM)。值得注意的是，全林分模型和啞變量模型中林木間的空間關(guān)系是一致的，即林木間具有一致的空間自相關(guān)結(jié)構(gòu)，因此，只需建立基于全林分基礎(chǔ)模型的SLM、SEM和SDM共3個(gè)SAR模型，并通過比較分析模型殘差空間自相關(guān)、模型擬合及模型參數(shù)估計(jì)效果，即可找到最適宜的SAR模型。最后，利用最適宜的SAR模型構(gòu)建啞變量SAR模型。

2.1 基礎(chǔ)模型的確定

利用3個(gè)擬合指標(biāo)R2、RMSE和AIC，分別比較全林分和啞變量模型的3個(gè)備選基礎(chǔ)模型。由表3可知，全林分模型的備選基礎(chǔ)模型C具有最大的R2、最小的RMSE和AIC；啞變量模型的備選模型A具有最大的R2、最小的RMSE和AIC，模型C略低于模型A，最差為模型B。由表4可知，啞變量模型的備選模型A中主要樹種蒙古櫟和楊樺的參數(shù)β0和β1均顯著(P<0.01)，主要樹種紅松的參數(shù)β0顯著(P<0.01)，而參數(shù)β0不顯著(P>0.1)。因此，本文選用模型C作為全林分和啞變量的基礎(chǔ)胸徑-樹高模型(BM)。

表3 3個(gè)備選基礎(chǔ)模型的擬合指標(biāo)及參數(shù)估計(jì)①Tab. 3 Fitting index and parameter estimates for three alternative base models

①R2：決定系數(shù)Coefficient of determination;RMSE:均方根誤差Root mean square error;AIC:Akaike信息準(zhǔn)則Akaike information criterion.下同The same below.

表4 啞變量模型中3個(gè)主要樹種的參數(shù)估計(jì)Tab. 4 Parameter estimates for three main species in dummy variable model

2.2 全林分模型殘差空間自相關(guān)

由表5可知，空間加權(quán)矩陣SV的BM殘差MI值大于1，說明SV是不合理的空間加權(quán)矩陣，因此，以下結(jié)果分析中將不再考慮SV。除了SV外，其他6個(gè)空間加權(quán)矩陣的BM殘差MI值均顯著(P< 0.01)大于期望值I0=-1/(752-1)=-0.001 3，說明BM殘差存在正空間自相關(guān)。6個(gè)空間加權(quán)矩陣的SLM殘差MI值均顯著(P< 0.01)大于I0，說明SLM沒有消除模型殘差的空間自相關(guān)，但是降低了空間自相關(guān)影響，因其MI值小于相對(duì)應(yīng)的BM殘差MI值。除了GV和ID1外，其他4個(gè)空間加權(quán)矩陣的SDM殘差MI值均與I0差異不顯著(P> 0.1)，說明選擇合適的空間加權(quán)矩陣，SDM可以消除模型殘差的空間自相關(guān)，同時(shí)也表明產(chǎn)生模型殘差的空間自相關(guān)可能是由滯后變量WY(即相鄰加權(quán)樹高)和滯后變量WX(即相鄰加權(quán)胸徑)共同作用產(chǎn)生的。除了ID1外，其他5個(gè)空間加權(quán)矩陣的SEM殘差MI值均與I0差異不顯著(P> 0.1)，說明SEM與SDM類似，在合適的空間加權(quán)矩陣條件下，SEM可以消除模型殘差的空間自相關(guān)，但空間自相關(guān)可能是由滯后變量Wξ(即相鄰加權(quán)未觀測(cè)的變量)產(chǎn)生的，因?yàn)镾EM的前提假設(shè)是空間自相關(guān)過程來源于模型誤差。

此外，由表5還可知，SDM和SEM 2個(gè)SAR模型均能消除模型殘差的空間自相關(guān)矩陣中，依據(jù)MI值和P值，GV是最好的空間加權(quán)矩陣，因其有最大P值，但只對(duì)SEM有效，其次為EV和ID5，但I(xiàn)D5略遜于EV，最后為ID2和DT。

表5 4個(gè)模型殘差的空間自相關(guān)①Tab. 5 Test for spatial autocorrelation of four model residuals

① W表示空間加權(quán)矩陣。下同。 W denotes spatial weight matrix.The same below.

2.3 全林分模型擬合

利用3個(gè)擬合指標(biāo)R2、RMSE和AIC，比較分析BM與3個(gè)SAR模型(SLM、SDM和SEM)的擬合效果。由表6和表3可知，3個(gè)SAR模型的R2(0.578，0.669)均大于BM，RMSE和AIC均小于BM，說明3個(gè)SAR模型的擬合效果均優(yōu)于BM。3個(gè)SAR模型中，SDM與SEM擬合效果不相上下，且均明顯優(yōu)于SLM，這是因?yàn)镾LM沒有消除模型殘差的空間自相關(guān)(詳見2.2)。另外，空間加權(quán)矩陣GV在SDM中的擬合效果明顯比SEM差，結(jié)合表5可知，這是因?yàn)镚V在SDM中沒有消除模型殘差的空間自相關(guān)。因此，消除模型殘差的空間自相關(guān)，有助于改善模型擬合效果。

由表6還可知，依據(jù)R2、RMSE和AIC，在SDM和SEM中的空間加權(quán)矩陣(不考慮GV)從優(yōu)至劣排序?yàn)椋?ID2> DT > ID > ID5 > EV。

表6 3個(gè)SAR模型的擬合Tab. 6 Three SAR models fitting

2.4 全林分模型參數(shù)

對(duì)3個(gè)SAR模型及BM的回歸參數(shù)進(jìn)行T檢驗(yàn)，對(duì)3個(gè)SAR模型的自回歸參數(shù)進(jìn)行似然比檢驗(yàn)(likelihood ratio test)。由表7可知，無論采用哪個(gè)空間加權(quán)矩陣，3個(gè)SAR模型的回歸參數(shù)β1均與BM中的β1相似，且均顯著不為零(P< 2.2e-16)。相比β1，SEM與BM中的β0相似，但SDM和SLM中的β0與BM中的β0不相似，并且隨著空間加權(quán)矩陣的變化而變化，這主要是由于SDM和SLM的模型形式中增加了空間滯后變量WY和WX。所有的β0中，ID1應(yīng)用于SDM中的β0與零差異不顯著(P>0.1)，說明ID1在SDM中是不理想的空間加權(quán)矩陣。同時(shí)由表7可知，利用ID1得出的自回歸參數(shù)ρ、γ和λ均明顯高于其他空間加權(quán)矩陣得出的值，且無論應(yīng)用于哪個(gè)SAR模型，這進(jìn)一步說明了ID1是一個(gè)不理想的空間加權(quán)矩陣。

另外，由表8可知，GV只有在SEM中才能使自回歸參數(shù)λ顯著(P<0.01)。除了GV外，利用其他5個(gè)空間加權(quán)矩陣得出的ρ、γ和λ均顯著(P<0.01)，且均大于零，說明空間滯后變量WY、WX及Wξ均能產(chǎn)生正的空間自相關(guān)。

表7 3個(gè)SAR模型的回歸參數(shù)①Tab. 7 Regressive parameters of three SAR models

①BM:β0=3.014,P< 2.2e-16;β1=-7.967,P< 2.2e-16.

表8 3個(gè)SAR模型的自回歸參數(shù)Tab. 8 Autoregressive parameters of three SAR models

2.5 啞變量同步自回歸模型

由基于全林分基礎(chǔ)模型的SLM、SDM和SEM的殘差MI檢驗(yàn)效果、建模擬合效果和參數(shù)檢驗(yàn)效果可知，以ID2為空間加權(quán)矩陣的SEM模型是清除基礎(chǔ)模型殘差空間自相關(guān)現(xiàn)象和提高模型精度最佳的同步自回歸模型，因此，本文采用該模型建立啞變量SAR模型，得到其R2、RMSE和AIC分別為0.766、0.237和10.728，3個(gè)擬合指標(biāo)均優(yōu)于啞變量基礎(chǔ)模型(表3)，其自回歸參數(shù)λ為0.412，且顯著不為零(P<0.01)。啞變量SAR模型中的蒙古櫟、楊樺和紅松的參數(shù)估計(jì)結(jié)果見表9。由表9可知，3個(gè)主要樹種的截距參數(shù)β0和β1均顯著不為零(P<0.01)，參數(shù)檢驗(yàn)效果理想。以ID2和SEM建立的全林分SAR模型殘差、啞變量SAR模型殘差以及3個(gè)主要樹種的模型殘差見圖2。

表9 啞變量模型中3個(gè)主要樹種的參數(shù)估計(jì)Tab. 9 Parameter estimates for three main species in dummy variable model

圖2 不同模型殘差Fig.2 Residuals of different models A表示全林分SAR模型； B表示啞變量SAR模型； C表示主要樹種蒙古櫟； D表示主要樹種楊樺； E表示主要樹種紅松。A denotes whole stand SAR model; B denotes Dummy variable SAR model; C denotes main species of Quercus mongolica; D denotes main species of Populus davidiana and Betula platyphylla; E denotes main species of Pinus koraiensis.

3 結(jié)論

本文以線性化林木胸徑-樹高模型為基礎(chǔ)模型(BM)，加入3個(gè)SAR模型(SDM、SEM和SLM)構(gòu)建了基于空間自相關(guān)的天然蒙古櫟闊葉混交林林木胸徑-樹高模型，得到以下結(jié)論：

1) SDM和SEM能很好地處理空間自相關(guān)，而SLM最差，這一結(jié)論與Kissling等(2008)、Meng等(2009)和Lu等(2011)研究結(jié)果類似。

2) 無論采用哪個(gè)空間加權(quán)矩陣，SLM均不能消除模型殘差空間自相關(guān)，但可降低空間自相關(guān)，在一定程度上提高了模型擬合效果。

3) 模型殘差MI值均與I0差異不顯著，SDM和SLM雖能消除模型殘差空間自相關(guān)，提高模型擬合效果，但需要選擇合適的空間加權(quán)矩陣，如可選擇ID2、DT、ID5和EV，其中ID2最好，能將全林分基礎(chǔ)模型的R2從0.574提高到0.669，啞變量基礎(chǔ)模型的R2從0.732提高到0.766。

4) 所有7個(gè)空間加權(quán)矩陣中，SV和ID1應(yīng)用于胸徑-樹高模型時(shí)，均為不合理的空間加權(quán)矩陣。

5) 蒙古櫟闊葉混交林林木胸徑-樹高模型殘差的空間自相關(guān)可能由滯后變量WY(即相鄰加權(quán)樹高)和WX(即相鄰加權(quán)胸徑)共同產(chǎn)生或者可能由滯后變量Wξ(即相鄰加權(quán)未觀測(cè)的變量)產(chǎn)生。

6) 利用ID2和SEM構(gòu)建了啞變量模型，從而得出蒙古櫟、楊樺和紅松的同步自回歸胸徑-樹高模型。

本文研究了基于空間自相關(guān)的天然蒙古櫟闊葉混交林林木胸徑-樹高模型，得出ID2和SEM是最適宜的空間自相關(guān)模型。對(duì)于其他林分類型，該模型是否最適宜是下一步需要探討的問題。另外，本文只構(gòu)建了基于線性化空間自相關(guān)模型的天然蒙古櫟闊葉混交林林木胸徑-樹高模型，如何研究非線性化空間自相關(guān)模型就目前來說是一個(gè)難題，也是今后的重點(diǎn)研究方向。

陳新美, 張會(huì)儒, 姜慧泉. 2010. 東北過伐林區(qū)蒙古櫟林空間結(jié)構(gòu)分析與評(píng)價(jià). 西南林學(xué)院學(xué)報(bào), 30(6): 20-24.

(Chen X M, Zhang H R, Jiang H Q. 2010. Analysis and evaluation on spatial structure ofQuercusmongolicaforests in over-logged region in northeast China. Journal of Southwest Forestry University, 30(6): 20-24. [in Chinese])

洪玲霞, 雷相東, 李永慈. 2012. 蒙古櫟林全林整體生長(zhǎng)模型及其應(yīng)用. 林業(yè)科學(xué)研究, 25(2): 201-206.

(Hong L X, Lei X D, Li Y C. 2012. Integrated stand growth model of mongolian oak and its application. Forest Research, 25(2): 201-206. [in Chinese])

李希菲, 唐守正, 袁國(guó)仁,等. 1994. 自動(dòng)調(diào)控樹高曲線和一元立木材積模型. 林業(yè)科學(xué)研究, 7(5): 512-518.

(Li X F, Tang S Z, Yuan G R,etal. 1994. Self-adjusted height-diameter curves and one entry volume model. Forest Research, 7(5): 512-518. [in Chinese])

馬 武. 2012. 蒙古櫟林單木生長(zhǎng)模型系統(tǒng). 北京: 中國(guó)林業(yè)科學(xué)研究院碩士學(xué)位論文.

(Ma W. 2012. Individual tree growth model system for mongolian oak forest. Beijing: MS thesis of Chinese Academy of Forestry. [in Chinese])

田 超, 楊新兵, 李 軍, 等. 2011. 冀北山地不同海拔蒙古櫟林枯落物和土壤水文效應(yīng). 水土保持學(xué)報(bào), 25(4): 221-226.

(Tian C, Yang X B, Li J,etal. 2011. Hydrological effects of forest litters and soil ofQuercusmongolicain the different altitudes of north mountain of Hebei province. Journal of Soil and Water Conservation, 25(4): 221-226. [in Chinese])

許中旗, 李文華, 劉文忠, 等. 2006. 我國(guó)東北地區(qū)蒙古櫟林生物量及生產(chǎn)力的研究. 中國(guó)生態(tài)農(nóng)業(yè)學(xué)報(bào), 14(3): 21-24.

(Xu Z Q, Li W H, Liu W Z,etal. Study on the biomass and productivity of mongolian oak forests in northeast region of China. Chinese Journal of Eco-Agriculture, 14(3): 21-24. [in Chinese])

于順利, 馬克平, 陳靈芝, 等. 2001. 黑龍江省不同地點(diǎn)蒙古櫟林生態(tài)特點(diǎn)研究. 生態(tài)學(xué)報(bào), 21(1): 41-46.

(Yu S L, Ma K P, Chen L Z,etal. 2001. The ecological characteristics ofQuercusmongolicaforest in Heilongjiang Province. Acta Ecologica Sinica, 21(1): 41-46. [in Chinese])

于順利, 馬克平, 徐存寶, 等. 2004. 環(huán)境梯度下蒙古櫟群落的物種多樣性特征. 生態(tài)學(xué)報(bào), 24(12): 2932-2939.

(Yu S L, Ma K P, Xu C B,etal. 2004. The species diversity characteristics comparison ofQuercusmongolicacommunity along environmental gradient factors. Acta Ecologica Sinica, 24(12): 2932-2939. [in Chinese])

Anselin L. 1988. Spatial econometrics: methods and models. Dordrecht: Kluwer Academic Publishers.

Anselin L. 2001.Spatial effects in econometric practice in environment and resource economics. American Journal of Agricultural Economics, 83(3): 705-710.

Calama R, Montero G. 2004. Interregional nonlinear height-diameter model with random coefficients for stone pine in Spain. Canadian Journal of Forest Research, 34(1): 150-163.

Chhetri D B K, Fowler G W. 1996. Prediction models for estimating total heights of trees from diameter at breast height measurements in Nepal’s lower temperate broad-leaved forests. Forest Ecology and Management, 84(1/3): 177-186.

Cliff A D, Ord J K. 1981. Spatial processes-models and applications. London: Pion Ltd.

Colbert K C, Larsen D R, Lootens J R. 2002. Height-diameter equations for thirteen midwestern bottomland hardwoods species. Northern Journal of Applied Forestry, 19(4): 171-176.

Cressie N A C. 1993. Statistics for spatial data. Wiley Series in Probability and Mathematical Statistics, New York: Wiley.

Curtis R O. 1967. Height-diameter and height-diameter-age equations for second-growth douglas-fir. Forest Science, 13(4): 365-375.

Dormann C F, McPherson J M, Araújo M B,etal. 2007. Methods to account for spatial autocorrelation in analysis of species distributional data: a review. Ecography, 30(5): 609-628.

Frod E D, Diggle P J. 1981. Competition for light in a plant monoculture modelled as a spatial stochastic process. Annals of Botany, 48(4): 481-500.

Fortin M J, Dale M R T. 2005. Spatial analysis-guide for ecologists. Cambridge: Cambridge University Press.

Fox J C, Ades P K, Bi H. 2001. Stochastic structure and individual-tree growth models. Forest Ecology and Management, 154(1): 261-276.

Getis A, Aldstadt J. 2004. Constructing the spatial weight matrix using a local statistic. Geographical Analysis, 36(2): 90-104.

Gregoire T G. 1987. Generalized error structure for forestry yield models. Forest Science, 33(2): 423-444.

Haining R. 2003. Spatial data analysis: theory and practice. Cambridge: Cambridge University Press.

Huang S, Titus S, Wiens D P. 1992. Comparison of nonlinear height-diameter functions for major alberta tree species. Canadian Journal of Forest Research, 22(9): 1297-1304.

Kissling W D, Carl G. 2008. Spatial autocorrelation and the selection of simultaneous autoregressive models. Global Ecology and Biogeography, 17(1): 59-71.

Kr?mer W. 1980. Finite sample efficiency of ordinary least squares in linear regression model with autocorrelated errors. Journal of the American Statistical Association, 75(372): 1005-1009.

Krisnawati H, Wang Y, Ades P K. 2010. Generalized height-diameter models forAcaciamangiumWilld. plantations in south Sumatra. Journal of Forestry Research, 7(1): 1-19.

Legendre P, Dale M R T, Fortin M J,etal. 2002. The consequences of spatial structure for the design and analysis of ecological field surveys. Ecography, 25(5): 601-615.

Legendre P, Legendre L. 1998. Numerical ecology. Amsterdam: Elsevier.

Legendre P. 1993. Spatial autocorrelation: trouble or new paradigm? Ecology, 74(6): 1659-1673.

Lei X D, Peng C H, Wang H Y,etal. 2009. Individual height-diameter models for young black spruce (Piceamariana) and jack pine (Pinusbanksiana) plantations in New Brunswick, Canada. The Forestry Chronicle, 85(1): 43-56.

Lei Y C, Parresol B R. 2001. Remarks on height-diameter modeling. Research Note SE-IO. U.S. Department of Agriculture, Forest Service, Southeastern Forest Experiment Station, Asheville, NC.

Lennon J J. 2000. Red-shifts and red herrings in geographical ecology. Ecography, 23(1): 101-113.

Lu J F, Zhang L J. 2011. Modeling and prediction of tree height-diameter relationship using spatial autoregressive models. Forest Science, 57(3): 252-264.

Magnussen S. 1990. Application and comparison of spatial models in analyzing tree-genetics fiele trials. Canadian Journal of Forest Research, 20(5): 536-546.

Mateu J, Uso J L, Montes F.1998. The spatial pattern of a forest ecosystem. Ecological Modelling, 108(1): 163-174.

Meht?talo L. 2004. A longitudinal height-diameter model for Norway spruce in Finland. Canadian Journal of Forest Research, 34(1): 131-140.

Meng Q M, Cieszewski C J, Strub M R,etal. 2009. Spatial regression modeling of tree height-diameter relationships. Canadian Journal of Forest Research, 39(12): 2283-2293.

Paradis E. 2014. Moran’s autocorrelation coefficient in comparative methods.[2015-03-14]. http://www.cran.r-project.org/web/packages/ape/vignettes/MoranI.pdf.

Peng C H. 2001. Developing and validating nonlinear height-diameter models for major tree species of Ontario’s boreal forest. Northern Journal of Applied Forestry, 18(3): 87-94.

Sharma R P, Breidenbach J. 2015. Modeling height-diameter relationships for norway spruce, scots pine, and downy birch using Norwegian national forest inventory data. Forest Science and Technology, 11(1): 44-53.

Smirnov O, Anselin L. 2001.Fast maxinum likehood estimation of very large spatial autoregressive models:a characteristic pdynomial approach.Computational Statistics and Data Analysis,35(3):301-319.

Tomppo E. 1986. Models and methods for analysing spatial patterns of trees. Helsinki: The Finnish Forest Research Institute.

Trincado G, Van der Schaaf C L, Burkhart H E. 2007. Regional mixed-effects height-diameter models for loblolly pine (PinustaedaL.) plantations. European Journal of Forest Research, 126(2): 253-262.

West P W, Ratkowsky D A, Davis A W. 1984. Problems of hypothesis testing of regression with multiple measurements from individual sampling units. Forest Ecology and Management, 7(3): 207-224.

Zhang L J, Bi H Q, Cheng P F,etal. 2004a. Modeling spatial variation in tree diameter-height relationships. Forest Ecology and Management, 189(1): 317-329.

Zhang L J, Shi H J. 2004b. Local modeling of tree growth by geographically weighted regression. Forest Science, 50(2): 225-244.

Zhang L J. Gove J F, Heath L S. 2005. Spatial residual analysis of six modeling techniques. Ecological Modelling, 186(2): 154-177.

(責(zé)任編輯 石紅青)

mongolicaBroadleaved Natural Stands Based onSpatial Autocorrelation

Lou Minghua1,2Zhang Huiru1Lei Xiangdong1Li Chunming1Zang Hao3

(1.ResearchInstituteofForestResourceInformationTechniques,CAFBeijing100091; 2.NingboAcademyofAgricultureSciences,ZhejiangProvinceNingbo315040; 3.CollegeofForestry,JiangxiAgriculturalUniversityNanchang330045)

【Objective】 Considering spatial autocorrelation among individuals, individual diameter-height models based on spatial autocorrelation were constructed. It may provide a theoretical basis for sustainable management of natural mixed forests. 【Method】Three simultaneous autoregressive (SAR) models, including spatial lag model (SLM), spatial error model (SEM) and spatial Durbin model (or called spatial mixed model) (SDM) within seven spatial weight matrices, including Delaunay triangulation (DT), inverse distance raised to one power (ID1), inverse distance raised to two powers (ID2), inverse distance raised to five powers (ID5), spherical variogram (SV), gaussian variogram (GV) and exponential variogram (EV), was used to construct individual diameter at breast height and height models of mixedQuercusmongolicabroadleaved natural stands in Northeast China, and treating linearization base model (BM) as a benchmark model. Model parameters of BM were estimated by ordinary least squares (OLS), model parameters of three SAR models were estimated by maximum likelihood. Model coefficientsβ0andβ1of four models were tested byT-test, the autoregressive parametersρ,γandλwere all tested by likelihood ratio test. Moran’s I (MI) was selected to compared autocorrelation of four model residuals. Three statistics, i.e. coefficient of determination (R2), root mean square error (RMSE) and Akaike information criterion (AIC), were regarded as the appropriate criteria to identify the model fitting among BM, SLM, SDM and SEM. 【Result】MI values of BM residuals were larger than 1, when applying SV into BM. Therefore, SV was the unreasonable spatial weight matrix and did not regard as a spatial weight matrix in the following result analysis. MI values of BM and SLM residuals were significantly larger than the expected valueI0of MI in the all spatial weight matrices (except SV). MI values of SLM residuals were smaller than those of BM using the same spatial weight matrix. The difference between MI values of SDM residuals andI0was not significant in other four spatial weight matrices, except GV and ID1. Similarly, the difference between MI values of SEM residuals andI0was not significant in other five spatial weight matrices, except ID1. Three criteria of three SAR models were all better than those of BM. Using the same spatial weight matrix, MI values of SDM were very similar to those of SEM, meanwhile, MI values of SDM and SEM were both larger than those of SLM. Different spatial weight matrices (except GV) in SDM and SEM were sorted from best to worst according three criteria and the ranking was: ID2 > DT > ID > ID5 > EV. Model coefficientsβ1of three SAR were very similar to those of BM, regardless of which spatial weight matrix was used. Compared withβ1, model coefficientsβ0of SEM were similar to those of BM, while model coefficientsβ0of SDM and SLM were different to those of BM, and were changed along with the different spatial weight matrix. Among all spatial weight matrices within three SAR models, the autoregressive parametersρ,γandλusing ID1 were larger higher than any other spatial weight matrix. GV only applied to SEM, rather than SDM, could make the autoregressive parameterλsignificant not equal to zero. The autoregressive parametersρ,γandλwere all not equal to zero using five spatial weight matrices (except GV).【Conclusion】 Among all spatial weight matrices applied in three SAR models, SV and ID1 are the unreasonable spatial weight matrices. SLM do not remove, but reduce the spatial autocorrelation of model residuals, and slightly improve the model fitting. Model fitting of SLM was worse than those of SDM and SEM. Selecting appropriate spatial weight matrices, SDM and SEM can remove the spatial dependence of model residuals and improve the model fitting. ID2 is the best one among these selected appropriate spatial weight matrices. The diameter-height models ofQuercusmongolica,Populus-Betula(PopulusdavidianaandBetulaplatyphylla) andPinuskoraiensiswere constructed by species dummy variables SAR models based on ID2 and SEM.

spatial autocorrelation; diameter at breast height; tree height; autoregressive model; spatial lag

10.11707/j.1001-7488.20170608

2015-12-07；

2016-04-13。

“十二五”國(guó)家科技支撐計(jì)劃課題(2012BAD22B02)。

S757.2

A

1001-7488(2017)06-0067-10

* 張會(huì)儒為通訊作者。