基于M-K檢驗和地統計分析的沂蒙山區降雨侵蝕力時空變化趨勢研究

齊 斐, 張春強, 劉 霞, 于海鵬, 趙傳普, 吳 傲, 賈 敏

(1.南京林業大學 南方現代林業協同創新中心, 江蘇省水土保持與生態修復重點實驗室, 江蘇 南京 210037; 2.淮河水利委員會 淮河流域水土保持監測中心站, 安徽 蚌埠 233001)

土壤侵蝕是世界范圍的生態問題，而降雨是土壤侵蝕的重要驅動因子[1]。降雨侵蝕力不僅可定量反映降雨引發土壤侵蝕的潛在能力，也是通用土壤流失方程(USLE)、中國土壤流失方程(CSLE)等土壤流失計算模型的重要因子之一[2-3]。其最早由Wischmeier[4]提出，之后英國、前蘇聯、日本等國家開始了相關研究[5]；從20世紀80年代起，中國在全國各地對EI30經典算法[4]進行驗證，并提出區域降雨侵蝕力計算模型。由于經典算法以次雨量數據為基礎，獲取難度大，處理繁瑣，基于日降雨、月降雨、年降雨等降雨侵蝕力簡易算法的研究逐漸成為R值研究的重心[6]。眾多學者在國內外不同區域開展了豐富的研究[7-16]。其中，國內以章文波等[15]的日雨量模型應用較為廣泛，后來又根據冷暖季雨型的變化，對參數進行調整[17]。受氣候變化影響，降雨量的時空演變會對降雨侵蝕力產生影響，進而影響區域土壤流失風險[18]。國內外學者采用M-K檢驗、小波分析、回歸分析等方法，針對降雨侵蝕力的年際變化趨勢做了諸多研究[16, 19-21]。肖蓓等[22]采用山東省34個氣象站點研究發現1961—2015年年降雨量和降雨侵蝕力總體呈波動下降趨勢，但未進行空間上年際變化分析；馮若昂等[23]采用山東省23個氣象站資料研究發現1966—2015年降雨侵蝕力總體呈不明顯增加趨勢，局部上升或下降趨勢明顯。因此，不同的降雨數據站點數量和空間位置、不同的降雨數據年限可能會對區域降雨侵蝕力時空分析產生重要影響。

沂蒙山國家級水土流失重點治理區位于山東省中南部，地形破碎，土壤松散，是淮河流域易發生土壤侵蝕的區域之一。因此，本研究以沂蒙山國家級水土流失重點治理區為研究區，基于71個雨量站點數據，采用冷暖季日雨量公式計算降雨侵蝕力，分析其年內、年際變化特征，并采用經驗貝葉斯克里金和徑向基函數插值探索空間分異特征及其變化方向，從而為區域土壤流失監測和水土保持工作提供數據支撐。

1 材料與方法

1.1 研究區概況

研究區沂蒙山國家級水土流失重點治理區位于山東省中部和南部(116.74°—119.66°E，36.76°—34.90°N)，面積2.63×104km2。該區域地貌以低山丘陵為主，占81.24%，地形破碎，易發生水土流失。氣候屬半濕潤暖溫帶季風氣候，四季分明，雨熱同期，年平均氣溫在12～14 ℃之間，年均降水量743.52 mm。河流主要有沂河、沭河、泗河等，建有田莊、跋山、岸堤、唐村、日照等大型水庫。主要土壤類型有潮土、粗骨土、褐土、棕壤、紅黏土等；地帶性植被屬暖溫帶闊葉林帶，主要為次生林。

1.2 數據來源與處理

現有相關研究中，多采用氣象站點進行降雨侵蝕力的研究，氣象站點時間序列長，但數量少；而水文雨量站點數量多，但部分站點僅汛期進行降雨量觀測。本研究搜集研究區及其周邊雨量站點101個1980—2018年逐日雨量數據，經統計篩選，剔除汛期雨量站點和觀測年限不足站點，選用71個雨量站點進行降雨侵蝕力研究，統計各站點年/月降雨量、年/月侵蝕性雨量(日雨量≥10 mm)、汛期降雨量(6—9月)等降雨特征數據。雨量站點空間位置見圖1。

圖1 沂蒙山區雨量站點分布

1.3 降雨侵蝕力計算方法

采用冷暖季日降雨量公式[17]進行降雨侵蝕力因子估算，該模型基于全國16個氣象站逐分鐘和逐日降雨資料建立，并針對暖季和冷季的雨型差異采用不同參數，公式如下：

