許功偉, 徐立榮, 張云蘋, 王 坤, 惠 莉, 徐 晶, 趙芹蕊

(1.濟南大學 水利與環(huán)境學院, 濟南 250022; 2.禹城市水利局, 山東 德州 251200; 3.濟南黃河河務局 章丘黃河河務局, 濟南 250200)

在世界范圍內，土壤侵蝕問題早已備受關注，我國土壤侵蝕問題也尤為突出[1]。據研究表明，降雨會破壞土壤的結構，其形成的匯流會帶走土壤顆粒，這就直接造成了土壤的流失，而雨滴大小分布和雨滴終端速度決定的降雨侵蝕力能反映出降雨造成土壤侵蝕的潛力[2-4]。降雨侵蝕力的計算方法從1958年到現在已經提出多種。最初是1958年Wischmeier等[5]提出的利用次降雨動能和降雨強度的乘積來計算，后來又由科學家提出用次降雨量和降雨強度乘積來計算[6]，盡管此方法的計算避免了降雨動能的計算，但同樣需要詳細的次降雨量數據，并且計算的工作量繁重。再后來，科學家們提出了簡易估算模型，從最初孫保平等[7]建立并在西吉縣黃土丘陵區(qū)試驗的年雨量模型到后來的黃炎和等[8]提出的在閩東南地區(qū)試驗的月雨量模型，再到現在使用廣泛的章文波等[9]提出的通過日雨量計算降雨侵蝕力的日雨量模型，這就使得降雨侵蝕力的計算量得到顯著的下降，也使得精度逐漸的提高。

大汶河流域地勢東西差異大，從東至西落差超300 m。大汶河流域受季風影響，夏季暴雨頻繁，其他三季的降雨相對較少。受地勢和氣候的雙重影響，大汶河流域自建國以來就洪澇不斷，土壤侵蝕問題也比較嚴重，不僅影響日常生活，也阻礙經濟發(fā)展。這些現象表明，想要讓大汶河流域突破種種發(fā)展阻礙，首先要解決其土壤侵蝕問題，從大汶河流域的背景來看，要解決土壤侵蝕問題，對大汶河流域的降雨侵蝕力進行分析就顯得尤為重要。

因此，本文基于1970—2019年大汶河流域6個市級雨量站的日雨量數據和日降雨侵蝕力模型來計算降雨侵蝕力，并借助ArcGIS，MATLAB等工具分析其時間和空間的變化，結合前人對大汶河流域的研究，尋求大汶河流域中降雨侵蝕力的演化規(guī)律，共同為緩解大汶河流域的土壤侵蝕問題提供科學數據依據，以此來有效地緩解水土侵蝕對大汶河流域的影響。

1 研究區(qū)概況

大汶河(35°45′—36°29′N，116°16′—117°53′E)位于山東省中部，全長209 km，流域面積8 633 km2。大汶河流域屬大陸性季風氣候，旱澇季分明，夏季降雨多，而春、秋、冬季則降雨較少。大汶河流域地形地貌主要分為平原、丘陵、山地，還有少部分的洼地、水面。北部及東部的上游地區(qū)是山地，而其下游地區(qū)的海拔較低，從東北至西南方向地勢逐漸降低，這種地勢分布就造成了大汶河流域河流呈自東向西走向的走勢。其植被覆蓋，受暖溫帶影響，山區(qū)植被覆蓋相對較少，以松、柏等植物為主，中游等丘陵地區(qū)則以灌木居多，下游主要種植梧桐、楊柳等植物。由于這種特殊的氣候和地形以及植被，就使得此流域洪澇災害、水土流失等現象極易發(fā)生。

2 材料與方法

2.1 數據來源

本文數據采用1970—2019年大汶河流域6個雨量站(圖1)的日降雨數據，數據來源于中國氣象數據共享服務系統(http:∥cdc.cma.gov.cn)。

圖1 大汶河流域地形及雨量站分布

2.2 研究方法

2.2.1 降雨侵蝕力的計算 因為采用降雨動能的方法來計算，數據很難獲取，并且處理十分繁瑣，而以月或年的降雨量數據來計算降雨侵蝕力則精度不高，因此本文采用章文波等[9]提出的方法來計算。此方法已被廣泛應用[10-12]。具體計算公式如下：

(1)

