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

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

土壤侵蚀会导致土壤肥力、土地利用率和生产力大幅下降[1-2]，甚至制约生态环境和经济的可持续发展[3-4]。青海省共有水土流失面积16.21万km2，约占全省土地总面积的23.28%，其中中度侵蚀以上面积占6.46万km2。为减缓水土流失，需要对其进行定量化的研究，有一些学者已经对青藏高原进行的初步研究，例如康琳琦[5]基于USLE模型模拟青藏高原的降雨侵蚀，林慧龙等[6]以137Cs示踪法为基础，结合RUSLE模型，以GIS反演为手段，综合分析三江源区2001—2012年土壤侵蚀影响因子的特征和土壤侵蚀空间分布规律，陈豪等[7]基于USLE模型，对祁连山国家公园2005—2019年水力侵蚀进行模拟和计算，分析土壤侵蚀模数时空分布的动态变化及影响土壤侵蚀的主导因素，因此RUSLE模型对海拔相差较大区域具有一定的适应性。

马尔科夫模型(CA-Markov)可以根据目前时间的事物状态来模拟预测未来某个时间该事物的状态，它被广泛应用于土地利用变化[8]、植被覆盖动态变化[9]以及人口等[10]各项地理学研究，但也有对于土壤侵蚀的预测研究，例如武国胜[11]及赵博轩等[12]利用此模型对福建长汀土壤侵蚀动态预测，赵明松等[13]利用CA-Markov模型对安徽省的土壤侵蚀进行预测。迪氏指数因子分解模型(LMDI)，是可以从像元尺度定量判断各因素变化对土壤侵蚀变化的具体贡献值[14]，该模型能够对结果进行完全分解，无不能分解的残余项，与已往的影响因素分析方法例如地理探测器[15]、主成分分析法、相关性分析法[16-17]及通径分析等[18]方法不同，是研究增长内因及作用强度中广为应用的一种方法[19-20]。QIAN H E等[21]利用LMDI模型对青藏高原从像素尺度计算植被覆盖因子(C因子)和降雨侵蚀力因子(R因子)对土壤侵蚀变化的贡献值。贺倩等[14]利用LMDI模型对三江源土壤侵蚀影响因素进行研究，植被对土壤侵蚀变化整体上具有积极作用，贡献值范围集中在-100～100 t/(hm2·a)。

祁连山南坡地形复杂，生态系统丰富，可以说祁连山南坡是整个祁连山地区水源涵养的核心区[22]。本文通过对研究区不同地形条件下土壤侵蚀的定量研究、未来6 a不同土壤侵蚀强度预测及对影响土壤侵蚀内因的定量化表达，为研究区土壤侵蚀的治理及决策提供一定的科学支撑。

1 研究区概况与数据来源

1.1 研究区概况

祁连山南坡地处青海省境内，山脉包括东北-西南走向的走廊南山、托勒山、托勒南山、大通山和冷龙岭，地形地貌复杂多样[23]，研究区地理位置为东经98°08′13″—102°38′16″，北纬37°03′ 17″—39° 05′56″，总面积约为2.4×104km2，区域内海拔高差大，土壤类型丰富[24]，主要山脉均为西北—东南走向，其间分布山间谷地，该区是典型的高原大陆性气候，气温日较差大，年平均气温为-5.9℃，四季不明显，年内无绝对无霜期，年降水量约400 mm，主要集中在5—9月，7月、8月最为集中[25]。受地形影响，水热条件垂直变换明显，区内生态环境差异较大[26]。

1.2 数据来源

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

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

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

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

2 研究方法

2.1 土壤侵蚀方程

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

A=RKLSCP

(1)

式中:A为土壤侵蚀模数[t/(hm2·a)];R为降雨侵蚀力因子[(MJ·mm)/(hm2·a)];K为土壤可蚀性因子(t·hm2·h)/(hm2·hm2);LS为坡长坡度因子,无量纲;C为植被覆盖管理因子,无量纲,值域范围0～1;P为水土保持措施因子,与土地利用类型相关,无量纲,值域范围0～1,P=0,表示无侵蚀地区;P=1表示未采取任何水保措施的地区,P值赋值见表1-2。

表1 不同土地利用类型P值

表2 耕地P值

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

(2)

(3)

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

