钱 峰，董林垚，黄介生，刘洪鹄，韩 培，孙 蓓



基于Hurst指数与相关系数的降雨侵蚀力变异识别与分级方法

钱 峰1,2，董林垚2，黄介生1※，刘洪鹄2，韩 培2，孙 蓓2

（1. 武汉大学水资源与水电工程科学国家重点实验室，武汉 430072；2. 长江科学院水土保持研究所，武汉 430010）

变化环境下区域降雨侵蚀力的时空变化问题对区域水土流失防治工作提出了新的挑战。降雨侵蚀力序列不再是纯随机序列，往往存在趋势、跳跃或者周期的变化，在对降雨侵蚀力序列分析与计算时，现有研究往往采用单一的检验方法，缺乏对降雨侵蚀力序列各类成分的综合比较，所得到的结果可信度及其程度如何无法判断。该研究提出了基于Hurst系数和相关系数的降雨侵蚀力序列联合分析方法。该方法首先计算降雨侵蚀力序列的Hurst系数，引用水文序列变异的概念，从统计学角度将降雨侵蚀力序列确定性成分分为三级（无变异、弱变异和强变异）。然后通过多种检验方法综合检验，将得到的结果与原序列进行相关性分析提取相关系数最大的确定性成分（趋势、跳跃和周期），对其进行剔除，重复上述步骤，将降雨侵蚀力序列中的确定性成分进行一一分解，最终得出的降雨侵蚀力序列将是一个随机序列与确定性序列的组合。实际应用中，根据长江流域174个气象站点1961—2014年逐日降雨资料，对流域内各气象站点年降雨侵蚀力序列进行确定性成分分析与分级结果表明：长江流域174个气象站点中有130个站点降雨侵蚀力序列无明显变异，有31个站点降雨侵蚀力序列出现弱变异，有13个站点降雨侵蚀力序列出现强变异。以重庆奉节站为例进行综合检验，分析结果为整体强变异，该站年降雨侵蚀力序列存在复合周期和跳跃成分，其中复合周期为5 a和16 a，向下的跳跃点为2011年。该研究为变化环境下区域降雨侵蚀力预测提供理论依据。

侵蚀；分级；Hurst指数；相关系数；降雨侵蚀力；时空变异

0 引 言

降雨侵蚀力反映了降雨引起土壤侵蚀的潜在能力，它受降雨量的控制。在全球气候变化的背景下，局部地区干旱持续范围增加，暴雨及降雨程度增加，从而引起区域侵蚀性降雨量的变化。学者们在世界各地针对变化环境下的降雨侵蚀力序列展开了不同尺度的研究，例如：国家尺度上有美国、澳大利亚、中国、瑞士和希腊等[1-5]；区域尺度上有德国Rhine-Westphalia北部地区、中国珠江流域、中国干旱地区和陕南地区等[6-9]。区域年降雨侵蚀力序列与降雨量序列一样可以分解为两种成分：确定性成分和随机性成分。确定性成分具有一定的物理意义，包含周期的和非周期的成分；随机性成分则由不规则的振荡和随机影响组成[10]。

降雨侵蚀力序列中的趋势是指序列随着时间发展变化的动向。序列中的各值平均来说，或是增加或是减少，而序列的统计参数则随着时间的推移呈现出系统且连续的增加或减少的变化。降雨侵蚀力序列的趋势性反映了降雨侵蚀力演变的总体规律和趋势。马良等[11]运用Mann-Kendall方法和径向基函数插值分析了江西省降雨侵蚀力时空分布特征，结果表明江西省1957～2008年降雨侵蚀力呈增长趋势，空间差异显著。刘斌涛等[12]采用趋势系数、气候倾向率和克里格插值分析了西南山区1960～2009年降雨侵蚀力时空变化趋势，结果表明降雨侵蚀力空间分布特征与年降雨量空间分布一致，降雨侵蚀力呈上升趋势，降雨侵蚀力变化趋势系数随海拔升高而增加。

跳跃是降雨侵蚀力序列急剧变化的一种形式，是非线性系统对外界干扰的一种响应，当序列从一种状态过渡到另一种状态时凸显出来。对一个降雨侵蚀力序列而言，跳跃前后样本常常服从不同的分布，而跳跃分析就是要找出序列发生跳跃变异的时间、次数以及变化幅度。赵平伟等[13]采用Mann-Kendall跳跃检验方法分析了滇西北地区1980～2013年降雨侵蚀序列的时空变化，发现在1983年降雨侵蚀力发生显著性跳跃。钟科元和郑粉莉[14]基于松花江流域61个气象站点1960～2014年的日降雨数据，采用标准化序列累计距平法对55 a的降雨侵蚀力序列进行跳跃点检验，并采用秩检验法对可能的跳跃点进行检验，研究表明1982和1998年为松花江流域降雨侵蚀力跳跃点。

