1901-2016年黃土高原降雨侵蝕力時空變化

2022-07-03 03:09:44陳劍南劉益麟李朋飛胡晉飛高健健黨恬敏

水土保持研究 2022年4期

陳劍南, 劉益麟, 李朋飛, 胡晉飛, 高健健, 黨恬敏

(1.西安科技大學測繪科學與技術學院, 西安 710054; 2.黃河水利委員會綏德水土保持科學試驗站, 陜西榆林 719000; 3.黃河流域水土保持生態環境監測中心, 西安 710021)

黃土高原是全球土壤侵蝕最嚴重的地區之一[1]。水力侵蝕是黃土高原最重要的侵蝕類型，而降水及其產生的徑流是水力侵蝕的主要驅動力。降雨侵蝕力用于表征降雨引起土壤侵蝕的潛在能力[2]，目前已廣泛用于氣候變化監測、土壤侵蝕模擬等研究中[3-4]。因此，厘清降雨侵蝕力時空特征是開展黃土高原土壤侵蝕研究與防治的基礎。

20世紀50年代以來，中外學者針對降雨侵蝕力時空變化開展了大量研究。例如，Panagos等[5]估算了希臘1965—1994年月平均降雨侵蝕力；Liu等[6]分析了1980—2017年全球降雨侵蝕力的時空變化；劉斌濤等[7]估算了1960—2009年中國590個氣象站的降雨侵蝕力，并通過插值得到中國降雨侵蝕力的時空分布；殷水清等[8]繪制了1961—2016年全國1 km降雨侵蝕力等值線圖，分析了中國降雨侵蝕力的時空分布及重現期。黃土高原是我國降雨侵蝕力研究的重點區域。例如，穆興民等[9]繪制了1956—2002年陜北黃土高原降雨侵蝕力的等值線并分析了其時空變化特征；孫從建等[10]研究了黃土高原塬面保護區內降雨侵蝕力的時空變化趨勢及影響因素；李靜等[11]研究了1971—2000年黃土高原不同地貌分區的降雨侵蝕力的時空特征，并得到降雨侵蝕力存在2.7 a的波動周期；KEO等[12]利用Kringing內插法得到黃土高原空間連續分布的多年平均降雨侵蝕力分布并討論其在1960—2011年的變化趨勢。然而，已有研究多關注20世紀50年代以后黃土高原降雨侵蝕力變化，對于更長時段(如橫跨20世紀百年尺度)的降雨侵蝕力研究依然缺乏，限制了對黃土高原土壤侵蝕長時段變化及其與氣候變化關系的理解。

綜上所述，本文基于高分辨率地表氣候格網數據集的逐月降水數據，評估1901—2016年黃土高原降雨侵蝕力時空變化特征。首先采用逐月實測降水量驗證CHELSAcrust數據集的精度，并基于該數據集估算黃土高原1901—2016年降雨侵蝕力；其次，利用Mann-Kendall非參數檢驗法、經驗模態分解法(Empirical Mode Decomposition， EMD)、累積距平、逐像元趨勢檢驗等方法分析黃土高原1901—2016年降雨侵蝕力的時空變化特征，為黃土高原土壤侵蝕侵蝕防治與預報提供參考。

1 材料與方法

1.1 研究區概況

黃土高原地處中國中部偏北，東起太行山，西至烏鞘嶺，南接秦嶺，北抵長城(100.8°—114.6°E，33.7°—41.3°N)，可分為沙地沙漠區、農灌區、丘陵溝壑區、高原溝壑區、河谷平原區、土石山區(圖1)[13]。黃土高原溝壑縱橫，地形復雜，總面積約64萬km2。該區域屬于半濕潤、干旱和半干旱區，年均氣溫4～12℃，年均降水量400～600 mm，降水以短歷時暴雨為主，且集中于7—9月[14]。區域內土壤結構松散，抗蝕能力較差，植被覆蓋率較低，加之汛期暴雨多發，造成了嚴重的土壤侵蝕。20世紀70年代以來，黃土高原實施了大規模水土保持措施，尤其自1999年以來，實施了退耕還林(草)工程，有效緩解了區域內的土壤侵蝕。然而，近年來極端降水事件多發[15-16]，為已有水土保持措施帶來了挑戰。

圖1 研究區位置

1.2 數據來源

1901—2016年降水數據來自高分辨率地表氣候格網數據集(Climatologies at high resolution for the earth′s land surface areas，CHELSAcrust)。該數據集基于英國東英吉利大學氣候研究中心生產的CRU TS 4.01數據集，采用增量更改算法生產而成[17]，提供20世紀以來的氣候數據資料，空間分辨率為1 km，數據集獲取網址為https:∥chelsa-climate.org/chelsacruts/。此外，由國家氣象信息中心編制的《中國地面氣候資料日值數據集》，獲取了黃土高原內國家氣象站1980—2016年的實測降水量數據，用于驗證CHELSAcrust數據集的精度，獲取網址為：http:∥data.cma.cn/。

