十大孔兑上游土壤侵蚀空间分布及动态变化

刘 璐, 郭月峰, 姚云峰, 刘晓宇, 祁 伟,2

(1.内蒙古农业大学 沙漠治理学院, 呼和浩特 010018; 2.内蒙古自治区水利水电勘探设计院, 呼和浩特 010020)

土壤侵蚀是土壤在风力、水力、重力等自然营力与人类活动综合作用下被剥蚀、搬运、沉积的过程，其土壤养分流失，土地生产力下降，土壤退化，引起泥沙淤积，加剧洪涝、干旱等灾害发生，严重危及人类生存与发展，是各国关注的环境热点问题之一[1]，为此各国学者采用了多种技术对土壤侵蚀进行了大量探索性试验和深入研究，其中土壤侵蚀定量评价是土壤侵蚀研究中的主要方法，尤其是近几十年来，遥感技术和地理信息系统的快速发展极大地推动了土壤侵蚀定量化研究[2]，为土壤侵蚀研究提供了有力的技术支持，已在世界许多国家得到广泛应用和研究[3-6]。

国内对GIS和RUSLE的综合应用始于20世纪80年代末，RUSLE是将降雨侵蚀力、土壤可蚀性、植被覆盖、坡长、坡度和水土保持措施作为6大因子对土壤侵蚀进行定量评价的水蚀预报经验模型。卜兆宏等[7]完善了模型因子算式算法，改进了土壤粒级换算方法，在水土流失定量遥感方法中取得新进展，并在太湖流域加以应用，还有其他学者利用USLE或RULSE模型估算了广东省[8]、江西省[9]的土壤侵蚀强度及中国的[10]年平均土壤侵蚀率，基于RULSE模型在不同省份和地区相继出现了一系列研究报道[11-14]。基于区域土地类型对其地面坡度、土地利用/覆被类型和植被覆盖数据构建土壤侵蚀遥感反演模型，从而进行土壤水力侵蚀强度定性分析，郭碧云等[15]以坡度、植被覆盖度和土地利用方式3个方面作为水力侵蚀的判别指标，结合野外调查和室内分析，对内蒙古农牧交错区土壤水蚀强度进行评价。基于土壤侵蚀经验模型运用GIS定量评估土壤侵蚀已成为合理规划、利用区域水土资源的有效手段。

十大孔兑上中游区域土壤侵蚀严重，水土流失面积约占上中游总面积的95%，多年平均输沙量0.27亿t，是黄河河道主要的产沙区之一。区域内土壤类型与地表物质具有较大差异，其水土流失特征及类型也各不相同。上游地处砒砂岩外露区，汛期暴雨集中时，在水力侵蚀的作用下[16]，风化物源源不断地输送到中游和下游，是河道粗砂重要的补给，也是下游防汛、取水安全的重要隐患，故水力侵蚀是该区域水土流失的重要原因。

本文选取侵蚀类型为水力侵蚀的十大孔兑上游作为土壤侵蚀演算的研究区，基于GIS采用两种方法分别构建了十大孔兑上游流域1980年、1990年、2000年、2010年及2018年5个典型年份的土壤侵蚀反演模型。第1种方法在土壤侵蚀分类分级标准基础上，基于十大孔兑上游地形、土地利用、植被覆盖度数据对土壤侵蚀强度进行定性分析。第2种方法基于RUSLE演算土壤侵蚀模数并进行土壤侵蚀强度分级，并通过对比两种方法的有效性为生态脆弱区流域治理提供理论依据。

1 研究区概况

十大孔兑位于内蒙古鄂尔多斯市境内，地理坐标为108°47′—110°58′E，39°47′—40°30′N，其发源于鄂尔多斯台地，流经库布齐沙漠，横穿下游冲积平原后汇入黄河的10条一级支流，从西向东依次为毛不拉孔兑、布尔色太沟、黑赖沟、西柳沟、罕台川、壕庆河、哈什拉川、母哈尔河、东柳沟、呼斯太河。研究区属于典型干旱大陆性气候，气候干燥，风大沙多，降水年际变化大，年内分布不均，7—9月的降水量占全年降水量的70%，由南向北逐渐减少，其中风沙区雨量最少。地貌类型以砾质丘陵、沙质丘陵、冲洪积平原为基础，由南向北可分为上游黄土丘陵沟壑区、中游库布齐沙漠风沙区和下游黄河南岸冲积平原区。上游为鄂尔多斯台地北缘，沟壑纵横，地表支离破碎，以水力侵蚀为主；中游为库布齐沙漠区，横贯东西，西宽东窄，风力侵蚀极为严重；下游为冲洪积平原区，地势平坦，多年平均径流量、输沙量较大。流域内暴雨洪水频繁，水土流失严重，植被稀疏，生态环境十分脆弱(图1)。

