鄂西犟河流域表层土壤有机碳时空变化及其影响因素

符雅盛, 张利华, 郑芸蔷, 朱志儒, 马永明, 崔 越

〔中国地质大学(武汉) 地理与信息工程学院地理系, 湖北 武汉 430074)〕

土壤是陆地生态系统中最大的碳库[1]，是陆地植被碳库(5.00×1011～6.00×1011t)的2～3倍，是全球大气碳库(7.50×1011t)的2倍多[2]。土壤有机碳(soil organic carbon， SOC)包括植物、动物及微生物的遗体、排泄物、分泌物及其部分分解产物和土壤腐殖质。土壤有机碳的储量是植物残体进入土壤并受土壤微生物作用分解损失后的动态平衡结果[3]。土壤碳是全球碳储量的重要组成，其固定及释出，不仅影响到陆地植被的养分供给，也影响了陆地生态系统的组成、结构、分布及功能，更对温室气体及全球气候变化具有重要影响[4]。其库容大小受到气候、植被、土壤理化性质以及人类活动等诸多自然和人为因素的影响。深入研究SOC在小型丘陵山地集水区的分异规律及其影响因素对小流域碳循环、碳分配以及防止森林SOC流失、土壤肥力退化具有重要意义。

中外学者在土壤有机碳储量、密度、空间分异、时间变化及其影响因子等方面做了大量研究，对全球土壤碳库和碳循环研究具有借鉴意义。全球主要生态系统及土壤剖面的土壤有机碳库研究已较成熟，不同模型和方法的计算结果无数量级差异。区域插值模型、回归拟合模型和机器学习等模型被较多的运用于土壤有机碳含量的估算[5]。空间变异性、依赖性以及各环境因子对土壤有机碳的影响可以被用来认识和研究其含量[6]。当前，随着研究程度的加深和研究对象的细化，中小尺度的土壤有机碳空间分异和动态变化及其影响因素的研究尚存在一些不足。小尺度的土壤有机碳研究集中在单一土地覆盖类型下土壤有机碳含量与海拔、地形之间的关系[7]。首先，传统的有机碳空间分布预测局限于土壤普查等历史数据，数据陈旧、更新慢、时间分辨率低。全国土壤普查在小流域尺度的实测剖面不足，空间分辨率低。其次，目前的研究多以土地利用、生态类型和地貌类型为研究对象，研究较多的有耕作模式、植被类型以及生境对土壤有机碳含量的影响。对于地表复杂、山体破碎的小流域尺度土壤SOC的影响因子的研究还不够充分。此外，人类活动的规模不断扩大，地表植被的非自然演替以及土地利用使得传统方法不再具有普适性。小流域尺度的研究区还有着地表形态复杂、人为活动强度大、土壤变异性强等特点，地形的复杂弱化了空间自相关，地统计方法的效果被削弱。

本文基于野外工作、实验室测定、遥感影像解译和地形因子计算，研究环境变量对土壤有机碳含量的影响；通过地统计分析研究区土壤有机碳变异特征，引入植被、地形等因子进行拟合，比较其相关性并作主成分分析，确定研究区土壤有机碳含量时空变化的影响因子和主导因素，以期为当地和类似区域的水土保持提供科学建议。

1 研究区概况

犟河位于丹江口水库上游100 km处，是堵河的一级支流，发源于十堰市张湾区大独岭，东西流向，在十堰市黄龙镇东湾村附近汇入堵河，干流全长25 km左右。犟河流域跨越东经110°31′42″—110°43′24″，北纬32°29′51″—32°42′32″，流域面积326 km2，是南水北调中线工程的重要水源地之一。流域地处秦巴山区东延余脉，位于大巴山和武当山之间，地势南高北低，自西南向东北倾斜，间有河谷平地，海拔集中在150 m到1 100 m之间。流域地层以中元古界武当山群为主，岩性为浅变质岩系，以石英片岩为主。流域土壤质地以壤质砂土为主，流域气候属亚热带季风性气候，多年平均气温为15.4 ℃，多年平均降水量为769.6 mm，年总平均降水日数为115～120 d。流域属于北亚热带和温带过渡区的自然生态系统，森林覆盖率近80%，既具有保存完整的原生林，也具有次生经济林和涵养水源林，植被类型以常绿针叶林、落叶阔叶灌木和针阔混交林为主。

