土壤侵蚀是土壤及其他地表物质在风力、水力、重力、冻融等自然外力作用下，出现分离、破坏、剥蚀、搬运和沉积的现象，具有降低土壤肥沃性、破坏土地资源、减少陆地生态系统生产力、加剧旱涝灾害、影响全球碳循环等多重危害，严重威胁着生境质量、生态系统可持续发展、世界粮食安全、生物多样性和社会经济发展，已成为国内外学者研究土壤退化的重要课题^[1]。土壤侵蚀是长期且持续的变化过程，在环境中的强弱变化具有一定滞后性^[2]。研究流域内土壤侵蚀的时空演变通常需要较长时间序列的观测数据。目前，用于衡量土壤侵蚀的经验统计模型主要包括通用土壤流失方程^[3]和修正通用土壤流失方程(RUSLE)^[4]。其中RUSLE模型以较少数据和计算需求兼顾考虑土壤侵蚀的易感因素，在大尺度水平上的模拟结果受到国内外学者认可^[5-6]。

黄河流域生态质量对国家生态安全和经济社会发展都具有至关重要的意义^[7]。山西省位于黄河流域中游河段，是流经省区中土壤侵蚀发生最为密集和严重的省份之一，也是黄河泥沙含量的主要源区^[8]。脆弱的生态环境以及特殊的地质构造致使该区域土壤侵蚀成因复杂，自然因素的潜在影响和人类不合理活动的过度干扰诱导流域范围内土壤侵蚀程度进一步加剧^[9]。迄今为止，诸多学者对黄河流域开展土壤侵蚀发生机制^[10-11]、发展趋势及其驱动因素研究^[12-13]，但黄河流域间相应的自然要素、资源禀赋和发展程度差异显著^[7]，不同时空尺度之间难以进行结果的验证^[14]。探究土壤侵蚀时空分布格局以及判别侵蚀强度较高的区域对因地制宜地制定环境管理政策和土壤保持措施更具指导意义。

本研究以黄河流域山西段为研究对象，利用RUSLE模型和地理检测器对该流域2000、2005、2010、2015和2020年5期土壤侵蚀状况进行分析，旨在揭示不同时期黄河流域山西段土壤侵蚀强度的区域分异规律，量化单因子及多因子交互作用对土壤侵蚀强度的贡献权重，评估土壤侵蚀驱动因素的解释力，为决策者和规划者优化土地管理配置和实施生态保护战略提供科学依据和理论支撑。

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

黄河流域山西段(34°33′N ~ 40°18′N，110°11′E ~ 113°37′E)属于黄河中游地区，流域面积涵盖山西省11个市(地)86个县(市、区)，占全省国土总面积的62%。研究区属于暖温带大陆性季风气候，地貌属于典型的黄土广泛覆盖的山地高原。全区土壤包括粗骨土、栗褐土、褐土等多种类型，土质疏松，抗侵蚀能力差。

1.2 数据与方法 1.2.1 数据来源与数据预处理

本研究使用的基础数据主要包括气象数据、土壤数据、数字高程模型数据(DEM)、植被覆盖指数(NDVI)、植被型组空间分布数据、人口密度数据、坡度坡长因子数据和土地利用/覆盖类型数据(表 1)。以上原始数据均利用ArcGIS10.8中的数据管理工具(Data management)统一坐标系为WGS_1984_UTM_Zone_49N，使用数据转换工具(Data conversion)对数据格式进行均一化处理，重采样为100 m分辨率的栅格数据。2020年的NDVI数据参考2019年，受COVID-19的影响，2020年的NDVI值被假设为不变。

1.2.2 RUSLE模型

本研究综合考量影响土壤侵蚀的主要因素，凭借降雨侵蚀力、土壤可蚀性、坡长、坡度反映土壤侵蚀原动力，以植被覆盖度和水土保持措施度量地表抗力，采用RUSLE模型进行信息集成与处理，评估黄河流域山西段土壤侵蚀状况，计算公式如下：

$ A{\text{a}} = R \times K \times LS \times C \times P $

(1)

$ {A_{\text{p}}} = R \times K \times LS $

(2)

$ {A_{\text{C}}} = {A_{\text{P}}} - {A_{\text{a}}} $

(3)