R反映降雨对土壤侵蚀的影响，是土壤侵蚀预报的重要因子，我们利用祁连山的气象站点降雨观测数据，整理得到研究区月平均降雨量和年平均降雨量。通过计算得到整个研究区R，经过随机森林插值得到研究区面上R，插值精度均在0.8以上。

(4)

式中:pi为月均降雨量(mm)；p为年平均降雨量(mm)。

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

SN=1-SAN/100

(5)

(6)

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

坡长坡度因子LS可由DEM计算获取，基于渭河流域的90 m DEM数据，采用从国家科技基础条件平台—国家地球系统科学数据共享服务平台—黄土高原科学数据中心(http:∥loess.geodata.cn)申请得到的Launch LS工具计算研究区LS因子[30]。

2.2 土壤侵蚀变化率

本文对2000—2019年土壤侵蚀模数变化率用最小二乘法进行提取，具体公式如下[31]：

(7)

式中:β为目标变量的变化率；X为年份、Y为因变量；n为研究时段总年份。

2.3 马尔科夫模型

马尔科夫模型(Markov model)是 1907年由俄国数学家 Markov在研究布朗运动时发现的，可以用来模拟事物状态的转移，具有状态转移的无后效性[32]，本文利用马尔可夫模型对祁连山的土壤侵蚀进行预测，具体公式如下[33]：

(8)

式中：Pij为土壤侵蚀等级类型i到j的转移概率。

元胞自动机(cellular automata，CA)是一种能够模拟系统时空演变过程，具有空间计算特征的动力学模型，具体公式如下：

S(t+1)=f〔S(t),N〕

(9)

2.4 LMDI模型

迪氏指数因子分解模型(LMDI模型)，经常被应用于驱动因素分析研究，研究的指标因子R，K，LS，C，P，本文主要研究植被与降雨对土壤侵蚀变化的影响。具体公式如下：

(10)

ΔVtot=VT-V0=ΔVx1+ΔVx2+…+ΔVxn

(11)

式中:ΔVtot为目标量变化;VT为计划期T时的目标量;V0为基期0时的目标量;i=1,2,…,n。式(11)反映了从基期0到计划期T时间段内目标量的变化[14，34]。

(12)

其中，

(13)

3 结果与分析

3.1 土壤侵蚀的时空变化

从图1中可以看出，研究区土壤侵蚀模数呈现出西北向东南递减的趋势，2000—2005年侵蚀变化明显，极强度侵蚀突显，2010—2015年微度、轻度侵蚀逐渐向中度侵蚀转变。2000—2019年研究区整体土壤侵蚀模数呈增加趋势，平均值为1.29 t/(km2·a)，变化速率为0.06/a，且土壤侵蚀强烈的区域土壤侵蚀变化率也相对较大。

图1 祁连山南坡土壤侵蚀强度分级

3.2 土壤侵蚀强度的时空动态变化及预测

(1) 土壤侵蚀强度动态变化。由表3中可以得出，2000—2019年祁连山南坡除了剧烈侵蚀外，其余土壤侵蚀等级向高一级土壤侵蚀强度转化的比例大于向低一级土壤侵蚀转化的比例，土壤侵蚀改善不明显，分阶段看，2000—2005年土壤侵蚀强度低一级向高一级转化幅度最大，其中60.14%的极强度侵蚀转化为剧烈侵蚀、58%的强度侵蚀转化为极强度侵蚀。2005—2019年土壤侵蚀等级由高一级转化低一级较多，土壤侵蚀得到明显改善。因此总体上看土壤侵蚀未改善可能是由于2000—2005年土壤侵蚀加剧，即使2005年之后政府加强土壤侵蚀防治工作，仍然大于2000年的土壤侵蚀。

(2) 土壤侵蚀预测。本研究用2015—2019年转移矩阵预测2023年及2027年土壤侵蚀分类等级，2019年预测土地利用类型与实际对比，Kappa系数为0.92，验证经度较高。从表4中看出，强度以上的侵蚀面积逐渐减少，轻度侵蚀与中度侵蚀面积增加，微度侵蚀面积减少。表明整体上土壤侵蚀强度有减轻趋势，但同时也应注意降低土壤侵蚀强度低一级向高一级转变的风险。

