# Multifractal Analysis of Hydrologic Data Using Wavelet Methods and Fluctuation Analysis

Authors: Tongzhou Zhao; Liang Wu; Dehua Li; Yiming Ding
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

We study the multifractal properties of water level with a high-frequency and massive time series using wavelet methods (estimation of Hurst exponents, multiscale diagram, and wavelet leaders for multifractal analysis (WLMF)) and multifractal detrended fluctuation analysis (MF-DFA). The dataset contains more than two million records from 10 observation sites at a northern China river. The multiscale behaviour is observed in this time series, which indicates the multifractality. This multifractality is detected via multiscale diagram. Then we focus on the multifractal analysis using MF-DFA and WLMF. The two methods give the same conclusion that at most sites the records satisfy the generalized binomial multifractal model, which is robust for different times (morning, afternoon, and evening). The variation in the detailed characteristic parameters of the multifractal model indicates that both human activities and tributaries influence the multifractality. Our work is useful for building simulation models of the water level of local rivers with many observation sites.


---

## Body

## 1. Introduction

Long-range dependence (LRD) and multifractality are the inherent characteristics of many natural phenomena. Long-range dependence (LRD) was found in the flood process by Hurst [1] when he was studying the flow of the Nile River. Fractional Brown motion (FBM), which is introduced by Mandelbrot and Van Ness [2], is the commonly used model for studying long-range dependence. Long-range dependence is characterized by the Hurst parameter, which describes the intensity of long memory phenomena. Many studies [3–6] for estimating the Hurst parameter (or Hurst exponent) H to judge the correlation effect have been made available in recent years. These studies reveal that, in many areas of applied sciences, such as climate change, stock markets, telecommunication network, and river flow, the phenomena appear as long-range behaviour, and the estimation of Hurst has been widely used in decisions and predictions.The concept of multifractality was originally introduced by Mandelbrot [7] in the context of turbulence. The formalism of multifractality was first introduced using box-counting algorithms for the complexity of measures [8–10] and further generalized to functions using the structure function method [11]. Muzy et al. [12] first developed the multifractal formalism using wavelets in the discussion of turbulence data. Schertzer and Lovejoy [13] concluded that a multifractal field could be regarded as the result of a multiplicative cascade process. Until now, the multifractality and its formalism have been applied in hydrology and many other contexts [14–20].Previous studies on the LRD and multifractality of hydrology (see, e.g., [14–20]) were based on relatively large sampling periods, such as 1 day, 1 month, or even 1 year, rarely on high-frequency hydrological data. Venugopal et al. [21–23] studied the multifractality of high-resolution temporal rainfall data using wavelet methods. We focus on studying the properties of high-frequency hydrological data and would like to further explore the typical LRD phenomena and multifractal behaviours in hydrology. Our records are the high-frequency water level records (every 6 minutes) of one river (a branch of the Haihe River in China) at its 10 observation sites from April 2011 to September 2013. We have not yet seen studies conducted for such short sampling intervals and such a wide range of observation sites (along the river from upstream to downstream).It is known that the hydrologic data are often affected by trends or other nonstationarities, for example, due to the seasonal cycle or a change in climate, which may lead to an overestimation of the Hurst exponent. Furthermore, under the influence of the trends, the uncorrelated data may behave as long-term correlated data. Thus, the methods applied in this paper are based on wavelet analysis (estimation of Hurst exponents, multiscale diagram, and wavelet leaders for multifractal analysis (WLMF)) and multifractal detrended fluctuation analysis (MF-DFA). The wavelet method can eliminate some trends as a result of the vanishing moment property, and MF-DFA can avoid spurious detection of correlations that are artifacts of nonstationarity in the records [24]. The wavelet-based estimator of the Hurst exponent was proposed by Veitch and Abry [25] and has been improved in recent years [26, 27]. This estimate was further extended to the situation of random fields [28, 29]. In addition to the estimation, the wavelet-based method can also be used to obtain a multiscale diagram for detecting the multifractal behaviour [8, 9, 26, 30]. The multifractality has been studied in several methods since it was proposed; the primary and advanced methods are wavelet analysis [10, 31, 32] and multifractal detrended fluctuation analysis [24, 33]. The MF-DFA is based on the standard fluctuation analysis (FA) and is a multifractal generalization of the detrended fluctuation analysis (DFA) [24, 33]. Compared with FA, the method can be used to conduct nonstationary data for its detrending function. In recent years, the well-established MF-DFA has been widely used in hydrologic data (see Koscielny-Bunde et al. [14, 16]). Lashermes et al. introduced a new multifractal formalism based on wavelet leaders [34], which allows the multifractal spectrum to be obtained over its entire range. Since then, many studies on multifractality and its applications have been continuously conducted [35]. The analysis method, namely, wavelet leaders multifractal analysis (WLMF), is a discrete version of wavelet transform modulus maxima (WTMM) method [36, 37]. The DFA and WTMM methods are compared in many literatures (see, e.g., [16, 24, 38]). Comparing with WTMM, DFA is easy to implement and has smaller MSE than WTMM for very short time series. But the DFA is lack of deep theoretical concepts underlying the multifractal formalism [38].In this paper, we first apply the wavelet-based method to estimate the Hurst exponent with different sampling and wavelet scales and observe both the LRD phenomenonH>0.5 and multiscale behaviour in the high-frequency water level records. We think that the multiscale behaviour may induce multifractality. This multifractality is detected via multiscale diagram. Then we apply MF-DFA and WLMF to the records to study the multifractality in detail. The two methods give the same conclusion that the water level records at most sites (except some abnormal sites, where the water levels are seriously affected by human activities) satisfy the generalized binomial multifractal model, which is robust for different times (morning, afternoon, and evening). Another interesting result is that, by using MF-DFA, the values of a+b vary from approximately 1 to less than 0.9; we conclude that this phenomenon occurs because there are four tributaries that flow into the river, changing the multifractality of the river. As a technical improvement in the estimation of the Hurst parameter H, we show that the sampling intervals and choice of scales of the discrete wavelet transformation should be matched. Our work is useful for those who want to build simulation models of the water level of local rivers with many observation sites, because the LRD and multifractal properties of the simulation data can be checked by the methods used in our paper.The logical structure and organization of the paper are described in Figure1.Figure 1
Data processing flow chart.The remaining parts of the paper are organized as follows. We introduce the data source and data extraction method for analysis in Section2. The methods employed in this paper are introduced briefly in Section 3. The detailed analysis and the results for the records are presented in Section 4. In Section 5, we compare the results obtained from Section 4 with the previous results. We conclude our work in Section 6.


## 2. Data

The records that are collected in real time from one river (a branch of the Haihe River in China) represent the water level. The river is located in northern China, which has a significant continental monsoon climate. It is cold and dry in the winter and warm and humid in the summer. From upstream to downstream, 10 artificial canals have been built to control the water level. Next to these canals, 10 water level observation sites (Site 1–Site 10) have been set up along the river (see Figure2). Four tributaries are located downstream of Site 7 that flow into this river. The data were obtained from these 10 sites and collected every 6 minutes over the past two and a half years from April 2011 to September 2013. Table 1 shows the sampling information at each site. Figure 3 shows the water level records at Site 1 denoted by Xt and the cumulative sum of Xt-X¯. X¯ is the mean of Xt. The MF-DFA and WLMF methods require that the records should be a measure and its cumulative sum should be a function. The cumulative sum in Figure 3 displays a function form. This confirms the use of the MF-DFA and WLMF methods. Partial records for validation are available at GitHub: https://github.com/tongzhouzhao/water-level-records.Table 1
Sampling information.
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Site 7
Site 8
Site 9
Site 10
Number of records
182855
152307
189967
190959
188074
185557
190840
190445
185822
111600
Sampling period (days)
904
837
904
904
904
904
904
904
903
698
Sampling frequency
6 min
6 min
6 min
6 min
6 min
6 min
6 min
6 min
6 min
6 min
Note. The records numbers are not equal in the same sampling period because of the missing data.Figure 2
Site 1–Site 10 locations.Figure 3
The water level recordsXt at Site 1 (a) and the cumulative sum of Xt-X¯ (b).
(a)
(b)The water level of the river is always influenced by the weather conditions and human activities. Every year, during the traditional rainy season, the upstream canals are always filled with much rainwater. During the dry season, the upstream canals discharge the stored water. Both the storage and discharge of water should be adjusted at any time according to the growing condition of crops along the bank of the river and the actual demand of human activities. According to our survey, water storage and discharge at Site 2 and Site 7 are conducted frequently for crop growth and human activities.The following three preprocessing steps were conducted before the multifractal analyses: (1) removing the outliers and (2) extracting the appropriate data. These preprocessing steps were conducted to exclude outliers, noise, and artifacts in the data. Firstly, we deleted the outliers caused by the complex environment where the data were acquired. For instance, at 12:42:00 7/18/2013, the water level measured at Site 8 was −1.89, which is obviously an outlier. Secondly, we adjust the artifact data records. The artificial control of the sluice may cause the water level to rise or fall in a short period of time. For example, at 13:12:00 7/5/2013, the record (Site 8) dropped from 27.3 to 0.28. The phenomenon of large fluctuations in a short time may be misdiagnosed, as there is noise in the dataset. Thirdly, we fill the missing data appropriately. There is a small portion of data missing during times when the data-collecting equipment did not work well. For instance, at Site 5, we received one day’s last record at 16:42:00 4/27/2011, which should have been followed by another record six minutes later at 16:48:00 on the same day. However, the next record time in our dataset is 8:48:00 4/29/2011, which means that almost one-day records are missing. A situation like this can occur at other sites. The missing record is directly relevant to the nearby records. For avoiding the effect of time in one day, we use the records at the same time of adjacent days to fill the missing data. The used method is linear interpolation. For example, at Site 5, the record at 16:18:00 4/28/2011 is missing. According to linear interpolation method, it can be filled with the mean of the records at 16:18:00 4/27/2011 and 16:18:00 4/29/2011. The preprocessing may ensure a uniformly spaced sampling interval.


