祁連山南坡土壤侵蝕定量研究與影響因素分析

童 珊, 曹廣超, 閆 欣, 刁二龍, 張 卓

(1.青海師范大學 地理科學學院, 西寧 810008; 2.青海師范大學 青海省自然地理與環(huán)境過程重點實驗室,西寧 810008; 3.青藏高原地表過程與生態(tài)保育教育部重點實驗室, 西寧 810008)

土壤侵蝕會導致土壤肥力、土地利用率和生產力大幅下降[1-2]，甚至制約生態(tài)環(huán)境和經(jīng)濟的可持續(xù)發(fā)展[3-4]。青海省共有水土流失面積16.21萬km2，約占全省土地總面積的23.28%，其中中度侵蝕以上面積占6.46萬km2。為減緩水土流失，需要對其進行定量化的研究，有一些學者已經(jīng)對青藏高原進行的初步研究，例如康琳琦[5]基于USLE模型模擬青藏高原的降雨侵蝕，林慧龍等[6]以137Cs示蹤法為基礎，結合RUSLE模型，以GIS反演為手段，綜合分析三江源區(qū)2001—2012年土壤侵蝕影響因子的特征和土壤侵蝕空間分布規(guī)律，陳豪等[7]基于USLE模型，對祁連山國家公園2005—2019年水力侵蝕進行模擬和計算，分析土壤侵蝕模數(shù)時空分布的動態(tài)變化及影響土壤侵蝕的主導因素，因此RUSLE模型對海拔相差較大區(qū)域具有一定的適應性。

馬爾科夫模型(CA-Markov)可以根據(jù)目前時間的事物狀態(tài)來模擬預測未來某個時間該事物的狀態(tài)，它被廣泛應用于土地利用變化[8]、植被覆蓋動態(tài)變化[9]以及人口等[10]各項地理學研究，但也有對于土壤侵蝕的預測研究，例如武國勝[11]及趙博軒等[12]利用此模型對福建長汀土壤侵蝕動態(tài)預測，趙明松等[13]利用CA-Markov模型對安徽省的土壤侵蝕進行預測。迪氏指數(shù)因子分解模型(LMDI)，是可以從像元尺度定量判斷各因素變化對土壤侵蝕變化的具體貢獻值[14]，該模型能夠對結果進行完全分解，無不能分解的殘余項，與已往的影響因素分析方法例如地理探測器[15]、主成分分析法、相關性分析法[16-17]及通徑分析等[18]方法不同，是研究增長內因及作用強度中廣為應用的一種方法[19-20]。QIAN H E等[21]利用LMDI模型對青藏高原從像素尺度計算植被覆蓋因子(C因子)和降雨侵蝕力因子(R因子)對土壤侵蝕變化的貢獻值。賀倩等[14]利用LMDI模型對三江源土壤侵蝕影響因素進行研究，植被對土壤侵蝕變化整體上具有積極作用，貢獻值范圍集中在-100～100 t/(hm2·a)。

祁連山南坡地形復雜，生態(tài)系統(tǒng)豐富，可以說祁連山南坡是整個祁連山地區(qū)水源涵養(yǎng)的核心區(qū)[22]。本文通過對研究區(qū)不同地形條件下土壤侵蝕的定量研究、未來6 a不同土壤侵蝕強度預測及對影響土壤侵蝕內因的定量化表達，為研究區(qū)土壤侵蝕的治理及決策提供一定的科學支撐。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

祁連山南坡地處青海省境內，山脈包括東北-西南走向的走廊南山、托勒山、托勒南山、大通山和冷龍嶺，地形地貌復雜多樣[23]，研究區(qū)地理位置為東經(jīng)98°08′13″—102°38′16″，北緯37°03′ 17″—39° 05′56″，總面積約為2.4×104km2，區(qū)域內海拔高差大，土壤類型豐富[24]，主要山脈均為西北—東南走向，其間分布山間谷地，該區(qū)是典型的高原大陸性氣候，氣溫日較差大，年平均氣溫為-5.9℃，四季不明顯，年內無絕對無霜期，年降水量約400 mm，主要集中在5—9月，7月、8月最為集中[25]。受地形影響，水熱條件垂直變換明顯，區(qū)內生態(tài)環(huán)境差異較大[26]。