式中:Ri為第i個半月的降雨侵蝕力〔MJ· mm/(hm2·h·a)〕;n為半月內的降雨量≥12 mm的日數(d);Pk為半月內第k天的降雨量≥12 mm的降雨量(mm)。α和γ是模型參數,其計算如下:

(2)

α=21.586γ-7.1891

(3)

式中:Pd12為日均侵蝕性雨量(mm);Py12是年均侵蝕性雨量(mm)。

2.2.2 時間變化分析 文中主要用了3種方法對時間變化進行分析，分別分析了時間序列的趨勢變化、突變的年份以及周期的變化，具體方法如下：

(1) Mann-kendall法。此法分析時間序列的變化趨勢時，是通過統計值的改變來分析，統計值Z為正時，表現為上升的趨勢，當超過置信區(qū)間±1.96(α=0.05)時，則為顯著上升，反之同理。此法還可用來判斷系列突變開始的時刻，主要通過兩個變量來判斷，按時間順序統計的UF和其逆序列UB，兩者的交點，就是突變發(fā)生的時刻[13-14]。

(2) 累積距平法。很多時候單用Mann-kendall法分析突變點時，圖中會出現多個交點，無法判斷，這時就需要通過累積距平法與其共同分析出突變點。累積距平法的原理是通過求出此序列的平均值，然后用逐年的數據與平均值做差，得出的便是距平值，然后逐年累加，就得出累積距平序列[15-17]。公式如下：

(4)

(3) 小波分析。分析時間序列的問題時，基礎的時域與頻域的分析往往不能滿足一些受多種因素影響的現象，因此就需要通過小波分析來進行分析，結合時域與頻域，在時間的變化上分析其周期的改變性來滿足分析的需要。其分析流程是通過MATLAB來計算小波系數，然后通過Excel計算小波系數的實部，最后通過Furfer 9繪制出小波等值線圖進行分析。通過分析圖中枯豐的交替來判斷周期的變化規(guī)律[18-19]。

2.2.3 空間變化分析 空間變化的分析主要是通過反距離加權插值等方法來進行的。反距離加權插值法依據的原理是，彼此相近的兩者相較于遠處的事物更具有相似性。在通過此法計算時，距離觀測者較近的觀測物會被賦予較大的權重，較遠則相反。因此權重隨距離的增加而減少，故名為反距離加權插值。分析過程則是通過ArcGIS 10.4對大汶河流域的年降雨量和年降雨侵蝕力進行空間插值，得到空間插值分布圖進行具體分析。

3 結果與分析

3.1 降雨侵蝕力時間變化分析

3.1.1 年際變化 通過數據分析，大汶河流域1970—2019年年均降雨侵蝕力范圍為1 310.84～6 721.53 MJ·mm/(hm2·h·a)，均值為3 808.83 MJ·mm/(hm2·h·a)；年均降雨量為735.24 mm，最大值為1990年出現的1 135.47 mm，最小值為2002年的415.88 mm；并且降雨侵蝕力最值出現的年份與降雨量相同，驗證了降雨量對降雨侵蝕力的作用。

據圖2分析，年降雨量的變化趨勢為：1989年前呈下降趨勢，1989年波動上升至1994年，從1994年又下降至2002年、2002年上升至2007年后，繼續(xù)波動下降至2019年。年降雨侵蝕力的變化與其大體一致。圖2A中，降雨量的5 a移動平均呈微弱上升趨勢，但在圖2B中，降雨侵蝕力的趨勢卻是存在著下降的狀態(tài)，表明雖然降雨量在逐漸增加，但其侵蝕性雨量卻存在一個逐步下降的趨向。由于人口增長及人類活動的影響，全球氣溫正在不斷變暖，尤其近10 a來增加更為劇烈，全球氣溫變暖使得水循環(huán)加劇，大氣中水分含量增加，進而積雨云出現的概率增加[20]，大汶河流域東、北部是山地，研究站點大部分位于雨影區(qū)與平原區(qū)，就使得研究區(qū)容易出現強對流天氣，進而降雨量增加。降雨侵蝕力受侵蝕性降雨影響，大汶河流域侵蝕性降雨大多集中在夏、秋兩季。目前，全球氣溫變暖，ENSO事件不斷發(fā)生，ENSO暖事件使得大汶河流域夏、秋季降雨量減少，導致大汶河流域夏、秋季降雨侵蝕力減小，從而使得整體流域降雨侵蝕力減少。

