Joint Multifractal Analysis and Source Testing of River Level Records Based on Multifractal Detrended Cross-Correlation Analysis

The joint multifractal analysis is usually conducted in two diﬀerent variables for their cross-correlations but rarely used for two records of one variable collected at two diﬀerent places. It is important for the detection of change in multifractality in space. Besides, the cross-correlations in two analyzed series make the analysis of sources of joint multifractality diﬃcult. There are few studies on the source of joint multifractality. We focus on the two issues for two level records at pairs of adjacent sites along one river and carry out an extension of our previous work which is about the single multifractality of one record with the same data set. The data set is collected from 10 observation stations of a northern China river and contains about two million high-frequency river level records. Results of joint multifractal analysis via multifractal detrended cross-correlation analysis show that the change in joint multifractality at pairs of adjacent sites caused by weak cross-correlations can be detected by comparing the single generalized Hurst exponent with the joint scaling exponent function and reveal the eﬀects of human activities on joint multifractality. This analysis provides an approach for detecting the change in multifractality. Following the idea of our previous work, two robust hypothesis tests via a set of pairs of surrogate series are proposed for the source testing of joint multifractality. The analysis of the eﬀects of cross-correlations is carried out via a proposed simultaneously half-shifting technique which can both minimize the cross-correlations between original series and make full use of records. Results of source analysis show not only the eﬀects of autocorrelations in series and probability distribution of river levels but also the eﬀects of cross-correlations between series.


Introduction
e multifractality has been studied in many areas of applied sciences since its concept was first introduced by Mandelbrot [1]. e structure function method [2] is a classical multifractal analysis method which was generalized from boxcounting algorithms [3][4][5] and once widely used for multifractal analysis of time series. Until now, many multifractal analysis methods were developed, such as multifractal detrended fluctuation analysis (MF-DFA) [6] and wavelet transform modulus maxima (WTMM) [7] and its discrete version, namely, wavelet leaders multifractal analysis (WLMF) [8][9][10][11]. ese methods have been widely applied to many areas of applied sciences instead of the structure function method [12][13][14][15][16][17][18][19][20][21][22]. e MF-DFA is a multifractal generalization of the detrended fluctuation analysis (DFA) that is used to estimate Hurst exponent of monofractal process [17,23,24]. Similarly, the monofractal versions of wavelet-based methods can be also used for Hurst exponent [25][26][27][28][29][30][31]. Besides, many complex systems usually have the joint multifractal nature which is exhibited by simultaneous records of two variables. ere are also many methods for this issue such as multifractal detrended cross-correlation analysis (MF-DCCA) [32], multifractal cross wavelet analysis (MF-X-WT) [33], and joint multifractal analysis based on wavelet leaders (MF-X-WL) [34]. e MF-DCCA which was first proposed by Zhou [32] is a generalization of the MF-DFA [6] for two nonstationary time series. As MF-DFA, the MF-DCCA is easy to implement and robust when time series is short. Besides, it can also eliminate some polynomial trends that may exist in nonstationary time series. e MF-DCCA has been a commonly used joint multifractal analysis method in many fields including hydrology [35][36][37][38]. erefore, we conduct the MF-DCCA for the joint multifractal analysis of river level records. e river levels we analyze are collected from the whole 10 observation stations of a tributary of Haihe River in North China and contain about two million records. e sampling frequency is 6 minutes. e multifractal analysis of these data and their source testing has been conducted in [39,40]. Such many and high-frequency records are rarely seen in previous literatures. Venugopal et al. [41][42][43] studied the multifractality of high-resolution temporal rainfall data using wavelet methods.
is paper follows up our previous work [39,40] with the same data set and focuses on the joint multifractality and its source testing of river levels at adjacent stations. e joint multifractality is exhibited by two simultaneous records and distinguished from the single multifractality of one record studied in [40]. e joint multifractal analysis is usually conducted in two different variables for their cross-correlations but rarely used for two records of one variable collected at two different places like the data we analyze in this paper. It is very important for the detection of change in multifractality in space. In our study, the joint multifractal analysis of two river level records at adjacent stations is conducted using MF-DCCA. e obtained joint multifractal results are used to analyze the change in multifractality along the river and can detect the change in multifractality affected by human activities. It provides an approach for detecting the change in multifractality. Besides, there are also few studies on the source of joint multifractality. Jiang et al. [33] studied the source of joint multifractality in financial series using MF-X-WT. It is known that the single multifractality of one record has two main sources: the fatness of the probability distribution of original time series and the different correlations in small and large fluctuations. e two sources can be usually distinguished by analyzing the corresponding surrogate series [44][45][46][47][48]. In addition to the two sources that cause single multifractality, the cross-correlations that exist in two analyzed records also affect the joint multifractality. e cross-correlations in two analyzed records make some surrogate series methods ineffective and the analysis of sources of joint multifractality difficult. Following the idea in [40], we propose two hypothesis tests for the source of joint multifractality, which are based on the empirical distributions of scaling exponent function estimated from some pairs of surrogate series. e pairs of surrogate series are generated via shuffling method and rank-ordered remapping technique which are still effective for joint multifractality. For further testing the effect of cross-correlations, we also propose a simultaneously half-shifting technique that can minimize cross-correlations. Testing results show that the joint multifractality in water levels is mainly caused by the different correlations within and between original series and is also related to the probability distribution of river levels. e simultaneously half-shifting technique further shows the effects of cross-correlations. e MF-DCCA and related comparison results not only show the effects of cross-correlations but also can detect the change of multifractality. e paper is organized as follows: in Section 2, we introduce the MF-DCCA method and propose the source testing method and the simultaneously half-shifting technique. e data and related preliminaries are described in Section 3. In Section 4, we report the results and give detailed analysis. Our work is concluded in Section 5.