(1)

(2)

(3)

(4)

1.4 降雨侵蝕力和降雨量統計分析方法

1.4.1 年內變化集中性分析 采用FI(fournier index)和CI(concentration index)指數[24]進行降雨侵蝕力和降雨量的集中性分析。其中，當FI<50%時，月份集中度低；50%～100%，月份集中度較高；200%以上，月份集中度高。當CI<15%時，輕微月度集中性；15%～20%時，中度月度集中性；20%～50%時，高度月度集中性；50%以上時，月度集中性非常高。

(5)

(6)

式中：Mmax為全年最大月降雨侵蝕力或最大月降雨量;Mi為第i月降雨侵蝕力或降雨量,i=1,2,…,12。

1.4.2 年際變化趨勢及突變分析 采用5 a滑動曲線、累積距平曲線和Mann-Kendall(M-K)檢驗進行降雨侵蝕力年際變化趨勢分析[21]。其中，M-K趨勢檢驗[21, 25]可能有效區分數據自然波動或確定的變化趨勢，常用于降水、干旱頻次檢查。本研究將其用于降雨侵蝕力的年際變化趨勢定量分析，其檢驗統計量z值為正則表示成增加趨勢，反之，則降低，若|z|≥1.96則表示在95%的區間內有顯著性變化。M-K突變檢驗[26]則通過計算每個樣本的秩，計算統計量UFk和UBk，若UFk和UBk曲線出現交點，且交點在臨界直線(±1.96)之間，那么交點對應的時刻就是突變點。

1.5 基于地統計分析的空間差異性

采用ArcGIS地統計分析中經驗貝葉斯克里金和徑向基函數插值進行空間變異性分析。其中，經驗貝葉斯克里金法用于降雨侵蝕力和降雨量的空間插值，徑向基函數用于M-K檢驗z值和變異系數的空間插值，并采用交叉驗證精度為標準選取最優的函數模型。

經驗貝葉斯克里金法[27](empirical Bayesian Kriging)可通過構造子集和模擬的過程，自行構建克里金模型參數的地統計插值方法，標準誤差優于其他克里金方法，且參數設置方便、效果較好。徑向基函數法(radial basis function)[28]是用一個通過所有屬性樣點的曲面來推算待插值點，可推算大于最大測量值和小于最小測量值的確定性插值方法，方法靈活、精度較高。

2 結果與分析

2.1 降雨時空變化特征

2.1.1 降雨年際、年內變化分析 根據1980—2018年研究區降雨量年際變化特征值統計表(表1)和年際變化曲線(圖2a)，多年平均降雨量為743.52 mm，其中多年平均侵蝕性降雨量為598.84 mm，占降雨量的80.54%。年降雨量和侵蝕性雨量變異系數均高于0.2，表示年度差異較大，結合圖2a和z值，雖年際降雨量上下波動較大，但仍有上升趨勢，尤其是年降雨量和侵蝕性雨量z值高于1.65，表示在90%區間顯著性上升。多年平均侵蝕性降雨頻次為21.07 d，z值2.01，在95%區間呈顯著性上漲趨勢，而多年平均降雨頻次為72.87 d，z值僅0.4，上升趨勢不明顯。由此可知，1980—2018年降雨頻次沒有明顯增加，但侵蝕性降雨頻次出現顯著性增長，單場次降雨量增加。

表1 沂蒙山區降雨量年際變化特征值統計

根據多年平均降雨量和侵蝕性降雨量年內變化曲線(圖2)，降雨量和侵蝕性降雨量呈明顯“幾”字形分布。其中，1月降雨量最低為9.22 mm，7月降雨最高199.88 mm，汛期(6—9月)降雨量占72.56%；多年侵蝕性降雨量均值598.84 mm，占年降雨量的80.54%，1月最低(3.33 mm)，7月最高(199.88 mm)，在24個半月中占各統計時段降雨量的14.08%～89.44%，其中在5月上旬至10月上旬，侵蝕性雨量占比均在75%以上，6月下旬至8月下旬均在85%以上。根據降雨量和侵蝕性降雨量FI指數(53.73，53.01)和CI指數(0.17，0.19)，均具有中度月集中性。

a 年際變化曲線 b 年內變化曲線