式中：A_a，实际土壤侵蚀模数(t/(hm²·a))；A_p，潜在土壤侵蚀模数(t/(hm²·a))；A_c，土壤保持强度(t/(hm²·a))；R，降雨侵蚀力因子((MJ·mm)/(hm²·h·a))；K，土壤可蚀性因子((t·h)/MJ·mm))；L，坡长因子(无量纲)；S，坡度因子(无量纲)；C，植被覆盖因子(无量纲)；P，水土保持措施因子(无量纲)。

1) 降雨侵蚀力因子(R)选择Wischmeier和Smith^[15]的月尺度计算模型，计算公式如下：

$ R = \sum\limits_{i = 1}^{12} {1.735 \times {{10}^{\left[ {1.5\lg \left( {\frac{{{P_i}^2}}{{{\text{AP}}}}} \right) - 0.8188} \right]}}} $

(4)

式中：P_i，栅格水平上第i月(i=1, 2, ..., 12)平均降水量(mm)；AP，年降水量(mm)。

2) 土壤可蚀性因子(K)依据Erosion Productivity Impact Calculator模型中的K值估算土壤可蚀性^[16]。计算公式如下：

$ K = ( - 0.01383 + 0.51575K_{\text{EPIC}}) \times 0.1317 $

(5)

$ \begin{array}{l} K_{\text{EPIC}} = \left\{ {0.2 + 0.3{\text{exp}}\left[ { - 0.0256S_{\text{SAND}}\left( {1 - \frac{{S_{\text{SILT}}}}{{100}}} \right)} \right]} \right\} \hfill \\ \times {\left[ {\frac{{S_{\text{SILT}}}}{{S_{\text{CLAY}} + S_{\text{SILT}}}}} \right]^{0.3}} \times \left[ {1 - \frac{{0.25{\text{OC}}}}{{{\text{OC}} + \exp \left( {3.72 - 2.95{\text{OC}}} \right)}}} \right] \hfill \\ \times \left[ {1 - \frac{{0.7{\text{SN}}}}{{{\text{SN}} + \exp \left( { - 5.51 + 22.9{\text{SN}}} \right)}}} \right] \hfill \end{array} $

(6)

式中：K，修正后的土壤可蚀性因子；K_EPIC，修正前的土壤可蚀性因子；S_CLAY，黏粒(< 0.002 mm)百分比含量(%)；S_SILT，粉粒(0.002 ~ 0.005 mm)百分比含量(%)；S_SAND，砂粒(0.005 ~ 2 mm)百分比含量(%)；OC，有机碳百分比含量(%)；SN=1－S_SAND/100。

3) 地形因子(LS)：本研究基于泛第三极(20国)坡度坡长因子数据集^[17]提取黄河流域山西段LS因子数据集，其值的范围为0 ~ 72.12。

4) 植被覆盖因子(C)参考蔡崇法等^[18]提出的FVC与C值的对应关系式，计算公式如下：

$ {\text{FVC}} = \frac{{{\text{NDVI}} - {\text{NDVI}}_{\min}}}{{{\text{NDVI}}_{\max} - {\text{NDVI}}_{\min}}} $

(7)

$ C = \left\{ \begin{array}{l} 1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{FVC}} \leqslant 0.095 \hfill \\ 0.6508 - 0.3436\lg {\text{FVC }}0.095 < {\text{FVC}} \leqslant 0.783 \hfill \\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{FVC}} > 0.783 \hfill \end{array} \right. $

(8)

式中：FVC，植被覆盖度；NDVI，归一化植被指数；NDVI_max，完全植被覆盖的归一化植被指数；NDVI_min，无植被覆盖的归一化植被指数。

5) 水土保持措施因子(P)参照前人研究^[19]并结合研究区实际情况，根据土地利用/覆盖类型将耕地、草地、林地分别赋值为0.31、0.16、0.05，其余类型均赋值为1.00。

1.2.3 地理探测器模型

地理探测器由风险因子探测、交互作用探测、风险区探测、生态探测4个模块组成。因子探测一般用q值(公式(9))衡量其大小，q值越大，因变量的空间分异性越显著。交互作用探测可以分析两变量间的交互作用；生态探测可比较影响因子对土壤侵蚀空间分布的影响是否有显著差异；风险探测可判断影响因子的层间土壤侵蚀量是否有显著差别，并识别土壤侵蚀高风险区域^[20]。

$ q = 1 - \frac{{\sum\limits_{h = 1}^l {N_h {{\sigma_h}^2}} }}{{N{\sigma ^2}}} = 1 - \frac{{{\text{SSW}}}}{{{\text{SST}}}} $

(9)