Multifractal Detrended Cross-Correlation Analysis.
In this section, we introduce the MF-DCCA [32,49] which is used for joint multifractal analysis of river level fluctuations in this work. e MF-DCCA is a generalization of the MF-DFA [6] for two nonstationary time series. e method is described as follows.
Let X(t) and Y(t) be the two processes to be analyzed. First, the profile is computed as follows: where L is the minimum length of X(t) and Y(t) and X and Y are the means of X(t) and Y(t), respectively.
Dividing the profiles X ′ (t) and Y ′ (t) into N s : � int(L/s) nonoverlapping segments with equal length s and then computing the detrended covariance for each segment, we get where x v (i) and y v (i) are, respectively, the best polynomial fit of X ′ [(v − 1)s + i] and Y ′ [(v − 1)s + i] in each of the N s segments v. In this step, the trends are removed via the polynomial fit of profile. e m-order polynomial used in the fitting procedure can eliminate the polynomial trend of order m − 1 in the original series. is ensures the fluctuation analysis available for data affected by trends or other nonstationarities. Note that the detrended covariance computed by equation (2) uses absolute values of fluctuations for fluctuation strength, which is different from that of 2 Complexity original paper [32] using absolute values of F XY (v, s) for fluctuation itself. Besides, the operation of absolute value here can also avoid no obvious power-law scaling when q < 0 [32,37], which exists in this study. Average the fluctuation function F XY (v, s) over all the segments, given by where q can be any real values except 0. For q � 0, e final step is calculating h XY (q) via the slope of the log-log plot of F XY (q, s) versus s, which is based on the power-law [32]: e function h XY (q) is called the scaling exponent function [32]. When q � 2, h XY (q) is the bivariate Hurst exponent [50].
When X � Y, MF-DCCA degenerates to MF-DFA and h XY (q) is called the generalized Hurst exponent [6]. Let h X (q) and h Y (q) be the generalized Hurst exponents of X(t) and Y(t) estimated using MF-DFA, respectively. For binomial multifractal measures and some multifractal random walks, it is validated numerically in the case of dependent pairs [32,51]: Furthermore, more general relationship is derived [37]: where k is a positive constant. When scale s ⟶ + ∞, Since the computation of F XY (q, s) in equation (3) is determined by F XY (v, s) from the segments with large fluctuation for q > 0, h XY (q) of q > 0 shows the scaling behavior of the segments with large fluctuations. For q < 0, h XY (q) shows that with small fluctuations [6,17].
e joint multifractal spectrum f XY (α) is obtained via the Legendre transform of Renyi scaling exponent τ XY (q): where α is the singularity strength and τ XY (q) can be obtained by calculating τ XY (q) � qh XY (q) − 1 [35,38]. e strength of joint multifractality of records can be characterized by Δα XY � α max − α min . e dynamics of multifractal system are described by a continuous spectrum of exponents (multifractal spectrum), rather than a single exponent (fractal dimension). is indicates that the exponent α (singularity strength or Hölder exponent) values continuously in one range. e width of this range is the width of multifractal spectrum Δα XY . It can reflect the complexity of exponent value and characterize the strength of multifractality. Inherently, the multifractal formalism is originally used to characterize the complexity of invariant measures of nonlinear dynamical systems.
If h(q) satisfies the generalized binomial multifractal model [18] then the τ(q) and Δα can be calculated via the following formulas:

