土壤侵蚀是全球范围内最严重的土壤退化问题之一, 其发生发展变化过程直接影响到水土资源开发利用, 是威胁人类生存、社会经济发展的全球性问题^[1-2]。近几十年来有关土壤侵蚀的研究涵盖了地形、降水、土壤、植被、放牧强度和土地利用变化(Land Use/Cover Change, LUCC)等诸多因子在不同时空尺度上与其相互影响机制^[3-11]。在微观尺度上Vaezi等^[4]分析了半干旱地区降水对水土流失影响, 指出了雨滴动量冲击引起地土壤物理性质改变是造成水土流失的主要原因；在宏观尺度Du等^[10]对宁夏-内蒙地区的研究结果表明, LUCC是影响区域水土流失最重要因子之一。无论在哪一尺度, 共同的结论指向了植被因子具备减弱水土流失的力学基础和多尺度效应功能^[12]。从控制水土流失角度来看, 大量研究表明人工植被恢复过程是控制和减少土壤侵蚀的有效途径之一^{[6, 8, 11, 13]}, 尤其是草地植被恢复在控制土壤侵蚀的生态建设中起着至关重要的作用^[14-15]。从植被状况对水土流失影响研究的主要共识来看, 主要有以下几个特点：(1)植被盖度增加是减少水土流失的主要因素^{[8-9, 16]}；(2)植被类型对减流减沙具有综合效果, 相比较而言, 其在影响变化幅度上弱于盖度^{[6-7, 14]}；(3)植被可以通过覆盖坡面有效降低雨滴能量来减少冲刷作用^[17-18]；(4)植被覆盖通过冠层延缓降水增加土壤水分入渗量, 同时根系物理固结作用也增强了土壤抗侵蚀力^{[4, 12, 19-20]}。

这些研究结果对水土流失过程与植被的关系做了比较清晰的解译, 但是对植被组成、结构及盖度等数量化指标与土壤侵蚀量(模数)关系研究相对较弱。在这些指标中植被盖度是可用各种遥感手段快速获取的因子, 其具备的多源、多时像、多分辨率特点可满足不同空间范围研究需要。因此, 在特定研究区内植被盖度与土壤侵蚀量的数量化关系呈现什么状态是备受水土保持研究人员关注的重要问题, 因为二者数量化关系确定可对研究区土壤侵蚀量预判提供相对便捷途径。

十大孔兑水土流失区位于鄂尔多斯北部以中度和重度侵蚀为主^[21-23], 是我国水土流失最为严重的砒砂岩黄土区之一。因此, 本文基于通用土壤流失评价模型USLE应用ANUSPLIN软件对降雨侵蚀力因子插值^[24-25]和遥感手段计算区域土壤侵蚀模数, 采用整体拟合和分段拟合的方法对十大孔兑南部植被盖度和土壤侵蚀模数的关系进行分析, 尝试识别植被盖度变化梯度上土壤侵蚀模数的阈值, 为深入了解植被盖度与土壤侵蚀量的关系提供参考依据。

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

十大孔兑砒砂岩区位于鄂尔多斯北部, 地理位置为108°6′—111°E, 39°50′—40°10′N, 境内有黄河内蒙古段的十个一级支流(图 1), 是黄河内蒙古段的主要产沙支流。由西向东依次为毛卜拉孔兑、布日嘎色太沟、黑赖沟、西柳沟、罕台川、壕庆河、哈什拉川、母哈日沟、东柳沟和呼斯太河。地势表现为西南高东北低；气候属典型大陆性气候, 平均气温6℃左右；多年平均降雨量从西北到东南变化于200—400 mm之间, 年蒸发量为2200 mm^[26]。

1.2 研究方法 1.2.1 数据