1.3 研究方法

1.3.1 降雨侵蝕力估算方法降雨侵蝕力計算方法較多，初期研究多基于次降水強度及降雨動能估算降雨侵蝕力[18-19]。然而，次降水信息較難獲取，限制了該類方法的應用。2003年章文波等[20-21]建立了基于不同類型降水資料(如日降水量、月降水量、年降水量)的降雨侵蝕力估算方法，為利用氣象站整編降水資料評估降雨侵蝕力提供了方法基礎，已在黃土高原上得到廣泛應用。本文采用逐月降雨侵蝕力模型計算區域的降雨侵蝕力，計算方法如下:

(1)

(2)

式中:Pi,j為第i年j月的降水量(mm)；N為年數;R為多年平均降雨侵蝕力[MJ·mm/(hm2·h)]；α4，β4為模型參數。模型參數選用逐月雨量的模型參數,即α4=0.1833,β4=1.9957。

1.3.2 降水數據集精度驗證方法為驗證CHELSA crust數據集精度，選取黃土高原內均勻分布的13個氣象站點(圖1)，利用1980—2016年逐月降水量觀測數據，以Nash效率系數(Nash-Sutcliffe efficiency coefficient，NSE)、決定系數(Coefficient of determination，R2)的結果來評定逐月降水數據集精度。R2取值范圍0～1，Nash取值范圍為-∞～1，其值越接近1表明模擬精度越高，降雨數據集越可靠，一般認為，R2>0.6，Nash>0.5時，結果可信[22]。

1.3.3 降雨時空分布特征分析方法

(1) 經驗模態分解法(EMD)。利用經驗模態分解法分析1901—2016年降雨侵蝕力年際周期性。經驗模態分解法(Empirical mode decomposition，EMD)是Huang等[23]提出的一種自適應信號時頻處理方法。該方法依據數據系列的瞬時特征，將原始時間序列信號按照頻率特征，從高到低分解為若干個本征模函數(Intrinsic model function，IMF)，所分解出來的IMF分量包含原信號不同時間尺度的局部特征信號，以反映原始時間序列的周期與振幅[24]。

(2) 非參數Mann-Kendall突變檢驗。Mann-Kendall突變檢驗法是檢驗氣候、水文等要素時間序列發生突變的非參數檢驗方法。該方法通過繪制UF，UB曲線，當UF>0，則表明序列呈現上升趨勢，反之同理。給定顯著性水平α，若|UF|>Uα，則表明序列存在顯著的趨勢變化。若UF和UB相交，且交點在臨界線Uα之內，則交點對應的時刻可能為突變開始的時間[25]。本文使用該法評估1901—2016年降雨侵蝕力的突變特征，置信水平設置為95%(Uα=±1.96)。

(3) 累積距平法。使用累積距平法探明1901—2016年黃土高原和各個分區降雨侵蝕力的年際變化的趨勢和階段。其計算方法如下：

(3)

(4) 逐像元趨勢性檢驗。通過逐像元一元線性回歸來判斷1901—2016年黃土高原地區降雨侵蝕力時空變化特征及顯著性。使用ArcGIS 10.3.1柵格計算器，計算一元線性回歸中各項參數，得到逐柵格降雨侵蝕力線性回歸方程，并進行顯著性F檢驗，基于給定的顯著性α=0.05，通過查表判斷線性變化的顯著性，通過ArcGIS 10.3.1重分類工具實現降雨侵蝕力空間變化分布和顯著性分布。

2 結果與分析

2.1 降水數據集精度驗證結果

基于1980—2016年各氣象站點的降水量實測值，驗證CHELSAcrust數據集逐月降水量的精度(表1)。結果表明，各氣象站點Nash系數最小為0.65，最大為0.91，均值為0.79，R2最大為0.89，最小為0.65，均值為0.82，說明該數據集模擬精度較高，能夠滿足長時間序列分析需求。

表1 CHELSAcrust數據集逐月降水量精度驗證結果

2.2 降雨侵蝕力空間分布特征分析

2.2.1 年均降雨侵蝕力空間分布 1901—2016年間黃土高原年均降雨侵蝕力差距較大，從東南向西北遞減(圖2)，東部大部分地區降雨侵蝕力超過3 000 MJ·mm/(hm2·h)，而西北地區降雨侵蝕力較小，均在1 000 MJ·mm/(hm2·h)以下。降雨侵蝕力大于4 000 MJ·mm/(hm2·h)的地區主要集中在土石山區和丘陵溝壑區東部。各分區的降雨侵蝕力差異明顯，土石山區降雨侵蝕力較大，平均降雨侵蝕力為2 414 MJ·mm/(hm2·h)。農灌區的降雨侵蝕力均值較小，為699.30 MJ·mm/(hm2·h)。不同分區降雨侵蝕力由大到小依次為：土石山區>河谷平原區>丘陵溝壑區>高原溝壑區>沙地沙漠區>農灌區。

