A Modified Multifractal Detrended Fluctuation Analysis ( MFDFA ) Approach for Multifractal Analysis of Precipitation in Dongting Lake Basin , China

Multifractal detrended fluctuation analysis (MFDFA) method can examine higher-dimensional fractal and multifractal characteristics hidden in time series. However, removal of local trends in MFDFA is based on discontinuous polynomial fitting, resulting in pseudo-fluctuation errors. In this paper, we propose a two-stage modified MFDFA for multifractal analysis. First, an overlap moving window (OMW) algorithm is introduced to divide time series of the classic MFDFA method. Second, detrending by polynomial fitting local trend in traditional MFDFA is replaced by ensemble empirical mode decomposition (EEMD)-based local trends. The modified MFDFA is named OMW-EEMD-MFDFA. Then, the performance of the OMW-EEMD-MFDFA method is assessed by extensive numeric simulation experiments based on a p-model of multiplicative cascading process. The results show that the modified OMW-EEMD-MFDFA method performs better than conventional MFDFA and OMW-MFDFA methods. Lastly, the modified OMW-EEMD-MFDFA method is applied to explore multifractal characteristics and multifractal sources of daily precipitation time series data at the Mapoling and Zhijiang stations in Dongting Lake Basin. Our results showed that the scaling properties of the daily precipitation time series at the two stations presented a long-range correlation, showing a long-term persistence of the previous state. The strong q-dependence of H(q) and τ(q) indicated strong multifractal characteristics in daily precipitation time series data at the two stations. Positive ∆ f values demonstrate that precipitation may have a local increasing trend. Comparing the generalized Hurst exponent and the multifractal strength of the original precipitation time series data with its shuffled and surrogate time series data, we found that the multifractal characteristics of the daily precipitation time series data were caused by both long-range correlations between small and large fluctuations and broad probability density function, but the broad probability density function was dominant. This study may be of practical and scientific importance in regional precipitation forecasting, extreme precipitation regulation, and water resource management in Dongting Lake Basin.


Introduction
Precipitation is arguably the most important variable in water cycle and land-atmosphere systems accountable for shaping the climatic state and variability of water on the land surface, in soil, and in the atmosphere [1].Hydrological processes and the water cycle have been greatly affected by human activities during the past century, such as deforestation, emission of carbon dioxide, exploitation of water resources, and so forth, which has also led to evident global climate change [2].Therefore, increasing attention has been paid to the regulation of water resources and hydro-environmental protection [3].A better understanding of the evolution law of precipitation is crucial to further grasp land surface hydrological processes for hydrometeorological modeling, flood control, water resources management, and ecological conservation.
Precipitation variations are extremely complex because precipitation in hydrometeorological processes has been severely affected by global climate change and human activities [4].An important way to understand the dynamic evolution law of precipitation is to investigate the inner mechanisms of precipitation time series data.Long-range correlation of precipitation time series data is a useful characteristic for future forecasting.In 1951, an England scholar named Hurst first studied long-range statistical dependencies [5].Later, similar long-term correlated characteristics have also been reported for precipitation, temperature, runoff, and other geophysical variables [6][7][8][9][10].Earlier approaches exclusively focused on either the absolute values or the variances of the full distribution of fluctuations, which can be regarded as the first moment F 1 (s) and the second moment F 2 (s) of the fluctuations, respectively, in giving time segments of length s [11].After that, many scholars studied scaling properties and long-range correlations of time series data [12][13][14][15][16][17][18][19][20][21].However, the original Hurst exponent calculation method was rescaled range (R/S) analysis and spectral analysis, which have both been criticized by many scholars and may led to spurious results in presence trends [22][23][24].Then, Peng et al. [25] proposed detrended fluctuation analysis (DFA) for scaling property estimations when predicting DNA sequences.DFA was common and widely employed to study long-range correlation and fractal scaling properties of hydrometeorological variable time series data.Matsoukas et al. [20] utilized DFA to quantify the correlation properties of rainfall time series from various hydrologic and climatic regions in the United States; their results showed a common scaling behavior.Kantelhardt et al. [11] applied DFA to detect long-term correlations of precipitation.They found that precipitation has scaling behavior corresponding to a rapid decay of the autocorrelation function.To determine the scaling behavior in the presence of trends of daily precipitation records, Rybski et al. [26] applied DFA.They found that daily runoffs were characterized by a power-law decay of the autocorrelation function above some crossover time that was usually several weeks.Previous research on multifractal characteristics of precipitation time series data focused on analysis of individual sites, small regions, or a short time period [27][28][29][30][31][32].Inspired by the DFA method, Kantelhardt et al. [33] proposed a simple and powerful multifractal analysis method by considering the qth-order of the fluctuation to give a full description of a more complicated time series; this was named multifractal detrended fluctuation analysis (MFDFA) [34].Kantelhardt et al. [35] utilized MFDFA to study multifractal temporal scaling properties of river discharge and precipitation records.Zhang et al. [36] used the MFDFA technique to explore scaling and multifractal properties of long daily mean streamflow series from Cuntan, Yichang, Hankou, and Datong stations at the Yangtze River basin.They found that the streamflow series of the upper Yangtze River basin was characterized by anti-persistence, while the streamflow series of the lower Yangtze River basin was characterized by long memory.Yuan et al. [16] proposed an improved MFDFA to investigate the long-term monthly runoff records of runoff in the Yangtze River.They found that the monthly runoff series was characterized by long-range correlation and multifractality, and the multifractal source was dominated by a broad probability density function.Zhou et al. [37] proposed the multifractal moving-window detrended fluctuation analysis (MF-MWDFA) and the multifractal temporally weighted detrended fluctuation analysis (MF-TWDFA) methods, and applied these two modified methods to scaling behavior analyses in temperature series.However, dividing and detrending are two key steps in DFA and MFDFA.In classic DFA and MFDFA, the time series is divided into I s = int(n/s) non-overlapping segments of equal length s.Since the total length n of time series data is often not an integer multiple of the segmentation length s, it causes the tail of the time series data to be discarded.To avoid this, the above segmentation is repeated in reverse order of the original time series data.However, the problem is that reverse segmentation disrupts the order of the original time series data and also produces corresponding errors.Another problem is that if there is a trend in a real-time series, the functional form of this trend is usually unknown.A conventional strategy is to assume that local trends are in the form of a polynomial.Discontinuity at the segmentation points using polynomial fitting results in new pseudo-fluctuation errors and estimation errors in the qth-order generalized Hurst exponent [37,38].In 2009, Wu et al. [39] proposed a new model-free, timescale-adaptive detrending approach, which was called empirical ensemble mode decomposition (EEMD).The EEMD method can decompose time series data into a limited number of intrinsic mode functions (IMFs) and a residual component.The trend component obtained by EEMD can be considered a combination of the monotonic residue and IMFs, which are significantly distinguishable from pure white noise, and the timescale of an EMED-based trend is the averaged timescale of the decomposed IMFs.In order to solve the two problems, we propose a two-stage modified MFDFA method for multifractal analysis in this study.First, we introduce the overlap moving window (OMW) algorithm to divide the time series data.Second, we employ EEMD to obtain the local trend component to replace the polynomial local trend component.Thus, the modified MFDFA method is called OMW-EEMD-MFDFA.Then, numeric simulation experiments will be conducted to measure the performance of the OMW-EEMD-MFDFA.Lastly, we apply the proposed OWM-EEMD-MFDFA method to explore complex multifractal characteristics and also understand the multifractal source of daily precipitation time series in the Dongting Lake Basin.
The paper is organized as follows: Section 2 describes the MFDFA, OMW, EEMD, and the modified OMW-EEMD-MFDFA methods; Section 3 provides the performance of the OMW-EEMD-MFDFA method with extensive numerical experiments; the OMW-EEMD-MFDFA method is applied in Section 4 to investigate multifractal characteristics and the multifractal source of daily precipitation time series data in the Dongting Lake Basin; and finally, Section 5 presents the conclusion of the paper.