周期作为降雨侵蚀力序列中确定性成分之一，往往是普遍存在的。李静等[15]对黄土高原降雨侵蚀力时空分布研究发现，降雨侵蚀力存在一个2.7 a的波动周期。张家其等[16]利用Mann-Kendall方法和小波分析对湖北省1957～2008年的降雨侵蚀力分析发现，1978年以来湖北省降雨侵蚀力呈上升趋势，变化主周期为2～3 a。钟科元和郑粉莉[14]采用小波周期分析发现松花江流域降雨侵蚀力序列存在15.2 a的主周期和4.7 a的小周期。

综上所述，目前关于降雨侵蚀力序列确定性成分的研究取得了丰富的研究成果，但大多研究是采用单一的检验方法对某一确定性成分进行检测往往存在很大的不确定性，缺乏对序列各类成分的综合比较，如常用的Mann-Kendall方法计算虽然直观简便，但其不适用于多个跳跃点的情形，也不适用于等级数据序列[17]。谢平等[18]根据R/S分析和分数布朗运动理论，提出了基于Hurst系数的时间序列分析方法，该方法可以从整体上判断时间序列是否存在确定性成分。本研究在基于Hurst系数对降雨侵蚀力序列整体识别的基础上，提出的基于相关系数的降雨侵蚀力时间序列分析方法，即先对降雨侵蚀力时间序列进行整体判断，是否存在确定性成分，若存在，则采用多种检验方法综合判断其存在哪些确定性成分（周期，跳跃以及趋势），若存在多种确定性成分，则进一步通过拟合分级判断哪一种确定性成分是最主要的确定性成分，对其进行剔除，重复上述工作，将降雨侵蚀力序列中的确定性成分进行分解，最终得出的降雨侵蚀力序列将是一个随机序列与确定性序列的组合。通过对降雨侵蚀力时间序列的分析和分解，以期为变化环境下区域降雨侵蚀力预测研究提供理论依据。

1 降雨侵蚀力序列分析原理与方法

1.1 基于Hurst指数的整体识别

1951年，Hurst[19]首次采用Hurst指数值来定量表征时间序列的长期相关性。偏离0.5的程度，决定了序列持续效应的程度，越接近于0，则反持续性效应越强，即未来的变化趋势与过去变化趋势相反；越接近于1，正持续性效应越强，即未来变化趋势与过去变化趋势相同。由此可根据Hurst指数识别降雨侵蚀力序列是否含有确定性成分。谢平等[18]将水文序列中的趋势、跳跃、周期等确定性成分称为水文变异，并采用类似相关系数检验法对分数布朗运动增量的相关系数进行了检验，给出了判断序列变异程度的5个等级，为了与相关系数分级方法统一，本研究将其重新划分为3个等级（表1），即将原方法的中变异、强变异和巨变异都归为强变异。实际应用中，假设检验选用显著性水平时，通常选取2个水平、(>)，当假设检验通过了显著性水平为的检验，但未通过显著性水平为的检验，即相关函数满足r≤()<r（r、r为显著水平为、的相关系数），表明在显著性水平下序列变异显著，而在显著性水平下序列变异不显著，因此本研究将其划分为弱变异；当()<r时，本研究认为其变异程度为无变异；当r≤()≤1时，本研究认为其变异程度为强变异。

表1 基于Hurst指数的确定性成分程度分级

1.2 基于相关系数的确定性成分分析

1.2.1 趋势成分分析

设有某一区域年降雨侵蚀力序列{R, t=12, 3,}，假定其存在趋势，则可以采用线性相关方程描述序列趋势，其表达式为：

式中R为年降雨侵蚀序列中的第项；、为方程的参数，η为残余项。由最小二乘法可以求出、的估计值：

表2 基于相关系数的趋势成分程度分级

1.2.2 跳跃成分分析

求出原降雨侵蚀力序列和跳跃成分序列的相关系数：

同理可求得r2, r3, …,r，本研究取相关系数r值最大点作为原降雨侵蚀力序列的跳跃点。与Hurst系数法一样，以显著性水平和为界限(>)，将跳跃成分分为3个等级，如表3所示，如果样本相关系数r满足r≥r，则拒绝原假设r=0，认为序列存在跳跃成分；否则接受原假设，认为不存在跳跃成分。

