基于PESERA的三江源土壤侵蚀时空分布及其驱动因素

Spatiotemporal Pattern of Soil Erosion Rates and Its Driving Factors in the Three-River Headwater Region Based on the PESERA Model

摘要

目的 为深入剖析三江源地区土壤侵蚀空间分布和主导因素，保障区域生态安全。 方法 利用区域尺度土壤侵蚀过程模型PESERA模拟1983—2022年4个时段（每10 a为1个时段）三江源地区的平均侵蚀速率，并采用最优地理探测器分析单因子对侵蚀速率的解释力及多因子间的交互作用。 结果 PESERA模型在三江源地区验证结果较好（NSE=0.53，R²=0.69）。三江源地区的侵蚀速率呈由南向北逐渐降低的空间分布特征，剧烈侵蚀［>150 t∕（hm²·a）］、轻度侵蚀［<5 t∕（hm²·a）］和微度侵蚀［5~25 t∕（hm²·a）］分别集中在研究区南部的澜沧江源地区、西北部的长江源地区和东部黄河源地区，其余侵蚀强度主要在三江源南部地区。单因子探测结果表明，地形起伏度和降水量对侵蚀速率的解释力最高，分别为0.43和0.35。在交互分析中，降水量与植被覆盖度的交互作用解释力最高为0.73，地形起伏度与其他因子的共同作用显著提高对侵蚀速率的解释力。 结论 1983—2022年三江源土壤侵蚀以轻度侵蚀为主，侵蚀速率产生先上升后下降的时间变化，地形起伏度、降水量、植被覆盖度为影响该区域侵蚀的主要因素。研究结果为三江源地区的土壤侵蚀防治和生态保护提供有益参考。

Abstract

Objective This study aimed to comprehensively analyze the spatial distribution and dominant influencing factors of soil erosion in the Three-River Headwater Region for ensuring regional ecological security. Methods The average erosion rate for four decadal periods from 1983 to 2022 was simulated using a regional-scale process-based erosion model， PESERA. The optimal Parameters-Based Geographical Detector （OPGD） method was employed to evaluate the explanatory power of individual factors on erosion rates and the interactions between multiple factors. Results The accuracy of PESERA modelling results was acceptable in the Three-River Headwater region （NSE = 0.53， R² = 0.69）. The spatial distribution of erosion rates showed a gradual decrease from south to north. Severe erosion ［> 150 t∕（hm²·a）］， mild erosion ［< 5 t∕（hm²·a）］， and slight erosion ［5–25 t∕（hm²·a）］ were mainly concentrated in the headwater area of Lancang River （southern part of the study region）， the headwater area of Yangtze River （northwestern part of the study region）， and the headwater area of the Yellow River （eastern part of the study region）， respectively. Other erosion intensities were mainly distributed across the southern parts of the Three-River Headwater Region. The analysis of individual factor detection revealed that terrain relief （0.43） and precipitation （0.35） were the most influential factors. Interaction analysis indicated that the combination of precipitation and vegetation cover had the highest explanatory power （0.73）， while the interaction between terrain relief and other factors considerably enhanced their effects on erosion rates. Conclusions Between 1983 and 2022， mild erosion was the dominant type of soil erosion in the Three-River Headwater Region. Temporally， erosion rates exhibited an initial increasing trend followed by a subsequent decline. The main factors affecting soil erosion in the region were terrain relief， precipitation， and vegetation cover. The findings can provide valuable insights for soil conservation and ecological protection in the Three-River Headwater Region.

Graphical abstract

关键词

三江源 / 土壤侵蚀 / PESERA模型 / 最优地理探测器 / 驱动因素

Key words

Three-River Headwater Region / soil erosion / PESERA model / OPGD / driving factors

土壤侵蚀为全球性环境问题，严重威胁土壤质量、农业生产力和生态系统稳定性，增加泥石流、洪水等自然灾害风险^［1］。然而，传统的土壤侵蚀研究方法如同位素示踪、水文站点和径流小区观测、野外试验虽然可以提供精确的数据，但主要集中在坡面尺度和小流域尺度，对于区域尺度的侵蚀量化和预测等方面存在不足^［2］。相比之下，模型模拟通过考虑气候、地形、植被等多种影响因素，能够实现跨尺度侵蚀过程的动态模拟和预测，尤其在区域尺度的侵蚀模拟和风险评估中展现出明显优势^［3］。因此，模型模拟方法成为解决区域尺度土壤侵蚀研究的重要工具，为实现全面的侵蚀防治和土地管理提供坚实的基础。

区域尺度模拟中，目前使用较为广泛的模型为Revised Universal Soil Loss Equation（RUSLE）等经验模型，但其以经验公式为基础，依赖于气候、地形、植被等因素，缺乏明确的物理过程描述。而Pan-European Soil Erosion Risk Assessment（PESERA）模型作为基于物理过程的区域尺度模型，采用分布式计算框架，能够有效处理区域尺度的空间异质性。同时PESERA能够模拟降雨-产流过程和植被-水文之间的相互关系，提供更为详细的径流、植被与侵蚀模拟结果^［4］。PESERA模型在欧洲土壤侵蚀风险评估中已被广泛应用和改进，如COOPER等^［5］将其运用于英国水蚀风险评估；KARAMESOUTI等^［6］利用该模型模拟希腊部分区域火灾对于土壤侵蚀的响应；LI等^［7］为PESERA模型开发泥炭地冻融和干旱侵蚀模块，形成PESERA模型的泥炭地版本PESERA-PEAT。然而，PESERA在中国应用较少，仅在黄土高原的皇甫川流域有初步应用，因此PESERA模型在中国不同地区的适应性仍有待进一步探讨^［8］。