## 3. Methods

### 3.1. Wavelet-Based Estimation of the Hurst Exponent

Veitch et al. [25, 26] have proposed the wavelet-based method for many years. Recently, they improved it for non-Gaussian long-range dependent processes [27]. We use this method to estimate the Hurst exponent. The wavelet-based method can detrend the records with polynomial trend because of the property of vanishing moments of wavelet, can avoid the spurious result caused by trend, and is suitable for water level records.LetXt,t∈0,n, denote a second-order stationary process with a spectrum labelled ΓXν, which satisfies(1)ΓXν~CXν-2H-1,ν⟶0,0.5<H<1.The process Xt is said to be long-range dependent.Letψt be a wavelet function with a vanishing moment N. The wavelet coefficients of Xt are(2)dXj,k=∫Rψj,ktXtdt,whereψj,kt=2-j/2ψ2-jt-k.Because of the stationary property ofXt, we have(3)EdXj,k2=∫RΓXνΨ2jνdν,and considering the property of the spectrum density, (3) becomes(4)EdXj,k2~CX2j2H-1,2j⟶∞.As a result,(5)log2⁡EdXj,k2~j2H-1+log2⁡CX,2j⟶∞.Put(6)Sj≔1nj∑dXj,k2,where nj is the number of dXj,k actually available at scale 2j. EdXj,k2 can be estimated by Sj. Then, the estimation of the Hurst parameter can be realized with a linear regression.Because of the different variances oflog2⁡Sj at different values of j, the weighted least squares for this regression model are needed. The weight for this regression is the reciprocal of the variance of log2⁡Sj. Moreover,(7)varlog2⁡Sj≈log2⁡e2Cj,where Cj=vardX2j,·/njEdX2j,·2. The Cj term can be estimated for each octave by using the sample moment estimators of the fourth and second moments of dXj,·.Abry et al. [27] studied the statistical performance of this estimator for non-Gaussian long-range dependent processes obtained from point transformations of Gaussian processes. Their results suggest that the vanishing moment should satisfy N≥2. Furthermore, the linear regression should be performed over scales of 2j, ranging from the finest to the coarsest available. To obtain the estimators at different scales, we do not comply with this suggestion in our experiments. It is also shown that the reduction of the dependence range of wavelet coefficients from long to short by increasing N above 2, as observed for Gaussian processes, does not hold for non-Gaussian processes. This may lead to the fact that increasing N above 2 does not result in any decrease in the variance of the estimator.


### 3.2. Wavelet-Based Multiscale Diagram

Wavelet procedures provide us with a way to estimate the Hurst parameter. Sometimes, the Hurst index depends on timet (i.e., multiple exponents), such as multifractional process and multifractal process. For the multifractal process, its sample path has local regularity:(8)Xω,t+τ-Xω,t≈τHω,t,τ⟶0,where ω denotes an element of the probability space [26]. Abry et al. [26] introduced a wavelet-based multiscale diagram for detecting multifractality of process. This method abandons the idea of estimating H that changes with t and studies the process statistically by Legendre multifractal spectrum, which can be obtained by calculating the structure function (or partition functions) from the sample path:(9)Sqτ=∫Xω,t+τ-Xω,tqdt,and for the multifractal case, Sqτ≈τζq,τ→0,q∈R [26, 39]. Moreover, the Legendre transform of the function ζq is the Legendre multifractal spectrum.Wavelet processes are an effective tool for obtainingζq because of the fact that the increments of the sample path can be seen as simple examples of the wavelet coefficients. The wavelet coefficients are used to replace the increments in the partition function [9, 26]. Furthermore, the wavelet-based partition function exhibits power-law behaviour for small scales [8, 26, 30, 39]:(10)∫TXa,tqdt≈aζq+q/2,a⟶0,whereTXa,t=∫1aXuψu-tadu.According to the power-law, the statistics Sj,q≔1/nj∑kdXj,kq≈2jζq+q/2,q∈R, for small j. So ζq can be estimated using linear regression [26].Ifζq=qH (linear form), the process is often called “monofractal,” whose multifractal spectrum is determined by one H, such as self-similar process (e.g., FBM) and LRD process. Conversely, if ζq is not linear, the process exhibits multiscale property. The multifractality of a process can be detected by measuring the deviation of ζq from the linear form [26]. Let ζ^q be the estimation of ζq. The graph of ζ^q against q, together with the confidence intervals related to ζ^q, is designated as “multiscale diagram (MD)” [26]. The confidence intervals are obtained according to Delbeke and Van Assche [40]. The detection of multifractality can be solved by examining whether there is a lack of alignment in the multiscale diagram.


### 3.3. Multifractal Detrended Fluctuation Analysis

The wavelet-based multiscale diagram can detect whether the processXt is multifractal. We study the multifractality by using MF-DFA [24], which is based on the standard fluctuation analysis (FA).In MF-DFA, one investigates the series:(11)Yi≔∑k=1iXk-X¯,i=1,…,L,where L is the length of Xt and X¯ is the mean of Xt.LetNs≔intL/s. Yt is divided into Ns nonoverlapping segments with equal length s. The fluctuation F2v,s is calculated using(12)F2v,s=1s∑i=1sYv-1s+i-yvi2,v=1,…,Ns,where yvi is the best polynomial fit of Yv-1s+i in each of the Ns segments v.In FA, the fluctuationF2v,s=Yvs-Yv-1s2, the trend in the original series is not considered. In MF-DFA, the time series is detrended by subtracting yvi. The m-order polynomial used in the fitting procedure can eliminate the polynomial trend of order m-1 in the original series. This ensures the fluctuation analysis available for data affected by trends or other nonstationarities.For a multifractal description of records, one must consider all momentsFqs:(13)Fqs≔1Ns∑v=1NsF2v,sq/21/q,v=1,…,Ns.If Xt is long-range dependent, Fqs satisfies [24](14)Fqs~shq,for  larges.For very larges, that is, s>N/4, Fqs becomes inaccurate due to the statistical errors produced by the small number of segments Ns. So we will concentrate on the scales s lower than N/4. The minimum scale chosen is equal to 25 to avoid spurious results produced by s<30.h
q can be obtained by analyzing log-log plots Fqs versus s. If q=2, h2 is the fluctuation exponent H that is related to the correlation exponent γ by H=1-γ/2 [24]. For a stationary time series, h2 is the well-known Hurst exponent. The function hq is called the generalized Hurst exponent [24].It is known that [14, 24], for q>0, hq shows the scaling behaviour of the segments with large fluctuation, and, for q<0, hq shows that with small fluctuation.One can also get the Renyi scaling exponentτq by calculating τq=qhq-1. The way to describe the multifractality is the singularity spectrum fα (or multifractal spectrum), that is, the Legendre transform of τq [14, 16, 24]:(15)α=τ′q,fα=qα-τq,where α is singularity strength (or Hölder exponent). The strength of the multifractality of records can be characterized by Δα=αmax-αmin.Ifhq satisfies the following equation,(16)hq=1q-ln⁡aq+bqqln⁡2,which can be obtained from a generalized binomial multifractal model [16], τq=-ln⁡aq+bq/ln⁡2, Δα=h-∞-h∞. In this case,(17)Δα=ln⁡b-ln⁡aln⁡2.The three parametersa,b, and Δα can then be regarded as multifractal finger prints for records [14].


### 3.4. Wavelet Leaders for Multifractal Analysis

For comparison with MF-DFA, we introduce the method of wavelet leaders for multifractal analysis (WLMF) proposed by Wendt et al. (see [41]).Assume thatψt is a wavelet function with compact time support and N vanishing moment. Let λj,k=k2j,k+12j be dyadic intervals and let 3λj,k denote the intervals 3λj,k≔λj,k-1∪λj,k∪λj,k+1. The wavelet leaders are defined as(18)LXj,k≔supk′2j′,k′+12j′⊂3λj,k⁡dXj′,k′.The structure function is calculated by the averages of LXj,k:(19)SqLj=1nj∑k=1njLXj,kq.It can been shown [41] that, under a mild uniform Hölder regularity condition on Xt,(20)SqLj~Fq2jζq,2j⟶0.So the estimation for ζq can be realized by a linear regression in j versus log2⁡SqLj:(21)ζ^q=∑j=j1j2ωjlog2⁡SqLj.The weightsωj satisfy ∑j1j2jωj≡1 and ∑j1j2ωj≡0.To estimate the multifractal spectrum, Wendt et al. [41] suggest a parametric formulation that avoids the computation of the Legendre transform:(22)f^q=∑j1j2ωjULj,q,h^q=∑j1j2ωjVLj,q,where ULj,q=∑k=1njRXqj,klog2⁡RXqj,k+log2⁡nj, VLj,q=∑k=1njRXqj,klog2⁡LXj,k, and RXqj,k=LXj,kq/∑k=1njLXj,kq.We emphasize that the wavelet transform withN vanishing moments can eliminate the polynomial trend of order N-1 in the original series.


