胸径和树高是森林调查及经营管理的2个重要因子(Lei et al., 2009)。相对于胸径，树高的测量比较困难，存在成本高、耗时长、观测误差大、视觉阻碍(如树梢顶部容易遮挡，特别是天然混交林)等不利因素(Chhetri et al., 1996; Colbert et al., 2002; Mehtätalo, 2004; Trincado et al., 2007; Krisnawati et al., 2010; Sharma et al., 2015)，因此，学者们经常利用林木胸径-树高模型估计树高，进而估计林木材积、立地质量、林分动态以及其他有关森林生长与收获(Lei et al., 2001; Peng, 2001; Calama et al., 2004; Lei et al., 2009)。

传统林木胸径-树高模型均假定林木间是相互独立的，通常采用普通最小二乘法(ordinary least squares, OLS)进行估计。然而，随着研究的逐步深入，学者们逐渐意识到林木间的关系并不是相互独立的(Frod et al., 1981; Tomppo, 1986; Mateu et al., 1998)，这是因为林木间的竞争和微环境结合作用，导致林木间存在空间自相关(spatial autocorrelaton)(Magnussen, 1990; Fox et al., 2001)。空间自相关是指一些变量在同一分布区内的观测数据之间潜在的相互依赖性，一般认为，相邻林木间的空间自相关程度强于相离林木间的空间自相关程度(Legendre, 1993; Legendre et al., 1998), 因此，忽略空间自相关将会使OLS回归违背残差独立分布假设，导致犯第一类错误(原假设为真，假设检验拒绝了原假设)的可能性变大(Legendre, 1993; Lennon, 2000; Legendre et al., 2002; Dormann et al., 2007; Kissling et al., 2008)、模型参数标准差的有偏估计及回归模型估计有效性降低(Krämer, 1980; West et al., 1984; Gregoire, 1987; Fox et al., 2001; Zhang et al., 2004b)等。

鉴于林木间存在空间自相关，目前已有学者将空间自相关模型加入到林木胸径-树高模型中，如地理加权回归(geographically weighted regression)模型(Zhang et al., 2004a; 2004b; 2005)、同步自回归(simultaneous autoregressive, SAR)模型(Meng et al., 2009; Lu et al., 2011)等，其中，SAR模型全面考虑了空间自相关过程发生在模型应变量、模型解释变量及模型误差中的可能性，是常用的空间自相关模型(Kissling et al., 2008; Meng et al., 2009; Lu et al., 2011)。具体来说，SAR模型包括空间滞后模型(spatial lag model, SLM)、空间误差模型(spatial error model, SEM)和空间Durbin模型(spatial Durbin model, SDM)，SLM假设空间自回归过程只发生在应变量中，SEM假设该过程只发生在模型误差项中，SDM假设该过程同时发生在应变量和解释变量中(Cressie, 1993; Haining, 2003)。

蒙古栎(Quercus mongolica)阔叶混交林是我国东北地区主要的天然次生林，蒙古栎主干坚硬耐腐、纹理美观，是良好的造材原料，枝条发热量高，可作为良好的薪炭材，叶可饲蚕和动物，根系发达且适应性强，具有涵养水源、防风固沙、减少径流、保持水土等作用，果实富含淀粉，经加工提炼可转化燃料乙醇及生物化工产品(田超等，2011；马武，2012)，对国民经济和生态环境均具有很高的价值(于顺利等，2001；2004；许中旗等，2006；陈新美等，2010；洪玲霞等，2012)。为此，本文以天然蒙古栎阔叶混交林为研究对象，考虑林木间的空间自相关，构建基于空间自相关的林木胸径-树高模型，以期为可持续经营天然混交林提供理论依据。

研究区位于吉林省汪清林业局塔子沟林场，地理坐标为129.971°—130.222°E，43.327°—43.492°N。属长白山系老爷岭山脉雪岭支脉，低山丘陵地貌，海拔300~1 200 m，坡度5°~25°。研究区属季风性气候，年均气温3.9 ℃左右，年均降水量600~700 mm，且集中在夏季，占全年总降水量的80%。土壤主要以暗棕壤为主，平均厚度40 cm左右。植被属长白山植物区系，群落结构复杂，植物种类较多，主要树种有蒙古栎、白桦(Betula platyphylla)、红松(Pinus koraiensis)、山杨(Populus davidiana)、榆(Ulmus pumila)、椴(Tilia tuan)、黑桦(Betula dahurica)、长白落叶松(Larix olgensis)、色木槭(Acer mono)、红皮云杉(Picea koraiensis)、冷杉(Abies fabri)、水曲柳(Fraxinus mandshurica)、黄檗(Phellodendron amurense)等。

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 固定样地林分特征^① Tab.1 Stand characteristics of permanent plots

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

