土壤侵蚀严重威胁人类社会的安全和可持续发展，目前全球土壤侵蚀依旧严峻。2019年世界土壤日把土壤侵蚀认定为土壤的头号威胁。中国作为世界上土壤侵蚀最严重的国家之一，据统计依然有超过50%的国土面积仍在遭受着土壤侵蚀的危害^[1-2]。西藏自治区“一江两河”(雅鲁藏布江、年楚河、拉萨河)地区作为中国西部典型的开发区，同时也是中国典型的生态脆弱区之一，其主要范围包括雅鲁藏布江及其两条主要支流拉萨河、年楚河的中部河谷地区。该区域地形落差较大，植被覆盖率低，加之全球气候变暖以及人类活动的影响，生态系统极为脆弱，水土流失严重^[3-5]。区域内侵蚀类型多样，以水力侵蚀、风力侵蚀和冻融侵蚀为主。其中水力侵蚀作用主要发生在拉萨河流域及雅鲁藏布江中下游山南地区，风力侵蚀主要发生在以辫状及乱流状水系极为发育的雅鲁藏布江宽谷地段，冻融侵蚀主要集中分于雅鲁藏布江北岸高寒山区^[5-6]。严重的水土流失已经造成“一江两河”地区土地生产力退化、水利设施淤积及生态环境负担加重等危害^[3]。随着“一江两河”地区的快速发展，对于生态环境质量的要求也愈来愈高，严重的土壤侵蚀所造成的生态退化已经成为制约该地区经济发展的重要因素之一^{[5, 7]}。因此，“一江两河”地区水土流失的有效治理，对于该地区可持续发展具有重要意义。目前，针对“一江两河”地区的土壤侵蚀研究已逐渐受到关注。同时由于高寒环境下土壤侵蚀监测较为困难导致实测资料相对较少，利用RUSLE, CSLE等土壤侵蚀模型对土壤侵蚀量进行计算，已经成为评估该区域土壤侵蚀状况的主要工具。Wang等^[8]利用RUSLE模型对雅鲁藏布江中游的土壤侵蚀进行了估算，研究发现，虽然2001—2015年土壤侵蚀呈下降趋势，但由于降水量和气温的增加，到2050年土壤侵蚀风险可能会加剧。方广玲等^[9]以RUSLE与空间信息技术(GIS和RS)相结合的方法对拉萨河进行了土壤侵蚀特征分析，该流域中度侵蚀面积比例达59.0%，强烈以上侵蚀面积很小，但侵蚀量所占比例较大，达到14.3%。虽然已有部分学者对“一江两河”地区的土壤侵蚀状况进行了评估，但多集中于对单年或多年平均土壤侵蚀状况的评价，而针对不同时期的土壤侵蚀时空变化规律等研究还较为缺乏。同时由于该地区降水集中，雨强及山地径流较大，其水力侵蚀较为严重。并且由于近年气候变化的作用，该地区的植被以及降雨等发生了较大的变化，对于该地区水力侵蚀产生了较大的影响。目前关于气候变化条件下该地区水力侵蚀的时空变化特征研究较少。在全球变化和经济快速发展的背景下，深入研究和理解“一江两河”地区土壤侵蚀时空分异规律，可有力支撑该区水土流失治理和可持续发展的需求。为此，本文以“一江两河”地区为研究对象，基于RUSLE模型和GIS平台，分析1995—2020年土壤侵蚀模数、强度等级及水土流失面积的时空变化特征，以期为该区土壤侵蚀防治和生态治理提供科学依据。

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