土壤侵蚀的主要驱动因素包括地形地貌、气候变化、植被覆盖和人类活动^［9］。其中，坡度是关键地形因素之一，较大的坡度加快径流流速，从而增加侵蚀量；降雨直接影响侵蚀过程；温度通过影响地表水的蒸散发间接作用于侵蚀；植被在一定程度上截留地表径流同时加固土壤，从而影响径流量和侵蚀量^［8］。在驱动因素定量分析方面，已有学者^［10-11］利用敏感性分析、相关分析、回归分析、随机森林等方法开展研究。刘悦俊等^［10］利用随机森林回归模型研究了青藏高原土壤侵蚀对气候变化和人类活动的响应，指出年降水量是青藏高原土壤侵蚀的主要驱动因素；AN等^［11］通过逐步回归分析计算土壤侵蚀因子的贡献，评估土壤侵蚀因子对全球四大代表性流域土壤保持服务空间分布的影响。这些方法能够有效评估单一因素的影响，但在定量分析多个变量间的交互作用方面仍有不足。地理探测器通过因子探测和交互探测弥补上述方法的不足，其能够量化不同因素的单因子解释力及其交互作用。因此，地理探测器为土壤侵蚀驱动因素的深入解析提供强有力的工具，被广泛应用于土壤侵蚀相关的研究^［12］。其中SONG等^［13］提出的最优参数的地理探测器（optimal parameters-based geographical detector，OPGD）能够克服传统地理探测器离散化过程的主观性和随机性，从而获得最佳的分类组合。其改进使得多因素交互作用分析更加精确，为土壤侵蚀驱动因素的全面解析提供新的可能。

三江源是青藏高原的核心区域，蕴藏丰富的资源，是我国重要的生态保护区。该区域地貌多样、气候复杂，长期遭受水力、风力、冻融等侵蚀作用的影响，水土流失不断加剧。作为黄河、长江和澜沧江的源头，三江源地区的水土流失不仅威胁到本地生态安全，也对下游省份的生态环境、粮食安全和社会经济产生重要影响^［14］。尽管已有研究^［14-15］采用经验模型分析三江源的土壤侵蚀空间格局变化，但经验模型难以描述土壤侵蚀相关过程，其结果精度和可靠性存在局限性。因此，为保障各水系的生态安全，促进区域高质量发展，亟需引入区域尺度过程模型来揭示侵蚀的时空分布特征及其驱动机制，以弥补经验模型的不足，以便为科学施策提供依据，推动水资源的可持续管理和利用。

鉴于三江源地区土壤侵蚀的复杂性与生态保护的迫切需求，本研究首先校准与验证PESERA模型，评估其在该区域的适用性。基于模型模拟结果，系统地揭示1983—2022年三江源地区土壤侵蚀的时空变化规律。采用最优地理探测器方法深入探究土壤侵蚀的驱动机制，定量分析主要影响因子的解释力与交互作用，揭示区域侵蚀过程的关键控制因素及其动态关系，为三江源地区未来的水土保持规划提供科学依据，助力区域生态环境的可持续发展。

1 材料与方法

1.1 研究区概况

三江源（32°22′36″~36°47′53″N，89°50′57″~99°14′57″E）位于中国青海省南部，是长江、黄河和澜沧江的发源地，总面积约3.03×10⁵ km²（图1）。黄河源区、长江源区和澜沧江源区分别约占总面积的46%、44%和10%。地貌以山地为主，呈西高东低的趋势，海拔1 991~5 786 m。三江源区域河流、湖泊、冰川、沼泽广泛分布，水资源储量丰富、地表径流充足。气候属于典型的高原大陆性气候，年降水量为262~772 mm，年平均气温5.6~3.8 °C^［15］。土地利用类型以草地为主，平均植被覆盖度约为70%，主要植被类型为高寒草甸和高寒草原^［15］。随着气候变化和人类活动的增加，该区域生态环境日益恶化，土地荒漠化、草地退化和水土流失等生态问题严重。自2005年以来，中国政府启动三江源生态保护和建设工程，陆续实施退牧还草、禁牧封育、设立自然保护区等措施，对该区域进行全面的生态保护^［16］。

1.2 研究方法

1.2.1 PESERA模型原理

PESERA模型是具有可靠的物理基础的区域尺度的过程模型，分为水文、植物生长和侵蚀3个模块。其将地表概化为一系列坡面，并假设每个坡面均直接连接至河道（不考虑河道内的土壤侵蚀、输沙和泥沙淤积等过程），便于在区域尺度的应用^［8］。

水文模块主要根据水量平衡原则，将降水量分解为地表径流、土壤水下渗和存储、植被截留和水分蒸散发等过程^［4］。地表径流（f_{\mathrm{0}}）被认为是总降水量中超出径流阈值的部分，计算方法为：

f_{\mathrm{0}}=p\left(\theta -\mu \right)（1）

式中：θ为总降水量，mm；\mu为地表产流阈值，mm；p为降雨量中除去植被截留及土壤吸收下渗后转换为地表径流的比例。

植物生长模块基于生物量碳平衡，每月更新植被覆盖度、植被生物量和土壤有机质。植物生长的过程主要受土壤水分控制，同时也对土壤含水量以及地表产流产生影响。因此，植物生长模块与水文模块紧密联系，通过与水文模块的交互响应，进一步影响侵蚀模块的结果。

在该模型中，地表径流是侵蚀的主要驱动力，土壤侵蚀速率主要是坡面上受降水和径流冲刷携带到河道的泥沙量，其并未对河道侵蚀、沟道侵蚀和河道输沙过程进行估算^［4］，土壤侵蚀速率等同于径流运输泥沙能力的大小。最终的月土壤侵蚀速率结合植被覆盖度的影响，通过将裸地部分和植被覆盖部分的土壤侵蚀速率进行加权求和得到最终的侵蚀速率，计算方法为：

