1961-2020年沂蒙山区降雨侵蚀力时空变化特征

徐 晶, 时延锋, 徐征和, 徐立荣

(1.济南大学 水利与环境学院, 济南 250022; 2.山东建筑大学 市政与环境工程学院, 济南 250101)

土壤侵蚀是全球性的环境问题之一，不但导致土壤退化、土地生产力降低，影响农业生产和食物安全，同时也会对相邻地区的人类生存、生态环境和社会经济发展带来严重影响[1-2]。土壤侵蚀的影响因素有很多，其中降雨是最主要的影响因素之一，一定强度的降雨量降落地表后会形成地表径流，进而造成土壤侵蚀，其侵蚀能力的大小用降雨侵蚀力来表示[3-5]。准确计算研究区的降雨侵蚀力，对定量评估土壤侵蚀、流域综合治理等具有重要意义。

降雨侵蚀力的计算最初是由Wischmeier等[6]提出的，计算公式为EI30，其中E代表降雨动能，I30为最大30 min降雨强度，该方法的优势在于计算的精度较高，但同时对数据的要求高，需要长序列的场次降雨，国内外许多地方都很难满足这样的数据要求，并且数据处理起来非常的繁琐、耗时间。因此，许多学者开始探索并提出了利用日、月或年降雨量计算降雨侵蚀力的模型[7-10]。相对于其他模型，基于日降雨量计算降雨侵蚀力的模型精度高，得到了广泛的应用[11-15]。曾瑜等[16]利用章文波的日降雨模型分析了1961—2014年鄱阳湖流域降雨侵蚀力的时空变化特征，得到鄱阳湖流域年均降雨侵蚀力为9 537. 9 (MJ·mm)/(hm2·h·a)，年际变化呈上升趋势，但趋势不显著；刘慧等[17]则利用1961—2015年雅鲁藏布江流域 8个气象站日降雨量气象资料，得到年降雨侵蚀力平均值为 758.1 (MJ·mm)/(hm2·h·a)，总体呈波动上升趋势，空间上东部地区的降雨侵蚀力高于西部地区；索笑颖等[18]分析了河北山区2000—2018年降雨侵蚀力时空变化特征，该地区夏季的水土流失最严重，尤其是燕山山区的部分地区水土流失问题最为突出。

沂蒙山区属于典型的土石山区，山地和丘陵面积广，人口密集，加上表层的土壤较为疏松，土壤水分涵养能力较低，土地利用不合理，使得当地的植被遭到严重的破坏，生态环境脆弱，因此土壤侵蚀问题尤为突出[19-20]。沂蒙山区是山东省内水土流失最为严重的区域之一，因此本文以该区域为研究区，基于日降雨数据，利用日降雨侵蚀力模型，结合Mann-Kendall趋势/突变检验法、累积距平法、小波分析、反距离加权插值(IDW)等方法系统分析1961—2020年沂蒙山区降雨侵蚀力的时空变化演替特征及变化规律，以期为区域的水土流失监测与防治、生态环境保护等提供科学参考。

1 材料与方法

1.1 研究区概况

沂蒙山区位于山东省中南部(116°34′—119°39′E，34°26′—36°23′N)，总面积约3.26×104km2，共辖6个地市，27个县(市、区)(图1)。研究区的气候类型属于半湿润暖温带季风气候，年均降水量为 700～900 mm，降水较为丰富，年平均气温为12～14℃。土壤以棕壤和褐土为主，地貌复杂，以山地、丘陵为主。自然植被以暖温带落叶阔叶林为主，但由于人类活动的影响及破坏，现在自然植被比较少，多以人工植被或次生林为主。

图1 沂蒙山区高程及气象站点分布

1.2 数据来源

本研究所用的气象数据来源于中国气象数据网(http:∥data. cma. cn)，包括20个国家气象站点1961—2020年的逐日降雨数据，各站点的基本信息见表1。

