土壤侵蚀是土体在特定时间和空间条件下的迁移过程，是最重要的土地退化问题，它直接影响水、土资源的利用与保护，危害区域生态环境，妨碍区域社会经济和可持续发展，是区域水土保持、荒漠化防治以及生态保护领域的关键科学问题之一^[1]。在全球范围内，因水力剥蚀、搬运作用形成的水蚀面积达1.09×10⁸ km²，因气流冲击、风力剥蚀形成的风蚀侵蚀面积近5.78×10⁷ km²，不考虑冻融侵蚀以及其他形式的侵蚀，全球土壤侵蚀面积已达1.64×10⁸ km²，而中国因各种因素造成不同形式的土壤侵蚀遍布全国、态势严重^[2]。湟水河流域作为黄河源区重要的生态屏障，地处高寒生态脆弱区^[3]，随着全球气候及环境变化，流域内土壤出现了不同程度的退化及侵蚀现象。因此，开展湟水河流域尺度的土壤侵蚀时空变化研究，对黄河上游地区水土流失治理具有重要的意义。

近几十年来，国内外学者在土壤侵蚀方面发展了通用土壤流失方程USLE(universal soil loss equation)、分布式土壤水力侵蚀预报模型WEPP(water erosion prediction project)等区域土壤侵蚀定量分析方法，在不同地区开展了土壤侵蚀调查、评估、预测等方面的研究^[4-6]；同时，在土壤侵蚀动力机制与过程模拟，土壤侵蚀与物质迁移，土壤侵蚀与气候变化等方面取得了进展^[7-9]，发现气候变化可通过改变降水径流、地表植被覆盖度以及人类活动影响侵蚀过程，未来土壤侵蚀的变化可通过气候变化模式与土壤侵蚀模型的耦合来进行预测^[10-11]，而在激光雷达和高时空分辨率遥感等技术的支持下，能够将地表信息精细化表达，进一步增强土壤侵蚀的模拟^[12]。另外，为适应不同区域、不同尺度的土壤侵蚀预报，构建了中国土壤流失方程CSLE(Chinese soil loss equation)^[13]、AGNPS(agricultural non-point source)模型^[14]和坡面尺度的水蚀预报模型^[15]，等。修正的通用土壤流失方程RUSLE(revised universal soil loss equation)在不同地区土壤侵蚀研究中也有应用，开展了诸如土壤侵蚀量估算、土壤侵蚀强度分析、土壤侵蚀时空演变及其空间分异特征等方面的研究^[16-18]。但是，在土壤侵蚀模型模拟结果的检验、土壤侵蚀动力机制等方面，特别是应用先进的对地观测技术开展实地调查与验证工作有待深入研究。本研究基于2000年和2018年的MODIS，Landsat，降雨、人口密度、经济等多源数据，利用修正的通用土壤流失方程RUSLE，分析湟水河流域土壤侵蚀的特点与时空变化规律，并基于无人机遥感技术对结果进行实地验证与评价，以进一步明确湟水河流域近20 a来的土壤侵蚀时空变化特征及驱动因素，以期为黄河上游地区的水土保持与防治工作提供数据基础与决策依据。

1 研究区概况

湟水河是黄河流域重要的水源地之一，也是黄河上游最大的一级支流，处于青藏高原与黄土高原的交错地带，位于北纬36°05′—38°30′，东经99°00′— 103°43′，海拔1 545~5 232 m，年均气温5.2~6.8 ℃，年降水量350~550 mm，属高海拔半干旱地区，地貌过渡性和复杂性明显，流域水系呈树叶状分布，干流南北两岸支流发育^[19]。湟水河源区多沼泽湿地，流经地区多峡谷、盆地相间分布，存在高寒草甸地貌，分布优良牧场和良好草地；中游一带平坦开阔，存在有较大面积的灌溉农田和非灌溉农田分布；中下游主要为石林山区、土石山区、黄土丘陵沟壑区；流域东南部多分布寒钙土、新积土和冷钙土。流域西部的山地以灰褐土和黑钙土为主，流域中部的中高山地以黑钙土、黄土和壤土为主。研究区总人口约为3.30×10⁷人，农田耕地超过3.00×10⁶ hm²，是青海省主要的人口居住区和产粮区^[20-21]。

