水土流失不仅破坏土壤肥力，导致土地生产力的下降，威胁国家粮食生产安全^[1]，而且携带了大量的泥沙淤积河道，造成严重的危害^[2]。土壤保持是指生态系统对土壤侵蚀的防控能力以及对泥沙存储的保持能力^[3]，在维护区域侵蚀控制和生态安全方面起到了独特的作用，为区域降低灾害提供了重要保障。1965年，Wischmeier和Smith^[4]首次提出了通用水土流失方程（Universal Soil Loss Equation，USLE），对该方程各个因子的含义和算法进行不断的修正和调整，最终形成了修订后的通用土壤流失方程（RUSLE）^[5]。2002年Liu等^[6]提出了更适合中国区域的中国土壤流失方程（CSLE）。随着学者们对土壤侵蚀模型研究的不断完善，生态系统服务和交易的综合评估模型（Integrated Valuation of Ecosystem Services and Trade-offs，InVEST）应运而生，此模型弥补了USLE和RUSLE方程缺少地块拦截上游泥沙模块这一缺陷，是量化大区域/国家层面上土壤侵蚀较为合适的工具。近年来，InVEST应用较为广泛。基于500 m空间分辨率的土壤数据，Marques等^[7]研究了葡萄牙土壤侵蚀的时空分布。Xiao等^[8]采用空间分辨率为90 m的DEM数据和1 km的NDVI数据对三峡库区土壤保持功能进行了分析，并探讨了影响土壤保持功能的因素。Guo等^[9]分析了淮河地区土壤侵蚀和土壤保持的时空变化，并通过随机森林和地理探测器对驱动因素进行评价和量化。赵欣悦等^[10]采用空间分辨率为90 m的DEM数据对北三河流域的土壤保持情况进行研究。贾婉琳等^[11]基于空间分辨率为1 km的土壤数据分析赤水河流域的土壤保持情况。总体而言，国内外有关区域尺度土壤侵蚀的研究采用的栅格数据空间分辨率相对较粗。1 km的世界土壤数据库（HWSD），1 km的NDVI数据和90 m DEM数据。然而，无论是基于何种方程，在大尺度上获取准确的土壤侵蚀参数对研究结果均至关重要^[12]。玉院和王金亮等^[13]通过分析不同空间分辨率（30 m、60 m、90 m、120 m）下的地形因子的差异性，表明在土壤侵蚀等模型中对地形因子进行估算时，应该充分考虑空间分辨率的影响。同样，土壤侵蚀格局的精细表征也需要高分辨率土壤可蚀性因子数据支撑^[14]。提高研究区土壤侵蚀模型输入数据的空间分辨率可以有效地提高该区域土壤侵蚀结果精度。

川滇区是我国重要的生态功能区，对构建国家整体生态安全格局的战略地位十分突出。然而长期以来，川滇区生态环境“家底”不清，生态环境变化的系统性和完整性调查工作尚属空白。现有的土壤侵蚀或土壤保持数据集产品的空间分辨率较低（1 km），且时效性较差（2000—2015年）。数据源空间分辨率的高低一定程度上会影响土壤侵蚀结果的分析，研究的时间范围也会对土壤侵蚀趋势有一定的影响。因此，本文采用空间分辨率为90 m的土壤数据、30 m的DEM数据和NDVI数据对川滇生态屏障区2000年、2010年、2015年、2020年的土壤侵蚀进行研究，并运用地理探测器对土壤侵蚀的主控因子进行分析。研究结果将为川滇区的水土保持综合治理和生态环境保护政策提供理论依据和科学决策参考。

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

川滇生态屏障区位于川滇西部及四川盆地的过渡地带（图 1），在24°78′~33°73′N，98°72′~106°77′E之间，横跨我国长江、黄河、澜沧江和怒江等四大江河，地跨四川、云南两省，由川西73个县（区）和滇北17个县（区）构成，幅员面积2.37×10⁵ km²。位于青藏高原东缘地形急变带的川滇区，地形西北高、东南低，平均海拔约4 000 m。以高原山地温带、寒温带季风气候为主，属于高原山地垂直体系。区域内地震、泥石流、滑坡等山地灾害类型多样、分布广泛、活动频繁，并具有规模巨大、灾害链效应明显、破坏力强、容易形成巨灾等鲜明特点^[15]。