由表 1可知，天然蒙古栎阔叶混交林中有4个主要树种，即蒙古栎、白桦、红松和山杨。考虑到山杨和白桦具有相近的林木生长特性，本文将山杨和白桦称为杨桦，作为同一建模树种，并建立全林分、蒙古栎、杨桦和红松4个胸径-树高模型。为不改变林木间的空间关系即林木间的空间自相关结构，采用哑变量方法即以树种为哑变量构建分树种的胸径-树高模型(简称哑变量模型)。

当前SAR模型均为线性的(Meng et al., 2009)，因此本文只考虑线性化的基础林木胸径-树高模型。根据图 1所示的胸径与树高散点图，选用3个备选线性化基础模型，具体见表 2。

图 1 不同树种的胸径-树高散点示意 Fig.1 Scatter plot of diameter at breast height against tree height for different species

表 2 备选基础模型^① Tab.2 Alternative base models

表 2 备选基础模型① Tab.2 Alternative base models

SAR模型包括空间滞后模型(SLM)、空间误差模型(SEM)和空间Durbin模型(SDM)，其表达式(Cliff et al., 1981; Anselin, 1988; Haining, 2003; Fortin et al., 2005; Dormann et al., 2007)如下：

(1)

(2)

(3)

式中：Y表示应变量矩阵；X表示解释变量矩阵；W表示空间加权矩阵；β表示回归系数；ε表示模型残差；ρ、λ和γ表示空间自回归系数；ξ表示空间相关残差。

WY、WX和Wξ为3个空间滞后变量。当ρ=0时，式(1) 变为线性回归模型；当λ=0时，式(2) 变为线性回归模型；当ρ和γ同时为零时，式(3) 变为线性回归模型。采用极大似然法(maximum likelihood)估计3个SAR模型的参数，通过R统计语言spdep包的lagsarlm和errorsarlm函数实现。

空间加权矩阵(spatial weight matrix)是SAR模型的重要组成部分，本文选取常用的空间加权矩阵，即相接邻近(contiguous neighbors)矩阵、逆距离幂(inverse distances raised to some power)矩阵和地统计(geostatistical)矩阵(Getis et al., 2004; Lu et al., 2011)。

