祁连山南坡土壤侵蚀时空变化及迁移特征分析

童 珊, 曹广超, 张 卓, 闫 欣, 刁二龙

(1.青海师范大学 地理科学学院, 西宁 810008; 2.青海省自然地理与环境过程重点实验室,西宁 810008; 3.青藏高原地表过程与生态保育教育部重点实验室, 西宁 810008)

土壤侵蚀是世界性的环境灾害之一[1-2]，不仅会引起土壤退化，土壤肥力下降，更能将土壤中的氮、磷、钾等营养元素及有害物质进入河道，进而造成水体污染[3]，且土壤一旦遭受破坏在短时间内难以修复，因此加强地区水土保持工作，降低土壤侵蚀对人类生产生活造成危害，成为相关学者关注的重要问题。修正后的水土流失方程(RUSLE)是现今应用最广泛的土壤侵蚀模型，模型结构清晰、简单易懂。随着GIS与RS技术的发展，国内外许多学者将土壤侵蚀模型与地理信息系统相结合，对土壤侵蚀进行定量化研究[4-6]。Borrelli等[7]利用RUSLE模型对意大利山间盆地的土壤侵蚀进行研究；陈峰等[8]基于RUSLE模型和GIS/RS空间信息技术，对滇南山区土壤侵蚀时空变化进项研究；林锦阔[9]基于RUSLE和TLSD模型，构建适用于河西地区的水蚀过程模型；康琳琦等[10]采用RUSLE模型分析青藏高原1984—2013年的土壤侵蚀情况；林慧龙等[11]采用RUSLE模型结合137Cs示踪法对三江源土壤侵蚀空间分布及影响因子进行研究。

随机森林对异常值的容忍性较强[12]，且具备较强的非线性挖掘能力[13]，同时解决了因子共线性问题，可以较好地刻画变量之间的非线性关系，因此被广泛应用。郑利林等[14]利用随机森林模型预测虉草和南荻植被分布，预测精度达89.6%，89.3%；包青岭等[15]利用随机森林模型对土壤有机质进行反演，得出随机森林在预测上相较于灰色关联分析表现出更优的可靠性与稳定性；卢宏亮等[16]利用随机森林模型对安徽省的土壤属性分布进行预测，得出利用该模型能够反映大尺度区域上土壤属性的空间分布；重心迁移模型最早应用于人口及经济方面，指事物的集中点或平衡点[17]，或指某事物的发展方向，本文利用重心迁移模型探究土壤侵蚀的变化及移动方向。赵珍珍等[18]利用重心迁移模型对科尔沁沙地土地利用变化进行研究，范俊甫等[19]利用重心迁移模型对黄河流域1980—2015年降雨侵蚀力变化进行研究。

祁连山南坡气候独特，生态环境脆弱且主要以山地地貌为主，是青海省主要的原始森林分布区和水源涵养区，维系着青藏高原生态平衡和西部区域生态安全[20-21]。然而根据青海省2019年水土保持公报，祁连山—黑河国家级水土流失预防区水土流面积为4 422.12 km2，约占全省的土壤总面积的21.68%，全部为水力侵蚀，水土流失较为严重，因此对水土流失的保护已刻不容缓。本文通过对研究区近20 a来土壤侵蚀时空变化及侵蚀重心迁移的定量研究，为研究区生态环境保护提供理论基础。

1 材料与方法

1.1 研究区概况

祁连山南坡研究区地处祁连山中段腹地，位于青藏高原东北缘，黄土高原西缘，为黄土高原向青藏高原的过渡地带，是青藏高原特征和黄土高原特征的综合体。地理位置处于东经 98°08′13″—102°38′16″，北纬37°03′17″—39°05′56″，总面积2.4万km2，研究区分布着黑河、大通河、石羊河等河流，山地由丘陵、低山、中山和高山4大类型共同组成。研究区行政区分为祁连县、门源县、刚察县、天峻县及海晏县。年均温为-5.9℃，年极端最高气温30.5℃，极端最低气温-37.1℃，年平均降水在400 mm左右，研究区主要以草地为主，林地与耕地面积较小，不同土地利用类型空间分布见图1。

图1 土地利用类型

1.2 数据与方法

1.2.1 数据来源 气象数据来源于中国气象科学数据共享服务网(http:∥cdc.cma.gov.cn)下载国家台站监测数据，包括祁连山及周边气象站点的日降雨资料，共33个站点，整理后得到站点的月降雨量及年降雨量数据。

本文使用2000—2019年的植被遥感数据MOD13Q1，来源于NASA网站(https:∥ladsweb.modaps.eosdis.nasa.gov/)，空间分辨率为250 m×250 m，时间分辨率为16 d，经过拼接、裁剪及最大合成法形成年NDVI数据。