表3 基于相关系数的跳跃成分程度分级

1.2.3 周期成分分析

假定年降雨侵蚀力序列R存在周期成分，采用正弦函数来描述其周期成分，其表达式为：

通过最小二乘法计算，的值：

通过选用不同方法确定出不同的周期长度，以周期分量序列与原降雨侵蚀力序列之间的相关系数r值最大的作为主周期。

与Hurst系数法一样，以显著性水平和为界限(>)，将周期成分分为3个等级，如表4所示，如果样本相关系数r满足r≥r，则拒绝原假设r=0，认为序列存在周期成分；否则接受原假设，认为不存在周期成分。

表4 基于相关系数的周期成分程度分级表

2 实例应用

2.1 研究区概况及数据获取

长江流域位于90°33′～122°25′E、24°30′～35°45′N之间（图1），流域面积约180万平方千米，涉及19个省级行政区，流域自西而东横贯，呈多级阶梯性地形，上游为深切割高原区，中上游为中切割山地区，中下游以低山丘陵与平原区为主。气候类型属于亚热带季风气候。流域年降水量分布极不均匀，大部分地区多年平均降雨量在800～1 600 mm之间，江源地区年降雨量小于400 mm，四川、江西、湖南和湖北部分地区年降雨量大于1 600 mm，青海、甘肃以及川西高原年降雨量在400～800 mm之间。长江流域水土流失主要集中在上中游地区，其中四川、贵州、湖北、云南、湖南、重庆、青海及陕西是水土流失面积较大的地区[20]。

本研究采用长江流域1961—2013/2014年的174个气象站点的日降雨量（北京时间20～20时，24 h降雨量）数据。数据来源于中国气象科学数据共享服务网。本研究所取站点数据符合欧洲气候评估数据集采用标准：数据不少于40 a；单个站点缺失数据少于10%；每一年缺失数据低于20%且小于3个月。

2.2 降雨侵蚀力计算

本研究采用章文波等[21]参照RUSLE以半月为计算时段来估算降雨侵蚀力季节分布的方法，该方法已在中国大部分地区成功运用，测算精度较高，其计算公式如下：

图1 长江流域位置及气象站分布图

式中R为第个半月时段的侵蚀力值（MJ·mm/(hm2·h)）；为第个半月时段内的天数，将每个月的前15 d划分为上半月，其余天数即为下半月；D为半月时段内第天的日雨量，当日雨量小于12 mm时，认为不发生侵蚀，D为0；、为模型参数，可由下式计算：

式中P12是一年内侵蚀性降雨日平均雨量（mm）；P12是侵蚀性降雨年总量的多年均值（mm）。

利用式（15）计算得到逐年各半月的降雨侵蚀力，累加即可得到年降雨侵蚀力。在降雨侵蚀力空间格局分析中，采用目前地统计学使用较为广泛的Kriging插值法对各气象站点的多年平均降雨侵蚀力进行空间插值，得到空间连续分布的长江流域多年平均降雨侵蚀力分布图。

本研究选取常用的肯德尔（Kendall）秩次相关检验法和斯波曼（Spearman）秩次相关检验法分析趋势成分[22-23]；选取滑动游程检验法、Mann-Kendall检验法和贝叶斯（BSYES）检验法分析跳跃成分[10,24]；选取傅里叶（累积图）法、功率谱密度法和简单分波法分析周期成分[25-26]。所选的显著性水平为=0.05、=0.01，相应的第一Hurst系数置信限1=0.673，第二Hurst系数置信限2=0.717。

3 结果与分析

3.1 降雨侵蚀力空间分布特征

根据式（14）和（15）计算得到各气象站点的逐年各半月的降雨侵蚀力，经累加得到年降雨侵蚀力，最后计算降雨侵蚀力多年平均值，通过Kriging插值得到长江流域多年平均降雨侵蚀力分布图（图2）。长江流域多年平均降雨侵蚀力值的分布范围在110.7～15 799.9 MJ·mm/(hm2·h)之间，平均值为6 041.2 MJ·mm/(hm2·h)。流域降雨侵蚀力总体上从东部到西部呈递减趋势，其中下游的浙江、江苏、江西等地降雨侵蚀力值较大，而上游青海省则较小。174个气象站点中有16个站点年均降雨侵蚀力超过10 000 MJ·mm/(hm2·h)，这些站点主要分布在江西和浙江，除此之外还有四川雅安站、峨眉山站、湖南南岳站和安徽黄山站。从图2中可以看出，上述4个站点正位于四川盆地西南部、雪峰山东麓和皖赣交界地区降雨侵蚀力高值区域的中心点