2 材料与方法

2.1 数据来源与处理

基于ArcGIS的水文分析模块对十堰市ALOS World 3 D-DSM(空间分辨率为30 m，下载于http:∥www.eorc.jaxa.jp/ALOS/en/aw3 d30/)进行水文分析，获得犟河流域的流域边界，以犟河流域的ALOS World 3 D-DSM和2017年3月的Landsat8-OLI影像(下载于http:∥www.gscloud.cn/)为基础资料，提取犟河流域的主要水系(见封2，附图4)和土地覆盖类型(见封2，附图5)。使用ENVI的地形模块计算坡度、坡向、剖面曲率、横向曲率、纵向曲率和水平曲率。提取植被指数包括经典的归一化植被指数(NDVI)和比值植被指数(RVI)。将犟河干流及其支流分为上游、中游、下游3个部分进行土壤样品采集(见封2，附图4)：样品采集时间为2017年1，6月和10月以及2018年1月，共采集表层土壤样品100个；土壤采集样区的土地覆盖方式以针叶林、阔叶灌木和农田为主(见封2，附图5)，采集深度均为去除枯枝落叶层后10—15 cm，每个样区面积为5 m×5 m，均采集四角及对角线中心点的样品进行混合，然后采用四分法采集总量均为1 000 g左右的混合样品装入样袋并进行记录，回室内后自然风干研磨过筛(0.25 mm)备用。土样SOC含量由中国地质大学(武汉)地球科学学院的总有机碳分析仪(德国Elementar公司的vario TOC)进行分析和测定。

2.2 研究方法

用地统计的克里格模型预测研究区土壤有机碳含量的空间分布并作误差分析，用SPSS软件作影响因子的非线性拟合以及主成分分析。

3 结果与分析

3.1 描述性统计和地统计分析

3.1.1 土壤有机碳含量的描述性统计 土壤有机碳含量的描述性统计是一种能够反映土壤有机碳含量基本特征的统计学方法。土壤样品分为2017年冬(2017年1月)，2017年夏(2017年6月)，2017秋(2017年10月)，2018冬(2018年1月)共4期，次年春季的土样因实验室仪器损坏而损毁。根据SPSS描述性统计及正态检验结果(见表1)，土壤有机碳含量的变异系数大于水体有机碳，因为土壤质地、覆盖类型和地形因素等都对土壤有机碳有影响。总体而言，变异系数由大到小依次为：冬季、夏季、次年冬季、秋季。

表1 4期土壤有机碳含量描述性统计特征

Kolmogorov-Smirnov检验(K-S检验)是一个比较频率分布f(x)与理论分布g(x)或者两个观测值分布的检验方法。K-S检验是一种适用于小样本的非参数检验方法。根据K-S检验结果，有3期样本的显著性大于0.05，满足正态分布，为后续地统计分析提供参考。

3.1.2 土壤有机碳的地统计空间分布模拟 描述性统计可以描述数据的总体分布特征，但无法准确反映土壤有机碳的空间变异结构。地统计学基于参数估计和假设检验的思想，既考虑样本值的大小，又重视样本空间位置及样本间的距离，弥补了经典统计学忽略空间方位的缺陷。因此采用地统计学方法对研究区SOC含量的结构性与随机性进行定量描述。2018年1月样本分布最均匀，且样本数量较大，因此选择该期样本进行克里金插值，拟合全流域的土壤有机碳分布状况。

分别用高斯模型和指数模型对半变异函数点云进行拟合并制图，并比较两种模型生成的预测图的优劣，比较的统计指标见表2。对比表2中的数据，确定高斯模型作为最适模型。插值结果(见封2，附图6)示，犟河流域土壤有机碳含量的低值出现在阳南沟上游—鲍花沟中游—安沟下游—安沟一线以北以及财神沟中游地区；土壤有机碳含量的中等水平分布在大西沟下游、鲍花沟下游、犟河中下游低坡地区；土壤有机碳含量的高值分布在大西沟上游及财神沟上游以南的山地地区。由附图7可知：犟河中游以西—阳南沟鲍花沟中游以南—大西沟和财神沟中游以北地区预测误差小，预测误差最大处主要在大西沟和枧堰沟上游及以南山区。当模型非常适合用于生产和制图时，预测误差的平均值应趋于0，标准均方根趋于1，均方根以及平均标准误差应尽可能小。高斯模型的参数如下：块金值C0=36.762，偏基台C=0，变程=0.217，方向=9°。空间异质比为块金值/基台值，用于衡量随机因素和结构性变异对土壤空间分异性的影响。C0/(C+C0)=1>0.75，说明土壤有机碳的空间相关性极弱，即距离不是土壤有机碳含量发生空间变异的主导因子，其变异主要由结构性因素和随机性因素导致。