\mathrm{S}{\mathrm{Y}}_{\mathrm{d}}=k{f}_{\mathrm{0}}^m{\mathrm{\Delta }}^n（2）

E=\mathrm{S}{\mathrm{Y}}_{\mathrm{d}}(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{v})+\frac{\mathrm{S}{\mathrm{Y}}_{\mathrm{d}}}{\mathrm{10}}\mathrm{c}\mathrm{o}\mathrm{v}（3）

式中：E为土壤侵蚀量，t∕（hm²·a）；k为土壤可侵蚀性因子；∆为地形起伏度；m、n为经验参数，本研究中，m取值2，n取值1。

1.2.2 模型输入

PESERA模型输入参数主要包括气象、土地利用/覆盖、地形和土壤4个类型，共需128层输入数据。依据数据情况，本文研究时段确定为1983—2022年。将研究时段以10 a为单位划分为4个时段，分别为1983—1992年、1993—2002年、2003—2012年、2013—2022年开展模拟，以便于探究气候变化的周期性和人类活动的长期影响，以及后续驱动力分析。模型输入数据来源见表1。

PESERA模型输入数据通过重采样和降尺度的方法，均统一为1 km×1 km空间分辨率。所需6个气象输入数据中月总降雨量、月平均温度、月潜在蒸散量获取自数据产品，月平均温差由月最高气温和月最低气温的差值计算所得，月降雨天数降雨量根据月总降雨量和降雨天数计算所得、月降雨天数变异系数（\mathrm{c}\mathrm{v}\mathrm{r}{\mathrm{f}}_{\mathrm{2}}）根据公式（4）计算所得，降雨天数降尺度处理根据公式（5）计算。

\mathrm{c}\mathrm{v}\mathrm{r}{\mathrm{f}}_{\mathrm{2}}=\mathrm{1}+\mathrm{18}\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\mathrm{D}\mathrm{i}\mathrm{s}\mathrm{t}(\mathrm{w}\mathrm{e}\mathrm{t},\mathrm{1.7,11},\mathrm{f}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e})（4）

\mathrm{w}\mathrm{e}\mathrm{t}={(\frac{\mathrm{P}\mathrm{R}\mathrm{E}\times \mathrm{W}\mathrm{E}\mathrm{T}\_\mathrm{N}\mathrm{O}\mathrm{R}{\mathrm{M}}^{\frac{\mathrm{1}}{\mathrm{0.45}}}}{\mathrm{P}\mathrm{R}\mathrm{E}\_\mathrm{N}\mathrm{O}\mathrm{R}\mathrm{M}})}^{\mathrm{0.45}}（5）

式中：cvrf₂为降雨天数变异系数；wet为降雨天数，d；PRE_NORM为1 km分辨率降雨量各研究时段逐月多年平均值；WET_NORM为根据CRU TS降雨天数数据集（分辨率0.5°）采用双线性内插方法降尺度为1 km×1 km的各研究时段年逐月多年平均降雨天数；PRE为逐月1 km×1 km分辨率降雨数据。

PESERA模型的地形输入为地形起伏度，根据DEM（digital elevation model）计算为高程的标准差。土壤数据和土地利用数据根据PESERA模型手册进行计算和重分类^［17］。植被覆盖度数据根据像元二分模型进行计算，其中2000年后利用MODIS NDVI数据计算。2000年之前的数据，参考MA等^［18］提出的NDVI降尺度算法，采用AVHRR NDVI数据进行降尺度处理，该算法捕捉AVHRR NDVI和MODIS NDVI之间的统计关系，考虑二者之间的变异系数，并融合这2个数据的时空信息，生成1 km×1 km分辨率的NDVI数据，再基于Google Earth Engine云平台计算植被覆盖度。

1.2.3 模型输出

PESERA模型的输出数据同样包括水文、植被和侵蚀3类。本文主要对模型生成的侵蚀结果进行分析，根据水利部颁布的《土壤侵蚀分类分级标准》（SL 190—2007）^［19］进行分类，分别为微度［<5 t∕（hm²·a）］、轻度［5~25 t∕（hm²·a）］、中度［25~50 t∕（hm²·a）］、强烈［50~80 t∕（hm²·a）］、极强烈［80~150 t∕（hm²·a）］和剧烈［>150 t∕（hm²·a）］。

1.2.4 模型校准和验证

模型校准参数最主要的为径流阈值p（范围在0~1），采用试错法调整p值对模型进行率定。PESERA仅预测坡面的土壤侵蚀率，并未考虑冲沟和河道中的泥沙产生和输送，因此选取坡面的土壤侵蚀率来验证更为准确。本文搜集三江源范围内的多个¹³⁷Cs同位素示踪所测得的侵蚀数据进行验证^［20-23］。根据决定系数（R²）和纳什系数（nash-sutcliffe efficiency coefficient，NSE）评价模型预测值的准确度，采用均方根误差（root mean squared error，RMSE）和平均绝对误差（MAE，mean absolute error）评价模拟值与实测值的数值差异，计算方法为：

R^{\mathrm{2}}=\frac{{\left[{\sum }_i\left(O_i-\overline{O}\right)\left(S_i-\overline{S}\right)\right]}^{\mathrm{2}}}{{\sum }_i{\left(O_i-\overline{O}\right)}^{\mathrm{2}}{\sum }_i{\left(S_i-\overline{S}\right)}^{\mathrm{2}}}（6）