1.2 數(shù)據(jù)來源

氣象數(shù)據(jù):數(shù)據(jù)來源于中國氣象科學數(shù)據(jù)共享服務網(wǎng)(http:∥cdc.cma.gov.cn)下載國家臺站監(jiān)測數(shù)據(jù)，包括祁連山及周邊氣象33個站點的降雨資料，整理后得到站點的月降雨量及年降雨量數(shù)據(jù)，時間為2000年、2005年、2010年、2015年及2019年。

NDVI數(shù)據(jù)：本文使用的植被遙感數(shù)據(jù)MOD13Q1，來源于NASA網(wǎng)站(https:∥ladsweb.modaps.eosdis.nasa.gov/)，空間分辨率為250 m×250 m，時間分辨率為16 d，經(jīng)過拼接、裁剪及最大合成法形成年NDVI數(shù)據(jù)，時間為2000年、2005年、2010年、2015年及2019年。

DEM數(shù)據(jù)：來源于地理空間數(shù)據(jù)云平臺(http:∥www.gscloud.cn)，空間分辨率為90 m×90 m。

不同土地利用類型數(shù)據(jù)：來源于中國科學院地理科學與資源研究所(http:∥www.resdc.cn/)，空間分辨率為30 m×30 m，時間為2000年、2005年、2010年、2015年及2020年。

2 研究方法

2.1 土壤侵蝕方程

本文使用RUSLE土壤流失方程，對研究區(qū)長時期土壤侵蝕進行量化評估，表征降雨和下墊面共同作用下的單位面積潛在土壤流失速率[27]:

A=RKLSCP

(1)

式中:A為土壤侵蝕模數(shù)[t/(hm2·a)];R為降雨侵蝕力因子[(MJ·mm)/(hm2·a)];K為土壤可蝕性因子(t·hm2·h)/(hm2·hm2);LS為坡長坡度因子,無量綱;C為植被覆蓋管理因子,無量綱,值域范圍0～1;P為水土保持措施因子,與土地利用類型相關,無量綱,值域范圍0～1,P=0,表示無侵蝕地區(qū);P=1表示未采取任何水保措施的地區(qū),P值賦值見表1-2。

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

表2 耕地P值

本文運用蔡崇法等人建立的植被覆蓋管理因子(C)與植被覆蓋度的關系計算公式得到各個時期的植被覆蓋管理因子柵格數(shù)據(jù)。修正后的公式如下[28]：

(2)

(3)

式中:C為植被蓋度因子；c為植被覆蓋度。

R反映降雨對土壤侵蝕的影響，是土壤侵蝕預報的重要因子，我們利用祁連山的氣象站點降雨觀測數(shù)據(jù)，整理得到研究區(qū)月平均降雨量和年平均降雨量。通過計算得到整個研究區(qū)R，經(jīng)過隨機森林插值得到研究區(qū)面上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數(shù)據(jù)，采用從國家科技基礎條件平臺—國家地球系統(tǒng)科學數(shù)據(jù)共享服務平臺—黃土高原科學數(shù)據(jù)中心(http:∥loess.geodata.cn)申請得到的Launch LS工具計算研究區(qū)LS因子[30]。

2.2 土壤侵蝕變化率

本文對2000—2019年土壤侵蝕模數(shù)變化率用最小二乘法進行提取，具體公式如下[31]：

(7)

式中:β為目標變量的變化率；X為年份、Y為因變量；n為研究時段總年份。

2.3 馬爾科夫模型