表2 克里金交叉验证参数

3.2 影响土壤有机碳的因素分析

影响土壤有机碳含量的结构性因素和随机因素主要有土壤质地、土地覆盖、植被等环境因子以及海拔、坡度、坡向、地表曲率等地形因子。

3.2.1 土地覆盖 根据由犟河流域采集得到的46个表层土壤样品的SOC含量数据以及流域土地利用数据，采用描述性统计分析方法得到如表3所示的统计结果。

表3 不同土地利用方式下土壤有机碳含量统计特征

犟河流域SOC含量为14.32±8.62 g/kg，不同土地利用下SOC含量差异不显著(p>0.01)。其中，农田SOC含量的均值最大，园地、混交林和针叶林次之，灌木最小。不同的土地利用方式会直接导致地表覆被的改变，植被覆盖越高，枯枝落叶堆积越多，有利于表层土壤有机碳的富集。土壤SOC的含量存在显著差异。园地受人为因素影响严重，生活污水的倾倒和排放也间接给作物提供肥力。作物不同，对农田进行施肥、播种、松土等活动的频次和强度不同，使得农田的土壤有机碳含量的差异也较大。由此可见，农田和园地受人为因素影响最大。比农田的土壤有机碳含量稍低的是混交林，混交林整体的土壤有机碳水平较针叶林高。混交林中的树种有杉树、竹、杨树、泡桐树等，由于生物群落复杂，生态系统稳定性高，有利于土壤有机碳的生成与积累[8]。

3.2.2 土壤类型 成土过程、发育程度及耕作管理措施的差异形成了不同的土质类型，不同土质的理化性质也存在不同。石灰性冲积土的土壤有机碳含量相对较高，简育高活性淋溶土和不饱和雏形土的土壤有机碳含量接近(见表4)。

表4 不同土壤类型的土壤有机碳含量统计特征

3.2.3 坡度 分析4期土壤样品有机碳含量与坡度的Pearson线性相关性，发现仅2017年夏季样本相关系数为0.618(p<0.05)，即土壤有机碳含量随坡度的增加而增加，其他几期并无明显规律。其中，线性模型为最优拟合模型，模型表达式为：Y=6.067+0.299X，表明其有机碳含量随坡度的增加而增加(见表5)。

表5 2017年夏季土壤有机碳含量与坡度回归拟合模型

3.2.4 土层深度 课题组于2018年在流域东部的头堰水库采集样本，样点位于人迹罕至的山坡坡顶，排除了人为干扰，植被覆盖为针叶阔叶混交林，坡度为42°。去掉枯枝落叶层后的土壤厚度为25 cm。以5 cm为一个单位厚度采集同一剖面的土样5个，测得其土壤有机碳的含量分别为5.22 g/kg(0—5 cm)，15.57 g/kg(5—10 cm)，38.51 g/kg(10—15 cm)，19.17 g/kg(15—20 cm)，8.93 g/kg(20—25 cm)。该剖面土壤有机碳含量的平均值为17.48 g/kg，且随土壤厚度的加深有先增后减的趋势。

3.2.5 坡位 在头堰水库的42°坡的坡中采集土样，覆盖为针阔叶混交林，其土壤有机碳含量为14.34 g/kg；在坡脚采集土样，土地覆盖为灌丛，测得其土壤有机碳平均值为13.99 g/kg。坡顶、坡中、坡底的土壤有机碳含量接近，有不明显的递减趋势。

