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

婁明華 張會儒 雷相東 李春明 臧 顥

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

?

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

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

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

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

Individual Diameter-Height Models for MixedQuercus

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

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

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

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

1 材料與方法

1.1 研究區概況與樣地調查

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

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

1.2 基礎林木胸徑-樹高模型

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

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

1.3 同步自回歸模型

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

(1)

(2)

(3)

式中：Y表示應變量矩陣；X表示解釋變量矩陣；W表示空間加權矩陣；β表示回歸系數；ε表示模型殘差；ρ、λ和γ表示空間自回歸系數；ξ表示空間相關殘差。

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

①N表示林木株數； 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 不同樹種的胸徑-樹高散點示意Fig.1 Scatter plot of diameter at breast height against tree height for different species

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

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

WY、WX和Wξ為3個空間滯后變量。當ρ=0時，式(1)變為線性回歸模型； 當λ=0時，式(2)變為線性回歸模型； 當ρ和γ同時為零時，式(3)變為線性回歸模型。采用極大似然法(maximum likelihood)估計3個SAR模型的參數，通過R統計語言spdep包的lagsarlm和errorsarlm函數實現。

1.4 空間加權矩陣

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

相接鄰近矩陣包括Rook相接矩陣、Queen相接矩陣、Bishop相接矩陣(Anselin, 2001)和Delaunay三角網(Delaunay triangulation, DT)矩陣(Smirnovetal., 2001)等。其中，Rook相接矩陣、Queen相接矩陣和Bishop相接矩陣一般用于規則的空間相接關系，DT矩陣可處理不規則的空間相接關系，因為林木間的空間關系一般是不規則的，因此本文采用DT矩陣。DT矩陣中，林木間為相接鄰近關系時，其空間加權值為1；否則為0。同時，利用R統計語言spdep包的trib2nb函數計算空間加權值。

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

(4)

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

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

(5)

(6)

(7)

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

利用R統計語言nlme包的corExp、corGaus和corSpher函數分別計算SV、GV和EV矩陣的空間加權值。

本文采用DT、ID1、ID2、ID5、SV、GV和EV共7個空間加權矩陣，應用于3個SAR模型中。

1.5 空間自相關檢驗指標

選擇常用的Moran’s I(MI)指數檢驗模型殘差的空間自相關，其表達式(Paradis, 2014)如下：

(8)

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

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

2 模型比較分析

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

2.1 基礎模型的確定

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

表3 3個備選基礎模型的擬合指標及參數估計①Tab. 3 Fitting index and parameter estimates for three alternative base models

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

表4 啞變量模型中3個主要樹種的參數估計Tab. 4 Parameter estimates for three main species in dummy variable model

2.2 全林分模型殘差空間自相關

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

此外，由表5還可知，SDM和SEM 2個SAR模型均能消除模型殘差的空間自相關矩陣中，依據MI值和P值，GV是最好的空間加權矩陣，因其有最大P值，但只對SEM有效，其次為EV和ID5，但ID5略遜于EV，最后為ID2和DT。

表5 4個模型殘差的空間自相關①Tab. 5 Test for spatial autocorrelation of four model residuals

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

2.3 全林分模型擬合

利用3個擬合指標R2、RMSE和AIC，比較分析BM與3個SAR模型(SLM、SDM和SEM)的擬合效果。由表6和表3可知，3個SAR模型的R2(0.578，0.669)均大于BM，RMSE和AIC均小于BM，說明3個SAR模型的擬合效果均優于BM。3個SAR模型中，SDM與SEM擬合效果不相上下，且均明顯優于SLM，這是因為SLM沒有消除模型殘差的空間自相關(詳見2.2)。另外，空間加權矩陣GV在SDM中的擬合效果明顯比SEM差，結合表5可知，這是因為GV在SDM中沒有消除模型殘差的空間自相關。因此，消除模型殘差的空間自相關，有助于改善模型擬合效果。

由表6還可知，依據R2、RMSE和AIC，在SDM和SEM中的空間加權矩陣(不考慮GV)從優至劣排序為： ID2> DT > ID > ID5 > EV。

表6 3個SAR模型的擬合Tab. 6 Three SAR models fitting

2.4 全林分模型參數

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

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

表7 3個SAR模型的回歸參數①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個SAR模型的自回歸參數Tab. 8 Autoregressive parameters of three SAR models

2.5 啞變量同步自回歸模型

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

表9 啞變量模型中3個主要樹種的參數估計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 結論

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

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

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

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

4) 所有7個空間加權矩陣中，SV和ID1應用于胸徑-樹高模型時，均為不合理的空間加權矩陣。

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

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

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

陳新美, 張會儒, 姜慧泉. 2010. 東北過伐林區蒙古櫟林空間結構分析與評價. 西南林學院學報, 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. 蒙古櫟林全林整體生長模型及其應用. 林業科學研究, 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])

李希菲, 唐守正, 袁國仁,等. 1994. 自動調控樹高曲線和一元立木材積模型. 林業科學研究, 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. 蒙古櫟林單木生長模型系統. 北京: 中國林業科學研究院碩士學位論文.

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

田 超, 楊新兵, 李 軍, 等. 2011. 冀北山地不同海拔蒙古櫟林枯落物和土壤水文效應. 水土保持學報, 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. 我國東北地區蒙古櫟林生物量及生產力的研究. 中國生態農業學報, 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. 黑龍江省不同地點蒙古櫟林生態特點研究. 生態學報, 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. 環境梯度下蒙古櫟群落的物種多樣性特征. 生態學報, 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.

(責任編輯 石紅青)

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。

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

S757.2

A

1001-7488(2017)06-0067-10

* 張會儒為通訊作者。