1901—2016年黃土高原土壤侵蝕格局演變及其驅動機制

穆興民，李朋飛，劉斌濤，趙廣舉，高 鵬，孫文義

(1.西北農林科技大學 黃土高原土壤侵蝕與旱地農業國家重點實驗室，陜西 楊凌 712100;2.中國科學院 水利部 水土保持研究所，陜西 楊凌 712100;3.西安科技大學 測繪科學與技術學院，陜西 西安 710054;4.中國科學院 水利部 成都山地災害與環境研究所，四川 成都 610041）

土壤侵蝕對陸地生態系統安全構成極大威脅[1-2］,是我國頭號生態環境問題。黃土高原地形復雜、溝壑縱橫、土質松散、降雨集中,是我國乃至世界水土流失危害最嚴重的地區之一[3-5］。氣候變化和人類活動從正反兩個方面影響土壤侵蝕的發生和發展。全新世中期以前,黃土高原土壤侵蝕主要受氣候變化影響,而近1000 a來,人類活動是黃土高原土壤侵蝕不斷加劇之主導因素[6］,人口增長導致農作物種植面積擴大與土地利用強度提高,進而加劇水土流失(面積增加、強度提高)[7-8］；近100 a來,在氣候變化和人類活動雙重作用下,黃土高原土地利用及地表覆蓋顯著變化,土壤侵蝕模數劇烈波動[9-11］。探明黃土高原土壤侵蝕模數時空變化特征,對黃土高原土壤侵蝕演變及其驅動機制的認識意義重大。水土保持作為新興學科,歷史觀測數據不足制約了對歷史時期黃土高原土壤侵蝕時空動態變化及其驅動機制認識的深化。

土壤侵蝕模數獲取方法有控制模擬試驗[12-14］、野外定位觀測[15-16］、遙感影像分析[17-18］、土壤侵蝕模型模擬[19-21］等。控制模擬試驗、野外定位觀測費時費力且成本較高,適用于小尺度研究；遙感影像分析僅能獲取20世紀80年代以后的土壤侵蝕變化特征[22］；土壤侵蝕模型模擬雖因目前對土壤侵蝕發生發展的機理尚不非常清楚而存在諸多問題,但仍不失為土壤侵蝕時空演變評估的有效方法。已有較多學者采用上述方法對黃土高原進行了不同時間尺度的土壤侵蝕調查、監測與模擬,但研究時段多集中于20世紀80年代以后[23-24］,鮮有針對百年尺度土壤侵蝕模數的模擬評估。

本文基于修正通用土壤流失方程RUSLE(Revised Universal Soil Loss Eguation)構建了1901—2016年逐年黃土高原土壤侵蝕模數柵格數據集,并采用人類活動干擾較小的22個黃河支流水文站的輸沙量實測數據對其進行驗證,基于雙累積曲線法分析了不同時期土壤侵蝕模數變化的主要驅動因素,以期解決黃土高原土壤侵蝕模數歷史數據缺乏問題,為黃土高原土壤侵蝕演變過程研究及水土保持生態建設規劃實施等提供方法與數據支撐。

1 資料與方法

1.1 黃土高原概況

黃土高原地處我國西北部,總面積約64萬km2,地勢西北高、東南低,大致分為高塬溝壑區、農灌區、沙地沙漠區、丘陵溝壑區、土石山區、河谷平原區等6個土壤侵蝕地貌類型區(見圖1)。黃土高原屬干旱和半干旱氣候區[25］,多年平均降水量150～700 mm[8］,降水量遠低于潛在蒸散發能力[26］,土壤含水率低,植被類型從南部的森林過渡到北部的荒漠草原[27］。由于黃土高原地區降雨主要發生在6—9月且具有降雨歷時短、強度大的特點[28］,加之地形復雜、溝壑縱橫、土質疏松,極易發生嚴重的水土流失[4,10-11］。20世紀以來,黃土高原人類活動日益加劇,如農業生產、戰爭、森林砍伐、植被恢復、水土保持措施與水利工程建設、城市化等,對區域土壤侵蝕的發生發展產生重要的影響[5］,加之氣候變化導致的極端天氣事件,使黃土高原土壤侵蝕表現出顯著的時空分異特征。

圖1 黃土高原土壤侵蝕地貌類型分區

1.2 RUSLE

RUSLE根據降雨、土壤、地形、植被覆蓋、水土保持措施等因子計算土壤侵蝕模數,其概念清晰,所需數據較少且易于獲取,計算方便[29］,是應用廣泛的土壤侵蝕模型之一。RUSLE形式為

式中:A為土壤侵蝕模數,t/(hm2·a)；R為降雨侵蝕力,MJ·mm/(hm2·h·a)；K為土壤可蝕性因子,t·hm2·h/(hm2·MJ·mm)；LS為地形因子(其中L為坡長因子、S為坡度因子)；C為植被覆蓋因子；P為水土保持措施因子。

(1)降雨侵蝕力的確定。已有大量研究提出了我國不同區域的降雨侵蝕力計算方法,本研究采用以月平均雨量估算侵蝕力的簡易算法[30］計算黃土高原地區降雨侵蝕力:

其中

式中:F為與年降雨量的季節分布有關的指標,mm；p、pi分別為年降雨量、第i個月份降雨量(降雨量數據源自高分辨率地表氣候格網數據集,網址為https://chelsa-climate.org/chelsacruts/),mm。

(2)土壤可蝕性因子的確定。基于EPIC模型估算土壤可蝕性因子,并采用Zhang等[31］提出的方法對其進行修正:

式中:KEPIC為基于EPIC模型估算的土壤可蝕性因子；San、Sil、Cla分別為土壤的沙粒、粉粒、黏粒含量,%；OC為土壤有機碳含量,%；Sn為常數,Sn＝1-San/100。土壤類型和屬性數據從中國土壤科學數據庫(http://vdb3.soil.csdb.cn/extend/jsp/introduction)獲取。

(3)坡長和坡度因子的確定。坡長和坡度因子基于數字高程模型(DEM)進行計算,DEM由分辨率為30 m的全球數字高程模型ASTER GDEM(https://nordpil.com/blog/astergdem/)經重采樣得到。計算坡長和坡度因子的傳統方法適用于坡度≤18%的區域[32］,而黃土高原大部分區域地表坡度＞18%[33］,Liu等[34］發現坡度為9%～55%的坡面土壤侵蝕強度與坡面坡度的正弦呈線性關系,江忠善等[35］根據我國不同地區土壤侵蝕研究成果提出了冪函數形式的坡度坡長因子計算方法,本研究綜合考慮上述研究成果,采用如下公式計算坡度和坡長因子:

其中

式中:θ為坡度,%；r為坡長,m；m為指數。

(4)植被覆蓋因子的確定。植被覆蓋因子根據土地利用類型和植被特征來確定。對源于全新世土地利用數據集HYDE3.2.1(https://archaeology.datastations.nl/dataset.xhtml？persistentId＝doi:10.17026/dans-25g-gez3)的土地利用數據進行重采樣,獲得分辨率為1 km的柵格數據[36］,將土地利用類型分為森林、草地、耕地、居民地、水體,把水體、居民地的植被覆蓋因子分別設定為0.001、0.1,森林、草地、耕地的植被覆蓋因子根據2011年全國土壤侵蝕調查數據[37］與1 km分辨率NDVI之間的關系進行計算:

根據遙感影像可以獲取1982年之后的NDVI數據,但1982年以前因遙感影像資料缺乏而難以確定NDVI。筆者分析NDVI與月降水量的關系并構建了二者的線性回歸方程(其相關系數為0.77、顯著性水平為0.001),據此方程,結合由CHELSA cruts數據集得出的月降水量構建了1901—2016年NDVI數據集,并依此估算歷史時期的C值。

(5)水土保持措施因子的確定。水土保持措施因子根據HYDE3.2.1和土地利用數據集CCI＿LC(https://www.esa-landcover-cci.org/)來確定,其中CCI＿LC提供了灌溉農田和雨養農田的空間分布[38］。將HYDE3.2.1數據重采樣至300 m,與CCI＿LC數據進行疊加,提取各項水土保持措施和耕地(把耕地分為梯田和坡耕地,其中梯田包括灌溉農田和緩坡雨養農田)。參照已有研究成果[24］,把梯田、居民地的P值分別設定為0.2、0.01,把森林、草地、水體、坡耕地的P值均設定為1。

1.3 RUSLE模擬計算精度驗證

區域土壤侵蝕模數是難以測度的。本研究基于RUSLE計算土壤侵蝕模數的空間分辨率為1 km,無法根據小區實測數據進行驗證,故采用中小河流實測輸沙量對RUSLE計算結果進行驗證。黃土高原地區在無壩庫攔蓄情況下中小河流泥沙輸移比接近于1[8］,在1970年以前水土保持措施尤其淤地壩數量很少,中小河流輸沙量接近坡面侵蝕產沙量,因此本研究選取有實測資料的22個黃河支流水文站1919—1969年實測輸沙量(換算為輸沙模數)對RUSLE模擬計算的對應支流土壤侵蝕模數進行驗證,驗證結果用二者的相關系數、Nash系數及顯著性來反映。

1.4 土壤侵蝕強度分級

根據水利部發布的《土壤侵蝕分類分級標準》(SL 190—2007),依據土壤侵蝕模數將土壤侵蝕強度劃分為6級:微度,侵蝕模數＜1000 t/(km2·a)；輕度,侵蝕模數為1000～2500 t/(km2·a)；中度,侵蝕模數為2500～5000 t/(km2·a)；強烈,侵蝕模數為5000～8000 t/(km2·a)；極強烈,侵蝕模數為8000～15000 t/(km2·a)；劇烈,侵蝕模數＞15000 t/(km2·a)。

1.5 驅動因素貢獻評估

雙累積曲線是有效評估氣候和人類活動對河川徑流量及輸沙量變化影響的有效方法之一[39］。本研究采用土壤侵蝕模數與降水量雙累積曲線確定土壤侵蝕模數變化的各驅動因素的貢獻(如圖2所示),在不受人類活動影響的情況下,某一河流的多年累積土壤侵蝕模數和累積降水量二者之比為常數,雙累積曲線為線段OT,TT′為其延長線；當流域土壤侵蝕受人類活動影響劇烈時,雙累積曲線變為TB1或TB2(TB1為人類活動加劇侵蝕,TB2為人類活動減緩侵蝕),其與直線OT′的偏差為人類活動對土壤侵蝕模數變化的貢獻量。

圖2 土壤侵蝕模數與降水量雙累積曲線示意

2 結 果

2.1 RUSLE模擬計算精度驗證結果

分別計算22個支流水文站實測輸沙量與RUSLE模擬計算的土壤侵蝕模數的相關系數和Nash系數,以驗證RUSLE模擬計算結果的精度。分析表1可知:1957—1969年選取的16個支流水文站實測輸沙量與計算的土壤侵蝕模數相關系數為0.48～0.80,其中有10個支流的相關系數大于0.60、有9個支流的相關性達到極顯著水平；除白家川和舊縣外,其余支流模擬結果的Nash系數接近1,其中有9個支流的Nash系數大于等于0.4。此外,1919—1953年享堂、靜樂、蘭村、趙城、頭、南河川等水文站共有實測輸沙資料22 a,與對應年份RUSLE計算的土壤侵蝕模數相關系數為0.75、相關性達到極顯著水平。驗證結果說明,RUSLE模擬計算結果總體可信,對部分支流模擬計算的可信度較高。

表1 1957—1969年RUSLE模擬計算精度驗證結果

2.2 黃土高原土壤侵蝕模數時空變化情況

黃土高原地區1901—2016年多年平均土壤侵蝕模數為5056.86 t/(km2·a),不同區域土壤侵蝕模數存在明顯差異(見圖3),其中:黃土高原西北部沙地沙漠區與農灌區土壤侵蝕模數大都小于1000 t/(km2·a)；黃土高原中部的丘陵溝壑區、高塬溝壑區平均侵蝕模數分別達8570.06、5781.82 t/(km2·a),黃土高原土壤侵蝕模數大于8000 t/(km2·a)的極強烈和劇烈侵蝕主要集中在這2個類型區,因此本研究把這2個類型區作為分析的重點區域(圖3中紅線所圍區域)。

圖3 黃土高原年平均土壤侵蝕模數空間變化情況

1901—2016年黃土高原及重點區域土壤侵蝕模數年際變化情況見圖4。1901—1909年,黃土高原平均土壤侵蝕模數為4160.55 t/(km2·a),維持在相對較低水平；1910—1919年,黃土高原平均土壤侵蝕模數與1901—1909年相比有一定的上升,達到4928.70 t/(km2·a)；1920—1929年平均土壤侵蝕模數上升至5663.00 t/(km2·a)；1930—1979年,黃土高原平均土壤侵蝕模數為6195.69 t/(km2·a),強烈以上等級的土壤侵蝕面積急劇增長,有13 a平均侵蝕模數超過8000 t/(km2·a)；1980—2009年,黃土高原的平均土壤侵蝕模數顯著下降,平均侵蝕模數為3401.30 t/(km2·a),比侵蝕嚴重的1930—1979年平均值下降了45.1%；2010—2016年,黃土高原平均土壤侵蝕模數有所反彈。

圖4 黃土高原及其重點區域土壤侵蝕模數年際變化過程

表3 1901—2016年丘陵溝壑區各級土壤侵蝕面積占比變化情況 %

1901—2016年黃土高原及2個重點區域各級土壤侵蝕面積占比變化情況見表2～表4。黃土高原中度及以下土壤侵蝕面積占比較大(為54.56%～81.31%),尤其1930年以前與1980年以后中度及以下土壤侵蝕面積占比較大、強烈及以上土壤侵蝕面積占比較小,而1930—1979年強烈及以上土壤侵蝕面積占比相對較大(為43.80%～45.44%)；丘陵溝壑區強烈及以上土壤侵蝕面積占比1979年以前為56.21%～72.16%、1980年以后下降到36.24%～55.96%,高塬溝壑區強烈及以上土壤侵蝕面積占比1920—1979年為50.24%～53.40%、1980年后為24.38%～36.05%；2010—2016年黃土高原及2個重點區域強烈及以上土壤侵蝕面積占比均有所反彈。

表2 1901—2016年黃土高原各級土壤侵蝕面積占比變化情況 %

表4 1901—2016年高塬溝壑區各級土壤侵蝕面積占比變化情況 %

根據黃土高原氣候變化及人類活動情況及其對土壤侵蝕的影響,可將研究時段(1901—2016年)大致劃分為自然狀態時期(1901—1929年)、災害破壞時期(1930—1949年)、初步治理時期(1950—1969年)、全面治理時期(1970—1998年)及高質量生態恢復時期(1999—2016年)等5個時期(見圖5),各時期黃土高原年均土壤侵蝕模數分別為4943.49、6415.55、6022.78、4575.70、3431.79 t/(km2·a)。從自然狀態時期到災害破壞時期,黃土高原土壤侵蝕模數呈上升態勢,在災害破壞時期以后隨著人類活動(各類水土保持措施的實施)土壤侵蝕模數呈降低趨勢。

圖5 黃土高原各時期土壤侵蝕模數變化趨勢

黃土高原土壤侵蝕模數與汛期(6—9月)降雨量關系密切,二者波動變化趨勢基本一致,尤其1920—1980年丘陵溝壑區二者的波動一致性更為明顯,見圖6。極端暴雨往往造成極端侵蝕,把黃土高原平均土壤侵蝕模數超過10000 t/(km2·a)的年份稱為極端侵蝕年份,則研究時段(1901—2016年)有4 a為極端侵蝕年份(見圖5),分別為1917年、1925年、1949年、2013年,這4個年份土壤侵蝕模數分別為11717、15202、10853、10968 t/(km2·a)。

圖6 土壤侵蝕模數與汛期降雨量逐年變化過程

2.3 不同時期土壤侵蝕模數變化情況及驅動機制

黃土高原土壤侵蝕以水蝕為主,其主要影響因素可分為氣候(主要是降水)和人類活動影響下的土地利用、地表覆蓋及微地形改變兩類。采用土壤侵蝕模數和降水量雙累積曲線(見圖7)分析黃土高原及重點區域土壤侵蝕模數變化的階段性特征及驅動機制:黃土高原及2個重點區域土壤侵蝕模數均發生了2次突變,其中黃土高原及高塬溝壑區土壤侵蝕模數2次突變年份均為1924年和1981年、丘陵溝壑區土壤侵蝕模數2次突變年份為1931年和1981年。按突變年份可把研究時段(1901—2016年)分為3個階段,把第一次突變前、后分別稱為基準期、變化期(分為第一變化期和第二變化期),黃土高原和高塬溝壑區的基準期為1901—1923年、丘陵溝壑區的基準期為1901—1930年,按照前述方法基于雙累積曲線計算降水和人類活動造成變化期土壤侵蝕模數的變化量及其占比,結果見表5。

圖7 土壤侵蝕模數與降水量雙累積曲線

分析表5可知:第一變化期,黃土高原及2個重點區域土壤侵蝕模數較基準期均有所上升,黃土高原地區土壤侵蝕模數階段平均值由4367.6 t/(km2·a)上升至6317.9 t/(km2·a),其中人類活動是土壤侵蝕模數上升的主要驅動因素(造成的變化量為1652.7 t/(km2·a),貢獻占比為84.7%)；第二變化期,黃土高原及2個重點區域土壤侵蝕模數較基準期均顯著降低,黃土高原土壤侵蝕模數階段平均值降至3475.9 t/(km2·a),相對于基準期、第一變化期的降幅分別為20.4%、45.0%,人類活動是土壤侵蝕模數降低的驅動因素,降水對土壤侵蝕模數降低的貢獻占比為負值,盡管該階段降水量較之前有所增大、理論上應增大土壤侵蝕模數,但大規模的水土保持措施和水利工程建設等人類活動仍然使土壤侵蝕模數大幅度降低。

表5 降水與人類活動對侵蝕模數的影響

3 討 論

3.1 近百年來黃土高原土壤侵蝕模數演變特征

我國區域土壤侵蝕評估研究多聚焦于20世紀80年代以后,而河流水沙變化研究始于20世紀20年代。本研究通過區域土壤侵蝕模擬計算,首次構建了黃土高原長時段逐年土壤侵蝕模數數據集,并探討了百年時間尺度的黃土高原土壤侵蝕模數變化情況,以期為黃土高原歷史時期侵蝕模數研究提供方法與數據支撐。黃土高原20世紀初土壤侵蝕模數相對較低、30—70年代土壤侵蝕模數明顯上升、80年代后呈下降趨勢,21世紀初處于歷史最低水平,2010年以后土壤侵蝕模數有所反彈。與黃土高原土壤侵蝕模數的年代際變化相對應,黃河輸沙量變化趨勢為在20世紀30年代之前相對較小、30—60年代持續增大、80年代后急劇減小[40-41］。

3.2 人類活動對土壤侵蝕模數變化的影響

人類活動是黃土高原百年尺度土壤侵蝕模數變化的主導因素,這與以往研究結論基本一致[42-44］,上述第一變化期不合理的人類活動導致黃土高原土壤侵蝕模數上升,第二變化期大規模水土保持措施的實施有效降低了土壤侵蝕模數。

20世紀30—40年代我國處于戰爭時期,人口遷移導致黃土高原中部人口大量增加[41］、森林和草地被大面積砍伐和墾殖,加之當時人們對生態環境保護認識不足,因而使土壤侵蝕迅速加劇。

新中國成立后的20世紀50—70年代,國家實施“以糧為綱”的方針,黃土高原地區人口快速增長,經濟發展仍以種植農業為主,為了增加耕地面積和解決農村薪柴問題,森林被進一步砍伐、荒草地被開墾種植[45］,導致土壤侵蝕進一步加劇。盡管這一時期已開始加強水土保持工作,但其效果并不理想[8］。

20世紀80年代以后國家進一步重視水土保持工作,黃土高原地區推廣旱地農業增產技術和以小流域為單元的山水林田路綜合治理,在保障農業生產的同時,使土壤侵蝕模數有所降低、入黃泥沙量顯著減少。生態系統是人與自然生命共同體之基,是人類經濟社會存在和發展之基。隨著改革開放的推進,全國社會經濟快速發展,黃土高原乃至西北地區水土流失嚴重、生態環境脆弱問題與經濟社會發展不適應的矛盾日益突出,為了再造一個山川秀美的西北地區,從20世紀90年代末大規模退耕還林(草)工程的實施,以及黃土高原化石能源開發利用根本上改變了農村能源結構,黃土高原植被覆蓋率明顯提高,積極貫徹落實生態文明建設基本國策、科學實施各類水土保持項目使水土保持工作逐步向高質量發展轉變,同時城市化吸引農村人口大規模向城市遷移,因而有效緩解了黃土高原地區的人地矛盾[46］,使黃土高原地區土壤侵蝕環境得到系統性改變,土壤侵蝕模數顯著降低(尤其丘陵溝壑區和高塬溝壑區降低更為顯著)[47］、入黃泥沙量顯著減少。深思細想,各項水土保持措施是黃土高原土壤侵蝕模數降低的驅動因素,而退耕還林(草)等正確的水土保持政策是控制黃土高原土壤侵蝕的根本驅動力。

3.3 極端降雨及其侵蝕災害的積極應對

黃土高原時段平均土壤侵蝕模數在2000—2009年處于近百年最低水平的情況下,于2010—2016年出現小幅反彈,其與2013年7月降雨異常偏多造成極端侵蝕有關。隨著社會經濟發展,在全球極端暴雨頻發的背景下,未來黃土高原水土保持工作須重點考慮強監管、重維護、提標準,以提高防御極端降雨引發極端侵蝕的能力。長期以來,對水土保持措施存在重建輕管的問題,使得各項水土保持措施因缺乏維護維修資金而老化失修,制約了其效應的持續有效發揮,如大量老舊的小型淤地壩幾近淤滿,淤地壩和梯田的管理維護不到位,存在垮塌致災風險,這是極端暴雨導致土壤侵蝕模數及河流輸沙量突增的主要原因之一,因此應加強“后水土保持”監管工作。經過20余a來退耕還林(草)政策引導和持續治理,黃土高原已基本呈現溝坡林灌成蔭、梁峁芳草鋪地的景觀特點,水土流失得到有效遏制,但如果對剛剛恢復的仍顯脆弱的林草“資源”進行不合理開發利用,那么將很快使林草植被遭到破壞、入黃泥沙量“回歸”高位,因此必須鞏固水土流失防治和退耕還林(草)成果。此外,植被恢復中重林輕草、盲目還“林”導致的人工林退化以及土壤干燥化問題比較突出,亟須在鞏固和提升退耕還林(草)成果基礎上進一步優化各項水土保持措施配置,提升生態系統保持水土的功能；現有水土保持措施投資標準及建設標準規范多產生于20世紀80年代,據此建設的水土保持措施難以抵御全球氣候變化背景下的極端暴雨及極端侵蝕事件,應盡快修訂和完善水土保持措施設計、施工、驗收等標準規范體系,以滿足新時代水土保持高質量發展的需求。

4 結論

基于RUSLE建立1901—2016年黃土高原土壤侵蝕模數逐年序列1 km分辨率柵格數據集,評估百年尺度的黃土高原土壤侵蝕模數時空變化情況,并探討了不同時期黃土高原土壤侵蝕模數變化的驅動機制,結論如下:

(1)用22個黃河支流水文站在人類活動影響較小時期的實測輸沙量對RUSLE模擬計算的土壤侵蝕模數進行驗證,結果表明二者具有較好的一致性,采用RUSLE模擬計算黃土高原土壤侵蝕模數的精度較高。

(2)百余年來黃土高原土壤侵蝕模數在1924年左右和1981年左右發生了2次突變,20世紀30—70年代土壤侵蝕面積持續增大、侵蝕模數持續上升,20世紀80年代—21世紀初黃土高原土壤侵蝕模數持續降低,至21世紀初降為百余年來的最低值,2010—2016年土壤侵蝕模數出現小幅反彈。

(3)人類活動是黃土高原土壤侵蝕模數變化的主要驅動因素。人口增長導致森林過度砍伐、耕地擴張及過度放牧是20世紀30—70年代土壤侵蝕加劇的主要原因,20世紀80年代及以后各項水土保持措施建設及退耕還林(草)政策的實施使黃土高原土壤侵蝕模數大幅降低。

(4)丘陵溝壑區和高塬溝壑區是黃土高原水土流失嚴重區域和治理的重點區域。

(5)極端降雨事件是2010—2016年黃土高原土壤侵蝕模數發生反彈的主要原因,鑒于全球氣候變化背景下未來極端暴雨事件或將更加頻繁,亟須創新黃土高原土壤侵蝕防治策略,提高防御極端暴雨造成嚴重水土流失災害的能力。