\mathrm{N}\mathrm{S}\mathrm{E}=\mathrm{1}-\frac{{\sum }_i{\left(O_i-S_i\right)}^{\mathrm{2}}}{{\sum }_i{\left(O_i-\overline{O}\right)}^{\mathrm{2}}}（7）

\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt[]{\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^n}{\left(O_i-S_i\right)}^{\mathrm{2}}}（8）

\mathrm{M}\mathrm{A}\mathrm{E}=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^n}\left|O_i-S_i\left|\right.\right.（9）

式中：O_i、S_i 分别为实测值和模拟值；\overline{O}、\overline{S}分别为实测值和模拟值的平均值；N为样本数量。R²与NSE越接近1，表明模型模拟值与实测值的拟合度越高，RMSE和MAE越小，模拟值越接近于实测值。

1.2.5 基于最优地理探测器的土壤侵蚀因子分析

地理探测器作为一种揭示空间异质性的分析模型，可以定量评估自变量对因变量的解释力，探测自变量两两之间的交互作用及其强弱关系^［24］。最优参数地理探测器（OPGD）在传统探测器的基础上克服离散方法具有主观性的缺陷^［13］，进行更加准确的空间分析，其具体步骤为：

根据PESERA模型的输入数据选取植被覆盖度（FVC）、年总降雨量（PRE）、年平均温度（TEM）、高程（DEM）、地形起伏度（STD）作为自变量（X₁~X₅），在研究区构建不同大小（1 km×1 km、5 km×5 km，10 km×10 km）的规则渔网数据，并统计渔网内各自变量的均值，作为采样数据，结果发现5 km×5 km分辨率的渔网数据能够更好地兼顾数据的空间差异及模型运行效率。借助R语言中的GD包实现对以上空间连续数据的最优离散化处理，其分级运用自然间断点法、等间隔法、标准差法、几何间隔分类法及分位数分类法^［13］，将类别设置为5~9类。根据各因子所得最高解释度（q值）对应的离散化分级方法作为最优分级方法，分别进行因子探测、交互探测、风险探测。

因子探测器探测自变量因子（X₁~X₅）对因变量侵蚀速率（Y）的解释力度。交互探测器用于评估自变量之间是否具有交互作用及二者共同作用时对因变量Y的解释力是否具有增强或减弱的作用。二者均用q值（0~1）表示，q值越大，表明解释力度越强。风险探测则用于判断2个子区域间的属性平均值是否有显著的差别。

q=\mathrm{1}-\frac{{\sum }_{N=\mathrm{1}}^L{N}_h{\delta }_h^{\mathrm{2}}}{N{\delta }^{\mathrm{2}}}（10）

式中：q为自变量因子的解释；h为自变量X的分级数；N_h 为层h的单元数；N为全区单元数；{\delta }^{\mathrm{2}}、{\delta }_h^{\mathrm{2}}为全区的侵蚀速率和层h的离散方差。

2 结果与分析

2.1 模型精度验证

根据才北山小流域的2018—2020年实测侵蚀和径流数据对模型参数进行率定（表2），其径流阈值p率定结果为0.75，径流和侵蚀的相对误差分别为0.80%和0.69%，校准结果良好。验证结果在兴海县等区域的相对误差较高在50%以上，其余相对误差在30%左右（表3）。经计算侵蚀验证结果的R²、NSE、RMSE和MAE分别为0.69、0.53、0.024和2.73（图2）。基于以上评估参数，说明基于PESERA模型模拟的结果基本可靠，可在三江源地区进行初步模拟。

2.2 模拟结果及驱动因素的时空分布变化

2.2.1 驱动因子的时空变化

由图3可知，三江源1983—2022年气候变化显著，年平均降水量、年平均气温均呈显著上升趋势。该地区年平均降水量为472 mm，并在1989年达到峰值532 mm；年平均气温为-4.5 °C，在2006年时年平均气温最高为-3.44 ℃，1997年最低为-5.47 ℃。由图4可知，在空间分布上，降雨量呈南高北低的空间分布趋势，黄河源南部为降水量最高的区域。温度呈西高东低的空间分布状况，澜沧江源南部及黄河源东部温度较高。三江源植被覆盖度平均值约为34%，1983—2022年呈显著上升趋势，2000年之后的平均植被覆盖度高于2000年之前的植被覆盖度，且在2010年时植被覆盖度达到峰值，为36.7%。在空间分布上，长江源区域的植被覆盖度较低，黄河源东部地区的植被覆盖度较高。气候和植被覆盖的时空分布变化进一步反映三江源地区气候向“暖湿化”方向发展的特征，以及“退牧还草”等生态工程政策的干预对生态恢复具有一定积极影响。三江源地区地势东高西低，由长江源到黄河源逐渐降低，地形起伏度高值集中在澜沧江源区域及黄河源中部，长江源地区地形较为平坦，澜沧江源地区由于其气候、植被和地形的共同影响，该区域为侵蚀监测的重点区域。

2.2.2 三江源土壤侵蚀速率的时空分布变化

在空间分布上，三江源土壤侵蚀呈南高北低的空间分布，其中微度侵蚀主要分布在西北部，剧烈侵蚀则集中在南部（图5）。长江源地区，西部以微度侵蚀为主，东南部分布少量极强烈及剧烈侵蚀；黄河源地区，南部主要为轻度侵蚀，剧烈侵蚀呈分散点状分布；澜沧江源地区侵蚀程度整体较高，主要侵蚀类型为强烈及以上。对于不同侵蚀程度的面积占比，其中轻度侵蚀占比最高，并在2003—2012年达到峰值，约占总面积的40.4%，随后出现下降。微度侵蚀面积占比产生先下降后上升的变化，在2013—2022年达到峰值，为39.3%（图6a）。强烈侵蚀及以上等级面积占比较低且波动较小，均呈先上升后下降的变化过程，其中剧烈侵蚀占比最低，约占总面积的5%以下。综上所述，对于南部地形陡峭且降水高的区域，侵蚀更为严重；而西北部地形较为平坦且降水较少的区域，侵蚀程度较低。