气象数据源于中国气象数据网(http://data.cma.cn/), 为十大孔兑周围站点(达拉特旗、杭锦旗、东胜区、包头、呼和浩特、临河、鄂托克旗、右玉和河曲)2000—2017年逐日降水数据。

MODIS NDVI数据源于美国地质调查局(USGS), 空间分辨率为250 m, 产品类型为MOD13Q1 (https://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.005), 时间为2000—2017年生长季(5—8月)的16 d合成数据, 原始数据为HDF格式, 通过MRT(ModisReprojection Tool)进行转换, 将NDVI数据由HDF格式转换为TIFF格式, 投影方式由等面积正弦曲线投影(Sinusoidal Projection)转换为通用横轴墨卡托投影(Universal Transverse Mercator Projection, UTM), 基准面为WGS-84；最后用resample工具形成30 m分辨率数据。DEM数据源于地理空间数据云(http://www.gscloud.cn/), 空间分辨率为30 m。

裸地NDVI采用野外实测数据, 数据获取时间为2017年8月7—18日, 研究区采样位置为24个地点(图 1), 其中每个采样地点设置30个样圆, 选用美国ASD公司的FieldSpecHandHeld手持便携式光谱仪, 波长范围为325—1075 nm。测定时, 探头垂直向下, 视场角为25°, 与实测目标地物的距离为0.67 m, 对应地面样本为随机选取的直径30 cm的样圆。光谱测定时间为10:00—14:00之间, 对同一个样圆裸地反射光谱测定5次, 并用数据处理软件ViewSpec Pro导出所测光谱值, 对应MODIS NDVI数据的红光(620—670 nm)和近红外波段(841—876 nm), 计算平均值。为保证光谱测定的准确性, 每隔10—15min利用标准白板进行校光。测定期间, 光照条件较稳定、风力弱, 对光谱反射值的干扰极小。

土地利用数据基础数据(2000、2005、2010、2015年)源于中科院地理与资源科学研究所, 为了更准确的反映研究区现状对研究区内林地用鄂尔多斯市林业局提供的同期森林资源分布图进行了修订。同时, 为避免植被类型、水保工程治理区等造成的干扰, 在ArcGIS里面仅选取剔除了林地、农田、水保工程治理区的草原区进行研究。

1.2.2 评价模型

采用通用土壤流失USLE评价模型, 方程式如下^[24]：

(1)

式(1)中, A为土壤流失量(t hm^-2a^-1)。R为降雨侵蚀力因子(MJ mm hm^-2 h^-1 a^-1), 筛选日降水量大于12 mm的降水数据依据文献计算^[27]。然后对得到的降雨侵蚀力因子通过ANUSPLIN软件插值形成grid格式的数据, 通过ArcGIS定义坐标体系并转换成栅格数据

K为土壤可蚀性因子(无量纲), 依据圪坨店国家水土保持监测站提供的实测数据和水利部最新标水土流失动态监测技术规定(办水保[2018]189号)计算完成。L为坡长因子, S为坡度因子, 根据GDEMDEM数据提取(分辨率30 m)；C为植被覆盖与管理因子(无量纲), 采用文献所应用方法^[21],

(2)

(3)

在(2)和(3)式中, fc为植被盖度, NDVI_veg为植被完全覆盖时的NDVI值, NDVI_soil为裸地NDVI实测值。

P为水土保持措施因子(无量纲), 同样根据水利部(办水保[2018]189号)规定对各种水土保持工程措施进行赋值, 其中梯田、水平沟和鱼鳞坑分别为0.414、0.335和0.249生成P因子图层。P值范围在0—1, 0代表基本不发生侵蚀, 水保措施很有成效, 而1代表未采取任何控制措施的土地。

1.2.3 趋势分析

在USLE模型中植被盖度与其他因子是以独立变量形式参与计算土壤侵蚀量, 所以植被盖度和土壤侵蚀量之间可能存在的是非线性关系(等级理论中的包容型等级结构)^[28-32], 即：植被盖度和土壤侵蚀模数关系随分析梯度变化产生多种线性规律变化趋势。因此, 本文用整体拟合和分段拟合两种趋势分析方法对二者之间关系进行分析, 整体拟合分析过程中阈值识别可通过对函数求导获得, 分段拟合分析中通过两个方程交点求得阈值。

研究区坡度是在ArcGIS空间分析模块中完成, 坡度分为 < 5°、5—10°和>10°作为控制条件。为了简化分析过程, 分别以侵蚀模数总体(多年平均值)、最大、最小和最接近平均水平特征等四个方面分别进行分析。

2 结果 2.1 土壤侵蚀模数总体变化趋势

从土壤侵蚀模数来看(图 2), 多年土壤侵蚀模数平均值为29.31 t hm^-2 a^-1, 最接近平均值的年份为2010年, 平均土壤侵蚀模数为28.81 t hm^-2 a^-1；2016年土壤侵蚀模数最大, 平均为65.6 t hm^-2 a^-1, 高于平均值36.29 t hm^-2 a^-1；2011年土壤侵蚀模数最小, 为10.95 t hm^-2 a^-1, 低于平均值18.35 t hm^-2 a^-1。从趋势分析来看, 研究期间土壤侵蚀模数在整体上呈增加趋势, 变化率为0.35 t hm^-2 a^-1。不同年份侵蚀模数差异较大, 在2003、2012和2016年分别出现了明显大于平均值的土壤侵蚀模数, 分别为54.95 t hm^-2 a^-1、63.6 t hm^-2 a^-1和65.6 t hm^-2 a^-1, 最低也高于平均值25.65 t hm^-2 a^-1。在2000、2005、2011、2015和2017年出现了明显低于平均值的土壤侵蚀模数, 最高也低于平均值14.54 t hm^-2 a^-1。这意味着土壤侵蚀模数在时间尺度上呈现多种变化趋势。

2.2 植被盖度和土壤侵蚀模数的关系 2.2.1 总体特征

从多年平均特征来看(图 3), 在不同坡度等级随着植被盖度增加土壤侵蚀模数变化在整体上都呈典型的抛物线型(P < 0.001)。

坡度 < 5°时, 土壤侵蚀模数在植被盖度 < 14%时呈增加趋势(图 3), 土壤侵蚀模数变化于12.72—22.56 t hm^-2 a^-1之间；当盖度>14%时, 土壤侵蚀模数呈减少趋势, 变化于2.02—18.66 t hm^-2 a^-1之间。对植被盖度和土壤侵蚀模数关系函数进行求导, 土壤侵蚀模数阈值为18.18 t hm^-2a^-1, 对应植被盖度为11.42%。

坡度为5—10°时, 土壤侵蚀模数在植被盖度 < 18%时呈增加趋势(图 3), 土壤侵蚀模数变化于26.91—35.85 t hm^-2 a^-1之间；当盖度>18%时, 土壤侵蚀模数呈减少趋势, 变化于6.76—33.73 t hm^-2 a^-1之间。对植被盖度和土壤侵蚀模数函数进行求导, 土壤侵蚀模数阈值为34.29 t hm^-2 a^-1, 对应植被盖度为16.51%。

坡度>10°时, 土壤侵蚀模数在植被盖度 < 18%时呈增加趋势(图 3), 土壤侵蚀模数变化于52.87—75.59 t hm^-2 a^-1之间；当盖度>18%时, 土壤侵蚀模数呈减少趋势, 变化于12.23—72.61 t hm^-2 a^-1之间。对植被盖度和土壤侵蚀模数函数进行求导, 土壤侵蚀模数阈值为74.56 t hm^-2 a^-1, 对应植被盖度为16.5%。

上述结果表明, 无论坡度 < 5°、5—10°还是坡度>10°, 植被盖度和土壤侵蚀模数变化趋势一致(抛物型), 但是在不同坡度条件下, 土壤侵蚀模数阈值对应植被盖度不同, 坡度 < 5°时阈值对应植被盖度比坡度为5—10°和坡度>10°时阈值对应植被盖度小。

2.2.2 典型年份特征

从三种典型年份植被盖度和土壤侵蚀模数关系来看, 最接近平均年(2010)平均土壤侵蚀模数随着植被盖度增加在整体上也呈抛物线型变化(P < 0.001)(图 4)。在坡度 < 5°、5—10°和坡度>10°时三个级别上, 土壤侵蚀模数阈值对应植被盖度分别为8.73%、11.64%和15.37%。

最小年(2011)平均植被盖度和土壤侵蚀模数关系较复杂(图 4), 坡度 < 5°和>10°时, 植被盖度和土壤侵蚀模数关系呈多峰值不显著变化(P> 0.05)。坡度 < 5°时, 土壤侵蚀模数阈值对应植被盖度分别为6.27%和27.57%；坡度为5—10°时, 植被盖度和土壤侵蚀模数关系呈抛物线型变化, 土壤侵蚀模数阈值对应植被盖度为11.24%；坡度>10°时, 土壤侵蚀模数阈值对应植被盖度分别为6.7%和21.67%。

最大年(2016)年植被盖度和土壤侵蚀模数的关系在整体上呈典型的等级结构(图 4), 可以分为土壤侵蚀模数减少区、相对稳定区和减少区三个变化区。坡度 < 5°和5—10°时, 土壤侵蚀模数随着盖度增加呈极显著单调递减趋势(P < 0.001), 需采用分段拟合分析。在坡度 < 5°级别, 盖度 < 12%时, 线性方程为y=-2.8797x+73.295(P < 0.01)；盖度为14%—28%时, 线性方程为y=-0.4578x+51.452(P < 0.01)；两条线的交点对应植被盖度为9.02%。盖度>28%时, 线性方程为y=-1.4315x+83.062(P < 0.01), 与植被盖度为14%—28%时的趋势线的交点对应植被盖度为32.46%。坡度为5—10°的区域, 盖度 < 8%时, 线性方程为y=-3.106x+105.56(P> 0.05), 盖度为10%—24%时, 线性方程为y=0.1628x+78.022(P < 0.01), 交点对应植被盖度为8.42%；盖度>24%时, 线性方程为y=-1.3477x+114.95(P < 0.001), 与植被盖度为10%—24%时趋势线交点对应植被盖度为24.45%。坡度>10°时, 通过对函数求导得到土壤侵蚀模数阈值对应植被盖度为12.98%和25.74%。

上述结果表明, 最大和最小年份与总体平均水平变化趋势相差较大, 而在研究周期(18年)中, 较大偏离平均值出现了3次, 较小偏离出现了5次(图 2), 说明区域土壤侵蚀模数变化受最大和最小年份影响较大。

3 讨论

土壤侵蚀状况与许多因素有关, 如气候条件、植被状况、土壤性质和人类活动等^{[9, 16, 33-34]}。从与本文有关的研究结论来看, 土地利用类型转换(如草地转为人工林)和状况改善(低覆盖草地向高覆盖草地转换)等导致的植被盖度变化都会影响地表径流进而影响土壤侵蚀变化^{[4-5, 11, 18, 33]}。

从研究区土壤侵蚀模数时间变化趋势来看(图 2), 最大降水年(2016)和最小降水年(2011)土壤侵蚀模数分别为最大(65.6 t hm^-2 a^-1)和最小(10.95 t hm^-2 a^-1), 说明降水变化是影响土壤侵蚀模数的最主要因子。从土壤侵蚀模数(图 2中数值)与同期年降水量相关关系来看, 呈极显著(y = 0.2032x - 15.666, R= 0.8599, P < 0.001), 其中2003年、2012年和2016年的较大降水量导致了较高土壤侵蚀模数(54.95 t hm^-2 a^-1、63.6 t hm^-2 a^-1和65.6 t hm^-2 a^-1), 拉高了研究时期的线性变化率。

从植被盖度与土壤侵蚀模数关系来看(图 3), 在总体上呈极显著抛物线型趋势(P < 0.001), 这一结果与朱冰冰等的结果一致, 但与张孝中等和曹梓豪等的植被盖度和土壤侵蚀模数呈幂函数的结果相反^{[6, 35-36]}。歧义原因一是随着地域不同, 植被盖度和土壤侵蚀模数会发生多种变化趋势^[28-32], 二是在植被盖度较小的区域, 由于多年面状侵蚀导致土层较薄且表面混杂大量砾石, 可侵蚀土壤少、土壤贫瘠无法提供足够养分维持植被生长, 从而导致植被盖度降低的恶性循环。本文研究区低植被盖度区几乎全属于近似裸露的砒砂岩区, 从而出现了图 3植被盖度较低且土壤侵蚀量也较少的情况。此外, 本研究结果与吴光艳等在南方红壤区和李斌等以整个黄土高原区所做的研究结果“植被盖度越高, 土壤侵蚀模数越小的变化趋势”相悖^{[8, 16]}, 产生这个原因同样是砒砂岩区特殊的地质和土壤条件所致。

从不同坡度级别土壤侵蚀特点来看, 坡度越小土壤侵蚀模数的阈值对应植被盖度越小。原因是坡度越大, 水流速度和径流系数越大, 植被阻力系数越小^[37-38], 2010年的阈值也证明了这一结果是与阻力系数有关(图 4)。与已有研究阈值为植被盖度25%和20%的结果相比^{[6, 39]}, 本研究结果阈值较小, 产生原因一是本文做了严格条件控制, 排除了乔灌林、水保工程措施和农田(梯田)等干扰；二是植被盖度计算方法不同, 本文采用野外实测裸地光谱数据和MODIS数据相结合来计算研究区植被盖度, 可更客观反应真实植被状况；而采用模拟控制实验所得结果为理论值^{[6, 39]}, 与土壤侵蚀发生过程可能会有差异。在植被盖度为32%—36%时, 土壤侵蚀模数出现了与变化趋势不符的异常高值(图 3), 原因可能是受植被阻力系数影响, 在植被盖度较小时, 土壤直接随水分流走, 而在较大植被盖度时, 植被阻碍作用产生土壤堆积, 当遇到较大降雨时由于可冲刷量较多发生了较大土壤侵蚀模数(即随降雨量增加土壤侵蚀量增加现象所致)。因此可以推断, 在黄土砒砂岩区, 土壤侵蚀模数不仅受到了USLE模型中诸因子的影响, 而且也受到了土壤可侵蚀量限制。

从典型年份植被盖度与土壤却侵蚀模数关系来看(图 4), 在最接近平均年植被盖度和土壤侵蚀模数变化趋势和总体年份变化趋势最为相似, 在最大年土壤侵蚀模数随着植被盖度增加逐渐减少的趋势最显著, 而最小年植被盖度和土壤侵蚀模数的关系和总体关系相差甚远, 说明了区域多年土壤侵蚀量变化主要是受最大土壤侵蚀模数事件影响, 而最小和中等年份对其影响较小。因此可以推断在黄土砒砂岩区多年平均土壤侵蚀模数主要是受较大土壤侵蚀模数年份的影响, 判断区域土壤侵蚀状况与植被盖度关系时应高度关注较大土壤侵蚀年份。

4 结论

2000—2017年十大孔兑砒砂岩区土壤侵蚀模数在时间尺度上呈现复杂的变化趋势, 多年土壤侵蚀模数平均值为29.31 t hm^-2 a^-1；最大和最小值分别为65.6 t hm^-2 a^-1和10.95 t hm^-2 a^-1。从总体来看, 土壤侵蚀模数随着植被盖度增加呈极显著抛物线型变化趋势(P < 0.001)；在坡度级别分别为 < 5°、5—10°和>10°时, 土壤侵蚀模数的阈值(18.18 t hm^-2a^-1、34.29 t hm^-2a^-1和74.56 t hm^-2 a^-1)对应的植被盖度分别为11.42%、16.51%和16.5%。结果解释了在砒砂岩区土壤侵蚀模数不仅受到了USLE模型中诸因子的影响, 而且也受到了土壤可侵蚀量限制。此外, 判断区域土壤侵蚀与植被盖度关系时应高度关注较大侵蚀量较大年份。