## 3.1. Wavelet-Based Estimation of the Hurst Exponent

Veitch et al. [25, 26] have proposed the wavelet-based method for many years. Recently, they improved it for non-Gaussian long-range dependent processes [27]. We use this method to estimate the Hurst exponent. The wavelet-based method can detrend the records with polynomial trend because of the property of vanishing moments of wavelet, can avoid the spurious result caused by trend, and is suitable for water level records.LetXt,t∈0,n, denote a second-order stationary process with a spectrum labelled ΓXν, which satisfies(1)ΓXν~CXν-2H-1,ν⟶0,0.5<H<1.The process Xt is said to be long-range dependent.Letψt be a wavelet function with a vanishing moment N. The wavelet coefficients of Xt are(2)dXj,k=∫Rψj,ktXtdt,whereψj,kt=2-j/2ψ2-jt-k.Because of the stationary property ofXt, we have(3)EdXj,k2=∫RΓXνΨ2jνdν,and considering the property of the spectrum density, (3) becomes(4)EdXj,k2~CX2j2H-1,2j⟶∞.As a result,(5)log2⁡EdXj,k2~j2H-1+log2⁡CX,2j⟶∞.Put(6)Sj≔1nj∑dXj,k2,where nj is the number of dXj,k actually available at scale 2j. EdXj,k2 can be estimated by Sj. Then, the estimation of the Hurst parameter can be realized with a linear regression.Because of the different variances oflog2⁡Sj at different values of j, the weighted least squares for this regression model are needed. The weight for this regression is the reciprocal of the variance of log2⁡Sj. Moreover,(7)varlog2⁡Sj≈log2⁡e2Cj,where Cj=vardX2j,·/njEdX2j,·2. The Cj term can be estimated for each octave by using the sample moment estimators of the fourth and second moments of dXj,·.Abry et al. [27] studied the statistical performance of this estimator for non-Gaussian long-range dependent processes obtained from point transformations of Gaussian processes. Their results suggest that the vanishing moment should satisfy N≥2. Furthermore, the linear regression should be performed over scales of 2j, ranging from the finest to the coarsest available. To obtain the estimators at different scales, we do not comply with this suggestion in our experiments. It is also shown that the reduction of the dependence range of wavelet coefficients from long to short by increasing N above 2, as observed for Gaussian processes, does not hold for non-Gaussian processes. This may lead to the fact that increasing N above 2 does not result in any decrease in the variance of the estimator.


## 3.2. Wavelet-Based Multiscale Diagram

Wavelet procedures provide us with a way to estimate the Hurst parameter. Sometimes, the Hurst index depends on timet (i.e., multiple exponents), such as multifractional process and multifractal process. For the multifractal process, its sample path has local regularity:(8)Xω,t+τ-Xω,t≈τHω,t,τ⟶0,where ω denotes an element of the probability space [26]. Abry et al. [26] introduced a wavelet-based multiscale diagram for detecting multifractality of process. This method abandons the idea of estimating H that changes with t and studies the process statistically by Legendre multifractal spectrum, which can be obtained by calculating the structure function (or partition functions) from the sample path:(9)Sqτ=∫Xω,t+τ-Xω,tqdt,and for the multifractal case, Sqτ≈τζq,τ→0,q∈R [26, 39]. Moreover, the Legendre transform of the function ζq is the Legendre multifractal spectrum.Wavelet processes are an effective tool for obtainingζq because of the fact that the increments of the sample path can be seen as simple examples of the wavelet coefficients. The wavelet coefficients are used to replace the increments in the partition function [9, 26]. Furthermore, the wavelet-based partition function exhibits power-law behaviour for small scales [8, 26, 30, 39]:(10)∫TXa,tqdt≈aζq+q/2,a⟶0,whereTXa,t=∫1aXuψu-tadu.According to the power-law, the statistics Sj,q≔1/nj∑kdXj,kq≈2jζq+q/2,q∈R, for small j. So ζq can be estimated using linear regression [26].Ifζq=qH (linear form), the process is often called “monofractal,” whose multifractal spectrum is determined by one H, such as self-similar process (e.g., FBM) and LRD process. Conversely, if ζq is not linear, the process exhibits multiscale property. The multifractality of a process can be detected by measuring the deviation of ζq from the linear form [26]. Let ζ^q be the estimation of ζq. The graph of ζ^q against q, together with the confidence intervals related to ζ^q, is designated as “multiscale diagram (MD)” [26]. The confidence intervals are obtained according to Delbeke and Van Assche [40]. The detection of multifractality can be solved by examining whether there is a lack of alignment in the multiscale diagram.


## 3.3. Multifractal Detrended Fluctuation Analysis

The wavelet-based multiscale diagram can detect whether the processXt is multifractal. We study the multifractality by using MF-DFA [24], which is based on the standard fluctuation analysis (FA).In MF-DFA, one investigates the series:(11)Yi≔∑k=1iXk-X¯,i=1,…,L,where L is the length of Xt and X¯ is the mean of Xt.LetNs≔intL/s. Yt is divided into Ns nonoverlapping segments with equal length s. The fluctuation F2v,s is calculated using(12)F2v,s=1s∑i=1sYv-1s+i-yvi2,v=1,…,Ns,where yvi is the best polynomial fit of Yv-1s+i in each of the Ns segments v.In FA, the fluctuationF2v,s=Yvs-Yv-1s2, the trend in the original series is not considered. In MF-DFA, the time series is detrended by subtracting yvi. The m-order polynomial used in the fitting procedure can eliminate the polynomial trend of order m-1 in the original series. This ensures the fluctuation analysis available for data affected by trends or other nonstationarities.For a multifractal description of records, one must consider all momentsFqs:(13)Fqs≔1Ns∑v=1NsF2v,sq/21/q,v=1,…,Ns.If Xt is long-range dependent, Fqs satisfies [24](14)Fqs~shq,for  larges.For very larges, that is, s>N/4, Fqs becomes inaccurate due to the statistical errors produced by the small number of segments Ns. So we will concentrate on the scales s lower than N/4. The minimum scale chosen is equal to 25 to avoid spurious results produced by s<30.h
q can be obtained by analyzing log-log plots Fqs versus s. If q=2, h2 is the fluctuation exponent H that is related to the correlation exponent γ by H=1-γ/2 [24]. For a stationary time series, h2 is the well-known Hurst exponent. The function hq is called the generalized Hurst exponent [24].It is known that [14, 24], for q>0, hq shows the scaling behaviour of the segments with large fluctuation, and, for q<0, hq shows that with small fluctuation.One can also get the Renyi scaling exponentτq by calculating τq=qhq-1. The way to describe the multifractality is the singularity spectrum fα (or multifractal spectrum), that is, the Legendre transform of τq [14, 16, 24]:(15)α=τ′q,fα=qα-τq,where α is singularity strength (or Hölder exponent). The strength of the multifractality of records can be characterized by Δα=αmax-αmin.Ifhq satisfies the following equation,(16)hq=1q-ln⁡aq+bqqln⁡2,which can be obtained from a generalized binomial multifractal model [16], τq=-ln⁡aq+bq/ln⁡2, Δα=h-∞-h∞. In this case,(17)Δα=ln⁡b-ln⁡aln⁡2.The three parametersa,b, and Δα can then be regarded as multifractal finger prints for records [14].


## 3.4. Wavelet Leaders for Multifractal Analysis

For comparison with MF-DFA, we introduce the method of wavelet leaders for multifractal analysis (WLMF) proposed by Wendt et al. (see [41]).Assume thatψt is a wavelet function with compact time support and N vanishing moment. Let λj,k=k2j,k+12j be dyadic intervals and let 3λj,k denote the intervals 3λj,k≔λj,k-1∪λj,k∪λj,k+1. The wavelet leaders are defined as(18)LXj,k≔supk′2j′,k′+12j′⊂3λj,k⁡dXj′,k′.The structure function is calculated by the averages of LXj,k:(19)SqLj=1nj∑k=1njLXj,kq.It can been shown [41] that, under a mild uniform Hölder regularity condition on Xt,(20)SqLj~Fq2jζq,2j⟶0.So the estimation for ζq can be realized by a linear regression in j versus log2⁡SqLj:(21)ζ^q=∑j=j1j2ωjlog2⁡SqLj.The weightsωj satisfy ∑j1j2jωj≡1 and ∑j1j2ωj≡0.To estimate the multifractal spectrum, Wendt et al. [41] suggest a parametric formulation that avoids the computation of the Legendre transform:(22)f^q=∑j1j2ωjULj,q,h^q=∑j1j2ωjVLj,q,where ULj,q=∑k=1njRXqj,klog2⁡RXqj,k+log2⁡nj, VLj,q=∑k=1njRXqj,klog2⁡LXj,k, and RXqj,k=LXj,kq/∑k=1njLXj,kq.We emphasize that the wavelet transform withN vanishing moments can eliminate the polynomial trend of order N-1 in the original series.


