基于改进RUSLE模型的西南土石山区水土流失评价——以湖南省龙山县为例

龙思佳，汤媛媛,2*，戴亮亮，乔 双，樊旺东，佘 雄，孔巍巍

1.中国地质调查局长沙自然资源综合调查中心，湖南长沙 410600；2.中国地质大学（北京）地球科学与资源学院，北京 100083

湘西土家族苗族自治州地处武陵山区，位于湘、鄂、黔、渝4省(市)交界处，土壤侵蚀类型分区属于西南土石山区，是我国西南典型岩溶山区，也是长江中上游地区的岩溶区（潘佑堂和杨胜刚，1995；吴兆清和苏绍眉，2004）。受自然条件和人为因素影响，湘西土家族苗族自治州的水土流失问题十分严重。该地区位于长江中游地区的沅水流域，是湖南省水土流失问题相对集中和严重的地区之一。2010年，龙山县成为湘西土家族苗族自治州水土流失最为严重的县市之一，整体水土流失状况以轻度和中度侵蚀为主，强烈侵蚀次之（尹黎明等，2012）。水土流失导致当地生态环境破坏严重，土壤流失和退化，土壤蓄水抗旱能力降低，对人类的生存发展造成极大的影响。当地百姓流传一句话：“水在地下流，人在岸上愁”。防治水土流失、加大水土保持力度对于生态文明建设和推动湘西土家族苗族自治州经济社会可持续发展有着重大而深远的意义。

土壤侵蚀模型是检测和预报土壤流失的重要工具。美国的修正通用土壤流失方程RUSLE 是目前较为通用的土壤侵蚀模型（陈朝良等，2021；张素等，2021；高凡洁等，2022；王丰等，2022），已应用于我国多地土壤侵蚀时空变化分析和水土保持研究。RUSLE 土壤流失方程共包含5 个土壤侵蚀因子：降雨侵蚀力因子（R）、土壤可蚀性因子（K）、坡度坡长因子（LS）、植被覆盖管理因子（C）、水土保持措施因子（P）。植被覆盖管理因子（C）作为抑制水土流失的重要因素，以植被类型和植被覆盖度对其的影响最大（陈学兄，2013）。植被覆盖度一般由植被指数反演而来，是描述地表植被状况的重要参数，是土壤侵蚀模型的重要输入参数，也是当今研究植被因子与水土流失关系中用得最多的一个参数（张灿等，2015；Jia K et al.，2016），其估算精度影响土壤侵蚀预测的精度（吴志杰和徐涵秋，2011）。然而，植被覆盖度遥感估算精度与植被指数等紧密相关，且大多数的植被指数都没考虑地形的影响（刘亚迪等，2015）。目前植被覆盖管理因子的估算均使用均一化植被指数（NDVI）来提取植被覆盖度，能够较好反映区域植被情况。但是在地形起伏的山区，遥感影像中有着大量的阴影分布，利用NDVI 计算的结果会影响土壤侵蚀模型的精度，存在一定误差（孙桂凯等，2021；罗杰等，2022）。龙山县多山丘，地形起伏大，其引起的光谱信息变化影响山地植被覆盖度的提取精度，从而成为植被覆盖管理因子准确估算的主要障碍。

目前用于解决山地植被信息提取的方法有：对影像地形校正后计算植被指数、复合植被指数、地形调节植被指数和归一化差值山地植被指数（NDMVI）（吴志杰和徐涵秋，2011；吴志杰等，2016，2017），其中归一化差值山地植被指数可有效提高山地植被信息提取的精度，该指数已被广泛应用于南方丘陵山区植被覆盖度的遥感估算（吴志杰等，2016，2017；陈学兄等，2020）。而这些植被指数多应用于植被覆盖度的提取，较少应用于水土流失监测与评价。本研究以Landsat5 TM 和Landsat8 OLI 为遥感数据源，尝试将NDMVI 应用于西南土石山区植被覆盖管理因子的提取，探讨改进的RUSLE 模型在龙山县水土流失监测与评价中的应用，并对龙山县2000 年和2020 年土壤侵蚀进行估算，分析龙山县土壤侵蚀空间分布特征，实现对龙山县土壤侵蚀动态变化的快速定量监测，旨在为今后以龙山县为代表的西南土石山区土壤侵蚀监测和治理、水土保持规划等提供科学依据。

1 研究区概况及数据获取

1.1 研究区概况