三江源地区1983—1992年、1993—2002年、2003—2012年和2013—2022年的平均侵蚀速率分别为25.75、24.05、32.9、30.69 t∕（hm²·a），此变化过程呈先上升后下降的特征。黄河源、长江源、澜沧江源1983—2022年的平均侵蚀速率同样呈先上升后下降的变化模式，均在2003—2012年时段内达到峰值，其中澜沧江源的侵蚀速率远高于其他区域（图6b）。其可能是因为2003—2012年更多的人类活动（如铁路建设等）导致侵蚀加剧，而随着生态工程（如退耕还草、天然林保护等）的实施可能促进侵蚀的恢复。不同侵蚀等级在各区域的面积占比存在明显差异，其中约有70%的微度侵蚀分布在长江源地区，1%分布在澜沧江源地区。轻度和微度侵蚀在黄河源区的面积占比最大，约有50%的轻度侵蚀和45%微度侵蚀分布在黄河源地区。对于极强烈及以上的侵蚀等级，约有30%以上的极强烈侵蚀和50%以上的剧烈侵蚀，分布在澜沧江源地区（图7）。针对澜沧江源地区剧烈侵蚀严重的现象，未来需优先实施坡面固结工程，减少该区域的水土流失，以促进生态环境的恢复。

2.3 基于最优参数地理探测器（OPGD）的土壤侵蚀影响因素分析

2.3.1 因子探测分析

由图8可知，各因子对三江源整体区域的土壤侵蚀速率的影响程度的排序为STD>PRE>TEM>DEM>FVC。对于长江源区域，地形起伏度（STD）和降雨量（PRE）为该区域侵蚀变化的主导因素，其解释力平均值分别为0.43和0.35。相比长江源地区，黄河源区域的降雨量（PRE）和地形起伏度（STD）的解释力较低，温度（TEM）和高程（DEM）的解释力则更高。在1983—1992年黄河源地区的解释力STD>PRE>TEM；在1993—2022年3个时段内，解释力STD>TEM>PRE。澜沧江源区域，植被覆盖度（FVC）和地形起伏度（STD）的解释力分别最高和最低，其中1983—2002年2个时段内植被覆盖度对侵蚀的解释力最高，高程次之；在2003—2022年2个时段内植被覆盖度对侵蚀的解释力最高；降雨量次之。黄河源地区温度（TEM）的解释力增强说明气候变暖可能加剧该区域的侵蚀（如冻融作用增强）。在澜沧江源降雨量（PRE）解释力的增加可能因极端降雨事件增多或人类活动（如工程建设）改变地表径流条件，从而导致侵蚀加剧。因此，土壤侵蚀驱动因子在三江源区域具有显著的区域异质性和时间动态性，其治理策略需因地制宜。

2.3.2 交互探测分析

由图9可知，在三江源整体区域，双因子探测的解释力普遍高于单因子。在4个时段内，不同区域的交互探测结果表明，三江源整体区域FVC ∩ PRE的解释力最高，为0.730。同时FVC在与PRE和STD的交互中显著增强对侵蚀的解释力，表明FVC具有重要辅助作用。长江源地区的交互探测分析结果与三江源整体类似，其中FVC ∩ PRE在该区域的解释力最高为0.672，其次为PRE ∩ TEM，其解释力为0.622。在黄河源地区，交互探测的结果表明，FVC ∩ PRE、TEM ∩ STD 、FVC ∩ STD的解释力均较高，分别为0.573、0.553、0.523。在澜沧江源区域，各因子之间的交互探测结果相较于其他区域均较高，其中FVC ∩ STD的对该区域解释力最高，为0.777。所有区域的交互探测结果均表明，地形起伏度（STD）在与其他因子的交互中具有显著的增强作用，说明地形起伏度（STD）不仅直接影响侵蚀空间分异，还通过与其他因子的共同作用显著提高对侵蚀空间分异的解释力。

2.3.3 风险探测分析

由图10可知，在长江源区域植被覆盖度、降雨量、温度、高程、地形起伏度的最大风险范围分别为21%~26.6%、>532 mm、-11~-8 °C、5 300~5 700 m、>104 m。在黄河源区域植被覆盖度、降雨量、温度、高程、地形起伏度的最大风险范围分别为24.7%~33.3%、>675 mm、-11.2~-7.5 °C、>4 650 m、113 m。而在澜沧江流域植被覆盖度<25.5%、年总降雨量>644 mm、温度-9.8~-6.2 °C、高程>4 900 m、地形起伏>142 m的区域侵蚀风险最大。综合各区域分析，在降雨量较高、植被覆盖较低、地形起伏度较高的区域侵蚀容易发生；反之侵蚀风险较低。

3 讨 论

3.1 PESERA模型精度验证

本研究模拟结果与基于¹³⁷Cs的实测侵蚀数据存在显著相关性（R²=0.69，p<0.05）、模型精度较好（NSE=0.53）。基于¹³⁷Cs的实测侵蚀数据解决以往研究在使用土壤侵蚀空间分布结果验证时，由于缺乏与实测数据的对比而产生的精度受限问题^［25］。此外，水文站点输沙数据也是侵蚀模型验证的重要数据源。然而，由于水文站点控制面积较大，泥沙的传输和沉积过程对结果影响较大，难以准确反映实际侵蚀量^［16］，基于¹³⁷Cs侵蚀数据则能够直接反映土壤在局部尺度上的侵蚀特征，有效提高验证的可靠性。