馬爾科夫模型(Markov model)是 1907年由俄國數(shù)學家 Markov在研究布朗運動時發(fā)現(xiàn)的，可以用來模擬事物狀態(tài)的轉移，具有狀態(tài)轉移的無后效性[32]，本文利用馬爾可夫模型對祁連山的土壤侵蝕進行預測，具體公式如下[33]：

(8)

式中：Pij為土壤侵蝕等級類型i到j的轉移概率。

元胞自動機(cellular automata，CA)是一種能夠模擬系統(tǒng)時空演變過程，具有空間計算特征的動力學模型，具體公式如下：

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

(9)

2.4 LMDI模型

迪氏指數(shù)因子分解模型(LMDI模型)，經(jīng)常被應用于驅動因素分析研究，研究的指標因子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中可以看出，研究區(qū)土壤侵蝕模數(shù)呈現(xiàn)出西北向東南遞減的趨勢，2000—2005年侵蝕變化明顯，極強度侵蝕突顯，2010—2015年微度、輕度侵蝕逐漸向中度侵蝕轉變。2000—2019年研究區(qū)整體土壤侵蝕模數(shù)呈增加趨勢，平均值為1.29 t/(km2·a)，變化速率為0.06/a，且土壤侵蝕強烈的區(qū)域土壤侵蝕變化率也相對較大。

圖1 祁連山南坡土壤侵蝕強度分級

3.2 土壤侵蝕強度的時空動態(tài)變化及預測

(1) 土壤侵蝕強度動態(tài)變化。由表3中可以得出，2000—2019年祁連山南坡除了劇烈侵蝕外，其余土壤侵蝕等級向高一級土壤侵蝕強度轉化的比例大于向低一級土壤侵蝕轉化的比例，土壤侵蝕改善不明顯，分階段看，2000—2005年土壤侵蝕強度低一級向高一級轉化幅度最大，其中60.14%的極強度侵蝕轉化為劇烈侵蝕、58%的強度侵蝕轉化為極強度侵蝕。2005—2019年土壤侵蝕等級由高一級轉化低一級較多，土壤侵蝕得到明顯改善。因此總體上看土壤侵蝕未改善可能是由于2000—2005年土壤侵蝕加劇，即使2005年之后政府加強土壤侵蝕防治工作，仍然大于2000年的土壤侵蝕。

(2) 土壤侵蝕預測。本研究用2015—2019年轉移矩陣預測2023年及2027年土壤侵蝕分類等級，2019年預測土地利用類型與實際對比，Kappa系數(shù)為0.92，驗證經(jīng)度較高。從表4中看出，強度以上的侵蝕面積逐漸減少，輕度侵蝕與中度侵蝕面積增加，微度侵蝕面積減少。表明整體上土壤侵蝕強度有減輕趨勢，但同時也應注意降低土壤侵蝕強度低一級向高一級轉變的風險。

3.3 土壤侵蝕在不同地形條件下的變化特征

(1) 不同海拔范圍內土壤侵蝕變化特征。從表5中可以看出，2 700 m以上，海拔越高土壤侵蝕模數(shù)越大，海拔在4 700～5 200 m的土壤侵蝕模數(shù)最大，達到1.05×104t/(km2·a)，土壤侵蝕模數(shù)最低值在海拔2 700～3 200 m，為514.61 t/(km2·a)，而海拔3 700～4 200 m的土壤侵蝕量最大，高達1.26×106t/a，最低值出現(xiàn)在海拔2 200～2 700 m，為1.47×104t/a，海拔在4 700～5 200 m的土壤侵蝕量不大，但單位面積侵蝕量最大。由于祁連山南坡隨著海拔的升高引起水熱條件的變化，土壤植被類型及土地利用類型也發(fā)生著巨大的變化，在海拔4 200 m以上分布著大量的裸巖、永久冰川及積雪等，植被覆蓋較少，易發(fā)生土壤侵蝕，而在3 200～3 700 m人類活動減少，植被覆蓋較好，主要為草地、高寒草甸等，海拔2 700 m以下土壤侵蝕主要受居民點建設、交通等基礎設施的修建等人類活動影響[35]。