图2 长江流域多年平均降雨侵蚀力分布

3.2 年降雨侵蚀力序列整体识别

受全球气候变化和人类活动的影响，长江流域的降水、气温等气候要素发生了显著变化，尤其是90年代后，大量的水电站开发建设对长江流域的水文过程造成了较大的干扰，导致降雨序列发生了变化。本研究将长江流域174个气象站点的年降雨侵蚀力序列进行基于Hurst指数的确定性成分整体识别与分级，结果表明长江流域174个气象站点中有130个站点年降雨侵蚀力序列变异程度不显著，即未受变化环境的影响，有31个站点降雨侵蚀力序列出现确定性成分，基于Hurst指数分级确定为弱变异，有13个站点降雨侵蚀力序列确定性成分分级为强变异，变异站点主要分布在长江中下游地区（图3）。

图3 长江流域降雨侵蚀力序列变异站点分布

降水时空变化是改变降雨侵蚀力时空分布的主要原因[7]，采用线性回归法构建多年平均降雨侵蚀力与多年平均降雨量之间的相关关系（见图4）。长江流域多年平均降雨侵蚀力值与多年平均降雨量之间呈线性相关(2=0.89)。

为分析地理位置对降雨侵蚀力的影响，采用线性回归方法分析各气象站点的经度、纬度分别与多年降雨侵蚀力的相关关系。由图5可知，降雨侵蚀力值总体上随着经度的增加而增加（2=0.61，<0.01）；降雨侵蚀力值与纬度无明显线性关系（2=0.08）。

图4 长江流域多年平均降雨侵蚀力与多年平均降雨量相关关系

关颖慧[27]计算长江流域1960—2012年均降雨量的变异系数发现，长江流域降雨量的变异系数呈现由西南向东北增加的趋势，这与本研究中长江流域降雨侵蚀力变异站点由西南向东北增加的趋势一致。许多学者研究发现，长江中下游地区受季风气候影响显著，降雨强度、极端降雨量、连续5 d最大降雨量等指标均显著增加[27-29]，这与本研究中降雨侵蚀力变异站点集中在中下游地区的结论一致。另一方面，随着经济和社会发展，近几十年长江中下游地区大量的耕地和林草地向城镇用地转变，土地利用变化改变了水文循环的下垫面条件，进而影响了降雨径流过程。

3.3 代表性气象站点年降雨侵蚀力序列确定性成分分析

重庆奉节站地处于长江上游的末端，是长江流域生态屏障的咽喉地带，受到东亚季风和印度季风的双重影响，同时又有三峡水库的反馈作用，特殊的地理位置使该区域降雨侵蚀力呈现复杂多样的变化。因此，本研究选取重庆奉节站1961—2014年的年降雨侵蚀力序列为例进行系统分析，分别对其跳跃、趋势及周期成分进行识别，并对于识别存在的确定性成分进行分级，其结果如表5所示。

图5 长江流域多年平均降雨侵蚀力与经纬度关系

表5 降雨侵蚀力序列确定性成分识别与分级表

注：跳跃识别中（+）代表在所取信度水平下该点跳跃显著，（–）表示跳跃不显著；趋势识别中（+）代表在所取信度水平下趋势显著，（–）代表该方法检测的趋势不显著；跳跃综合中“↑”代表跳跃向上，“↓”则代表跳跃向下.下同。

Note：(+) maked variation are significant at the selected level; (–) maked variation are not significant at the selected level; “↑” maked jump variation are the upward jump point; “↓” maked jump variation are the downward jump point. The same as below.

从表5中可以看出，重庆奉节站年降雨侵蚀力序列Hurst系数为0.722，大于Hurst系数上限值0.717，表明序列间存在显著变异，变异程度为强变异。本研究选取了滑动游程检验、Mann-kendall检验和BSYES检验对跳跃成分进行了识别，得到的跳跃点分别为2007年、1963年和2004年，由此可知不同方法识别出的结果并不相同；同理针对周期成分，由傅里叶累计图法、功率谱密度法和简单分波法识别得到的结果也不相同。因此，采用单一的检验方法对某一确定性成分进行检测时，往往存在很大的不确定性。进一步识别结果表明，重庆奉节的年降雨侵蚀力序列的周期成分相关系数（0.356）高于跳跃（0.262），趋势成分不显著，此时确定性成分中存在跳跃和周期性成分。通过相关系数的分析，先对周期成分进行提取，其中傅里叶（累积方差图）检验法得出的5 a的强周期，功率谱密度得出的周期4 a弱周期，简单分波得出的周期为不显著，因此选取5 a作为值序列的主周期，采用正弦曲线对该序列周期成分进行拟合得：