在ZHAO等^［26］的验证区域中（三江源地区），本研究的模拟侵蚀速率比RUSLE的结果更接近实测值，在达日县、玛沁县军功乡及秦珍乡RUSLE模拟结果与实测值的相对误差分别为44.85%、38.75%、10.50%，高于PESERA模拟的相对误差（31.33%、28.13%、8.91%）（表3）。总体来看，PESERA模型在空间异质性和过程模拟方面优于传统经验模型，能够为三江源地区的土壤侵蚀精准评估及科学管理提供更为可靠的理论支撑。

然而，PESERA模拟结果依然存在误差。尤其在兴海县、野牛沟乡等区域仍存在相对误差较大的问题，其相对误差超过50%，存在低估现象。可能是由于该地区沟道变化较为复杂，多种侵蚀现象共存，而PESERA模型模拟过程并未考虑沟道侵蚀、重力侵蚀等，仅考虑细沟和细沟间侵蚀，从而产生低估^［8］。此外，多种侵蚀营力的作用、过度放牧导致牲畜的过度啃食对土壤质地造成影响及输入数据的分辨率等客观因素均可能导致模拟结果产生一定误差^［16］。

对PESERA模拟结果的时空变化格局分析表明，在时间变化趋势上，各区域侵蚀速率均在1983—2012年有所增加，2013—2022年有所下降，三江源整体区域相较于峰值期2003—2012年，其平均侵蚀速率下降约6%，其中长江源区的下降最为显著，相较于峰值期约下降13%。不同侵蚀程度的面积占比在2013—2022年有明显变化，其中强烈及以上等级侵蚀面积占比减少，微度侵蚀面积占比有大幅增加。此变化主要受益于2005年以来，对三江源地区的生态保护，如建立自然保护区、退牧还草等生态工程的实施有效遏制侵蚀加剧的趋势。但目前侵蚀速率仍然高于1982—2002年，其原因可能为生态保护措施（如退牧还草、植被恢复）虽然在一定程度上减缓侵蚀速率，但由于滞后效应和局部效果差异，生态环境尚未完全恢复^［16］。三江源地区的气温升高、降水量增加、极端气候事件增加是导致侵蚀速率仍然较高的主要原因^［27］。气候变化导致冻融侵蚀的加剧，从而使土壤结构破坏、土壤松散度增加，进一步放大水力侵蚀的作用。然而，需要注意的是，PESERA模型并未考虑冻融循环对土壤结构和可侵蚀性等的影响，未来可结合LI等^［7］研究，进一步优化对冻融侵蚀的表述。

空间分布上，该研究结果与RUSLE模拟结果基本一致^［27］，均表现为侵蚀速率在地形起伏较大的丘陵沟壑区较高，平坦区域较低。但相比RUSLE模拟结果，PESERA结果中强烈及以上等级侵蚀的面积占比更高，微度侵蚀的面积占比则较低，表明PESERA在捕捉高强度侵蚀方面更具敏感性，更适合反映三江源地区剧烈侵蚀的空间分布特征。

在侵蚀量级上，PESERA模拟结果略低于贺倩等^［25］模拟结果［37.96 t∕（hm²·a）］但高于曹巍等^［16］利用RUSLE模拟结果［11.4 t∕（hm²·a）］。造成此差异的原因除数据来源、计算方法及分辨率差异外，模型原理的不同也对结果产生影响。三江源地区的土壤侵蚀机制较为复杂，而RUSLE作为经验模型，未能充分考虑这些复杂因素。但PESERA作为过程模型，通过分布式模拟计算地表径流量，并将植被模块与水文模块耦合，能够同时考虑植被对降雨的截留作用、蒸散、入渗等关键过程，从而更合理地再现土壤侵蚀的发生机制^［8］。因此，相比RUSLE依据经验公式直接计算侵蚀速率，其理论依据更为严谨。

3.2 三江源土壤侵蚀的主要驱动力

对于单因子探测分析，不同区域的主导因素具有一定差异，长江源和黄河源地区地形较为平坦，局部地区较为明显的地形变化导致侵蚀速率的变化显著，因此，该区域地形起伏度（STD）对侵蚀速率的解释力最高。而在地形本身较为复杂的澜沧江源地区，较高的植被覆盖度能够对侵蚀产生抑制作用，因此，澜沧江源地区植被覆盖度（FVC）的单因子解释力最高，地形起伏度（STD）的解释力最低。其与WANG等^［28］所得出的坡度对土壤侵蚀的解释力在相对平坦的地区显著高于地势较高的地区的结论类似。

双因子交互的解释力均高于单因子探测结果，说明侵蚀现象是多因素共同作用的结果，影响机制复杂多变，与廖洪圣等^［29］所得结论一致。降雨量（PRE）和植被覆盖度（FVC）的交互作用对三江源侵蚀速率空间分异解释力均最高，是影响侵蚀变化的主要因素，与LU等^［30］和吴成勇等^［31］分别利用地理探测器和地理加权回归在青藏高原和黄河上游地区探究影响的土壤保持服务主导因子的研究结果类似。三江源属于高原大陆性气候，植被以高寒草甸和高寒草原为主^［9］，降水量的变化直接影响植被的生长状况。适度的降水可以促进植被生长，提升植被覆盖度，缓解侵蚀速率的升高。但过多的降水量增加降水侵蚀力，从而导致侵蚀加剧。地形起伏度与其他因素的交互均增强对侵蚀速率的解释力，由于较高的地形起伏度使得土壤容易受到水流冲刷等外力的作用^［28］，若降水量增加则导致径流增大，流速加快，土壤水下渗和存储减少，植被截留的水分减少，泥沙迁移增多，进而导致侵蚀增加^［26］。如长江源地区的东南部，地形起伏度和降水量均较高，二者的协同作用共同导致该地区侵蚀速率较高。因而，处于复杂地形、稀疏植被及较高降雨量的区域受到的侵蚀风险程度最高，故未来应因地制宜实施植被恢复、修建梯田、减少放牧等水土保持措施，以缓解土壤侵蚀对生态安全造成的负面影响。

