森林生态系统作为最大的陆地碳库, 其植被和土壤分别约占全球植被碳库的86%以上和土壤碳库的73%以上, 在减缓气候变化和维护全球碳循环等方面发挥着重要的作用^[1-2]。因此, 森林生态系统的固碳能力受到广泛关注。碳密度是衡量森林生态系统固碳能力的重要指标之一^[3-4], 也是准确估算森林碳储量的基础。由于森林生态系统是一个具有较强的空间异质性的复杂系统, 其碳密度受气候、地形地貌及周围其他环境等因素的影响^[5], 而这些因素在空间上也往往存在相关性, 在空间上相互影响、相互作用。因此, 从区域尺度上探明影响森林碳密度的主要因素及构建估测模型, 可为森林碳储量估算及碳汇林业的经营管理等提供科学依据。

在区域尺度上, 森林碳密度/碳储量的估测主要是通过构建碳密度/碳储量与其影响因子之间的关系模型。由于传统的回归模型在构建中忽略了各变量之间的空间异质性, 因此, 空间误差模型^[6]、空间滞后模型^[7]、空间滤波模型^[8]、空间Durbin模型^[9]、地理加权回归模型^[10]等空间模型受到许多学者的关注。如李月等^[11]、刘畅等^[12]采用空间误差模型构建碳储量与影响因素关系模型时均得出其拟合精度明显优于最小二乘模型；Lin等^[13]、欧光龙等^[14]、郭含茹等^[15]基于地理加权回归模型进行研究, 均表明地理加权回归模型比最小二乘模型在碳密度/碳储量拟合及预估等方面具有明显优势。然而, 由于森林生态系统的自然地理要素及森林类型等的不同可能会造成不同区域相同森林生态系统或同一区域不同森林生态系统碳密度空间异质性的差异, 从而可能导致其碳密度估测最优空间模型的不同, 而目前多数学者主要集中在某一空间模型与传统回归模型在碳密度估测中的比较分析, 较为缺乏对不同空间模型之间的比较研究, 加强这方面的研究将有利于更加全面了解和掌握空间模型在森林生态系统碳密度估测方面的应用。

马尾松(Pinus massoniana)是我国南方分布广泛的乡土针叶树种, 具有适应性强、耐干旱与贫瘠的特点。江西省马尾松林总面积为239.4万hm², 占乔木林总面积的28.9%, 是江西省主要的森林类型之一, 在区域碳循环及应对气候变化中扮演着重要角色^[16]。本研究以江西省马尾松林生态系统为研究对象, 分别利用最小二乘模型、空间误差模型、空间滞后模型和地理加权回归模型构建其碳密度与影响因子之间的关系模型, 比较分析选出最优拟合模型, 以丰富空间模型在森林碳密度估测中的应用及为江西省马尾松林碳汇林业的经营规划与管理提供参考依据。

1 材料与方法 1.1 研究区概况

江西省位于我国东南部, 地理坐标为24°29′—30°04′N, 113°34′—118°28′E, 面积16.69×10⁴ km²。该区属于亚热带温暖湿润季风气候, 气候温和, 雨量充沛, 年均气温16.3—19.5℃, 年均降雨量1675 mm。地貌以山地、丘陵为主, 林地的土壤类型主要有红壤、黄红壤、黄壤等, 成土母岩主要有花岗岩、页岩、砂岩、板岩等。全省森林资源丰富, 森林覆盖率为63.1%, 主要森林类型有常绿阔叶林、针叶林、针阔叶混交林、竹林等。

1.2 数据与方法 1.2.1 数据来源

样地数据。根据研究区的森林分布图及森林资源统计表, 综合考虑马尾松林在全省的地域分布、林龄结构、地形地貌特征等因素, 于2011—2016年期间设置样地调查, 并收集2016年江西省森林资源连续清查马尾松林样地数据(对部分缺失林下植被、凋落物等相关数据的样地进行补充调查)。样地面积为900 m²和800 m², 样地数共计611个, 其样地分布见图 1。乔木层每木调查起测胸径为≥5 cm(胸径＜5 cm的乔木视为灌木), 测定其胸径、树高、林分密度等林分因子及地理坐标、海拔、坡度、坡向、坡位、土层厚度、腐殖质层厚度等立地因子, 同时选取标准木取样测定其干、枝、叶、根的碳含量。在样地的上、中、下部分别布设2 m×2 m的灌木样方, 在所选灌木样方中各设置1个1 m×1 m的草本样方和凋落物样方。采用“样方收获法”分别收集灌木(分枝、叶、根)、草本(分地上、地下部分)及凋落物并称量鲜重, 取样用于测定其含水率和碳含量。在样地内选择代表性地块挖取100 cm深度(未达100 cm的挖至母岩)的土壤剖面, 根据土层厚度按0—10 cm、10—20 cm、20—30 cm、30—50 cm、50—100 cm五个层次, 对各层土壤采用环刀法(体积为100 cm³)取样来测定土壤容重, 并分别取约1 kg的土壤样品带回实验室用于土壤碳含量的测定。