DEM数据来源于地理空间数据云平台(http:∥www.gscloud.cn)，空间分辨率为90 m×90 m。

土地利用类型数据来源于中国科学院地理科学与资源研究所(http:∥www.resdc.cn/)，时间为2000年、2005年、2010年、2015年及2020年，空间分辨率为30 m×30 m。

1.2.2 研究方法

(1) 本文使用RUSLE土壤流失方程，对研究区长时期土壤侵蚀进行量化评估，表征降雨和下垫面共同作用下的单位面积潜在土壤流失速率[22]:

A=R·K·LS·C·P

(1)

式中:A为土壤侵蚀模数[t/(hm2·a)]；R为降雨侵蚀力因子[(MJ·mm)/(hm2·a)]；K为土壤可蚀性因子[(t·hm2·h)/(hm2·hm2)]；LS为坡长坡度因子,无量纲；C为植被覆盖管理因子,无量纲,值域范围0～1；水土保持措施因子P与土地利用类型相关，无量纲，值域范围0～1[2,10,23]，最终将结果划分为5个等级[23]。

本文运用蔡崇法等建立的植被覆盖管理因子(C)与植被覆盖度的关系计算公式得到各个时期的植被覆盖管理因子栅格数据。修正后的公式如下[24]：

(2)

(3)

式中:C为植被盖度因子;c为植被覆盖度。

R反映降雨对土壤侵蚀的影响，利用青海省分布的气象站点降雨观测数据，整理得到研究区月平均降雨量和年平均降雨量，计算出各站点降雨侵蚀力后再通过随机森林插值得到整个研究区降雨侵蚀力因子。

(4)

式中：R为降雨侵蚀力；pi为月均降雨量(mm)；p为年平均降雨量(mm)。

土壤可蚀性K值大小表示土壤是否受侵蚀破坏的性能，是控制土壤承受降雨和径流分离及输移等过程的综合效应[25]，其计算公式如下：

SN=1-SAN/100

(5)

(6)

式中：SAN,SIL,CLA和C分别为砂粒含量(%)、粉粒含量(%)、黏粒含量(%)和有机碳含量(%)。

坡长坡度因子LS以DEM为基础数据，采用从国家科技基础条件平台—国家地球系统科学数据共享服务平台—黄土高原科学数据中心(http:∥loess. geodata.cn)申请得到的LaunchLS工具计算研究区LS因子[26]。

(2) 土壤侵蚀量计算。

(7)

式中:A为单位面积上土壤侵蚀量;mi为第i个像元的面积；i为像元。

(3) 重心迁移模型。为了解土壤侵蚀的空间聚集及迁移特征，本研究采用重心迁移模型对土壤侵蚀进行定量化表达，其公式如下[27-28]：

(8)

(9)

式中:X和Y分别为重心经、纬度;n为研究区栅格数量;i为栅格序号;Xi,Yi分别为第i个栅格的中心经、纬度;Ai为第i个栅格的面积;Wi为第i个栅格的植被覆盖度值。

重心移动距离公式如下：

(10)

式中:d为第t+1年相对于第t年重心移动的距离;常数C为由地球表面的坐标(°)转化为平面距离(km)的系数,其值为111.111[29-30]。

(4) 随机森林。随机森林(RF)是一种统计学习理论，该算法是由Breiman[31]在2001年提出的用于分类和回归的一种机器学习方法，大多数应用于处理非线性问题[32]。RF算法通过R语言中的Random Forest包实现。RF是由决策树构成的森林[33]。有2个关键参数：ntree和mtry，ntree代表决策树数量，mtry代表随机变量的数量。利用RF进行回归运算时，最终的结果为所有决策树预测结果的平均值。

2 结果与分析

2.1 降雨侵蚀力插值结果验证

运用随机森林方法，以海拔作为影响因子，对青海省33个气象站点的降雨侵蚀力进行插值，得到的插值验证结果如图2所示，精度均在0.8以上，插值结果较好。

图2 降雨侵蚀力插值结果验证

2.2 祁连山南坡土壤侵蚀时空变化特征