4 结论

1） PESERA模型在三江源地区验证结果良好，其R²和NSE分别为0.69和0.53，但由于多种侵蚀营力的作用和过度放牧等因素，在兴海县等局部地区存在低估现象。

2） 1983—2022年的平均侵蚀速率为28.36 t∕（hm²·a），各区域平均侵蚀速率平均值呈先上升后下降的变化特征，均在2003—2012年达到峰值。侵蚀空间分布呈南高北低的分布格局。微度侵蚀、轻度侵蚀和剧烈侵蚀分别主要分布在长江源地区、黄河源地区和澜沧江源地区。轻度侵蚀和微度侵蚀面积占比均较高，分别占总面积35%和30%以上，剧烈侵蚀面积占比最低，约为总面积的3%。

3）在降雨量较高、植被覆盖度较低且地形起伏度较高的区域土壤侵蚀风险较高。在地形较为平坦（长江源和黄河源地区）和复杂（澜沧江源地区）的地区，地形起伏度和植被覆盖度的解释力分别最高。降雨量和植被覆盖度的交互作用在三江源地区具有最高解释力，地形起伏度与其他因子的交互作用，具有增强解释力作用，植被覆盖度在与其他因子的交互过程中具有重要辅助作用。

参考文献

[1] BORRELLI P， ROBINSON D A， PANAGOS P， et al. Land use and climate change impacts on global soil erosion by water （2015—2070）［J］.Proceedings of the National Academy of Sciences of the United States of America，2020， 117（36）：21994-22001.

[2] LI P F， LI D， HU J F， et al. Improving the application of UAV-LiDAR for erosion monitoring through accounting for uncertainty in DEM of difference［J］.Catena，2024，234：e107534.

[3] MAQSOOM A， ASLAM B， HASSAN U， et al.Geospatial assessment of soil erosion intensity and sediment yield using the revised universal soil loss equation （RUSLE） model［J］.ISPRS International Journal of Geo-Information，2020，9（6）：e356.

[4] LI P F， MU X M， HOLDEN J， et al. Comparison of soil erosion models used to study the Chinese Loess Plateau［J］.Earth-Science Reviews，2017，170：17-30.

[5] COOPER D， FOSTER C， Goodey R， et al. Use of ＇UKCIP08 Scenarios＇ to determine the potential impact of climate change on the pressures/threats to soils in England and Wales［M］.Final Report to Defra， Project SP0571. NERC，2010.

[6] KARAMESOUTI M， PETROPOULOS GP， PAPANIKOLAOU ID， et al. Erosion rate predictions from PESERA and RUSLE at a Mediterranean site before and after a wildfire： Comparison and implications［J］.Geoderma， 2016，261：44-58.

[7] LI P F， HOLDEN J， IRVINE B ．Prediction of blanket peat erosion across Great Britain under environmental change［J］.Climatic Change，2016，134（1/2）：177-191．

[8] LI P F， ZANG Y Z， MA D D， et al. Soil erosion rates assessed by RUSLE and PESERA for a Chinese Loess Plateau catchment under land-cover changes［J］.Earth Surface Processes and Landforms，2020，45（3）：707-722.

[9] GUO Y R， PENG C H， ZHU Q A， et al. Modelling the impacts of climate and land use changes on soil water erosion： Model applications， limitations and future challenges［J］.Journal of Environmental Management，2019，250：e109403.

[10] 刘悦俊，程江浩，张燕杰.青藏高原土壤侵蚀对气候变化和生态政策的非线性响应［J］.水土保持研究，2024，31（4）：126-134.

[11] LIU Y J， CHENG J H， ZHANG Y J. Nonlinear response of soil erosion in the Tibetan Plateau to cliamte change and ecological policies［J］.Research of Soil and Water Conservation，2024，31（4）：126-134.

[12] AN Y M， ZHAO W W， LI C J， et al. Temporal changes on soil conservation services in large basins across the world［J］.Catena，2022，209：e105793.

[13] 彭大力，张斌，吴林蓬，等.龙川江流域近20年土壤侵蚀时空变化及驱动因素分析［J］.水土保持学报，2024，38（4）：29-37.

[14] PENG D L， ZHANG B， WU L P，et al.Analysis of spatial and temporal changes and driving factors of soil erosion in Longchuan River basin in recent 20 years［J］.Journal of Soil and Water Conservation，2024，38（4）：29-37.

[15] SONG Y Z， WANG J F， GE Y， et al. An optimal parameters-based geographical detector model enhances geographic characteristics of explanatory variables for spatial heterogeneity analysis： Cases with different types of spatial data［J］.GIScience and Remote Sensing，2020，57（5）：593-610.

[16] 黄婷婷，赵辉，赵院，等.三江源国家公园土壤侵蚀及其分布特征［J］.水土保持通报，2023，43（5）：95-103.

[17] HUANG T T， ZHAO H， ZHAO Y， et al. Soil erosion and its spatial distribution characteristics in three-river-source National Park［J］.Bulletin of Soil and Water Conservation，2023，43（5）：95-103.