龙山县位于湘西土家族苗族自治州西北部，地处武陵山脉腹地，座落于湘、鄂、渝三省（市）交界处，地理坐标为109°10′～109°53′E，28°45′～29°30′N（图1）。全县地形以山地为主，地势北高南低，东陡西缓，海拔218.2～1736.5 m。龙山县属亚热带大陆性湿润季风气候，年平均温度15.7 ℃，极端最高温度41.8 ℃，极端最低温度-14 ℃，无霜期238～333天，年均降雨量1046.2～1740 mm。龙山县土壤主要由灰岩、白云岩、板页岩、河流冲积物（包括古河流和近代河流冲积物）、紫色砂页岩等富含矿物质营养元素的母质风化发育而成。

图1 研究区地理位置图Fig.1 Geographical location of the study area

1.2 数据来源

研究区基础数据包括：来自美国地质调查局的2000年Landsat 5、2020年Landsat 8 多光谱影像；数据来源于NASA EARTHDATA 的数字高程模型ALOS DEM，空间分辨率为12.5 m；来自中国土壤数据库的研究区土壤类型矢量图；来源于Terra Climate 的2000年和2020年研究区降雨数据，精度4 km；中国科学院资源环境科学与数据中心提供的研究区土地利用类型图。遥感数据预处理包括：对两期影像数据进行镶嵌、裁剪、正射校正和辐射校正处理，将影像数据各波段像元灰度值转换为表观反射率值。

2 研究方法

2.1 研究原理

2.1.1 山地植被指数（NDMVI）提取

吴志杰和徐涵秋（2011）提出一种归一化差值山地植被指数(NDMVI)，该植被指数由NDVI 变换而来，它不需要DEM的支持，仅依据植被光谱特征和遥感影像数据。该指数能有效减弱地形效应，可用于解决复杂地形对山地植被信息的影响问题。NDMVI的计算公式为：

式中：ρnir和ρred为近红外波段和红光波段的表观反射率；Rmin为红光波段表观反射率的最小值；NIRmin为近红外波段表观反射率的最小值。为避免其取值的主观性以及不确定性，不剔除建成区、水体等信息。

2.1.2 土壤侵蚀模型RUSLE

采用修正通用土壤流失方程(RUSLE)对研究区的土壤侵蚀进行定量评估，RUSLE 模型可用如下公式表达：

式中：A为年土壤侵蚀模数，量纲为t·hm-2·a-1；R为降雨侵蚀力因子，量纲为MJ·mm·hm-2·h-1·a-1；K 为土壤可蚀性因子，量纲为t·h·MJ-1·mm-1；LS 为坡长坡度因子，无量纲；C为植被覆盖与管理因子，无量纲；P为水土保持措施因子，无量纲。

2.2 RUSLE 模型因子计算

2.2.1 降雨侵蚀力因子(R)

降雨侵蚀力因子反映降雨对土壤的潜在剥蚀能力，是一项评估降雨引起的土壤分离和搬运的动力指标（李晓松等，2011）。目前R 的计算方法有很多种，本研究利用龙山县逐月降雨数据，计算研究区域降雨侵蚀力因子R值，计算公式如下：

式中：Pi为各月降雨量，单位mm；P 为年降雨量，单位mm。

2.2.2 土壤可蚀性因子(K)

土壤可蚀性因子反映土壤对侵蚀外营力剥蚀和搬运的敏感程度，是评价土壤遭受侵蚀敏感程度的指标（高峰等，2014；陈正发等，2021）。本研究将收集到的龙山县土壤类型分布图按照制定的坐标位置匹配后，把这些土壤进行空间矢量化分类，并按照我国土壤分类标准（GB/T 17296-2009）（中华人民共和国国家质量监督检验检疫总局和中国国家标准化管理委员会，2009）将其分为以下土壤亚类：沼泽土、水稻土、石灰（岩）土、黄壤土、紫色土、黄棕壤土、红壤土。参考前人研究（梁音和史学正，1999；江莉佳，2015），根据研究区不同位置土壤物质的组成成分，对土壤类型的K值进行赋值（表1）。

表1 研究区不同土壤类型的K值Table 1 K values of different soil types in the study area

2.2.3 坡长坡度因子(LS)

坡长坡度因子反映地形地貌特征对土壤侵蚀的影响，是降雨侵蚀动力的加速因子（邓辉等，2013）。坡度越大，土壤的重力势能越大，越容易被剥蚀；坡长越大，坡面水流沿程能量积累越大，土壤剥蚀量越大（付兴涛和张丽萍，2015）。本文利用12.5 m分辨率DEM，采用基于累积流量的单位汇水面积法计算LS 值（Wischmeir and Smith,1978），公式如下：