Testing Method for Source of Joint Multifractality.
is section focuses on the testing for the source of joint multifractality. As shown in previous articles [33,44,47,48,51], the source of multifractality can be clarified by analyzing the corresponding surrogate series. Unlike the source testing of single series, we need to consider cross-correlations between series which also play an important role in generating joint mulitfractality [33,34]. In view of the above, we propose two hypothesis tests using pairs of surrogate series for the following two different sources: (1) e different correlations within and between original series. (2) e fatness of probability distribution of original series.
e following two surrogate methods [51] can generate pairs of surrogate series which keep one of the two sources and destroy another. Let X(t) and Y(t) be the two series to be analyzed (t � 1, . . . , L).
(2) For testing Type (2) source, the pair of surrogate series (X sg2 (t), Y sg2 (t))|tt � n1, 2q, h . . . , xL is generated using rank-ordered remapping technique: let Z(t)|tt � n1, 2q, h . . . , xL be a sequence of random numbers generated from the standard normal distribution. We rearrange Z(t) { } such that the rearranged series Z X (t)|tt � n1, 2q, h . . . , xL has the same rank ordering as the original series }. e series Z Y (t)|tt � n1, 2q, h . . . , xL is generated using the same technique for the original series , where σ X and σ Y are the sample standard deviations of X(t) and Y(t), respectively. e pair of surrogate series (X sg1 (t), Y sg1 (t)) preserves the distribution of original series and removes the correlations within and between original series. If the joint multifractality is caused in part by different correlations within and between original series, their joint multifractality will be significantly different from that of two original series. e pair of surrogate series (X sg2 (t), Y sg2 (t)) preserves the correlations within and between original series, and both have the Gaussian distribution. If the joint multifractality is caused in part by the fat-tailed distribution of original series, they will show a weaker joint multifractality than original series. en, the two types of joint multifractality can be distinguished by comparing h XY (q) of two original series with that of corresponding pair of surrogate series.
Following the idea in [40], we propose two hypothesis tests for the source of joint multifractality, which are based on the empirical distribution of scaling exponent function estimated from 1000 pairs of surrogate series. It can not only achieve more robust results than the previous clarification method using one pair of surrogate series but also show more details about the source of joint multifractality such as the effects of large and small fluctuations. e null hypotheses for the two hypothesis tests are described as follows: (i) H 0sg1 : the joint multifractality is not due to the correlations within and between series. (ii) H 0sg2 : the joint multifractality is not due to the fatness of probability distribution.
Monte Carlo simulation is adopted for the acceptance region of null hypotheses [52]. e proposed procedure for the two hypothesis tests is described as follows: (i) Assume that the significance level of test is α.
Generate 1000 pairs of surrogate series and calculate the scaling exponent function h sg1 (q) (or h sg2 (q)) using MF-DCCA to obtain the empirical probability distribution of the scaling exponent function, where h sg1 (q) and h sg2 (q) denote the scaling exponent for null hypothesis at the significance level α i � (α/N), based on the empirical probability distribution. In detail, h l (q i ) is estimated by the (α i /2)-quantile of the empirical probability distribution and where N is the number of the analyzed qs.
(iii) Connect h l (q i )s and h u (q i )s, respectively, using straight line. Since h l (q i )s or h u (q i )s varies continuously, it is suitable to connect using straight line. e area surrounded by these connection lines from q 1 to q N is the acceptance region.
(iv) Calculate the scaling exponent function h XY (q) of original series using MF-DCCA. (v) Accomplish the hypothesis testing via comparing the h XY (q) with the acceptance region.
Take testing for the null hypothesis H 0sg1 as an example. If h(q) falls into the acceptance region for H 0sg1 , we accept that the joint multifractality is not mainly due to different correlations within and between series. If h(q) falls out of the acceptance region for H 0sg1 , we believe that the joint multifractality is caused in part by different correlations within and between series.
Note that the significance level for each q is α i � (α/N) because this procedure is a multiple-testing procedure so that Bonferroni correction is conducted. In this paper, N � 31 and we use 5% as its significance level. So the significance level for each q is α i ≈ 0.16%.