图1 研究区示意图

2 数据与方法

2.1 基础数据来源与处理

由于十大孔兑在降水、植被、土地利用等方面资料较为匮乏，本研究的基础数据主要依靠遥感技术获取。十大孔兑流域1980年、1990年、2010年、2018年降水数据选用地面测站数据，基于泰森多边形法将测站降水数据插值成面状数据，因2000年地面测站降水数据短缺，选用TRMM降雨数据替代，通过对2000年1—12月的TRMM3B43降雨数据统计计算，获取2000年的年降雨数据。收集覆盖十大孔兑流域1980年、1990年、2000年、2010年、2018年的6—10月Landset MSS/TM/ETM/landsat8影像数据25期，选取影像质量较好的8景影像作为5个典型年份的遥感数据源。收集覆盖十大孔兑流域的DEM数据1套，空间分辨率约为30 m×30 m。收集覆盖十大孔兑流域的土壤类型数据1套，比例尺为1∶100万。

2.2 研究方法

2.2.1 基于土地类型土壤侵蚀研究 基于《土壤侵蚀分类分级标准》(SL190—2007)，利用地面坡度、土地利用/覆被类型和植被覆盖数据，构建土壤侵蚀遥感反演模型，按照土地利用类型、植被覆盖度、地面坡度3种参数在不同类型或区间内的归属不同来划分土壤侵蚀强度，利用ArcGIS软件的Model Builder模块，构建土壤侵蚀模型，针对栅格图像建立3种参数的类型或区间的判断矩阵，满足条件的像元，通过类型合并，获取规定的6类土壤侵蚀强度数据[17-18]。

2.2.2 基于RUSLE土壤侵蚀研究 采用RUSLE法对十大孔兑上游近40 a来土壤侵蚀模数变化展开定量研究，并结合水土保持措施实施情况分析其产生的水土保持效益。研究利用ArcGIS软件的Model Builder模块，分别构建各土壤侵蚀因子计算模型，然后将以上6个土壤侵蚀因子和转换常量相乘后获取土壤侵蚀模数结果[19]。

水土流失定量测算采用RUSLE方程，其表达式如下：

A=f·R·K·LS·CP

(1)

式中:A为年土壤流失量[t/(km2·a)]；f为使A的单位转换为国际通用单位的常数,其值为224.2[20]；R为降雨侵蚀力因子;K为土壤可侵蚀性因子;LS为坡长坡度因子的乘积;CP为植被覆盖与管理因子、水土保持因子的乘积。R,K量纲单位与美国通用水土流失方程相同；LS,CP无量纲单位。

(1) 降雨侵蚀力因子(R)。降雨侵蚀力是降水对土壤造成侵蚀的动力指标[21]，是侵蚀的主要驱动力[22]，本文R因子的计算公式为:

R=4.1×F-152

(2)

(3)

式中:F为修正参数;Pj为第j月的降雨量(mm)；p为年平均降雨量(mm)。将计算所得的各气象站点的年降雨侵蚀力在ArcGIS中进行克吕格空间插值[23-24],得到研究区年降雨侵蚀力图。

(2) 土壤可蚀性因子(K)。土壤可蚀性反映了土壤对侵蚀营力分离和搬运作用的敏感性[25-26]，K值越高表示土壤抗蚀能力越弱，反之K值越低表示土壤抗蚀能力越弱。计算公式如下[27]：

(4)

Sn=1-Sa/100

(5)

式中:Sa为砂粒含量百分比(%);Si为粉粒含量百分比(%);Cl为黏粒含量百分比(%);C0为有机碳含量百分比(%)。在ArcGIS中,以区域土壤类型图为基础,根据K值对应的土壤类型进行赋值,最后以K值作为字段进行数据转换。

(3) 坡长坡度因子(LS)。LS因子是坡长因子和坡度因子叠加相乘的结果，借助于GIS技术，对填洼处理后的数字高程模型提取坡长和坡度值[27]，利用ArcGIS的表面分析来提取坡度，同时运用非累积流量的算法确定坡长，并按格网单元的非累计坡长再进行累加[28]。通过运算分别生成坡长、坡度因子图，在ArcGIS中使用栅格计算器运算后最后生成LS空间分布图。计算公式为：

(6)

L=(λ/22.1)m

(7)

(8)

式中:L和S分别为坡长和坡度因子;θ为利用数字高程模型(digital elevation model,DEM)数据提取的坡度值(°);m为坡长指数;λ为坡长(m)。