MFDFA Algorithm
Based on the DFA, the MFDFA was proposed by Kantelhardt et al. [33].This method is considered to be an effective tool for measuring whether multifractal characteristics exist in time series data such as hydrometeorological time series.
For time series X = {x τ } n τ=1 , n is the length of time series data, and the detailed computational procedure can be described as follows [34,40]: Step 1: calculate the profile time series Y = y τ n τ=1 .
where x is the mean value of X = {x τ } n τ=1 .

•
Step 2: divide the profile time series Y = y τ n τ=1 into I s = int(n/s) non-overlapping segments of equal length s.In segmentation, the total length n of time series data is often not an integer multiple of the segmentation length s, which causes the tail of the profile time series data to be discarded.To avoid this, the above segmentation process was repeated in reverse order of the profile time series data, and a total of 2I s segments were obtained.

•
Step 3: Calculate the trend of each 2I s segment, and the detrended time series z t (i) can be obtained by the following Equation: where v t (i) and p t (i) represent the segment time series and trend time series at each segment t, respectively.

•
Step 4: Determine the variance F 2 (s, t) by the following Equation: and then the qth-order fluctuation function of all 2I s segments are obtained by the following Equation:

•
Step 5: There is a power-law relationship between the fluctuation function F q (s) and scale s.
The scaling behavior of the fluctuation function could be determined by analyzing log-log plots of F q (s) versus scale s for each value of q [41].Then, the least squares function was employed to fit the ln F q (s) and ln s, and the slope was considered as the generalized Hurst exponent H(q).If H(q) and q were independent to each other, H(q) had constant values.This meant that the time series data did not have multifractal characteristics.On the contrary, if H (q) was significantly dependent on q, then the original time series data had multifractal characteristics.For positive values of q, H(q) described the scaling behavior of segments with large fluctuations, while for negative values of q, H(q) depicted the scaling behavior of segments with small fluctuations [14,42].In addition, if q was equal to 2, H(2) was known as the Hurst exponent, which was proposed by Hurst, a British hydrologist, in 1951 to study whether time series data had long-range correlation characteristics [36].A scaling exponent H equal to 0.5 indicated that the time series data were uncorrelated.The range 0.5 < H < 1 denoted long-range correlation (long memory, or long-term persistence), while the range 0 < H < 0.5 meant short-range correlation (anti-memory, or anti-persistence) [36].∆H = H(q min ) − H(q max ), measured the multifractal strength of the time series data.