式中，为数据年限，取值范围为：1961≤≤2014。如果降雨侵蚀力序列无周期成分，且保持一致的话，序列的均值将是过直线X,2()第一点（=1961）的一条水平线，其方程为X,1=4 613.396 0，它反映了值序列出现周期成分前的平均情况。因此，上述值序列出现周期成分前后的差值即为重庆奉节年降雨侵蚀力序列的周期成分，如式（17）所示：

式中X为奉节年降雨侵蚀力序列的周期成分。同样假设值序列由周期性成分X和剩余成分Y组成，则剩余成分Y=R–X：

在对重庆奉节年降雨侵蚀力序列进行5 a的周期成分提取之后，对得到的新序列Y重新进行确定性成分检验，其结果如表6所示。

Y序列Hurst指数为0.742，大于Hurst指数上限值0.717，表明序列间存在显著变异，变异程度为强变异。进一步识别结果表明，Y序列的周期成分相关系数（0.302）高于跳跃（0.297），趋势成分不显著，此时确定性成分中存在跳跃和周期性成分，通过相关系数的分析，先对周期成分进行提取，其中傅里叶（累积方差图）和功率谱密度检验法得出的3 a的周期，相关系数为0.0304，不显著；简单分波法计算得出的周期为16 a的弱周期，因此选取16 a作为Y序列的主周期，并采用正弦曲线对该序列周期成分进行拟合得：

表6 降雨侵蚀力序列确定性成分识别与分级（Yt）

如果Y序列不存在周期成分，且保持一致的话，序列的均值将是过直线2()第一点（=1961）的一条水平线，其方程为=3 926.033，它反映了Y序列出现周期成分前的平均情况。由此可知，Y序列出现周期成分前后的周期差值,2-1，即为Y序列的周期成分，如式（20）所示：

根据时间序列的线性叠加原理，即假设Y序列由周期性成分和剩余成分组成，可知剩余成分Y′ t= Y-X′ t：

由此得到Y序列去除16 a周期成分后的序列如图7所示。对得到的新序列重新进行确定性成分分析，其结果如表7所示。

图7 重庆奉节降雨侵蚀力序列周期成分（序列）

从表7中可以看出，序列Hurst系数为0.709，大于Hurst系数下限值0.673，小于Hurst系数上限值0.717，表明序列间存在变异，变异程度为弱变异。进一步识别结果表明，序列存在跳跃成分，且通过相关系数的分析，Bayes方法检测出向上跳跃的强跳跃点2011年相关系数大于Mann-kendall法得出的向上跳跃的跳跃点1970年，因此选取2011年作为跳跃点对序列中的跳跃性成分进行提取。

序列1961—2014年的年降雨侵蚀力序列在0=2011年前后可分成2个系列，即1961～2011年为第一个系列（过去），2012—2014年为第2个系列（现状），如图7所示。第一个系列的均值为4 068.8 MJ·mm/(hm2·h)，第2个系列的均值为1 498.3 MJ·mm/(hm2·h)，2个系列的均值差为2 570.6 MJ·mm/(hm2·h)，因此，可以得出序列的跳跃成分J如式（22）所示：

据此求得序列去除跳跃成分后的序列S如图8所示。

S=Y′ t-J，其中J为Y序列的跳跃成分，S为剩余成分；对得到的新序列S再次进行确定性成分检验，结果表明，其Hurst系数为0.646，小于于Hurst系数下限值0.673，表明序列间不存在，变异，变异程度不显著。

表7 降雨侵蚀力序列确定性成分识别与分级表（序列）

综上所述，提取出周期以及跳跃成分之后的重庆奉节年降雨侵蚀力序列，不存在任何形式的确定性成分，满足平稳随机的要求。参照随机水文学的理论，将重庆奉节降雨侵蚀力序列采用线性叠加的方式表示：

4 结 论

本研究提出了基于相关系数的确定性成分识别方法，首先基于Hurst指数的确定性成分初步识别与分级，同时引用水文序列变异的概念，从统计学角度将降雨侵蚀力序列确定性成分分为3级（无变异、弱变异和强变异）。通过多种检验方法综合检验，将得到的结果与原序列进行相关性分析得出最可能的趋势成分、跳跃点以及周期成分，并对不同形式的确定性成分（趋势、跳跃和周期）进行分级。根据长江流域174个气象站点1961—2014年逐日降雨资料，对流域内各气象站点年降雨侵蚀力序列进行确定性成分分析与分级，主要结论如下：