表1 气象站基本信息

2 研究方法

2.1 降雨侵蚀力的计算

利用日降雨量数据计算降雨侵蚀力的模型有很多，其中章文波模型在国内被广泛应用[21-24]，本文主要利用该模型来计算降雨侵蚀力[25-26]，其计算公式为：

(1)

(2)

α=21.586γ-7.1891

(3)

式中:Ri为第i个半月的降雨侵蚀力[(MJ·mm)/(hm2·h·a)];n为半月内侵蚀性降雨量(≥12 mm)的天数(d);Pk为半月内第k天的侵蚀性降雨量(≥12 mm),mm;Pd12和Py12是日雨量≥12 mm的日平均降雨量和年平均降雨量(mm);α和γ均为模型的参数,利用日降雨资料进行估算,能够反映不同区域的降雨特征。年降雨侵蚀力主要通过累加半月降雨侵蚀力得到。

2.2 时空变化分析方法

本文利用泰森多边形法[27-28]计算出各气象站的权重，然后计算各气象站降雨侵蚀力等的加权平均值，并将其作为整个研究区的平均值；利用Mann-Kendall趋势分析法和气候倾向率法分析降雨侵蚀力等序列的时间变化趋势[29-30]；利用累积距平法和Mann-Kendall法共同进行数据序列的突变检验[31-32]；利用小波分析来分析降雨侵蚀力的周期变化，通过MATLAB计算小波系数，并通过Excel计算小波系数的实部，最后利用Surfer9绘制出小波等值线图进行分析，通过分析图中丰枯的交替来判断周期的变化规律[33]。

空间变化特征主要利用空间插值来进行分析，常见的空间插值法包括样条函数插值、克里金(Kriging)插值和反距离加权插值(IDW)等方法。本文主要是通过反距离加权插值方法来进行空间变化特征分析[34]。通过ArcGIS 10.4对沂蒙山区的年降雨量、年侵蚀性降雨量及年降雨侵蚀力进行空间插值，得到插值分布图进行空间分析。

3 结果与分析

3.1 降雨侵蚀力时间变化特征

3.1.1 年际变化 由于降雨侵蚀力与降雨量、侵蚀性降雨量密切相关，因此本文在分析降雨侵蚀力的同时，对研究区降雨量和侵蚀性降雨量进行相关分析。研究区1961—2020年年均降雨量、侵蚀性降雨量和降雨侵蚀力分别为785.88 mm，603.55 mm及5 081.59 (MJ·mm)/(hm2·h·a)，其中侵蚀性降雨量占年降雨量的76.80%。三者的年际变化特征见表2，其极值比分别为2.4，2.8，3.7，变异系数Cv分别为0.20，0.24，0.29，由此可以看出，降雨侵蚀力的年际波动程度大于降雨量和侵蚀性降雨量。

表2 沂蒙山区降雨量、侵蚀性降雨量及降雨侵蚀力年际变化特征

由图2分析结果可知，年降雨侵蚀力和年降雨量、年侵蚀性降雨量均存在幂函数关系，且相关系数均在0.9以上，年降雨侵蚀力与年侵蚀性降雨量更相关一些。

图2 年降雨侵蚀力与年降雨量及年侵蚀性降雨量的关系

采用气候倾向率法分析研究区各气象站点降雨侵蚀力、降雨量和侵蚀性降雨量可知(表3)，1961—2020年期间，大部分气象站的降雨侵蚀力和侵蚀性降雨表现为正气候倾向率，变化最大的均是台儿庄站和泗水站，每10 a增加271.92 (MJ·mm)/(hm2·h·a)，143.05 (MJ·mm)/(hm2·h·a)，28.95 mm，11.58 mm；降雨侵蚀力变化最小的是蒙阴站和薛城站，每10 a增加2.17 (MJ·mm)/(hm2·h·a)，17.98 (MJ·mm)/(hm2·h·a)。只有日照站、莒县站、莒南站和邹城站的降雨侵蚀力表现为负气候倾向率，变化最大的日照站，每10 a减少181.33 (MJ·mm)/(hm2·h·a)；莒县站、邹城站、莒南站、费县站、平邑站和曲阜站的侵蚀性降雨量表现为负气候倾向率，变化最大的是莒县站，每10 a减少12.70 mm。降雨量的变化趋势与降雨侵蚀力和侵蚀性降雨有所不同，大部分气象站点表现为负气候倾向率，但变化最大的站点与降雨侵蚀力相同。