式中：h=1, 2，...，l；l，因变量或自变量的分区；N_h和N，层h和全区的数据单元数；σ_h²和σ²，层h和全区因变量的方差；SSW，层内方差之和；SST，全区总方差。

鉴于自变量因素和土壤侵蚀强度(因变量)在空间统计尺度上存在显著差异，将海拔(X1)、坡度(X2)、植被覆盖度(X3)、年均降水量(X4)、植被类型(X5)、人口密度(X6)、土地利用类型(X7)和土壤侵蚀强度(Y)分别进行离散化处理，其中植被类型和土地利用/覆盖类型基于其数据类别直接离散化分级；植被覆盖度(FVC)根据 < 10%、10% ~ 30%、30% ~ 45%、45% ~ 60%、60% ~ 75%、> 75% 离散化为6级；坡度根据 < 5°、5° ~ 8°、8° ~ 15°、15° ~ 25°、25° ~ 35°、> 35°离散化为6级；年均降水量、海拔和人口密度选用自然断点法离散化为6级。

2 结果与分析 2.1 土壤侵蚀因子

本研究发现5个时间节点R值具有自西北向东南逐渐递增的特征，表明研究区域东南方位的降雨侵蚀力大于中部和西北地区(图 1)。此外，R值呈现出随海拔升高而增大的趋势。从全区平均来看，年际间的降水情况具有高度偶然性，2000年、2005年、2010年、2015年和2020年黄河流域山西段降雨侵蚀力因子R值分别为1 484.43、1 323.17、1 595.23、1 488.35和1 748.74(MJ·mm)/(hm²·h·a)，呈波动上升趋势，年均上升速率为69.38(MJ·mm)/(hm²·h·a)。

黄河流域山西段土壤可蚀性因子K的值域范围为0.004 ~ 0.021(t·h)/(MJ·mm)，在空间上呈现明显的斑块化分布特征(图 2)。全区土壤侵蚀风险低，高值主要出现在粉砂、黏土含量较低，沙砾占比较多的区域。

黄河流域山西段植被覆盖因子C的值域范围为0.0 ~ 1.0，5个时间节点植被覆盖因子空间格局基本一致，整体表现为西高东低、北高南低、中部地区波动的空间异质性特征(图 3)。C的峰值多出现在大同、忻定、太原、潞安、临汾、运城、晋城盆地的城市核心区，呈零星分布状态，这些区域城市化进程较快、生产要素密集、社会经济发展水平高，但存在人口聚集度高、工业污染严重等弊端，植被覆盖稀疏。从全区平均来看，2000—2020年间，黄河流域山西段植被覆盖因子C值总体呈现降低趋势，表明研究区域内植被覆盖情况有所改善，因覆被稀疏产生的土壤侵蚀呈减少态势。

2.2 土壤侵蚀强度时空变化特征

黄河流域山西段土壤侵蚀强度总共划分为微度、轻度、中度、强烈、极强烈和剧烈侵蚀6个侵蚀级别(图 4)，5个时间节点土壤侵蚀强度空间分布格局基本一致，均表现出东弱西强、南弱北强、流水地貌区弱黄土地貌区强的分异特征。

5个时间节点黄河流域山西段实际土壤侵蚀模数均值分别为31.29、25.07、29.32、22.12和25.67 t/(hm²·a)，呈波动显著下降趋势(P < 0.01)。结合表 2可知，研究时段内黄河流域山西段以微、轻度侵蚀为主，占总面积的57.05% 以上；中度侵蚀所占面积比例保持稳步下降态势；强烈、极强烈、剧烈侵蚀面积占比均呈现“减少-增加-减少-增加”的波动变化趋势。2000—2020年黄河流域山西段土壤侵蚀状况总体上得到明显控制，局部地区仍存在较高土壤侵蚀风险。

将实际土壤侵蚀模数与潜在土壤侵蚀模数进行空间叠置差值分析可以得到土壤保持强度^[21]。从图 5可以看出，2000—2020年黄河流域山西段土壤保持强度A_C值基本呈现出西部和东南部高、北部和中东部至西南部低、地形起伏且林草覆盖的丘陵和山地高、人类活动强度较高的盆地和平原低的空间分布特征，其中低值区多集中连片分布，中、高值区以交错分布为主。

2.3 土壤侵蚀格局的驱动因素分析 2.3.1 成因单因子探测