Simultaneously Half-Shifting
Technique. As shown in [33], the shift of two series does not change the correlations within each series but can weaken cross-correlations between them. For further analyzing the effects of crosscorrelations on joint multifractality, we propose the following technique in MF-DCCA for both minimizing the cross-correlations between original series and making full use of records: . , X(L) . Two constructed pairs of shifting series are (X 0 (t), Y 1 (t)) and (X 1 (t), Y 0 (t)) . ey keep the autocorrelations in each series and minimizing the cross-correlations between series. (iii) For making full use of records, the joint multifractality without cross-correlations is obtained by simultaneously conducting two pairs of shifting series using MF-DCCA. For details, in the step of computation of F XY (q, s), for each s, the F XY (q, s)s is computed via averaging the fluctuations of both . e other steps are the same.
is technique simultaneously uses two pairs of series that are shifted half their length relative to each other. We thus call it simultaneously half-shifting technique (SHST). e scaling exponent function h XYS (q) obtained via simultaneously half-shifting technique is compared with the scaling exponent function estimated from two original series and the acceptance region for H 0sg1 .
For MF-DFA of single series X(t), there are no crosscorrelations. In this case, the simultaneously half-shifting technique degenerates the operation of half-dividing; that is, the X(t) is divided into two series X 0 (t) and X 1 (t) which are defined above. Further for multifractal detrended fluctuation analysis of the two divided series, F X (q, s)s is computed via averaging the fluctuations of both X 0 (t) and X 1 (t). e 4 Complexity operation of half-dividing should not change the multifractality since it can keep the autocorrelations in series and the distribution of series. We can analyze the change in multifractality of each series (X(t) and Y(t)) after the operation of half-dividing for validating the simultaneously half-shifting technique in keeping autocorrelations.