## 4. Computational Results

### 4.1. Wavelet-Based Results

This section describes the wavelet-based estimation of Hurst parameters for the water level records of the 10 sites. We estimate the Hurst exponents at different scales and different sampling interval times. The results indicate that there is an interesting connection between the wavelet scale and the sampling period. For comparison, we repeat the procedure for a synthetic exponential long-range dependent process (ELRD) with length of 180000, which is close to the actual sampling data length. The vanishing moment of the used Daubechies wavelet isN=3 (db3). The simulation results for ELRD are listed in Tables 2 and 3. The values in brackets are errors bars. The errors bars are based on the weighted linear least square fit which is used to estimate H and reveal 95% confidence interval.Table 2
Simulation data analysis forH (H = 0.6, data length = 180000).
Interval
j
1~5
2~6
3~7
4~8
5~9
1
0.61 (0.009)
0.61 (0.013)
0.60 (0.018)
0.59 (0.025)
0.59 (0.036)
2
0.61 (0.013)
0.61 (0.018)
0.60 (0.025)
0.59 (0.036)
0.60 (0.051)
4
0.61 (0.018)
0.61 (0.025)
0.60 (0.036)
0.59 (0.051)
0.61 (0.071)
8
0.60 (0.025)
0.60 (0.036)
0.60 (0.051)
0.61 (0.071)
0.62 (0.106)Table 3
Simulation data analysis forH (H = 0.8, data length = 180000).
Interval
j
1~5
2~6
3~7
4~8
5~9
1
0.83 (0.009)
0.82 (0.013)
0.81 (0.018)
0.80 (0.025)
0.80 (0.034)
2
0.83 (0.013)
0.82 (0.018)
0.81 (0.025)
0.80 (0.034)
0.82 (0.049)
4
0.83 (0.018)
0.82 (0.025)
0.81 (0.034)
0.81 (0.049)
0.80 (0.072)
8
0.82 (0.025)
0.82 (0.034)
0.82 (0.049)
0.80 (0.072)
0.79 (0.100)The simulation results show that the simulations for bothH=0.6 and H=0.8 are stable under changes in the interval length and estimation scale. One can see that the wavelet-based estimation of ELRD for the Hurst exponent is valid and reliable, and the selection of wavelet scale and sampling cycle is independent of the estimation.Tables4–6 display the results of the practical water level. There exists an expression aij≈ai+1j-1 (aij is the element of the tables). For example, in Table 4,H = 0.82 (12 min, scales (1~5)) andH = 0.84 (6 min, scales (2~6)) are approximately equal;H = 0.70 (24 min, scales (1~5)),H = 0.68 (12 min, scales (2~6)) andH = 0.68 (6 min, scales (3~7)) are almost equal, and so forth. The results represent the consistent relationship between the wavelet scale and the sampling period. We present a mathematical proof in the Appendix.Table 4
Multiscales and multisampling forH of Site 1 (data length = 182855).
Interval
j
1~5
2~6
3~7
4~8
6 min
0.84 (0.062)
0.84 (0.047)
0.68 (0.046)
0.53 (0.054)
12 min
0.82 (0.067)
0.68 (0.054)
0.53 (0.054)
—
24 min
0.70 (0.046)
0.55 (0.056)
—
—
48 min
0.57 (0.053)
—
—
—Table 5
Multiscales and multisampling forH of Site 2 (data length = 152307).
Interval
j
1~5
2~6
3~7
4~8
6 min
0.80 (0.101)
0.86 (0.086)
0.81 (0.077)
0.61 (0.086)
12 min
0.89 (0.082)
0.80 (0.079)
0.61 (0.086)
0.50 (0.111)
24 min
0.81 (0.079)
0.63 (0.087)
0.50 (0.112)
—
48 min
0.65 (0.085)
0.52 (0.116)
—
—Table 6
Multiscales and multisampling forH of Site 3 (data length = 189967).
Interval
j
1~5
2~6
3~7
4~8
6 min
0.79 (0.047)
0.97 (0.041)
0.87 (0.044)
0.69 (0.048)
12 min
0.95 (0.046)
0.89 (0.041)
0.69 (0.047)
0.50 (0.059)
24 min
0.91 (0.044)
0.71 (0.048)
0.50 (0.059)
—
48 min
0.74 (0.046)
0.50 (0.060)
—
—We can also see that the estimations for the Hurst exponent are not stable under changes in the interval length and estimation scale, while the estimations along the secondary diagonal direction of the tables are still stable. We conclude that the massive hydrologic data have a multiscale consistent property, while the simulation data do not. This indicates that there may exist multifractality in the records.


### 4.2. Multiscale Diagram for Hydrological Data

In this section, we use the multiscale diagram to test the multifractality that may exist in the records. The multiscale diagram is the graph ofζ^q against q, together with the confidence intervals related to ζ^q. The multifractality can be detected by measuring the deviation of ζ^q from the linear form [26].We use a linear multiscale diagram to examine the alignment in the multiscale diagram, in whichHq=ζq/q is drawn against q. The multifractality can be determined by examining whether there is a horizontal alignment in the linear multiscale diagram. The horizontal alignment indicates that the Hurst parameter H can control all of the scaling property. The value of H can be received by estimating the level. The linear multiscale diagram and multiscale diagram are statistically equivalent.To determine the alignment of hydrological data, we impose the synthesized FBM (H=0.8) known as the “linear form” for comparison. The regressions for estimating ζq for each of the 10 different q are performed for small scales: j1,j2=2,6 according to the power-law behaviour of Sj,q discussed in Section 3.2. For Site 1–Site 10, as well as the FBM, the multiscale diagrams are shown in Figure 4, and the values of ζq for all 10 sites are also shown in Table 7 (the sampling interval is 6 minutes). We choose db3 as the wavelet of the wavelet transform, which is orthogonal to quadratic trends.Table 7
ζ
(
q
) for hydrological data of 10 sites; ζ(q) is estimated with a sampling interval of 6 min and (j1,j2)=(2,6).
q
0.5
1
1.5
2
2.5
3
4
5
6
8
Site 1
0.52
0.86
1.13
1.33
1.49
1.61
1.80
1.97
2.15
2.52
Errors bars
0.008
0.015
0.025
0.041
0.064
0.095
0.170
0.250
0.319
0.423
Site 2
0.51
0.97
1.34
1.56
1.68
1.76
1.87
1.97
2.09
2.39
Errors bars
0.010
0.023
0.046
0.083
0.123
0.160
0.233
0.308
0.383
0.526
Site 3
0.52
1.00
1.39
1.68
1.85
1.96
2.07
2.12
2.16
2.26
Errors bars
0.007
0.014
0.024
0.045
0.077
0.117
0.201
0.282
0.357
0.471
Site 4
0.59
1.16
1.61
1.70
1.59
1.47
1.27
1.15
1.07
1.00
Errors bars
0.007
0.017
0.064
0.156
0.222
0.272
0.362
0.441
0.511
0.618
Site 5
0.55
0.94
1.10
0.98
0.76
0.49
−0.06
−0.61
−1.13
−2.10
Errors bars
0.013
0.031
0.086
0.161
0.228
0.288
0.394
0.481
0.552
0.655
Site 6
0.42
0.80
1.12
1.37
1.58
1.73
1.91
2.02
2.10
2.26
Errors bars
0.008
0.015
0.023
0.037
0.060
0.090
0.161
0.227
0.285
0.386
Site 7
0.48
0.77
1.02
1.23
1.42
1.58
1.86
2.08
2.25
2.56
Errors bars
0.016
0.031
0.050
0.075
0.110
0.157
0.282
0.420
0.537
0.684
Site 8
0.58
1.07
1.52
1.77
1.82
1.80
1.74
1.71
1.71
1.76
Errors bars
0.015
0.027
0.067
0.188
0.334
0.437
0.542
0.589
0.615
0.653
Site 9
0.46
0.88
1.25
1.51
1.66
1.76
1.98
2.22
2.46
2.94
Errors bars
0.007
0.015
0.035
0.079
0.140
0.202
0.304
0.373
0.423
0.495
Site 10
0.59
1.12
1.59
1.86
1.95
1.98
2.03
2.11
2.22
2.50
Errors bars
0.012
0.024
0.057
0.135
0.220
0.286
0.383
0.455
0.513
0.599
Note. The errors bars are based on the weighted linear least square fit which is used to estimate ζ(q) and reveal 95% confidence interval.Figure 4
Testing for multiscaling. (a) Multiscale diagram. Compared with the alignment for FBM, hydrological data with no alignment suggests multiscale behaviour. (b) Linear multiscale diagram. The FBM series exhibits horizontal alignment. The errors bars are based on the weighted linear least square fit which is used to estimateζq and reveal 95% confidence interval.
(a)
(b)In Figure4, the multiscale diagram in (a) shows a lack of alignment for the hydrologic data compared to the clear alignment of FBM for its self-similarity, which indicates that it has multifractal behaviour. Thus, the multifractal model may be relevant to the hydrological data. The linear multiscale diagram in (b) confirms this conclusion. The Hurst parameter H of the synthesized FBM series can be estimated by the horizontal alignment, and its estimator is approximately 0.8. This is not the case for the hydrologic data. In summary, the hydrological data have multifractal behaviour and do not comply with the equation ζq=Hq. The errors bars in Figure 4 are based on the weighted linear least square fit, which is used to estimate ζq. They are calculated as 1.96×Stdζq. Under the assumption of normality, this reveals a 95% confidence interval. The detailed values of ζq for other sites can be seen in Table 7.