如图 6A所示，除人口密度以外，其他6个因子对土壤侵蚀强度的解释能力均具有显著影响(P < 0.01)。5个时段内各因子呈现出不同的变化趋势，植被覆盖度的解释力由2000年的0.123 2波动上升至2020年的0.143 4，海拔的解释力由0.024 7波动上升至0.046 4，两者20 a间上升趋势明显；土地利用/覆盖类型的解释力由2000年的0.066 5上升至2010年的0.087 1后出现下降；坡度的解释力变化最大，呈现先降后增又持续下降的明显下降趋势；年均降水量、植被类型的解释力总体轻微上升，q值均在2015年达到最大。

2000—2020年植被覆盖度和土地利用/覆盖类型是对土壤侵蚀空间分异性影响程度最大的两个因子，其中植被覆盖度的解释力最强($\bar q$=0.161 6)；坡度、年均降水量和植被类型的解释力位次变化显著；海拔的解释力位次呈轻微变化状态，在第3、4位次间交替往复；人口密度的解释力十分有限($\bar q$=0.002 2)，20 a间始终处于最末位(图 6B)。

2.3.2 成因交互作用探测

交互探测结果表明(表 3)，2000—2020年植被覆盖度∩坡度的交互协同作用表现为最强($\bar q$=0.220 0)，是影响土壤侵蚀空间分布的第1主导交互因子；植被覆盖度∩年均降水量($\bar q$=0.217 2)、植被覆盖度∩海拔($\bar q$=0.194 6)和植被覆盖度∩土地利用/覆盖类型($\bar q$=0.183 6)的交互影响力次之。

2000年坡度∩植被覆盖度为第1主导交互因子(q=0.183 3)，2005年坡度∩土地利用类型跃居为第1主导交互因子(q=0.208 7)，2010—2015年，坡度∩植被覆盖度又回升至第1位次(q=0.248 0、q=0.302 5)，2020年植被覆盖度∩年均降水量成为第1主导交互因子(q=0.200 9)。综上发现，在单因子探测中坡度(X2)的解释能力较为有限，但与植被覆盖度(X3)交互后表现最为突出，整体呈上升趋势且均为双因子增强，综合解释力较单因子解释力的显著增加表明20 a间黄河流域山西段土壤侵蚀的发生更大程度取决于双因子交互作用，驱动机理复杂。

2.3.3 风险区分析

由表 4可知，2000—2020年黄河流域山西段在海拔629 ~ 932 m、坡度5° ~ 8°、植被覆盖度30% ~ 45% 区间范围内出现土壤侵蚀现象的风险性最高，植被类型为栽培植被、土地利用/覆盖类型为草地或未利用地时土壤侵蚀最易发生，年均降水量和人口密度未表现出统一规律。

3 讨论

本文基于RUSLE模型定量评估了2000—2020年间黄河流域山西段土壤侵蚀的时空演变状况，20 a间实际土壤侵蚀模数范围为22.12 ~ 31.29 t/(hm²·a)，从2000年的31.29 t/(hm²·a)波动下降至2020年的25.67 t/(hm²·a)；空间上呈“局部加剧，整体改善”的变化趋势，全区土壤侵蚀强度逐渐降低，与张卓佩等^[19]关于黄河中游多沙粗沙区土壤水蚀时空变化特征研究较为一致。贾磊等^[22]使用RUSLE模型对黄河支流渭河流域的土壤侵蚀进行研究，发现2000年、2005年、2010年和2015年土壤侵蚀模数分别为50.72、42.27、39.72和35.94 t/(hm²·a)，总体表现为“减蚀降级”的变化特点；Yan等^[23]采用RUSLE和SEDD模型得到黄河二级支流北洛河流域的平均土壤侵蚀模数由1990年的181.90 t/(hm²·a)下降到2000年的74.09 t/(hm²·a)再降至2010年的28.58 t/(hm²·a)，土壤侵蚀程度明显缓解。黄河流域土壤侵蚀过程复杂、类型多样，即使同属于黄河流域，土壤侵蚀模数的估算结果也会因为研究时段、研究区域、研究尺度和研究方法的不同存在较大差异。山西省作为黄河流域重要省份之一，是保障国家能源安全和粮食安全的重要区域，也是拱卫京津冀和黄河生态安全的重要屏障，划定黄河流域山西段土壤侵蚀重点区域给予精准防治在推动黄河流域生态保护和高质量发展国家战略中具有举足轻重的地位。本研究发现黄河流域山西段57.05% 以上的土壤为微度或轻度侵蚀，对土壤侵蚀总量的贡献率不足20%。然而，半数以上的土壤侵蚀量源自仅占研究区面积20% 左右的强烈、极强烈和剧烈侵蚀地带，具有高强度侵蚀发生集中、危害程度严重的特点，应以黄土高原丘陵沟壑区和残塬沟壑区为重点治理区域，加快实施黄土高原水土流失综合治理、黄河重点生态区历史遗留矿山生态修复工程，全面推进“两山七河一流域”生态保护与修复。