西藏自治区“一江两河”地处西藏自治区的中南部(28°60′—30°29′N, 87°4′—92°37′E)，主要包括雅鲁藏布江及其两条主要支流拉萨河、年楚河的中部河谷地区，其中雅鲁藏布江贯穿地区东西，年楚河位于地区西部、雅鲁藏布江以南，拉萨河位于地区东部、雅鲁藏布江以北。该区东西长约520 km，南北宽约220 km，区域面积约6.67万km²，涵盖拉萨市、日喀则市、山南市下辖18个区县。区域海拔介于3 284~7 168 m之间且落差较大，地形空间分异性强，南北海拔高且多山地，年楚河、拉萨河以及雅鲁藏布江周边的河谷地区海拔较低^[10]。气候属高原温带半干旱季风气候，多年平均气温为6.82 ℃，多年平均降水量约为344.24 mm，降水主要集中在5—9月，占全年总降水量90%以上。植被类型主要以草地以及高原草甸为主^[11]，耕地主要分布于河谷地带。该区总人口占西藏自治区总人口33%，在西藏自治区的经济发展中有重要地位^[12]。当地农业与畜牧业较为发达，但其易受自然灾害等影响，抵抗外部风险的能力较弱^[13]。在较为匮乏的自然资源背景、经济发展及人口增长的压力之下，出现了过度开垦、过度放牧等问题，进一步加剧了当地的土壤侵蚀^[5]。

1.2 数据来源与处理

本研究采用的数据主要有降雨、数字高程模型DEM、土壤、归一化植被指数(NDVI)及土地覆被等(表 1)。NDVI数据为美国国家航空航天局(National Aeronautics and Space Administration, NASA)全球监测与模型研究组(Global Inventor Modeling and Mapping Studies, GIMMS)发布的NOAA AVHRR数据集和NASA陆地过程分布式数据档案中心(The Land Processes Distributed Active Archive Center, LP DAAC/NASA)提供的EOS MODIS NDVI数据集。本研究采用的AVHRR GIMMS NDVI 3g V1半月时段为198201—201512，空间分辨率为8 km; MODIS MOD13Q1的16 d数据时段为200002—202112，空间分辨率250 m。每月和每年的NDVI采用了最大值合成法(Maximum value composites, MVC)，以进一步消除云雾、大气等因素的干扰，提高数据的精度。GIMMS NDVI与MODIS NDVI数据的重叠时段为200002—201512，基于以上时段建立逐月数据的线性关系，并将此线性关系应用到GIMMS NDVI 198201-200001时段，得到新的时段数据，以完成上述2种数据的集成^[14-15]。利用两个数据间的线性拟合关系调整了198201—200003的GIMMS NDVI数据，并与200004—202112 MODIS NDVI数据结合(NDVI_MODIS=1.035 7 NDVI_GIMMS-0.023 1, R²=0.770)，形成新的1982—2021年时间序列NDVI数据。

降雨数据采用1986—2020年逐月降水数据，来自于WorldClim的Global climate and weather data数据(1960—2018年)和ERA5_Land逐月平均降水量及平均气温数据(1986—2020年)。WorldClim逐月降水量及气温数据集由东英吉利大学气候研究室的CRU-TS-4.03降尺度而来，并使用WorldClim 2.1数据进行了偏差校正^[16-17]。ERA5_Land数据集是陆面应用中最先进的全球再分析数据集之一，分辨率为9 km(1950至今)。本研究中主要采用WorldClim的CRU降尺度数据，其分辨率为4.5 km，时间序列为198601—201812，但缺乏201901—202012数据。而ERA5_land降水数据与CRU, CMFD等其他数据间在时间趋势上表现出高度的一致性，但是整体偏高^[18]。因此，可借鉴NDVI数据融合的方法，将WorldClim 198601—201812与ERA5-land 201901—202012调整为198601—202012时间序列降水数据(PRE_CRU=0.466 PRE_ERA5, R²=0.971)。

土壤数据为国际土壤参考和资料中心(International Soil Reference and Information Centre, ISRIC)的SoilGrids土壤剖面数据，空间分辨率为250 m。与HWSD全球土壤数据库、中国土壤普查数据等相比，该数据库中土壤属性数据更全、分辨率更高，因此更适用于K因子的计算^[19]。数字高程数据DEM为SRTM 90 m Digital Elevation Database v4.1，由美国太空总署(NASA)和国防部国家测绘局(NIMA)联合测量，是计算坡度坡长因子的主要数据来源之一。土地覆被数据采用了欧洲航天局1992—2020年The Land Cover CCI土地覆被数据(空间分辨率300 m)。土地利用分类系统参照欧洲航天局提供的分类标准将研究区地类重分类为农地、林地、草地、水域、灌木林地、稀疏植被、湿地、居民地、裸地、水域和冰川/雪地等。