(4) 植被覆盖与管理因子(C)。植被覆盖与管理因子表示植被覆盖和管理措施对土壤侵蚀的作用[29]。计算公式为：

(9)

(10)

(11)

式中:NDVI为归一化植被指数;ρNIR，ρR分别为近红外波段和红波段反射率;f为植被覆盖度(%);NDVImax，NDVImin分别为NDVI的最大和最小值。

(5) 水土保持措施因子(P)。水土保持措施因子指在特定水土保持措施的土壤流失与起伏地耕作的相应土壤流失之比[30]，在RUSLE模型中，P值介于0～1，其中未采取任何水土保持措施的地块P值为1，不发生土壤侵蚀的地块则P值为0[31](表1)。

表1 不同土地利用/覆被类型P值

3 结果与分析

3.1 基于土地类型的土壤侵蚀分析

十大孔兑上游1980—2018年土壤侵蚀强度面积统计结果见表2、图2。1980年十大孔兑上游轻度侵蚀等级及以上侵蚀类型面积2 436.42 km2，1990年为2 722.28 km2，2000年为2 717.23 km2，2010年为2 638.21 km2，2018年为2 462.51 km2。1980—1990年十大孔兑上游轻度侵蚀等级及以上侵蚀类型面积增加285.86 km2；1990—2000年侵蚀类型面积减少5.05 km2；2000—2010年侵蚀类型面积减少79.02 km2；2010—2018年侵蚀类型面积减少175.70 km2。

表2 基于土地类型十大孔兑上游1980-2018年土壤侵蚀强度面积 km2

图2 基于土地类型十大孔兑上游1980-2018年土壤侵蚀强度空间分布

3.2 基于RUSLE的土壤侵蚀研究

3.2.1 十大孔兑上游水土流失情况分析 十大孔兑上游1980—2018年土壤侵蚀强度面积统计信息见表3、图3。1980年十大孔兑上游轻度侵蚀等级及以上侵蚀类型面积为1 721.59 km2，1990年为2 889.63 km2，2000年为2 622.81 km2，2010年为2 783.77 km2，2018年为2 686.10 km2。1980—1990年十大孔兑上游轻度侵蚀等级及以上侵蚀类型面积增加了1 168.04 km2，1990—2000年减少了266.82 km2，2000—2010年增加了160.96 km2，2010—2018年减少了97.67 km2。

表3 基于RUSLE十大孔兑上游1980-2018年土壤侵蚀强度面积 km2

图3 不同年份十大孔兑上游土壤侵蚀模数及侵蚀强度空间分布

1980年十大孔兑上游平均土壤侵蚀模数为2 744.88 t/(km2·a)；1990年十大孔兑上游平均土壤侵蚀模数为14 557.51 t/(km2·a)；2000年十大孔兑上游平均土壤侵蚀模数为3 058.06 t/(km2·a)；2010年十大孔兑上游平均土壤侵蚀模数为4 884.66 t/(km2·a)，2018年十大孔兑上游平均土壤侵蚀模数为2 997.03 t/(km2·a)。5个典型年份土壤侵蚀模数由高到低依次为1990年、2010年、2000年、2018年、1980年。1980—1990年十大孔兑上游平均土壤侵蚀模数增加了11 812.63 t/(km2·a)，1990—2000年减少了11 499.45 t/(km2·a)，2000—2010年增加了1 826.60 t/(km2·a)，2010—2018年减少了1 887.63 t/(km2·a)。

综上，1990—2018年十大孔兑上游轻度及以上等级土壤侵蚀面积由2 889.63 km2下降至2 686.10 km2，土壤侵蚀模数由14 557.51 t/(km2·a)下降至2 997.03 t/(km2·a)，该区域土壤侵蚀面积减小，土壤侵蚀强度下降。

3.2.2 十大孔兑上游各孔兑水土流失情况分析 十大孔兑流域内的各孔兑上游平均土壤侵蚀模数统计信息见表4。1980年十大孔兑上游各孔兑土壤侵蚀模数由高到低依次为哈什拉川、罕台川、母哈尔河、布尔色太沟、黑赖沟、西柳沟、毛不拉孔兑、壕庆河、东柳沟、呼斯太河。1990年十大孔兑上游各孔兑土壤侵蚀模数由高到低依次为母哈尔河、哈什拉川、罕台川、壕庆河、东柳沟、西柳沟、黑赖沟、呼斯太河、布尔色太沟、毛不拉孔兑。2018年十大孔兑上游各孔兑土壤侵蚀模数由高到低依次为西柳沟、母哈尔河、毛不拉孔兑、呼斯太河、黑赖沟、壕庆河、罕台川、哈什拉川、东柳沟、布尔色太沟。