Data and Preliminaries
e river level records we consider have been analyzed for single multifractality in previous papers [39,40]. e records are collected in real time from 10 observations along a tributary of Haihe River in North China. e locations of 10 observations are shown in Figure 1. We mark these locations from upstream to downstream with Site 1-Site 10. e records are of high-frequency, collected every 6 minutes at 10 water level observations of this river from April 2011 to September 2013. e location of river has a significant continental monsoon climate, which is cold and dry in the winter and warm and humid in the summer. More information about records can be seen in Section 3 of [40] and Table 1. e sudden change in records in Table 1 at Site 2 and Site 7 is led by water storage and discharge at the two sites. Besides, there are some tributaries located downstream of Site 7 that flows into this river. e main tributary of this river is located at Site 9.
e preprocessing of records considers three issues: (1) the noise effects led by high-frequency data; (2) the missing data led by abnormal work of data-collecting equipment; and (3) seasonal trend in the records led by monsoon climate. Following the literature [40], we conduct the preaverage of one day's record and the linear interpolation, respectively, for the first two issues. e preaverage of one day's record allows us to use the daily average river levels for analysis. More detail about the preprocessing can be seen in Section 3 of [40].
For the third issue, it is known that MF-DCCA cannot remove the seasonal periodic trends in the records, which cause cross-overs in fluctuation function and further give incorrect results [53][54][55][56][57][58]. ere are some robust methods for the time series with periodic trends such as Fourier detrended fluctuation analysis (F-DFA) [54,55], singular value decomposition (SVD) [56,57,59], and the polynomial of varying order l [60]. e method of SVD can achieve the same goal as F-DFA [36]. A useful review on this issue can refer to the literature [58]. In this work, we use F-DFA for removing the seasonal periodic trends. After truncating the first 38 coefficients of low Fourier frequencies, we obtain the preprocessed daily average river levels for analysis, which show no cross-overs in fluctuation function (see Figure 2).
Note that since we calculate joint multifractality of two records, the two records used for analysis should keep the same in time.
is results in less available data than the sampling days.

Multifractal Results.
In this section, we study the joint multifractality between the preprocessed daily average river levels using MF-DCCA and mainly focus on the change in joint multifractality along the river. So the joint multifractal properties of two records at Site k and Site k + 1 are obtained, where k � 1, . . . , 9. Parameters for MF-DCCA are selected following [39,40]. ey have been proved to be robust and useful for the same records. e third-order polynomial is adopted to calculate the best polynomial fit of the profile in each segment, which can eliminate the secondorder polynomial trend in the original series. Limited by the series length and following [39,40], the scaling range of s is selected from 2 5.1 to 2 7.2 every 0.1 power. e length of river levels is short and only about 800. For avoiding inaccurate results at large q caused by the finite-size effects [61][62][63][64], the q-range cannot be large. e range of q is chosen − 6 ≤ q ≤ 6 carefully. e range is large enough to contain accurate h XY (q)s for our analysis. Figure 2 shows the log-scaling plots of F XY (q, s) versus s for two preprocessed daily average river levels at Site k and Site k + 1. It can be seen that although there is some volatility in log 2 F XY (q, s) which is led by short series, log 2 F XY (q, s) increases approximately linearly with log 2 s for all qs at all adjacent sites. erefore, there exists the cross power-law behavior of two preprocessed daily average river levels at adjacent sites. e joint multifractal results between two preprocessed daily average river levels at Site k and Site k + 1 are shown in Figure 3. e errors bars of h XY (q) are based on the linear least square fit. It is shown in Section 2.1 that there exist some relationships between h XY (q) and (h X (q) + h Y (q))/2. We also plot the values of (h X (q) + h Y (q))/2 in Figure 3 for comparison. From Figure 3, we can find the values of (h X (q) + h Y (q))/2 deviate obviously from the values of h XY (q) and fall out of the errors bars for q < 0 at all adjacent sites except Site 3-Site 4 and Site 9-Site 10. e deviations of all pairs of adjacent sites become smaller as they get away from Site 2 and Site 7. As mentioned in Section 2.1, it has been validated numerically in the case of dependent pairs [32,51]: For dependent pairs, the values of (h X (q) + h Y (q))/2 should not differ much from that of h XY (q). So the deviations of h XY (q) from (h X (q) + h Y (q))/2 can reflect the strength of cross-correlations. e small deviations correspond to the strong cross-correlations. e results shown in Figure 3 are consistent with this conclusion and indicate weak cross-correlations between two daily average river levels of these adjacent sites except Site 3-Site 4 and Site 9-Site 10. In fact, we note that there are water storage and discharge activities at Site 2 and Site 7, which can lead to the weak cross-correlations between them and their upstream or downstream (Site 1, Site 3, Site 6, and Site 8). Site 4 and Site 5 are close to Site 7 (see Figure 1) and also affected by this activity. e h XY (q) is usually used to characterize the joint multifractality. us, we conclude that the human activities (water storage and discharge) at Site 2 and Site 7 can change the corresponding joint multifractality of river levels (the part caused by cross-correlations) by weakening cross-   correlations. e significant deviations at Site 8-Site 9 may be due to the lateral inflow from the main tributary at Site 9 which can also weaken cross-correlations. It can be also seen in Figure 3 that the shapes of h XY (q) for all pairs of adjacent sites suggest the two-parameter binomial model. e fitting results are shown in Figure 3 and reported in Table 2. e results show that the h XY (q) given by MF-DCCA is fitted well with two-parameter binomial model (see equation (10)), and the values of R 2 are all larger than 0.95. For two-parameter binomial model, the strength of joint multifractality can be characterized by Δα � ((ln b − ln a)/ln 2). From Table 2, the values of Δα are larger than 0.6, which indicates that there is strong joint multifractality at all adjacent sites.