Multifractal Strength Analysis
The relationship between the generalized Hurst exponent H(q) and the multifractal scaling exponent τ(q) can be calculated by the following Equation: If τ(q) and q had a linear relationship, then the time series data was monofractal.On the contrary, if τ(q) and q had a nonlinear relationship, the time series data was multifractal; the stronger their nonlinear relationship was, the stronger the multifractal characteristics were.
For multifractal time series data, another way to describe its multifractal characteristics is the singularity spectrum f (α), which can be determined from τ(q) via the first-order Legendre transform.
where α is the Hölder exponent or singularity index, which was utilized to describe the singularity degree of the time series.The physical meaning of α was: if α = 1, it meant that the distribution of the time series data was uniform; if α was larger, the singularity degree was smaller.On the contrary, if α was smaller, the singularity degree was larger.f (α) represented the dimensions of the segments of the time series data that were characterized by α.The width of spectrum f (α) represented the strength of multifractal effects in the time series data.For monofractal time series data, the spectrum f (α) was a single point, and both functions τ(q) and H(q) were linear.
In order to analyze and understand the multifractal characteristics of time series data, two multifractal parameters were used to deduce from an approximate α ∼ f (α) spectrum.The width of the spectrum can be defined as: ∆α = α max − α min (11) where α max and α min are the maximum and minimum of the α.∆α represents the inhomogeneity of the distribution of the probability measured on the overall fractal structure.To be specific, the larger ∆α was, the larger the multifractal characteristics intensity and complexity of the time series data were; the distribution was more uneven and the fluctuation was more severe in the time series.
The difference in fractal dimensions between the minimum probability and maximum probability segments can be obtained by the following Equation: where f (α min ) and f (α max ) indicate the values of f (α) at the specific α min and α max , respectively.∆ f describes the proportion of the number of elements at the minimum and maximum values in the segments, indicating the proportion of the small and large peak values of the fluctuation time series data.
In order to measure quantitative characterizations, the multifractal spectrum can be fitted to a second-order equation around the position of its maximum at α 0 [43,44] by the following function: where the coefficients A, B, and C are obtained by the least squares method.α 0 is the corresponding α to the maximum f (α).Parameter B represents the asymmetry of the multifractal spectrum.A B of zero represented symmetric shape of the spectrum, and positive and negative B values represented right-skewed and left-skewed shapes, respectively.

Overlap Moving Window (OMW)
In order to solve the problem of discontinuity of data segmentation, in this study we introduced the overlap moving window (OMW) algorithm.The workflow of OWM algorithm is shown in Figure 1.For time series data with a total length of N, first determine the size s of the moving window.Then determine the step size k of each forward movement, and the number of movement times can be obtained by subtracting the size s of the moving window from the total length N of time series data.

Overlap Moving Window (OMW)
In order to solve the problem of discontinuity of data segmentation, in this study we introduced the overlap moving window (OMW) algorithm.The workflow of OWM algorithm is shown in Figure 1.For time series data with a total length of N, first determine the size s of the moving window.Then determine the step size k of each forward movement, and the number of movement times can be obtained by subtracting the size s of the moving window from the total length N of time series data.

Ensemble Empirical Mode Decomposition (EEMD)
Ensemble empirical mode decomposition (EEMD) is a noise-assisted data analysis method based on empirical mode decomposition (EMD), which can empirically and adaptively process nonlinear and nonstationary time series data with high efficiency [39].The EMD method can decompose a time series data into a limited number of intrinsic mode functions (IMFs) and a residual component.More details about the EMD can be found in reference [45].Although EMD has obvious advantages in time series data processing, there are also some unavoidable weaknesses.The major weakness is mode-mixing.Mode-mixing means that the same IMF contains different frequency components, or the frequency of the same and similar scales is distributed in different IMFs.Thus, mode-mixing not only causes mixing of various scale vibration modes, but may even change the physical meaning of the individual IMF.In order to overcome these weaknesses, the ensemble EMD method was proposed as a powerful improvement to the EMD with the idea of ensemble mean.The main steps are expressed as follows: • Step 1: Add white noise w i (t) to the original time series data x(t). • Step 2: Decompose the new time series data x i (t) to many IMFs using the EMD method.

•
Step 3: Repeat steps 1 and 2 with a different white noise each time.

•
Step 4: Obtain the ensemble means of the corresponding IMFs of the decompositions as the final results.
For the EEMD method, the first important step is to determine the ensemble times and the amplitude of adding noise.But how to select the best ensemble number and amplitude of adding noise is still an open question.The effect of adding white noise should obey the following statistical rule: where ε n is the final standard deviation of the error, which is described as the difference between the input signal and the corresponding IMFs.ε denotes the amplitude of the added noise, and N is the number of ensemble numbers.In this work, the ensemble number was set to 100, and the amplitude of added white noise was set to 0.2 times the standard deviation of the corresponding data [40].