### 4.3. Multifractal Analysis and Results

In this section, we analyze the daily water level records at the 10 observation sites using MF-DFA analysis for the multifractal behaviour. In order to ensure that there are enough fluctuations in each segment so that MF-DFA can be applied, the high-frequency collected daily records are used. For sufficient samples, we get 100 sets of daily records randomly at different times of one day. The number of records of each sets is about six hundred. We also apply wavelet leaders for multifractal analysis (WLMF) based on the daily water level records for a comparison with MF-DFA.Statistical Tests. It is well known that finite-size effects and lack of statistical convergence will introduce biases on the estimate of hq when applying the method of MF-DFA and WLMF. Also the robustness of the results obtained by different orders of polynomials in MF-DFA and different vanishing moments of wavelet in WLMF should be tested. In this part, we perform some statistical tests for these issues.Following the test method reported by Khalil et al. [42], we consider the problem of statistical convergence of hq through the statistical convergence of Fqqs and SqLj in the two methods. Recall that Fqqs and SqLj can be transformed into a continuous integral over M,(23)SqLj=1nj∑k=1njLXj,kq=∫MqPjMdM,Fqqs=1Ns∑v=1NsF2v,sq/2=∫MqPsMdM.M means the modulus of LXj,k or F2v,s1/2. PjM means the probability density function (pdf) of LXj,k. PsM means the probability density function (pdf) of F2v,s1/2. The well-defined shape of the pdf’s plots indicates a proper convergence [42, 43]. Figures 5(b), 5(f), 6(b), and 6(f) show PjM and PsM computed at three different scales at Site 4 and Site 7. It can be seen that there are not enough statistics to give a well-defined shape of the pdf. However, the statistical convergence of Fqqs and SqLj does not need so many statistics that can give a well-defined shape of the pdf. Then other statistical tests of convergence are conducted. We need to check the plots of Fqqs or SqLj versus the number of datasets used (from 1 to 100, implicating nj or Ns). If Fqqs or SqLj in plots is convergent with the increase of number of datasets used, the statistical convergence of Fqqs or SqLj is confirmed. Figures 5(b), 5(f), 6(b), and 6(f) also show that the shapes of pdf are close to 0 with the decrease of scale. So we investigate the convergence for q=10 with a larger scale and q=-10 with a smaller scale (see Figures 5(c), 5(d), 5(g), 5(h), 6(c), 6(d), 6(g), and 6(h) for Site 4 and Site 7). All the convergence plots except for Site 7 show that Fqqs and SqLj are convergent for the range of q from −10 to 10, the range of s from 25 to 27.2, and the range of j from 3 to 6. From Figures 6(d) and 6(h), it can be seen that Fqqs and SqLj at Site 7 are not convergent for q<0 with some small scales. This can be confirmed with the log-log plots in Figures 6(a) and 6(e).Figure 5
Statistical tests at Site 4. (a) and (e) show the log-log plots of MF-DFA and WLMF; different colors indicate different orders of polynomial (or different vanishing moments of wavelet). (a) 2-order (red), 3-order (blue), and 4-order (green). (e) 2-moment (red), 3-moment (blue), and 4-moment (green). (b) shows the pdf ofF2v,s1/2. (c) and (d) show the convergence of Fqqs. The results in (b), (c), and (d) are based on 3-order. (f) shows the pdf of LXj,k. (g) and (h) show the convergence of SqLj. The results in (f), (g), and (h) are based on 3-moment. These results correspond to an averaging over 100 sets of daily records.Figure 6
Statistical tests at Site 7. (a) and (e) show the log-log plots of MF-DFA and WLMF; different colors indicate different orders of polynomial (or different vanishing moments of wavelet). (a) 2-order (red), 3-order (blue), and 4-order (green). (e) 2-moment (red), 3-moment (blue), and 4-moment (green). (b) shows the pdf ofF2v,s1/2. (c) and (d) show the convergence of Fqqs. The results in (b), (c), and (d) are based on 3-order. (f) shows the pdf of LXj,k. (g) and (h) show the convergence of SqLj. The results in (f), (g), and (h) are based on 3-moment. These results correspond to an averaging over 100 sets of daily records.Figures5(a), 5(e), 6(a), and 6(e) show the log-log plots at Site 4 and Site 7 obtained by three different orders of polynomials in MF-DFA and three different vanishing moments of wavelet in WLMF. It can be seen that three different orders of polynomials in MF-DFA give the consistent results, as three different vanishing moments in WLMF do. The robustness of the results obtained by different orders of polynomials in MF-DFA and different vanishing moments of wavelet in WLMF is confirmed.Selection of Parameters. Before comparing MF-DFA and WLMF, we should select parameters to be used in the two methods. In MF-DFA, for detrending, following [16], we use a third-order polynomial to calculate the best polynomial fit of the profile in each segment, which can eliminate the second-order polynomial trend in the original series and is efficient for the hydrologic data in [16]. According to Section 3.3, the scales s should be lower than N/4, and the values of q should be both positive and negative. For all of the sites, the number of daily water level records is approximately 600, and the range of time scale s is selected from 25 to 27.2. The upper scale is selected for sufficient statistics. The lower scale is selected for avoiding abnormal points caused by small scales (see Figure 7 for details). According to the statistical tests above and following [14, 16], we choose the range -10≤q≤10. In WLMF, compared with the results of MF-DFA, we appropriately choose 3 as the number of the vanishing moment, which can eliminate the second-order polynomial trend in the original series. The scale j is chosen to be from 3 to 6.Figure 7
The log-scaling plots of the water level dynamics of all of the sites. The different symbols and colors indicate different moments:q=-10 (blue circles), q=-1 (red squares), and q=10 (green triangles). The colored lines indicate the regression lines of the corresponding colors from scale 25 to scale 27.2. These results correspond to an averaging over 100 sets of daily records. The errors bars are obtained from these sets and reveal 95% confidence interval. Some abnormal points at small scales are not shown in the figure because of their very small values such as points of scale 24.4 and scale 24.8 at Site 8.Multifractal Model. After using MF-DFA and WLMF, we apply the least-squares method to fit the generalized binomial multifractal model (see (16)). Some other multifractal models may be available such as LS ansatz model with the functional form(24)hq=1+H′-Cqα′-1-1α′-1,q>0,with the three parameters H′, C, and α′, which has been used successfully to describe the multifractal behaviour of rainfall and runoff records for q>0 [16, 44, 45]. However, the LS ansatz cannot be used for q<0. The two-parameter binomial model is suitable for both positive and negative values of q and has fewer parameters. So the two-parameter binomial model is relevant.Computational Results. By using the results of the fit of the generalized binomial multifractal model, we calculate Δα via (17). We can also calculate the fluctuation exponent H=h2 with MF-DFA; see Table 8 for details. H can be used to determine the long-range dependence because of its relationship with the correlation exponent γ. For a stationary time series, H is the well-known Hurst exponent.Table 8
The results of MF-DFA and WLMF for 10 sites.
MF-DFA
WLBMF
H
(
h
(
2
)
)
h
(
2
)’s errors bars
a
b
Δ
α
a
b
Δ
α
Site 1
0.88
0.069
0.40
0.70
0.81
0.48
0.97
1.00
Site 2
1.10
0.043
0.39
0.56
0.52
0.39
0.66
0.76
Site 3
0.90
0.040
0.34
0.66
0.95
0.26
0.83
1.66
Site 4
0.89
0.028
0.30
0.67
1.17
0.36
0.87
1.27
Site 5
0.92
0.014
0.26
0.70
1.44
0.26
0.83
1.66
Site 6
0.95
0.029
0.26
0.63
1.27
0.32
0.77
1.26
Site 7
0.88
0.039
0.37
0.66
0.82
0.52
0.86
0.73
Site 8
0.99
0.018
0.16
0.63
1.99
0.19
0.95
2.28
Site 9
1.20
0.027
0.29
0.54
0.88
0.35
0.64
0.87
Site 10
1.37
0.017
0.25
0.46
0.85
0.30
0.53
0.81
Note. Since the MF-DFA and WLMF fail to estimate h(q) for q<0 at Site 7, the fitting results at Site 7 are based on values of h(q) when q>0.Figure7 shows the log-scaling plots of the water level dynamics of all of the sites. It can be seen that except at Site 7, log2⁡Fqs increases approximately with log2⁡s for q=-10,-1,10 and can be fitted well by a straight line. hq can be estimated by the slope of the regression line. Moreover, MF-DFA fails to estimate hq for q<0 at Site 7.Figure8 shows the comparison of the shape of hq assessed, respectively, by MF-DFA and WLMF, and the fit of the two-parameter binomial model. From the comparison of the shape of hq, although the values of hq estimated by the two methods are different, the shapes assessed by MF-DFA and WLMF at all sites (except Site 7) are similar and suggest the two-parameter binomial model. Especially at Site 1, Site 4, Site 6, and Site 10, the shapes of hq assessed by MF-DFA and WLMF can be consistent by up- and downtranslation. This indicates that the multifractality of the water level records assessed by MF-DFA is credible. Since the MF-DFA and WLMF fail to estimate hq for q<0 at Site 7, the fitting results at Site 7 are based on values of hq when q>0. From Figure 8, it can be also seen that there are large errors bars for q<0 at Site 2, which also indicates poor estimation of hq for q<0.Figure 8
Comparison of the shape ofhq assessed by MF-DFA and WLMF. The curves are obtained by fits of the two-parameter binomial model, yielding (16). The different colors indicate the corresponding different fitting results. The resulting model parameters are reported in Table 8. These results correspond to an averaging over 100 sets of daily records. The errors bars are obtained from these sets and reveal 95% confidence interval.To further confirm the results of MF-DFA, we apply MF-DFA to the daily records at different times (morning, afternoon, and night) of a day. As seen in Figure9, the shapes of hq at different times are similar. Then, through the mutual verification of the results at different times, we can further confirm the results of MF-DFA that are shown in Figure 8.Figure 9
The generalized Hurst exponentshq for ten representative daily water level records at different times of a day. The errors bars are based on the linear least square fit which is used to estimate hq and reveal 95% confidence interval.From Figure8, we can also find some details of the results. As mentioned above, the estimation of hq for q<0 by MF-DFA and WLMF at Site 2 and Site 7 is poor, especially at Site 7. The adjacent sites are also affected by these two sites. It can be seen that at Site 6 and Site 8 hq near 0 deviate from their respective fitting curves. Excluding these sites, hq of most other sites can be fitted well with the two parameters (a, b). The records of most sites satisfy the two-parameter binomial model (or hq satisfy (16)). For the interpretation of the special behaviour at Site 2 and Site 7, we note that the frequent water storage and discharge activities at these two sites lead to a different water level property when compared to other sites. Moreover, apart from Site 2, Site 7, and those that are seriously affected by these two sites (i.e., Site 1, Site 6, and Site 8), the two parameters a and b exhibit a relationship such that a+b≈1, whereas in the standard binomial multifractal model [24], the relationship is close to a+b=1. However, for all of the sites downstream of Site 7, the sum of a and b is smaller than 0.9. The decrease in a+b can also be seen in the figures for WLMF. From Figure 2, we think that the four tributaries downstream of Site 7 have a significant impact on the river; that is, they change the multifractality of this river’s water level.From Table8, we can see that all of the Hurst exponents are larger than 0.8. This indicates that the water level records of the river have the long-range dependence property, which is the same as the result given by the wavelet-based method. The values of Δα are larger than 0.5, which indicates that there is strong multifractality at all sites.