表3 祁連山南坡土壤侵蝕強度面積轉移矩陣 %

表4 祁連山南坡2019-2027年土壤侵蝕強度面積占比 %

(2) 不同坡度條件下土壤侵蝕變化特征。土壤侵蝕模數(shù)隨著坡度的增加而增加，坡度>30°的土壤侵蝕模數(shù)最大達到7 256.32 t/(km2·a)，坡度<5°的土壤侵蝕模數(shù)最小，為872.10 t/(km2·a)，而土壤侵蝕量隨著坡度升高而降低，坡度在8°～11°時，土壤侵蝕量達到最大為4.24×105t/a，坡度>30°時土壤侵蝕量最小為1.56萬 t/a。坡度越大，植被活動愈弱[36-37]，遇到強降雨時更容易發(fā)生土壤侵蝕(表6)。

(3) 不同坡向條件下土壤侵蝕變化特征?？傮w上可以看出(表7)，土壤侵蝕大小排序為坡向西>北>南>東>水平方向，土壤侵蝕模數(shù)最大在西北方向，3.37×103t/(km2·a)，土壤侵蝕量最大值出現(xiàn)在西南方向，為3.72×105t/a，無坡向區(qū)土壤侵蝕最弱。由于無坡向植被覆蓋最好[35]，因此土壤侵蝕最小，陳紅等[38]認為各坡向土壤侵蝕差異是由于降雨和季風性刮風等因素造成的，但祁連山南坡地形復雜，降雨量整體上呈現(xiàn)出由東南向西北遞減的趨勢[35]，與侵蝕模數(shù)趨勢正好相反，反而受植被覆蓋影響較大。

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

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

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

3.4 土壤侵蝕影響因素定量研究

(1) 植被覆蓋因子對土壤侵蝕影響定量研究。統(tǒng)計可知植被對土壤侵蝕貢獻值的減少量范圍為0～80 000 t/(km2·a)，面積占研究區(qū)的23.50%，主要分布在研究區(qū)西北部，植被對土壤侵蝕貢獻值的增加量范圍為0～80 000 t/(km2·a)，約占26.97%，研究區(qū)不受植被影響區(qū)域約占49.53%，主要分布在研究區(qū)東南部微度侵蝕區(qū)域(門源縣內)，植被覆蓋較好。2005—2010年植被覆蓋因子對土壤侵蝕貢獻值的減少量范圍為0～80 000 t/(km2·a)，面積增加至26.40%，植被對土壤侵蝕貢獻值的增加量范圍為0～80 000 t/(km2·a)，約占19.92%，研究區(qū)不受植被影響區(qū)域約占53.60%；2010—2015年植被對土壤侵蝕貢獻值的減少量范圍為0～80 000 t/(km2·a)，面積為15.93%，植被對土壤侵蝕貢獻值的增加量范圍為0～80 000 t/(km2·a)，約占36.22%，研究區(qū)不受植被影響區(qū)域約占47.84%；2015—2019年植被對土壤侵蝕貢獻值的減少量范圍為0～80 000 t/(km2·a)，面積為31.53%，植被對土壤侵蝕貢獻值的增加量范圍為0～80 000 t/(km2·a)，約占17.78%，研究區(qū)不受植被影響區(qū)域約占50.68%(圖2)。

研究區(qū)整體上植被對土壤侵蝕變化無影響區(qū)域面積呈增加趨勢，這與研究區(qū)植被覆蓋變化趨勢相符，植被對土壤侵蝕貢獻值的減量效應面積總體上呈擴張趨勢，而增量效應面積整體上呈萎縮趨勢，且除了2010—2015年植被對土壤侵蝕貢獻值的減少量面積小于增加量面積，其余均是減少量面積大于增加量，除說明植被覆蓋對土壤侵蝕具有一定的緩解作用。