采用Mann-Kendall趋势检验法对研究区各气象站的降雨侵蚀力、降雨量和侵蚀性降雨量进行分析，见(表3)。Z为Mann-Kendall趋势检验法中的检验数，当|Z|≥1.96时，可判断该序列呈现显著的变化趋势。由分析结果可知，大部分气象站的降雨侵蚀力和侵蚀性降雨量均表现为上升趋势，只有小部分气象站表现为下降趋势，而大部分气象站的降雨量表现为下降趋势，并且所有气象站的降雨侵蚀力、侵蚀性降雨和降雨量均未通过α=0.05水平的检验，即变化均不显著。综合气候倾向率法和Mann-Kendall趋势检验法的结果可知，两种方法的分析结果基本一致。

表3 沂蒙山区各气象站点降雨侵蚀力和降雨量趋势检验结果

由图3可以看出，整个研究区的年降雨侵蚀力与侵蚀性降雨的变化趋势相同，呈现波动上升的变化趋势，而降雨量则呈现波动下降的变化趋势。降雨量出现了负气候倾向率，表现为每10 a减少3.2 mm，降雨侵蚀力和侵蚀性降雨为正增长趋势，每10 a增加的31.7 (MJ·mm)/(hm2·h·a)，5.8 mm。研究区降雨侵蚀力的变化趋势总体表现为：1969年以前呈现下降趋势，之后上升至1974年，随后波动下降，1990年之后波动上升和下降交替进行。

图3 沂蒙山区降雨量、侵蚀性降雨量及降雨侵蚀力年际变化

由沂蒙山区降雨侵蚀力的小波分析结果可知(图4)，研究区多年平均降雨侵蚀力存在显著的周期变化，在整个时间序列上存在5～10 a，16～24 a两类尺度的周期变化，且在1961—2020年均稳定分布，具有全域性。计算降雨侵蚀力系列小波方差并绘制小波方差图，图中峰值记为降雨侵蚀力序列变化过程中存在的主周期。由图4可以看出，存在2个峰值，分别为7 a和22 a，第一峰值是22 a，说明年均降雨侵蚀力序列在22 a左右周期震荡最强，是降雨侵蚀力变化的第一主周期，7 a尺度为第二主周期。

由Mann-Kendall突变检验结果可知(图5)，降雨侵蚀力的正序列UF和逆序列UB相交于1963年、1965年、2002年、2007年，且均未超过0.05显著水平，均有可能为突变年份。为了进一步验证突变年份，结合累积距平法来具体分析。由图6可以看出，降雨侵蚀力在2002年明显由下降趋势转换为增长趋势，因此判断2002年为降雨侵蚀力的突变年份。

图5 年降雨侵蚀力Mann-Kendall突变分析 图6 年降雨侵蚀力累积距平曲线

3.1.2 年内及季节变化 降雨侵蚀力年内分布规律见图7，多集中在汛期6—9月份，占全年降雨侵蚀力的84.15%；降雨侵蚀力最大月份出现在7月，为1 706.44 (MJ·mm)/(hm2·h·a)，占全年的33.58%；最小月份出现在1月份，为14.35 (MJ·mm)/(hm2·h·a)，仅占全年的0.41%，月际差异明显。

图7 沂蒙山区降雨侵蚀力年内分布