3.3 土壤侵蚀在不同地形条件下的变化特征

(1) 不同海拔范围内土壤侵蚀变化特征。从表5中可以看出，2 700 m以上，海拔越高土壤侵蚀模数越大，海拔在4 700～5 200 m的土壤侵蚀模数最大，达到1.05×104t/(km2·a)，土壤侵蚀模数最低值在海拔2 700～3 200 m，为514.61 t/(km2·a)，而海拔3 700～4 200 m的土壤侵蚀量最大，高达1.26×106t/a，最低值出现在海拔2 200～2 700 m，为1.47×104t/a，海拔在4 700～5 200 m的土壤侵蚀量不大，但单位面积侵蚀量最大。由于祁连山南坡随着海拔的升高引起水热条件的变化，土壤植被类型及土地利用类型也发生着巨大的变化，在海拔4 200 m以上分布着大量的裸岩、永久冰川及积雪等，植被覆盖较少，易发生土壤侵蚀，而在3 200～3 700 m人类活动减少，植被覆盖较好，主要为草地、高寒草甸等，海拔2 700 m以下土壤侵蚀主要受居民点建设、交通等基础设施的修建等人类活动影响[35]。

表3 祁连山南坡土壤侵蚀强度面积转移矩阵 %

表4 祁连山南坡2019-2027年土壤侵蚀强度面积占比 %

(2) 不同坡度条件下土壤侵蚀变化特征。土壤侵蚀模数随着坡度的增加而增加，坡度>30°的土壤侵蚀模数最大达到7 256.32 t/(km2·a)，坡度<5°的土壤侵蚀模数最小，为872.10 t/(km2·a)，而土壤侵蚀量随着坡度升高而降低，坡度在8°～11°时，土壤侵蚀量达到最大为4.24×105t/a，坡度>30°时土壤侵蚀量最小为1.56万 t/a。坡度越大，植被活动愈弱[36-37]，遇到强降雨时更容易发生土壤侵蚀(表6)。

(3) 不同坡向条件下土壤侵蚀变化特征。总体上可以看出(表7)，土壤侵蚀大小排序为坡向西>北>南>东>水平方向，土壤侵蚀模数最大在西北方向，3.37×103t/(km2·a)，土壤侵蚀量最大值出现在西南方向，为3.72×105t/a，无坡向区土壤侵蚀最弱。由于无坡向植被覆盖最好[35]，因此土壤侵蚀最小，陈红等[38]认为各坡向土壤侵蚀差异是由于降雨和季风性刮风等因素造成的，但祁连山南坡地形复杂，降雨量整体上呈现出由东南向西北递减的趋势[35]，与侵蚀模数趋势正好相反，反而受植被覆盖影响较大。

表5 不同海拔梯度下的土壤侵蚀量

表6 不同坡度下的土壤侵蚀量

表7 不同坡向下的土壤侵蚀量

3.4 土壤侵蚀影响因素定量研究

(1) 植被覆盖因子对土壤侵蚀影响定量研究。统计可知植被对土壤侵蚀贡献值的减少量范围为0～80 000 t/(km2·a)，面积占研究区的23.50%，主要分布在研究区西北部，植被对土壤侵蚀贡献值的增加量范围为0～80 000 t/(km2·a)，约占26.97%，研究区不受植被影响区域约占49.53%，主要分布在研究区东南部微度侵蚀区域(门源县内)，植被覆盖较好。2005—2010年植被覆盖因子对土壤侵蚀贡献值的减少量范围为0～80 000 t/(km2·a)，面积增加至26.40%，植被对土壤侵蚀贡献值的增加量范围为0～80 000 t/(km2·a)，约占19.92%，研究区不受植被影响区域约占53.60%；2010—2015年植被对土壤侵蚀贡献值的减少量范围为0～80 000 t/(km2·a)，面积为15.93%，植被对土壤侵蚀贡献值的增加量范围为0～80 000 t/(km2·a)，约占36.22%，研究区不受植被影响区域约占47.84%；2015—2019年植被对土壤侵蚀贡献值的减少量范围为0～80 000 t/(km2·a)，面积为31.53%，植被对土壤侵蚀贡献值的增加量范围为0～80 000 t/(km2·a)，约占17.78%，研究区不受植被影响区域约占50.68%(图2)。