流域年均降雨侵蚀力总体上从东部到西部呈递减趋势，174个气象站点中有16个站点年均降雨侵蚀力超过10 000 MJ·mm/(hm2·h)，这些站点主要分布在江西和浙江。近52年来，长江流域174个气象站点中有130个站点降雨侵蚀力序列无明显变异，有31个站点降雨侵蚀力序列出现弱变异，有13个站点降雨侵蚀力序列出现强变异，其中变异站点主要分布在长江中下游地区。

选取代表性站点重庆奉节站年降雨侵蚀力序列进行基于Hurst系数的初步识别之后显示存在确定性成分，且不同检验方法得到的结果不尽相同，针对存在多种确定性成分的序列，通过拟合分级判断哪一种确定性成分是最主要的确定性成分，对其进行剔除，重复上述工作，将降雨侵蚀力序列中的确定性成分进行一一分解，得出的降雨侵蚀力序列将是一个随机序列与确定性序列的组合。重庆奉节站为整体强变异，存在复合周期和跳跃成分，其中复合周期为5 a和16 a，向下的跳跃点为2011年。结果表明，采用单一检验方法，而不将确定性成分一一提取，无法判断检验结果的准确性，容易产生较大的误差。

[1] Nearing M A. Potential changes in rainfall erosivity in the U.S. with climate change during the 21st century.[J]. Journal of Soil & Water Conservation, 2001, 56(3): 229－232.

[2] Lu H, Yu B. Spatial and seasonal distribution of rainfall erosivity in Australia[J]. Soil Research, 2002, 40(6): 887－901.

[3] Sun G, McNulty SG, Moore J, et al. Potential impacts on climate change on rainfall erosivity and water availability in China in the next 100 years[C]//Proceedings of the 12th International Soil Conservation Conference, Beijing, China, 2002.

[4] Meusburger K, Steel A, Panagos P, et al. Spatial and temporal variability of rainfall erosivity factor for Switzerland[J]. Hydrol. Earth Syst. Sci. 2012,16: 167－177.

[5] Panagos P, Ballabio C, Borrelli P, et al. Spatio-temporal analysis of rainfall erosivity and erosivity density in Greece [J]. Catena, 2016, 137: 161－172.

[6] Sauerborn P, Klein A, Botschek J, et al. Future rainfall erosivity derived from large-scale climate models: Methods and scenarios for a humid region [J]. Geoderma, 1999, 93 (3－4): 269－276.