春(3—5月)、夏(6—8月)、秋(9—11月)、冬(12-翌年2月)四季降雨侵蚀力的变化特征见图8。除秋季外，春季、夏季和冬季的降雨侵蚀力均呈现上升的变化趋势。春季和秋季的平均降雨侵蚀力相差不大，分别为512.96 (MJ·mm)/(hm2·h·a)，720.03 (MJ·mm)/(hm2·h·a)，占全年的降雨侵蚀力的10.09%和14.17%。夏季降雨侵蚀力最大，平均为3 787.59 (MJ·mm)/(hm2·h·a)，占全年的74.54%。冬季由于降雨较少，降雨侵蚀力最小，仅为61.02 (MJ·mm)/(hm2·h·a)，占全年的1%左右。

图8 沂蒙山区降雨侵蚀力四季变化

3.2 降雨侵蚀力空间变化特征

3.2.1 降雨侵蚀力空间分布特征 各气象站点年均降雨侵蚀力存在差异，变化范围为4 081.71～5 982.92 (MJ·mm)/(hm2·h·a)，降雨侵蚀力最大的是临沂站，最小的是沂源站，最大值约为最小值的1.5倍(图9)。

图9 沂蒙山区各气象站点多年平均降雨侵蚀力

为了进一步对比降雨侵蚀力与降雨量、侵蚀性降雨量空间分布的特征，利用IDW法进行空间插值，得到插值分布见图10所示。由结果可以看出，三者的空间分布特征基本上是一致的，均是由东南向西北呈带状逐渐递减。根据降雨侵蚀力由大到小可以将研究区大体分为3个区域：东南部地区包括日照、莒南、临沂、费县、枣庄、峄城及台儿庄，降雨侵蚀力基本达到5 200 (MJ·mm)/(hm2·h·a)以上；中西部地区包括微山、薛城、蒙阴、沂南、莒县及五莲等地，降雨侵蚀力4 900～5 200 (MJ·mm)/(hm2·h·a)；北部及西北地区包括沂源、沂水、滕州、平邑、曲阜、邹城和泗水，降雨侵蚀力在4 900 (MJ·mm)/(hm2·h·a)以下。

图10 沂蒙山区降雨侵蚀力、降雨量及侵蚀性降雨量空间分布

由于季节的变化对降雨侵蚀力也存在一定的影响，因此对沂蒙山区四季降雨侵蚀力进行空间分布特征分析，结果见图11。各季节降雨侵蚀力的空间变化整体也是由东南向西北呈带状逐渐递减，其中夏季降雨侵蚀力占全年的比例最大，各气象站点变化范围为3 144.03～4 327.68 (MJ·mm)/(hm2·h·a)，空间分布特征也与年均降雨侵蚀力的分布最为相似；秋季除日照附近的降雨侵蚀力相对较大之外，其余地区降雨侵蚀力的分布相对较为均匀，基本在600～800 (MJ·mm)/(hm2·h·a)。

图11 沂蒙山区季节降雨侵蚀力空间分布

3.2.2 降雨侵蚀力空间变化分析 利用变异系数Cv和Mann-Kendall统计量Z值来分析研究区降雨侵蚀力的空间变化(图12)。各气象站变异系数的范围是0.32～0.53，且地区差异比较明显，西部地区的变异系数相对较大，南部地区的变异系数相对较小，主要是由于南部地区的侵蚀性降雨量较大，并且年际变化相对稳定，因此变异系数基本小于0.40。由降雨侵蚀力的统计量Z值可知，除邹城、滕州、莒南、费县、莒县和日照站的降雨侵蚀力表现出降低的趋势外，其余站点均呈现增长的趋势。但无论是降低还是增长的变化趋势，统计量|Z|<1.96，均未通过显著性检验，均表现为不显著的变化趋势。

图12 沂蒙山区降雨侵蚀力变异系数及Z值空间分布

4 讨 论

