土壤侵蚀可致土地退化和生产力下降，对生态、经济和社会发展造成威胁，是全球性环境问题^[1]。近几十年，随着气候变化和人类活动升级，区域土壤侵蚀的形成机制愈发复杂，时空异质性不断提升^[2]。一方面，气候变化如降水增加，能够促发更多土壤侵蚀^[3]；另一方面，垦殖等土地资源开发，导致地表扰动，增大了土壤侵蚀发生的概率^[4]。土壤侵蚀量化及其驱动机制的深入研究，对区域水土保持和生态系统健康维系等，具有理论和实践方面的重要意义^[5-6]。

各类侵蚀模型是定量评估土壤流失的主流手段^[7]，常见的有基于物理的如WEPP，GUESS和EUROSEM等，也有基于经验的如HSPF，LASCAM和MMF等^[8-9]。大部分物理模型支持复杂参变量体系搭建，能很好地刻画土壤侵蚀过程，但率定和验证工作量大，较难应用于区域尺度；经验模型如通用土壤流失方程修订版(revised universal soil loss equation, RUSLE)，模型结构相对简单，数据获取方便，参数化本地化和模型率定容易实现，得到迅速推广，在侵蚀预报和水土保持方面发挥着重要作用^[10-11]，并被延伸应用至土壤侵蚀时空变化机制分析及基于地类的侵蚀强度分级和防控规划等领域^[12-13]。

甘肃省黄河流域覆被稀疏，下垫面脆弱；特别是陇中和陇东地区，降雨多呈“短历时，高强度”模式，容易激发土壤侵蚀，是甘肃省乃至全中国严重的水力侵蚀区^[14]。过去几十年，甘肃省黄河流域开展了大量水保工作，取得了显著成效^[15]；但由于降水总体稀少，人工植被成活率低，亦有很多沟道坝系等水保措施存在年久失修问题，水土流失防治的群体性、系统性和全局性仍面临挑战^[16]。当前，甘肃省黄河流域面临产业经济转型、乡村振兴和区域协调发展等历史新机遇，积极开展全域性土壤侵蚀量化评估，可为新时期水土流失防治提供理论储备，也是甘肃省水土保持发展战略实施的重要内容^[17]。但相关内容较少见于报道。

鉴于以上，本研究基于RUSLE模型量化评估甘肃省黄河流域土壤侵蚀时空分布及变化，探究不同下垫面情形下的土壤侵蚀规律，对比分析其驱动机制，为甘肃省黄河流域土壤侵蚀防控提供基础数据支撑和方法借鉴。

1 研究区概况

甘肃省黄河流域(33°6′—37°29′N，100°45′—108°43′E)系指甘肃省境内除内陆河流域(河西)和陇南长江流域以外的区域，大致位于秦岭以北，乌鞘岭以东，总面积14.59万km²，占甘肃省总面积的34.30%(图 1)，占黄河流域面积的19.00%。海拔介于739~4 865 m之间，地势西高东低；西部地处青藏高原东北缘，东部位于黄土高原西端，北部毗邻内蒙高原。按地理—生态空间格局及行政区划，甘肃省黄河流域可划分为甘南高原、陇中和陇东黄土高原等三大地貌类型^[18]。其中，甘南高原是黄河上游重要的水源涵养区，发育有黄河上游重要支流洮河、大夏河和长江二级支流白龙江等水系；陇中属于典型的黄土丘陵沟壑区，是甘肃省重要的经济文化核心区，发育有祖厉河和葫芦河等水系，源自青海的湟水、大通河和庄浪河等河流于该区汇入黄河干流；陇东属于典型的黄土台塬沟壑区，发育有泾河、蒲河和马莲河等水系，上述水系与陇中的葫芦河，同属黄河一级支流渭河水系；渭河于陕西潼关汇入黄河干流。陇中和陇东，是黄河干流重要的泥沙源区。

2 数据和方法 2.1 数据来源