黄河流域山西段土壤保持强度的空间格局特征明显，高值点主要位于以吕梁山为主干的黄土高原和太岳山、中条山、王屋山等块状山地，地表隆起形势强烈且植被覆盖度高，土壤保持能力较强，与《山西省黄河流域生态保护和高质量发展规划》中构筑以吕梁山区域、太岳山-中条山区域为主体的两大生态安全屏障范围基本一致。两大生态安全屏障是国家黄土高原-川滇生态屏障的重要构成，在保障京津冀生态安全中发挥重要作用；低值区则广泛分布在由大同、忻定、太原、临汾、运城等组成的“多”字形断陷盆地地区，基本呈东北-西南的条带状走向，景观格局相对破碎，土壤保持能力弱，土壤保持强度高低值分布与宁婷等^[24]和吴树荣等^[25]的研究结果较为接近。

由于土壤侵蚀会受到人地相互作用影响^[26]，揭示自然环境因素(坡度陡度、降水侵蚀力、土壤易蚀性)和人类活动因素(人口密度、地表覆被类型、植被覆盖度和土地利用/覆盖类型)与土壤侵蚀的关系对后续采取差别化的管理措施以提高防治效率尤为重要。研究区域内地形以山地丘陵为主，广泛发育黄土地貌，土壤纬度地带性鲜明，以半淋溶土和钙层土为主要类型，质地粗、砂性大致使土体松弛、自身持水性差、土壤侵蚀的成因和运动方式复杂多变。本研究通过地理探测器模型定量刻画各要素及两两要素间的交互作用对土壤侵蚀强度分布格局的贡献量大小^[27]，发现黄河流域山西段土壤侵蚀强度格局主要是坡度、植被覆盖度、年均降水量和土地利用/覆盖类型综合作用的结果，植被覆盖度始终是侵蚀发生的首要驱动因素，与尹作堂和常军^[26]、贾磊等^[22]和张艳等^[28]得出植被覆盖度在控制土壤侵蚀的各因素中占据主导地位的结论一致。植被覆盖度高的下垫面地下根量大、根系发达，拦蓄径流、抗冲抗蚀性能相应提升，有利于土壤保持^[29]；覆盖度低的生境承载力差，生态系统结构功能退化，地表植被的水源涵养能力有限，超过土壤蓄渗限度的降水会形成细沟侵蚀，严重时还会出现下切和溯源侵蚀^[30]。

4 结论

本研究以黄河流域山西段为例，基于2000—2020年多源数据通过RUSLE模型和GIS空间分析技术揭示土壤侵蚀强度的时空分布格局及变化情况，剖析侵蚀发生高风险区域所在的区位特征，利用地理探测器模型量化影响土壤侵蚀强度的主导因素及交互影响程度。发现黄河流域山西段2000—2020年实际土壤侵蚀模数的空间分异特征显著，侵蚀主要发生在研究区西部的黄土高原丘陵沟壑区和中、北部煤矿集中分布区域。近20 a该流域大部分地区土壤侵蚀强度向低强度侵蚀转化，但因局部水土流失依然严重，全区常年处于中度侵蚀级别；对各因子进行土壤侵蚀风险探测发现：栽培植被是土壤侵蚀的主要植被类型；草地和未利用地较其他地类而言更易发生土壤侵蚀；5° ~ 8°坡度等级时实际侵蚀量出现最大值，土壤保持能力较弱；海拔为629 ~ 932 m的区域需注重土壤侵蚀的防治工作；植被覆盖度在30% ~ 45% 区间范围内土壤侵蚀强度最大，应作为重点关注的梯度带；影响土壤侵蚀的众多因子中，植被覆盖度和土地利用/覆盖类型始终是引起土壤侵蚀变化的主导因素；对因子交互作用分析可知，因子间的交互作用大于单因子的影响，其他因子与植被覆盖度交互后对侵蚀强度的贡献量明显增大，20 a间黄河流域山西段土壤侵蚀格局演变经历由坡度与植被覆盖度、坡度与土地利用类型或植被覆盖度与年均降水量交替主导的过程。