The Modified OMW-EEMD-MFDFA Method
In step 2, the profile time series Y = y τ n τ=1 was divided into I s = int(n/s) non-overlapping segments of equal length s.The problem was that discontinuity at the segmentation points using polynomial fitting resulted in a new pseudo-fluctuation error and estimation error in the qth-order generalized Hurst exponent.In order to avoid the above error, the OWM method was introduced to improve segmentation of the time series data in step 2. Specifically, modification of the MFDFA method was only in steps 2 and 4, keeping the other steps unchanged.
Step 2: divide the profile time series into n − s + 1 overlapping segments of equal length s by the OWM algorithm.
Step 4: The original time series data were divided into n − s + 1 segments by the OMW algorithm.Thus, the qth-order fluctuation function F q (s) needed to be modified as follows: For real-time series data, it is usually unknown whether there is a trend component and, if so, what the functional form of this trend is.The conventional strategy in the MFDFA method is to assume that the local trends are in the form of a polynomial.Now, EEMD is embedded into MFDFA to modify Step 3 of the method while keeping all other steps unchanged.
Step 3: For each segment v t (i), the local trend component r n (i) is obtained by the EEMD method replaced by conventional polynomial fitting.
Note that the trend r n (i) should be determined for each 2I s segment on each timescale s.

Experimental Analyses
The multiplicative cascading method, as the widely utilized method to model multifractal tests in many complex systems, was introduced to test the performance of the modified OWM-EEMD-MFDFA method [46,47].The p-model algorithm was the easiest and simplest one in the multiplicative cascading method, and was successfully applied to simulate the energy dissipation field in turbulent flows [48].In this study, the p-model algorithm was employed to generate multifractal simulation data to test the stability and performance of the modified OMW-EEMD-MFDFA method.The p-model algorithm is defined as follows: here, 0.5 < p < 1, k = 1, 2, 3, . . ., N, N = 2 n max is the total length of the generated multifractal simulation data; n max is a power of 2; and n(k) is a value of 1 in the binary representation of index k (e.g., n(10) = 2, since 10 corresponds to binary 1010).If n max is infinite, the scaling exponent τ(q) can be expressed as follows: The Hurst exponent H(q) can be determined by Equations ( 6) and (19).
Then, the singularity exponent α and the singularity spectrum f (α) can be obtained via Equations ( 7) and (8).
In order to measure the performance of the proposed OMW-EEMD-MFDFA method, the OMW-EEMD-MFDFA method was compared with the MFDFA and OMW-MFDFA methods by adopting multifractal simulation data generated by the p-model.The experiment was carried out in Matlab 2018a on the graphic workstation machine.Before the experiment, some parameters needed to be set.In the p-model, the p was set to 0.75 and n max was set to 12.The length of the generated multifractal simulation data was N = 2 12 = 4096.In the modified method, the size s of the overlap moving window ranged from 10~int(N/10) with an increment of 1.The step k was set to 1.The qth-order ranged from −10~10 with an increment of 0.1.
The scaling exponent τ(q), Hurst exponent H(q), singular exponent α, and spectrum f (α) of multifractal simulation data were calculated by MFDFA, OMW-MFDFA, and OMW-EEMD-MFDFA methods.Comparison of the three modes with the theoretical value is shown in Figure 2. It can be seen from Figure 2 that the results of the three methods were very asymptotic to the theoretical values.Their performance was very good, and the accuracy was also high.In general, the three methods had good results for negative q and deviated from the theoretical value for positive q.For positive q, the three methods underestimated τ(q), but the OWM-EEMD-MFDFA method was more asymptotic to theoretical values.From the multifractal spectrum in Figure 2c, the OWM-EEMD-MFDFA performed slightly better than the other two methods for small and large singular exponent α.
Water 2019, 11, x FOR PEER REVIEW 8 of 17 methods underestimated () , but the OWM-EEMD-MFDFA method was more asymptotic to theoretical values.From the multifractal spectrum in Figure 2c, the OWM-EEMD-MFDFA performed slightly better than the other two methods for small and large singular exponent .In addition, in order to quantitatively test the performance of the OMW-EEMD-MFDFA method with different lengths of simulated data, we repeated the experiments with  setting of 10, 11, 12, 13, and 14, and we computed the sum of squares for error (SSE) of the generalized Hurst exponent () and the scaled index () of the OMW-EEMD-MFDFA, OMW-MFDFA, and MFDFA methods with the theoretical values.It can be seen from Table 1 that as the  increased, the SSE of () and () to the theoretical values for the three methods became smaller and smaller, but the OMW-EEMD-MFDFA hybrid method had smaller SSE than the OMW-MFDFA and MFDFA methods.In summary, the performance of the OMW-EEMD-MFDFA method was better than the other two methods, and it could more accurately describe the multifractal characteristics of time series data.In addition, in order to quantitatively test the performance of the OMW-EEMD-MFDFA method with different lengths of simulated data, we repeated the experiments with n max setting of 10, 11, 12, 13, and 14, and we computed the sum of squares for error (SSE) of the generalized Hurst exponent H(q) and the scaled index τ(q) of the OMW-EEMD-MFDFA, OMW-MFDFA, and MFDFA methods with the theoretical values.It can be seen from Table 1 that as the n max increased, the SSE of H(q) and τ(q) to the theoretical values for the three methods became smaller and smaller, but the OMW-EEMD-MFDFA hybrid method had smaller SSE than the OMW-MFDFA and MFDFA methods.In summary, the performance of the OMW-EEMD-MFDFA method was better than the other two methods, and it could more accurately describe the multifractal characteristics of time series data.
Table 1.The sum of squares for error (SSE) of H(q) and τ(q) of the three methods with the theoretical values.