3.2.6 植被指数 秋季是利用遥感影像光谱信息评估不同植被生长状况的最佳观测时期。经预处理后计算出归一化植被指数(NDVI)和比值植被指数(RVI)(见图1)。NDVI和RVI都是反映植被长势和营养信息的重要指标之一，前者在低植被区的敏感性强，在植被覆盖度高(50%)的地区敏感性降低，而RVI对高值被覆盖度下的植被差异响应明显。归一化植被指数图中亮色为低值，暗色为高值；图1中显示犟河流域除了犟河干流、流域靠近市区的东北部以及南北数条沟内为非植被，其余地区都覆有植被。比值植被指数图中颜色越暗，植被越密、长势越好。由图1可知RVI以犟河干流为轴向两侧逐渐增加，而南北两侧又沿着沟向两侧山坡逐渐增加。由于犟河南部的面积比犟河北侧部分的面积大，人类活动少，海拔差异也越大，植被表现出一定程度的差异：犟河南侧东北部头堰水库附近有较多人工种植的次生针叶林，财神沟海拔中游和上游海拔700 m以上为原生针叶林，中游海拔500～700 m以下为原生针叶林灌木混交林，下游与头堰水库一致，为次生针叶林。

图1 犟河流域植被指数分布

(1) NDVI值。秋季土壤样本与NDVI存在中等程度的正相关(相关系数为0.49，p<0.01)，即土壤有机碳含量随NDVI的增加而增加(见表6)。

表6 2017年秋季土壤有机碳含量与坡度回归拟合模型

(2) RVI值。秋季土壤样本与RVI存在强正相关(相关系数为0.494，p<0.05)，即土壤有机碳含量随RVI的增加而增加(见表7)。

表7 2017年秋季土壤有机碳含量与RVI回归拟合模型

3.2.7 地表曲率 地面曲率是度量地形表面单位点的扭曲变化程度的重要因子，是地面复杂程度以及地面切割破碎程度的体现。垂直方向上的分量称为剖面曲率，水平方向上的水平曲率。

水平曲率是等高线方向的变化率，即坡向的变化率，它是一个反映等高线弯曲程度的指标，反映了坡向的变化。剖面曲率是指地面上任意一点位地表坡度的变化率，即垂直方向上地形的凹凸变化。经过线性和非线性拟合，在2017年1月和6月的样本中均发现土壤有机碳含量随地表曲率绝对值的增大而增加的规律。使用线性、指数、对数等模型进行拟合分析，结果表明指数函数的拟合效果优于其他函数，最优拟合模型的参数见表8。而秋季和次年冬季的样本并无明显规律。2017年秋季发生持续性降水，犟河干流水位达到1986年以来的最高值，土壤有机碳流失加重；2018年冬季采样主要目的是分析不同土地利用方式对土壤有机碳的影响，因此土壤有机碳含量更多的受土地利用方式控制。

3.2.8 主成分分析 土壤有机碳含量的影响因子颇多，且各因子之间相互关联，相互作用。用“降维”思想对各影响因子做主成分分析(见表9)。主成分对应特征值大于1是主成分分析的原则。由表9可知，最终主成分的特征值分别为3.997，1.818，即原始变量的信息量分别被扩大了3.997，1.181倍，其中前两个主成分指标的贡献率分别达到66.228%和30.294%，累计能够解释所有原始变量信息的96.582%，且特征值均大于1，达到主成分分析的要求。

表8 犟河流域表层土壤有机碳含量与地表曲率回归拟合模型

表9 土壤有机碳含量影响因子主成分的特征值和贡献率

根据表10可知，第一主成分由剖面曲率、横向曲率、水平曲率和纵向曲率决定，它们在主因子上的载荷分别为0.982，-0.978，0.963，0.944，是犟河流域土壤有机碳的主导影响因子，而曲率是地形因子的反映；第二主成分由NDVI和RVI构成，它们在主因子上的载荷分别为0.926，0.904，是犟河流域土壤有机碳的次要解释变量，代表了植被因子。

表10 地形和植被因素主成分因子载荷

4 讨论与结论

4.1 讨 论

犟河流域属于鄂西典型的丘陵山地小流域，山体破碎、地形起伏大。流域内土壤样本有机碳的平均含量为11.39 g/kg，低于全国平均值(22.28 g/kg)；变幅为0.6～63.9 g/kg，变异系数为76.97%，低于全国水平(96.10%)，最大值是最小值的106.5倍，变异幅度较大，具有强空间异质性。块基比为1，表明随机部分引起的空间异质性占系统总变异的比例高，即样本间的变异更多的是由随机因素引起，而非距离，因此SOC的空间自相关极弱[9]。地统计结果表明土壤有机碳含量在水平空间分布上无明显趋势，建立在地理学第一定律基础上的地统计学方法在山体破碎、地表形态复杂地区的适用性不强。