在20世紀20年代之前，黃土高原平均降雨侵蝕力為1 233.35 MJ·mm/(hm2·h)，其中土石山區的降雨侵蝕力最大，為5 253.75 MJ·mm/(hm2·h)。降雨侵蝕力<2 000 MJ·mm/(hm2·h)的面積占整個黃土高原面積的80.92%～88.44%，而降雨侵蝕力>4 000 MJ·mm/(hm2·h)的面積僅占總面積的0.34%～1.43%；表明該時期黃土高原的降雨侵蝕力總體維持在較低水平；20世紀30年代開始，黃土高原的降雨侵蝕力激增，尤其是1930—1939年，土石山區以及丘陵溝壑區東部的極強降雨侵蝕力[>5 000 MJ·mm/(hm2·h)]顯著增加，面積占比達到3.10%(圖2)。1930—1979年，降雨侵蝕力在1 000～3 000 MJ·mm/(hm2·h)的面積占比達到60%以上，降雨侵蝕力>4 000 MJ·mm/(hm2·h)的面積占比為1.99%～6.53%，比1901—1929年增加4.5倍以上，各分區在這一時期的降雨侵蝕力達到峰值，且增大趨勢向西北地區的高原溝壑區和沙地沙漠區延伸。20世紀80年代—21世紀初，降雨侵蝕力逐漸下降，降雨侵蝕力<2 000 MJ·mm/(hm2·h)的面積占到整個黃土高原面積的79.04%～88.59%，而降雨侵蝕力>4 000 MJ·mm/(hm2·h)的面積僅占0.57%～0.16%，且降雨侵蝕力較大的地區逐漸退至黃土高原東部的土石山區局部地區。2010—2016年，降雨侵蝕力>4 000 MJ·mm/(hm2·h)面積占比增至3.06%，表明黃土高原降雨侵蝕力高值區域再次增加。

圖2 1901-2016年黃土高原降雨侵蝕力時空分布

如圖3所示，1901—2016年，黃土高原降雨侵蝕力增加的區域主要集中在黃土高原中部及西部部分區域，減少的區域主要分布在黃土高原西部和東北部。降雨侵蝕力的顯著性變化在空間分布上呈現出邊緣地區多為非顯著變化而中部地區多為顯著變化的規律。其中，丘陵溝壑區中部、沙地沙漠區中部、高原溝壑區西部和東部等區域的降雨侵蝕力顯著增長；沙地沙漠區東部和高原溝壑區部分地區的降水侵蝕力顯著減少。

圖3 1901-2016年黃土高原降雨侵蝕力變化趨勢

2.3 降雨侵蝕力時間演變特征

2.3.1 降雨侵蝕力年際變化 1901—2016年，黃土高原平均降雨侵蝕力呈非顯著下降趨勢(圖4)。年均降雨侵蝕力為1 462.17 MJ·mm/(hm2·h)，其中，1925年降雨侵蝕力為研究時間序列最大值，為3 497.3 MJ·mm/(hm2·h)，較平均值高出約139.2%。1991年降雨侵蝕力出現最小值，為545.4 MJ·mm/(hm2·h)，相較于平均值減少了約62.7%。對各分區而言，土石山區降雨侵蝕力明顯高于其他五大分區，沙地沙漠區及農灌區的降雨侵蝕力相比之下較小。沙地沙漠區和河谷平原區的降雨侵蝕力非顯著上升，其他4個分區降雨侵蝕力非顯著下降。各分區降雨侵蝕力隨時間的變化相似，在1901—1916年期間，各區降雨侵蝕力維持在較低水平[<4 000 MJ·mm/(hm2·h)]，1917年降雨侵蝕力激增，丘陵溝壑區達到峰值5 242.12 MJ·mm/(hm2·h)；1917—1960年，降雨侵蝕力明顯增大且各分區較大降雨侵蝕力集中出現，1939年土石山區達到峰值9 001.02 MJ·mm/(hm2·h)，1940年沙地沙漠區達到峰值2 037.20 MJ·mm/(hm2·h)，1949年河谷平原區和高原溝壑區均達到峰值，分別為5 072.44，2 497.41 MJ·mm/(hm2·h)，1959年農灌區達到峰值2 030.51 MJ·mm/(hm2·h)。1961—2000年間，降雨侵蝕力波動下降，除河谷平原區其他各區均達到谷值，沙地沙漠區、農灌區、丘陵溝壑區谷值在1965年，高原溝壑區出現在1972年，土石山區出現在1997年。2010年后降雨侵蝕力出現反彈，2013年出現峰值。