n m ax
MFDFA OMW-MFDFA OMW-EEMD-MFDFA H(q) τ(q) H(q) τ(q) H(q) τ(q) In this study, daily precipitation of the Dongting Lake Basin was analyzed using data from Mapoling and Zhijiang weather gauging stations, which were obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn).Daily precipitation for each station was measured from 1 January 2006 to 31 December 2016.Data had undergone a series of quality control by the China Meteorological Administration (CMA), including extreme values and internal consistency checks.The accuracy rate of daily precipitation data was generally more than 99%.Mapoling station is located on the lower reaches of the Xiangjiang river in Changsha city, near the Dongting Lake, while the Zhijiang station is located in the upper mountain areas of Yuanshui river in Zhijiang county.The locations of the Mapoling and Zhijiang weather gauging stations at Dongting Lake basin are shown in Figure 3.The daily precipitation of the two stations is shown in Figure 4 (Top row).

Data Preprocessing
Some scholars found that periodicity fluctuations existed in meteorological and hydrological time series data, which significantly affected the results of MFDFA analyses [2,16].So, before MFDFA analysis, annual periodicity components of the meteorological and hydrological time series data were eliminated.Cleveland et al. [49] proposed a versatile and robust method for decomposing time series, which was called Seasonal and Trend decomposition using Loess (STL), while the Loess is a method for estimating nonlinear relationships.STL can decompose time series data into seasonal (S t ), trend (T t ), and remainder (R t ) components.In this study, STL was employed to decompose precipitation into seasonal, trend, and remainder components (Figure 4).Then, the periodicity component of daily precipitation at the Mapoling and Zhijiang stations was removed from the decomposed results.The remaining data were used as the new daily precipitation time series data for MFDFA analysis.

Segmentation using the OWM Algorithm
In this paper, the OMW method was used to segment the daily precipitation time series data.Before the segmentation, some parameters needed to be set.The size s of the overlap moving window ranged from 12~int(N/12) with an increment of 1, where N was the length of the daily precipitation time series data.The step k was set to 1.We obtained N -s + 1 segments by using the OMW algorithm for each s.

Detrending using the EEMD Method
The EEMD method was used to decompose the N − s + 1 segments obtained by the above OMW algorithm for each s.Then, the local trend was determined.At last, we removed the local trend of each segment obtained by EEMD method.

Selection of the Range of q
Multifractal analysis was to study the characteristics of fractals under different scales by the q.Thus, the selection range of q played a very important role in computing multifractal characteristics.Different studies had different range of q.This meant that different ranges of q had some influence on the multifractal spectrum.While in the study, if the selected range of q was too large, the calculation time would increase exponentially with the increase of q, and the fluctuation of the multifractal spectrum would be stable; if the selected range of q was too small, the fluctuation of the multifractal spectrum would be relatively large with the increase of q, which would make it impossible to comprehend and accurately describe multifractal characteristics in the data.At present, there is no uniform standard for the selection of q.Many scholars suggest that the selected range of q is set to [−10~10], which can generally meet research needs [2,16,35].s, in this study, a q was chosen from -5 to 5 with an increment of 0.5.