圖2 大汶河流域年降雨量和年降雨侵蝕力的變化趨勢

對大汶河流域年降雨侵蝕力和年降雨量的離散程度進行分析，得出前者的變異系數CV值為0.34，后者為0.21，此二者同屬中等變異。但相較之，前者的變異程度大于后者，說明降雨侵蝕力在此分析時段內變化的更強烈。兩者通過Mann-kendall法求得的Z值分別為-1.18，0.05，表明在此分析時段內年降雨侵蝕力的變化趨勢為下降，而年降雨量的變化趨勢為上升，這和5 a移動平均的分析結果一致。盡管與年降雨量相比，年降雨侵蝕力體現出較大的變化趨勢，但均小于1.96，變化均為不顯著。

圖3為周期分析圖，從圖中可以分析出，在演化過程中，降雨侵蝕力出現了3類周期，分別為4～7 a，8～12 a和17～25 a。其中，8～12 a周期上存在著“枯—豐”交替的準7次振蕩，17～25 a存在4次。同時還能分析出，8～12 a和17～25 a這兩類周期在全局都比較穩(wěn)定，而4～7 a這類的周期只在1988—2008年的變化較為穩(wěn)定。

圖3 基于小波分析的年降雨侵蝕力的周期變化

圖4和圖5為突變分析的結果，在圖4中，1976年之前UF>0，1976—1978年UF<0，1978—1980年UF>0，1979年之后UF<0，說明1976年之前都呈上升趨勢，1976—1980年發(fā)生短暫的先下降后上升趨勢，1980年直至最后都呈下降趨勢。UF與UB相交于1976年和1979年，且兩個交點都位于±1.96(α=0.05)內，為了準確確定突變年份，結合累積距平法來具體分析。發(fā)現，圖5曲線的變化趨勢與圖4基本相同，且在1976年曲線沒有發(fā)生顯著的突變變化，而在1979年以后曲線明顯呈下降趨勢，這與圖4中1979年之后曲線明顯下降，并且在1984年以后超出置信區(qū)間(變化顯著)的分析結果恰好相同。通過上述兩種方法共同證明了突變發(fā)生在1979年。

圖4 年降雨侵蝕力M-K突變分析

圖5 年降雨侵蝕力累積距平曲線

3.1.2 年內變化 從圖6中可以發(fā)現降雨侵蝕力和降雨量的年內分布大致相同，主要分布在6—9月。其中，7月、8月尤為突出。這種現象主要是由于受東南季風的影響，夏季陸地溫度高于海洋，導致陸地空氣的密度小于海洋，從而使得空氣從氣壓相對較高的海洋吹向氣壓相對較低陸地，形成季風，造成集中降水。此時間段內降雨量和降雨侵蝕力分別占全年的82.78%，73.58%，前者的最大值出現在7月，為208.10 mm，占全年的28.25%，同樣，后者的最大值出現的時間與之相同，其值為1 288.72 MJ·mm/(hm2·h·a)，占全年的33.82%；其次是8月，前后二者的值分別為166.67 mm和1 000.55 MJ·mm/(hm2·h·a)，分別占全年的22.62%，26.26%；降雨侵蝕力和降雨量的最低值都出現在1月，分別為3.77 MJ·mm/(hm2·h·a)和6.47 mm，各占全年的0.10%，0.88%；7月、8月的降雨侵蝕力和降雨量所占全年比例分別為60.08%，50.87%，都占到全年的半數以上。通過年內變化的分析表明在7月、8月加大水土保持的力度對緩解大汶河流域水土侵蝕問題顯得尤為重要。

圖6 大汶河流域月降雨量和降雨侵蝕力分布