相接邻近矩阵包括Rook相接矩阵、Queen相接矩阵、Bishop相接矩阵(Anselin, 2001)和Delaunay三角网(Delaunay triangulation, DT)矩阵(Smirnov et al., 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种形式(Getis et al., 2004)，其表达式如下：

(4)

式中：d_ij表示林木i与其相邻木j的距离，w_ij表示林木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表示变程，w_ij和d_ij意义同上。

利用R统计语言nlme包的corExp、corGaus和corSpher函数分别计算SV、GV和EV矩阵的空间加权值。

本文采用DT、ID1、ID2、ID5、SV、GV和EV共7个空间加权矩阵，应用于3个SAR模型中。

选择常用的Moran’s Ⅰ(MI)指数检验模型残差的空间自相关，其表达式(Paradis, 2014)如下：

(8)

式中:I表示模型残差的MI; n表示林木株数; y_i和y_j分别表示林木i和林木j的模型残差值；y表示平均模型残差；w_ij表示林木i和林木j的空间加权值。

其原假设为I等于期望值I₀=-1/(n-1)，表明不存在空间自相关。如果I显著(α=0.01) 大于I₀，则模型残差为正自相关，说明相邻林木间的相似程度高；反之如果I显著(α=0.01) 小于I₀，则模型残差为负自相关，说明相邻林木间的不相似程度高。利用R统计语言spdep包的moran.test函数进行检验。

为便于比较模型之间的拟合程度，本文选用决定系数(coefficient of determination, R²)、均方根误差(root mean square error, RMSE)和Akaike信息准则(Akaike information criterion, AIC)3个拟合指标比较分析模型的拟合效果。

从3个备选模型中选择全林分模型和哑变量模型建模效果均最优的线性化基础模型作为全林分和哑变量的基础模型(base model, BM)。值得注意的是，全林分模型和哑变量模型中林木间的空间关系是一致的，即林木间具有一致的空间自相关结构，因此，只需建立基于全林分基础模型的SLM、SEM和SDM共3个SAR模型，并通过比较分析模型残差空间自相关、模型拟合及模型参数估计效果，即可找到最适宜的SAR模型。最后，利用最适宜的SAR模型构建哑变量SAR模型。

利用3个拟合指标R²、RMSE和AIC，分别比较全林分和哑变量模型的3个备选基础模型。由表 3可知，全林分模型的备选基础模型C具有最大的R²、最小的RMSE和AIC；哑变量模型的备选模型A具有最大的R²、最小的RMSE和AIC，模型C略低于模型A，最差为模型B。由表 4可知，哑变量模型的备选模型A中主要树种蒙古栎和杨桦的参数β₀和β₁均显著(P＜0.01)，主要树种红松的参数β₀显著(P＜0.01)，而参数β₀不显著(P＞0.1)。因此，本文选用模型C作为全林分和哑变量的基础胸径-树高模型(BM)。

表 3 3个备选基础模型的拟合指标及参数估计^① Tab.3 Fitting index and parameter estimates for three alternative base models

表 3 3个备选基础模型的拟合指标及参数估计① Tab.3 Fitting index and parameter estimates for three alternative base models

表 4 哑变量模型中3个主要树种的参数估计 Tab.4 Parameter estimates for three main species in dummy variable model

表 4 哑变量模型中3个主要树种的参数估计 Tab.4 Parameter estimates for three main species in dummy variable model

由表 5可知，空间加权矩阵SV的BM残差MI值大于1，说明SV是不合理的空间加权矩阵，因此，以下结果分析中将不再考虑SV。除了SV外，其他6个空间加权矩阵的BM残差MI值均显著(P＜0.01) 大于期望值I₀=-1/(752-1)=-0.001 3，说明BM残差存在正空间自相关。6个空间加权矩阵的SLM残差MI值均显著(P＜0.01) 大于I₀，说明SLM没有消除模型残差的空间自相关，但是降低了空间自相关影响，因其MI值小于相对应的BM残差MI值。除了GV和ID1外，其他4个空间加权矩阵的SDM残差MI值均与I₀差异不显著(P＞0.1)，说明选择合适的空间加权矩阵，SDM可以消除模型残差的空间自相关，同时也表明产生模型残差的空间自相关可能是由滞后变量WY(即相邻加权树高)和滞后变量WX(即相邻加权胸径)共同作用产生的。除了ID1外，其他5个空间加权矩阵的SEM残差MI值均与I₀差异不显著(P＞0.1)，说明SEM与SDM类似，在合适的空间加权矩阵条件下，SEM可以消除模型残差的空间自相关，但空间自相关可能是由滞后变量Wξ(即相邻加权未观测的变量)产生的，因为SEM的前提假设是空间自相关过程来源于模型误差。

表 5 4个模型残差的空间自相关^① Tab.5 Test for spatial autocorrelation of four model residuals

表 5 4个模型残差的空间自相关① Tab.5 Test for spatial autocorrelation of four model residuals

此外，由表 5还可知，SDM和SEM 2个SAR模型均能消除模型残差的空间自相关矩阵中，依据MI值和P值，GV是最好的空间加权矩阵，因其有最大P值，但只对SEM有效，其次为EV和ID5，但ID5略逊于EV，最后为ID2和DT。

利用3个拟合指标R²、RMSE和AIC，比较分析BM与3个SAR模型(SLM、SDM和SEM)的拟合效果。由表 6和表 3可知，3个SAR模型的R²(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 3个SAR模型的拟合 Tab.6 Three SAR models fitting

表 6 3个SAR模型的拟合 Tab.6 Three SAR models fitting

由表 6还可知，依据R²、RMSE和AIC，在SDM和SEM中的空间加权矩阵(不考虑GV)从优至劣排序为：ID2＞DT＞ID＞ID5＞EV。

对3个SAR模型及BM的回归参数进行T检验，对3个SAR模型的自回归参数进行似然比检验(likelihood ratio test)。由表 7可知，无论采用哪个空间加权矩阵，3个SAR模型的回归参数β₁均与BM中的β₁相似，且均显著不为零(P＜2.2e-16)。相比β₁，SEM与BM中的β₀相似，但SDM和SLM中的β₀与BM中的β₀不相似，并且随着空间加权矩阵的变化而变化，这主要是由于SDM和SLM的模型形式中增加了空间滞后变量WY和WX。所有的β₀中，ID1应用于SDM中的β₀与零差异不显著(P＞0.1)，说明ID1在SDM中是不理想的空间加权矩阵。同时由表 7可知，利用ID1得出的自回归参数ρ、γ和λ均明显高于其他空间加权矩阵得出的值，且无论应用于哪个SAR模型，这进一步说明了ID1是一个不理想的空间加权矩阵。

表 7 3个SAR模型的回归参数^① Tab.7 Regressive parameters of three SAR models

表 7 3个SAR模型的回归参数① Tab.7 Regressive parameters of three SAR models

另外，由表 8可知，GV只有在SEM中才能使自回归参数λ显著(P＜0.01)。除了GV外，利用其他5个空间加权矩阵得出的ρ、γ和λ均显著(P＜0.01)，且均大于零，说明空间滞后变量WY、WX及Wξ均能产生正的空间自相关。

表 8 3个SAR模型的自回归参数 Tab.8 Autoregressive parameters of three SAR models

表 8 3个SAR模型的自回归参数 Tab.8 Autoregressive parameters of three SAR models

由基于全林分基础模型的SLM、SDM和SEM的残差MI检验效果、建模拟合效果和参数检验效果可知，以ID2为空间加权矩阵的SEM模型是清除基础模型残差空间自相关现象和提高模型精度最佳的同步自回归模型，因此，本文采用该模型建立哑变量SAR模型，得到其R²、RMSE和AIC分别为0.766、0.237和10.728，3个拟合指标均优于哑变量基础模型(表 3)，其自回归参数λ为0.412，且显著不为零(P＜0.01)。哑变量SAR模型中的蒙古栎、杨桦和红松的参数估计结果见表 9。由表 9可知，3个主要树种的截距参数β₀和β₁均显著不为零(P＜0.01)，参数检验效果理想。以ID2和SEM建立的全林分SAR模型残差、哑变量SAR模型残差以及3个主要树种的模型残差见图 2。

表 9 哑变量模型中3个主要树种的参数估计 Tab.9 Parameter estimates for three main species in dummy variable model

表 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.

本文以线性化林木胸径-树高模型为基础模型(BM)，加入3个SAR模型(SDM、SEM和SLM)构建了基于空间自相关的天然蒙古栎阔叶混交林林木胸径-树高模型，得到以下结论：

1) SDM和SEM能很好地处理空间自相关，而SLM最差，这一结论与Kissling等(2008)、Meng等(2009)和Lu等(2011)研究结果类似。

2) 无论采用哪个空间加权矩阵，SLM均不能消除模型残差空间自相关，但可降低空间自相关，在一定程度上提高了模型拟合效果。

