降雨侵蚀力反映了降雨对土壤侵蚀的潜在能力, 是通用土壤流失方程USLE^[1]和修正通用土壤流失方程^[2]的主要影响因素之一。因此, 分析区域尺度降雨侵蚀力时空分布特征, 对流域土壤侵蚀监测、评估、预报, 制定合理的防治措施, 农业耕种和开发建设项目进度安排都具有重要意义。自次降雨总动能E与30 min最大雨强I₃₀的乘积EI₃₀作为降雨侵蚀能力指标并被应用于USLE以来, 有关降雨侵蚀力R值的研究越来越被人们所重视^[3]。EI_7.5、EI₁₀、EI₆₀等计算R值的指标^[4-5]都基于小时或分钟, 雨量信息资料难以获取。为了克服这些限制, 学者们开始应用基于年、月、日尺度的R值计算方法^[6-7]。与基于年、月、日尺度的R值计算方法相比, 日降雨资料更易获得, 利用日降雨量建立的降雨侵蚀力估算模型计算的R值精度更高, 含有更多降雨特征信息并综合反映了降雨侵蚀力的特点^[8]。国内计算R值一般以月为尺度, 但由于地表植被覆盖在1个月内可能变化极大, 故以半月为尺度分析降雨侵蚀力的空间分布^[9]。章文波等^[10]提出的日降雨侵蚀力估算方法使用最为广泛^[11-13]。

降雨侵蚀力具有年内分布不均、年际年代变化明显、时空分布与降雨量基本一致等特征。钟科元等^[14]分析了松花江流域降雨侵蚀力的时空分布, 显示降雨侵蚀力的空间分布与降雨量基本一致, 自东南向西北递减并且年际变化趋势呈现明显的区域分异。邱美娟等^[15]同样证明了吉林省降雨侵蚀力的空间分布与降雨量基本一致, 时间分布上与多年平均降水量的时间分布特征具有高度一致性。戴金梅等^[16]研究结果表明多年平均降雨侵蚀力与降雨量呈极显著正相关, 空间上西高东低, 与降雨量分布基本一致, 研究期内降雨侵蚀力呈不明显上升趋势。

连云港赣榆区低山丘陵较多, 土层浅薄, 属于北方土石山区, 水力侵蚀是该区土壤侵蚀主要类型, 开展降雨侵蚀力相关研究十分必要。已有研究大多集中在全国、大流域和省区, 小区域降雨侵蚀力研究很少^[16], 赣榆区缺乏降雨侵蚀力相关研究。笔者基于连云港赣榆区12个雨量站1954—2014年日降雨资料, 通过章文波等提出的日降雨侵蚀力估算方法^[17]计算降雨侵蚀力, 探讨其时空变化规律和特点，以期为赣榆区的水土流失治理和水土保持规划提供科学依据, 为北方土石山区的水土流失预测、水土保持、生态环境保护以及灾害防治等提供参考。

1 研究区概况

赣榆区位于连云港西北部(图 1), 东滨黄海, 地理位置为北纬34°41'~35°07', 东经118°45'~119°18'。该区年均气温22.8 ℃, 属暖温带海洋性季风气候区, 四季分明, 气候温和。地势由西北向东南倾斜^[18], 总面积1 514 km², 低山丘陵区约占总面积的29.2%, 平原区约占61.8%, 其余区域为沿海滩涂。自然植被以马尾松(Pinus massoniana)、侧柏(Platycladus orientalis)、阔叶树、灌木丛及荒草等次生植被为主, 但多遭到破坏, 大多为人工植被。

2 数据来源与分析方法 2.1 数据来源 2.1.1 赣榆区日降雨资料

收集连云港赣榆区12个雨量站1954—2014年的日降雨资料, 依次是夹谷山、岳庄、石桥、朱汪、清水涧、黑林、汪子头、小塔山水库、青口、班庄、朱堵集、范河闸, 表 1介绍了12个站点的详细情况。

2.1.2 径流小区监测资料