气象数据。本研究首先根据各样地的地理坐标从Wang等^[17]编写的提取亚太地区气候数据的ClimateAP软件中获取1980—2016年期间各样地每年的年降水量与年温度数据, 然后计算得出各样地的年均降水量与年均温度。

1.2.2 碳含量测定与碳密度计算

所有样品研磨粉碎后过0.25 mm筛, 采用重铬酸钾氧化-外加热法来测定其碳含量^[18]。乔木层、林下植被层及凋落物层碳密度为各组分生物量乘以对应的碳含量, 其中乔木层生物量采用异速生长方程W=a×D^b, 应用我国立木生物量模型及碳计量参数系列行业标准计算^[19-22]。土壤层碳密度的具体计算公式如下^[16]：

(1)

式中, D_soc为土壤层碳密度(t/hm²), C_i为土壤有机碳含量(g/kg), D_i为土壤容重(g/cm³), E_i为土层厚度(cm), G_i为大于2 mm的石砾所占的体积百分比(%)。

乔木层、林下植被层、凋落物层与土壤层的碳密度之和为生态系统碳密度。

1.2.3 模型选取与评价

模型选取。(1)最小二乘模型(OLS)：是建立在整体区域上的一种因变量偏差的平方和最小的模型估计；(2)空间误差模型(SEM)：是指模型的误差项导致了空间变量之间的相关性, 变量之间的空间相互作用存在于误差项, 研究的主要是相邻样地点的误差冲击对其他样地点的影响；(3)空间滞后模型(SLM)：是指被解释变量之间的空间依赖性对模型非常重要从而导致了空间相关, 它研究的是因变量在邻接地区的行为对整个区域其他地区行为的影响；(4)地理加权回归模型(GWR)：是一种局部参数估计的模型, 其本质是为数据集中的每一个要素都建立独立的方程^[23], 通过将空间结构嵌入线性回归模型中, 以探测空间关系的非平稳性, 该方法的估计结果有明确的解析表示, 而且还能对其结果的参数估计进行统计检验, 其计算公式为：

(2)

其中, (y_i; x_i1, x_i2, …, x_ip)为在地理位置(u_i, v_i)处的因变量y_i和自变量x_i1, x_i2, …, x_ip的观测值(i= 1, 2, …, n)。β_k(u_i, v_i)(k=1, 2, …, p)为第i个观测点(u_i, v_i)处的未知参数, 是(u_i, v_i)的任意函数, ε_i(i= 1, 2, …, n)为独立同分布的误差项, 通常假定其服从N(0, σ²)。

模型评价。选用决定系数(R²)、均方误差(MSE)和赤池信息量准则(AIC)3个统计量对模型拟合效果进行评价。R²主要是用来评价模型的拟合优度的指标, 其值越高则表明模型的拟合效果越好, 即拟合出来的模型稳定性越好；MSE主要是用来衡量因变量实测值与模型预估值之间偏差的指标, 其值越低表明模型的预估结果越接近于实际测量值, 即模型的拟合能力越高；AIC可以用来衡量模型的实用性和复杂性, 其值越小表明模型的拟合效果越好。

模型残差检验。选择Moran′s Ⅰ指数检验模型残差的空间自相关性, 其计算公式如下：

(3)

式中, x_i和x_j分别表示研究区域内样点i和j上的观测值；x表示研究区域内所有样点观测值的平均值；n表示研究区域内的样点总数；w_ij(d)表示空间权重, 即样点i和j之间的距离小于选定尺度值时权重为1, 超过尺度值时权重为0。

在R软件中, 采用多元线性逐步回归方法对生态系统碳密度的影响因子进行筛选, 并构建最小二乘模型；利用Geoda软件建立空间误差模型和空间滞后模型；在GWR 4.0软件中建立地理加权回归模型。

2 结果与分析 2.1 影响因子筛选