2.2.4 植被覆盖管理因子(C)

植被覆盖管理因子反映植被对土壤侵蚀的消减作用，一般来说植被覆盖度越大，地表植被越多，土壤在植物根系的固定作用下越不易被剥蚀。C值的取值范围为0～1，值越小，表示植被对土壤侵蚀的抑制作用越大。本文采用蔡崇法等（2000）的方法，根据NDMVI指数提取植被覆盖度，进而对C 值进行估算。

（1）计算研究区植被覆盖度FVC

式中：FVC 为植被覆盖度；NDMVI 为像元的山地植被指数；NDMVImin为全裸土地表的植被指数；NDMVImax为完全由植被覆盖地表的植被指数；取置信度为2%和98%所对应的NDMVI值。

使用2000 年11 月Landsat5 TM、2020 年11月的Landsat8 OLI 表观反射率数据计算各期的NDMVI 和NDVI，统计各期NDMVI 和NDVI 的最小值、最大值，分析植被指数变化的范围（表2），统计各期NDMVI 和NDVI 的频数、频率，分析植被指数各组数值频率频数变化（表3，图2）。

表2 研究区NDMVI与NDVI数值范围比较Table 2 Comparison of NDMVI and NDVI numerical ranges in the study area

表3 研究区NDMVI与NDVI数值频数频率比较Table 3 Comparison of NDMVI and NDVI numerical frequency comparison in the study area

图2 研究区2020年、2022年NDMVI与NDVI频率分布图Fig.2 Frequency distribution of NDMVI and NDVI of 2020 and 2022 in the study area

图3 基于NDVI和NDMVI提取的龙山县2020年C因子对比图Fig.3 Comparison chart of 2020 C-factor based on NDVI and NDMVI extraction in Longshan County

由表2、表3及图2可以看出，NDMVI的数值及频率变化范围较NDVI 更宽，在-1～1 范围内，2000年的NDMVI 数值范围较NDVI 增加了0.3158，2020年的NDMVI 数值范围较NDVI 增加了0.2076，增加幅度均较大。2000 年和2020 年的NDVI 在[-1，-0.7)和[0.7，1]区间的频数，均小于2000 年和2020 年的NDMVI。这说明NDMVI区分地物的能力更强，具有较强地消除复杂地形影响的能力。因此，使用该指数计算植被覆盖度和植被覆盖管理因子C，比使用NDVI计算更为精确。

（2）计算植被覆盖管理因子C

当植被覆盖度大于78.3%时，发生土壤侵蚀的概率很小，C 值取0；当植被覆盖度为0 时，最易发生土壤侵蚀，C 值取1（蔡崇法等，2000）。研究区植被覆盖管理因子2000 年平均值为0.59，标准差为0.31，2020年平均值为0.39，标准差为0.38。

为了验证NDMVI 提取C 因子的科学性，利用2020 年Landsat 影像分别提取NDVI 和NDMVI 指数，并计算对应的C因子。通过影像对比，可以看出NDMVI 区分地物的能力要优于NDVI，提取城镇用地、水体等地物的精度更高，尤其在地形起伏地区以及山坡的阴影地区，能更好地反演植被覆盖管理因子。

2.2.5 水土保持措施因子(P)

水土保持措施因子反映水土保持措施对土壤侵蚀的抑制作用（杨冉冉等，2013；陈红等，2021）。P 值的取值范围为0～1，值越小，表示水土保持能力越好，土壤越不易被侵蚀。已有研究中经常根据土地利用类型确定P值，本文参考已有的相关研究（江莉佳，2015），根据研究区土地利用/地表覆盖数据，确定研究区不同土地利用类型的P值（表4）。

表4 研究区不同土地利用类型的P值选择Table 4 P value selection of different land use types in the study area

3 结果与分析

3.1 土壤侵蚀结果

将两期RUSLE模型各因子计算结果进行叠加计算，得到龙山县土壤侵蚀模数分布图，龙山县2000 年平均土壤侵蚀模数为2116.18 t·km-2·a-1，2020 年平均土壤侵蚀模数为1275.84 t·km-2·a-1。根据我国《土壤侵蚀分类分级标准》(SL190-2007)（中华人民共和国水利部，2007），龙山县属于西南土石山区，容许土壤侵蚀模数为500 t·km-2·a-1，其2000年平均土壤侵蚀模数属于轻度接近中度水平。按照该标准分类，得到龙山县两期各等级土壤侵蚀面积如表5所示，各等级土壤侵蚀空间分布如图4所示。