1.2 数据来源

本文所需数据包括：（1）川滇区边界和行政区划的矢量数据，源于中国科学院资源环境科学与数据中心；（2）DEM数据，采用空间分辨率为30 m的SRTM1 V3.0的中国区域数据；（3）土壤数据，数据来源于国家科技资源共享服务平台—国家地球系统科学数据中心-土壤分中心（http://soil.geodata.cn）^[16-17]；（4）逐月降雨数据集，采用国家青藏高原科学数据中心，空间分辨率为1 km的逐月降水量数据集^[18]；（5）NDVI数据，源于中国科学院地理科学与资源研究所土地利用与全球变化遥感团队制作的中国30 m年最大NDVI数据集^[19]；（6）土地利用数据，采用了在GEE平台上制作的源自Landsat的年度中国土地覆盖数据集^[20]，空间分辨率为30 m。上述数据中，逐月降雨数据集、NDVI数据和土地利用数据的时间范围为4期：2000年、2010年、2015年和2020年。

1.3 InVEST模型原理

InVEST模型以GIS作为应用平台，以空间数据为参考基础，通过输入较少的数据就可以输出较多的数据，打破了传统方法的局限性，使复杂问题简易化，运行结果可以在ArcGIS上以地图的方式呈现出来，使结果可视化^[21]，较纯文字的表达方式更加直观^[22]。目前，此模型已非常成熟，适用于全球范围各尺度上的生态系统服务功能评估，推广性强。其中土壤保持模块的具体公式如下：

${\text{USL}}{{\text{E}}_i} = {R_i} \times {K_i} \times {\rm{L}}{{\rm{S}}_i} \times {{\rm{C}}_i} \times {{\rm{P}}_i}$

(1)

$RKL{S_i} = {R_i} \times {K_i} \times L{S_i}$

(2)

${\rm{SD}}{{\rm{R}}_i} = {\rm{RKL}}{{\rm{S}}_i} - {\rm{USL}}{{\rm{E}}_i}$

(3)

式中，${\rm{RKL}}{{\rm{S}}_i}$为第$i$个栅格单元的年潜在土壤侵蚀量，${\text{USL}}{{\text{E}}_i}$为第$i$个栅格单元的年实际土壤侵蚀量，${\rm{SD}}{{\rm{R}}_i}$为第$i$个栅格单元的土壤保持量，年潜在土壤侵蚀量与年实际土壤侵量之差为年土壤保持量。${R_i}$、${K_i}$、${\rm{L}}{{\rm{S}}_i}$、${{\rm{C}}_i}$、${{\rm{P}}_i}$分别为第$i$个栅格单元的降雨侵蚀力因子、土壤可蚀性因子、坡长因子、地表覆盖管理因子、土壤保持管理因子。

1.4 土壤侵蚀强度分级

土壤侵蚀强度分级标准采用中华人民共和国水利部颁布的《土壤侵蚀分类分级标准》（SL 190—2007）。以年平均侵蚀模数为判别指标，具体分类分级标准如表 1所示。

1.5 模型参数处理

（1）降雨侵蚀力因子R，采用四期月降雨量的栅格数据，计算公式具体见参考文献[23]。（2）土壤可蚀性因子K，采用Williams等^[24]提出的EPIC模型，模型中的粒径包括砂粒、粉粒、黏粒，具体公式见参考文献[25]。（3）地形因子LS，可由InVEST模型根据填洼后的DEM数据自动计算出地形因子LS^[26]。（4）地表覆盖管理因子C，采用目前适用较广的蔡崇法模型^[27]。（5）土壤保持措施因子P，本文参考郎燕等^[23]对川滇区土壤保持措施因子的P取值，耕地为0.4，林地、草地、水域、建设用地、未利用地和其他用地为1。