3.1.3 季節(jié)變化 圖7為大汶河流域降雨侵蝕力四季變化圖，據分析，春季到冬季的變化趨勢為，從顯著上升到微弱下降再到顯著下降最后呈顯著上升趨勢，除了夏季變化不顯著，春、秋、冬季均為顯著變化，夏、秋兩季呈現出減少趨勢(這與3.1.1年際變化中的分析結果相同)，主要是因為全球氣溫逐漸升高，使得ENSO事件發(fā)生的頻率增加，而經馮若昂等[21]發(fā)現ENSO暖事件會對會山東省內的降雨侵蝕力造成減少的影響。春季到冬季的平均R值依次為398.13，2 744.96，635.98，29.75 MJ·mm/(hm2·h·a)，各占全年降雨侵蝕力的10.45%，72.07%，16.70%，0.78%。夏季明顯大于其他三季，而冬季明顯低于其他三季，這與夏季降雨多、冬季降雨少有直接關系。通過對季節(jié)的變化的分析，可以看出夏季降雨侵蝕力較為集中，還可以看出夏季降雨侵蝕力的變化相較于其他季節(jié)較為穩(wěn)定，集中且穩(wěn)定的變化就可針對性的在每年的7月、8月加大水土保持工作的力度，這就大大減少了水土保持工作的強度。

圖7 大汶河流域降雨侵蝕力四季變化

3.2 降雨侵蝕力空間變化分析

3.2.1 降雨量和降雨侵蝕力的空間分布分析 由上述時間分析發(fā)現，降雨量的改變確實在某種程度上影響著降雨侵蝕力，因此在空間上分析降雨侵蝕力，也首先對降雨量進行分析。通過圖8A發(fā)現泰山周圍的降雨量最高，基本都在900 mm之上，其次是萊蕪、新泰，二者降雨量在730 mm以上，再次是肥城、泰安，也都超越650 mm，最低是東平，不足650 mm。大汶河流域降雨整體分布為中部最多，明顯高于東、西部，再者是東部，西部最低。通過圖8B可以看出，降雨侵蝕力的分布與前面分析的降雨量的分布有著很高的相似度，R值方面，中部地區(qū)最高，其值可達4 000 MJ· mm/(hm2·h·a)，以泰山為中心的地區(qū)甚至達到5 000 MJ·mm/(hm2·h·a)，東部整體高于西部。這種現象與大汶河流域的地勢分布有關，其北部及東部是山地，地勢較高，這就使得氣流抬升，然后形成降雨，造成降雨較多，從而產生相對較大的降雨侵蝕力。因此在大汶河流域的水土保持工作中，強化以泰山為中心的山地地區(qū)的水土保護措施，會使其工作事半功倍。

圖8 大汶河流域降雨量和降雨侵蝕力的空間分布

3.2.2 年降雨侵蝕力空間變化分析 通過對Z值和CV值的空間插值圖(圖9)進行分析，可以看出大汶河流域大多數雨量站的Z值都小于0，說明分析區(qū)站點的降雨侵蝕力都存在下降的趨向，盡管東平、萊蕪等東、西部地區(qū)的下降趨勢比泰山、泰安等中部地區(qū)的下降趨勢大一些，但都沒超過±1.96，說明均為不顯著下降，肥城的Z值恰好為0，這就表明肥城地區(qū)的降雨侵蝕力在近50 a內的變化都較為穩(wěn)定。變異系數的分布與Z值分布相似，都是東平最大，其次是新泰、泰安和萊蕪，最低為肥城和泰山，整體上看中部較小，東、西部較大。說明降雨侵蝕力在東平、萊蕪等東、西部地區(qū)不僅下降趨勢較大，而且離散程度也較大。

圖9 大汶河流域降雨侵蝕力Z值和變異系數的空間分布

3.3 不同地形區(qū)降雨侵蝕力分析

根據大汶河流域地形特點將大汶河流域劃分為東、北部山地地區(qū)和中、西部平原地區(qū)兩種地形區(qū)(以200 m海拔為界)，進而分析降雨侵蝕力與Z值、變異系數等的關系(表1)。

通過表1可以看出，東部、北部山地地區(qū)的降雨量和降雨侵蝕力明顯高于中、西部平原地區(qū)，這是因為從海洋吹向陸地的暖濕氣團遇到較高地勢會造成氣流抬升，從而使降雨增多，侵蝕性降雨也隨之增加，進而表現出較大的降雨侵蝕力。還可以看出東、北部山地地區(qū)的|Z|比中、西部平原地區(qū)大，說明其下降趨勢比中、西部平原地區(qū)明顯。兩類地形區(qū)和整體流域的變異系數依次為0.33，0.36，0.34，盡管中、西部平原地區(qū)的變異系數最大，但三者整體相比，差異不大，離散程度相當。在相關性方面，降雨侵蝕力與降雨量，除了中、西部為0.64，其余均超0.90，呈顯著正相關，其與海拔的相關性同樣如此，且都超過0.90，呈強相關性。

