姚 昆, 周 兵, 何 磊, 劉 斌， 羅 涵， 劉敦龍， 李玉霞

(1.西昌學院 資源與環(huán)境學院, 四川 西昌 615000;2.國家氣候中心 氣候服務室, 北京 100081; 3.成都信息工程大學 軟件工程學院, 成都 610225;4.四川省信息化應用支撐軟件工程技術研究中心, 成都 610225; 5電子科技大學 自動化工程學院, 成都 611731)

土壤侵蝕會造成河道淤塞、耕地質量下降和促進泥石流頻發(fā)等環(huán)境問題[1-3]。區(qū)域的地形地貌、氣候環(huán)境和土壤類型等條件是促進土壤侵蝕形成和改變的先天因素，而后天不合理的人類活動，又會驅動土壤侵蝕強度和分布面積的改變[4-6]。截至2020年，四川省仍有12.1萬km2的地區(qū)存在水土流失現(xiàn)象，其中水力侵蝕就占據(jù)侵蝕總面積的94.21%[7]。因此，加強區(qū)域土壤侵蝕(特別是水力侵蝕)的實時監(jiān)測，掌握其動態(tài)變化規(guī)律并分析背后的驅動因素是實現(xiàn)侵蝕綜合治理重要基礎。

為實現(xiàn)土壤侵蝕狀況的定量評價，學者們將RS和GIS技術應用于該領域，并且提出了通用土壤流失方程(USLE)[8-10]、修正的土壤侵蝕模型(RUSLE)[11-12]、中國土壤流失方程(CSLE)[13-15]和四川省土壤流失方程(SCSLE)等[16]眾多模型。田宇等[5]利用RUSLE模型完成1990—2015年三峽庫區(qū)土壤侵蝕定量評價后，還對其時空變化規(guī)律和驅動因素進行了探討；康琳琦等[12]以RUSLE模型為基礎，分析了青藏高原近30 a來土壤侵蝕的時空變化規(guī)律；姜琳等[17]基于該模型完成了岷江上游流域2000—2010年土壤侵蝕時空格局的動態(tài)變化分析，并探討了侵蝕狀況與高程、坡度和土地利用類型的關系；龔雪梅等[18]基于RUSLE模型完成了1995—2014年岷江上游地區(qū)土壤侵蝕的敏感性評價，并分析了“汶川地震”和土地利用類型改變對其變化的影響。

雖然以上成果均為區(qū)域土壤侵蝕綜合治理提供了重要的參考價值，但是也存在部分不足需要進一步完善。如以往成果在進行驅動力探索時，大都以定性分析或數(shù)理統(tǒng)計(相關性或回歸分析)[5,17-18]為主要方法；前者存在較強的主觀性未能實現(xiàn)驅動因素較準確的定量分析，而后者實施的前提條件是假設在整個時間序列中，土壤侵蝕與驅動因子間存在顯著線性關系，而現(xiàn)實中此關系卻不一定存在[19-20]；同時，以上兩種方法都未將各因子交互產生的協(xié)同作用對侵蝕改變的驅動作用進行分析；而針對以上兩點不足地理探測器卻較好解決[3,20]。綜上，研究以RUSLE模型為基礎，完成岷江上游地區(qū)2000—2018年土壤侵蝕狀況定量評價；借助斜率變化模型完成其時空變化規(guī)律分析；以地理探測器為工具對影響區(qū)域土壤侵蝕變化的驅動力進行探索。期望能為該地區(qū)土壤侵蝕綜合治理提供科學理論參考。

1 研究區(qū)概況

岷江上游地區(qū)地處103°32′—104°15′E，30°45′—33°09′N，包括松潘、黑水和茂縣等5個行政區(qū)，幅員面積約2.4萬km2，海拔表現(xiàn)為東部比西部相對較高，地形起伏差異相對明顯，降水主要集中于夏季(5—10月)，季節(jié)性變化明顯，全域植被生態(tài)系統(tǒng)相對復雜，闊葉林、針葉闊葉林和灌叢等為主要植被類型，土壤則主要包括棕壤、粗骨土和草氈土等，區(qū)域地質構造復雜，地質災害相對頻發(fā)。