[18] 陶洁怡，章锦河，郭丽佳，等.三江源国家公园水源涵养服务的时空演变特征及驱动因素研究［J］.生态学报，2024，45（3）：1226-1238.

[19] TAO J Y， ZHANG J H， GUO L J， et al. Temporal and spatial evolution characteristics and driving factors of water conservation service in Sanjiangyuan National Park［J］.China Industrial Economics，2024，45（3）：1226-1238.

[20] 曹巍，刘璐璐，吴丹.三江源区土壤侵蚀变化及驱动因素分析［J］.草业学报，2018，27（6）：10-22.

[21] CAO W， LIU L L， WU D. Soil erosion changes and driving factors in the Three-River Headwaters region［J］.Acta Prataculturae Sinica，2018，27（6）：10-22.

[22] KIRKBY M J， IRVINE A B J， JONES A R J A， et al. The PESERA coarse scale erosion model for Europe. I.-Model rationale and implementation［J］.European Journal of Soil Science，2008，59（6）：1293-1306.

[23] MA Z M， DONG C Y， LIN K R， et al. A global 250-m downscaled NDVI product from 1982 to 2018［J］.Remote Sensing，2022，14（15）：e3639.

[24] 中华人民共和国水利部. 土壤侵蚀分类分级标准：SL 190-2007 ［S］.北京：中国水利水电出版社，2008.

[25] Ministry of Water Resources of the People's Republic of China. Standards for classification and gradation of soil erosion：SL 190-2007 ［S］.Beijing： China Water and Power Press，2008.

[26] 李元寿，王根绪，王军德，等. ¹³⁷Cs示踪法研究青藏高原草甸土的土壤侵蚀［J］.山地学报，2007，25（1）：114-121.

[27] LI Y S， WANG G X， WANG J D， et al. ¹³⁷Cs trace technique to study soil erosion at alpine meadow of Tibetan Plateau［J］.Journal of Mountain Science，2007，25（1）：114-121.

[28] 李俊杰.应用于土壤侵蚀的环境放射性核素示踪技术方法研究［D］.北京：中国农业科学院，2008.

[29] LI J J. Study for tracing techniques and methods using enviroment radionuclides apllied to soil erosion［D］.Beijing： Chinese Academy of Agricultural Sciences，2008.

[30] 邵全琴，肖桐，刘纪远，等.三江源区典型高寒草甸土壤侵蚀的¹³⁷Cs定量分析［J］.科学通报，2011，56（13）：1019-1025.

[31] SHAO Q Q， XIAO T， LIU J Y， et al. Quantitative analysis of soil erosion of typical alpine meadow in the source region of Three Rivers by ¹³⁷Cs［J］.Chinese Science Bulletin，2011，56（13）：1019-1025.

[32] LI Y， LI J J， ARE K S， et al. Livestock grazing significantly accelerates soil erosion more than climate change in Qinghai-Tibet Plateau： Evidenced from ¹³⁷Cs and ²¹⁰Pbex measurements［J］.Agriculture， Ecosystems and Environment，2019，285：e106643.

[33] 王劲峰，徐成东.地理探测器：原理与展望［J］.地理学报，2017，72（1）：116-134.

[34] WANG J F， XU C D. Geo detector： Principle and prospective［J］.Acta Geographica Sinica，2017，72（1）：116-134.

[35] 贺倩，戴晓爱.基于LMDI模型的三江源区植被对土壤侵蚀变化影响的定量分析［J］.长江科学院院报，2020，37（7）：61-67.

[36] HE Q， DAI X A. Quantitative analysis of the effects of vegetation on soil erosion in the Three-River Headwaters Region using logarithmic mean divisia index model［J］.Journal of Yangtze River Scientific Research Institute，2020，37（7）：61-67.

[37] ZHAO W T， CHENG Y Z， JIAN J S， et al. Water erosion changes on the Qinghai-Tibet Plateau and its response to climate variability and human activities during 1982—2015［J］.Catena，2023，229：e107207.

[38] WANG Z， WANG J B. Changes of soil erosion and possible impacts from ecosystem recovery in the Three-River Headwaters Region， Qinghai， China from 2000 to 2015［J］.Journal of Resources and Ecology，2019，10（5）：e461.

[39] WANG H， GAO J B， HOU W J. Quantitative attribution analysis of soil erosion in different geomorphological types in karst areas： Based on the geodetector method［J］.Journal of Geographical Sciences，2019，29（2）：271-286.

[40] 廖洪圣，卫伟，石宇.黄土丘陵区典型流域土壤侵蚀时空演变特征及其驱动机制：以祖厉河为例［J］.生态环境学报，2024，33（6）：908-918.

[41] LIAO H S， WEI W， SHI Y. Characteristics of spatial and temporal evolution of soil erosion in typical watersheds in loess hilly areas and its driving mechanisms：A case study of Zuli River［J］.Ecology and Environmental Sciences，2024，33（6）：908-918.

[42] LU R R， DAI E F， WU C S. Spatial and temporal evolution characteristics and driving factors of soil conservation services on the Qinghai-Tibet Plateau［J］.Catena，2023，221：e106766

[43] 吴成永，曹广超，陈克龙，等.黄河上游地区土壤保持服务时空变化及归因［J］.水土保持学报，2022，36（4）：143-150.

[44] WU C Y， CAO G C， CHEN K L， et al. Spatio-temporal variation in soil conservation service and its influencing factors in the upper reaches of the Yellow River［J］.Journal of Soil and Water Conservation，2022，36（4）：143-150.

基金资助

国家自然科学基金项目(42207407)

国家自然科学基金项目(U2243211)

榆林高新区“科学家+工程师”项目(YGXKG-2023-105)

中国科学院A类战略性先导科技专项(XD20040202)