1.6 土壤侵蚀主控因子识别

采用地理探测器进行土壤侵蚀主控因子的识别，其为探测空间异质性的工具，该方法的原理是如果某个自变量对某个因变量有重要影响，那么自变量和因变量的空间分布应该具有相似性^[28]。本文使用地理探测器的“分异及因子探测”和“交互作用探测”模块对研究区土壤保持服务的主要影响因子和交互作用进行定量探测。将影响因子坡度、高程、植被覆盖度、土地利用类型、降雨侵蚀力因子、土壤可蚀性因子定义为自变量，并对其进行编码X1、X2、X3、X4、X5、X6，土壤侵蚀数据为因变量Y。除土地利用类型数据外，输入的其他数据均为连续变量，本文按照自然断点法对自变量进行重分类，得到类别数据。再根据研究区矢量数据生成10 km×10 km的规则网格，网格的中心点为采样点，按照采样点对X、Y进行提取分析，去除异常值后将数据导入地理探测器模型中，分析得到各因子对土壤侵蚀的影响力以及双因子对土壤侵蚀的交互作用。

1.7 数据处理与绘图

在InVEST模型中运算土壤侵蚀和土壤保持数据，模型所需参数R、K、LS、C、P（空间分辨率为30 m）以及后续的分析、出图均采用ArcGIS 10.8.1软件进行处理；在GeoDetector 2018软件中实现土壤侵蚀主控因子的识别；数据统计和绘图采用Microsoft Excel 2021和Origin 2018软件。

2 结果 2.1 土地利用变化特征

图 2为川滇区不同时期的土地利用情况和空间分布图，分析可得川滇区土地利用类型以林地为主，约占62%，所占面积呈增长趋势；其次为耕地和草地，分别约占17.5%和16.6%，所占面积呈减少趋势；建设用地面积约占0.8%，且呈增长趋势，其他土地利用类型的面积占比较少，且四期无明显变化。耕地和建设用地主要分布在川滇区的东北地区，其中耕地面积较大，呈片状分布；建设用地呈零星线状由中心城区向四周散开，如图 2中的成都市。川滇区的西部和南部多为山地，土地利用方式以林地为主。

2.2 降雨侵蚀力变化特征

根据公式计算得到四期降雨侵蚀力因子R栅格图层，如图 3所示。川滇区降雨侵蚀力因子R总体上呈现从西北至东南递增的趋势，从2000年至2020年，降雨侵蚀力的空间分布没有明显差异，降雨侵蚀力较高的区域主要集中在邛崃市、蒲江县、彭山区、东坡区、丹棱县、名山区、雨城区、洪雅县等区域。2000年至2020年四期的降雨侵蚀力平均值分别为2.7×10⁵、2.3×10⁵、2.6×10⁵、3.0×10⁵ MJ⋅mm⋅km^–2⋅hm^–2⋅a^–1，最低值和最高值均出现在2015年，分别为1.1×10⁵、9.6×10⁵ MJ⋅mm⋅km^–2⋅hm^–2⋅a^–1。

2.3 土壤侵蚀

根据InVEST模型输出结果，在ArcGIS中对其进行分析，得到2000年、2010年、2015年、2020年川滇区土壤侵蚀情况。总体上看，川滇区土壤侵蚀呈先减少后增长的趋势，土壤侵蚀模数从2000年的144.8 t⋅km^–2先降低至2010年的137.5 t⋅km^–2，然后增长至2020年的231.4 t⋅km^–2。土壤侵蚀总量从2000年的3.4×10⁷ t先降低至2010年的3.3×10⁷ t，然后增长至2020年的5.5×10⁷ t，与2000年相比，2020年增加了2.1×10⁷ t。对川滇区土壤侵蚀强度进行划分，得到图 4（空间分辨率为30 m）。