1.3 研究方法 1.3.1 土壤侵蚀模型

采用修正的通用土壤侵蚀模型(RUSLE)模拟计算该区1995，2005，2015和2020年土壤侵蚀模数，以期探讨该区土壤侵蚀潜在风险及时空分异规律。RUSLE模型的具体形式为:

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

(1)

式中：A为土壤流失模数〔t/(hm²·a)〕; R为降雨侵蚀力因子〔MJ·mm/(hm²·h·a)〕; K为土壤可蚀性因子〔t·hm²·h/(hm²·MJ·mm)〕; L为坡长因子，无量纲; S为坡度因子, 无量纲; C为作物与覆盖因子, 无量纲; P为水土保持措施因子, 无量纲。

1.3.2 降雨侵蚀力因子R值

根据逐月降雨量估算降雨侵蚀力R因子，采用章文波提出的简易模型进行计算。该方法针对逐月降雨数据进行了参数优化，以包含西藏自治区拉萨气象站的66个气象站资料进行验证，其模型决定系数达到0.861，计算结果精度较高^[20]。

$ R=0.183\;3 \times\left(\sum\limits_{i=1}^{12} P_i^2 / P\right)^{1.995\;7} $

(2)

式中：R为降雨侵蚀力因子〔MJ·mm/(hm²·h·a)〕; P_i为第i月的降水量(mm); P为年降水量(mm)。基于1986—1995年、1996—2005年、2006—2015年和2011—2020年平均月降水量分别计算得到1995，2005，2015和2020年降雨侵蚀力(图 1)。

1.3.3 土壤可蚀性因子K值

基于土壤属性数据，采用EPIC模型中的土壤可蚀性因子计算方法，该方法在国内外应用均较为广泛。

$ \begin{array}{l} K= \left\{0.2+0.3 \times \exp 〔 -0.025\;6 \times S_{\mathrm{ND}} \times\left(1-\frac{S_{\mathrm{LT}}}{100}\right)〕\right\} \times\left(\frac{S_{\mathrm{LT}}}{C_{\mathrm{LY}}+S_{\mathrm{LT}}}\right)^{0.3} \times 〔 1.0-\frac{0.25 \times C}{C+\exp (3.72-2.95 \times C)}〕\times \\ \;\; \;\;{〔 1.0-\frac{0.7 \times \mathrm{SN}_1}{\mathrm{SN}_1+\exp \left(-5.51+22.9 \mathrm{SN}_1\right)}〕 \times 0.131\;7 \;\;\;\; \; \quad \mathrm{SN}_1=1-\frac{S_{\mathrm{AN}}}{100} } \end{array} $

(3)

式中：K为土壤可蚀性因子〔t·hm²·h/(hm²·MJ·mm)〕; S_AN为黏粒含量(%); S_LT为粉粒含量(%); S_ND为砂粒含量(%); C为土壤有机碳(%); 0.131 7为美制单位值转换公制单位的系数。本研究基于SoilGrids土壤数据, 采用0—30 cm土层的平均砂粒、粉粒、黏粒及有机碳含量进行土壤可蚀性K值的计算(图 2)。

1.3.4 坡度坡长因子LS值

坡度因子S依据不同的坡度范围选用相应的公式计算，其中10°以下采用McCool等^[21]的公式计算，10°以上采用Liu等^[22]的公式计算; 坡长因子L采用Foster等^[23]提出的坡长因子公式计算。

基于DEM数据，采用符素华等^[24]设计的《坡度坡长因子计算工具》来实现(图 3)。

$ L=\left(\lambda_{\text {out }}^{m+1}-\lambda_{\text {in }}^{m+1}\right) /〔 \left(\lambda_{\text {out }}-\lambda_{\text {in }}\right) \times 22.13^m〕 $

(4)