2 數(shù)據(jù)來源與研究方法

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

基礎數(shù)據(jù)如下：(1) 1980—2018年岷江上游地區(qū)及周邊氣象站逐日累計降水量數(shù)據(jù)，來源于國家氣候中心；(2) 各年份1∶10萬土地類型矢量成果，由中國科學院資源環(huán)境科學中心(http:∥www.resdc.cn)提供；(3) 90 m分辨率DEM數(shù)據(jù)，來源于地理空間數(shù)據(jù)云平臺；(4) 研究區(qū)30 m土壤可侵蝕性(K)因子，由中科院山地所劉斌濤課題組提供；(5) 1∶100萬中國土壤類型矢量數(shù)據(jù)，來自于中國土壤數(shù)據(jù)庫；(6) 各年份NDVI為MOD13Q1成果，來自USGS。

2.2 研究方法

2.2.1 RUSLE模型 通用土壤流失方式表達式如下：

A=R×LS×K×C×P

(1)

式中:A為侵蝕模數(shù)[t/(hm2·a)];R為降雨侵蝕力[(MJ·mm)/(hm2·h·a)];K為土壤可侵蝕[(t·hm2·h)/(MJ·mm·hm2)];LS為坡度坡長；C,P分別為植被管理與水土保持因子,均無單位;侵蝕模數(shù)乘以100單位轉換為t/(hm2·a)。

其中，R因子以章文波等[21]提出的日降雨量侵蝕力模型計算；K因子為基于EPIC方程[22]計算，結合張科利等[23]提出的修正模型進行校正，最后采用GIS空間插值與平滑優(yōu)化等算法對數(shù)據(jù)進一步優(yōu)化以提高精度的成果(圖1)；L因子則以Liu等[24]提出模型計算；S因子則采用劉斌濤等[25]修正的方程(圖1)；C因子以蔡崇法等[26]基于植被覆蓋度的模型計算(圖2)，而植被覆蓋度則采用像元二分模型[27]；P因子參考已有成果[28-29]結合區(qū)域實際分別賦水田0.15，旱地0.35，草地0.8，林地和未利用為1及其他為0。

2.2.2 斜率變化 斜率可以實現(xiàn)土壤侵蝕與時間的回歸擬合，對其未來發(fā)展趨勢進行預測[30]，本文以最小二乘法實現(xiàn)逐像元侵蝕模數(shù)隨時間的變化分析。表達式如下：

(2)

式中:K為斜率;n為年份;Ai為第i年侵蝕模數(shù)。K<0代表土壤侵蝕狀況有良好的發(fā)展趨勢;相反,K>0則侵蝕狀況有加重的發(fā)展傾向。同時,采用F檢驗,以p為0.05為分界點實現(xiàn)斜率變化的顯著性檢驗。

2.2.3 地理探測器 地理探測器是一種能揭示地理現(xiàn)象空間分異差異和分析自變量與因變量之間相互作用的數(shù)學模型。相比傳統(tǒng)的相關系數(shù)模型，其不僅能實現(xiàn)定量數(shù)據(jù)分析，也能完成定性數(shù)據(jù)的處理，同時又能對各變量之間的交互作用進行分析[20,31]。

圖1 坡度坡長和K因子

圖2 研究區(qū)C因子

因子探測器主要用于探測各個自變量X對因變量Y空間分異的影響力度[20,31]，采用q值進行表征：

(3)

(4)

式中:q為因子解釋力,數(shù)值在[0,1];h為變量X或Y的分類區(qū)間;Nh,N為各層或者全域的分區(qū)數(shù)量；SSW,SST為層內和全區(qū)總方差；q數(shù)值越大表明X對Y的影響力越明顯。

交互探測器可用于分析當各不同自變量X因子存在交互作用時，其交互作用對因變量Y的驅動解釋力是增強或者減弱[20,32-33]。其先分別單獨計算自變量因子X1和X2分別對Y的解釋力q(X1)和q(X2)，其次計算q(X1∩X2)，最后將這3個變量的q值進行比較以此完成交互影響作用的判定(表1)。