夹谷山坡面径流场位于水土保持重点防护区, 是连云港市的主要监测站点。该坡面径流小区主要土壤为花岗片麻岩发育的砂壤土。径流小区平行布设在坡度为8°的同一坡面上, 共设立6个5 m×20 m小区, 分别为乔木小区、经济林小区、牧草小区、顺坡耕作区、标准小区和梯田小区。选取2015年夹谷山坡面径流场基于30 min雨强I₃₀计算得到的降雨侵蚀力资料。

2.2 研究方法 2.2.1 降雨侵蚀力计算方法

章文波等^[10]修正的Richardson日降雨侵蚀力估算模型为

$ {R_{半月}} = \alpha \sum\limits_{k = 1}^m {P_k^\beta } , $

(1)

$ \alpha = 21.586{\beta ^{ - 7.1891}}, $

(2)

$ \beta = 0.8363 + \frac{{18.144}}{{{P_{{\rm{d12}}}}}} + \frac{{24.455}}{{{P_{y{\rm{12}}}}}}。$

(3)

式(1)~(3)中, k=1, 2, 3, …, m, 为某半月内降雨时间, d; p_k为半月内第k天的日雨量(要求日雨量≥12 mm, 否则以0计), mm; α和β为模型参数; R_半月为某半月内的降雨侵蚀力, MJ·mm·hm^-2·h^-1·a^-1; P_d12为日雨量≥12 mm的日平均雨量, 即1 a中≥12 mm日雨量总和与相应时间的比值, mm; P_y12为日雨量≥12 mm的年平均雨量, 即1 a中≥12 mm日雨量年累加值的多年平均值, mm。利用式(1)~(3)计算每年各半月的降雨侵蚀力, 经过汇总得月降雨侵蚀力、年降雨侵蚀力、多年平均降雨侵蚀力。该文中侵蚀性降雨量即为日雨量≥12 mm的日平均雨量。

2.2.2 Mann-Kendall检验

Mann-Kendall检验是一种基于秩的非参数统计检验方法, 统计值S和标准化检验统计值Z的计算公式为

$ S = \sum\limits_{i = 1}^{n - 1} {\sum\limits_{j = i + 1}^n {{\mathop{\rm sgn}} \left( {{x_j} - {x_i}} \right)} } , $

(4)