[7] 赖成光，陈晓宏，王兆礼，等. 珠江流域1960-2012年降雨侵蚀力时空变化特征[J]. 农业工程学报，2015，31(8)：159－167. Lai Chengguang, Chen Xiaohong, Wang Zhaoli, et al. Spatial and temporal variation of rainfall erosivity in the Pearl river basin during 1960-2012[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(8): 159－167. (in Chinese with English abstract)

[8] Yang F, Lu C. Spatiotemporal variation and trends in rainfall erosivity in China's dryland region during 1961－2012[J]. Catena, 2015, 133: 362－372.

[9] 程金文，岳大鹏，达兴，等. 陕南地区1960-2014年降雨侵蚀力变化研究[J]. 山地学报，2017，35(1)：48－56. Cheng Jinwen, Yue Dapeng, Da Xing, et al. Research on rainfall erosivity changes during 1960-2014 in the southern Shaanxi[J]. Mountain Research, 2017, 35(1): 48－56. (in Chinese with English abstract)

[10] 丁晶，邓育仁. 随机水文学[M]. 成都：成都科技大学出版社. 1988.

[11] 马良，姜广辉，左长清，等. 江西省50余年来降雨侵蚀力变化的时空分布特征[J]. 农业工程学报，2009，25(10)：61－68. Ma Liang, Jiang Guanghui, Zuo Changqing, et al. Spatial and temporal distribution characteristics of rainfall erosivity changes in Jiangxi province over more than 50 years[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2009, 25(10): 61－68. (in Chinese with English abstract)

[12] 刘斌涛，陶和平，宋春风，等. 我国西南山区降雨侵蚀力时空变化趋势研究[J]. 地球科学进展，2012，27(5)：499－509. Liu Bintao, Tao Heping, Song Chunfeng, et al. Temporal and spatial variations of rainfall erosivity in southwest China from 1960 to 2009[J]. Advance in Earth Science, 2012, 27(5): 499－509. (in Chinese with English abstract)

[13] 赵平伟，郭萍. 1980-2013年滇西北地区降雨侵蚀力变化特征[J]. 地理科学，2015，35(10)：1306－1311. Zhao Pingwei, Guo Ping. Spatial and temporal variations of rainfall erosivity in northwest Yunnan province[J]. Scientia Geographica Sinca, 2015, 35(10): 1306－1311. (in Chinese with English abstract)

[14] 钟科元，郑粉莉. 1960-2014年松花江流域降雨侵蚀力时空变化研究[J]. 自然资源学报，2017，32(2)：278－291. Zhong Keyuan, Zheng Fenli. Spatial and temporal variation characteristics of rainfall erosivity in the Songhua River Basin from 1960 to 2014[J]. Journal of natural resources, 2017, 32(2): 278－291. (in Chinese with English abstract)

[15] 李静，刘志红，李锐. 黄土高原不同地貌类型区降雨侵蚀力时空特征研究[J]. 水土保持通报，2008，28(3)：124－127. Li Jing, Liu Zhihong, Li Rui. Temporal and spatial characteristics of rainfall erosivity with different geomorphology types on the Loess Plateau[J]. Bulletin of Soil and Water Conservation, 2008, 28(3): 124－127. (in Chinese with English abstract)

[16] 张家其，龚箭，吴宜进. 基于日降雨数据的湖北省降雨侵蚀力初步分析[J]. 长江流域资源与环境，2014，23(2)：274－280. Zhang Jiaqi, Gong Jian, Wu Yijin. A preliminary study on rainfall erosion force in Hubei based on the data of daily rainfall[J]. Resources and Environment in the Yangtze Basin, 2014, 23(2): 274－280. (in Chinese with English abstract)

[17] 项静恬，史久恩. 非线性系统中数据处理的统计方法[M]. 北京：科学出版社，1997.

[18] 谢平，陈广才，雷红富. 基于Hurst系数的水文变异分析方法[J]. 应用基础与工程科学学报，2009，17(1)：32－39. Xie Ping, Chen Guangcai, Lei Hongfu. Hydrological alteration analysis method based on hurst coefficient [J]. Journal of Basic Science and Engineering, 2009, 17(1): 32－39. (in Chinese with English abstract)

[19] Hurst H E. Long-term storage capacity of reservoirs[J]. Trans ASCE, 1951, 116: 770－808

[20] 水利部长江水利委员会. 2006-2015长江流域水土保持公报[R]. 武汉：长江水利委员会，2017.

[21] 章文波，谢云，刘宝元. 利用日雨量计算降雨侵蚀力的方法研究[J]. 地理科学，2002，22(6)：705－711. Zhang Wenbo, Xie Yun, Liu Baoyuan. Rainfall erosivity estimation using daily rainfall amounts[J]. Scientia Geographica Sinica, 2002, 22(6): 705－711. (in Chinese with English abstract)

[22] Mann H B. Nonparametric test against trend[J]. Econometrica, 1945, 13(3): 245－259.

[23] Gautheir T D. Detecting trends using Spearman's rank correlation coefficient[J]. Environmental forensics, 2001, 2(4): 359－362.

[24] 周园园，师长兴，范小黎，等. 国内水文序列变异点分析方法及在各流域应用研究进展[J]. 地理科学进展，2011，30(11)：1361－1369. Zhou Yuanyuan, Shi Changxing, Fan Xiaoli, et al. Advances in the research methods of abrupt changes of hydrologic sequences and their applications in drainage basins in China[J]. Advances in Geography, 2011, 30(11): 1361－1369. (in Chinese with English abstract)

[25] 肖志国. 几种水文时间序列周期分析方法的比较研究[D]. 南京：河海大学，2006. Xiao Zhiguo. Comparison of Different Analysis Methods for the Periodicity of Hydrological Time Series[D]. Nanjing: Hohai University, 2006. (in Chinese with English abstract)

[26] 范钟秀. 中长期水文预报[M]. 南京：河海大学出版社，1999.

[27] 关颖慧. 长江流域极端气候变化及其未来趋势预测[D].杨凌：西北农林科技大学，2015. Guan Yinghui. Extreme Climate Change and Its Trend Prediction in the Yangtze River Basin[D]. Yangling: Shaanxi Northwest A&F University, 2015. (in Chinese with English abstract)

[28] Zhao G J, Hörmann G, Fohrer N, et al. Spatial and temporal characteristics of wet spells in the Yangtze River Basin from 1961 to 2003[J]. Theoretical & Applied Climatology, 2009, 98(1－2): 107－117.

[29] Chen J, Wu X, Finlayson B L, et al. Variability and trend in the hydrology of the Yangtze River, China: Annual precipitation and runoff[J]. Journal of Hydrology, 2014, 513(5): 403－412.

Identification and classification of rainfall erosivity variation based on Hurst and correlation coefficient

Qian Feng1,2, Dong Linyao2, Huang Jiesheng1※, Liu Honghu2, Han Pei2, Sun Bei2

(1.430072,; 2.430010,)

The spatial distributions and temporal trends of rainfall erosivity are critical for accurately assessing soil erosion rates, especially under the circumstances of climate change. Temporal trends of rainfall erosivity have been noted by researchers. However, reports on the methods for temporal changes of rainfall erosivity, especially the comprehensive comparison of its components (trend, jump, periodicity, and so on), are still lacking, which reduces the accuracy of assessing soil erosion risk. The single test method of rainfall erosivity series showed large uncertainties. Through the comprehensive test methods, the most reliable components could be extracted, which was an effective way to reduce the uncertainty. In this study, a joint analysis method for rainfall erosivity series based on Hurst and correlation coefficient was proposed. Firstly, the Hurst coefficient of rainfall erosivity series was calculated, and the variation was divided into 3 intervals: no variation, weak variation and strong variation. The variation components were analyzed by a variety of test methods, and correlation analysis was conducted between the variation components and the original rainfall erosivity series to extract the variation component with the largest correlation coefficient. Then this component was eliminated, and the above steps were repeated, until all the variation components were removed from the series. Finally, the original rainfall erosivity series would be a combination of random series and the variation components. In practical applications, long-term daily rainfall data from 1961 to 2013 or 2014 in 174 national weather stations were assembled to characterize the spatial and temporal patterns of annual rainfall erosivity across the Yangtze River basin. Kendall rank correlation test and Spearman rank correlation test were employed to detect the temporal trends. Sliding run test, Mann-Kendall test and Bayes test were employed to detect the jump variations. Fourier (cumulative variance chart), power spectral density and simple partial wave method were employed to detect the periodic variations. The results showed that: 1) A total of 130 stations in the 174 meteorological stations were not affected by human activities, and there was no significant variation. There were 31 stations with weak variation, and 13 stations with strong variation. Stations with rainfall erosivity variation increased from southwest to northeast, which was consistent with the trend of precipitation, and these stations were mainly located in middle and lower reaches of the Yangtze River. 2) The average annual rainfall erosivity in the Yangtze River basin was 6041.2 MJ·mm/(hm2·h). Long-term average annual rainfall erosivity decreased from east to west, ranging from 110.7 to 15 799.9 MJ·mm/(hm2·h). The value of average annual rainfall erosivity increased with the increase of longitude. There was no significant relationship between rainfall erosivity and latitude. 3) A representative weather station (Fengjie, Chongqing) was selected for a comprehensive test. Results of the test verified the feasibility of the proposed method, and also showed that the results calculated from the single test method were uncertain. Based on Hurst and correlation coefficient analysis, the variation degree of annual rainfall erosivity series in Fengjie was strong, and the variation forms were periodic and jump variations, in which the compound period was 5 and 16 a, and the downward jump point was in 2011. This method was derived from Hurst coefficient and the relationship between correlation coefficient and variation components, and could grade the levels of variation in rainfall erosivity series. The results provide valuable information for soil erosion prediction.

erosion; classification; Hurst coefficient; correlation coefficient; rainfall erosivity; spatial and temporal variation

10.11975/j.issn.1002-6819.2018.14.018

S157.1

A

1002-6819(2018)-14-0140-09

2018-01-29

2018-06-06

国家自然科学基金（51409007）；国家重点研发计划项目（2017YFC050530302）

钱 峰，男，浙江建德人，博士生，主要从事土壤侵蚀研究。Email：qianfeng@whu.edu.cn.

黄介生，男，湖北天门人，博士生导师，教授，主要研究方向为水土资源高效利用。Email：jshuanga@public.wh.hb.cn

钱 峰，董林垚，黄介生，刘洪鹄，韩 培，孙 蓓. 基于Hurst指数与相关系数的降雨侵蚀力变异识别与分级方法[J]. 农业工程学报，2018，34(14)：140－148. doi：10.11975/j.issn.1002-6819.2018.14.018 http://www.tcsae.org

Qian Feng, Dong Linyao, Huang Jiesheng, Liu Honghu, Han Pei, Sun Bei. Identification and classification of rainfall erosivity variation based on Hurst and correlation coefficient[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(14): 140－148. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2018.14.018 http://www.tcsae.org