1970-2019年大汶河流域降雨侵蚀力时空变化分析

许功伟, 徐立荣, 张云苹, 王 坤, 惠 莉, 徐 晶, 赵芹蕊

(1.济南大学 水利与环境学院, 济南 250022; 2.禹城市水利局, 山东 德州 251200; 3.济南黄河河务局 章丘黄河河务局, 济南 250200)

在世界范围内，土壤侵蚀问题早已备受关注，我国土壤侵蚀问题也尤为突出[1]。据研究表明，降雨会破坏土壤的结构，其形成的汇流会带走土壤颗粒，这就直接造成了土壤的流失，而雨滴大小分布和雨滴终端速度决定的降雨侵蚀力能反映出降雨造成土壤侵蚀的潜力[2-4]。降雨侵蚀力的计算方法从1958年到现在已经提出多种。最初是1958年Wischmeier等[5]提出的利用次降雨动能和降雨强度的乘积来计算，后来又由科学家提出用次降雨量和降雨强度乘积来计算[6]，尽管此方法的计算避免了降雨动能的计算，但同样需要详细的次降雨量数据，并且计算的工作量繁重。再后来，科学家们提出了简易估算模型，从最初孙保平等[7]建立并在西吉县黄土丘陵区试验的年雨量模型到后来的黄炎和等[8]提出的在闽东南地区试验的月雨量模型，再到现在使用广泛的章文波等[9]提出的通过日雨量计算降雨侵蚀力的日雨量模型，这就使得降雨侵蚀力的计算量得到显著的下降，也使得精度逐渐的提高。

大汶河流域地势东西差异大，从东至西落差超300 m。大汶河流域受季风影响，夏季暴雨频繁，其他三季的降雨相对较少。受地势和气候的双重影响，大汶河流域自建国以来就洪涝不断，土壤侵蚀问题也比较严重，不仅影响日常生活，也阻碍经济发展。这些现象表明，想要让大汶河流域突破种种发展阻碍，首先要解决其土壤侵蚀问题，从大汶河流域的背景来看，要解决土壤侵蚀问题，对大汶河流域的降雨侵蚀力进行分析就显得尤为重要。

因此，本文基于1970—2019年大汶河流域6个市级雨量站的日雨量数据和日降雨侵蚀力模型来计算降雨侵蚀力，并借助ArcGIS，MATLAB等工具分析其时间和空间的变化，结合前人对大汶河流域的研究，寻求大汶河流域中降雨侵蚀力的演化规律，共同为缓解大汶河流域的土壤侵蚀问题提供科学数据依据，以此来有效地缓解水土侵蚀对大汶河流域的影响。

1 研究区概况

大汶河(35°45′—36°29′N，116°16′—117°53′E)位于山东省中部，全长209 km，流域面积8 633 km2。大汶河流域属大陆性季风气候，旱涝季分明，夏季降雨多，而春、秋、冬季则降雨较少。大汶河流域地形地貌主要分为平原、丘陵、山地，还有少部分的洼地、水面。北部及东部的上游地区是山地，而其下游地区的海拔较低，从东北至西南方向地势逐渐降低，这种地势分布就造成了大汶河流域河流呈自东向西走向的走势。其植被覆盖，受暖温带影响，山区植被覆盖相对较少，以松、柏等植物为主，中游等丘陵地区则以灌木居多，下游主要种植梧桐、杨柳等植物。由于这种特殊的气候和地形以及植被，就使得此流域洪涝灾害、水土流失等现象极易发生。

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为正时，表现为上升的趋势，当超过置信区间±1.96(α=0.05)时，则为显著上升，反之同理。此法还可用来判断系列突变开始的时刻，主要通过两个变量来判断，按时间顺序统计的UF和其逆序列UB，两者的交点，就是突变发生的时刻[13-14]。

(2) 累积距平法。很多时候单用Mann-kendall法分析突变点时，图中会出现多个交点，无法判断，这时就需要通过累积距平法与其共同分析出突变点。累积距平法的原理是通过求出此序列的平均值，然后用逐年的数据与平均值做差，得出的便是距平值，然后逐年累加，就得出累积距平序列[15-17]。公式如下：

(4)