黃土高原和各分區的年降雨侵蝕力累積距平結果變化趨勢類似(圖5)，可將降雨侵蝕力變化劃分為1901—1930年、1931—1980年、1980—2016年3個階段，1930年前累積距平值為負，說明降雨侵蝕力低于平均值。自1930年起，降雨侵蝕力與多年均值差值為正，累積距平值呈現持續上升趨勢，到1957年，累積距平值達到0，之后持續上升，直到1979年達到峰值。1980年開始，累積距平值逐步下降，表明降雨侵蝕力低于多年平均水平。各分區的降雨侵蝕力距平與黃土高原降雨侵蝕力累積距平的變化趨勢基本一致。值得注意的是，黃土高原及各分區的降雨侵蝕力累積距平變化于2010年出現較大反復，說明2010年之后的降雨侵蝕力有增強跡象。

利用Mann-Kendall法對黃土高原及6個分區進行突變診斷(圖6)，整個黃土高原降雨侵蝕力UF和UB序列均存在多個交叉點，6個分區的降雨侵蝕力與黃土高原的UF和UB曲線波動基本一致，交點大多分布在20世紀20年代以前和90年代以后。結合年降雨侵蝕力變化趨勢(圖4)及年降雨侵蝕力累積距平(圖5)，說明黃土高原的降雨侵蝕力不存在明顯突變。黃土高原及6個分區降雨侵蝕力的正序列在1930年之前均表現為先升后降，UF值在0刻度線上下浮動，無持續性趨勢；1930年之后，黃土高原、高原溝壑區及河谷平原區的正序列曲線表現為先升后降，但UF值始終>0，表明降雨侵蝕力在增加，在1980年之后增加幅度減小，且均在60—80年代左右突破0.05顯著性水平線，說明在1960—1980年降雨侵蝕力增加較顯著。土石山區與以上3個區域相似，但在2007年之后UF值在0值上下波動，說明2007年之后該區降雨侵蝕力反復變化但不顯著。農灌區、沙地沙漠區、丘陵溝壑區在2000年之后UF<0且未超過置信區間，表明該時段內降雨侵蝕力非顯著減小。

2.3.2 降雨侵蝕力周期性分析利用EMD法對黃土高原地區1901—2016年的116 a降雨侵蝕力進行分析，得到了5個本征模分量(IMF)及最低頻分量(RES)(圖7)。表2計算了各個IMF分量、方差以及方差貢獻率。可以看出5個本征模分量(IMF)累計方差貢獻達到96.03%，而其中IMF1貢獻率最大，為56.53%，其次為IMF2分量，為32.77%。IMF5分量的貢獻率最小，為1.48%。由此得出，黃土高原地區降雨侵蝕力變化存在周期性規律，2.75 a變化周期最顯著，其次為5.33 a，再其次為10.29，24.5，44.5 a，這3個周期的方差貢獻率差異不明顯，均為弱周期。降雨侵蝕力的周期變化受制于該區域氣候要素的波動周期。太陽黑子和大氣環流的異常活動可以影響水熱平衡，從而改變降水時序和空間分布。研究表明，太陽黑子存在11 a的活動周期[26]，這與本研究得出的本征模分量IMF3(10.29 a)的周期基本吻合。大氣環流異常因子厄爾尼諾現象也存在2～7 a的變化周期[27]，本研究所得出的本征模分量IMF1(2.75 a)，IMF2(5.33 a)的周期與其基本一致。

注：△為降雨侵蝕力最大值；▽為降雨侵蝕力量小值。

圖5 1901-2016年黃土高原及各分區年均

3 結論

(1) CHELSAcrust數據集降水量數據精度較高，可用于黃土高原長時段研究需求。

(2) 黃土高原降雨侵蝕力總體為東部大于西部。黃土高原地區各個分區降雨侵蝕力差距較大，不同分區降雨侵蝕力由大到小依次為：土石山區>河谷平原區>丘陵溝壑區>高原溝壑區>沙地沙漠區>農灌區。1901—2016年，降雨侵蝕力增加的區域主要集中在黃土高原中部及西部部分區域，減少的區域主要分布于黃土高原西部和東北部。降雨侵蝕力的顯著性變化在空間分布上呈現出邊緣多為非顯著變化而中部多為顯著變化的規律。

(3) 1901—2016年黃土高原地區年平均降雨侵蝕力約為1 462.17 MJ·mm/(hm2·h)。時間序列上6個分區與黃土高原的變化趨勢基本一致。研究時段內，黃土高原降雨侵蝕力變化不顯著且無明顯突變點，降雨侵蝕力變化總體可分為1901—1930年、1930—1980年和1980—2016年3個階段。