表1 大汶河流域不同地形區(qū)降雨侵蝕力分析

4 討 論

大汶河流域年際降雨量與降雨侵蝕力分別表現出了微弱上漲和不顯著下降的趨勢。降雨侵蝕力的變化與馮若昂等[21]在分析山東省降雨侵蝕力的變化中提到以泰山為中心一帶的年降雨侵蝕力的變化趨勢以及肖蓓等[22]提到的魯中南地區(qū)降雨侵蝕力變化趨勢相符；降雨量的變化趨勢與朱祎真等[23]分析的相同，也與徐晶等[13]分析其鄰近濟南市的變化相一致，但與夏繼勇等[24]分析的大汶河流域近58 a降雨量的變化趨勢不同，主要是由研究中采取不同長度的降雨序列(朱祎真等研究序列區(qū)間為1971—2010年，徐晶等研究序列區(qū)間為1970—2016年，夏繼勇等研究序列區(qū)間為1956—2013年，本研究序列區(qū)間為1970—2019年)所造成的。

ENSO冷事件，會增加山東省夏秋季降水，從而使降雨侵蝕力增加，ENSO暖事件則相反[21]。ENSO冷事件一般伴隨在ENSO暖事件之后發(fā)生，但由于全球氣溫正在不斷變暖，ENSO暖事件的發(fā)生次數遠多于ENSO冷事件。因此，大汶河流域受ENSO事件影響，主要表現為暖事件。由于地理位置的原因，大汶河流域的降水主要受東南季風的影響，降雨集中在夏秋季，而ENSO暖事件會減弱東南季風，減少對大陸的水汽輸送，從而減少流域夏秋季的降雨，這就使得流域夏秋季的侵蝕性降雨也隨之減少，進而流域夏秋季降雨侵蝕力減小，而夏秋季降雨侵蝕力又是整年降雨侵蝕力的主要來源，因此，大汶河流域整體降雨侵蝕力年際變化呈下降趨勢。

5 結 論

(1) 大汶河流域1970—2019年年均降雨侵蝕力范圍為1 310.84～6 721.53 MJ·mm/(hm2·h·a)，均值為3 808.83 MJ·mm/(hm2·h·a)。降雨侵蝕力的變化存在著4～7 a，8～12 a，17～25 a共3類周期，5 a移動平均顯示其變化趨勢為微弱下降，雖然整體有著減少的趨向，但年際變化比較強烈，并且在1979年發(fā)生突變。

(2) 按月分析，降雨侵蝕力和降雨量在7月、8月的分布尤為突出。一年中，7月、8月份的降雨侵蝕力和降雨量占比高達60.08%，50.87%。四季的年際趨勢變化為，春季顯著上升，夏季微弱下降，秋季顯著下降，冬季顯著上升，除了夏季變化不顯著，春、秋、冬季均為顯著變化。7月、8月兩月的降雨侵蝕力，不僅集中而且穩(wěn)定，加大這兩個月的防控力度是整個水土保持工作的重中之重。

(3) 大汶河流域降雨整體分布為中部最多，明顯高于東、西部，再者是東部，西部最低，降雨侵蝕力的分布與之類似；大汶河流域整體雨量站的Z值都小于0，說明分析區(qū)內站點的降雨侵蝕力都存在下降的趨向，盡管東平、萊蕪等東、西部地區(qū)的下降趨勢比泰山、泰安等中部地區(qū)的下降趨勢大一些，但都沒超過±1.96，說明均為不顯著下降；大汶河流域離散程度的分布與Z值類似，東、西部離散程度較大，中部較小。

(4) 東部、北部山地地區(qū)的降雨量和降雨侵蝕力高于中、西部平原地區(qū)，且東、北部降雨侵蝕力的下降趨勢也比中、西部明顯，但離散程度相當，差異不大；降雨侵蝕力與降雨量、海拔在不同地區(qū)都有較強的關聯性，這表明在東、北部地勢較高的地區(qū)要注重水土流失的防控，也再次驗證了在7月、8月雨季集中的時候要加大水土保持的力度的結論。