本文選取多年平均NDVI(X1)、高程(X2)、多年平均降雨量(X3)、土壤類型(X4)、土地利用類型(X5)和坡度(X6)共計6個因子為自變量，而將多年平均土壤侵蝕模數(shù)為因變量。

表1 交互作用判定

基于王勁峰等[20]提出的數(shù)據(jù)離散化和先驗知識，在參考已有成果[4,19,30]的基礎上利用“自然間斷點分級”將高程、坡度和多年平均降雨量均分為9類，而多年平均NDVI分為13類；土壤類型參照《1∶100萬中華人民共和國土壤圖》將其分為17類，土地利用類型按照耕地、林地和草地等分為6個一級類。同時，以1 km×1 km格網完成該地區(qū)中心樣點的規(guī)則取樣。

3 結果與分析

3.1 土壤侵蝕評價

研究基于USLE模型完成2000—2018年各年份土壤侵蝕模數(shù)的計算，并按《土壤侵蝕分類分級標準(2008)》文件[34]完成各年份侵蝕程度的等級劃分(圖3)。

結合圖3分析，該地區(qū)土壤侵蝕強度空間分布變化規(guī)律與姜琳等[17]基本吻合，間接證明研究成果基本可靠。以2018年岷江上游地區(qū)侵蝕強度分級結果為例，該地區(qū)侵蝕以微度和輕度為主，它們占據(jù)了全域總面積的77.59%，這些土壤侵蝕現(xiàn)象并不明顯的地區(qū)在整個研究區(qū)的分布范圍最廣，分布在全域的大部分地區(qū)。結合土地利用類型和地形地貌數(shù)據(jù)分析可知，微度和輕度區(qū)地貌大多以中—高海拔的山地為主，而這些地區(qū)土地類型則又主要為有林地、灌木林地和高度蓋度草地等，這些區(qū)域植被生態(tài)環(huán)境相對良好且覆蓋度相對較高，良好的植被生態(tài)環(huán)境對調節(jié)氣候變化和生態(tài)環(huán)境好轉又起到了有利的促進作用，微度侵蝕全域分布范圍廣可能就是得益于此。

圖3 2000-2018年土壤侵蝕強度分級

相反，中度及以上程度在全域分布范圍僅占據(jù)全域總面積的22.41%，它們絕大部分分布于高海拔且地形起伏差異明顯的高山和極高山地區(qū)，這些地區(qū)呈現(xiàn)出植被覆蓋度相對較低，地形陡峭，特殊的氣候易引發(fā)巖石風化和土壤剝離搬運等現(xiàn)象，因此土壤侵蝕程度相對明顯。此次，它們還有極小部分分布于岷江、黑水河和雜谷腦河等江河水系兩側，這些地區(qū)雖然地勢相對平緩，但它們是社會經濟發(fā)展的核心地帶，也是受人類活動影響相對明顯的區(qū)域，全域主要的農牧生產均主要集中于此，促成了這些地區(qū)土壤侵蝕現(xiàn)象相對明顯現(xiàn)象的產生。

3.2 趨勢變化分析

結合岷江上游地區(qū)土壤侵蝕斜率變化(圖4)計算結果分析可知：2000—2018年內，土壤侵蝕呈現(xiàn)好轉(K<0)的地區(qū)全域面積占比為75.60%；而整個岷江上游地區(qū)僅有25.40%的侵蝕呈現(xiàn)出嚴重化發(fā)展趨勢。結合顯著性檢驗結果又可知：整個時段內，有21.14%的地區(qū)土壤侵蝕呈現(xiàn)明顯好轉的發(fā)展狀態(tài)，它們主要分布于松潘縣東部地區(qū)、岷江、黑水河以及雜谷腦河等江河水系兩側地形相對平緩的區(qū)域；侵蝕狀況未發(fā)生明顯改變(好轉或嚴重化)的地區(qū)全域面積占比最大為76.99%，在全域的大部分范圍均有分布；僅有1.87%的地區(qū)呈現(xiàn)出土壤侵蝕明顯嚴重化的發(fā)展傾向，結合地形地貌數(shù)據(jù)分析可以發(fā)現(xiàn)，它們主要分布在海拔4 000 m以上的高山地區(qū)。