计算降雨侵蚀力的模型有很多，其中利用降雨动能和降雨强度计算降雨侵蚀力时精度相对最高，但由于该方法对数据量、数据处理等要求极高，在很多地方均难以实现。由于本次研究数据的限制，缺乏长序列的次降雨数据，因此本文选取章文波的日降雨侵蚀力模型，得出沂蒙山区的年均降雨侵蚀力为5 081.59 (MJ·mm)/(hm2·h·a)，相对于年降雨侵蚀力及月降雨侵蚀力模型，该模型能够较为准确地反映研究区降雨侵蚀力的时空分布特征。但在今后的研究中，可以考虑通过加强对各类数据的收集与整理，尝试利用降雨强度、降雨动能等数据对研究区的降雨侵蚀力进行更深入的探讨。降雨侵蚀力与降雨量及侵蚀性降雨均呈现幂函数关系，其相关系数均大于0.9，且幂指数均大于1，由此说明降雨侵蚀力与降雨的关系非常密切。年内降雨侵蚀力集中分布在汛期6—9月份，尤其以7月、8月份分布最为集中，因此要重视对研究区7月、8月份的水土流失灾害防治。

研究区降雨侵蚀力由东南向西北呈带状逐渐递减，其中降雨侵蚀力最大的几个气象站是日照、莒南、临沂、费县、枣庄、峄城及台儿庄，根据张芷温[35]、杨韶洋[19]的研究，沂蒙山区强烈侵蚀主要集中分布在海拔相对较高、坡度较大以及植被覆盖差的山丘区，如沂源县、沂水县、蒙阴县、莒县、五莲县等。坡度大的地方地表径流的流速就越快，对土壤的冲刷侵蚀力就越强；有植被覆盖的地方可以较好地涵蓄水分，对土壤起到保护作用；同时人类也会通过破坏植被和稳定地形，增加水土流失。这也从一定程度上说明降雨侵蚀力的强弱除了与降雨有密切的关系之外，还与海拔、坡度、植被覆盖及人类活动等都有一定的关系。因此无论是沂蒙山区还是其他相类的地区，在防治水土流失过程中，要考虑各因素的影响，进行综合防治。

5 结 论

(1) 沂蒙山区年均降雨量、侵蚀性降雨量和降雨侵蚀力分别为785.88 mm，603.55 mm及5 081.59 (MJ·mm)/(hm2·h·a)，年降雨侵蚀力与年降雨量及侵蚀性降雨量存在幂函数关系，且相关系数均在0.9以上；年降雨侵蚀力与侵蚀性降雨的变化趋势相同，呈现波动上升的变化趋势，而降雨量则呈现波动下降的变化趋势，且降雨侵蚀力的年际波动程度大于降雨量和侵蚀性降雨量。各气象站点均表现为不显著的变化趋势，大部分气象站的降雨侵蚀力和侵蚀性降雨量表现为上升趋势，而大部分气象站的降雨量表现为下降趋势。

(2) 沂蒙山区降雨侵蚀力在整个时间序列上存在5～10 a，16～24 a两类尺度的周期变化，且在1961—2020年均稳定分布，具有全域性，其中22 a左右为第一主周期，7 a为第二主周期；降雨侵蚀力在2002年发生突变。

(3) 降雨侵蚀力年内分布多集中在汛期6—9月份，占全年的84.15%；除秋季外，春季、夏季和冬季的降雨侵蚀力均呈现上升的变化趋势。

(4) 降雨侵蚀力、降雨量与侵蚀性降雨的空间分布均呈现由东南向西北呈带状逐渐递减，夏季降雨侵蚀力的空间变化与年均降雨侵蚀力最为相似；各气象站变异系数的范围是0.32～0.53，且地区差异比较明显，西部地区的变异系数相对较大，南部地区的变异系数相对较小；除邹城、滕州、莒南、费县、莒县和日照站的降雨侵蚀力呈现降低趋势外，其余站点均呈现不显著增长的趋势。