Source of Joint Multifractality.
In this section, we test the source of joint multifractality obtained from river level records at pairs of adjacent sites in the section above. e testing method is described in Section 2.2. All the testing results are based on 1000 pairs of surrogate series. Figure 4 shows the testing results for the null hypothesis H 0sg1 : the joint multifractality is not due to the correlations within and between series. e shadow area indicates acceptance region of h XY (q) at 5% significance level estimated from 1000 pairs of surrogate series (X sg1 (t), Y sg1 (t)) . For testing of the null hypothesis H 0sg1 , if the joint multifractality is only due to the correlations within and between series, the pair of surrogate series (X sg1 (t), Y sg1 (t)) should show monofractality with h XY (q) � 0.5, and the acceptance region should be around 0.5 for all qs. In this case, our testing method will be invalid for h XY (q) near 0.5 since h XY (q) near 0.5 will always fall into the acceptance region and accept H 0sg1 . is suggests that we should focus on the h XY (q) away from 0.5 when conducting this testing. From Figure 4, triangles indicate h XY (q)s of two preprocessed daily average river levels at adjacent sites. It can be seen that the acceptance region is always around 0.5 and most values of h XY (q) fall out of the shadow area at all pairs of adjacent sites except those near 0.5; that is, we should reject H 0sg . So we can conclude that the joint multifractality is mainly due to the correlations within and between series. Note that most and ((h X (q) + h Y (q))/2) for two preprocessed daily average river levels at adjacent sites. Different colors and symbols indicate the corresponding different results: blue triangles for values of h XY (q) and red circles for values of ((h X (q) + h Y (q))/2). e curves are obtained by fitting the two-parameter binomial model (equation (10)) to points of h XY (q). e resulting model parameters are reported in Table 2. e errors bars of h XY (q) are based on the linear least square fit. All the errors bars reveal 95% confidence interval.
Complexity 7 values of h XY (q) are smaller than the acceptance region, which indicates strong negative correlations within and between records. For further analyzing the effects of cross-correlations on joint multifractality, we obtain the h XYS (q) (see circles in Figure 4) via simultaneously half-shifting technique which is described in Section 2.3 for minimizing the cross-correlations between original series. It can be seen that the difference    Figure 5: Comparison of values of the generalized Hurst exponent h X (q), respectively, obtained from the preprocessed daily average river levels and its half-dividing series. Different symbols indicate the corresponding different results: triangles for results of the preprocessed daily average river levels and circles for its half-dividing series. e errors bars for the preprocessed daily average river levels are based on the linear least square fit. All the errors bars reveal 95% confidence interval.
mainly caused by the autocorrelations within each series; For other pairs of adjacent sites, both autocorrelations and crosscorrelations have significant effects on the joint multifractality. Besides, we also compare the values of the generalized Hurst exponent h X (q), respectively, obtained from the preprocessed daily average river levels and its half-dividing series at all sites for validating the simultaneously halfshifting technique in keeping autocorrelations. e results are shown in Figure 5. It shows that the values of h X (q) change little after the operation of half-dividing. is validates the effectiveness of simultaneously half-shifting technique in keeping autocorrelations. Figure 6 shows the testing results for the null hypothesis H 0sg2 : the joint multifractality is not due to the fatness of probability distribution. e shadow area indicates acceptance region of h XY (q) at 5% significance level estimated from 1000 pairs of surrogate series (X sg2 (t), Y sg2 (t)) . It can be seen that almost all values of h XY (q) fall out of the shadow area; that is, we should reject H 0sg2 . But the difference between h XY (q) and acceptance region is small. is confirms that the probability distribution of river levels also has effects on the joint multifractality, but the effects are small.