圖4 斜率及顯著檢驗

3.3 侵蝕驅動因子探測

3.3.1 侵蝕驅動因子探測 結合(表2)對土壤侵蝕驅動因子的關系進行分析可知：各個因子對區(qū)域土壤侵蝕空間分布格局變化的驅動作用存在不同差異；其影響力大小關系為多年平均NDVI>多年平均降水>高程>土壤類型>土地利用類型>坡度；6個因子中多年平均NDVI的q值最大為0.331 4，說明其對土壤侵蝕空間分布格局形成和改變的影響力最大，坡度的q值最小僅為0.017 2，和其他5個因子相比坡度對侵蝕強度空間分布格局變化的驅動作用最不明顯。

表2 驅動因子探測結果

結合q值變化與實際分析可知：岷江上游地區(qū)森林和中高覆蓋度草地分布范圍相對較廣，全域95%以上為林草區(qū)(其中，中高植被覆蓋區(qū)又占據(jù)林草區(qū)總面積的80%以上)，全域植被覆蓋狀況整體相對良好，植被對土壤侵蝕的形成與變化起到相對明顯的抑制作用；區(qū)域海拔差異明顯，地形起伏相對較大，降雨又主要集中于下半年(5—10月)和海拔2 500 m以下的河谷地帶，夏季暴雨相對頻發(fā)，降雨強度也相對較大，是促進岷江上游地區(qū)土壤侵蝕面積擴張和程度加深的主要因素；此外，多年平均植被覆蓋度、高程和多年平均降雨量這3個因子的解釋力均明顯遠大于土地利用類型；綜上，客觀程度可以判定，自然因子是影響區(qū)域土壤侵蝕空間分布格局變化的主要驅動因素。

3.3.2 交互探測 交互探測主要用于判定，若兩個驅動因子產生交互關系時，其交互作用對土壤侵蝕變化的影響作用為增強、減弱或者相互獨立。通過交互探測器完成岷江上游地區(qū)土壤侵蝕與各驅動因子的交互探測，結果表明各驅動因子間不存在相互獨立作用，因子交互作用表現(xiàn)為非線性增強和雙因子增強兩種。

結合(表3)對各因子之間的交互影響作用進行分析：多年平均NDVI分別與高程、多年平均降雨量或土地利用類型；高程分別與多年平均降雨量、土壤類型、土地利用類型、多年平均降雨量或坡度。以上幾類驅動因子產生交互關系時，它們的交互協(xié)同作用對土壤侵蝕分布格局變化產生雙因子增強影響；其他幾類因子產生交互關系時，交互協(xié)同作用對土壤侵蝕分布格局變化產生的影響作用為非線性增強；當多年平均NDVI分別與多年平均降雨量、高程或土壤類型產生交互關系時，它們的協(xié)同作用對土壤侵蝕分布格局變化的影響力均達到40%以上，客觀程度表征了3類組合產生的交互作用對土壤侵蝕分布格局變化影響的貢獻作用最大；相反，土壤類型和坡度產生交互關系時，其協(xié)同作用對土壤侵蝕分布格局變化的影響力最小，僅產生5%左右的驅動作用；多年平均NDVI與其他驅動因子產生交互關系時，它們產生的協(xié)同作用對土壤侵蝕分布格局變化的影響力都在30%以上，均比其他因子間產生交互的影響力大，這也意味著植被覆蓋度差異較大的地區(qū)土壤侵蝕程度通常也呈現(xiàn)出較明顯的差異。

表3 因子交互探測結果

4 討論與結論

4.1 討 論

岷江上游地區(qū)既是四川省西部生態(tài)高原建設的重要功能區(qū)，也是土壤侵蝕問題相對典型的地區(qū)之一。因此，較全面掌握其空間分布特征、時間變化規(guī)律和變化驅動因素對科學實現(xiàn)該地區(qū)土壤侵蝕綜合治理有著重要意義。研究以RUSLE為基礎，完成該地區(qū)2000—2018年土壤侵蝕狀況的定量評價與侵蝕動態(tài)變化規(guī)律的探索，結合地理探測器對影響岷江上游地區(qū)土壤侵蝕空間分布格局和變化的驅動因素進行分析。對研究結果分析可知，本文采用的方法可行，預期目標也基本實現(xiàn)。