表4 十大孔兑上游各孔兑流域土壤侵蚀模数 t/(km2·a)

1980—1990年毛不拉孔兑、布尔色太沟、黑赖沟、西柳沟、罕台川、壕庆河、哈什拉川、母哈尔河、东柳沟、呼斯太河的土壤侵蚀模数分别增加了1 881，2 797，4 812，6 182，22 449，18 932，21 576，36 725，11 944，6 048，11 813 t/(km2·a)。其中母哈尔河土壤侵蚀模数增量最大，毛不拉孔兑相对较小。

1990—2018年，毛不拉孔兑、布尔色太沟、黑赖沟、西柳沟、罕台川、壕庆河、哈什拉川、母哈尔河、东柳沟、呼斯太河土壤侵蚀模数分别减少了2 130，2 532，3 160，3 857，23 319，17 371，24 101，36 833，11 910，12 003，6 236 t/(km2·a)。其中母哈尔河土壤侵蚀模数减少量最为显著，毛不拉孔兑减少量最小。

综合5个典型年份土壤侵蚀情况，1980—2018年母哈尔河、哈什拉川、罕台川三条孔兑上游土壤侵蚀最为严重，东柳沟、呼斯太河两条孔兑上游土壤侵蚀相对较轻。其中1990年土壤侵蚀程度尤为显著，主要是由于20世纪90年代还未对母哈尔河、哈什拉川、罕台川孔兑开展大规模水土保持综合治理等措施，且该区域所发生的“89·7·21”特大暴雨，其最大雨强达119.1 mm/h，洪峰流量达3 090 m3/s，引发特大洪水造成罕台川水土流失严重。

3.3 基于两种方法对比分析

基于土地类型与RUSLE方程的1980—2018年十大孔兑上游土壤侵蚀强度反演结果如图4所示。1990年、2000年、2010年、2018年两种方法的轻度侵蚀强度及以上等级的土壤侵蚀反演结果基本相同，1980年两种方法的轻度侵蚀强度及以上等级的土壤侵蚀反演结果差别较大。基于土壤侵蚀分类分级标准遥感反演的5个典型年土壤侵蚀强度以无明显侵蚀和中度侵蚀为主。基于RUSLE方程遥感反演的5个典型年无明显侵蚀等级面积最多，且除1990年外，其余典型年无明显侵蚀、轻度侵蚀、中度侵蚀、强度侵蚀、极强度侵蚀、剧烈侵蚀所占面积均呈现逐渐减少趋势。

注：方法一为基于土地类型；方法二为基于RUSLE方程。

1980—1990年十大孔兑上游植被覆盖度呈降低趋势，轻度侵蚀等级及以上侵蚀类型面积、平均土壤侵蚀模数呈增加趋势；1990—2000年植被覆盖度呈增加趋势，轻度侵蚀等级及以上侵蚀类型面积、土壤侵蚀模数呈较少趋势；2000—2010年植被覆盖度呈增加趋势，轻度侵蚀等级及以上侵蚀类型面积、土壤侵蚀模数呈增加趋势；2010—2018年植被覆盖度呈增加趋势，轻度侵蚀等级及以上侵蚀类型面积、土壤侵蚀模数呈减少趋势。

土壤侵蚀和植被覆盖度呈负相关关系，即相同条件下，植被覆盖度越高土壤侵蚀越轻，反之越高。1980—2000年十大孔兑上游植被覆盖度与轻度侵蚀等级及以上侵蚀类型面积、平均土壤侵蚀模数呈负相关关系，但2000—2018年没有呈负相关关系，这是因为土壤侵蚀除受植被覆盖度影响外，还与降雨、土壤类型、地形、实施水土保持措施等因素有关。两种方法对比发现，RUSLE方程考虑因子更为全面，模拟更接近实际，同时适用于水土流失定量研究。

4 讨 论

十大孔兑上游从20世纪50，60年代开始开展了水土保持生态建设，中游开展了防风固沙林带和引洪淤地工程建设；60年代陆续开展种树种草项目，70年代实施“逐步退耕还林还牧，林牧为主，多种经营”的治理方针，80年代贯彻“三种五小”和“资源开发与环境保护并重”两个治理方针，90年代以后进一步加大林草治理力度，在利用柠条改良草场、沙棘治理砒砂岩沟壑等方面取得了显著成就。但在2000年以前，受到人力、财力及管理权限的限制，鄂尔多斯林草的破坏速度大于恢复速度，治理成效甚微，生态环境一直处在局部治理、整体恶化状态，1990年十大孔兑上游土壤侵蚀模数为14 557.51 t/(km2·a)，处于5个典型年份中土壤侵蚀模数最大值，侵蚀等级主要以极强度侵蚀与剧烈侵蚀为主。