## 4.1. Wavelet-Based Results

This section describes the wavelet-based estimation of Hurst parameters for the water level records of the 10 sites. We estimate the Hurst exponents at different scales and different sampling interval times. The results indicate that there is an interesting connection between the wavelet scale and the sampling period. For comparison, we repeat the procedure for a synthetic exponential long-range dependent process (ELRD) with length of 180000, which is close to the actual sampling data length. The vanishing moment of the used Daubechies wavelet isN=3 (db3). The simulation results for ELRD are listed in Tables 2 and 3. The values in brackets are errors bars. The errors bars are based on the weighted linear least square fit which is used to estimate H and reveal 95% confidence interval.Table 2
Simulation data analysis forH (H = 0.6, data length = 180000).
Interval
j
1~5
2~6
3~7
4~8
5~9
1
0.61 (0.009)
0.61 (0.013)
0.60 (0.018)
0.59 (0.025)
0.59 (0.036)
2
0.61 (0.013)
0.61 (0.018)
0.60 (0.025)
0.59 (0.036)
0.60 (0.051)
4
0.61 (0.018)
0.61 (0.025)
0.60 (0.036)
0.59 (0.051)
0.61 (0.071)
8
0.60 (0.025)
0.60 (0.036)
0.60 (0.051)
0.61 (0.071)
0.62 (0.106)Table 3
Simulation data analysis forH (H = 0.8, data length = 180000).
Interval
j
1~5
2~6
3~7
4~8
5~9
1
0.83 (0.009)
0.82 (0.013)
0.81 (0.018)
0.80 (0.025)
0.80 (0.034)
2
0.83 (0.013)
0.82 (0.018)
0.81 (0.025)
0.80 (0.034)
0.82 (0.049)
4
0.83 (0.018)
0.82 (0.025)
0.81 (0.034)
0.81 (0.049)
0.80 (0.072)
8
0.82 (0.025)
0.82 (0.034)
0.82 (0.049)
0.80 (0.072)
0.79 (0.100)The simulation results show that the simulations for bothH=0.6 and H=0.8 are stable under changes in the interval length and estimation scale. One can see that the wavelet-based estimation of ELRD for the Hurst exponent is valid and reliable, and the selection of wavelet scale and sampling cycle is independent of the estimation.Tables4–6 display the results of the practical water level. There exists an expression aij≈ai+1j-1 (aij is the element of the tables). For example, in Table 4,H = 0.82 (12 min, scales (1~5)) andH = 0.84 (6 min, scales (2~6)) are approximately equal;H = 0.70 (24 min, scales (1~5)),H = 0.68 (12 min, scales (2~6)) andH = 0.68 (6 min, scales (3~7)) are almost equal, and so forth. The results represent the consistent relationship between the wavelet scale and the sampling period. We present a mathematical proof in the Appendix.Table 4
Multiscales and multisampling forH of Site 1 (data length = 182855).
Interval
j
1~5
2~6
3~7
4~8
6 min
0.84 (0.062)
0.84 (0.047)
0.68 (0.046)
0.53 (0.054)
12 min
0.82 (0.067)
0.68 (0.054)
0.53 (0.054)
—
24 min
0.70 (0.046)
0.55 (0.056)
—
—
48 min
0.57 (0.053)
—
—
—Table 5
Multiscales and multisampling forH of Site 2 (data length = 152307).
Interval
j
1~5
2~6
3~7
4~8
6 min
0.80 (0.101)
0.86 (0.086)
0.81 (0.077)
0.61 (0.086)
12 min
0.89 (0.082)
0.80 (0.079)
0.61 (0.086)
0.50 (0.111)
24 min
0.81 (0.079)
0.63 (0.087)
0.50 (0.112)
—
48 min
0.65 (0.085)
0.52 (0.116)
—
—Table 6
Multiscales and multisampling forH of Site 3 (data length = 189967).
Interval
j
1~5
2~6
3~7
4~8
6 min
0.79 (0.047)
0.97 (0.041)
0.87 (0.044)
0.69 (0.048)
12 min
0.95 (0.046)
0.89 (0.041)
0.69 (0.047)
0.50 (0.059)
24 min
0.91 (0.044)
0.71 (0.048)
0.50 (0.059)
—
48 min
0.74 (0.046)
0.50 (0.060)
—
—We can also see that the estimations for the Hurst exponent are not stable under changes in the interval length and estimation scale, while the estimations along the secondary diagonal direction of the tables are still stable. We conclude that the massive hydrologic data have a multiscale consistent property, while the simulation data do not. This indicates that there may exist multifractality in the records.


## 4.2. Multiscale Diagram for Hydrological Data