相比以往成果，本研究引入用于分析長時間序列數(shù)據(jù)空間變化規(guī)律分析的線性回歸模型，其能較好避免由于某時期數(shù)據(jù)突變而造成整體趨勢變化規(guī)律出現(xiàn)異?，F(xiàn)象的發(fā)生，這對較全面掌握區(qū)域土壤侵蝕真實的變化規(guī)律更具價值；采用地理探測器對影響區(qū)域土壤侵蝕空間分布格局變化因素進行探索，其不僅能較有效探測地理現(xiàn)象空間差異性，揭示其背后驅動力，還能分析因子交互協(xié)同作用對侵蝕強度空間分布格局變化的影響作用；其較有效彌補了以往成果僅從單因子對驅動力進行分析的不足，進一步促進了對驅動力的較全面掌握。

植被覆蓋度(NDVI)、海拔和降水是影響岷江上游地區(qū)土壤侵蝕分布及程度改變的主要驅動因素。特別是植被覆蓋度，單因子探測結果顯示其影響力幾乎是海拔和降水的兩倍。雙因子交互探測結果也顯示，植被覆蓋度(NDVI)與其他因子產生交互協(xié)同作用的驅動力也比其他因子兩兩交互產生的影響力大。以上兩點均表明，植被覆蓋度(NDVI)是驅動該地區(qū)土壤侵蝕形成與變化的最主要因素，其影響力遠遠大于其他因子的驅動作用。在整個研究時段內，岷江上游地區(qū)土壤侵蝕狀況整體得到有效遏制，對祝聰?shù)萚30]研究成果分析可知，整個階段內岷江上游地區(qū)植被覆蓋度整體處于良好發(fā)展的趨勢，植被生態(tài)環(huán)境的改善對區(qū)域土壤侵蝕的改變起到了有效的抑制作用。因此，今后對該地區(qū)土壤侵蝕進行綜合治理時，加強區(qū)域植被的保護與努力提高植被覆蓋度仍然是一項重要的工作。

由于影響土壤侵蝕變化的因素相對較多，同時受數(shù)據(jù)收集的限制，本研究仍然存在一定的不足尚需完善。如在計算土壤侵蝕模數(shù)時，本文以250 m的柵格尺寸為尺度，雖然利用遙感技術進行時空變化規(guī)律分析時大都更側重宏觀角度，該地區(qū)范圍也相對較大，該柵格尺度未能對動態(tài)變化規(guī)律分析產生較明顯干擾，但若能采用更高分辨率數(shù)據(jù)進行分析，那研究結果的可靠性會進一步提高；同時，研究在進行驅動力因子選取時，主要選擇了部分典型因子，今后的研究中若加入地層巖性、人口和GDP等更多因子進行探討，那么驅動力探討的結果也更具價值。

4.2 結 論

(1) 岷江上游地區(qū)土壤侵蝕整體狀況相對良好，中度以上侵蝕區(qū)全域的分布范圍相對較少，并且它們主要集中于地形起伏較大且植被覆蓋度相對較低的高海拔山地區(qū)；近20 a內，整個地區(qū)土壤侵蝕呈現(xiàn)出良好的發(fā)展趨勢。

(2) 因子探測結果表明，在眾多驅動因素中，植被覆蓋度(NDVI)、降雨量和海拔三者占據(jù)了促進土壤侵蝕形成與變化因素的前3位，其中植被覆蓋度(NDVI)的驅動作用最明顯；交互探測結果顯示，當影響土壤侵蝕變化的各因子產生交互作用時，它們的協(xié)同驅動作用均比單因子的影響力明顯。

(3) 對該地區(qū)土壤侵蝕時空變化規(guī)律和驅動力分析可知，加強該地區(qū)植被的恢復建設，努力提高區(qū)域植被覆蓋度是有效治理地區(qū)土壤侵蝕的重要方法，也是今后治理工作的重點努力方向。