Multifractal Characteristics Analysis
The fluctuation function F q (s) with different q for the daily precipitation time series data, removing the periodicity at the Mapoling and Zhijiang stations, was calculated by using the OMW-EEMD-MFDFA method.Figure 5 shows their log-log plots of the qth fluctuation F q (s) versus the time scale s and their coefficients of determination (R 2 ) for individual q values.It can be seen from the figure that the daily precipitation of the two stations had a better power-law relationship.That is to say, the daily precipitation of the multifractal characteristics of the two stations and the fluctuation behavior of its internal evolution were not random, but caused by the orderly process of long-range correlation features.Their fluctuation behavior was not random, but caused by long-range correlation (long-term persistence).the figure that the daily precipitation of the two stations had a better power-law relationship.That is to say, the daily precipitation of the multifractal characteristics of the two stations and the fluctuation behavior of its internal evolution were not random, but caused by the orderly process of long-range correlation features.Their fluctuation behavior was not random, but caused by long-range correlation (long-term persistence).Figure 6 shows the curves of ()~ and ()~, and the relationship between singularity spectrum () and singularity exponents  for the daily precipitation of the Mapoling and Zhijiang stations.It was clearly found that with the increase of q, the () nonlinearly decreased.This indicated that the daily precipitation time series data at the two stations had multifractal characteristics, which could not be characterized by monofractality.
Figure 6b shows the variation of scaling exponents () with q for the daily precipitation time Figure 6 shows the curves of H(q) ∼ q and τ(q) ∼ q, and the relationship between singularity spectrum f (α) and singularity exponents α for the daily precipitation of the Mapoling and Zhijiang stations.It was clearly found that with the increase of q, the H(q) nonlinearly decreased.This indicated that the daily precipitation time series data at the two stations had multifractal characteristics, which could not be characterized by monofractality.
Water 2019, 11, x FOR PEER REVIEW 12 of 17 data.Moreover, the nonlinear characteristics of ()~ at Zhijiang station were more significant, reflecting stronger multifractal characteristics.Figure 6c shows the singularity spectrum ()~ of daily precipitation at the Mapoling and Zhijiang stations.As shown in Figure 6c, the shape of the singularity spectrum ()~ curve had left truncation.This result demonstrated that the daily precipitation time series data at the two stations was characterized by multifractality, reflecting the internal dynamic characteristics of the daily precipitation time series data.The multifractal spectrum structures of the daily precipitation time series data at the two stations were very similar, indicating a similar evolution law.In order to measure the multifractal strength of daily precipitation time series data at the two stations, five multifractal parameters (e.g., (2), ∆(), ∆α, ∆, and B) were calculated and listed in Table 2.When q was equal to 2, the Hurst exponents (2) of the daily precipitation time series data at the Mapoling and Zhijiang stations were 0.5016 and 0.5056, respectively, indicating that daily precipitation was characterized by long-range correlation showing the long-term persistence of the previous state.The overall trend of future precipitation variability was the same as in the past.It can be seen from the table that the (2) of the daily precipitation time series data at the Zhijiang station was larger than Mapoling station, indicating a stronger long-term persistence.In addition, according to the range of the generalized Hurst exponent ∆(), the ∆() of daily precipitation time series data at the Mapoling station was larger than Zhijiang station, suggesting stronger multifractal characteristics of daily precipitation at Mapoling station.Moreover, the multifractal spectrum width ∆α was used to measure the degree of multifractality, which denoted a nonuniform clustering structure for daily precipitation.The bigger the ∆α was, the larger the variety of daily precipitation was; the stronger the degree of multifractality that was obtained, the more complex the structure of daily precipitation became.The ∆α of daily precipitation at the Mapoling and Zhijiang stations were 1.0145 and 0.9870, respectively.The multifractal spectrum width of Mapoling station was bigger than Zhijiang station.The difference between the daily maximum and minimum precipitation at the Figure 6.The curves of (a) the generalized Hurst exponents H(q) and order q; (b) the scaling exponents τ(q) and order q; (c) the singularity spectrum f (α) and singularity exponents α for the Mapoling and Zhijiang stations.
Figure 6b shows the variation of scaling exponents τ(q) with q for the daily precipitation time series at the Mapoling and Zhijiang stations.It was clear that the τ(q) increased nonlinearly with the increase of q, indicating significant multifractal characteristics of the daily precipitation time series data.Moreover, the nonlinear characteristics of τ(q) ∼ q at Zhijiang station were more significant, reflecting stronger multifractal characteristics.
Figure 6c shows the singularity spectrum f (α) ∼ α of daily precipitation at the Mapoling and Zhijiang stations.As shown in Figure 6c, the shape of the singularity spectrum f (α) ∼ α curve had left truncation.This result demonstrated that the daily precipitation time series data at the two stations was characterized by multifractality, reflecting the internal dynamic characteristics of the daily precipitation time series data.The multifractal spectrum structures of the daily precipitation time series data at the two stations were very similar, indicating a similar evolution law.
In order to measure the multifractal strength of daily precipitation time series data at the two stations, five multifractal parameters (e.g., H(2), ∆H(q), ∆α, ∆ f , and B) were calculated and listed in Table 2.When q was equal to 2, the Hurst exponents H(2) of the daily precipitation time series data at the Mapoling and Zhijiang stations were 0.5016 and 0.5056, respectively, indicating that daily precipitation was characterized by long-range correlation showing the long-term persistence of the previous state.The overall trend of future precipitation variability was the same as in the past.It can be seen from the table that the H(2) of the daily precipitation time series data at the Zhijiang station was larger than Mapoling station, indicating a stronger long-term persistence.In addition, according to the range of the generalized Hurst exponent ∆H(q), the ∆H(q) of daily precipitation time series data at the Mapoling station was larger than Zhijiang station, suggesting stronger multifractal characteristics of daily precipitation at Mapoling station.Moreover, the multifractal spectrum width ∆α was used to measure the degree of multifractality, which denoted a nonuniform clustering structure for daily precipitation.The bigger the ∆α was, the larger the variety of daily precipitation was; the stronger the degree of multifractality that was obtained, the more complex the structure of daily precipitation became.The ∆α of daily precipitation at the Mapoling and Zhijiang stations were 1.0145 and 0.9870, respectively.The multifractal spectrum width of Mapoling station was bigger than Zhijiang station.The difference between the daily maximum and minimum precipitation at the Mapoling station was relatively larger than Zhijiang station.This result indicated that the distribution of daily precipitation at Mapoling station was relatively uneven, and the multifractal strength of daily precipitation was bigger than Zhijiang station.Moreover, referring to Figure 6c, the multifractal spectrum f (α) ∼ α of the daily precipitation at the two stations was asymmetric and all had left truncation, indicating that precipitation at the two stations was insensitive to larger magnitudes of local fluctuations.All the ∆ f values were greater than zero, demonstrating that the probability of fluctuation at a high-level of daily precipitation at the two stations was dominant, and the internal variables of the system were obvious.Precipitation may have a local increasing trend.Furthermore, regarding the asymmetry parameter B, it was apparent that the values of Mapoling and Zhijiang stations were 0.4731 and 0.4695, respectively, indicating right-skewed spectra and relatively strongly weighted high fractal exponents.The results can fully prove the complexity of the precipitation evolution process at the Mapoling and Zhijiang stations.From the above experiment and analysis, we can see that there were multifractal characteristics of daily precipitation time series data at Mapoling and Zhijiang stations.However, the source that caused multifractality was unknown.Generally, many scholars believe there are two main causes that led to multifractal characteristics in time series data: (1) multifractal characteristics are caused by long-range correlations between small and large fluctuations in original time series data and (2) multifractal characteristics are caused by broad probability density function in original time series data [30,31].For the first cause, all correlations were destroyed by shuffling the time series data.Shuffling the time series data did not destroy the fluctuation characteristics and probability density function of the original time series.Thus, if multifractal characteristics only were in different correlations between small and large fluctuations, then the shuffled daily precipitation time series data showed random characteristics, and the generalized Hurst exponent was a constant value (H shu f (q) = 0.5).In addition, for the second cause, multifractal characteristics could not be removed by shuffling the original time series data.In order to obtain multifractal characteristics caused by the broad probability density function, the phases of a discrete Fourier transform coefficient of the original time series data were replaced with a set of pseudo-independent distributed, uniform (0, 2π) quantities in the surrogate method [50].The surrogate time series data did not destroy the correlation characteristics of the original time series data, it only weakened the Gaussian distribution of the original time series data [51].Therefore, if the multifractal characteristic was only caused by a broad probability density function, then the Hurst exponent H surr (q) determined by surrogate time series data would be independent of q and a constant value.If the multifractal characteristics of the daily precipitation time series data were caused by two kinds of reasons, the multifractal characteristics of shuffled and surrogate daily precipitation time series data would be weaker than those of the original time series data.
The daily precipitation time series data at the Mapoling and Zhijiang stations were shuffled and randomized, and the generalized Hurst exponent H(q), singularity exponent α, and singularity spectrum f (α) of the corresponding shuffled and surrogate daily precipitation time series data were calculated by using the OMW-EEMD-MFDFA method.The comparative results of H(q) ∼ q and f (α) ∼ α at the two stations are shown in Figure 7.The multifractal parameters of the shuffled and surrogate time series data were also calculated and listed in Table 3.It was clearly found that the generalized Hurst exponent H shu f (q) was not equal to 0.5.Moreover, the ∆H(q) values of shuffled and surrogate daily precipitation time series data at the two stations were smaller than the original time series data, demonstrating that the multifractal strength of the shuffled and surrogate time series data were weaker than the original time series data.Thus, from the above analysis, the multifractal strength of the daily precipitation time series data was caused both by long-range correlation between small and large fluctuations and a broad probability density function.Furthermore, the width of the multifractal spectrum ∆α values of the shuffled and surrogate daily precipitation time series at the two stations were also smaller than the original time series data.This result also supported the finding that shuffled and surrogate time series data had a weaker multifractality than the original time data series.As seen from Figure 7 and Table 3, the ∆H(q) and ∆α values of surrogate daily precipitation time series data were smaller than that of the shuffled time series data.This meant that the impact of the broad probability density function was greater than the long-range correlations between small and large fluctuations.Therefore, the long-range correlations between small and large fluctuations also contributed to the multifractal formation of daily precipitation time series data.The results obtained by Baranowski et al. [52] showed that the main source of multifractality was dominated by a broad probability density function.Our result was similar to theirs.
Water 2019, 11, x FOR PEER REVIEW 14 of 17 time series data, demonstrating that the multifractal strength of the shuffled and surrogate time series data were weaker than the original time series data.Thus, from the above analysis, the multifractal strength of the daily precipitation time series data was caused both by long-range correlation between small and large fluctuations and a broad probability density function.Furthermore, the width of the multifractal spectrum ∆α values of the shuffled and surrogate daily precipitation time series at the two stations were also smaller than the original time series data.This result also supported the finding that shuffled and surrogate time series data had a weaker multifractality than the original time data series.As seen from Figure 7 and Table 3, the ∆() and ∆α values of surrogate daily precipitation time series data were smaller than that of the shuffled time series data.This meant that the impact of the broad probability density function was greater than the long-range correlations between small and large fluctuations.Therefore, the long-range correlations between small and large fluctuations also contributed to the multifractal formation of daily precipitation time series data.The results obtained by Baranowski et al. [52] showed that the main source of multifractality was dominated by a broad probability density function.Our result was similar to theirs.

