A Method for Detecting Abrupt Change of Sediment Discharge in the Loess Plateau, China

Abstract

In recent decades, the streamflow and sediment of the Yellow River has decreased sharply, especially the sediment discharge. The factors that lead to this phenomenon have become a widely concerned problem of the whole society. The analysis of abrupt change points of hydrological series is the key to divide datum period, so it is an important work in the research of streamflow and sediment reduction cause. So far, many methods have been proposed to detect abrupt change. However, most methods have great uncertainty due to the deficiencies of irrational structure of test statistics, ideal hypothesis or subjectivity. In this paper, a new method called moving average difference method for abrupt change points detection is proposed. It is proved to be effective through comparison with four commonly used methods via both synthetic series and real data case study. The results show that the proposed method has four distinct advantages: (1) The test statistic structure of the method has physical significance and is intuitive to understand; (2) It is more accurate in abrupt change detection; (3) It can detect all of the abrupt change points at one time; (4) It can detect the abrupt changes and calculate the corresponding mutation intensity simultaneously.

1. Introduction

The Catastrophe Theory can be traced back to the late 1960s. It was first proposed by French mathematician R. Thom [1], and was developed and improved rapidly in the following decade in both theory and practical application due to its practicability. Then, the Catastrophe Theory has been applied extensively in the fields of mathematics, biology, astronomy, meteorology, social sciences, etc. The Catastrophe Theory is based on ordinary differential equations. The essence of the theory is the singularity, and the main point of it is to examine the mutation of a system or a process from a stable state to another stable state [2]. Statistically speaking, abrupt change in hydro-climate can be defined as a significant change in a statistical characteristic, e.g., in the mean, variance or trend of a time series. Currently, in hydrology and meteorology, researches are mostly focus on the abrupt change of the mean of a series [3].

So far, many methods for abrupt change detection have been developed in the case of mean mutation [4,5,6,7], including the filter test method, moving t-test, the Cramer method, the Yamamot method and order cluster analysis [8,9]. In 1945, H.B. Mann and M.G. Kendall proposed a nonparametric test method [10]. Initially, it was only used to analyze the trend of time series. Later, Gerstengarbe and Werner [11] developed this method to the widely used Mann–Kendall test, which can analyze the trend of a series as well as detect the abrupt change points. The Pettitt method is a nonparametric test method firstly proposed and used by A.N. Pettitt [12] to test abrupt change points in 1979. The method is similar to Mann–Kendall test and is also extensively used in change point test for the time being [13]. In view of the characteristics of nonstationary time series, which make it even harder to detect the abrupt change, Bernaola-Galvan in Boston University proposed a heuristic segmentation algorithm (BG), and introduced it to probe the temporal organization and heterogeneities of human heart rate interval time series [14]. Feng [15] verified the validity of the algorithm through ideal time series, and applied BG algorithm to detect and analyze the mutation points in different scales of the tree annual ring width sequence in the northern hemisphere. Methods, such as moving t-test, Cramer method and Yamamot method, are often considered to have great uncertainty, because the division and length of subsequences in the algorithms are determined arbitrarily by people. Additionally, the random selection of subsequence length will cause the drift of detected change point. Although the other prevalent methods, such as ordered clustering analysis (OC), Mann–Kendall test [3,16], Pettitt and BG segmentation method, are not affected by human factors in calculation, the rationality of those methods’ statistical structure needs to be discussed [3,17]. Furthermore, most of prevalent methods for abrupt change analysis are based on an implicit hypothesis that the series has only one change point. However, either climate or hydrologic sequence is a multi-scale system, and the size and fluctuation (anomaly) has the characteristics of periodic changes. Therefore, these methods are difficult to find all change points on different scales and levels, which is contrary to the analysis of the inducers of abrupt change [18,19].

The Yellow River is the second largest river in China, with a total length of 5464 km. It is also well-known as the sandiest river in the world, because it runs across the Loess Plateau of North China, which is characterized by broken terrain, deep soil, heavy rain, sparse vegetation and serious soil erosion. The Loess Plateau contributes about 80% of the sediment of the Yellow River (Figure 1). The construction of landscape projects, terraces, dams and reservoirs from 1970s to 1990s, and comprehensive soil and water conservations since 1990s have made great changes in the erosion environment of the underlying surface of the Loess Plateau, which lead to a dramatic decrease in runoff and sediment in the Yellow River [20,21,22,23,24,25,26]. According to statistics, the annual sediment discharge of Tongguan gauging station is reduced from 1.6 billion tons to 260 million tons, and even reduced to less than 1 million tons by 2014. At the same time, the annual runoff of Tongguan gauging station decreased from 56 billion m³ to 44 billion m³. The sharp decrease of water and sediment in the Yellow River has aroused wide attention from government and all walks of life [27]. For a long time, many scholars and research institutes have carried out researches on this issue and have obtained many achievements [18,19,28,29,30,31,32,33,34,35]. The analysis of the abrupt points of water and sediment series is the key to divide the datum period of the water and sediment problem research, which is an important content of the research of water and sediment influence factors. However, the methods used to analyze and predict the change of water and sediment are limited and uncertain, and the data adopted are not uniform. The results obtained from the same data by different methods are not exactly the same.