此外，犟河流域属于鄂西北丹江口水库周边山地丘陵水质维护保土区，全区水土流失面积占区域国土面积的29.10%，研究区北部以中度侵蚀—强烈侵蚀—极强烈侵蚀为主，南部沟谷也有条带状轻度侵蚀—中度侵蚀，集中在城镇建设用地以及大西沟和枧堰沟内，中部犟河干流的侵蚀强度由东至西降低[10]。仅有7.77%的土样的有机碳含量高于全国平均值，土壤有机碳储存能力较弱。通常土壤中有机碳主要集中在土壤表层，并随土层深度的增加而减少。人为耕作和管理可以提高土壤有机碳含量，但不合理的土壤资源利用也会破坏土壤碳储量。犟河下游以石灰性冲积土为主，土壤类型样本均值和地统计插值都显示了高值。除了土壤类型因素，下游较多的农田和园地，即土地覆盖类型，亦是区域高值出现的原因。季节尺度的SOC变化体现了低温下微生物活动分解作用减弱、SOC得以保存，而温度上升和降水增加都能加重SOC的消耗。

相关分析表明：水平曲率、剖面曲率等地表曲率的绝对值以及RVI和NDVI两个植被指数都和土壤有机碳呈显著正相关关系。地表曲率能影响流水的侵蚀能力，而流水是地表重力作用的主要营力之一，地表曲率直接影响着降水对地表土壤的侵蚀强度，进而影响土壤有机碳的流失，但与此同时也促使泥沙中有机碳的富集[11]。中等坡度与凹面剖面曲率和剖面曲率的组合导致了更高的SOC含量[12]。鲍伟佳[13]在安徽省舒城县龙潭小流域的表层土壤有机碳研究中也发现了地表曲率与有机碳含量的相关性，龙潭小流域属于山地丘陵小流域，但海拔差和地形复杂程度均不如犟河流域，研究发现平面曲率和剖面曲率与SOC含量的相关性分别达到0.46，0.56(p<0.01)，即SOC含量与平面曲率和剖面曲率存在极显著的正相关，但SOC含量与NDVI和坡向的关系不显著；卢宏亮等[14]利用随机森林模型，预测了安徽省全省的土壤属性，其中以水平曲率为代表的地形因子在随机森林预测模型中的重要性位列第三，仅次于高程和NDVI，考虑安徽省山地和丘陵地貌类型的占比较大，其规律和本文结论接近。此外，在水文学和农学领域的研究和经验也能证明地表曲率在水土流失过程中的作用，例如凹形坡(收敛形山坡)的含水量较凸形坡多，通常凹陷地貌的土壤碳储量更高[15]。通过改造地表曲率，可以提升土壤的固碳能力，例如“顺坡沟垄”和“横坡沟垄”耕作方式可以减轻土壤有机碳的流失，且后者效果更好[16]。估算相邻研究区丹江口某村庄的土壤有机碳库的空间变化时，考虑建立在区域地形基础上发展的地形农业，可以提高预测的准确性和可比性，研究发现水道或支流的分水岭附近的SOC含量通常较高[17]。

因此，自然状态下地表形态复杂的山地丘陵地区的土壤有机碳主要受地表曲率影响，地表的凸处和凹处分别具有侵蚀强度弱和堆积能力强的储碳优势。植被指数代表了植被生长状况，而植被又是土壤有机碳主要的重要来源之一[18]，因此植被指数作为次要解释变量，可以评价相邻或地表形态相近研究区的土壤有机碳含量。

4.2 结 论

(1) 丘陵山地小流域山体破碎、地形起伏大，其土壤有机碳含量变异幅度较大、空间异质性强且空间自相关极弱。样本间的变异更多的是由随机因素引起，而非距离。建立在地理学第一定律基础上的地统计学方法对类似地貌的研究区的适用性不强。

(2) 研究区不同季节土壤有机碳的控制因子并不单一，夏秋两季土壤有机碳主要受地表曲率影响，意味着地表的凸处和凹处分别具有侵蚀强度弱和堆积能力强的储碳优势。植被指数作为次要解释变量，可以评价相邻或地表形态相近研究区的土壤有机碳含量。