(2) 降雨對土壤侵蝕影響定量研究。由圖3可知，降雨對土壤侵蝕模數(shù)的影響主要在0～50 000 t/(km2·a)與0～50 000 t/(km2·a)兩個范圍內。統(tǒng)計可知，2000—2005年降雨對土壤侵蝕貢獻值的減少量范圍為0～50 000 t/(km2·a)，面積占研究區(qū)的8.52%，整個研究區(qū)均有涉及，降雨對土壤侵蝕貢獻值的增加量范圍為0～50 000 t/(km2·a)，約占32.02%，研究區(qū)不受降雨影響區(qū)域約占59.26%，整體分布與植被對土壤侵蝕變化貢獻值分布大致相同。2005—2010年降雨對土壤侵蝕貢獻值的減少量范圍為0～50 000 t/(km2·a)，面積增加至25.29%，降雨對土壤侵蝕貢獻值的增加量范圍為0～50 000 t/(km2·a)，減少至14.45%，研究區(qū)不受降雨影響區(qū)域約占60.24%；2010—2015年降雨對土壤侵蝕貢獻值的減少量范圍為0～50 000 t/(km2·a)，面積為38.26%，降雨對土壤侵蝕貢獻值的增加量范圍為0～50 000 t/(km2·a)，約占2.12%，研究區(qū)不受降雨影響區(qū)域約占59.38%；2015—2019年降雨對土壤侵蝕貢獻值的減少量范圍為0～50 000 t/(km2·a)，面積為15.20%，降雨對土壤侵蝕貢獻值的增加量范圍為0～50 000 t/(km2·a)，約占25.88%，研究區(qū)不受降雨影響區(qū)域約占58.77%。

研究區(qū)整體上降雨對土壤侵蝕模數(shù)變化的貢獻值無影響區(qū)域面積變化較小，降雨對土壤侵蝕貢獻值的減量效應面積總體上呈擴張趨勢，而增量效應面積整體上呈萎縮趨勢。由于2000—2005年及2015—2019年降雨對土壤侵蝕貢獻值的減少量面積小于增加量面積，而2005—2010年及2010—2015年剛好相反，說明降雨量的增加并不一定導致土壤侵蝕的增加[39]。

圖2 植被對土壤侵蝕變化貢獻值分布

圖3 降雨對土壤侵蝕變化貢獻值分布

4 結 論

(1) 土壤侵蝕模數(shù)呈現(xiàn)出西北向東南遞減的趨勢，變化速率為0.06/ a；從土壤侵蝕強度動態(tài)變化來看，2005—2019年土壤侵蝕等級由高一級轉化低一級較多，土壤侵蝕得到明顯改善，且2019—2027年，土壤侵蝕有減輕的趨勢，但也要防止極強度以下的侵蝕低級向高轉變。

(2) 從地形對土壤侵蝕影響來看，土壤侵蝕模數(shù)隨著海拔及坡度的增加而增加，土壤侵蝕量隨著坡度的增加而減小。在海拔4 700～5 200 m及坡度>30°的區(qū)域土壤侵蝕模數(shù)達到最大，分別為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；不同坡向下的土壤侵蝕排序為西>北>南>東>水平方向，土壤侵蝕模數(shù)最大在西北方向，為3 372.58 t/(km2·a)，土壤侵蝕量最大值出現(xiàn)在西南方向，為3.72×105t/a。

(3) 從影響土壤侵蝕變化的內因來看，土壤侵蝕受植被及降雨影響較小區(qū)域主要分布在門源縣；且兩者對土壤侵蝕貢獻值的減量效應面積總體上呈擴張趨勢，而增量效應面積整體上呈萎縮趨勢，說明植被對土壤侵蝕的影響是積極的，降雨量的增加不一定導致土壤侵蝕量的增加。