For this problem, a new method for detecting abrupt change of time series was proposed in this paper. It was investigated via both synthetic and real-data cases. Furthermore, the advantages and effectiveness of the new method were validated through comparison with ordered clustering analysis (OC), Mann–Kendall test, Pettitt test and BG segmentation methods. Then, this method is applied to analyze the data of runoff and sediment discharge in the Loess Plateau of the Yellow River, and investigate the influence factors of the abrupt change of runoff and sediment.

2. Materials and Methods

2.1. Moving Average Difference Method Combined with Wavelet Analysis

The detection of time series abrupt change points is to calculate the detection statistics by observed sample sequence. The process is usually based on some certain statistical hypothesis, and would be affected by three factors. The first one refers to the random error, which generally displays a zero mean and constant variance distribution, including observation error of basic sample sequence, sample selection error and so on. The second factor refers to the variation of the periodicity of sample. The variation is driven by the physical properties of the time series, which is fluctuant around a constant. When one or more physical driving factors of the time variable mutate, leading to a significant change in the mean value of the periodicity, it is the abrupt change point to be detected in the time series. The third one is effects of multiple change points. It is difficult to accurately and effectively detect all the change points, when time series has continuous multiple mutations.

Considering the three factors talked above, we should minimize or counteracting the influence of observation errors and interannual periodicity variation on abrupt changes detecting. An effective approach to eliminate the influence of random error is average calculation, and the effective means to solve the second and third problems is to take the average periodicity as an investigation variable. Taking all the above into account, it is desirable to select the physical periodicity of time series as the moving average period. Therefore, the selection of the moving average period is crucial. If the period is too short, it cannot offset observation error and interannual fluctuation error, which will increase the uncertainty of the results. Conversely, if the period is too long, it can lead to the excessive smoothing of the sequence, which will submerge the catastrophe point. In this study, the physical periodicity of time series analyzed by Wavelet Analysis is selected as the befitting moving average period.

Wavelet analysis was proposed by Molet in the early 1980s [36]. It can reveal the variety of periodicity in time series and forecast the developing tendency of the system in different time scales. The wavelet analysis has been widely applied to characterize the frequency, intensity, time position and duration of the variations for hydro-climatic time series. The wavelet analysis is to approximate a signal or function with a series of wavelet function [37]. In this study, a Morlet wavelet function $\phi \left(t\right)$, with a mean value of zero and localized in time and frequency space is defined as:

$$\phi \left(t\right)={\pi }^{-1/4}{e}^{i{\omega }_{0}t}{e}^{-t^2/2}$$

(1)

where ${\omega }_0$ denotes the dimensionless frequency, and t is the non-dimensional time frequency. The continuous wavelet transform of a discrete sequence X_i is defined as the convolution of X_i with a scaled and translated version of $\phi \left(t\right)$:

$$w_f\left(a\right)={\left|a\right|}^{-1/2}\Delta t{\displaystyle \sum _{i=1}^n{X}_i\overline{\phi }\left(\frac{i\Delta t-b}{a}\right)}$$

(2)

where the $\overline{\phi }$ represents the complex conjugate of $\phi$, a is the wavelet scale, and b is the time index. $\Delta t$ is the time interval, which is one year in our study.

The wavelet variance is defined as:

$$\mathrm{var}\left(a\right)={\displaystyle \sum _{i=1}^n{w}_f{\left(a\right)}^2}$$

(3)

The wavelet variance reflects the distribution of periodicity on the wavelet scale. The main time scale of a hydrological sequence, namely the primary periodicity, can be determined by the wavelet variance graph.

If the physical periodicity of time series $X_i\left(i=1,2,\cdots ,n\right)$ is P, then the forward moving average sequence can be constructed as:

$$M{U}_i=\frac{1}{k}{\displaystyle \sum _{j=1}^k{X}_{i-j}}i=2,3,\cdots ,n;k=\min (p,i-1);p<\frac{1}{2}n$$

(4)

Similarly, the reverse moving average sequence can be described as:

$$M{D}_i=\frac{1}{k}{\displaystyle \sum _{j=1}^k{X}_{i+j-1}}i=2,3,\cdots ,n;k=\min (p,n-i+1);p<\frac{1}{2}n$$

(5)

MU_i and MD_i denote the two subsequences before and after the sample point X_i.

We can investigate the simplest ideal time series with one abrupt change point (without observation errors and the value before and after the mutation point is constant), for instance,