研究区整体上土壤侵蚀模数呈现出西北向东南递减的趋势，平均值为2 598.53 t/(km2·a)，从不同年份看，土壤侵蚀模数的平均值由2000年的1 966.63 t/(km2·a)增加至2005年3 228.51 t/(km2·a)，最后下降至2019年的2 299.06 t/(km2·a)。根据《土壤侵蚀强度分类分级标准(SL190—2007)》，将祁连山南坡土壤侵蚀强度划分为微度、轻度、中度、强度、极强度及剧烈6个级别，见图3。2000—2005年侵蚀变化明显，极强度侵蚀突显，2010—2015年微度、轻度侵蚀逐渐向中度侵蚀转变。剧烈侵蚀与极强度侵蚀零散分布在研究区海拔4 500 m左右的裸土裸岩区，强度、中度及轻度侵蚀主要分布在祁连县、门源县西北部及天峻县西部，微度侵蚀主要分布在门源县、海晏县及刚察县海拔3 500 m以下的区域内。

研究区2000—2005年期间土壤侵蚀变化明显，微度、轻度及中度侵蚀面积均有不同程度的减少，微度侵蚀面积减少最大，面积减少662.69 km2，占比减少2.76%，剧烈侵蚀面积增加最大，面积增加544.99 km2，占比增加2.27%；2005—2010年及2015—2019年，除微度侵蚀面积增加，其余侵蚀等级面积均是减少；而2010—2015年，微度、轻度及中度侵蚀面积增加，其余侵蚀等级面积减少。从强度、极强度及剧烈侵蚀在2005年侵蚀面积增大外，其余年份均在减少，而微度、轻度及中度侵蚀面积变化不稳定(表1)。

2.3 不同行政单元土壤侵蚀变化特征

总体上看，2015—2019年研究区土壤侵蚀程度相对减轻，2005年土壤侵蚀最为严重，祁连县、天峻县、刚察县、海晏县及门源县平均侵蚀模数分别为35.22，18.95，21.53，25.09，26.92 t/(km2·a)，其次是2010年、2015年、2019年，最后为2000年，整体上呈现出增加再减少的趋势；从表2—3来看，2000年与2015年平均侵蚀模数各县大小为祁连县>门源县>天峻县>刚察县>海晏县；2005年各县土壤侵蚀量较2000年均有不同程度的增加，2005年、2010年与2019年平均侵蚀模数大小为祁连县>门源县>海晏县>刚察县>天峻县。土壤侵蚀量的排序为祁连县>门源县>刚察县>天峻县>海晏县，相较于2005年，研究区土壤侵蚀改善明显(表2—3)。

图3 土壤侵蚀分级强度

表1 祁连山南坡 2000-2019年各侵蚀强度等级的侵蚀面积及比例

表2 祁连山南坡不同行政单元2000-2019年土壤侵蚀模数 t/(km2·a)

表3 祁连山南坡不同行政单元2000-2019年土壤侵蚀量 t/a

从图4来看，5县整体上主要以微度及轻度侵蚀为主，除剧烈侵蚀外，随着侵蚀程度的增加侵蚀模数逐渐下降，5县中海晏县土壤侵蚀面积波动最大，微度侵蚀在不同侵蚀程度中面积变化最大，由2005年的57.67%上升至2015年71.67%，其余侵蚀面积均在下降，说明海晏县政府在土壤侵蚀防治方面取得一定成绩；刚察县2005年土壤侵蚀达到最大值，是土壤侵蚀最为严重一年、2010年土壤侵蚀得到有效改善，海晏县与刚察县2015年之后土壤侵蚀有加重趋势；天峻县土壤侵蚀在2005年与2010年相对较大，2010年之后有所减小，但各等级土壤侵蚀面积变化较小，面积变化在1%以下，说明天峻县在土壤侵蚀预防方面实施了有效的措施，但在治理上稍有欠缺；祁连县的土壤侵蚀在2005年达到最大，但其强度、极强度及强烈侵蚀在2010年、2015年及2019年土壤侵蚀面积的基本为连续减少，说明土壤侵蚀得到了一定程度的改善；门源县2005年各等级土壤侵蚀达到最大，微度侵蚀面积有2005年的76.84%增加至2019年的78.22%，强度，极强及剧烈侵蚀面积也由4.11%，4.57%及4.68%降至3.87%，4.57%及2.76%，门源县政府在水土保持方面做了巨大的贡献。

2.4 主要土地利用类型土壤侵蚀变化特征

不同土地利用类型的土壤侵蚀模数存在差异，而2015年与其余年份不同的是裸岩的土壤侵蚀模数小于草地的土壤侵蚀模数，其余年份土壤侵蚀模数大小排序为永久冰川>裸岩>草地>林地>盐碱地、沙地>裸土>沼泽>耕地；土壤侵蚀量排序为：草地>裸岩>林地>永久冰川>沼泽>盐碱地、沙地>耕地>裸土(表4—5)。

图4 不同行政单元各侵蚀等级面积变化

表4 祁连山南坡不同土地利用类型2000-2019年土壤侵蚀模数 t/(km2·a)