$ S= \begin{cases}10.8 \sin \theta+0.03 & \theta<5^{\circ} \\ 16.8 \sin \theta-0.50 & 5^{\circ} \leqslant \theta<10^{\circ} \\ 21.9 \sin \theta-0.96 & \theta \geqslant 10^{\circ}\end{cases} $

(5)

式中：L为坡长因子，无量纲; S为坡度因子，无量纲; θ为坡度(°); λ_out，λ_in分别为栅格出口及入口的坡长(m); m为坡长指数, 根据刘宝元等^[25]的研究结果进行取值: m=0.2, θ＜0.5°; m=0.3, 0.5°≤θ＜1.5°; m=0.4, 1.5°≤θ＜3°; m=0.5, θ＞3°。

1.3.5 植被覆盖和管理因子C值

由于NDVI与植被覆盖和管理C因子具有良好的非线性关系，采用Knijff提出的C因子模型进行计算^[26-27](图 4)，该方法在大尺度区域计算C值是最有效、最实用的方法^[28-30]。

$ C=\exp 〔 -\alpha \frac{\mathrm{NDVI}}{\beta-\mathrm{NDVI}}〕 $

(6)

式中：C为作物与覆盖因子，无量纲; α，β均为常数, 本研究分别取值2.5, 1^{[8, 31]}。

1.3.6 土壤保持措施因子p值

依据青藏高原土壤侵蚀相关研究^[32-34]，结合土地利用图对研究区各地类进行水土保持措施因子赋值(图 5)。其中农地和林地分别赋值为0.4和0.7，居民地和水域赋值为0，草地、灌木林地、稀疏植被地、裸地和冰川/雪地等地类赋值为1。

在ArcGIS平台中将计算得出的RUSLE模型各因子的空间分辨率统一重采样至250 m，相乘得出各年土壤侵蚀模数。根据《土壤侵蚀分类分级标准(SL190—2007)》，将研究区的土壤侵蚀强度划分为5个等级: 微度〔＜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)〕，得到该区域土壤侵蚀强度分级图。

2 结果与分析 2.1 土壤侵蚀模数变化特征

“一江两河”地区1995，2005，2015和2020年的平均土壤侵蚀模数分别为30.32, 31.66, 29.79, 29.64 t/(hm²·a)，呈先增加后降低的趋势(图 6), 26 a间土壤侵蚀模数降低了2.24%。1995—2005年时段内，土壤侵蚀模数增幅为4.40%;2005年以后逐渐降低，2005—2015年时段土壤侵蚀模数降幅为5.10%，而2015—2020年降幅仅为0.48%，与前一时段差异较小。拉萨河地区1995，2005，2015和2020年平均土壤侵蚀模数分别为24.89, 25.80, 23.86, 23.25 t/(hm²·a)，整体趋于稳定，26 a间土壤侵蚀模数减少了6.56%。其中，1995—2005年时段内也呈增加趋势，增幅为3.65%;2005—2015年和2015—2020年时段呈降低趋势，降幅分别为7.49%和2.55%。年楚河地区分别为45.11, 52.11, 47.88, 41.70 t/(hm²·a), 26 a间土壤侵蚀模数降低了7.56%。其年际间变化幅度较大，1995—2005年时段内土壤侵蚀模数增加了15.53%;2005—2015年和2015—2020年时段则分别降低了8.13%和12.91%。总体上年楚河土壤侵蚀模数高于拉萨河地区，且年际间波动也较拉萨河地区剧烈。

土壤侵蚀模数较高的区域主要集中在日喀则市桑珠孜区和南木林县，1995—2020年研究时段内平均侵蚀模数分别为56.41, 44.99 t/(hm²·a); 而土壤侵蚀模数降低的区域主要分布在达孜县和桑日县，1995—2020年研究时段内平均侵蚀模数分别为12.69, 17.10 t/(hm²·a)。