In this section, we use the multiscale diagram to test the multifractality that may exist in the records. The multiscale diagram is the graph ofζ^q against q, together with the confidence intervals related to ζ^q. The multifractality can be detected by measuring the deviation of ζ^q from the linear form [26].We use a linear multiscale diagram to examine the alignment in the multiscale diagram, in whichHq=ζq/q is drawn against q. The multifractality can be determined by examining whether there is a horizontal alignment in the linear multiscale diagram. The horizontal alignment indicates that the Hurst parameter H can control all of the scaling property. The value of H can be received by estimating the level. The linear multiscale diagram and multiscale diagram are statistically equivalent.To determine the alignment of hydrological data, we impose the synthesized FBM (H=0.8) known as the “linear form” for comparison. The regressions for estimating ζq for each of the 10 different q are performed for small scales: j1,j2=2,6 according to the power-law behaviour of Sj,q discussed in Section 3.2. For Site 1–Site 10, as well as the FBM, the multiscale diagrams are shown in Figure 4, and the values of ζq for all 10 sites are also shown in Table 7 (the sampling interval is 6 minutes). We choose db3 as the wavelet of the wavelet transform, which is orthogonal to quadratic trends.Table 7
ζ
(
q
) for hydrological data of 10 sites; ζ(q) is estimated with a sampling interval of 6 min and (j1,j2)=(2,6).
q
0.5
1
1.5
2
2.5
3
4
5
6
8
Site 1
0.52
0.86
1.13
1.33
1.49
1.61
1.80
1.97
2.15
2.52
Errors bars
0.008
0.015
0.025
0.041
0.064
0.095
0.170
0.250
0.319
0.423
Site 2
0.51
0.97
1.34
1.56
1.68
1.76
1.87
1.97
2.09
2.39
Errors bars
0.010
0.023
0.046
0.083
0.123
0.160
0.233
0.308
0.383
0.526
Site 3
0.52
1.00
1.39
1.68
1.85
1.96
2.07
2.12
2.16
2.26
Errors bars
0.007
0.014
0.024
0.045
0.077
0.117
0.201
0.282
0.357
0.471
Site 4
0.59
1.16
1.61
1.70
1.59
1.47
1.27
1.15
1.07
1.00
Errors bars
0.007
0.017
0.064
0.156
0.222
0.272
0.362
0.441
0.511
0.618
Site 5
0.55
0.94
1.10
0.98
0.76
0.49
−0.06
−0.61
−1.13
−2.10
Errors bars
0.013
0.031
0.086
0.161
0.228
0.288
0.394
0.481
0.552
0.655
Site 6
0.42
0.80
1.12
1.37
1.58
1.73
1.91
2.02
2.10
2.26
Errors bars
0.008
0.015
0.023
0.037
0.060
0.090
0.161
0.227
0.285
0.386
Site 7
0.48
0.77
1.02
1.23
1.42
1.58
1.86
2.08
2.25
2.56
Errors bars
0.016
0.031
0.050
0.075
0.110
0.157
0.282
0.420
0.537
0.684
Site 8
0.58
1.07
1.52
1.77
1.82
1.80
1.74
1.71
1.71
1.76
Errors bars
0.015
0.027
0.067
0.188
0.334
0.437
0.542
0.589
0.615
0.653
Site 9
0.46
0.88
1.25
1.51
1.66
1.76
1.98
2.22
2.46
2.94
Errors bars
0.007
0.015
0.035
0.079
0.140
0.202
0.304
0.373
0.423
0.495
Site 10
0.59
1.12
1.59
1.86
1.95
1.98
2.03
2.11
2.22
2.50
Errors bars
0.012
0.024
0.057
0.135
0.220
0.286
0.383
0.455
0.513
0.599
Note. The errors bars are based on the weighted linear least square fit which is used to estimate ζ(q) and reveal 95% confidence interval.Figure 4
Testing for multiscaling. (a) Multiscale diagram. Compared with the alignment for FBM, hydrological data with no alignment suggests multiscale behaviour. (b) Linear multiscale diagram. The FBM series exhibits horizontal alignment. The errors bars are based on the weighted linear least square fit which is used to estimateζq and reveal 95% confidence interval.
(a)
(b)In Figure4, the multiscale diagram in (a) shows a lack of alignment for the hydrologic data compared to the clear alignment of FBM for its self-similarity, which indicates that it has multifractal behaviour. Thus, the multifractal model may be relevant to the hydrological data. The linear multiscale diagram in (b) confirms this conclusion. The Hurst parameter H of the synthesized FBM series can be estimated by the horizontal alignment, and its estimator is approximately 0.8. This is not the case for the hydrologic data. In summary, the hydrological data have multifractal behaviour and do not comply with the equation ζq=Hq. The errors bars in Figure 4 are based on the weighted linear least square fit, which is used to estimate ζq. They are calculated as 1.96×Stdζq. Under the assumption of normality, this reveals a 95% confidence interval. The detailed values of ζq for other sites can be seen in Table 7.


## 4.3. Multifractal Analysis and Results

In this section, we analyze the daily water level records at the 10 observation sites using MF-DFA analysis for the multifractal behaviour. In order to ensure that there are enough fluctuations in each segment so that MF-DFA can be applied, the high-frequency collected daily records are used. For sufficient samples, we get 100 sets of daily records randomly at different times of one day. The number of records of each sets is about six hundred. We also apply wavelet leaders for multifractal analysis (WLMF) based on the daily water level records for a comparison with MF-DFA.Statistical Tests. It is well known that finite-size effects and lack of statistical convergence will introduce biases on the estimate of hq when applying the method of MF-DFA and WLMF. Also the robustness of the results obtained by different orders of polynomials in MF-DFA and different vanishing moments of wavelet in WLMF should be tested. In this part, we perform some statistical tests for these issues.Following the test method reported by Khalil et al. [42], we consider the problem of statistical convergence of hq through the statistical convergence of Fqqs and SqLj in the two methods. Recall that Fqqs and SqLj can be transformed into a continuous integral over M,(23)SqLj=1nj∑k=1njLXj,kq=∫MqPjMdM,Fqqs=1Ns∑v=1NsF2v,sq/2=∫MqPsMdM.M means the modulus of LXj,k or F2v,s1/2. PjM means the probability density function (pdf) of LXj,k. PsM means the probability density function (pdf) of F2v,s1/2. The well-defined shape of the pdf’s plots indicates a proper convergence [42, 43]. Figures 5(b), 5(f), 6(b), and 6(f) show PjM and PsM computed at three different scales at Site 4 and Site 7. It can be seen that there are not enough statistics to give a well-defined shape of the pdf. However, the statistical convergence of Fqqs and SqLj does not need so many statistics that can give a well-defined shape of the pdf. Then other statistical tests of convergence are conducted. We need to check the plots of Fqqs or SqLj versus the number of datasets used (from 1 to 100, implicating nj or Ns). If Fqqs or SqLj in plots is convergent with the increase of number of datasets used, the statistical convergence of Fqqs or SqLj is confirmed. Figures 5(b), 5(f), 6(b), and 6(f) also show that the shapes of pdf are close to 0 with the decrease of scale. So we investigate the convergence for q=10 with a larger scale and q=-10 with a smaller scale (see Figures 5(c), 5(d), 5(g), 5(h), 6(c), 6(d), 6(g), and 6(h) for Site 4 and Site 7). All the convergence plots except for Site 7 show that Fqqs and SqLj are convergent for the range of q from −10 to 10, the range of s from 25 to 27.2, and the range of j from 3 to 6. From Figures 6(d) and 6(h), it can be seen that Fqqs and SqLj at Site 7 are not convergent for q<0 with some small scales. This can be confirmed with the log-log plots in Figures 6(a) and 6(e).Figure 5
Statistical tests at Site 4. (a) and (e) show the log-log plots of MF-DFA and WLMF; different colors indicate different orders of polynomial (or different vanishing moments of wavelet). (a) 2-order (red), 3-order (blue), and 4-order (green). (e) 2-moment (red), 3-moment (blue), and 4-moment (green). (b) shows the pdf ofF2v,s1/2. (c) and (d) show the convergence of Fqqs. The results in (b), (c), and (d) are based on 3-order. (f) shows the pdf of LXj,k. (g) and (h) show the convergence of SqLj. The results in (f), (g), and (h) are based on 3-moment. These results correspond to an averaging over 100 sets of daily records.Figure 6
Statistical tests at Site 7. (a) and (e) show the log-log plots of MF-DFA and WLMF; different colors indicate different orders of polynomial (or different vanishing moments of wavelet). (a) 2-order (red), 3-order (blue), and 4-order (green). (e) 2-moment (red), 3-moment (blue), and 4-moment (green). (b) shows the pdf ofF2v,s1/2. (c) and (d) show the convergence of Fqqs. The results in (b), (c), and (d) are based on 3-order. (f) shows the pdf of LXj,k. (g) and (h) show the convergence of SqLj. The results in (f), (g), and (h) are based on 3-moment. These results correspond to an averaging over 100 sets of daily records.Figures5(a), 5(e), 6(a), and 6(e) show the log-log plots at Site 4 and Site 7 obtained by three different orders of polynomials in MF-DFA and three different vanishing moments of wavelet in WLMF. It can be seen that three different orders of polynomials in MF-DFA give the consistent results, as three different vanishing moments in WLMF do. The robustness of the results obtained by different orders of polynomials in MF-DFA and different vanishing moments of wavelet in WLMF is confirmed.Selection of Parameters. Before comparing MF-DFA and WLMF, we should select parameters to be used in the two methods. In MF-DFA, for detrending, following [16], we use a third-order polynomial to calculate the best polynomial fit of the profile in each segment, which can eliminate the second-order polynomial trend in the original series and is efficient for the hydrologic data in [16]. According to Section 3.3, the scales s should be lower than N/4, and the values of q should be both positive and negative. For all of the sites, the number of daily water level records is approximately 600, and the range of time scale s is selected from 25 to 27.2. The upper scale is selected for sufficient statistics. The lower scale is selected for avoiding abnormal points caused by small scales (see Figure 7 for details). According to the statistical tests above and following [14, 16], we choose the range -10≤q≤10. In WLMF, compared with the results of MF-DFA, we appropriately choose 3 as the number of the vanishing moment, which can eliminate the second-order polynomial trend in the original series. The scale j is chosen to be from 3 to 6.Figure 7
The log-scaling plots of the water level dynamics of all of the sites. The different symbols and colors indicate different moments:q=-10 (blue circles), q=-1 (red squares), and q=10 (green triangles). The colored lines indicate the regression lines of the corresponding colors from scale 25 to scale 27.2. These results correspond to an averaging over 100 sets of daily records. The errors bars are obtained from these sets and reveal 95% confidence interval. Some abnormal points at small scales are not shown in the figure because of their very small values such as points of scale 24.4 and scale 24.8 at Site 8.Multifractal Model. After using MF-DFA and WLMF, we apply the least-squares method to fit the generalized binomial multifractal model (see (16)). Some other multifractal models may be available such as LS ansatz model with the functional form(24)hq=1+H′-Cqα′-1-1α′-1,q>0,with the three parameters H′, C, and α′, which has been used successfully to describe the multifractal behaviour of rainfall and runoff records for q>0 [16, 44, 45]. However, the LS ansatz cannot be used for q<0. The two-parameter binomial model is suitable for both positive and negative values of q and has fewer parameters. So the two-parameter binomial model is relevant.Computational Results. By using the results of the fit of the generalized binomial multifractal model, we calculate Δα via (17). We can also calculate the fluctuation exponent H=h2 with MF-DFA; see Table 8 for details. H can be used to determine the long-range dependence because of its relationship with the correlation exponent γ. For a stationary time series, H is the well-known Hurst exponent.Table 8
The results of MF-DFA and WLMF for 10 sites.
MF-DFA
WLBMF
H
(
h
(
2
)
)
h
(
2
)’s errors bars
a
b
Δ
α
a
b
Δ
α
Site 1
0.88
0.069
0.40
0.70
0.81
0.48
0.97
1.00
Site 2
1.10
0.043
0.39
0.56
0.52
0.39
0.66
0.76
Site 3
0.90
0.040
0.34
0.66
0.95
0.26
0.83
1.66
Site 4
0.89
0.028
0.30
0.67
1.17
0.36
0.87
1.27
Site 5
0.92
0.014
0.26
0.70
1.44
0.26
0.83
1.66
Site 6
0.95
0.029
0.26
0.63
1.27
0.32
0.77
1.26
Site 7
0.88
0.039
0.37
0.66
0.82
0.52
0.86
0.73
Site 8
0.99
0.018
0.16
0.63
1.99
0.19
0.95
2.28
Site 9
1.20
0.027
0.29
0.54
0.88
0.35
0.64
0.87
Site 10
1.37
0.017
0.25
0.46
0.85
0.30
0.53
0.81
Note. Since the MF-DFA and WLMF fail to estimate h(q) for q<0 at Site 7, the fitting results at Site 7 are based on values of h(q) when q>0.Figure7 shows the log-scaling plots of the water level dynamics of all of the sites. It can be seen that except at Site 7, log2⁡Fqs increases approximately with log2⁡s for q=-10,-1,10 and can be fitted well by a straight line. hq can be estimated by the slope of the regression line. Moreover, MF-DFA fails to estimate hq for q<0 at Site 7.Figure8 shows the comparison of the shape of hq assessed, respectively, by MF-DFA and WLMF, and the fit of the two-parameter binomial model. From the comparison of the shape of hq, although the values of hq estimated by the two methods are different, the shapes assessed by MF-DFA and WLMF at all sites (except Site 7) are similar and suggest the two-parameter binomial model. Especially at Site 1, Site 4, Site 6, and Site 10, the shapes of hq assessed by MF-DFA and WLMF can be consistent by up- and downtranslation. This indicates that the multifractality of the water level records assessed by MF-DFA is credible. Since the MF-DFA and WLMF fail to estimate hq for q<0 at Site 7, the fitting results at Site 7 are based on values of hq when q>0. From Figure 8, it can be also seen that there are large errors bars for q<0 at Site 2, which also indicates poor estimation of hq for q<0.Figure 8
Comparison of the shape ofhq assessed by MF-DFA and WLMF. The curves are obtained by fits of the two-parameter binomial model, yielding (16). The different colors indicate the corresponding different fitting results. The resulting model parameters are reported in Table 8. These results correspond to an averaging over 100 sets of daily records. The errors bars are obtained from these sets and reveal 95% confidence interval.To further confirm the results of MF-DFA, we apply MF-DFA to the daily records at different times (morning, afternoon, and night) of a day. As seen in Figure9, the shapes of hq at different times are similar. Then, through the mutual verification of the results at different times, we can further confirm the results of MF-DFA that are shown in Figure 8.Figure 9
The generalized Hurst exponentshq for ten representative daily water level records at different times of a day. The errors bars are based on the linear least square fit which is used to estimate hq and reveal 95% confidence interval.From Figure8, we can also find some details of the results. As mentioned above, the estimation of hq for q<0 by MF-DFA and WLMF at Site 2 and Site 7 is poor, especially at Site 7. The adjacent sites are also affected by these two sites. It can be seen that at Site 6 and Site 8 hq near 0 deviate from their respective fitting curves. Excluding these sites, hq of most other sites can be fitted well with the two parameters (a, b). The records of most sites satisfy the two-parameter binomial model (or hq satisfy (16)). For the interpretation of the special behaviour at Site 2 and Site 7, we note that the frequent water storage and discharge activities at these two sites lead to a different water level property when compared to other sites. Moreover, apart from Site 2, Site 7, and those that are seriously affected by these two sites (i.e., Site 1, Site 6, and Site 8), the two parameters a and b exhibit a relationship such that a+b≈1, whereas in the standard binomial multifractal model [24], the relationship is close to a+b=1. However, for all of the sites downstream of Site 7, the sum of a and b is smaller than 0.9. The decrease in a+b can also be seen in the figures for WLMF. From Figure 2, we think that the four tributaries downstream of Site 7 have a significant impact on the river; that is, they change the multifractality of this river’s water level.From Table8, we can see that all of the Hurst exponents are larger than 0.8. This indicates that the water level records of the river have the long-range dependence property, which is the same as the result given by the wavelet-based method. The values of Δα are larger than 0.5, which indicates that there is strong multifractality at all sites.