川滇区土壤侵蚀强度以微度侵蚀为主，主要分布在川滇区东部地势较为平坦的低海拔区域，2000年占据总面积的86.1%，2010年后微度侵蚀面积开始有所减少，2015年和2020年微度侵蚀所占比例分别为85.9%和81.6%。其中变化最为明显的区域为甘洛县、峨边彝族自治县、马边彝族自治县、美姑县和越西县，由2000年少量区域为轻度侵蚀发展为2020年大部分区域为轻度侵蚀。轻度及以上侵蚀强度多分布于川滇区的西部地区。如康定市，该区域以轻度侵蚀为主，局部区域出现剧烈侵蚀，与泸定县接壤处剧烈侵蚀情况相对明显，整个川滇区的剧烈侵蚀大多分布于此。随着年份的增长，轻度及以上侵蚀强度所占面积也在增长，轻度侵蚀所占面积由2000年的2.7×10⁴ km²增长为2020年的3.5×10⁴ km²，中度侵蚀所占面积由2000年的3.6×10³ km²增长为2020年的5.5×10³ km²。与2000年相比，2020年强度侵蚀、极强度侵蚀和剧烈侵蚀所占面积均有不同程度的增加，分别增长了649.2、363.7、185.2 km²，其中汶川县和理县、小金县的接壤处变化最为明显。

2.4 土壤保持功能

通过InVEST模型计算得到川滇区土壤保持情况，结果如图 5所示。土壤保持呈先下降后增长的趋势，从2000年至2010年土壤保持模数开始降低，从4.0×10³ t⋅km^–2降至3.5×10³ t⋅km^–2。2010年至2015年土壤保持模数仍呈下降趋势，但仅降低了38.6 t⋅km^–2，趋于好转。从2015年至2020年，土壤保持模数呈增长趋势，从3.5×10³ t⋅km^–2增长至4.5×10³ t⋅km^–2。土壤保持总量的变化趋势与土壤保持模数保持一致，2000年至2015年为下降趋势，2015年至2020年为增长趋势。由图 5可知，土壤保持强度低于200 t⋅km^–2的区域多为海拔较高的山顶、水域和海拔较低的居民聚集地，如位于康定市和泸定县接壤处的燕子沟和海螺沟冰川森林公园；途经石棉县、汉源县、甘洛县的大渡河流域；川滇区东北部的东坡区、新津县、广汉市等地。从2000年至2020年，土壤保持强度大于5 000 t⋅km^–2的面积有一定的增加，中部和北部地区变化相对明显。

图 6a为不同土地利用类型下的土壤保持模数，其中冰川和裸地的土壤保持模数太小可忽略不记；从2000年至2020年耕地、林地、灌木、草地、建设用地和湿地的土壤保持模数呈先减少后增多的趋势，节点出现在2010年；水域的土壤保持模数逐年升高。林地和灌木的土壤保持模数远远大于其他土地利用类型，且二者很接近，其次为草地、耕地、水域、湿地、建设用地、裸地和冰川。图 6b为不同土地利用类型下的土壤保持总量，水域、建设用地、湿地、冰川和裸地的土壤保持总量相对较少，林地的土壤保持总量占主要部分，占比约为75%以上，2000年至2020年4期的土壤保持总量分别为1.0×10⁹、9.1×10⁸、9.9×10⁸、1.2×10⁹ t，与对应的土壤侵蚀模数波动趋势一致；其次为耕地、草地和灌木。

2.5 土壤侵蚀主控因子

通过分析土壤侵蚀各影响因子的解释力q值，发现2000年、2010年、2015年和2020年的主控因素均为土地利用类型，其解释力q值逐年降低，由2000年的0.2418降至为2020年的0.1850。高程、植被覆盖度、土壤可蚀性、降雨侵蚀力因子和坡度的解释力q值依次降低，对土壤侵蚀的影响依次减弱。采用地理探测的交互作用探测模块对双因子的交互作用进行定量探测，可知任何两种变量对土壤侵蚀空间分布的交互作用均要大于第一种变量的独自作用。土地利用类型与其余5个影响因子的交互协同作用对解释力q值提升均较为明显，其中土地利用类型与土壤可蚀性因子的交互协同作用对土壤侵蚀的影响最大，2000年、2010年、2015年和2020年的q值分别为0.51、0.51、0.38和0.5。