图 7为不同时段的土壤侵蚀模数差值图。“一江两河”地区1995—2005年、2005—2015年和2015—2020年土壤侵蚀模数降低的区域面积比例分别为65.44%, 40.86%和54.03%，侵蚀模数增加的区域面积比例分别为34.56%, 59.14%和45.97%。由此可知，研究时段内土壤侵蚀模数增加的面积呈先增加后减少的趋势，而侵蚀模数降低的面积呈现先降低后增加的趋势。1995—2005年，土壤侵蚀模数增加区域主要集中于研究区西部日喀则地区和东部雅鲁藏布江、拉萨河等河谷地带，而减小区域主要集中于东部的拉萨和山南地区。侵蚀模数增加大于10 t/(hm²·a)的区域面积比例为21.19%，而侵蚀模数降幅大于10 t/(hm²·a) 的区域面积比例为31.60%;两者均集中分布于研究区西部的谢通门县、拉孜县、南木林县和桑珠孜区，说明该区域土壤侵蚀模数在此10 a期间变化较为剧烈。此外，雅鲁藏布江、年楚河及拉萨河南段的河谷地带土壤侵蚀模数增加也较为剧烈，其中年楚河河谷侵蚀模数增加最为剧烈。2005—2015年，土壤侵蚀模数增减区域的空间分布特征与1995—2005年期间相比产生了较大变化。增降幅小于2 t/(hm²·a)的区域面积比例最大，为27.53%，且集中分布于雅鲁藏布江、拉萨河和年楚河河谷地带。而侵蚀模数增幅较大的〔增幅>10 t/(hm²·a)〕区域面积比例为8.42%，主要分布于谢通门县东部和南木林县。侵蚀模数降幅较大〔降幅>10 t/(hm²·a)〕的区域面积比例13.58%，主要分布于拉孜县、桑珠孜区、江孜县和白朗县; 除此之外，雅鲁藏布江、年楚河及拉萨河河谷地带的土壤侵蚀模数减少也较为明显。2015—2020年，侵蚀模数变化幅度较小的〔降幅＜2 t/(hm²·a)〕区域面积比例28.17%，在全域均有分布。降幅较大的区域面积比例10.85%，主要分布于研究区西部，尤以雅鲁藏布江和年楚河河谷区最甚。侵蚀模数增幅较大的区域面积比例12.04%，主要分布于谢通门县和南木林县北部。侵蚀模数降幅较大的区域面积比例11.68%，主要分布于南木林县、桑珠孜区和白朗县北部。

总体而言，研究区西部的日喀则地区土壤侵蚀模数在26 a间得到了显著降低，拉萨和山南地区土壤侵蚀模数在2005年之后变化幅度较小。

2.2 土壤侵蚀强度变化特征

1995—2020年土壤侵蚀各强度等级情况如图 8所示。不同年份的土壤侵蚀均以微度侵蚀和轻度侵蚀为主，两者面积比例总和可达70%以上。其中，微度侵蚀面积比例由1995年的21.80%增加至2020年的42.91%，轻度侵蚀面积比例则由35.52%减少至26.16%。对于微度侵蚀而言，4个研究年份2~5 t/(hm²·a) 区间微度侵蚀面积占微度侵蚀总面积比例分别为36.02%, 26.85%, 28.51%和27.31%，且随时间呈降低趋势。但是强烈及以上土壤侵蚀面积比例总和也较高，且年际间差异较小，分别为19.75%, 19.66%, 18.62%和18.85%。1995—2020年中度侵蚀和强烈侵蚀面积比例均呈降低趋势，分别由1995年的22.93%和12.01%降低至2020年的12.09%和7.25%。极强烈侵蚀面积比例年际间较为接近，但呈先增加后降低的总体趋势，由1995年的6.46%增加至2005年的7.74%，随后又降低至2020年的7.41%。剧烈侵蚀面积在1995—2005年变化较大，其面积比例由1.28%增加至4.89%, 2005—2020年呈缓慢减少趋势，至2020年面积比例降为4.19%。