Conclusions
In this study, in order to overcome the weakness of the classic MFDFA method, we proposed a two-stage modified MFDFA method for multifractal analysis.In the first stage we introduced the OMW algorithm to divide daily precipitation time series data of the classic MFDFA method; in the second stage we used the EEMD method to obtain the local trend component, which replaced the polynomial local trend component in the classic MFDFA method.The modified MFDFA method was called OMW-EEMD-MFDFA.Then, the performance of the proposed OMW-EEMD-MFDFA method was tested with simulation experiments based on the p-model of the multiplicative cascading method.In all simulation experiments, the OMW-EEMD-MFDFA method performed better than classic MDFA and OMW-MFDFA methods, and it could give better results for studying multifractal characteristics of time series data.The only shortcoming of the OMW-EEMD-MFDFA was that detrending in EEMD was a time-consuming process.
Then, we applied this proposed OMW-EEMD-MFDFA method to daily precipitation time series data at Mapoling and Zhijiang stations to investigate the multifractal behaviors.The main conclusions are summarized as follows: (1) The scaling properties of the daily precipitation time series at the two stations exhibited a long-range correlation, showing a long-term persistence of the previous state.That is, the overall trend of future precipitation variations was the same as in the past.
(2) The order q-dependence of H(q) and τ(q) indicate that the daily precipitation time series data at the two stations have multifractal characteristics.In addition, according to the range of the generalized Hurst exponent ∆H(q) and the width of the singularity spectrum ∆α, the ∆H(q) and ∆α of Zhijiang station are the biggest, indicating stronger multifractal characteristics than Mapoling station.The ∆ f value of the Mapoling and Zhijiang stations is larger than 0, suggesting that future precipitation may show a local increasing trend.
(3) In order to investigate the multifractal source of the daily precipitation time series data at the two stations, the generalized Hurst exponent H(q), singularity exponent α, and singularity spectrum f (α) of the corresponding shuffled and surrogate daily precipitation time series data were calculated.We compared the generalized Hurst exponent and the multifractal strength of the original precipitation time series data with its shuffled and surrogate time series data.We found that the multifractal characteristics of the daily precipitation time series data were caused by different correlations between small and large fluctuations as well as a broad probability density function, but the broad probability density function was dominant.
In this study, we introduced an advanced and complex statistical analysis by scaling exponent and multifractal spectrums to study daily precipitation time series data.We can obtain a better understanding of statistical characteristics from multifractal parameters of precipitation time series data, which has great importance for precipitation forecasting and extreme precipitation regulation.Therefore, the results of this study can be helpful in water resource management.

