土壤作为陆地生态系统中最大的碳库^[1-2], 在全球碳循环中占据重要作用^[3-5]。土壤中的碳素主要以有机碳的形式存在^[6], 土壤有机碳含量及其储量的变化都会对生态系统的稳定及功能发挥产生重要影响。土壤有机碳含量受土地利用方式^[7-8]、管理措施^[9-10]、群落结构^[11-12]、地形地势^[13-14]等诸多因素的综合影响, 在空间上具有明显的差异性^{[13, 15]}。土壤有机碳含量的微小变化会对生态系统和全球气候变化产生重要影响^[16-17]。对土壤有机碳空间分布特征及影响因素的研究有助于认清土壤有机碳变化过程及区域化反应。同时, 能够加强对土壤碳库的科学认识, 进一步了解土壤碳库"源"和"汇"的关系, 在土壤有机碳的提升和土地资源管理方面有重大意义。

传统方法采用空间地统计分析对土壤有机碳空间分布进行研究, 该方法只考虑样点之间距离的空间关联性, 而忽视了各影响因素的贡献率^[18-19]。探究土壤有机碳影响因素的研究多采用最小二乘法(Ordinary least squares，OLS)和逐步回归等方法^[20-21]。这些线性回归模型都是基于全局回归模型, 来反映区域总体状况, 其认为各要素对空间上所有区域的影响是相同的, 忽略了各影响因子的局部性, 导致模拟结果与实际情况相差较大^[22-23]。地理加权回归模型(Geographically weighted regression, GWR)是一种基于样点地理位置的局域空间分析方法, 能够将局部范围的因变量和解释变量进行合并, 从而得到每个样点各解释变量的回归系数。与传统回归模型相比, 该模型在研究土壤有机碳的空间建模中具有更高的精度。因此该方法在研究土壤有机碳分布及与环境因素、土壤性质的关系方面得到了广泛应用^{[13, 24-25]}。受研究区范围影响, 在空间大尺度上侧重环境因素和土壤有机碳的关系^[13], 在小尺度上侧重土壤性质对有机碳空间特征的影响^[26], 综合考虑环境因素和土壤理化因子对有机碳的共同作用的成果较少。

流域是社会-经济-自然要素综合作用的复合体, 其组成要素的多样性决定了流域的复杂性。在流域尺度上探讨各因素对土壤有机碳的影响, 是对开展区域生态学问题研究、维持各生态系统良性发展、综合考虑各类要素间利用和协调发展问题进行的有益探索。本文以伊河流域土壤为例, 利用GWR模型探索土壤有机碳与环境因子和其它土壤性质之间的关系, 揭示伊河流域土壤有机碳的空间分布特征及影响因素, 以期为伊河流域土地利用和管理提供一定的科学依据。

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

伊河流域位于河南省西部山区(33°39′—34°41′N、111°19′—112°54′E), 地处我国二、三阶梯的过渡地带。伊河发源于熊耳山, 流域总面积约6100 km^{2 [27]}, 海拔88—2128 m。地势西南高, 东北低, 地貌类型丰富, 主要包括山地、丘陵和平原。该区域处于北亚热带向暖温带过渡区, 属季风性气候, 年均温12.4—15.1℃, 年均降水量700—900 mm^[28], 四季变化明显。伊河流域雨热同期, 降水多集中于夏季, 且降雨强度较大, 土壤易受侵蚀；流域内上游降水较为充足, 其年降水量约是下游的2倍。伊河流域土壤类型较多, 分布最广的是褐土, 同时还存在棕壤、红黏土、沙土等。该流域植被类型以暖温带落叶阔叶林为主, 在海拔较高地区存在针叶、落叶阔叶混交林, 平原区和低山丘陵区是人类活动的主要场所, 大面积的低山丘陵被开发为农田, 自然植被分布较少。

伊河流域总人口237.2万, 城镇化率42.68%, 人口集中分布在河流沿岸, 尤其是中下游的丘陵和平原区^[28]。伊河流域土地利用类型以耕地和林地为主(图 2), 城镇和村落等建设用地面积约325.10 km², 主要分布在河流沿岸的丘陵和平原区；耕地面积约2458.87 km², 主要分布在丘陵和平原区, 在低山区也有少量分布；林地面积约2714.57 km², 主要分布在低山和中山区, 在丘陵和平原区有零星分布；草地面积约360.82 km², 主要分布在低山区, 在其它地区有零星分布。

1.2 样品采集和指标测定

根据研究区面积和交通可达性, 对研究区采用6 km×6 km网格取样, 采集表层(0—20 cm)土壤样品141个, 其中农田样品76个, 森林样品59个, 草地样品6个。将采集的土样带回室内自然风干, 取适量土样研磨过筛(2 mm、0.25 mm、0.15 mm)备用。土壤理化性质测定方法参照《土壤农业化学分析方法》^[29], 具体方法如下(表 1)：