1995—2020年土壤侵蚀强度均呈现西高东低的空间分布特征，总体而言日喀则地区土壤侵蚀强度远高于拉萨市和山南市(图 9)。“一江两河”地区微度侵蚀面积比例最高的区域分布于研究区达孜县、桑日县和墨竹工卡县，且面积随时间均呈增加趋势，1995—2020年分别由30.03%, 32.35%和25.93%增加至73.19%, 66.68%和66.51%。极强度侵蚀面积比例最大的区域主要分布于桑珠孜区、南木林县、江孜县和白朗县。其中极强烈侵蚀面积平均比例分别为13.39%, 11.66%, 10.43%和10.08%，且随时间呈降低趋势，分别由1995年的15.93%, 12.77%, 10.04%和10.27%降低至2020年12.29%, 10.70%, 9.18%和9.16%。剧烈侵蚀面积比例较大的区域主要分布于桑珠孜区和南木林县，平均比例分别为11.39%和7.45%;其中，桑珠孜区随时间呈增加后降低的趋势，南木林县随时间呈增加趋势。桑珠孜区1995—2005年剧烈侵蚀面积比例由9.84%增加至13.64%，至2020年降低至9.95%;南木林县剧烈侵蚀面积比例则由1995年的2.04%逐渐增加至2020年的9.33%。总体而言，“一江两河”地区强烈及以上的土壤侵蚀主要发生在西部，轻度及以下的土壤侵蚀主要发生在东部。剧烈侵蚀主要集中于桑珠孜区，其次为南木林县。而东部的达孜县、墨竹工卡县和桑日县土壤侵蚀主要以微度和轻度侵蚀为主，土壤侵蚀风险较小。

拉萨河地区土壤侵蚀强度以微度和轻度侵蚀为主。其中，微度侵蚀面积随时间呈显著增加趋势，由1995年的24.32%上升至2020年的61.46%。中度侵蚀面积比例由23.30%降低至6.81%。强烈及以上侵蚀面积比例较小，面积合计比例由1995年的14.33%降低至2020年的14.19%。年楚河地区以微度、轻度和中度侵蚀为主。其中，微度和轻度侵蚀面积呈增加趋势，分别由1995年的17.95%和21.83%增加至2020年的29.83%和27.38%; 但是中度侵蚀面积呈降低趋势，面积比例由23.39%减少至15.09%。强烈及以上侵蚀面积比例较大且呈降低趋势，由1995年36.83%增加至2020年的27.71%，主要以极强烈侵蚀面积增加为主。

1995—2020年各等级土壤侵蚀面积转移流向如图 10所示。1995—2005年各等级土壤侵蚀面积转移的主体为轻度侵蚀→微度侵蚀，其次为轻度侵蚀→中度侵蚀、中度侵蚀→轻度侵蚀/强烈侵蚀、强烈侵蚀→极强烈侵蚀/剧烈侵蚀。上述10 a为“一江两河”地区土壤侵蚀变化最复杂剧烈的时期。2005—2015年主要表现为微度侵蚀→轻度侵蚀、轻度侵蚀→微度/中度侵蚀、中度侵蚀→轻度侵蚀。2015—2020年主要表现为轻度侵蚀→微度侵蚀、中度侵蚀→轻度侵蚀。

2.3 土壤流失面积变化特征

“一江两河”地区1995，2005，2015和2020年水土流失面积比例分别为78.20%, 57.52%, 58.83%和57.09%，总体上呈逐渐减小的趋势(图 11)。其中拉萨河地区水土流失面积比例分别为75.68%, 38.91%, 40.68%和38.54%; 年楚河地区分别为82.05%, 74.55%, 73.66%和70.17%。拉萨河和年楚河地区水土流失面积也随时间呈减小趋势，且年楚河地区的水土流失面积高于拉萨河地区。在空间分布上，白朗县、江孜县、南木林县、谢通门县和桑珠孜区的水土流失面积最高，是侵蚀风险较高的区域; 26 a间平均水土流失面积比例分别为73.88%, 73.34%, 73.15%, 72.13%和72.02%。5个区县的水土流失面积比例随时间均呈降低趋势，南木林县降幅最大(23.28%)，谢通门县降幅最小(8.10%)。位于研究区东部的达孜县、桑日县和墨竹工卡县的水土流失面积较小，土壤侵蚀风险较小。1995—2020年平均水土流失面积比例分别为38.50%, 42.32%和43.61%，且随时间均呈降低趋势，其降幅分别为60.68%, 50.74%和54.79%，是土壤侵蚀改善最为显著的区县。

