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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)來源與處理

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

圖1 沂蒙山區(qū)雨量站點分布

1.3 降雨侵蝕力計算方法

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

(1)

(2)

(3)

(4)

1.4 降雨侵蝕力和降雨量統(tǒng)計分析方法

1.4.1 年內(nèi)變化集中性分析 采用FI(fournier index)和CI(concentration index)指數(shù)[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]可能有效區(qū)分數(shù)據(jù)自然波動或確定的變化趨勢，常用于降水、干旱頻次檢查。本研究將其用于降雨侵蝕力的年際變化趨勢定量分析，其檢驗統(tǒng)計量z值為正則表示成增加趨勢，反之，則降低，若|z|≥1.96則表示在95%的區(qū)間內(nèi)有顯著性變化。M-K突變檢驗[26]則通過計算每個樣本的秩，計算統(tǒng)計量UFk和UBk，若UFk和UBk曲線出現(xiàn)交點，且交點在臨界直線(±1.96)之間，那么交點對應的時刻就是突變點。

1.5 基于地統(tǒng)計分析的空間差異性

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

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

2 結(jié)果與分析

2.1 降雨時空變化特征

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

表1 沂蒙山區(qū)降雨量年際變化特征值統(tǒng)計

根據(jù)多年平均降雨量和侵蝕性降雨量年內(nèi)變化曲線(圖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個半月中占各統(tǒng)計時段降雨量的14.08%～89.44%，其中在5月上旬至10月上旬，侵蝕性雨量占比均在75%以上，6月下旬至8月下旬均在85%以上。根據(jù)降雨量和侵蝕性降雨量FI指數(shù)(53.73，53.01)和CI指數(shù)(0.17，0.19)，均具有中度月集中性。

a 年際變化曲線 b 年內(nèi)變化曲線

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

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

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

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

2.2.1 年際變化趨勢及突變分析 降雨侵蝕力和降雨量變化趨勢存在一致性，年降雨侵蝕力與年降雨量、侵蝕性雨量相關系數(shù)分別為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)，年際變化規(guī)律與年降雨侵蝕力基本一致，平均每年增長21.70 MJ·mm/(hm2·h·a)，平均漲幅0.67%。因此，年降雨侵蝕力的上升主要由汛期降雨侵蝕力的增長組成，非汛期降雨侵蝕力增長量較少，但平均漲幅達1.01%。

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

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

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

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

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

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

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

圖8 沂蒙山區(qū)降雨侵蝕力年內(nèi)變化特征

2.3 降雨侵蝕力空間變異性

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

表2 沂蒙山區(qū)降雨侵蝕力面積統(tǒng)計

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

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

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

3 討論與結(jié)論

3.1 討 論

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

a M-K趨勢檢驗z值 b 變異系數(shù)

3.2 結(jié) 論

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

(2) 區(qū)域多年平均降雨侵蝕力介于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%；總體上呈現(xiàn)南高北低、東高西低的趨勢。

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