3 讨论

位于青藏高原东缘地形急变带的川滇区，是我国乃至全球范围内地质灾害最为严重的区域之一，正因如此，该区的水土保持情况也备受关注。在以往的研究中，有学者采用CSLE和RUSLE方程对川滇区土壤侵蚀和土壤保持情况做了研究^[29-30]。在此基础上本研究采用InVEST模型进一步对川滇区土壤保持情况进行分析，结果表明，从2000年至2020年川滇区土壤侵蚀模数和土壤保持模数均呈先降低后升高的趋势。其中2015年至2020年侵蚀总量呈增大趋势，与徐庆等^[29]对川滇区2015年至2016年的研究结果基本一致。但与陈江玲^[30]的研究结果有一定的差别，主要原因是在研究期限、研究方法和研究区域等方面存在差异。陈江铃采用CSLE方程，研究了滇北4个县和川西27个县2010年和2015年的土壤侵蚀情况，结果表明，川滇区从2010年至2015年土壤侵蚀模数有减少的趋势，但侵蚀面积仍然有扩展的趋势。与前人研究相比，本研究生成的数据产品的空间分辨率更高（30 m）、结果更直观、可视化更强。虽然川滇区以山地为主，地形陡峭、海拔较高、自然灾害较频繁，导致其潜在土壤侵蚀较大，但是该区的植被覆盖较广，土地利用类型以林地和草地为主。而该用地类型对土壤保持的能力最强，效果最好，这与曹泰彰^[31]的研究结果保持一致。土壤保持能力较低的区域多为地势平坦的居民聚集地，这与张佳惠^[32]的研究结果保持一致，住宅区一般都是开发较为完善的区域，地表缺少自然保护，土壤保持能力偏低。土壤侵蚀和土壤保持不仅与自然因素密切相关，人为因素对其同样有很大的影响，本研究采用地理探测器对其驱动因子进行分析，结果表明川滇区土壤侵蚀的主控因子为土地利用类型，与黄鑫等^[33]对沂蒙山区土壤侵蚀影响因素的研究结果一致。

本研究还存在以下不足：（1）InVEST模型打破了传统方法的局限性，使复杂问题简易化，将复杂机理模型拆分成五个因子，同时也导致了模型输出结果对模型输入参数非常敏感^[34]。因此，在以后的研究中，考虑对各输入参数进行敏感性分析和精度验证。（2）本研究使用的是空间分辨率为90 m的土壤数据、30 m的DEM数据和NDVI数据。今后研究可以考虑进一步提高基础数据的空间分辨率，例如，采用空间分辨率为12.5 m的DEM数据，土地利用类型数据和NDVI数据可以通过哨兵2数据（空间分辨率为10 m）来加工生成，土壤数据可以通过数字制图方法来生产空间分辨率为10 m的研究区基本土壤属性数据。

4 结论

本文采用InVEST模型对川滇区2000年、2010年、2015年和2020年的土壤侵蚀情况进行分析，对不同土地利用类型、不同植被覆盖度、不同坡度以及不同海拔的土壤保持功能进行了分析，并运用地理探测器对其土壤侵蚀主控因子进行研究，主要结论如下：（1）土壤侵蚀呈先减少后增长的趋势，土壤侵蚀模数和土壤侵蚀总量分别从2000年的144.8 t⋅km^–2和3.4×10⁷ t增长至2020年的231.4 t⋅km^–2和5.5×10⁷ t，侵蚀强度以微度侵蚀为主。（2）土壤保持呈先下降后增长的趋势，2000年至2015年土壤保持模数和土壤保持总量呈下降趋势，从2015年开始逐渐上升，2015年出现最低值，分别为3.5×10³ t⋅km^–2和8.2×10⁸ t。（3）各影响因子对土壤侵蚀的影响程度从高到低依次为：土地利用类型、高程、植被覆盖度、土壤可蚀性因子、降雨侵蚀力因子、坡度，任何两个影响因子的交互作用对土壤侵蚀的影响均大于单因子的影响。