高程数据采用中国科学院地理空间数据云平台(http://www.gscloud.cn)发布的数字高程模型(Digital Elevation Model，DEM)，空间分辨率30 m，主要用于水系勾勒和基础地形因子计算。降水数据(1980—2020)采用中国区域高时空分辨率地面气象要素驱动数据集(http://data.tpdc.ac.cn/)，该数据的适用性已有大量对比和验证^[19-20]，空间分辨率为0.1°。NDVI数据采用中国科学院地理科学与资源研究所发布的中国5 km逐月NDVI数据集(http://www.geodata.cn/)。土壤数据采用世界土壤数据库1∶100万土壤空间分布以及属性数据(http://www.ncdc.ac.cn/portal)。土地利用/覆被数据采用中国科学院地理科学与资源研究所发布的中国土地利用现状遥感监测数据(CNLUCC)数据(https://www.resdc.cn)，空间分辨率为1 km。为满足叠合分析，所有数据空间分辨率经重采样统一至1 km。

2.2 研究方法 2.2.1 RUSLE

RUSLE是通用土壤流失方程USLE的修订升级版，由美国农业部开发，被广泛用于土壤侵蚀和规划保护等领域的量化评估^[21]。该方程基于降水营力和下垫面因素的综合考量，以5类共6个因子乘积形式，估计单位面积土壤侵蚀量，公式为

$ A=R \times K \times L S \times C \times P $

(1)

式中：A为年土壤侵蚀模数, 即单位面积年土壤侵蚀量〔t/(hm²·a)〕; R为降雨侵蚀因子〔(MJ·mm)/(hm²·h·a)〕; K为土壤可蚀性因子〔(t·hm²·h)/(hm²·MJ·mm)〕; LS为地形因子, 其中L为坡长因子; S为坡度因子; C为植被覆盖和管理因子, 取值范围0~1; P为水土保持措施因子, 取值范围0~1; L, S, C和P均为无量纲因子。RUSLE模型应用较多, 各因子计算公式可参阅文献[22]。

2.2.2 Theil-Sen斜率

Theil-Sen斜率法是一种非参数统计的趋势计算方法，通过计算样本在不同长度下的变化率构造序列，取变率序列的中值反映时间序列的变化趋势及幅度。该方法的优点是能一定程度避免缺失或异常值对统计结果的影响^[23]。计算公式为：

$ \operatorname{Sen}_{i j}=\operatorname{median} \frac{\left(X_j-X_i\right)}{(j-i)} $

(2)

式中：Sen_ij为Theil-Sen斜率; X_i和X_j分别为第i和第j时刻的要素序列值, 1 < i < j < n; n为序列长度。

2.2.3 地理探测器

两个或多个地理要素之间存在机制性影响时，它们具有类似的空间分布；由此，可根据要素空间相关性，辨析它们之间的数据关联机制及影响程度，上述思想启发下，诞生了地理探测器模型^[24]。其中，因子探测器用于探测因变量的空间分异性，以及某个自变量对因变量的空间分布的解释能力，大小用q值来度量，公式为：

$ q=1-\frac{\sum\limits_{h=1}^L N_h \sigma_h^2}{N \sigma^2} $

(3)

式中：h=1, 2, 3…; L为变量的层次; N_h和σ_h²分别为h层样本的单位数和方差, N和σ²表示研究区的总单位数和方差; q为影响因子对土壤侵蚀空间差异的解释力, q∈[0, 1], q越大表示某因子的解释力越强。

交互探测器用于识别不同因子交互作用影响，即评价两个因子共同作用时的因变量解释能力；风险探测器用于判断两个子区域间的属性均值是否有显著的差别以识别风险区，用t统计量来检验(公式4)，如侵蚀模数的高值区^[25]。

$ S t=\frac{\bar{Y}_{h=1}-\bar{Y}_{h=2}}{\left[\frac{\operatorname{var}\left(Y_{h=1}\right)}{n_{h=1}}+\frac{\operatorname{var}\left(Y_{h=2}\right)}{n_{h=2}}\right]} $

(4)

式中：Y_h为子区域h内年土壤侵蚀模数A的均值; n_h为子区域h内样本数量; var表示方差。

3 结果分析 3.1 土壤侵蚀量化 3.1.1 RUSLE

基于多源遥感数据对RUSLE模型各因子进行本地化，结果见图 2。降雨侵蚀力因子R值介于313.60~2 248.24 (MJ·mm)/(hm²·h·a)之间，高值主要分布在甘南黄河玛曲段、洮河上游以及陇东，低值主要分布在陇中。土壤可蚀性因子K值介于0.01~0.04 (t·hm²·h)/(hm²·MJ·mm)之间，高值分布在陇中和陇东大部分地区，低值主要分布在陇东马莲河流域上游，甘南高原土壤可蚀性居中。坡度坡长因子LS值介于0.06~34.00之间，空间分布与地形相关，高值分布在坡度较陡及海拔较高分水岭地带，低值集中分布在西南部玛曲及陇中北部地势相对平坦地带。植被覆盖及管理因子C值介于0~1之间，高值主要分布在陇中北部植被覆盖较差地带，低值主要分布在甘南及陇中、陇东南部植被覆盖较好地带；水保措施因子P值介于0~0.85之间，高值主要分布于陇中、陇东农耕区，低值多位于甘南高原等植被生态良好区域。

3.1.2 时空变化

基于RUSLE模型的本地化，计算得到甘肃省黄河流域土壤侵蚀模数A，其多年均值介于0~215.00 t/(hm²·a)之间，全域平均为37.38 t/(hm²·a)。从分区来看，甘南、陇中北部、陇东北缘和东缘等处的土壤侵蚀较小，侵蚀模数一般在10 t/(hm²·a)以下；陇中中部A值较高，大致介于10~80 t/(hm²·a)之间；陇东大部土壤侵蚀剧烈，多地高于80 t/(hm²·a)(图 3a)。统计期间，甘南玛曲西北部、陇中西南部以及陇东南部和东部地区，土壤侵蚀有所减少，减率最大值为-3.54 〔t/(hm²·a)〕/a；陇中东南部、北部和陇东北部土壤侵蚀呈增加趋势，以马莲河流域上游增率最为显著，最大值达到1.86 〔t/(hm²·a)〕/a (图 3b)。甘肃省黄河流域及各区侵蚀模数A值的年际变化见图 4。总体来看，甘南A值最小，陇中居中，陇东最大，三区平均A值分别为11.45，38.13，63.44 t/(hm²·a)。1980—2020年，甘南A值微弱减小，减率为-0.02 〔t/(hm²·a)〕/a；陇中和陇东A值发生增加，增率分别为0.13和0.08 〔t/(hm²·a)〕/a；甘肃省黄河流域A值总体增加，增率0.12 〔t/(hm²·a)〕/a。

3.1.3 强度分级

根据水利部颁布实施的《土壤侵蚀强度分级标准》(SL190-2007)，按侵蚀模数＜5，5~25，25~50，50~80，80~150，＞150 t/hm²定义为微度、轻度、中度、强度、极强度和剧烈6个强度等级类别。基于RUSLE模拟结果，对甘肃省黄河流域土壤侵蚀强度进行评估(图 5)。结果表明，甘南高原以微度侵蚀为主，面积比例达到68%，轻度及以下面积比例85%，中度及以上面积的比例为15%，未出现剧烈侵蚀。陇中轻度及微度侵蚀面积的比例分别为29%和19%，中度及以上面积的比例为52%，其中包括11%的极强度侵蚀和1%的剧烈侵蚀区。陇东轻度及以下侵蚀面积的比例为33%，中度及以上侵蚀面积的比例67%，其中包括24%的极强度侵蚀和9%的剧烈侵蚀区。甘肃省黄河流域而言，微度侵蚀面积的比例为28%，轻度侵蚀面积比例为24%，中度及以上侵蚀面积的比例一般小于20%，其中，极强度侵蚀和剧烈侵蚀面积的比例分别为13%和3%，主要分布在陇中和陇东地区，尤以陇东为甚。

3.2 基于下垫面类型的A值统计 3.2.1 土地利用/覆被

统计表明，甘南、陇中和陇东均以农田A值最大，密林地A值最小，其余地类土壤侵蚀模数介于其间，从全域方面的土地利用A值表现为：农田>中覆盖草地>低覆盖草地>未利用地>灌丛>疏林地>建设用地>水域>高覆盖草地>密林地>沼泽湿地(图 6a—6d)，表明各区域及全域农田土壤侵蚀强度最大，密林地和沼泽湿地土壤侵蚀强度小。地类A值和对应面积的乘积得到侵蚀量，统计可知由于甘南草地分布较大，中覆盖草地土壤侵蚀量最大，其次为高覆盖草地、灌丛和农田，三者侵蚀量均超过4.00×10⁶ t/a (图 6e)；陇中地区农田和中低覆盖草地面积较大，因此3种土地利用类型土壤侵蚀量较高，而其余地类占地少，侵蚀量较小(图 6f)。陇东以农田和中覆草地的土壤侵蚀最为显著，两者侵蚀量之和占陇东总侵蚀量的88.46%(图 6g)。就甘肃省黄河流域而言，农田侵蚀量最大，中、低覆草地侵蚀量次之，其余地类侵蚀量相对较小。灌丛和高覆草地侵蚀量在1.00×10⁷~2.00×10⁷ t/a之间，密林地侵蚀量最小, 约为3.00×10⁶ t/a (图 6h)。

3.2.2 地形

高程、坡度和坡向等地形因子对土壤侵蚀的控制及影响，在一定范围内其规律较为相似^[26-27]。因甘肃省黄河流域幅员辽阔，分别以1 000 m，10°和45°为海拔、坡度和坡向等阶，进行地形因子区间划分和A值统计。结果表明，海拔1 000 m以下区域土壤侵蚀显著，A值较大，平均为59.11 t/(hm²·a)；海拔3 000~4 000 m区间侵蚀最少，A值较小，平均为6.73 t/(hm²·a)。总体来看，甘肃省黄河流域土壤侵蚀在低海拔区较大，在高海拔区较小(图 7a)。坡度来看，10°~20°区间侵蚀最显著，平均为46.21 t/(hm²·a)；大于40°区间侵蚀较少，A值平均为26.81 t/hm²·a)(图 7b)。坡向分布相对均匀，该因子对侵蚀的统计差异较小，总体观之，南坡地带A值较北坡略大，即阳坡侵蚀较阴坡显著(图 7c)。前述分析基于单因子统计得出，其中，A值在海拔4 000 m以上地带再度升高，随坡度增加呈先增后减态势，等等，说明不同立地条件下土壤侵蚀的发生发展，受多因子的共同控制和影响。

