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

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

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

A=f·R·K·LS·CP(1)

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

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

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

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

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

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

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

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

式中:NDVI为归一化植被指数; ρ_NIR,ρ_R分别为近红外波段和红波段反射率; f为植被覆盖度(%); NDVI_max,NDVI_min分别为NDVI的最大和最小值。

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

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

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

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

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

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

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

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

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

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/(km²·a)。其中母哈尔河土壤侵蚀模数增量最大,毛不拉孔兑相对较小。

表4 十大孔兑上游各孔兑流域土壤侵蚀模数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/(km²·a)。其中母哈尔河土壤侵蚀模数减少量最为显著,毛不拉孔兑减少量最小。

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

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

图4 基于两种方法的十大孔兑上游土壤侵蚀强度面积对照

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

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