1.3 研究方法 1.3.1 环境变量提取

本文以伊河流域DEM数据为基础, 通过ArcGIS 10.2空间分析模块对DEM数据进行地形指标的提取计算。共计算了11个地形因子, 包括：海拔、坡度、坡向、曲率、平面曲率、剖面曲率、地形起伏度、地表粗糙度、复合地形指数(CTI)、汇流动力指数(SPI)和沉积物运移指数(STI)。其中CTI、SPI和STI分别由公式(1)、(2)、和(3)计算获得。研究区DEM来源于91卫图的GoogleEarth高程数据, 空间分辨率为15 m。选取的2个气象因子(年均温、年均降雨量)通过逐日气象数据计算获得, 数据来源于中国气象局(http://data.cma.cn/)。另外选取与土壤样品采集同时期的归一化植被指数(NDVI), 来反映研究区植被覆盖情况, 数据来源于美国NASA网站(https://www.nasa.gov/), 空间分辨率为250 m。将上述14个因子作为探讨影响土壤有机碳的环境因子。

(1)

(2)

(3)

式中, AC为垂直于水流方向的特定汇流面积, β为坡度。

1.3.2 GWR模型

最小二乘法模型(OLS)是随机变量(y)与确定性变量(x₁, x₂、x_i)的多元线性函数, 这是一种基于全局回归的函数(式4)。GWR模型^[30]是对OLS模型的拓展, 是一种局部回归模型, 将数据的地理位置嵌入到回归参数之中, 可以实现对参数的局部估计(式5)。

(4)

(5)

式中, y_i为样点i的因变量, x_ik为第i个点上第k个变量的观测值, (μ_i, ν_i)为第i个点的位置坐标, β₀(μ_i, ν_i)为截距, β_k(μ_i, ν_i)为第i个的回归系数, ε_i为误差项。

GWR模型中参数的估算关键在于空间权重函数的选取, 本文比较4种常用空间权重函数, 以确定适合本区域的最优空间函数, 分别是固定高斯函数(Fixed Gaussian)(式6)、自适应高斯函数(Adaptive Gaussian)(式7)、固定截尾型函数(Fixed bi-square)(式8)和自适应截尾型函数(Adaptive bi-square)(式9)。

(6)

(7)

(8)

(9)

式中, ij代表给定研究区域的任意一点, θ>0代表窗宽或者光滑参数。

带宽的选择是影响GWR模型分析结果的关键因素, 带宽的大小会直接影响回归参数的估计。本文采用黄金分割搜索来选择带宽, 以AICc作为带宽选择准则。

(10)

(11)

其中, K是参数数量, L是似然函数。

1.4 数据处理

实验数据在Excel中记录整理；地形因子通过ArcGIS 10.2提取、计算；GWR模型在GWR 4.0中运行；文中空间分布图均采用WGS_1984地理坐标系；相关分析、多元线性回归和共线性诊断在R软件中实现, 相关分析采用斯皮尔曼(Spearman)法(P < 0.05)。

2 结果与分析 2.1 伊河流域土壤有机碳的描述性统计特征及与环境变量和土壤性质的相关性

伊河流域表层土壤有机碳范围在3.37—38.34 g/kg, 平均含量为12.23 g/kg, 变异系数为0.47, 属于中等变异水平(表 2)。流域上、中、下游有机碳平均含量分别为12.43 g/kg、12.48 g/kg和11.72 g/kg, 变异系数分别为0.59、0.46、0.29, 均属于中等变异水平, 变异系数自上游到下游逐渐降低, 表明伊河流域土壤有机碳存在空间分布差异, 其中上游差异最大, 下游差异最小。

将土壤有机碳与环境因子和其它土壤性质进行相关分析(图 3)。环境因子中的年均气温与土壤有机碳呈显著负相关(P < 0.05), 其它环境因子与有机碳的相关性较弱；土壤因子中的容重、铵态氮与土壤有机碳表现出显著负相关(P < 0.01)的趋势；50—250 μm和250—1000 μm粒级含量、活性有机碳、全氮、硝态氮和全磷与有机碳呈正相关, 且相关性显著(P < 0.01)；pH和速效磷与有机碳相关性不显著。

2.2 GWR模型构建 2.2.1 模型变量选择

相关性作为模型解释变量选择的初步条件, 通过相关分析, 初步筛选出年均气温、活性有机碳、全氮、硝态氮、铵态氮、容重、全磷、50—250 μm和250—1000 μm粒级含量9个指标作为土壤有机碳的解释变量。为减少解释变量存在共线性引起GWR预测产生偏差, 需对解释变量进行共线性诊断, 并筛选剔除存在共线性的变量。本文选取方差因子(Variance inflation factor, VIF)、特征值、条件索引和方差比例来判定各解释变量之间是否存在共线性。由表 3发现, 各解释变量VIF值均小于10, 在所有主成分中, 只有9和10的条件索引大于30, 且每一主成分各解释变量的方差比例只存在单一变量的方差比例大于0.5, 加之最后一个主成分不能提取出来信息, 特征值接近0, 条件索引相应会很大。综上所述, 年均气温、活性有机碳、全氮、硝态氮、铵态氮、容重、全磷、50—250 μm和250—1000 μm粒级含量均为土壤有机碳的解释变量。

2.2.2 模型诊断与选择

以土壤有机碳为因变量, 年均气温、活性有机碳、全氮、硝态氮、铵态氮、容重、全磷、50—250 μm和250—1000 μm粒级含量为自变量, 选择不同的空间权函数进行建模。在模型运行时需对模型的拟合效果进行诊断和选择, 选择最优模型。诊断参数包括残差平方和、AICc、R²和调整后R^2[31]。不同权函数GWR模型残差平方和和AICc值均小于OLS, R²和调整后R²均高于OLS(表 4), 表明对于本研究区来说, GWR模型优于OLS模型。高斯空间权函数的残差平方和、AICc值均高于截尾型函数, R²和调整后R²均低于截尾型函数。自适应截尾型空间权函数的残差平方和和AICc值最低, 分别为1945.76和825.20, R²与调整后R²最大, 分别为0.58和0.49, 自适应截尾型空间权函数是本研究区的最优空间权函数, 其模型精准度最高。

2.3 GWR模型土壤有机碳解释因子系数空间分布

GWR模型预测伊河流域土壤有机碳空间分布整体结果较好, 局部决定系数在0.49—0.64之间, 自下游到上游, 决定系数逐步提高, 对上游的预测精度更高(图 4)。土壤有机碳各解释因子回归系数具有一定差异。回归系数最大的是AOC, 范围在1.28—2.98, 表明土壤有机碳受AOC影响最大, 容重的回归系数绝对值最小, 在-0.45—-0.01之间, 表明有机碳受容重影响较小。各解释因子空间分布具有局部性和非均匀性, 表明同一因子对不同地点有机碳的影响程度不同。年平均气温对土壤有机碳的影响整体均为负效应, 在伊河流域年平均气温越高, 土壤有机碳含量越低, 其系数绝对值中、上游最大, 下游最小, 年平均气温对中、上游有机碳限制性影响最大, 下游影响微弱。全氮对有机碳的影响为正效应, 全氮含量越高, 有机碳含量越高, 全氮系数自上游到下游逐步增加, 下游有机碳受全氮影响最大, 上游全氮系数在0.173—0.607, 上游有机碳受全氮影响较弱。活性有机碳对有机碳是正效应, 系数自上游到下游逐步减小, 活性有机碳对上游有机碳的影响最大, 中游次之, 下游最小。硝态氮的系数在0.18—1.61之间, 自下游到上游逐步增加, 上游系数较小, 硝态氮对上游有机碳影响较弱, 对下游影响较大。铵态氮系数在-1.13—0.58之间, 铵态氮对有机碳是一个负效应, 在中游地区最为明显, 上游和下游系数较小, 影响较弱。容重整体系数较小, 对有机碳影响较弱, 也存在空间差异, 自上游到下游, 对有机碳的负效应逐步加大。全磷对有机碳的影响在空间上有较大差异, 系数范围在-0.18—1.29之间, 在下游地区, 全磷系数为负值, 全磷对有机碳影响为负效应, 中游和上游其系数为正值, 对有机碳影响为正效应, 在上游地区全磷系数较大, 上游全磷对有机碳的正效应影响最大。土壤颗粒组成对土壤有机碳的影响在空间上分异较大。其中50—250 μm粒级系数在中上游为负值, 表明对有机碳的影响为负效应, 但绝对值较小, 说明在对中上游有机碳影响较弱, 其系数下游为正值, 表明对下游有机碳的影响为正效应。250—1000 μm粒级系数范围在0.83—1.68, 下游系数较小, 中游最大, 表明250—1000 μm粒级对下游有机碳影响较小, 对中游和上游影响较大。

2.4 土壤有机碳空间分布特征

GWR模型预测伊河流域土壤有机碳含量在4.31—26.93 g/kg之间(图 5), 土壤有机碳含量的最大值和最小值都在上游地区, 且低值区占的面积较大, 说明上游地区有机碳含量差别较大, 空间变异性最大。下游地区有机碳含量较接近, 大多在10.95—22.03 g/kg之间, 空间变异性较小。中游地区高值区、中值区和低值区范围接近, 空间变异性居中。总体来说, 土壤有机碳空间分布的高值区大多集中在海拔较高的中山区和地势平坦的平原区, 低值区主要分布在低山、丘陵和河流沿岸。

3 讨论

环境因子影响土壤有机碳的空间分布^{[13-14, 24, 32]}, 研究发现, 地形因子、NDVI和气象因子中的年均降水量与有机碳相关性不显著, 与罗梅等^[13]研究年均降水量、归一化植被指数、高程以及地形粗糙指数与有机碳含量呈极显著相关关系有差别, 可能是因为研究区人类活动剧烈, 从而干扰了环境因子的作用, 加之研究区范围较小, 环境因子作用不明显。从土壤有机碳的空间分布来看, 与地形有较为一致的规律, 有机碳分布的高值区大多在海拔较高地区。高海拔地区气温较低, 植被覆盖以灌木和草本为主, 大量枯枝落叶堆积, 对降水的缓冲截留, 避免了土壤有机碳的流失。但从相关分析来看, 有机碳的空间分布与海拔关系较弱, 这方面的原因可能是由于采样造成的, 海拔较高的地方难以到达, 采样点较少, 从而影响了海拔与有机碳的相关关系。研究发现, 年平均温度对土壤有机碳是负相关, 这与吴春生等^[33]研究一致, 气温降低, 土壤有机碳释放速率降低, 碳氮矿化速率减慢, 同时环境微生物活性降低, 分解动植物残体速度减慢, 致使土壤有机碳和全氮含量积累量显著^[34-35]。

自然因素和人类活动综合影响土壤有机碳的空间分布特征, 且具有明显的地域性。在海拔较高的中上游区域, 土壤有机碳含量主要受年平均气温、活性有机碳和全磷含量的影响。活性有机碳主要来源于土壤有机碳的水解, 植物根系的分泌物、及其凋落物的分解, 和土壤微生物本身及其代谢产物^[36-37]；而全磷主要来源于成土母质^[38]。这表明在自然状态下, 土壤有机碳主要受立地环境、成土母质和地表覆盖的影响。在中上游的低山丘陵区, 土壤有机碳主要受容重、250 μm以上粒级含量和速效氮的影响。低山丘陵区大多被改造为耕地, 受季风性气候的影响, 降水集中, 多暴雨, 土壤侵蚀严重, 细颗粒物、速效氮易被淋溶、冲刷, 是引起土壤中碳含量减少的原因, 加之低山丘陵区耕作交通不便, 化肥和有机肥难以大量施用, 收获的秸秆主要被用来当作燃料, 使得有机碳归还较少, 也就导致该区有机碳含量较低。这表明在低山丘陵区, 人类活动和环境因素, 尤其是在地形和气候的共同作用影响了土壤有机碳的含量。平原区有机碳含量较高, 主要受全氮和硝态氮的影响, 可能是因为平原区主要是耕地, 地势平坦, 交通便利, 大量化肥和有机肥的投入, 从而使得土壤有机碳含量较高^[9]。

GWR模型预测土壤有机碳可以恰当地反映出因子在局部区域的影响, 并且其预测效果在自然状态下或人为干扰较弱的地区更好。通过对比四种空间权重矩阵函数结果, GWR模型在空间权重矩阵函数的选择上, 自适应型函数可以根据采样点之间的距离调整权重, 更具有合理性^[39]。该模型在土壤方面应用的最大的问题之一就是解释变量的局部共线性^[40], 土壤性质间相互作用, 有较强的相关性^[41], 很容易产生局部共线性问题, 这是现阶段GWR模型在土壤方面应该被重视的问题之一, 也是今后该模型需要解决的主要问题。

4 结论

伊河流域表层土壤有机碳范围在3.37—38.34 g/kg, 流域上、中、下游有机碳平均含量分别为12.43 g/kg、12.48 g/kg和11.72 g/kg, 伊河流域土壤有机碳存在空间分布差异, 其中上游差异最大, 下游差异最小。相关分析表明, 土壤有机碳含量受年平均气温以外的环境因子影响较小, 与土壤理化性质相关性显著。

GWR模型预测土壤有机碳可以恰当地反映出各因素在局部区域的影响。结果表明, 在海拔较高的中上游区域, 土壤有机碳含量较高, 主要受立地环境、成土母质和地表覆盖的影响。在中上游的低山丘陵区, 人类活动和环境因素共同影响下土壤有机碳含量较低。中下游平原区农业活动和化肥投入是影响土壤有机碳含量较高的原因。这为伊河流域各亚区生态系统的合理发展和管理提供了依据。

GWR模型预测伊河流域土壤有机碳空间分布整体结果较好, 局部决定系数在0.49—0.64之间, 自下游到上游, 决定系数逐步提高, 对上游的预测精度最高。该模型能够较好预测伊河流域土壤有机碳的空间分布特征, 对人为干扰较弱的区域预测效果更好。