3.3 基于地理探测器的因子驱动分析

影响区域土壤侵蚀的因素众多，根据文献[27]和甘肃黄河流域实际情况，选择降水、气温、NDVI、坡度、高程、土地利用和土壤类型等7个因子，采用地理探测器就单因子和土壤侵蚀的空间关联进行区间对比分析。其中，土地利用和土壤类型按其分类体系进行统计，其余因子采用自然间断分级法，按7个数阶进行统计。

3.3.1 因子探测

因子探测器用于单因子驱动分析，结果表明(表 1)，坡度对甘南土壤侵蚀的解释力最强，q值为0.31；其次是降水和高程，q值均为0.20；气温、NDVI、土地利用和土壤等因子解释力较低，q值均小于0.1。土壤类型对陇中土壤侵蚀的解释力最强，q值为0.19；其次为坡度和降水，q值分别为0.17，0.13；其余因子q值均小于0.1。陇东而言，仍是坡度的解释力最强，q值为0.19；降水q值达到0.10，其余各因子的q值均小于0.1。上述情形表明，单因子层面，甘肃省黄河流域坡度、降水、土壤类型和高程的空间分异明显，对土壤侵蚀的发生发展较具解释力，相对而言，其余因子的单向解释力有限。

3.3.2 交互探测

交互探测器用于因子交互作用下的解释力分析。分析结果表明，任意两个因子交互对土壤侵蚀的解释力都大于单因子(表 2)，表明区域土壤侵蚀影响因素的多元性。甘南高原坡度和土地利用交互影响下的解释力最强，q值为0.41，说明不同坡度背景下，各土地利用方式的土壤侵蚀相差悬殊；此外，坡度与高程、NDVI、土壤及降水的交互解释力，均有明显提升。陇中地区土壤与坡度、气温及NDVI的交互解释力最强，q值均达到0.27；此外，土壤与土地利用的交互解释力也较强，q值为0.26。陇东与甘南相似，坡度与NDVI的交互解释力最强，q值为0.25；与其他因子的交互解释力也有所提升，q值接近或大于0.20。

