基于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)。