$$X_i=\{\begin{array}{cc}10\hfill & 1\le i\le 25\\ 0\hfill & 26\le i\le 50\end{array}$$

(6)

The forward and reverse moving average sequence of it is:

$$\begin{array}{cccccccccccccccc}M{U}_i=& 10& \cdots & 10& 10& 10& \cdots & 10& 10& 9& 8& \cdots & 1& 0& \cdots & 0\\ M{D}_i=& 10& \cdots & 10& 9& 8& \cdots & 1& 0& 0& 0& \cdots & 0& 0& \cdots & 0\end{array}$$

(7)

Thus, the moving average difference of MU_i and MD_i is:

$$\begin{array}{cccccccccccccccc}\Delta M_i=\left|M{U}_i-M{D}_i\right|=& 0& \cdots & 0& 1& 2& \cdots & 9& 10& 9& 8& \cdots 1& 0& \cdots & 0& \end{array}$$

(8)

Series (6) and (8) are drawn in Figure 2. It is found that the point of the maximum moving average difference is exactly the mutation point, which means that the difference between the mean value of the subsequences before and after the point reaches the maximum. Furthermore, the maximum difference is the mutation intensity of the change point. Therefore, the maximum points of statistics $\Delta M_i$ indicate an abrupt change in the mean value of the subsequence of the sample, that is, the maximum points may be points of abrupt change. Thus, we proposed the detection index of the moving average difference method:

$$I_{cr}=\underset{i}{{\displaystyle \max }}\left\{\Delta M_i\right\}$$

(9)

and the intensity index of it is:

$$M_{cr}=\frac{\underset{\Delta M}{{\displaystyle \max }}\left\{\Delta M_i\right\}}{M{U}_{Icr}}$$

(10)

2.2. Other Methods Adopted in the Study

2.2.1. The Ordered Clustering Analysis (OC)

The ordered clustering analysis is a method to estimate the most likely abrupt point of hydrological time series by statistical analysis. The basic idea of this method is to find a point to divide the series to make the sum of deviation squared of the subsequences minimum, and the sum of deviation squared between subsequences maximum.

If the possible abrupt change point of the sample X_i is τ, the deviation squared sum before and after the mutation can be respectively described as

$$V_{\tau }={\displaystyle \sum _{i=1}^{\tau }{\left(X_i-\overline{X_{\tau }}\right)}^2}$$

(11)

$$V_{n-\tau }={\displaystyle \sum _{i=\tau +1}^n{\left(X_i-\overline{X_{n-\tau }}\right)}^2}$$

(12)

where $\overline{X_{\tau }}$ and $\overline{X_{n-\tau }}$ are the mean values of the two subsequences.

The total deviation squared sum of the time series is

$$S_{OC}\left(\tau \right)=V_{\tau }+V_{n-\tau }$$

(13)

When $S=\min \left\{S_{OC}\left(\tau \right)\right\}\left(2\le \tau \le n-1\right)$, τ is the optimal bipartition point, which is concluded as the most likely abrupt change point.

Generally, if there is only one abrupt change in the series, the series of total deviation squared sum presents single minimum. If there are two or more change points, then the time series change of total deviation squared sum has two or more local minima. In this way, we can detect change points according to local minima in the series of test statistic $S_{OC}\left(\tau \right)$.

2.2.2. Mann–Kendall Test

In this study, the Mann–Kendall test procedure follows Gerstengarbe and Werner [11] who developed the method to test an assumption about the beginning of the development of trend within a sample X_i of the random variable X. The null hypothesis is formulated as follows: the sample under investigation shows no beginning of a developing trend, and a rank series S_k of the progressive and retrograde rows of this sample is constructed as:

$$S_k={\displaystyle \sum _{i=1}^k{r}_i},k=2,3,\cdots ,n$$

(14)

where

$$r_i=\{\begin{array}{l}\begin{array}{cc}1& X_i>X_j\end{array}\\ \begin{array}{cc}0& X_i\le X_j\end{array}\end{array}1\le j\le i$$

(15)

The test statistic is defined as

$$U{F}_k=\frac{S_k-E\left(S_k\right)}{\sqrt{\mathrm{var}\left(S_k\right)}},k=1,2,\cdots ,n$$

(16)

where $U{F}_1=0$ is given, $E\left(S_k\right)$ and $\mathrm{var}\left(S_k\right)$ are respectively the mean and variance of the rank series S_k, which can be calculated by Formulae (17) and (18).

$$E\left(S_k\right)=\frac{k\left(k-1\right)}{4}$$

(17)

$$\mathrm{var}\left(S_k\right)=\frac{k\left(k-1\right)\left(2k+5\right)}{72}$$

(18)

UF_i is a statistic sequence calculated by time series $X_1,X_2,\cdots ,X_n$. Arrange the time series X_i in reverse order, repeat the above calculation process and make

$$\begin{array}{l}U{B}_i=-U{F}_i\left(i=n,n-1,\cdots ,1\right)\\ U{B}_1=0\end{array}$$

(19)

If the two sequences of UF_i and UB_i have intersection points and the points are within the significance level interval, the intersection points are significant abrupt change points.

2.2.3. Pettitt

This method is a non-parametric approach developed by Pettitt [12] to detect significant changes in the mean of a time series. The test uses a version of the Mann–Whitney statistic U_t,n, which is given by:

$$U_{t,n}=U_{t-1,n}+{\displaystyle \sum _{j=1}^n\mathrm{sgn}\left(X_t-X_j\right)},2\le t\le n$$

(20)

$$\mathrm{sgn}\left(X_t-X_j\right)=\{\begin{array}{l}\begin{array}{cc}+1& X_t>X_j\end{array}\\ \begin{array}{cc}0& X_t=X_j\end{array}\\ \begin{array}{cc}-1& X_t<X_j\end{array}\end{array}$$

(21)

The test statistic counts the number of times a member of the first sample exceeds a member of the second sample. The null hypothesis of the Pettitt test is the absence of a change point. The test statistic K_N and the associated probability (P) used in the test are given as:

$$K_{t_0}=\underset{1\le t\le n}{\max }\left|U_{t,n}\right|$$

(22)

$$P\left(t_0\right)\cong 2\exp \left[-6K_{t_0}^2\left(n^3+n^2\right)\right]$$

(23)

t₀ is concluded as a significant change point when $P\left(t_0\right)\le 0.5$.

2.2.4. Bernaola-Galvan Segmentation Method

The BG segmentation method is proposed by Bernaola-Galvan et al. [14] to detect abrupt changes in non-stationary time series. It regards the abrupt change detection problem of non-stationary time series as a segmentation procedure, that is, the non-stationary time series is composed of multiple subsequences with different mean values, and the purpose of this method is to find the position where the mean value difference between each subsequences reaches maximum. The test statistic S_BG is defined as Formulae (24) and (25).

$$S_{BG}=\left|\frac{\overline{X_{left}}-\overline{X_{right}}}{s_d}\right|$$

(24)

$$s_d={\left(\frac{s_{left}^2+s_{right}^2}{n_{left}+n_{right}-2}\right)}^{\frac{1}{2}}\cdot {\left(\frac{1}{n_{left}}+\frac{1}{n_{right}}\right)}^{\frac{1}{2}}$$

(25)

where $\overline{X_{left}}$ and $\overline{X_{right}}$ are respectively the mean values of the left and right subsequences of the segmentation point. $s_{left}$ and $s_{right}$ are the standard deviations of the left and right subsequences. $n_{left}$ and $n_{right}$ are the number of sample of the left and right subsequences.

Calculate the statistic $S_{BG}\left(i\right)$ for each sample point of the time series X_i successively. The larger the $S_{BG}\left(i\right)$ is, the greater the difference of mean values of the left and right subsequences. The following relationships are applied:

$$S_{BG\max }=\max \left(S_{BG}\left(\tau \right)\right)$$

(26)

$$P\left(\tau \right)\approx {\left[1-I_{\left[v/\left(v+{\tau }^2\right)\right]}\left(\delta \nu ,\delta \right)\right]}^{\gamma }$$

(27)

where $v=n-2$ is the degree of freedom, $\delta =0.40$, $\gamma =4.19\ln n-11.54$, which can be simulated by Monte Carlo simulation. $I_{\gamma }\left(a,b\right)$ is incomplete β function.

A critical value of P₀ is set. If $P\left(\tau \right)\ge P_0$, then the time series X_i can be divided into two subsequences with significant differences in the mean values at sample point $\tau$; otherwise, no segmentation is conducted.

Repeat the above algorithm for the two newly segmented subsequences. This is repeated until the subsequence is no longer separable.

2.3. Study Area and Data Set

The structure of the proposed moving average difference method is very intuitive and easy to understand. We have demonstrated the validity of the method in theoretical deduction. In order to investigate the rationality and accuracy of the method, three synthetic series were designed as the ideal validation model in which the value, periodicity and abrupt change points were specified. The three synthetic cases are single mutation point series without errors, multi-mutation points series without errors and multi-mutation points series with errors.

2.3.1. Synthetic Series

Single Mutation Point Series without Errors

First, we consider the simplest situation, single mutation point series without errors, which is constructed as Formula (6).

Multi-Mutation Points Series without Errors

The second synthetic series is designed with five abrupt change points, the values before and after each mutation point is constant without any error, which can be described as:

$$X_i=\{\begin{array}{l}\begin{array}{cc}10& 1\le i\le 15\end{array}\\ \begin{array}{cc}5& 16\le i\le 28\end{array}\\ \begin{array}{cc}3& 29\le i\le 43\end{array}\\ \begin{array}{cc}7& 44\le i\le 59\end{array}\\ \begin{array}{cc}4& 60\le i\le 76\end{array}\\ \begin{array}{cc}10& 77\le i\le 93\end{array}\end{array}$$

(28)

Multi-Mutation Points Series with Errors

The former two ideal cases are synthetic series consisting of constant subsequences without any error. In order to investigate the validity of the method in the case with error, we added random error to the second synthetic series to form the multi-mutation points series with errors, which is the third synthetic case we studied.

2.3.2. Real Cases

The Yellow River originates from the Qinghai–Tibetan Plateau in western China, flows eastward through the Loess Plateau and drains into the Bohai Sea, with a drainage area of 795,000 km². The river can be divided into three segments according to the geographic feature and river characteristics: the upper, middle and lower reaches. The upper reaches extend from the river source to the Toudaoguai hydrological station in Togtoh County, Inner Mongolia. The middle reaches are located between Toudaoguai hydrological station and Taohuayu in Zhenzhou, and flows over 1206 km through the highly erodible Loess Plateau with an area of 34.46104 km². The Loess Plateau, containing more than 30 tributaries with a catchment area larger than 1000 km², contributes more than 80% of the sediment to the Yellow River.

A data set of annual sediment discharge from six key hydrological stations in the Loess Plateau (Toudaoguai, Longmen and Tongguan in the Yellow River mainstream, Wenjiachuan, Hejin and Huaxian in the primary tributary of the Yellow River) was selected in this study to verify the validity and advantages of the moving average difference method through comparison with the other four methods. The annual sediment discharge of the six stations were provided by the Yellow River Conservancy Commission (YRCC).

The basic information of the hydrological station is displayed in

Table 1. Information of the hydrological stations.

Stations	Drainage	Longitude (E)	Latitude (N)	Controlled Area (km2)	Length of Data
Toudaoguai	mainstream of the Yellow River	111°04′	40°15′	367,898	1964–2015
Wenjiachuan	Kuye River	110°45′	38°29′	8515	1960–2015
Longmen	mainstream of the Yellow River	110°35′	35°40′	497,552	1964–2015
Hejin	Fenhe River	110°48′	35°34′	39,728	1964–2015
Huaxian	Weihe River	109°46′	34°35′	106,498	1962–2015
Tongguan	mainstream of the Yellow River	110°18′	34°37′	682,166	1950–2015

, and the specific location of the study region and stations is shown in Figure 1.

3. Results

The moving average difference method was compared with the ordered clustering analysis (OC), Mann–Kendall test, Pettitt and BG segmentation method via synthetic cases and real cases to examine its effectiveness and advantages.

3.1. Synthetic Case

3.1.1. Synthetic Case 1

The detecting results of the single mutation point series without errors by the five methods chosen for comparative study are shown in

Table 2. Detection results of single mutation point series without errors.

Methods	Tested Series	MAD	BG	OC	Pettitt	Mann–Kendall
Change points	25	25	25	25	25	No
Statistic value	10	10	1.64 × 109	5.55 × 10−13	625	/

.

The results in Table 2 show that the moving average difference method, the OC method, Pettitt test and BG segmentation method all have detected the abrupt change points exactly. The detection statistic value of moving average difference method is precisely equal to the mutation intensity of the tested series, while the detection statistic values of the other three methods vary greatly. This illustrate that the moving average difference method has physical significance of mean value abrupt change, while other methods are only statistical.

The sequence of test statistics of the five methods is depicted in Figure 3 and Figure 4. In order to display the sequence of different detection indexes in one figure, the sample series and the detection statistics sequences of the methods (except Mann–Kendall test) are standardized by the following formula,

$$x_i^{\prime }=\frac{x_i-\min \left\{x_i\right\}}{\max \left\{x_i\right\}-\min \left\{x_i\right\}}$$

(29)

As shown in Figure 3, there is an obvious identical turning point in the detection statistics sequences of the methods adopted (except Mann–Kendall test), which is the same as the abrupt change point of the tested series. Figure 4 shows that the two sequences of Mann–Kendall test have no intersection point at all, so Mann–Kendall test failed to detect the abrupt change point of the sequence in this case.

3.1.2. Synthetic Case 2

The detecting results of multi-mutation points series without errors are listed in

Table 3. Detecting results of multi-mutation points series without errors.

Methods	Test Series	MAD	BG	OC	Pettitt	Mann–Kendall
Change point 1	15	15	15	15	15	No
Statistic value 1	5	5	54.6	503.2	915	/
Change point 2	28	28	no	no	no	No
Statistic value 2	2	2	/	/	/	/
Change point 3	43	43	no	no	no	No
Statistic value 3	4	4	/	/	/	/
Change point 4	59	59	no	no	no	No
Statistic value 4	3	3	/	/	/	/
Change point 5	76	76	76	76	76	No
Statistic value 5	6	6	60.3	468.6	1037	/

. The results show that only the moving average difference method has accurately detected all the abrupt change points and calculated their mutation intensity at one time without any deviation. The BG algorithm, OC method and Pettitt test have detected the first and the fifth abrupt change points, while the second, third and fourth mutation points are not detected. According to these methods’ algorithm, the remaining abrupt change points can be detected when the subsequence of detected abrupt points is removed successively. The Mann–Kendall test has detected none of the five abrupt change points in this case.

The sequences of the test statistics of the five methods in this case are plotted in Figure 5. As indicated in Figure 5, the sequence of the test statistics of the moving average difference method has five local maximum points, which are exactly the locations of five abrupt change points of the synthetic series. From the sequence of the detecting statistics of the BG algorithm and OC method, we can see two local extremum points, corresponding to the first and the fifth mutation points respectively. The sequence of the test statistics of Pettitt test showed three local extremum points of point 1, point 3 and point 5, but the third mutation point did not reach the significant level. Therefore, the Pettitt test also has only detected two abrupt change points of point 1 and point 5. For Mann–Kendall test, similar to the result in the first synthetic case, the two sequences of UF and UB have no intersection point at all, none of the five abrupt change points is detected.

Comparing the test statistic value of each method with the mutation intensity at each point in Table 3, it is found that the test statistic value of each method is in direct proportion to the mutation intensity of each change point. The relationship between the test statistic value of each method and the mutation intensity at each change point are plotted in Figure 6. In order to display the sequence of test statistic values of different methods in one figure, the test statistics values of the methods (except Mann–Kendall test) are transformed by the following formulas,

$$M_{BG}=6\times S_{BG}/60.3$$

(30)

$$M_{OC}=6\times 68.6/\left(S_{OC}-400\right)$$

(31)

$$M_P=6\times \left(U_{t,n}-400\right)/637$$

(32)

where S_BG, S_OC and U_t,n are the test statistic values of the BG algorithm, OC method and Pettitt test, respectively. M_BG, M_OC and M_P refer to the transformation of S_BG, S_OC and U_t,n.

As shown in Figure 6, the test statistic value of the moving average difference method is equal to the actual mutation intensity of each change point. The moving average difference method can detect all the change points and calculate their mutation intensity accurately at one time because the structure of the test statistic of this method does not lose effective information of the sample in the calculation process. However, the test statistics of the other three methods will inevitably lose effective information of the sample during the algorithm. The lower the mutation intensity is, the greater the loss. As we can see in Figure 5, the mutation intensity of change points 2 and 4 is lower than the others. Hence, the three methods failed to detect the two change points due to the information loss. This is a common shortcoming of the three methods.

3.1.3. Synthetic Case 3

The detection results of multi-mutation points series with errors are presented in

Table 4. Detection results of multi-mutation points series with errors.

Methods	Tested Series	MAD	BG	OC	Pettitt	Mann–Kendall
Change point 1	15–16	15–16	15–16	15–16	17–18	no
Statistic value 1	5	4.97	45	663	938	/
Change point 2	28–29	28–29	no	no	no	no
Statistic value 2	2	1.83	/	/	/	/
Change point 3	43–44	45–46	no	no	no	no
Statistic value 3	4	3.2	/	/	/	/
Change point 4	59–60	59–60	no	no	no	no
Statistic value 4	3	2.9	/	/	/	/
Change point 5	76–77	76–77	76–77	76–77	76–77	85–86
Statistic value 5	6	6.3	49.2	600	1006	/

.

Again, as shown in Table 4, the detection results are similar to that of the second case. Only the moving average difference method succeed to detect all the five abrupt change points. However, there is a small deviation in the calculation of mutation intensity due to the random error of the series. The BG algorithm, OC method and Pettitt test have detected the change point 1 and point 5, while the point 2, point 3 and point 4 are not detected. The results of Mann–Kendall test differ greatly. This method detects the abrupt change point of the sequence at the sample points 85–86, but the point is obviously not the actual abrupt change point of the series.

3.2. Real Case

The wavelet transform was applied to analyze the periodicities in annual sediment discharge at the six hydrological stations. Figure 7 and Figure 8 illustrate the wavelet analysis of the annual sediment discharge at Tongguan station. As indicated in Figure 7, there are multi-time scales in the annual sediment discharge at Tongguan station. In general, there are four types of time scales, namely, 45~60 years, 30~40 years, 10~20 years and 3~7 years. The former two periodicity scales are stable during the whole analysis period. But the later two periodicity scales are not full-scale. They weakened or even disappeared after 1990, because of the intensive comprehensive soil and water conservations.

There are four distinct peaks in the curves of wavelet variance of sediment discharge at Tongguan station (Figure 8), which are 54 years, 16 years, 35 years and 5 years, respectively. The maximum peak corresponds to the time scale of 54 years, indicating that the periodicity around 54 years is the strongest, which is the primary periodicity of the sediment discharge at Tongguan station. By parity of reasoning, the second peak, which corresponds to 16 years, is the second major periodicity of the series. Since the moving period p < $\frac{1}{2}$n, p = 16 was selected in this study.

The detecting results by five methods at Tongguan station are presented in Figure 9. As shown in Figure 9, two significant abrupt change points were detected by the moving average method and the corresponding mutation intensities of year 1979 and 2003 were 6.528 and 5.467, respectively. The OC method and the Pettitt test detected the abrupt change in 1979 and the abrupt change detected by the BG algorithm is in 1996. For Mann–Kendall test, the two test statistics UF and UB have an intersection in 2000, but the intersection point is not within the α = 0.05 confidence interval (u_0.05 = ±1.96), so the intersection is not a significant abrupt change. It can be clearly seen that the test statistics sequence of the BG algorithm, OC method and Pettitt test all have two obvious local extremum points in the curve. The first local extremum point is in 1979. According to the analysis series and the test results of the moving average difference method, the mutation intensity in 1979 is the largest, so 1979 is a significant abrupt change point of the analysis series. The second local extremum point is in 1996. It can be seen from the test statistics sequence curve of the moving average difference method that 1996 is indeed a change point. However, the mutation intensity of year 1996 is 4.868, which is comparatively smaller than 1979 and 2003. The abrupt change detected by the BG algorithm is not in 1979, when the mutation intensity is the strongest, but in 1996, when the mutation intensity is comparatively weak. Furthermore, in 2003, the mutation intensity should be second only to the abrupt change in 1979. But the BG algorithm, OC method and Pettitt test have not detected this abrupt change. This is caused by the implicit assumption in these methods that there is only one abrupt change point in the sequence. In summary, there are two significant abrupt change points in the annual sediment discharge at Tongguan station, which occurred in 1979 and 2003 respectively.

The detecting results by different methods of the annual sediment sequences at the other five hydrological stations (Toudaoguai, Longmen, Wen Jiachuan, Hejin and Huaxian) are shown in

Table 5. Detecting results of sediment discharge of other stations.

Stations	Abrupt Change	MAD	OC	Pettitt	BG	Mann–Kendall
Toudaoguai	Change point	1986	1985	1986	1985	1985
	Statistic value	0.76	13.77	561	29.44	/
Longmen	Change point 1	1979	1979	1979	1979	No
	Statistic value	5.76	711.68	468	27.44	/
	Change point 2	1996	1996	1996	1996	No
	Statistic value	4.31	840.43	591	29.05	/
Wenjiachuan	Change point 1	1979	1979	1979	1979	No
	Statistic value	0.66	27.86	460	21.94	/
	Change point 2	1996	1996	1996	1996	No
	Statistic value	0.69	26.76	663	29.83	/
Hejin	Change point 1	1979	1979	No	1979	No
	Statistic value	0.24	0.77	/	26.08	/
	Change point 2	1996	No	1996	No	No
	Statistic value	0.044	/	585	/	/
Huaxian	Change point 1	1979	1979	No	No	No
	Statistic value	1.75	204.85	/	/	/
	Change point 2	1996	1996	1996	1996	No
	Statistic value	1.48	204.34	503	23.96	/
	Change point 3	2003	No	No	2003	No
	Statistic value	1.87	/	/	24.14	/

and Figure 10.

Four conclusions can be drawn from Table 5 and Figure 10:

(1) When there is only one significant abrupt change point in the hydrological series (for example, the annual sediment discharge sequence at Toudaoguai station), all the five methods we investigate can detect the abrupt change, and the results are nearly the same.

(2) When there are two or more than two abrupt change points in the hydrological series, only the moving average difference method can accurately detect all abrupt change points and calculate the corresponding mutation intensity at one time. The abrupt change point detected by the BG algorithm, OC method and Pettitt test is one of the change points detected by the moving average difference method. But it is not necessarily the point of the maximum mutation intensity, and even may be the weakest mutation point among the abrupt change points. This illustrates that there is a problem of catastrophe points drift in the three methods.

(3) In many situations, the Mann–Kendall test cannot detect the abrupt change in the time series or the detecting result is wrong.

(4) The sediment discharge in the upper reaches of the Yellow River changed abruptly in 1986 (annual sediment discharge at Toudaoguai station). The water used for production and living in the upper reaches of the Yellow River increased greatly since 1980s, about 52 percent more than in the 1950s [27]. Besides, the construction of water conservancy projects in the upper reaches of the Yellow River after 1970s has greatly changed the annual and interannual distribution of water and sediment in the basin. For example, after the completion of Longyangxia reservoir in 1986, it was jointly operated with Liujiaxia reservoir, which impounded a large amount of runoff and sediment, greatly reducing the sediment discharge in the upstream [38]. The sediment discharge in the middle reaches of the Yellow River mainly has two abrupt changes, in 1979 and 1996, respectively. The middle reaches of the Yellow River run across the Loess Plateau of North China. Large-scale soil and water conservation measures in the Loess Plateau since 1970s have greatly accelerated the reduction of water and sand sediment in this region [35,39]. In the early 1990s, the second stage water and soil conservation management project was carried out, leading to the further increase of vegetation coverage and further decline in runoff and sediment. Therefore, there are two abrupt change points, respectively in 1979 and 1996, in the annual sediment discharge sequences at Longmen, Wenjiachuan and Hejin stations. There are many influence factors in the middle reaches, and the further downstream, the more complex the factors are. The Huaxian station is the control gauging station of the Weihe River, near the Tongguan station. Abrupt changes of the annual sediment discharge at the Huaxian station occurred in 1979, 2003 and 1996, which were consistent with the results of the Tongguan station. Heavy coal exploitation after 2000 and the increase of social water use in the middle reaches of the Yellow River are the dominating factors causing the abrupt change in 2003.

4. Discussion

The authors consider that the main reason for the effective information loss of the three methods (the BG algorithm, OC method and Pettitt) is the unreasonable statistical structure of the methods. There is an implied hypothesis in the algorithm of the three methods that there is only one abrupt change point and the subsequences before and after the change point are stationary sequences with different mean values. Therefore, when there are multiple abrupt change points, the points of low mutation intensity are concealed by those of high mutation intensity.

Such as the BG segmentation algorithm, as we can deduce from the Formulas (24) and (25), the value of the structure ${\left(\frac{1}{n_{left}}+\frac{1}{n_{right}}\right)}^{1/2}$ consequentially reaches the maximum at the two tails of the sequence (where n_left·n_right reaches the minimum) and reaches the minimum in the middle of the sequence (where n_left·n_right reaches the maximum). Therefore, this section in the formula has no relation with the abrupt change. It cannot excavate the effective information of the sample, on the contrary, it provides harmful information that interferes with the detection of catastrophe points.

As for Mann–Kendall test, this non-parametric test is based on the fact that, under the hypothesis of a stationary time series [1,10], which has no systematic change in mean value, variance and periodicity, the succession of climatological values must be independent and the probability distribution must remain always the same. Therefore, it is inapplicable for non-stationary series. In recent decades, large-scale human activities have led to the abrupt change in annual runoff and sediment discharge in the Yellow River, especially in mean value. Consequently, the sediment discharge series in the Yellow River basin is not a stationary sequence. Hence, the Mann–Kendall test is not applicable in such problems.

5. Conclusions

(1) A new method named moving average difference method for detecting abrupt change is proposed in this paper. It can detect the multiple abrupt change points of time series and calculate the corresponding mutation intensity at the same time. The method was investigated via both synthetic data and real data, showing that the method works as theoretically expected. Furthermore, it was also applied to validating the advantages and effectiveness through comparison with four widely used methods (the clustering analysis method, Mann–Kendall test, Pettitt test and BG segmentation method). It is found that the traditional four methods have the following deficiencies: First, there is an implied hypothesis in the algorithm of the four methods that there is only one abrupt change point in the series. When the real data does not conform to this assumption, the algorithm will lead to the loss of effective information, so that the abrupt change point with low mutation intensity cannot be detected. Second, the test statistics structure of the BG segmentation method is unreasonable and remains to be improved. In summary, the moving average difference method has four advantages: Firstly, the test statistic structure of the method has physical significance and is intuitive and simple to understand. Secondly, abrupt change detection is more accurate than other methods. Thirdly, it can detect all of the abrupt change points at one time. In addition, it can simultaneously detect the abrupt changes and calculate the corresponding mutation intensity. Therefore, the method proposed is worthy of being popularized.

(2) The annual sediment discharge series at six key hydrological stations in the Yellow River are analyzed. The results show that the sediment discharge in the upper reaches of the Yellow River changes abruptly in 1986. The sediment discharges in the middle reaches of the Yellow River (Longmen station), the Kuyehe River (Wenjiachuan station) and the Fenhe River (Hejin station) change abruptly in 1979 and 1996. Furthermore, 1979 was the most significant catastrophe point. There are three abrupt change points in the annual sediment discharge series at Huaxian station and Tongguan station, respectively in 1979, 2003 and 1996, with diminishing mutation intensities. Human activities, including the construction of water conservancy projects, large-scale soil and water conservation measures and coal mining are the primary factors that lead to the abrupt change of sediment discharge in the Yellow River.

(3) Little attention is paid to the variation of streamflow in the Yellow River as well as in other regions in this paper [40,41,42]. More work is needed to investigate the effectiveness and superiority of the proposed method. Therefore, we will apply the method to runoff discharge in the Yellow River as well as in other regions in the following research.