3.3.3 风险探测

风险探测器用于甄别土壤侵蚀高发(风险)区的影响因子值域，3个地理分区对比来看，可以发现一些异同(表 3)。一般而言，三区降水较丰、坡度较陡和NDVI较小区间更容易发生土壤侵蚀；甘南和陇东的农田、陇中的低覆草地发生侵蚀的风险更大；水分相对亏缺的钙层土、干旱土和初育土更易发生侵蚀；甘南高原的侵蚀高风险区多位于高海拔较低温地带，陇中和陇东地区则位于较低海拔的中温地带。上述因子阈值在空间上并不重叠，有利于区域土壤侵蚀防控及水土流失治理的潜在靶区遴选。

4 讨论

RUSLE模型基于降雨、土壤、地形、植被和水保等因子构建，保持了与USLE相同的机制表达和模型结构，利于区域尺度逐像元雨力侵蚀预测。针对甘肃省黄河流域土壤侵蚀研究表明，该区土壤侵蚀高值区主要位于陇中和陇东，甘南总体表现为微度侵蚀，与魏健美等^[22]的研究结果基本一致。基于地类的统计结果表明，农田以及中、低覆盖草地的土壤侵蚀相对剧烈，林地、湿地的土壤侵蚀较为轻微。梁占岐等^[28]对黄土高原LUCC和土壤侵蚀的研究也表明，土壤侵蚀加剧的地区主要为农田和稀疏草地，而造成这一现象的原因是农田、稀疏草地对降水的蓄涵能力较弱，降水动力学机制促发土壤侵蚀的能力相对较强，导致潜在的流失风险。林地生态系统通过冠层截留、枯落物蓄持等方式削弱降水击溅和径流对土壤的剥离和搬运^[29]。资料显示，1980年代以来甘肃省黄河流域实施了多种形式较具规模的水土保持工作，土壤环境和陆表生态得以改良。本研究基于RUSLE得出的潜在土壤侵蚀呈增加态势，一定程度表明变化环境下甘肃省黄河流域的水土流失风险仍然存在，而降水增加是造成这一结果的主要原因^[30]。总体来看，RUSLE模型在甘肃省黄河流域的适用性良好。

