Height-diameter model of Pinus massoniana and Schima superba mixed forest in central Guizhou Province

研究区位于贵州省贵阳市开阳县，该县西连安顺市与毕节市，东与南连接黔南布依族苗族自治州，北接遵义。位置为26°11′~27°22′N，106°07′~107°16′E，地势南北高、中间低、平均海拔为1 200 m，地形以山地、丘陵为主、喀斯特地貌分布广泛，占全区面积的71.8%。气候类型为亚热带高原季风湿润性气候，年均气温为15.3 ℃，年均降水量为1 196.9 mm。所选研究区森林类型为天然马尾松-木荷混交林，林内树种有马尾松、木荷、水青冈Fagus longipetiolata、光叶海桐Pittosporum glabratum、云贵鹅耳枥Carpinus pubescens和杉木Cunninghamia lanceolata等。

2017年8月，在研究区马尾松-木荷混交林内，根据CTFS(Centre for Tropical Forest Science)样地建设标准方法设置40 m×50 m的固定样地，并将样地划分为20个10 m×10 m的小样方，对每个样方内所有胸径≥1 cm的树木进行每木检尺，即记录树种名称，测量胸径、树高、冠幅、坐标等因子。样地内林木调查因子描述性统计见表1。

根据调查数据(表1)可知林内主要树种为马尾松和木荷。因此，将树种作为哑变量引入OLS、GAM、LMM、GWR和GWRK模型中，分别构建马尾松，木荷和其他树种树高-胸径模型。树种分类用定性代码0和1表示，第i个树种表示为S_i，定性数据S_i转换为(0, 1)形式；S_i为0或1，当x是第i个树种时为1，否者为0。其中i＝1、2、3；S₁、S₂、S₃分别是马尾松、木荷和其他树种的定性代码。将没有引入哑变量的模型作为全林木基础模型。

采用幂函数(Y=aX^b)的对数变换形式(lnY=a+blnX)作为本研究的基础模型。其模型表达式如下：

式(1)中，H_t为树高(m)，D_BH为胸径(cm)，a与b为模型参数，$ \varepsilon $为模型误差。

增加样方水平的单水平随机效应，在模型参数a和b的位置上分别增加随机参数$ {\mu }_{1}{\mathrm{和}\mu }_{2} $，如式(2)所示：

式(2)中，a和b为模型参数，$ {\mu }_{1} $和$ {\mu }_{2} $为随机参数，$ \varepsilon $为模型随机误差。

GAM是广义线性模型的半参数拓展，唯一的基本假设是平滑函数是可加的^[25]。将lnD_BH作为解释变量，lnH_t作为响应变量，平滑回归项采用平滑样条函数构建GAM，其表达式：

式(3)中，s为平滑样条函数，$ \varepsilon $为模型误差，$ \alpha $为截距。

GWR模型作为传统回归方法的拓展，将样本点的地理位置信息加入建模，使得模型既能描述响应变量与解释变量间的关系，又能消除模型残差的空间自相关^[26]。使用Gauss函数作为权函数，构建开阳县马尾松-木荷混交林树高-胸径模型。

式(4)中，$ \mathrm{l}\mathrm{n}{H}_{t} $为响应变量，$ \left({u}_{i},{v}_{i}\right) $为第i个样本点的坐标，β₀为第i个样本点的常数估计值，X_ki为第k个自变量在第i样本点的值，$ {\beta }_{k}\left({u}_{i},{v}_{i}\right) $为第i个样本点上的第k个回归参数，是关于地理位置的函数，p为用于建模的林木总数，$ \varepsilon $为模型误差。

GWRK模型是GWR模型的延伸与拓展，将GWR与克里格插值相结合的模型，通过GWR回归残差进行空间插值，然后将插值结果和GWR回归估计值相加，从而获得GWRK估计值^[27]。

分别采用全局与局域莫兰指数(Moran’s I)评价模型残差的空间分布^[9]。全局Moran’s I是用来描述所有的空间单元在整个区域上与周边地区的平均关联程度，取值为−1~1，多数统计检验在开始时都先确定零假设，在计算全局Moran’s I时通常会给出Z值，如果Z＞1.96或者＜−1.96，拒绝零假设，表示空间分布可能是聚集模式或分散模式；若−1.96≤Z≤1.96，接受零假设，表示空间分布可能是随机过程产生的结果，Z能够判断n个样点内的观测值之间是否存在空间相关性^[28]。而局域Moran’s I可以进一步研究局域空间聚集现象，评价模型残差是否存在局域聚集情况，分析模型残差局域空间自相关性。

