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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

圖1 研究區(qū)位置

1.2 數(shù)據(jù)來(lái)源

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

1.3 研究方法

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

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

1.3.3 降雨時(shí)空分布特征分析方法

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

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

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

(3)

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

2 結(jié)果與分析

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

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

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

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

2.2.1 年均降雨侵蝕力空間分布 1901—2016年間黃土高原年均降雨侵蝕力差距較大，從東南向西北遞減(圖2)，東部大部分地區(qū)降雨侵蝕力超過(guò)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世紀(jì)20年代之前，黃土高原平均降雨侵蝕力為1 233.35 MJ·mm/(hm2·h)，其中土石山區(qū)的降雨侵蝕力最大，為5 253.75 MJ·mm/(hm2·h)。降雨侵蝕力<2 000 MJ·mm/(hm2·h)的面積占整個(gè)黃土高原面積的80.92%～88.44%，而降雨侵蝕力>4 000 MJ·mm/(hm2·h)的面積僅占總面積的0.34%～1.43%；表明該時(shí)期黃土高原的降雨侵蝕力總體維持在較低水平；20世紀(jì)30年代開(kāi)始，黃土高原的降雨侵蝕力激增，尤其是1930—1939年，土石山區(qū)以及丘陵溝壑區(qū)東部的極強(qiáng)降雨侵蝕力[>5 000 MJ·mm/(hm2·h)]顯著增加，面積占比達(dá)到3.10%(圖2)。1930—1979年，降雨侵蝕力在1 000～3 000 MJ·mm/(hm2·h)的面積占比達(dá)到60%以上，降雨侵蝕力>4 000 MJ·mm/(hm2·h)的面積占比為1.99%～6.53%，比1901—1929年增加4.5倍以上，各分區(qū)在這一時(shí)期的降雨侵蝕力達(dá)到峰值，且增大趨勢(shì)向西北地區(qū)的高原溝壑區(qū)和沙地沙漠區(qū)延伸。20世紀(jì)80年代—21世紀(jì)初，降雨侵蝕力逐漸下降，降雨侵蝕力<2 000 MJ·mm/(hm2·h)的面積占到整個(gè)黃土高原面積的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年黃土高原降雨侵蝕力時(shí)空分布

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

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

2.3 降雨侵蝕力時(shí)間演變特征

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

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

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

2.3.2 降雨侵蝕力周期性分析 利用EMD法對(duì)黃土高原地區(qū)1901—2016年的116 a降雨侵蝕力進(jìn)行分析，得到了5個(gè)本征模分量(IMF)及最低頻分量(RES)(圖7)。表2計(jì)算了各個(gè)IMF分量、方差以及方差貢獻(xiàn)率?？梢钥闯?個(gè)本征模分量(IMF)累計(jì)方差貢獻(xiàn)達(dá)到96.03%，而其中IMF1貢獻(xiàn)率最大，為56.53%，其次為IMF2分量，為32.77%。IMF5分量的貢獻(xiàn)率最小，為1.48%。由此得出，黃土高原地區(qū)降雨侵蝕力變化存在周期性規(guī)律，2.75 a變化周期最顯著，其次為5.33 a，再其次為10.29，24.5，44.5 a，這3個(gè)周期的方差貢獻(xiàn)率差異不明顯，均為弱周期。降雨侵蝕力的周期變化受制于該區(qū)域氣候要素的波動(dòng)周期。太陽(yáng)黑子和大氣環(huán)流的異?；顒?dòng)可以影響水熱平衡，從而改變降水時(shí)序和空間分布。研究表明，太陽(yáng)黑子存在11 a的活動(dòng)周期[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ù)精度較高，可用于黃土高原長(zhǎng)時(shí)段研究需求。

(2) 黃土高原降雨侵蝕力總體為東部大于西部。黃土高原地區(qū)各個(gè)分區(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)。時(shí)間序列上6個(gè)分區(qū)與黃土高原的變化趨勢(shì)基本一致。研究時(shí)段內(nèi)，黃土高原降雨侵蝕力變化不顯著且無(wú)明顯突變點(diǎn)，降雨侵蝕力變化總體可分為1901—1930年、1930—1980年和1980—2016年3個(gè)階段。

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

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

表2 黃土高原降雨侵蝕力時(shí)間序列IMF分量及其方差貢獻(xiàn)率

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