表5 龙山县2000年、2020年土壤侵蚀情况Table 5 Soil erosion of Longshan County in 2000 and 2020

图4 龙山县2000年、2020年土壤侵蚀空间分布图Fig.4 Spatial distribution map of soil erosion of Longshan County in 2000 and 2020

3.2 土壤侵蚀空间特征分析

西南土石山区的土壤侵蚀允许量为500 t·km-2·a-1，据此，龙山县2000年的土壤侵蚀面积为1794.39 km2，占比57.75%。2020年的土壤侵蚀面积为1025.88 km2，占比33.02%，相比于2000 年减少了24.87%。其中，轻度、中度、强烈、极强烈侵蚀面积占比均比2000年有所下降。自1999 年以来，龙山县通过退耕还林、封山护林、坡改梯、实施节能工程等措施，治理水土流失（彭华，2006；江莉佳，2015；王海涛，2020），水土流失状况已得到很大改善。

从龙山县2000 年、2020 年土壤侵蚀空间分布（图4）来看，地形对土壤侵蚀具有显著的影响。坡度等级分级参照《土壤侵蚀分类分级标准》(SL190-2007)（中华人民共和国水利部，2007）的土壤水力侵蚀中对地面坡度分级标准，将龙山县水力侵蚀坡度分为[0°，5°]、(5°，8°]、(8°，15°]、(15°，25°]4 类。龙山县地形以山地为主，2000 年中度以上等级土壤侵蚀基本遍布整个区域，强烈及极强烈等级土壤侵蚀分布在水力侵蚀坡度在(8°，15°]的地区，剧烈侵蚀分布在水力侵蚀坡度在(15°，25°] 的地区。2020年中度、强烈及极强烈等级土壤侵蚀多发生坡度坡长较大的地区，即丘陵山壑地带，坡度为(8°，15°]。剧烈侵蚀土壤侵蚀区域在空间上呈条脉状分布，主要发生在山脚等陡坡附近，因为坡度较大，为(15°，25°]，加之降雨充沛，导致土壤侵蚀较为严重。

3.3 土壤侵蚀转移矩阵

引入转移矩阵（表6）对龙山县2000 年、2020年土壤侵蚀的演变状况进行分析，结合表5 进行分析可知：从2000 年到2020 年，研究区各等级间相互转换关系均相对显著；整个时段内，微度、轻度和剧烈侵蚀未发生改变的比例为89.98%、51.53%和91.87%，均具有较高的稳定性；除微度侵蚀外，各类别从高强度向低强度侵蚀转换的比例分别为47.32%、57.31%、56.49%、0%和59.36%，这说明从2000 年到2020 年研究区土壤侵蚀得到有效的治理和保护，整体向良好的方向发展。

表6 龙山县2000年、2020年土壤侵蚀强度转移矩阵Table 6 Soil erosion transfer matrix of Longshan County in 2000 and 2020

4 结论

本文尝试基于NDMVI山地植被指数反演植被覆盖管理因子C，利用改进RUSLE 模型对龙山县2000 年、2020 年土壤侵蚀的时空演变及发展规律进行分析，主要结论如下：

（1）利用2020 年Landsat 影像分别提取NDVI和NDMVI 指数并反演植被覆盖管理C 因子。通过影像对比，可以看出NDMVI 区分地物能力要优于NDVI，提取城镇用地、水体等地物的精度更高，尤其在地形起伏地区以及山坡的阴影地区，能更好地反演植被覆盖管理因子。验证了改进RUSLE 模型在西南土石山区的适用性，其可作为西南土石山区水土流失监测的技术手段。

（2）通过空间分析和数理统计，得出龙山县2000年平均土壤侵蚀模数为2116.18 t·km-2·a-1，土壤侵蚀面积为1794.39 km2；2020年平均土壤侵蚀模数为1275.84 t·km-2·a-1，土壤侵蚀面积为1025.88 km2。2020年相比2000年，龙山县各等级侵蚀面积均发生较明显的变化，微度侵蚀增长768.5 km2，轻度、中度、强烈、极强烈均呈不同程度的下降，除微度外，各类别高等级向低等级转换的比例分别为47.32%、57.31%、56.49%、0%和59.36%。龙山县土壤侵蚀得到有效的治理和保护，整体向良好的方向发展。