根据森林碳密度研究的相关文献, 采用频度统计法及样地调查内容选取影响生态系统碳密度的因子, 主要包括立地因子(海拔、坡度、坡向、坡位、土层厚度、腐殖质厚度)、植被因子(林龄、胸径、株数密度、树高、郁闭度、林下植被盖度、凋落物厚度)以及气象因子(年均温度、年均降水量)等15个因子。对坡向、坡位定性因子进行赋值, 将坡向划分为阳坡、半阳坡、半阴坡、阴坡、无坡向, 并依次赋值为1—5；将坡位划分为山脊、上部、中部、下部、山谷、平地、全坡, 依次赋值为1—7。利用多元线性逐步回归方法对各因子与生态系统碳密度进行分析, 从中筛选出与碳密度相关性显著且各自变量之间的方差膨胀因子均小于5(即不存在多重共线性)的因子, 用于数据建模。筛选出的因子分别为：海拔、坡度、土层厚度、胸径、年均温度、年均降水量, 其统计量见表 1。

2.2 模型构建 2.2.1 最小二乘模型、空间误差模型与空间滞后模型

最小二乘模型、空间误差模型与空间滞后模型的参数估计结果见表 2, 坡度与碳密度的系数为负值, 这表明碳密度随着坡度的增大呈减少的趋势；海拔、土层厚度、胸径、年均温度、年均降水量与碳密度的系数均为正值, 表明碳密度与这些因子之间存在正相关, 其中系数最大的为胸径, 最小的为年均降水量。海拔、坡度、土层厚度和胸径对碳密度的影响均达显著性水平(P < 0.05)。

2.2.2 地理加权回归模型

地理加权回归模型的参数估计结果见表 3, 海拔、土层厚度、胸径与碳密度的系数平均值分别为0.0817、0.4616、5.0531, 其第1四分位数至最大值的范围分别为0.0074—0.2920、0.3081—1.1967、3.5514—10.3628, 可看出这些因子75%的系数均为正值, 说明碳密度与海拔、土层厚度、胸径3个因子之间存在正相关关系；坡度与碳密度的系数平均值为-2.4527, 其最小值至第3四分位数的范围为-4.4625至-1.8330, 表明坡度与碳密度之间为负相关关系, 即碳密度随着坡度的增大呈逐渐减少的趋势。年均温度、年均降水量与碳密度的系数平均值为1.2671和0.0237, 但50%的系数为正值, 另外50%部分为负值, 说明气象因子对碳密度的影响比较多变。

从模型参数估计结果可知, 4种模型均表明碳密度与坡度呈负相关, 与海拔、土层厚度、胸径呈正相关。最小二乘模型、空间误差模型和空间滞后模型的结果均表明海拔、坡度、土层厚度、胸径对碳密度有显著影响(P < 0.05), 年均温度和年均降水量对碳密度的影响并不显著(P>0.05), 而地理加权回归模型的结果则反映出气象因子对碳密度的影响比较多变。

2.3 模型评价与检验 2.3.1 模型评价

据表 4可得, 4种模型的决定系数(R²)最大的为地理加权回归模型, 其次为空间误差模型与空间滞后模型, 最小的为最小二乘模型；地理加权回归模型的均方误差(MSE)与赤池信息准则(AIC)值均为4种模型中最小, 而最小二乘模型的AIC与MSE值为最大。由此可得, 模型拟合效果由高到低依次为地理加权回归模型>空间误差模型>空间滞后模型>最小二乘模型。

图 2为碳密度实测值与4种模型估计值之间的散点图。4种模型的估计值与碳密度实测值的配对点均离散分布在1:1的参考线附近, 当碳密度实测值较小时, 最小二乘模型、空间误差模型及空间滞后模型拟合得出的估计值普遍偏高, 而地理加权回归模型的估计值比其他3种模型更为接近实测值；当碳密度实测值较大时, 最小二乘模型、空间误差模型及空间滞后模型拟合得出的估计值均普遍偏低, 但地理加权回归模型的估计值与实测值偏低的程度明显小于其他3种模型。由此可得, 地理加权回归模型的拟合结果优于其他3种模型。

2.3.2 模型残差的空间自相关性

由表 5可得, 最小二乘模型、空间误差模型和空间滞后模型的Moran′s Ⅰ均为正值, 且P < 0.05；地理加权回归模型的Moran′s Ⅰ为负值, 且P>0.05。由此得出, 最小二乘模型、空间误差模型和空间滞后模型的模型残差均存在显著的空间自相关性, 而地理加权回归模型的模型残差的空间自相关性不显著, 这表明地理加权回归模型可有效降低模型残差的空间自相关性。

4种模型残差的Moran′s Ⅰ随距离的变化见图 3。最小二乘模型、空间误差模型和空间滞后模型残差的Moran′s Ⅰ均为正值, 即3种模型残差的空间相关性为正相关, 且模型残差的Moran′s Ⅰ在带宽较小时, 其降低幅度均较大, 随着距离的增加其降低幅度逐渐减弱, 逐渐趋近于0；地理加权回归模型残差的Moran′s Ⅰ为负值, 即其模型残差的空间相关性为负相关, 其模型残差的Moran′s Ⅰ在带宽较小时, 其增加幅度较大, 随着距离的增大其增加幅度逐渐减弱, 也逐渐趋近于0。模型残差的Moran′s Ⅰ由大到小表现为最小二乘模型>空间滞后模型>空间误差模型>地理加权回归模型。由此也可得出, 与其他模型相比, 地理加权回归模型有效降低了模型残差的空间自相关性。