空间异质性作为生态系统的主要属性之一，是产生空间格局的主要原因^[29]。使用模型残差组内方差评价模型残差空间异质性^{[28, 30]}。

模型评价指标选用决定系数(R²)、均方误差(MSE)和赤池信息准则(AIC)进行模型比较，除了决定系数值越大越好外，其他指标均是越小越好。综合以上评价指标选择最优模型。

使用R语言的minpack.lm包、mgcv包、lme4包和GWmodel包分别对OLS、GAM、LMM和GWR模型进行建模分析，随机选择581株(80%)林木作为建模数据，剩余的146株(20%)林木作为检验数据。用ArcMap 10.6对模型残差进行普通克里格插值，ROOKCASE (用VBA编写的Excel宏文件)计算全局和局域Moran’s I，作图在Excel中完成。

全林木基础模型拟合结果见表2，5种模型的R²为0.63~0.80，MSE为0.08~0.16，AIC为462.77~759.23，都能很好地反映树高与胸径之间的关系。将树种作为哑变量引入后，模型拟合精度整体高于全林木基础模型，5种回归模型的R²提高了0.01~0.06，MSE和AIC也有不同程度的降低，全林木基础模型和哑变量模型的3种拟合指标从劣到优均依次为OLS、GAM、LMM、GWR、GWRK。OLS和LMM哑变量模型参数估计值均在0.001水平上显著，GWR模型的参数估计值为一个范围，能很好地表现林木间的差异(表3)。

由表2可知：OLS和GAM模型残差全局Moran’s I在α=0.05的显著水平下，Z＞1.96，存在显著的空间自相关性。而LMM、GWR和GWRK模型残差具有更低的全局Moran’s I，且在α=0.05条件下，模型残差全局Moran’s I的Z在−1.96~1.96范围，模型残差空间自相关性不显著。以5 m为间隔，计算5~20 m滞后距离内模型残差全局Moran’s I (图1)。OLS和GAM模型残差全局Moran’s I随着滞后距离增加不断降低，而LMM、GWR和GWRK模型残差全局Moran’s I在各滞后距离处在0上下小幅度波动，且值相接近。

Figure 1.  Spatial correlation and local distribution of Moran’s I values of the five models’ residuals

OLS和GAM模型残差的局域Moran’s I整体分布不均匀且存在大而多的黑色气泡(正局域Moran’s I)(图1B~F)，说明这2种模型对热点的预测不足或过度。相比之下，LMM、GWR和GWRK模型残差产生更少的热点，局域Moran’s I多为负值，相邻点之间表现为相反特征。

图2为5~20 m滞后距离内各模型残差的组内方差。在该范围内，5种模型残差组内方差在5 m处最小，并随着滞后距离的增加而增加。GWR和GWRK模型残差组内方差在各滞后距离处均显著小于OLS、GAM和LMM模型残差组内方差。5种模型残差在各滞后距离处的组内方差从小到大排序为GWRK、GWR、LMM、GAM、OLS。

Figure 2.  Intra-block variance of residuals of five models

在全林木基础模型中，模型拟合精度由高到低排序为GWRK、GWR、LMM、GAM、OLS。OLS模型是对数据进行全局拟合，不考虑样本点之间空间自相关性的全局模型，即所得到的模型参数估计与数据的地理位置无关，估计结果是该点的最优无偏估计的同时也是研究区内所有样本点的最优无偏估计^[28]。GAM模型作为数据驱动模型，使用链接函数建立响应变量和解释变量的平滑函数之间的关系，为响应变量提供比OLS模型更好的预测^[31]。LMM模型拟合效果同样优于OLS模型，本研究通过增加样方水平的随机效应，将参数a和b同时作为随机效应参数，发现将样方作为随机效应可以消除样方差异对模型的影响^[32]，在拟合过程中展现出更好的效果。但在实际研究中，不同地理位置之间具有空间非平稳性，OLS、GAM和LMM模型属于全局模型，没有考虑样本点位置不同会造成不同参数预测值的问题，导致模型拟合效果不佳。GWR和GWRK模型作为一种加入空间权重函数的局域空间回归模型，考虑了不同地理位置之间的差异性，反映了局域参数的变化，比OLS模型拟合精度更高。GWRK模型对GWR模型残差进行普通克里格插值后，消除了残差空间相关对模型拟合的影响^[33]，模型拟合精度比GWR显著提高。将树种作为哑变量引入模型后，模型拟合精度比全林木基础模型均有提升，原因在于哑变量模型通常能更加精确地反映树木生长受到本身遗传性质的影响，即将不同树种存在的差异考虑进模型。娄明华等^[8]使用树种作为哑变量构建混交林林木树高-胸径模型时，同样发现哑变量模型拟合精度高于基础模型。有研究发现：不同树种所拟合的树高-胸径模型曲线不同^[20]，本研究构建的马尾松、木荷以及其他树种树高-胸径模型具有不同模型参数预测值，同样也说明了遗传特征是混交林树木生长的决定因素之一。整体而言，局域模型拟合效果优于全局模型，哑变量模型拟合精度高于全林木基础模型。