Conclusion
e joint multifractal analysis is usually conducted in two different variables for their cross-correlations but rarely used for two records of one variable collected at two different places. Besides, there are also few studies on the source of joint multifractality. In this paper, we focused on the two issues for river level records of a tributary of Haihe River in North China and proposed two hypothesis tests and simultaneously half-shifting technique for the source analysis of joint multifractality. is study is an extension of our previous work on the single multifractality of one record with the same data set. e joint multifractal analysis was conducted using MF-DCCA for two river level records at pairs of adjacent sites along the river. e obtained h XY (q)s was compared with ((h X (q) + h Y (q))/2) at each pair of adjacent sites. Results of comparison showed that human activities (water storage and discharge) can change the joint multifractality to deviate  Figure 6: Testing via pair of surrogate series (X sg2 (t), Y sg2 (t)) . X-axis denotes the range of q values. Triangles indicate h XYS (q)s of two preprocessed daily average river levels at adjacent sites. e shadow area indicates acceptance region of h XY (q) at 5% significance level estimated from 1000 pairs of surrogate series for the null hypothesis H 0sg2 : the joint multifractality is not due to the fatness of probability distribution.
from the average of two single multifractality. e change is due to the weak cross-correlations caused by human activities. It provides an approach for detecting the change in joint multifractality caused by cross-correlation. For the source of joint multifractality, we proposed two hypothesis tests, which are based on empirical distributions of the scaling exponent function estimated from 1000 pairs of surrogate series. It can achieve more robust results than the previous clarification method using one pair of surrogate series. Results of source testings showed that the joint multifractality in river level records is mainly caused by the correlations within and between series and is also related to the probability distribution of river levels. e simultaneously half-shifting technique was proposed for further analyzing the effects of the cross-correlations.
is technique can both minimize the cross-correlations between original series and make full use of records. e further analysis showed more details about the source of joint multifractality. For Site 1-Site 2, Site 2-Site 3, Site 5-Site 6, Site 6-Site 7, and Site 7-Site 8, the joint multifractality is mainly caused by the autocorrelations within each series. For other pairs of adjacent sites, both autocorrelations and cross-correlations have significant effects on the joint multifractality.
ese results are partly consistent with results of joint multifractal analysis. e multifractality of hydrologic dynamics and its source is an important topic in the fields of hydrology and meteorology. e hydrologic complex system has strong nonlinear correlations which cause the multifractal feature of dynamics. e source of multifractality is the direct influence factor of hydrologic dynamics. us, the join multifractal analysis of river levels in this study is useful for the theory and simulation of hydrological phenomena. Especially, the analysis of effects of human activities on joint multifractality can be a reference for the detection and control of human activities in hydrological and meteorological environment. From the research method, this study provides an approach for detecting the change in joint multifractality and some techniques for source analysis of joint multifractality. ese methods can be extended to other applications of joint multifractal analysis.

Data Availability
e data used to support the findings of this study have been analyzed in previous papers [39,40]. Partial records for validation are available at GitHub: https://github.com/ tongzhouzhao/water-level-records.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.