综上分析, 4种模型中, 地理加权回归模型决定系数(R²)最大, 均方误差(MSE)与赤池信息准则(AIC)均最小, 且能有效降低模型残差的空间自相关性, 表明地理加权回归模型的拟合效果优于其他3种模型。

3 讨论与结论 3.1 讨论

由于森林生态系统具有较强的空间异质性, 所以, 不少学者利用空间模型研究森林碳密度/碳储量/生物量时均得出空间模型比最小二乘模型(OLS)的拟合效果更好的结论。如李月等^[11]采用空间误差模型(SEM)构建碳储量与影响因素关系时得出其拟合精度明显优于OLS模型, 其主要是由于SEM减小了模型残差的空间自相关性；刘畅等^[24]利用空间滞后模型(SLM)预估林木含碳量得出SLM优于OLS模型, 认为SLM模型能更好地解决模型构建中的空间效应, 其稳定性更高, 对各因子的解释能力更强。

在空间模型的运用中, 更多学者采用了地理加权回归模型(GWR), 且其研究结果几乎都表明GWR优于其他模型。如Liu等^[25]利用OLS模型、线性混合模型和GWR模型拟合黑龙江省乔木层碳储量与影响因子间的关系时表明了GWR模型在拟合精度以及在预估方面都明显优于其他两种模型；Lin等^[13]利用OLS模型和GWR模型拟合福建省将乐县乔木层碳密度与影响因子之间的关系, 结果也表明了GWR模型在拟合精度及降低模型残差自相关性方面都明显优于OLS模型。此外, 欧光龙等^[14]、Zhang等^[26]的研究均表明GWR模型为最优模型。而本研究通过对OLS、SEM、SLM和GWR 4种模型的模型评价指标及残差检验的综合分析, 得出OLS模型的拟合效果最差, GWR模型最好。这可能主要是因为在区域尺度上, 森林碳密度本身存在明显空间效应, 而OLS模型是以忽略数据的空间效应作为前提下的无偏估计；SEM、SLM模型与OLS模型相比, 虽然减小了模型残差的空间自相关性, 提高了拟合效果, 但由于没有直接将这种空间自相关性纳入到拟合过程中, 难于解决空间异质性问题^[27]；而GWR模型适用于空间非平稳性的研究^[28], 它在对数据集所有位置点的参数进行估计时考虑了其空间非平稳性, 针对每一个坐标位置点都有相对应的参数^[29], 将模型残差的空间自相关性直接纳入到拟合过程中, 有效降低了模型残差的空间自相关性, 所以其拟合效果最优。

另外, 一些学者以森林碳储量/生物量为因变量, 以遥感影像数据为自变量构建空间模型时同样也得出了GWR模型优于其他空间模型的结论。如戚玉娇^[30]通过构建黑龙江省大兴安岭地区的森林植被地上碳储量与TM遥感影像数据之间的关系表明GWR模型的拟合效果最优；郭含茹等^[15]利用GWR模型、传统全局回归模型、协同克里格插值构建浙江省仙台县森林地上部分碳储量和Landsat TM遥感影像数据之间的关系也表明GWR模型是有效的, 闾妍宇等^[31]、Wang等^[32]的研究也得出同样的结论。可见, GWR模型能有效解决森林生态系统空间异质性问题。

在区域尺度上, 对森林碳密度/碳储量/生物量空间模型构建时, 多数学者仅采用某一空间模型与传统回归模型进行比较, 而本研究通过不同空间模型比较分析, 发现GWR模型拟合效果优于其他空间模型。然而, 基于GWR模型的扩展模型的相应出现, 如地理加权回归Kriging模型、地理加权回归Logistic模型、地理加权回归泊松模型等, 而这些扩展模型若用于本研究中其拟合精度是否会提高, 还有待今后进一步深入探讨。

3.2 结论

本研究以江西省为研究区, 通过多元线性逐步回归方法筛选出与马尾松林生态系统碳密度相关性显著且不存在多重共线性的影响因子为海拔、坡度、土层厚度、胸径、年均温度和年均降水量。采用OLS、SEM、SLM和GWR模型构建生态系统碳密度与其影响因子之间的关系模型, 通过对模型拟合效果的评价及残差检验的综合分析, 得出GWR模型拟合效果最优, 更适用于江西省马尾松林生态系统碳密度的估测。