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

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

圖1 研究區(qū)位置

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

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

1.3 研究方法

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

(1)

(2)

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

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

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

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

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

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

(3)

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

2 結(jié)果與分析

2.1 降水?dāng)?shù)據(jù)集精度驗證結(jié)果

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

表1 CHELSAcrust數(shù)據(jù)集逐月降水量精度驗證結(jié)果

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

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

在20世紀20年代之前，黃土高原平均降雨侵蝕力為1 233.35 MJ·mm/(hm2·h)，其中土石山區(qū)的降雨侵蝕力最大，為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年，土石山區(qū)以及丘陵溝壑區(qū)東部的極強降雨侵蝕力[>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倍以上，各分區(qū)在這一時期的降雨侵蝕力達到峰值，且增大趨勢向西北地區(qū)的高原溝壑區(qū)和沙地沙漠區(qū)延伸。20世紀80年代—21世紀初，降雨侵蝕力逐漸下降，降雨侵蝕力<2 000 MJ·mm/(hm2·h)的面積占到整個黃土高原面積的79.04%～88.59%，而降雨侵蝕力>4 000 MJ·mm/(hm2·h)的面積僅占0.57%～0.16%，且降雨侵蝕力較大的地區(qū)逐漸退至黃土高原東部的土石山區(qū)局部地區(qū)。2010—2016年，降雨侵蝕力>4 000 MJ·mm/(hm2·h)面積占比增至3.06%，表明黃土高原降雨侵蝕力高值區(qū)域再次增加。

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

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

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

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

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

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

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

圖5 1901-2016年黃土高原及各分區(qū)年均

3 結(jié) 論

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

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

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

圖6 1901-2016年黃土高原及其分區(qū)年降雨侵蝕力Mann-Kendall突變檢驗曲線

圖7 1901-2016年黃土高原降雨侵蝕力EMD分解結(jié)果

表2 黃土高原降雨侵蝕力時間序列IMF分量及其方差貢獻率

(4) 黃土高原地區(qū)降雨侵蝕力變化存在的周期性規(guī)律，2.62 a變化周期最顯著，其次為5.29 a，再其次為11.89，31.00，70.00 a，此3個周期都為弱周期。且變化周期與氣候要素的波動周期基本一致。