3) 模型残差MI值均与I₀差异不显著，SDM和SLM虽能消除模型残差空间自相关，提高模型拟合效果，但需要选择合适的空间加权矩阵，如可选择ID2、DT、ID5和EV，其中ID2最好，能将全林分基础模型的R²从0.574提高到0.669，哑变量基础模型的R²从0.732提高到0.766。

4) 所有7个空间加权矩阵中，SV和ID1应用于胸径-树高模型时，均为不合理的空间加权矩阵。

5) 蒙古栎阔叶混交林林木胸径-树高模型残差的空间自相关可能由滞后变量WY(即相邻加权树高)和WX(即相邻加权胸径)共同产生或者可能由滞后变量Wξ(即相邻加权未观测的变量)产生。

6) 利用ID2和SEM构建了哑变量模型，从而得出蒙古栎、杨桦和红松的同步自回归胸径-树高模型。

本文研究了基于空间自相关的天然蒙古栎阔叶混交林林木胸径-树高模型，得出ID2和SEM是最适宜的空间自相关模型。对于其他林分类型，该模型是否最适宜是下一步需要探讨的问题。另外，本文只构建了基于线性化空间自相关模型的天然蒙古栎阔叶混交林林木胸径-树高模型，如何研究非线性化空间自相关模型就目前来说是一个难题，也是今后的重点研究方向。