2 数据来源与研究方法 2.1 数据来源

本研究主要使用MODIS，Landsat，数字高程数据DEM、降雨、人口密度和经济等数据(产品)。

MODIS数据来源于美国航空航天局数据中心(https://www.earthdata.nasa.gov/)，选用2000，2018年行列号为h26 v05和h25 v05的MODIS MOD13Q1产品数据，空间分辨率为250 m，采用最大值合成法获得两个年度的NDVI数据；

Landsat数据来源于地理空间数据云(http://www.gscloud.cn/)，选取2000，2018年4月至10月轨道号为131—134，行编号为33—35，且云量小于5%的7幅影像，空间分辨率为30 m，对其进行辐射定标、几何纠正、镶嵌、分类等处理，获得两个年度的土地覆盖数据；DEM数据来源地理空间数据云(http://www.gscloud.cn/)，空间分辨率为30 m，对该数据进行流域分析，获得研究区流域边界与河网信息；2000，2018年研究区的降雨数据来源于中国气象数据网(http://www.data.cma.cn)，人口密度和经济数据来源于中国资源环境科学与数据中心(http://www.resdc.cn)，对上述数据进行投影转换，并利用流域边界数据对其进行裁切，获得研究区土壤侵蚀分析的数据集；

无人机遥感数据是利用大疆精灵4 RTK于2019年8月分别在研究区现场采集，设置航向重叠为80%，旁向重叠为70%，航高为120 m。将上述数据的空间分辨率统一重采样至250 m，均使用Albers Conical Equal Area投影。

2.2 研究方法 2.2.1 土壤侵蚀模型

本研究选用修正的通用土壤流失方程RUSLE计算研究区土壤侵蚀模数，公式如下^[22]：

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

(1)

式中：A为年均土壤侵蚀模数〔t/(hm²·a)〕; R为降雨侵蚀力因子〔(MJ·mm)/(hm²·h·a)〕; K为土壤可蚀性因子〔(t·h)/(MJ·mm)〕; LS为地形(坡长坡度)因子, 无量纲; C为植被覆盖和作物管理因子, 无量纲; P为水土保持措施因子, 无量纲。

2.2.2 降雨侵蚀力因子R值的确定

根据己有湟水河流域地区月度气象数据资料和模型的适用性，本研究采用Wischmeier等^[23]提出的R值计算方法，公式如下：

$ R = \sum\limits_{i = 1}^{12} 1 .735 \times {10^{〔 {1.5\lg \frac{{P{i^2}}}{P} - 0.818\;8} 〕}} $

(2)

式中：R为降雨侵蚀力值; P_i为月均降雨量(mm); P为年均降雨量。

根据青海省2000年和2018年22个气象站点的降雨数据，计算得到各个气象站的R值，经克里金空间插值计算得到研究区降雨侵蚀力因子R值分布，见图 1。

2.2.3 土壤可蚀性因子K值的确定

土壤可蚀性因子K与土壤粉粒、黏粒等含量有关，其计算方法如下^[24]：

$ \begin{aligned} K_{a}=&\left\{0.2+0.3 \exp \left[-0.0256 \mathrm{SAN} \times\left(1-\frac{\mathrm{SIL}}{100}\right)\right]\right\} \times\left(\frac{\mathrm{SIL}}{\mathrm{CLA}+\mathrm{SIL}}\right)^{0.3} \times \\ & {\left[1.0-\frac{0.25 C}{C+\exp (3.72-2.95 C)}\right] \times\left[1.0-\frac{0.7 \times \mathrm{SN}_{1}}{\mathrm{SN}_{1}+\exp \left(-5.51+22.9 \times \mathrm{SN}_{1}\right)}\right] } \end{aligned} $

(3)

其中：

$ \mathrm{SN}_{1}=1-\frac{\mathrm{SAN}}{100} $

(4)

式中：SAN为砂粒含量(%); SIL为粉粒含量(%); CLA为黏粒含量(%); C为土壤有机碳含量(%); K_a为土壤可蚀性因子。利用上述公式计算获得研究区土壤可蚀性因子的空间分布(图 2)。

2.2.4 坡度坡长因子LS值的确定

基于DEM数据计算研究区坡度、坡长数据，进而通过分段函数计算坡度因子S值和坡长因子L值，进而计算坡度坡长因子LS，公式如下^[2]：

$ L=\left(\frac{\lambda}{22.13}\right)^{m} $

(5)

$ m= \begin{cases}0.5 & (S \geqslant 5 \%) \\ 0.4 & (3 \%<S<5 \%) \\ 0.3 & (1 \%<S<3 \%) \\ 0.2 & (S \leqslant 1 \%)\end{cases} $

(6)

$ S=\left\{\begin{array}{cl} 10.8 \sin \theta+0.03 & \left(\theta<5^{\circ}\right) \\ 16.8 \sin \theta-0.5 & \left(5^{\circ} \leqslant \theta \leqslant 10^{\circ}\right) \\ 21.91 \sin \theta-0.96 & \left(\theta \geqslant 10^{\circ}\right) \end{array}\right. $

(7)

$ L S=L \times S $

(8)

式中：λ为坡长(m); θ为倾斜角(°); S为百分比坡度(%)。利用上述公式计算获得研究区坡度坡长因子LS值的空间分布(图 3)。

2.2.5 植被覆盖与作物管理因子C值的确定

以获得的研究区NDVI数据为基础，经处理得到植被覆盖度后^[22]，由蔡崇法等^[25]提出的植被覆盖度F_c与作物管理因子的转换公式计算C值：

$ C=\left\{\begin{array}{cl} 1 & \left(F_{c} \leqslant 10 \%\right) \\ 0.6508-0.3436 \lg F_{c} & \left(10 \%<F_{c}<78.3 \%\right) \\ 0 & \left(F_{c} \geqslant 78.3 \%\right) \end{array}\right. $

(9)

式中：C为植被覆盖与作物管理因子; F_c为植被覆盖度, 其计算方法如下:

$ F_{c}=\frac{N-N_{\mathrm{soil}}}{N_{\mathrm{veg}}-N_{\mathrm{soil}}} $

(10)

式中：F_c为年作物生长期的植被覆盖度; N为像元NDVI值; N_soil为裸地像元NDVI值; N_veg为纯植被覆盖像元NDVI值; N_veg和N_soil分别为NDVI统计结果在累积概率为95%和5%时的取值。利用上述公式计算获得研究区2000，2018年植被覆盖与作物管理因子C值的空间分布(图 4)。

2.2.6 水土保持措施因子p值的确定

采用Lufafa等^[6]提出的方法计算水土保持措施因子p值，公式如下：

$ P=0.2+0.03 \times \theta $

(11)

式中：θ为坡度, 由此计算获得研究区水土保持措施因子p值的空间分布(图 5)。

2.2.7 土壤侵蚀强度的验证

在土壤侵蚀典型区域选择实地调查样地进行模型模拟结果的验证，无人机数据空间分辨率优于0.05 m的影像数据，以250 m的格网为单元通过目视解译法解译调查样地的土地利用^[26]、水土保持措施，并计算植被覆盖度和坡度，依据《水土保持遥感监测技术规范(SL592-2012)》《土壤侵蚀分类分级标准(SL190-2007)》确定土壤侵蚀强度，计算总体精度^[2]：

$ P=\left[\left(1-\frac{\sum\limits_{i=1}^{n}\left|\mathrm{SE}_{i}-\mathrm{SE}_{i}^{*}\right|}{n}\right)\right] \times 100 \% $

(12)

式中：P为土壤侵蚀总体精度; SE_i为模型模拟土壤侵蚀强度等级; SE_i^*为无人机地面调查的土壤侵蚀强度等级; n为调查验证的总格网数。P值越接近100%表示模型计算的结果越准确。

3 结果与分析 3.1 2000年土壤侵蚀模型计算分析

基于2000年土壤侵蚀模型因子计算结果，利用式(1)计算获得研究区2000年土壤侵蚀模数，其最大值为15 706.50 t/(km²·a)，均值为477.81 t/(km²·a)。依据《土壤侵蚀分类分级标准》和祁连山地区土壤侵蚀的相关研究结果^[27]，利用表 1对计算获得的土壤侵蚀模数结果进行等级划分，结果见图 6。

由图 6可知，2000年研究区整体土壤侵蚀强度不大，微度侵蚀主要分布在湟水河源区及山间河谷地区的耕地和高植被覆盖度的草地、有林地；轻度、中度侵蚀主要分布在海拔较高、植被覆盖少的山地、荒地。

3.2 2018年土壤侵蚀模型计算分析

基于2018年土壤侵蚀模型因子计算结果，利用式(1)计算获得研究区2018年土壤侵蚀模数，其最大值为32 653.50 t/(km²·a)，均值为1 625.30 t/(km²·a)，利用表 1对计算获得的土壤侵蚀模数结果进行等级划分，结果见图 7。

由图 7可知，2018年研究区主要以微度侵蚀为主，分布在湟水河源区、高山之间的河谷盆地地区，轻度侵蚀主要分布在山区及其周围的耕地、坡耕地等区域，中度侵蚀主要分布在研究区东南部城镇居民聚集地带与河流滩地，强烈侵蚀和极强烈侵蚀零散分布于祁连山、达坂山等高山、秃岭裸地地区。

3.3 土壤侵蚀模型的验证分析

由以上获得的研究区土壤侵蚀各个模型因子计算结果可知，降水侵蚀力因子R值呈东北低西南高的特点，土壤可蚀性因子K值较大区域分布在湟水河源区和中部河流两侧地区，而坡度坡长因子LS值在中部及中部偏北地区较高，植被覆盖与作物管理因子C值东南高西北低，水土保持措施因子P值在湟水河源区普遍较低，在北部祁连山脉、中部达坂山和南部拉脊山一带较高。因此，分别在研究区北部的门源、中部的大通和南部的湟中选择总面积约3 km²的样地，采集的无人机高分辨率影像数据(见图 8)，建立解译标志。

依据解译标志和土壤侵蚀强度验证方法获得49个格网的侵蚀强度，并依据公式(12)计算土壤侵蚀验证总体精度P为91.8%，表明模型计算的研究区土壤侵蚀强度可靠。

3.4 土壤侵蚀强度时空变化分析

统计研究区2000年和2018年各个侵蚀强度等级的面积，绘制变化对比图(图 9)，发现微度侵蚀面积比例显著减少，中度、强烈、极强烈侵蚀面积显著增加，而轻度、剧烈侵蚀的面积虽有增加，但幅度不大，表明研究区土壤侵蚀面积有向中度、强烈和极强烈变化的趋势。对2000年和2018年研究区土壤侵蚀强度进行逐像元求差，得到0，1，2，4的运算结果，分别代表未发生变化、微度、中度和剧烈变化像元，空间分布如图 10所示。

由图 10可知，在湟水河源区、祁连山与达坂山之间的支流大通河谷地等区域土壤侵蚀强度未发生明显变化，在沿湟水河干流中下游的城镇居民聚集地带存在土壤侵蚀强度微度变化，面积约占28.12%；在祁连山、达坂山的高山裸地、秃岭地区存在土壤侵蚀强度中度变化，面积约占8.91%，而土壤侵蚀强度剧烈变化区域则零散分布在高山裸地及城镇地带。

4 讨论与结论 4.1 讨论

针对湟水河流域及其周边地区的土壤侵蚀过程与机制问题，已有学者开展了一定程度的研究。如王雪璐^[28]利用RUSLE模型估算青藏高原三江源区多年平均侵蚀量为3.10×10⁹ t/a；陈豪等^[24]利用USLE模型对祁连山公园的土壤侵蚀进行了模拟计算，得到2015年和2019年的平均侵蚀模数为1 063 t/(km²·a)和1 068 t/(km²·a)；徐进等^[29]对湟水河(西宁段)流域的土壤侵蚀危险性进行了综合评价，发现土壤侵蚀危险性较高的地方主要集中在城镇居民地聚集的市区及湟水河河谷地带；魏素娟等^[1]基于USLE模型对湟水流域土壤侵蚀进行了评估，得到湟水流域土壤侵蚀模数为523.8 t/(km²·a)，年侵蚀总量为8.37×10⁷ t/a。这些结果与本研究得到的湟水河流域土壤侵蚀模数和土壤侵蚀量基本一致，结果可靠。同时，本研究还利用低空无人机遥感技术对模型计算的土壤侵蚀强度结果进行了定量评价，具有一定的创新。

2000—2018年，研究区土壤侵蚀量增加，发生中、重度变化的区域主要集中在祁连山、达坂山的高山裸地、秃岭地区，沿湟水河干流的城镇居民地聚集地带的侵蚀强度也存在加剧的趋势。从本研究采用的RUSLE模型来看，土壤侵蚀过程是气候、植被覆盖等相对动态因子和地形、土壤等相对静态因子共同驱动的结果。相对静态因子主要控制土壤侵蚀的空间分布，不同的土壤类型孕育不同的植被，形成地域性特征；而相对动态因子控制土壤侵蚀的时空变化，即植被、气候等因子对土壤侵蚀起降低或加剧作用。根据研究区气象、人口、经济等数据，发现19 a来研究区平均降水量增加了113.12 mm，建设用地增加了5.25 km²，人口密度增加了30人/km²，而人口密度和建设用地的增加在一定程度上影响了地表植被覆盖度，再加上气候因子的持续性改变，导致相对动态因子逐步驱动研究区土壤侵蚀强度加剧。由此表明，研究区土壤侵蚀强度的增加，与气候暖湿化、人类活动强度增加存在一定的关系。

2000—2018年，湟水河源区、大通河谷地等区域土壤侵蚀强度未发生明显变化，在中下游城镇居民聚集地带发生微度变化，而高山裸地、秃岭地区土壤侵蚀强度发生中度、强烈的变化，存在显著的空间差异性。未发生明显变化的区域人口密度低、人类活动以放牧为主，植被覆盖稳定，植被控制由于降雨量增加侵蚀强度变大的作用显著；位于研究区中下游的城镇居民区聚集有青海省60%的人口，但在生态工程建设的有效驱动下，增加了地表的植被覆盖度和水土保持功能，在一定程度上阻滞了由于人类活动加剧和气候变化带来的侵蚀强度增大；而裸地、秃岭地区由于缺乏有针对性的水保措施，土壤与植被相互影响和作用，出现侵蚀加剧与覆盖度逐步降低的不良循环。

由于数据源时空分辨率的差异，导致模型计算结果的空间分辨率较低，且气象站数据较少造成降雨侵蚀力因子的精度也不高，下一步将结合土壤侵蚀定位观测试验进行更加精细的研究，并开展湟水河流域土壤侵蚀驱动机制等方面的研究。

4.2 结论

(1) 2000年研究区整体土壤侵蚀强度不大，以微度侵蚀为主，在祁连山、达坂山及拉脊山海拔较高、植被覆盖较少的地区存在中度侵蚀；2018年研究区中度侵蚀主要分布在研究区东南部城镇居民聚集地带与河流滩地，强烈侵蚀和极强烈侵蚀主要分布在祁连山、达坂山等高山、秃岭裸地地区。

(2) 2000—2018年，平均土壤侵蚀模数存在增加的趋势，轻度侵蚀量占比显著减少，强烈、极强烈和剧烈侵蚀量占比显著增加，这与研究区气候暖湿化、人类活动强度增加存在一定的关系。19 a间，湟水河源区、大通河谷地等区域土壤侵蚀强度未发生明显变化，祁连山、达坂山的高山裸地、秃岭地区的土壤侵蚀发生中、重度变化，沿湟水河干流的城镇居民地聚集地带的侵蚀强度也存在加剧的趋势，存在较为明显的空间差异性。这种差异性受区域人类活动特点、植被覆盖情况、生态工程建设、自然地理特性等因素的影响。