(3) 小波分析。分析时间序列的问题时，基础的时域与频域的分析往往不能满足一些受多种因素影响的现象，因此就需要通过小波分析来进行分析，结合时域与频域，在时间的变化上分析其周期的改变性来满足分析的需要。其分析流程是通过MATLAB来计算小波系数，然后通过Excel计算小波系数的实部，最后通过Furfer 9绘制出小波等值线图进行分析。通过分析图中枯丰的交替来判断周期的变化规律[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年后，继续波动下降至2019年。年降雨侵蚀力的变化与其大体一致。图2A中，降雨量的5 a移动平均呈微弱上升趋势，但在图2B中，降雨侵蚀力的趋势却是存在着下降的状态，表明虽然降雨量在逐渐增加，但其侵蚀性雨量却存在一个逐步下降的趋向。由于人口增长及人类活动的影响，全球气温正在不断变暖，尤其近10 a来增加更为剧烈，全球气温变暖使得水循环加剧，大气中水分含量增加，进而积雨云出现的概率增加[20]，大汶河流域东、北部是山地，研究站点大部分位于雨影区与平原区，就使得研究区容易出现强对流天气，进而降雨量增加。降雨侵蚀力受侵蚀性降雨影响，大汶河流域侵蚀性降雨大多集中在夏、秋两季。目前，全球气温变暖，ENSO事件不断发生，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这两类周期在全局都比较稳定，而4～7 a这类的周期只在1988—2008年的变化较为稳定。

图3 基于小波分析的年降雨侵蚀力的周期变化

图4和图5为突变分析的结果，在图4中，1976年之前UF>0，1976—1978年UF<0，1978—1980年UF>0，1979年之后UF<0，说明1976年之前都呈上升趋势，1976—1980年发生短暂的先下降后上升趋势，1980年直至最后都呈下降趋势。UF与UB相交于1976年和1979年，且两个交点都位于±1.96(α=0.05)内，为了准确确定突变年份，结合累积距平法来具体分析。发现，图5曲线的变化趋势与图4基本相同，且在1976年曲线没有发生显著的突变变化，而在1979年以后曲线明显呈下降趋势，这与图4中1979年之后曲线明显下降，并且在1984年以后超出置信区间(变化显著)的分析结果恰好相同。通过上述两种方法共同证明了突变发生在1979年。

图4 年降雨侵蚀力M-K突变分析

图5 年降雨侵蚀力累积距平曲线

3.1.2 年内变化 从图6中可以发现降雨侵蚀力和降雨量的年内分布大致相同，主要分布在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 季节变化 图7为大汶河流域降雨侵蚀力四季变化图，据分析，春季到冬季的变化趋势为，从显著上升到微弱下降再到显著下降最后呈显著上升趋势，除了夏季变化不显著，春、秋、冬季均为显著变化，夏、秋两季呈现出减少趋势(这与3.1.1年际变化中的分析结果相同)，主要是因为全球气温逐渐升高，使得ENSO事件发生的频率增加，而经冯若昂等[21]发现ENSO暖事件会对会山东省内的降雨侵蚀力造成减少的影响。春季到冬季的平均R值依次为398.13，2 744.96，635.98，29.75 MJ·mm/(hm2·h·a)，各占全年降雨侵蚀力的10.45%，72.07%，16.70%，0.78%。夏季明显大于其他三季，而冬季明显低于其他三季，这与夏季降雨多、冬季降雨少有直接关系。通过对季节的变化的分析，可以看出夏季降雨侵蚀力较为集中，还可以看出夏季降雨侵蚀力的变化相较于其他季节较为稳定，集中且稳定的变化就可针对性的在每年的7月、8月加大水土保持工作的力度，这就大大减少了水土保持工作的强度。

图7 大汶河流域降雨侵蚀力四季变化

3.2 降雨侵蚀力空间变化分析

3.2.1 降雨量和降雨侵蚀力的空间分布分析 由上述时间分析发现，降雨量的改变确实在某种程度上影响着降雨侵蚀力，因此在空间上分析降雨侵蚀力，也首先对降雨量进行分析。通过图8A发现泰山周围的降雨量最高，基本都在900 mm之上，其次是莱芜、新泰，二者降雨量在730 mm以上，再次是肥城、泰安，也都超越650 mm，最低是东平，不足650 mm。大汶河流域降雨整体分布为中部最多，明显高于东、西部，再者是东部，西部最低。通过图8B可以看出，降雨侵蚀力的分布与前面分析的降雨量的分布有着很高的相似度，R值方面，中部地区最高，其值可达4 000 MJ· mm/(hm2·h·a)，以泰山为中心的地区甚至达到5 000 MJ·mm/(hm2·h·a)，东部整体高于西部。这种现象与大汶河流域的地势分布有关，其北部及东部是山地，地势较高，这就使得气流抬升，然后形成降雨，造成降雨较多，从而产生相对较大的降雨侵蚀力。因此在大汶河流域的水土保持工作中，强化以泰山为中心的山地地区的水土保护措施，会使其工作事半功倍。

图8 大汶河流域降雨量和降雨侵蚀力的空间分布

3.2.2 年降雨侵蚀力空间变化分析 通过对Z值和CV值的空间插值图(图9)进行分析，可以看出大汶河流域大多数雨量站的Z值都小于0，说明分析区站点的降雨侵蚀力都存在下降的趋向，尽管东平、莱芜等东、西部地区的下降趋势比泰山、泰安等中部地区的下降趋势大一些，但都没超过±1.96，说明均为不显著下降，肥城的Z值恰好为0，这就表明肥城地区的降雨侵蚀力在近50 a内的变化都较为稳定。变异系数的分布与Z值分布相似，都是东平最大，其次是新泰、泰安和莱芜，最低为肥城和泰山，整体上看中部较小，东、西部较大。说明降雨侵蚀力在东平、莱芜等东、西部地区不仅下降趋势较大，而且离散程度也较大。

图9 大汶河流域降雨侵蚀力Z值和变异系数的空间分布

3.3 不同地形区降雨侵蚀力分析

根据大汶河流域地形特点将大汶河流域划分为东、北部山地地区和中、西部平原地区两种地形区(以200 m海拔为界)，进而分析降雨侵蚀力与Z值、变异系数等的关系(表1)。

通过表1可以看出，东部、北部山地地区的降雨量和降雨侵蚀力明显高于中、西部平原地区，这是因为从海洋吹向陆地的暖湿气团遇到较高地势会造成气流抬升，从而使降雨增多，侵蚀性降雨也随之增加，进而表现出较大的降雨侵蚀力。还可以看出东、北部山地地区的|Z|比中、西部平原地区大，说明其下降趋势比中、西部平原地区明显。两类地形区和整体流域的变异系数依次为0.33，0.36，0.34，尽管中、西部平原地区的变异系数最大，但三者整体相比，差异不大，离散程度相当。在相关性方面，降雨侵蚀力与降雨量，除了中、西部为0.64，其余均超0.90，呈显著正相关，其与海拔的相关性同样如此，且都超过0.90，呈强相关性。

表1 大汶河流域不同地形区降雨侵蚀力分析

4 讨 论

大汶河流域年际降雨量与降雨侵蚀力分别表现出了微弱上涨和不显著下降的趋势。降雨侵蚀力的变化与冯若昂等[21]在分析山东省降雨侵蚀力的变化中提到以泰山为中心一带的年降雨侵蚀力的变化趋势以及肖蓓等[22]提到的鲁中南地区降雨侵蚀力变化趋势相符；降雨量的变化趋势与朱祎真等[23]分析的相同，也与徐晶等[13]分析其邻近济南市的变化相一致，但与夏继勇等[24]分析的大汶河流域近58 a降雨量的变化趋势不同，主要是由研究中采取不同长度的降雨序列(朱祎真等研究序列区间为1971—2010年，徐晶等研究序列区间为1970—2016年，夏继勇等研究序列区间为1956—2013年，本研究序列区间为1970—2019年)所造成的。

ENSO冷事件，会增加山东省夏秋季降水，从而使降雨侵蚀力增加，ENSO暖事件则相反[21]。ENSO冷事件一般伴随在ENSO暖事件之后发生，但由于全球气温正在不断变暖，ENSO暖事件的发生次数远多于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年发生突变。

(2) 按月分析，降雨侵蚀力和降雨量在7月、8月的分布尤为突出。一年中，7月、8月份的降雨侵蚀力和降雨量占比高达60.08%，50.87%。四季的年际趋势变化为，春季显著上升，夏季微弱下降，秋季显著下降，冬季显著上升，除了夏季变化不显著，春、秋、冬季均为显著变化。7月、8月两月的降雨侵蚀力，不仅集中而且稳定，加大这两个月的防控力度是整个水土保持工作的重中之重。

(3) 大汶河流域降雨整体分布为中部最多，明显高于东、西部，再者是东部，西部最低，降雨侵蚀力的分布与之类似；大汶河流域整体雨量站的Z值都小于0，说明分析区内站点的降雨侵蚀力都存在下降的趋向，尽管东平、莱芜等东、西部地区的下降趋势比泰山、泰安等中部地区的下降趋势大一些，但都没超过±1.96，说明均为不显著下降；大汶河流域离散程度的分布与Z值类似，东、西部离散程度较大，中部较小。

(4) 东部、北部山地地区的降雨量和降雨侵蚀力高于中、西部平原地区，且东、北部降雨侵蚀力的下降趋势也比中、西部明显，但离散程度相当，差异不大；降雨侵蚀力与降雨量、海拔在不同地区都有较强的关联性，这表明在东、北部地势较高的地区要注重水土流失的防控，也再次验证了在7月、8月雨季集中的时候要加大水土保持的力度的结论。