Figure 1 .
Figure 1.The workflow of the overlap moving window (OWM) algorithm.

Figure 1 .
Figure 1.The workflow of the overlap moving window (OWM) algorithm.

Water 2019 ,
11, x FOR PEER REVIEW 9 of 17 was measured from 1 January 2006 to 31 December 2016.Data had undergone a series of quality control by the China Meteorological Administration (CMA), including extreme values and internal consistency checks.The accuracy rate of daily precipitation data was generally more than 99%.Mapoling station is located on the lower reaches of the Xiangjiang river in Changsha city, near the Dongting Lake, while the Zhijiang station is located in the upper mountain areas of Yuanshui river in Zhijiang county.The locations of the Mapoling and Zhijiang weather gauging stations at Dongting Lake basin are shown in Figure 3.The daily precipitation of the two stations is shown in Figure 4 (Top row).

Figure 3 .
Figure 3.The location of the Mapoling and Zhijiang weather stations.Figure 3. The location of the Mapoling and Zhijiang weather stations.

Figure 3 .
Figure 3.The location of the Mapoling and Zhijiang weather stations.Figure 3. The location of the Mapoling and Zhijiang weather stations.

Figure 3 .
Figure 3.The location of the Mapoling and Zhijiang weather stations.

Figure 4 .
Figure 4. Decomposition of daily precipitation using the Seasonal and Trend decomposition using Loess (STL) method at (a) Mapoling station and (b) Zhijiang station.From top to bottom, the plot shows the original observed daily precipitation, seasonal, trend, and the remainder component, respectively.
Water 2019, 11, x FOR PEER REVIEW 11 of 17

Figure 5 .
Figure 5.The log-log plots of qth fluctuation function  () versus time scale s of the daily precipitation time series data, removing the periodicity at the (a) Mapoling and (b) Zhijiang stations.

Figure 5 .
Figure 5.The log-log plots of qth fluctuation function F q (s) versus time scale s of the daily precipitation time series data, removing the periodicity at the (a) Mapoling and (b) Zhijiang stations.

Figure 6 .
Figure 6.The curves of (a) the generalized Hurst exponents () and order  ; (b) the scaling exponents () and order ; (c) the singularity spectrum () and singularity exponents  for the Mapoling and Zhijiang stations.

Figure 7 .
Figure 7.The ()~ curves and the singularity spectrum ()~ of the original, shuffled, and surrogate daily precipitation time series data for the (a) Mapoling station and (b) Zhijiang station.

Figure 7 .
Figure 7.The H(q) ∼ q curves and the singularity spectrum f (α) ∼ α of the original, shuffled, and surrogate daily precipitation time series data for the (a) Mapoling station and (b) Zhijiang station.

Table 1 .
The sum of squares for error (SSE) of () and () of the three methods with the theoretical values.

Table 2 .
Multifractal parameters of daily precipitation time series data at the Mapoling and Zhijiang stations.

Table 3 .
Multifractal parameters of original, shuffled, and surrogate daily precipitation time series data at the Mapoling and Zhijiang stations.

Table 3 .
Multifractal parameters of original, shuffled, and surrogate daily precipitation time series data at the Mapoling and Zhijiang stations.