考虑模型残差空间自相关时，LMM、GWR和GWRK模型比OLS和GAM模型得到显著改进。5种模型残差全局Moran’s I均表现在滞后距离为20 m时趋近于0，这表明树高之间的空间自相关性随着距离的增加逐渐减小，且当距离超过20 m时，空间相关性微弱且不显著。

OLS模型前提假设认为林木之间是相互独立的，然而大量研究发现林木间的关系并不是相互独立的^[34]，因此忽略空间自相关会使OLS违背独立性假设导致对模型参数的标准误差的偏见估计^[9]，导致犯第一类错误(原假设为真，假设检验拒绝了原假设)的可能性变大^[8]。与OLS模型相比，GAM模型残差全局Moran’s I降低，但空间自相关性仍表现显著(Z＞1.96，α=0.05)，表明GAM并没有显著降低残差空间自相关，因为GAM并没有将空间自相关纳入建模过程，尽管GAM模型由于其稳健性和灵活性而改进了模型拟合并产生了更好的预测，但是该模型本质上还是非空间的^[31]。同时，GAM和OLS模型残差局域Moran’s I存在较多聚集的热点，残差之间存在显著空间正相关性，即相邻点之间的表现为相同特征，造成对热点有较多过度预测或预测不足的值。然而，在α=0.05条件下，LMM、GWR和GWRK模型残差全局Moran’s I随着滞后距离的增加越来越接近于0，且−1.96≤Z≤1.96，残差空间自相关性不显著。LMM、GWR和GWRK模型残差局域Moran’s I在研究区内整体表现为较多冷点，研究区相邻点残差呈现相反的特征。LMM模型中协方差结构可以反映个体间的差异，类似GWR模型中的权函数，因此LMM模型在提高模型拟合效果的同时还能消除模型残差空间自相关性的影响^[35]。GWR和GWRK模型对空间数据具有较强的分析能力，在建模过程中解释变量的回归系数随着位置的变化而变化，因此能够有效地探测空间数据的空间非平稳性^[36]。张凌宇等^[37]在研究大兴安岭中部天然次生林更新分布时，比较全局模型和局域模型残差空间自相关时研究发现：局域模型残差自相关明显小于全局模型，这与本研究结果相同。陈科屹等^[38]利用GWR模型研究进界木空间分布，同样得到类似的结论。

5种回归模型在5~20 m滞后距离范围内，模型残差组内方差随着距离的增加逐渐增加。与OLS模型相比，LMM和GAM模型的组内方差在各滞后距离内变化不大，而GWR和GWRK模型却存在很大程度的降低，说明全局模型不能降低模型残差空间异质性，而局域模型不仅能消除模型残差空间自相关，还在降低空间异质性方面表现出明显的优势。ZHANG等^[31]对加拿大3种不同林分建立OLS、LMM和GWR树高-胸径模型，研究模型残差空间异质性时，得出结果与本研究结论相同。综上，OLS和GAM模型在模型拟合精度、降低空间自相关和空间异质性方面表现不佳，GWR和GWRK模型不仅能提高模型拟合精度，还能降低模型残差空间自相关和异质性，建议用于林木树高-胸径模型的建立。

本研究以开阳县马尾松-木荷混交林为研究对象，以OLS、GAM、LMM、GWR和GWRK模型为全林木树高-胸径基础模型，将树种作为哑变量引入基础模型，并比较全林木基础模型和哑变量模型的拟合效果，同时分析哑变量模型降低空间自相关和异质性程度，发现哑变量模型拟合效果整体优于基础模型，局域模型较全局模型不仅能显著提高模型拟合效果，还能降低模型残差空间自相关性和异质性。今后研究马尾松-木荷混交林树高-胸径模型时，建议应用局域模型得到精度更高的拟合结果，若考虑区域尺度问题，可以将混交比、立地质量等因子作为哑变量引入模型，进而构建更大尺度范围的混交林树高-胸径模型。