2.1.2 降雨空間分布特征 采用經驗貝葉斯法對多年平均降雨量和侵蝕性降雨量進行空間插值(圖3)。由圖3可知，區域內多年平均降雨量介于581.70～852.63 mm，最高值位于淄川中部，最低值位于泗水北部；主要集中于700～800 mm，占區域總面積的77.90%；800 mm以上的區域占7.75%。多年侵蝕性降雨量介于490.58～676.43 mm，最高值位于東港西部，最低值位于泗水北部；主要集中在550～650 mm，占區域總面積的82.93%；650 mm以上占9.44%。二者均呈現北低南高，西低東高的空間分布趨勢，在西北角出現高值旋渦。

根據多年平均侵蝕性降雨量占多年平均降雨量的比例(圖4)，研究區侵蝕性降雨量比例介于77.21%～85.24%，呈現出北低南高，西低東高的趨勢，與降雨量空間分布趨勢接近，但西北區域未出現高值旋渦。

a 多年平均降雨量 b 多年平均侵蝕性降雨量

2.2 降雨侵蝕力時間變化趨勢

2.2.1 年際變化趨勢及突變分析 降雨侵蝕力和降雨量變化趨勢存在一致性，年降雨侵蝕力與年降雨量、侵蝕性雨量相關系數分別為0.86，0.93(p<0.01)。由圖5可知，年降雨侵蝕力介于1 662.45～5 850.45 MJ·mm/(hm2·h·a)，平均為3 656.87 MJ·mm/(hm2·h·a)，呈波動上升趨勢，平均每年增長26.15 MJ·mm/(hm2·h·a)，平均漲幅0.72%；汛期降雨侵蝕力介于1 058.57～5 084.07 MJ·mm/(hm2·h·a)，平均3 220.68 MJ·mm/(hm2·h·a)，年際變化規律與年降雨侵蝕力基本一致，平均每年增長21.70 MJ·mm/(hm2·h·a)，平均漲幅0.67%。因此，年降雨侵蝕力的上升主要由汛期降雨侵蝕力的增長組成，非汛期降雨侵蝕力增長量較少，但平均漲幅達1.01%。

圖4 沂蒙山區多年平均侵蝕性降雨量占多年平均降雨量比例空間分布

圖5 沂蒙山區1980-2018年降雨侵蝕力年際變化特征

根據累積距平曲線和5 a滑動平均曲線(圖6)，沂蒙山區的降雨侵蝕力年際變化趨勢可分為6個階段：波動下降階段1980—1989，1995—2002，2012—2016，波動上升階段1989—1995，2002—2012，2016—2018年。因此，突變點可能位于1989，1995，2012，2016年。

圖6 沂蒙山區1980-2018年降雨侵蝕力累積距平曲線和5 a滑動平均值變化特征

利用M-K突變檢驗法對突變點做進一步檢驗(見圖7)。沂蒙山區的降雨侵蝕力的UBk和UFk在1995年和2012年均沒有交點，在1989年和2016年存在突變點。由于UFk在1990年之后普遍高于0，呈上漲趨勢，且在2008—2013年超過臨界值1.96，存在顯著性上漲。因此，年降雨侵蝕力和汛期降雨侵蝕力由1989年發生突變，開始上升，且存在5 a的顯著性上漲，將會增加區域水土流失風險，易形成山洪、泥石流等水土流失災害。

圖7 沂蒙山區降雨侵蝕力M-K突變檢驗曲線

2.2.2 年內變化特征及集中度分析 根據研究區24個半月多年平均降雨侵蝕力統計(圖8)，降雨侵蝕力主要集中在汛期(占全年85.92%)，其中7，8月降雨侵蝕力各占33.05%和30.21%；1—4月和11—12月降雨侵蝕力占比極少，僅占全年降雨侵蝕力的4.54%，其中1月份降雨侵蝕力最低，僅8.05 MJ·mm/(hm2·h·a)。根據FI指數(399.88)和CI指數(0.24)降雨侵蝕力具有高度月集中性特征。

圖8 沂蒙山區降雨侵蝕力年內變化特征

2.3 降雨侵蝕力空間變異性

2.3.1 多年平均降雨侵蝕力空間分布特征 采用經驗貝葉斯對多年平均降雨侵蝕力進行插值(圖9)，并統計不同梯度降雨侵蝕力面積(表2)。研究區多年平均降雨侵蝕力介于2 983.39～4 377.48 MJ·mm/(hm2·h·a)，總體呈現北低南高、西低東高的空間分布趨勢，最低值出現在安丘北部，最高值出現在山亭東南部。根據降雨侵蝕力梯度數據統計(表2)，降雨侵蝕力主要集中在3 250～4 000 MJ·mm/(hm2·h·a)，占區域總面積的79.40%；4 250 MJ·mm/(hm2·h·a)以上的區域僅占1.62%。