## 5. Discussion

The long-range dependence property and multifractality we detected in the water level records are identical to the previous results of hydrologic data (see, e.g., [14–16]), and our results are based on multiple observation sites and a high sampling frequency that we have not yet seen in other studies.In the previous results (see, e.g., [14–16]), we would like to note that the lateral inflows have an effect on the runoffs (or water level). Our results of the change of a+b seem to suggest that the contribution of lateral inflow (tributaries) leads to the change in the multifractality of the trunk stream’s water level.In light of the results obtained by MF-DFA, the values ofa are between 0.12 and 0.40 for all of the sites, and the values of b are between 0.54 and 0.70. From the previous results of 41 rivers’ runoffs studied by Koscielny-Bunde et al. [14], the values of a and b of the Würm river are in the range of our values. We would like to note that the Würm river is located in southern Germany (near Munich), which has a temperate continental climate. The river is also within a localized area, and its climate and area are similar to those of the analyzed river. The consistency of the two rivers in parameters validate the reliability of the obtained multifractal model.As previously mentioned, the spurious multifractality may be induced by intricate nonstationarities, which often exist in the hydrologic data, for example, due to the seasonal cycle or a change in climate. We consider that the obtained multifractality is reliable. First the method we used are based on wavelet analysis and MF-DFA which can eliminate some trends and avoid spurious detection of correlations. Second the multifractality is confirmed consistently via multiple methods (multiscale diagram, MF-DFA, and WLMF). The MF-DFA and WLMF also give the same multifractal model, and this multifractal model is robust in different location and different time. Last the results obtained from the analyzed river are consistent with the previous study of actual example.


## 6. Conclusions

In this paper, we proceed in an orderly manner to study the LRD and multifractal properties of high-frequency water level records of a river in northern China by using wavelet analysis and MF-DFA.We observe the long-range dependence property and multiscale behaviour in the time series, which indicates the multifractality. This multifractality is detected via multiscale diagram. Then we focus on the multifractality of records using MF-DFA and WLMF. At most sites, the generalized Hurst exponenthq can be described with the function hq=1/q-ln⁡aq+bq/(qln⁡2),a,b∈0,1 and the multifractal models are robust for different times (morning, afternoon, and evening). The validity of multifractality is (partially) supported with the same model obtained from the different methods (MF-DFA and WLMF) and by comparing previous study of actual example. The results also show that the poor estimation of hq for q<0 at Site 2 and Site 7 and the decreased change in a+b indicates human activities and tributaries have a great influence on the multifractality of the river.We obtain and list the values ofH, a, b, and Δα for all available records, including some abnormal records. At all available sites, MF-DFA yields values of Δα are larger than 0.5, whereas WLMF yields values are larger than 0.7, which shows a strong multifractality in the data.As a technical improvement in the estimation of the Hurst parameter, we show that the sampling intervals and choice of scales of the discrete wavelet transformation should be consistent (Tables4–6).In contrast to the previous studies, our studies are based on a short sampling time, which we have not seen before. This study contributes a new insight into previous research applications, and it may also be the extension of the study of multifractality. Further investigations are needed, especially for the massive data cases, for example, comparing the fitting ofhq that is obtained from the records partitioned for each year. As mentioned above, the human activities (frequent water storage and discharge activities) at Site 2 and Site 7 lead to poor estimation of hq for q<0 and the change in a+b. This is the key point for the future work on detection and prediction. A real time multifractal model may be built for hq,t, at, and bt. The detection and prediction of human activities and climate change in real time may be conducted by the change of hq,t, at, and bt. Besides, it would be also interesting to consider the water level dynamics as a multifractional process, study Ht at the different sites with different time lags, and provide predictive models of the water level records.



---

*Source: 3148257-2017-10-31.xml*