1999—2012年孔兑的草原植被得到明显的恢复，其平均高度、平均覆盖度呈现上涨趋势。根据鄂尔多斯市草原管理局测定，封禁前退化放牧草场植被平均高度为10～28 cm，覆盖度为15%～30%，每平方米内一般有4～6种植物，而在封禁后，经过一定时间的恢复期，其平均高度可达35～60 cm以上，平均覆盖度可达65%～75%以上，每1 m2内物种数可达到7～12种以上。此外气候因素(降水、气温、风速等)对林草植被恢复有一定的影响，如2010年由于降水较少，导致其草场返青日期、平均高度、平均覆盖度等各项指标低于2009年。

1990—2018年十大孔兑上游轻度及以上等级土壤侵蚀面积由2 889.63 km2下降至2 686.101 km2，土壤侵蚀模数由14 557.51 t/(km2·a)下降至2 997.03 t/(km2·a)。由于该区域20世纪80年代先后开展了黄河中游试点小流域、骨干坝工程、沙棘示范区、罕台川水土保持综合治理、黄土高原水土保持世行贷款项目、沙棘拦沙工程等水土保持项目。近年来，鄂尔多斯市推行了封山禁牧制度，生态修复成效显著。十大孔兑上中游水土保持治理措施包括工程措施、林草措施、防风固沙措施、其他水土保持措施等。工程措施包括淤地坝、谷坊、梯田、引洪淤地、沟头防护等；植物措施包括水土保持林(乔木林、灌木林、经济林)、人工种草、封育治理；防风固沙措施包括造林、种草和沙障等，防风固沙措施全部集中在十大孔兑中游的风沙区。截至2017年，流域上中游丘陵风沙区综合措施总面积5 826.15 km2，其中人工造林面积4 378.92 km2，建设淤地坝511座，控制面积984.36 km2，总库容2.05亿m3，拦沙库容1.11亿m3。由此可见，近些年十大孔兑上游土壤侵蚀得到了有效治理，水土保持治理工作取得了一定成效。如图5—6所示，2000—2018年C值与P值的空间分布有所增加，C值的确定与植被类型、植被覆盖度相关，P值即土地利用/覆被方式能够很好地表征人类活动对地表水土流失的作用效果。

图5 十大孔兑上游2000年与2018年C因子变化

图6 十大孔兑上游2000年与2018年P因子变化

基于两种方法演算不同侵蚀等级面积差异较大，一方面由于两类演算方法机理不同，考虑的影响因子亦不同。基于土地类型从坡度、植被覆盖度、地类定性划分，方法侧重于定性评价。基于RUSLE方程方法从水土流失机理出发，考虑多个影响因子(包括降雨、地形、土壤类型、植被覆盖、水保措施等)，且各影响因子有多种计算方法，各因子不同的计算方法计算结果也有差别，该方法侧重于定量研究。另一方面两种反演方法中使用的遥感影像不可避免会出现云雾、日照等影响，影像的灰度值会有差别，从而带来植被覆盖度反演误差。

5 结 论

(1) 基于土地类型法1980年十大孔兑上游轻度及以上等级土壤侵蚀面积为2 436.42 km2，1990年为2 722.28 km2，2000年为2 717.23 km2，2010年为2 638.21 km2，2018年为2 462.51 km2。基于RUSLE法1980年十大孔兑上游轻度及以上等级土壤侵蚀面积为1 721.59 km2，1990年为2 889.63 km2，2000年为2 622.81 km2，2010年为2 783.77 km2，2018年为2 686.10 km2。基于两种方法计算轻度及以上等级土壤侵蚀面积都分别呈现1980—1990年增加、1990—2018年减少趋势。

(2) 基于RUSLE法1980年十大孔兑上游平均土壤侵蚀模数为2 744.88 t/(km2·a)；1990年为14 557.51 t/(km2·a)；2000年为3 058.06 t/(km2·a)；2010年为4 884.66 t/(km2·a)，2018年为2 997.03 t/(km2·a)。十大孔兑上游平均土壤侵蚀模数呈现1980—1990年增加、1990—2000年降低、2000—2010年增加、2010—2018年减少趋势，其中母哈尔河、哈什拉川、罕台川三条孔兑上游土壤侵蚀最为严重，东柳沟、呼斯太河两条孔兑上游土壤侵蚀相对较轻。

(3) 两种方法对比发现，RUSLE方程考虑因子更为全面，模拟更接近实际，同时适用于水土流失定量研究。