$ {\mathop{\rm sgn}} \left( {{x_j} - {x_i}} \right) = \left\{ \begin{array}{l} 1,{x_j} > {x_i}\\ 0,{x_j} = {x_i}\\ - 1,{x_j} < {x_i} \end{array} \right., $

(5)

$ E\left( S \right) = 0, $

(6)

$ \begin{array}{l} V\left( S \right) = \left[ {n\left( {n - 1} \right)\left( {2n + 5} \right) - \sum\limits_{i = 1}^n {{t_i}\left( {i - 1} \right)} } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\left. {\left( {2t + 5} \right)} \right]/18, \end{array} $

(7)

$ Z = \left\{ \begin{array}{l} \frac{{S - 1}}{{\sqrt {{\rm{Var}}\left( S \right)} }},S > 0\\ 0,S = 0\\ \frac{{S + 1}}{{\sqrt {{\rm{Var}}\left( S \right)} }},S < 0 \end{array} \right.。$

(8)

式(4)~(8)中, sgn (x_j-x_i)为符号参数; x_i和x_j为统计时间序列中年份i和j所对应的数值; n为时间序列数据长度; t_i为第i个数对应的捆绑值。Mann-Kendall突变检验中, 正向序列曲线UF的变化趋势显示降雨侵蚀力的变化情况, ＞0表示呈上升趋势, ＜0表示呈下降趋势, 当UF值超过临界线则表示上升或下降趋势显著。Mann-Kendall统计量Z值可表示时间序列数据的变化趋势。Z为正, 表示呈增加趋势; Z为负, 表示呈下降趋势; Z的绝对值越大, 说明该序列的变化趋势越明显, 当Z的绝对值＞1.96时, 通过0.05置信水平的显著性检验^[19]。

2.2.3 重标极差法

赫斯特指数(H)基于重标极差法(rescaled range analysis, 简称R/S), 用于定量分析时间序列^[20]。时间序列{N_t}包含n个日观值, 极差R、标准差S、重标极差R/S和赫斯特指数H的计算公式为

$ {X_{k,n}} = \sum\limits_{t = 1}^k {\left( {{N_t} - {M_n}} \right)} ,k = 1,2, \cdots ,n, $

(9)

$ R = \max \left( {{X_{k,n}}} \right) - \min \left( {{X_{k,n}}} \right), $

(10)

$ S = \sqrt {1/n\sum\limits_{t = 1}^n {{{\left( {{N_t} - {M_n}} \right)}^2}} } , $

(11)

$ R/S = C{n^H}, $

(12)

$ \lg \left( {R/S} \right) = \lg C + H\lg n。$

(13)

式(9)~(13)中, X_k，n为n个时间序列期间的积累极差; N_t为第t日的观测值; M_n为n个时间序列期间N_t的算术平均值; C为常数。绘制lg (R/S)和lg n趋势图, 其回归模型斜率即为H的近似值。采用赫斯特指数H研究整个区域降雨侵蚀力R值的稳定度和发展趋势。当H＞0.5, 说明降雨侵蚀力的未来变化表现出持续性, 与过去趋势一致; 当H＜0.5, 说明降雨侵蚀力的未来变化与过去趋势表现出负相关性, 即与过去趋势相反。

2.2.4 变异系数

变异系数是原始数据标准差与原始数据平均数的比值。变异系数反映降雨侵蚀力的离散程度, 以0.5和1为界, 分为低等程度变异、中等程度变异和高等程度变异。

3 结果与讨论 3.1 日降雨侵蚀力模型可行性分析

以夹谷山径流小区基于I₃₀计算得到的降雨侵蚀力R₁为降雨侵蚀力真实值, 对基于章文波日降雨量模型计算的降雨侵蚀力R₂进行讨论并验证。表 2为赣榆区夹谷山径流小区2015年6—11月的降雨侵蚀力。

7月是R₁的高值期; 8和11月是中值期; 9月R₁占总降雨侵蚀力的比例最小。由于11月降雨多且雨强小, 导致11月R₁与R₂相比偏小。7月, R₁占总降雨侵蚀力的比例也远高于R₂, 这是由于章文波模型估算的半月侵蚀力对短历时高强度降雨可能会有一定程度的低估。因此夹谷山径流小区基于I₃₀计算得到的R₁对降雨极端情况反映灵敏; 基于日降雨量的章文波模型计算的降雨侵蚀力与真实值相比偏大, 且会低估高强度降雨的降雨侵蚀力。

虽然章文波日降雨量模型对集中且高强度降雨的估算存在偏差, 但连云港赣榆区为暖温带海洋性季风气候, 降雨均匀, 降雨强度相对平缓, 以蓄满产流为主, 同时考虑到资料可得性和应用的方便性, 基于日降雨量的章文波模型计算的降雨侵蚀力能反映当地的降雨及强度特征和长期趋势。

3.2 连云港赣榆区降雨侵蚀力的时间变化 3.2.1 年内变化

1954—2014年间, 连云港赣榆区降雨侵蚀力、降雨量和侵蚀性降雨量年内分布曲线一致且均为单峰型。三者年内分布主要集中于5—9月, 7月达最大, 其次是8月(图 2)。

由表 3可知，7月、7—8月及5—9月这3个时段中, 降雨侵蚀力占全年的比例均为最高, 其次为侵蚀性降雨量, 总降雨量占比最小。

结合赣榆区年内降雨侵蚀力分布特点, 将其分为4个时期:高值期、中值期、低值期和无侵蚀期。降雨侵蚀力高值期为7—8月, 7月降雨侵蚀力最大, 占年降雨侵蚀力的35.2%, 8月占比为26.3%。高值期7—8月侵蚀性降雨量占总降雨量的比例均为90%。降雨侵蚀力中值期为6和9月, 这2个月多年平均降雨侵蚀力占年降雨侵蚀力的比例差异性不大。中值期降雨量大概仅为8月的一半, 但月侵蚀性降雨量占总降雨量的比例均为80%。低值期为5月, 该月降雨侵蚀力占年降雨侵蚀力的6.5%。无侵蚀期为10月至翌年4月。因此, 连云港赣榆区土壤侵蚀预防工作的重点应该放在4—9月, 尤其是7—8月。

3.2.2 年际变化

1954—2014年, 赣榆区平均R值为6 154.93 MJ·mm·hm^-2·h^-1·a^-1, 1974年达最大值12 369.16 MJ·mm·hm^-2·h^-1·a^-1; 最小值出现在1966年, 为2 048.34 MJ·mm·hm^-2·h^-1·a^-1。降雨量均值为863.04 mm, 1974年出现最大值1 322.16 mm, 1954年出现最小值367.45 mm。研究区域内侵蚀性降雨量均值为684.41 mm, 同样在1974年出现最大值1 114.49 mm, 最小值302.65 mm也出现在1954年。上述数据表明, 连云港赣榆区降雨侵蚀力R值的年际变化比降雨量与侵蚀性降雨量的年际变化更为剧烈, 更为敏感。史东梅等^[21]也证实了这一结论。

由图 3可见，连云港赣榆区降雨侵蚀力年际变化为波浪形曲线, 除1961、1970、1974和1990年降雨侵蚀力较高外, 1954—2000年降雨侵蚀力波动较为平缓, 而21世纪极端降雨事件增多, 区域内年均降雨侵蚀力波动也明显变大。为分析降雨侵蚀力年际序列的突变情况, 采用Mann-Kendall突变检验方法对全区域61 a的降雨侵蚀力序列进行突变分析(图 4), 发现1954—1966年连云港赣榆区降雨侵蚀力呈显著上升趋势, 并在1957和1961年超过0.05置信水平, 但其余年间均未发生显著突变。1966—1970年间降雨侵蚀力呈下降趋势; 降雨侵蚀力在1970—1976年(除1973年外)开始上升; 1977—2006年, 降雨侵蚀力开始下降; 降雨侵蚀力在2007—2014年发生转折突变, 呈上升趋势。

3.2.3 年代变化

以1954—1959年为基准数据计算变化率, 从表 4可知, 1954—2014年降雨量、侵蚀性降雨量以及降雨侵蚀力R值的变化率基本为正值, 绝大部分(除20世纪80年代的侵蚀性降雨量、降雨侵蚀力外)相对于基准期增加。

其中2000—2009年降雨量、侵蚀性降雨量以及降雨侵蚀力R值的变化率均为最大值。降雨侵蚀力在20世纪50—70年代逐渐增加, 80年代减小为负值, 90年代又开始上升, 21世纪初呈减小趋势, 但不低于基准期。降雨侵蚀力在20世纪80年代数值为负, 这是由于侵蚀性降雨量减少, 非侵蚀性降雨量较多。这表明进入21世纪以来, 随着降雨侵蚀力增加, 水土流失可能会加剧。

3.3 连云港赣榆区降雨侵蚀力空间变化 3.3.1 空间分布

采用ArcGIS地统计模块的反距离权重插值法, 分析多年平均降雨量以及多年平均降雨侵蚀力的空间分布及变化趋势(图 5~6)。由图 6可知, 以岳庄—青口一线为界, 此线的西北部年均降雨侵蚀力呈现自东南向西北逐渐递减的趋势, 东南一侧年均降雨侵蚀力则由东南向西北逐渐递增。除岳庄—青口一线东南部外, 赣榆区降雨侵蚀力R值分布特征与降雨量(图 5)基本相似, 整体趋势为自东南向西北逐渐减少。这与刘斌涛等^[22]的研究结果一致。赣榆区受台风影响, 随台风势力的减弱, 降雨侵蚀力由沿海向内陆减少。在岳庄—青口一线东南侧, 多年平均降雨量的低值位于范河闸, 而降雨侵蚀力低值位于清水涧, 这说明降雨侵蚀力不仅受降雨量影响, 还与降雨历时、降雨强度等相关。

3.3.2 空间变化

1954—2014年, 连云港赣榆区各站点降雨侵蚀力R均值分布范围为5 358.1~6 792.6 MJ·mm·hm^-2·h^-1·a^-1, 在青口站最大, 清水涧站最小; 各站点多年降雨侵蚀力R值的最大值分布在11 059.0~17 443.6 MJ·mm·hm^-2·h^-1·a^-1之间(表 5), 2000年范河闸站达最大值; 各站点多年降雨侵蚀力R值的最小值分布在101.8~2 194.4 MJ·mm·hm^-2·h^-1·a^-1之间, 1976年清水涧站达最小值。

由图 7可见，各站点R值变异系数在0.41~0.52之间, 最大变异系数出现在范河闸站, 其次是朱汪站, 都属于中等程度变异, 其余站点R值变异系数均在0.41~0.48之间, 属于低等程度变异。沿海区域变异系数比内陆区域年际波动大。除夹谷山、石桥和黑林3个站点外, 连云港赣榆区其他站点Z值都＞0, 呈现上升趋势, 岳庄、清水涧和汪子头均表现为明显上升趋势(图 8)。意味着赣榆区未来水土保持压力可能会增大, 尤其是研究区南部的岳庄、清水涧和汪子头。赫斯特指数H值在0.43~0.76之间, 只有范河闸的H值＜0.5, 这意味着1954—2014年范河闸整体R值呈现不明显上升趋势, 但2014年之后降雨侵蚀力会开始下降。其他站点H值都＞0.5, 则会延续之前的趋势发展。

连云港赣榆区降雨侵蚀力R值变异系数为0.39, Z值为0.55, H值为0.6。总而言之, 整个区域呈低等程度变异且呈不明显的上升趋势。2014年之后, 连云港赣榆区将继续之前的上升趋势, 即由降雨侵蚀力引起的土壤流失风险增加。

4 结论

根据北方土石山区连云港赣榆区1954—2014年逐日降雨资料, 对区域年内、年际、年代及空间的降雨侵蚀力进行统计分析, 得出以下结论:

(1) 时间尺度上, 研究区降雨侵蚀力R值多年平均降雨侵蚀力为6 154.9 MJ·mm·hm^-2·h^-1·a^-1, 与降雨量和侵蚀性降雨量年内变化基本一致, 都为单峰型。该区89%的降雨侵蚀力集中于4—9月。R值的年际变化和年代变化均显示, 由于极端降雨的影响, 21世纪以后R值变化幅度会增大。

(2) 空间尺度上, 岳庄—青口一线西北部R值整体分布特征与降雨量基本一致, 呈现出自东南向西北逐渐递减的趋势; 东南一侧R值则自东南向西北逐渐递增; 沿海区域R值年际波动比内陆区域大, 且变化趋势具有明显的区域差异。

(3) 连云港赣榆区12个雨量站中, 夹谷山、石桥和黑林3个站点R值呈下降趋势, 占25%。其他站点R值均呈上升趋势, 其中岳庄、清水涧和汪子头表现为明显上升趋势; 整个区域中只有范河闸的赫斯特指数H值＜0.5, 表明2014年之后开始下降, 其他站点R值则会延续之前的发展趋势。

(4) 整个区域降雨侵蚀力R值变异系数为0.39, Z值为0.55, H值为0.6。该区R值呈低等程度变异且呈不明显的上升趋势, 在2014年之后, R值将继续呈上升趋势, 区域面临一定的水土流失潜在压力。

从连云港赣榆全区来看, 受极端气候环境的影响, 21世纪以来年降雨侵蚀力波动增大且呈现不明显上升趋势。未来应注重研究区4—9月水土流失工作的落实与开展, 特别是R值较大的岳庄—青口地区和R值增幅较大的南部及沿海区域。