3 讨论

“一江两河”地区是西藏自治区的经济、文化中心，也是重要的粮食产区，但该区地形地貌复杂多变、水土流失严重; 加之气候变化和人类活动的加强，使其成为土壤侵蚀潜在变化的核心关注区之一。本研究基于RUSLE模型评估了不同时期“一江两河”地区的土壤侵蚀状况，结果显示，该区的土壤侵蚀强度呈现西高东低的空间分布特征，且主要发生在雅鲁藏布江、拉萨河和年楚河等河谷地带，土壤侵蚀风险较高，这与张鹏等^[5]、Wang等^[8]的研究结果一致。Wang等^[8]估算的2001—2015年“一江两河”地区的平均土壤侵蚀模数为29.1 t/(hm²·a)，轻度侵蚀以上水土流失面积比例63.62%。本研究1995—2020年平均土壤侵蚀模数为30.35 t/(hm²·a)，轻度侵蚀以上水土流失面积比例62.91%，与Wang等^[8]的研究结果较为接近。对其子区域而言，方广玲等^[9]估算的拉萨河流域2010年土壤侵蚀模数为30.76 t/(hm²·a)，且以中度侵蚀为主。本研究针对拉萨河2005和2015年估算的土壤侵蚀模数分别为25.80和23.86 t/(hm²·a)，且以微度侵蚀为主。本研究结果略低于上述结果，其原因有可能是研究的年份不同和土壤侵蚀因子计算方法及数据源不同。由于高寒的特殊自然环境，青藏高原的气象、水文、土壤侵蚀等监测资料较为稀缺，因此本研究未对RUSLE模型估算结果进行验证，在后续的研究中还需对土壤侵蚀模型模拟结果进行必要的精度验证。与此同时，采用不同的土壤侵蚀模型和不同的数据源，对土壤侵蚀强度估算结果影响较大。因此有必要针对数据源精度、模型估算及结果校验等方面进行系统研究，为开发适用于青藏高原的土壤侵蚀预测模型提供技术支撑。

总体而言，“一江两河”地区水土流失面积呈降低趋势，但其侵蚀模数在26 a间却表现出先增加后减小的趋势。土壤侵蚀面积和侵蚀模数表现出降低趋势的原因可能是由于青藏高原地区实施了一系列生态环境保护工程，例如保护湿地、坡改梯、植树造林、引水灌溉、生态移民等，逐步恢复了生态系统的结构与功能，提高生态系统的稳定性，降低了土壤侵蚀模数和土壤侵蚀面积^[35-36]。其中1995—2005年土壤侵蚀模数表现出增加的趋势，这可能是1995—2005年人类活动强度增加所致。有研究表明，1995—2005年，该地区灯光指数增加迅速，而灯光指数能够在一定程度上反映人类活动的强弱。人类活动的强度快速增加可能是导致这一时期土壤侵蚀模数增加的原因^[7]。

4 结论

(1) “一江两河”地区土壤侵蚀强度呈西高东低的趋势，整体以微度和轻度侵蚀为主，但是研究区西部以强烈及以上侵蚀为主。

(2) 1995—2020年，“一江两河”地区土壤侵蚀模数呈先增大后降低的趋势。1995—2005年土壤侵蚀变化最剧烈，2005—2020年土壤侵蚀呈稳步降低趋势。

(3) “一江两河”地区水土流失面积在26 a间显著降低，土壤侵蚀状况得到明显改善。拉萨河地区土壤侵蚀状况整体优于年楚河地区，雅鲁藏布江和年楚河河谷区土壤侵蚀状况改善最为显著; 谢通门县、南木林县、白朗县和江孜县等日喀则市所辖区县的土壤侵蚀风险较高，土壤侵蚀防治亟需加强。