表2 沂蒙山區降雨侵蝕力面積統計

2.3.2 降雨侵蝕力空間變化趨勢 采用徑向基函數插值對71個站點M-K檢驗值z值和變異系數進行插值(圖10)。研究區z值介于-0.91～2.6，整體以>0為主，尤其是在西北的博山、淄川、臨朐存在明顯的上升趨勢(z值>1.96)，僅在費縣中部、沂南中部、東港西部等少量區域z值小于0，存在不顯著的下降趨勢(圖10a)。

根據圖10b可知，研究區年降雨侵蝕力變異系數介于0.35～0.79，均為中等變異，說明該區域降雨侵蝕力年際差異較大。其中，莒南整體變化相對較小，而東北部淄川—臨朐—安丘一帶變異系數較高。綜合考慮z值和變異系數，在西北部和北部z值和變異系數均較高，說明該區域年際降雨侵蝕力波動性較高，且上漲趨勢比較明顯。

圖9 沂蒙山區多年平均降雨侵蝕力空間分布

3 討論與結論

3.1 討 論

在本研究中，降雨侵蝕力總體呈不顯著上漲趨勢，局部呈現顯著性增長，與劉正佳等[29]、馮若昂等[23]、馬良等[28]降雨侵蝕力變化趨勢相似，但與肖蓓等[22]結果相反。肖蓓等[22]采用山東省34個氣象站點研究發現1961—2015年年降雨量和降雨侵蝕力總體呈波動下降趨勢，其中，魯中南山地丘陵區13個站點降雨量和降雨侵蝕力年際波動大，但總體降低趨勢輕微，平均每年降低0.43 mm和6.17 MJ·mm/(hm2·h·a)；劉正佳等[29]采用沂蒙山區及其周邊38個氣象站點1971—2008年數據，發現降雨侵蝕力最低的10 a出現在1980—1989年，而后10 a降雨侵蝕力最高；馮若昂等[23]采用山東省23個氣象站研究發現1966—2015年降雨侵蝕力總體呈不明顯增加趨勢，根據其降雨侵蝕力增長分布圖，發現在研究區內，淄川區周邊降雨侵蝕力增量較高，在臨沂市周邊呈下降趨勢，與本研究變化趨勢空間分布比較一致；馬良等[28]采用山東省22個氣象站點1951—2008年降雨數據發現，全省降雨侵蝕力總體未有顯著增減，但在在魯中南山區有明顯升高趨勢。綜合幾位學者研究與本研究結果來看，本研究區范圍內降雨侵蝕力高值年份出現在1964，1991，2003年，而2014年后降雨侵蝕力呈上升趨勢，但高值尚未出現，可能是肖蓓等[22]研究中出現輕微下降的原因。綜上所述，除了降雨數據站點的數量和空間分布位置，降雨的周期性變化也會對降雨侵蝕力的變化趨勢造成影響，可進一步延長降雨時間年限，分析降雨周期，并以此為單元進行降雨侵蝕力的變化趨勢分析。

a M-K趨勢檢驗z值 b 變異系數

3.2 結 論

(1) 降雨量和侵蝕性降雨量成中度月集中性，降雨侵蝕力具有高度月度集中性，均主要集中于7—8月份，呈單峰形。

(2) 區域多年平均降雨侵蝕力介于2 983.39～4 377.48 MJ·mm/(hm2·h·a)，平均3 656.87 MJ·mm/(hm2·h·a)，以3 250～4 000 MJ·mm/(hm2·h·a)分布范圍最廣，占79.40%；多年平均降雨量介于581.70～852.63 mm，平均743.52 mm，以700～800 mm分布范圍最廣，占77.90%；總體上呈現南高北低、東高西低的趨勢。

(3) 降雨侵蝕力和降雨量、侵蝕性雨量1980—2018年的年際波動規律相似，呈極顯著相關(p<0.01)。研究區總體年降雨侵蝕力、汛期降雨侵蝕力年際變化屬中度變異，具有不顯著的增長趨勢；年降雨量和侵蝕性降雨量年際變化為中度變異，且在90%區間顯著性增長。在空間分布上，降雨侵蝕力在西北的博山、淄川、臨朐存在明顯的波動上升趨勢(z>1.96)，僅在南部少量區域存在不顯著的下降趨勢(z<0)。