气候变化和下垫面共同影响下，区域土壤侵蚀变化机理复杂。一方面，降水增加致使降雨侵蚀力因子增大，土壤侵蚀和水土流失的风险提升。甘肃省黄河流域降水与土壤侵蚀模数总体呈正相关，即为印证。在甘南高原，降水与土壤侵蚀的数据关联被检验为负相关，原因可能在于该区林地、湿地面积分布广，所具有强大遮蔽和蓄调功能对侵蚀产生抑减效应，集总至区域尺度，降水对土壤侵蚀的正向促增效应未得体现。气温对侵蚀的影响是间接的，增温会促增区域蒸散发，导致土壤水文干旱，从而加剧土壤侵蚀，这种情况在陇中、陇东较为常见。

与甘南不同，陇中和陇东属半干旱气候，大部分地域植被稀疏，气温与土壤侵蚀被统计为负相关，一定程度说明降水稀少地区，增温不能有效促增植被生长，其所导致的干旱风险，大概率会加剧土壤侵蚀。这一情形，与Sokasmnang K.E.O.等^[31]和张雪才等^[32]的研究成果一致。总体来看，土壤侵蚀的影响因素多元，

甘肃省黄河流域的范例研究表明，坡度、降水和覆被等因子与区域侵蚀之间的驱动关联较强，但因子交互能更好解释土壤侵蚀的空间分布，该结果有望对水土流失治理和土地资源优化等提供思路和方向。

5 结论

研究基于RUSLE模型的本地化，对甘肃省黄河流域土壤侵蚀的时空分布进行模拟，使用Theil-Sen斜率方法分析了土壤侵蚀的年际变化，采用空间分析方法探究不同地类和地形条件下的土壤侵蚀分布规律，并基于地理探测器甄别和比较了区域土壤侵蚀的主要影响因素。

(1) 1980—2020年，甘肃省黄河流域土壤侵蚀均值为37.38 t/(hm²·a)；甘南中度以上侵蚀面积占比15%，陇中为52%，陇东为67%；甘南以微度侵蚀为主，陇中和陇东以强度以上侵蚀为主；甘南在统计期间土壤侵蚀呈减小趋势，陇中和陇东呈增加趋势，陇中、陇东地区面临较为严峻的土壤侵蚀问题。

(2) 三区均以农田侵蚀模数最大，农田是土壤侵蚀防护重点域；受面积影响，甘南草地侵蚀总量最高，陇中和陇东农田侵蚀量最高，全域以农田和中、低覆草地的侵蚀最显著，草地是主要的土壤侵蚀域；地形因子分析表明，土壤侵蚀随高程先减后增，随坡度则先增后减，高值主要分布在海拔低于1 000 m坡度介于10°~20°区间，阳坡侵蚀总体高于阴坡。

(3) 坡度对甘南和陇东地区土壤侵蚀的空间分布具有显著影响，土壤质地对陇中地区侵蚀的影响更大；因子交互对土壤侵蚀的解释力一般大于单因子，有利于区域土壤侵蚀的空间探析；甘肃省黄河流域侵蚀高风险区主要位于丰水、陡坡和覆被稀疏地区。