土壤侵蚀量与土地利用类型面积有很大关系，草地是研究区最大的土地利用类型，面积占比为50%以上，其次为裸岩，面积占比为19.40%，再次为林地面积占比为17.61%，沼泽占比为8.93%；最后耕地、裸土、未利用土地面积仅占比为2.81%，永久冰川面积仅占0.29%，通过土壤侵蚀量与面积的大小，本文主要选用林地、草地、永久冰川及裸岩进行深入研究。

如图5所示，草地与林地整体上呈增加—减少—增加—减少趋势，两者均主要以微度侵蚀为主，且强烈侵蚀所占面积大于极强度及强度侵蚀。2005—2015年土壤侵蚀强度及以上等级面积变化较小，起伏不定，可见对土壤侵蚀的治理效果不显著。2015—2019年剧烈侵蚀面积下降明显，强度及极强度侵蚀也有所下降，土壤侵蚀得到改善。

表5 祁连山南坡不同土地利用类型 2000-2019年土壤侵蚀量 t/a

以2010年为界，裸岩微度、轻度、中度侵蚀总体上先减少后增加，而强度、极强及剧烈侵蚀先增加后减少，说明后10 a裸岩土壤侵蚀得到了缓解，且裸岩是微度、轻度、中度侵蚀共同作用，不同侵蚀强度所占的面积随侵蚀等级强度的增加而减少，轻度侵蚀面积变化最大，由2005年的14.19%上升至2015年的19.93%，其次为极强度侵蚀，侵蚀面积变化由2005年的10.29%下降至2015年的5.01%，强度与强烈侵蚀面积也有一定程度的下降，说明在此时间段内裸岩的土壤侵蚀有所改善。

图5 主要土地利用类型各侵蚀等级面积变化

永久冰川在2005年土壤侵蚀得到有效控制，在2015年土壤侵蚀有所加重。主要以微度侵蚀及强烈侵蚀共同作用，两级分化较为严重，强烈侵蚀面积占50%以上，除微度侵蚀外，随侵蚀的强度增加，侵蚀面积增加。剧烈侵蚀面积变化最大，由2005年的最大值59.39%下降至2015年的53.89%，其余侵蚀面积均有不同程度的增加，说明永久冰川土壤剧烈侵蚀有所改善。

2.5 土壤侵蚀重心转移特征

2.5.1 祁连山南坡土壤侵蚀重心转移特征 基于物理中心迁移模型，以5 a为间隔统计土壤侵蚀，探究土壤侵蚀的空间迁移特征。从图6可知，2000—2019年整体土壤侵蚀重心总共迁移263.09 m，迁移方向为东南—西北—东南，迁移速率为13.15 m/a，2010年之前，土壤侵蚀重心整体向东南迁移，2010—2015年土壤侵蚀重心由东南向西北迁移，2015—2019年土壤侵蚀重心又由西北向东南迁移。

2.5.2 不同行政单元土壤侵蚀重心转移特征 从图7可以看出，祁连县土壤侵蚀重心位移最大，为351.89 m，迁移方向为东南—西北—东南，迁移速率为17.59 m/a，其次为刚察县与门源县，土壤侵蚀重心分别迁移25.28，21.19 m，迁移方向分别为东南—西南—西北及西北—西南—东南，迁移速率分别为1.26，1.06 m/a，最后为天峻县与海晏县，土壤侵蚀重心分别迁移13.05，3.84 m，迁移方向分别为东北—西南—东北及东南—西南—西北，迁移速率分别为0.65，0.19 m/a。相比之下得出海晏县土壤侵蚀得到有效的控制，祁连县的水土保持措施有待加强。

图6 祁连山南坡土壤侵蚀重心转移

2.5.3 主要土地利用类型土壤侵蚀重心转移特征 从图8中可以看出，主要土地利用类型中，永久冰川土壤侵蚀重心总迁移量最大，为367.78 m，尤其在2005—2010年，土壤侵蚀重心向西北迁移358.76 m，整体迁移方向为西北—东，迁移速率为18.39 m/a，其次为草地，土壤侵蚀重心迁移量为169.11 m，迁移方向为东南—西北—东南，迁移速率为8.46 m/a，最后为裸岩、林地，迁移速率分别为1.19，0.18 m/a。

3 讨 论