研究区整体上植被对土壤侵蚀变化无影响区域面积呈增加趋势，这与研究区植被覆盖变化趋势相符，植被对土壤侵蚀贡献值的减量效应面积总体上呈扩张趋势，而增量效应面积整体上呈萎缩趋势，且除了2010—2015年植被对土壤侵蚀贡献值的减少量面积小于增加量面积，其余均是减少量面积大于增加量，除说明植被覆盖对土壤侵蚀具有一定的缓解作用。

(2) 降雨对土壤侵蚀影响定量研究。由图3可知，降雨对土壤侵蚀模数的影响主要在0～50 000 t/(km2·a)与0～50 000 t/(km2·a)两个范围内。统计可知，2000—2005年降雨对土壤侵蚀贡献值的减少量范围为0～50 000 t/(km2·a)，面积占研究区的8.52%，整个研究区均有涉及，降雨对土壤侵蚀贡献值的增加量范围为0～50 000 t/(km2·a)，约占32.02%，研究区不受降雨影响区域约占59.26%，整体分布与植被对土壤侵蚀变化贡献值分布大致相同。2005—2010年降雨对土壤侵蚀贡献值的减少量范围为0～50 000 t/(km2·a)，面积增加至25.29%，降雨对土壤侵蚀贡献值的增加量范围为0～50 000 t/(km2·a)，减少至14.45%，研究区不受降雨影响区域约占60.24%；2010—2015年降雨对土壤侵蚀贡献值的减少量范围为0～50 000 t/(km2·a)，面积为38.26%，降雨对土壤侵蚀贡献值的增加量范围为0～50 000 t/(km2·a)，约占2.12%，研究区不受降雨影响区域约占59.38%；2015—2019年降雨对土壤侵蚀贡献值的减少量范围为0～50 000 t/(km2·a)，面积为15.20%，降雨对土壤侵蚀贡献值的增加量范围为0～50 000 t/(km2·a)，约占25.88%，研究区不受降雨影响区域约占58.77%。

研究区整体上降雨对土壤侵蚀模数变化的贡献值无影响区域面积变化较小，降雨对土壤侵蚀贡献值的减量效应面积总体上呈扩张趋势，而增量效应面积整体上呈萎缩趋势。由于2000—2005年及2015—2019年降雨对土壤侵蚀贡献值的减少量面积小于增加量面积，而2005—2010年及2010—2015年刚好相反，说明降雨量的增加并不一定导致土壤侵蚀的增加[39]。

图2 植被对土壤侵蚀变化贡献值分布

图3 降雨对土壤侵蚀变化贡献值分布

4 结 论

(1) 土壤侵蚀模数呈现出西北向东南递减的趋势，变化速率为0.06/ a；从土壤侵蚀强度动态变化来看，2005—2019年土壤侵蚀等级由高一级转化低一级较多，土壤侵蚀得到明显改善，且2019—2027年，土壤侵蚀有减轻的趋势，但也要防止极强度以下的侵蚀低级向高转变。

(2) 从地形对土壤侵蚀影响来看，土壤侵蚀模数随着海拔及坡度的增加而增加，土壤侵蚀量随着坡度的增加而减小。在海拔4 700～5 200 m及坡度>30°的区域土壤侵蚀模数达到最大，分别为10 460.72，7 256.32 t/(km2·a)。在海拔3 700～4 200 m的土壤侵蚀量最大，高达1.26×106t/a， 2 200～2 700 m土壤侵蚀最小，为1.47×104t/a，在坡度8°～11°时，土壤侵蚀量达到最大为4.24×105t/a，坡度>30°土壤侵蚀量最小为1.56万t/a；不同坡向下的土壤侵蚀排序为西>北>南>东>水平方向，土壤侵蚀模数最大在西北方向，为3 372.58 t/(km2·a)，土壤侵蚀量最大值出现在西南方向，为3.72×105t/a。

(3) 从影响土壤侵蚀变化的内因来看，土壤侵蚀受植被及降雨影响较小区域主要分布在门源县；且两者对土壤侵蚀贡献值的减量效应面积总体上呈扩张趋势，而增量效应面积整体上呈萎缩趋势，说明植被对土壤侵蚀的影响是积极的，降雨量的增加不一定导致土壤侵蚀量的增加。