祁连山南坡土壤侵蚀由西北向东南递减，这与降雨的趋势正好相反[34]，主要原因为降水量的增加与侵蚀强度的上升之间并没有显著的线性关系，雨滴溅蚀与径流剥蚀会受到其他因素的综合影响[35]，且研究区主要以微度侵蚀为主，强度及以上的侵蚀主要发生在高海拔区域，大都分布在水域及裸岩周围，发生强度降雨时，容易发生土壤侵蚀，这与康琳琦等[10]对青藏高原的研究结果一致。研究区土壤侵蚀在2005年达到最大值，可能是2000—2005年祁连山地区降雨量的增加，人类保护环境的观念较弱，导致祁连山在2005年土壤侵蚀最大。微度侵蚀主要在海拔3 500 m以下是由于低海拔区主要为建筑用地、耕地及林地，土壤侵蚀相对较小。近5 a祁连山南坡的土壤侵蚀得到改善与祁连山实施的封禁治理及修建坡改梯等生态保护与建设综合治理工程有关。不同土地利用类型及不同行政单元土壤侵蚀量的大小由各单元面积及侵蚀强度共同决定，因此与平均土壤侵蚀模数排序不一致。

图7 各县土壤侵蚀重心转移

图8 主要土地利用类型土壤侵蚀重心转移

从各县看祁连县土壤侵蚀量最大，这与祁连县面积及位置有一定的关系， 2005年、2010年、2019年海晏县的土壤侵蚀超过天峻县，这是由于海晏县降雨一直远大于天峻县(表6)，土壤侵蚀加重的风险大于天峻县。祁连县2000—2005年及2010—2015年土壤侵蚀重心迁移距离最大，分析原因可能是2000—2005年由于城镇化的发展，植被破坏土壤侵蚀增强，重心南移；2010—2015年由于人类活动稳定后，海拔高、植被覆盖低，侵蚀严重，进而重心向西北移动。其余4县土壤侵蚀重心移动相对较弱，说明了土壤侵蚀得到了一定的控制。

2015年裸岩，冰川的侵蚀量大幅减少，可能是降雨的减少及水土保持措施的实施相结合的结果，耕地土壤侵蚀最低，是由于耕地主要分布于门源县海拔较低区域，有田埂保护。草地与林地土壤侵蚀较大是因为其面积占比大以及高寒草原与草甸分布于研究区高海拔区域，降水变化敏感，且鼠害猖獗[36-37]。而2015—2019年政府调控力度及人们保护环境意识加强，土壤侵蚀得到改善。永久冰川与裸岩由于植被覆盖率小，地表径流及降雨量较大，导致土壤侵蚀最高。而冰川土壤侵蚀重心迁移距离最大，尤其是2005—2010年，向西北方向转移距离最大，原因可能是永久冰川与裸岩主要分布在高海拔区域内，高海拔降雨增加导致土壤侵蚀加重，重心北移(表7)。

表6 各行政单元降雨量 mm

表7 不同土地利用类型降雨量 mm

本文利用RUSLE模型计算土壤侵蚀，但还尚有一些不足，本文的研究工作主要借助遥感技术，结果缺乏实际的调查与监测数据，此外，土壤侵蚀量是土壤侵蚀模数与像元面积的乘积，土壤侵蚀时空分布及变化特征分析对于大尺度的土壤侵蚀研究仍然很重要。

4 结 论

(1) 祁连山南坡以微度侵蚀及轻度侵蚀为主，平均土壤侵蚀模数值为2 598.53 t/(km2·a)，整体呈现出由西北向东南递减的趋势，平均侵蚀模数由2000年的1 966.63 t/(km2·a)增加至2005年3 228.51 t/(km2·a)，最后下降至2019年的2 299.06 t/(km2·a)。

(2) 研究区2005年土壤侵蚀最为严重，2000年与2015年平均土壤侵蚀模数为祁连县>门源县>天峻县>刚察县>海晏县；2005年、2010年与2019年平均土壤侵蚀模数为祁连县>门源县>海晏县>刚察县>天峻县。土壤侵蚀量的大小排序为：祁连县>门源县>刚察县>天峻县>海晏县。

(3) 不同土地利用类型平均土壤侵蚀模数大小排序为永久冰川>裸岩、苔原>草地>林地>盐碱地、沙地>裸土>沼泽>耕地；土壤侵蚀量的大小排序为：草地>裸岩>林地>永久冰川>沼泽>盐碱地、沙地>耕地>裸土。

(4) 2000—2019年整体土壤侵蚀重心总共迁移263.09 m，不同县土壤侵蚀重心迁移速率排序为：祁连县(17.59 m/a)>刚察县(1.26 m/a)>门源县(1.06 m/a)>天峻县(0.65 m/a)>海晏县(0.19 m/a)；主要土地利用类型土壤侵蚀重心迁移速率排序为：永久冰川(18.39 m/a)>草地(8.46 m/a)>裸岩(1.19 m/a)>林地(0.18 m/a)。