1 INTRODUCTION

The term fractal, according to Mandelbrot, founder of fractal analysis, refers to "a shape made of parts similar to the whole in some way, " in which "similar... in some way" can be self-similarity, self-affinity, or statistical self-similarity, in terms of time or physical space. Fractal phenomena exist widely in nature (e.g., the shapes of coastlines, river networks, and the shapes of trees) and in physics and chemistry (such as fractal noise and soil particle distribution (Li and Zhang, 2005; Kou et al., 2015; Chen et al., 2016a, b). The nature of the fractal analysis is searching for innate law in seemingly unordered, high-dimension and uncertain data by exploring different fractal dimensions. By varying the fractal dimensions, one is enabled to model the dynamic changing law of the environment, and therefore optimize distribution model. Currently (Cai and Gouveia, 2013; Chen et al., 2016c)., the concept and idea of fractal analysis have been abstracted into a methodology, and there have been significant achievements made in data science based on the fractal analysis (Liu et al., 2015; Cao et al., 2016a, b ). Fractal analysis has become a new perspective to study complicated systems and bridges the parts and the overall system (Chen et al., 2008, 2017b; Li and Burgueño, 2010; Ponce-López et al., 2016; Chen and Wang., 2017).

Previously, the method for calculating wave heights of a certain return period was to choose a probability distribution model (Gumbel and MaxEntropy Distribution) as the probability distribution for annual extremes of wave heights (Wang et al., 2013a). The observed annual extreme wave heights were used to determine the parameters of the probability distribution, after which a cumulative rate determined the wave height that only occurs once in many years (a century or more) (Wang et al., 2016, 2017). Such a calculation method has three drawbacks: 1) it assumes that all annual extreme wave heights follow the same probability distribution while disregarding whether the time period is short term or long term (Wang et al., 2013b; Chen et al., 2017a); 2) the statistical characteristics of the data obtained over a short period of time and those obtained over a long period (a century or more) are both strictly selfsimilar(Wang et al., 2014); and 3) during the calculation, only annual extreme wave heights are used (only one data point each year) while most observational data are disregarded (Wang et al., 2015). All of these drawbacks lead to unconvincing calculation results. In fact, the long-term evolution of wave heights at an observation point is a very complex system, and only by conducting multifractal analysis to reveal the relationship between the short-term observational wave height sequences (parts) and the long-term observational wave height sequence (overall) can the latter be inferred from the former (Kantelhardt et al., 2002; Lu, 2004; Yang et al., 2018), thus providing a good foundation for more reasonable estimations of wave heights that only occur once in many years.

We apply multifractal analysis for the first time in the characterization of a wave height sequence observed in marine and hydrological stations, by conducting a comprehensive analysis on a part of the observed data using methods of statistical physics, from which the characteristics of the whole dataset can be calculated. Currently, two multifractal analysis methods have been proposed: the partition function method (Jiang and Zhou, 2008) and the multifractal detrended fluctuation analysis (MF-DFA) method (Hu and Song, 2003; Wang and Tian, 2006; Xu, 2011).

The detrending of fluctuations is an important step of the MF-DFA method; however, currently, there are some drawbacks associated with the calculations performed during this step: 1) the fitted polynomial function becomes discontinuous at the connection points of adjacent segments, which leads to the generation of new pseudo-fluctuation errors (Zhuang, 2012); 2) the selection of the polynomial order is highly subjective—low-order polynomials cannot accurately represent data fluctuations, whereas higher order polynomials can cause overfitting problems. These drawbacks will lead to erroneous analysis results on observed data.

In this study, we use the method of signal decomposition (Cai et al., 2011a, b, 2016; Xu and Wang, 2011; Chen and Wang, 2012; Hu et al., 2013; Ding and Fan, 2015, 2016; Escalante et al., 2016) to obtain a more reasonable fluctuation function. An improvement has been made to the detrending procedure of the MF-DFA method, which is described in Section 2. In Section 3, we use the MF-DFA method and its improved version to analyze the wave height sequence (significant wave heights) of Chaolian Island, Shandong, China, waves measured over a 27- year period. We make a comparison between the analysis results, including a detailed comparison between the fluctuation function obtained from polynomial fitting and that obtained from mode decomposition, which demonstrates the advantages of the improved method.

2 IMPROVEMENT ON THE CONVENTIONAL MF-DFA METHOD

The MF-DFA method is effective in verifying the multifractality of a non-stationary time series, which describes the multifractality of the object using the generalized Hurst exponent H(q). A key step of the MF-DFA method is to calculate the trend function in each segment, determining the fitting polynomial function using the least-squares method. The degree of the polynomial determines the type of detrended fluctuation, which leads to several problems: 1) thefi tted polynomial function is discontinuous at connection points between adjacent segments, which leads to the generation of new pseudo-fluctuation errors (Koag et al., 2014; Wei et al., 2017; Zhang et al., 2017b, 2018; Barrs and Chen, 2018; Bhimani et al., 2018); 2) the fitted polynomial can be linear, quadratic, or of higher orders (Liu et al., 2015; Zhang and Kleit, 2016; Xu et al., 2017a, b ; Zhang et al., 2017a; Zhe et al., 2017). Lower order polynomials cannot accurately represent the trend in data fluctuation, whereas higher order polynomials can cause overfitting problems. Selection of the specific type depends on the specific issues being studied, and this is a highly subjective procedure.

Based on the previous problems, this paper proposes using mode decomposition of signals to solve the trend function, which can better meet the requirements of solving fitted trends while avoiding the disadvantages of the conventional MF-DFA method.

2.1 Definition of multifractal and an improved MF-DFA method

Multifractal can be described by generalized fractal dimension or the spectral function; the definition is as following:

Definition 1: Suppose the function sequence of time{X(t)}is a stable stochastic process, if it satisfies:

We say that it has multifractal, of which t∈T, 0∈T; q∈Q, [0, 1]∈Q; T, Q is a real interval with positive length; τ(q), c(q) is a function on Q, Δt is the increment in time, τ(q) is the multifractal mass index.

The two key parameters of multifractal are the generalized Hurst index and the mass index respectively.

Definition 2: Generalized Hurst index H is defined as:

Definition 3: Definition of multifractal spectrum f(α): If the fractal object is divided into N intervals with ε as the linearity of each interval, the probability density of the i-th interval shall be:

Of which α is called singularity index. Let the number of intervals with the same α value be N α(ε), then

f(α) is the multifractal spectrum, indicating fractal characteristics of sub-intervals with different α values within the fractal object.

Definition 4: Mass index is defined as:

Of which χ_q(ε) is the partition function, indicating the weighted summation of probabilities of all intervals:

Zhang (2006) adapted the fast-bandpass signalfiltering method proposed by Xu and Wang (personal communication; published in 2011) to the empirical mode decomposition of signals, establishing a fast filtering method for signal mode decomposition that provided good results. We use the latter method to solve the trend function signals. A detailed description of this method can be found in the references (Zhang 2006, 2007; Xu and Wang, 2011) and a brief introduction is given here.

Discrete Fourier transform is conducted on the wave height sequence {x_n}, resulting in the sequence {x_m}. The pre-determined upper and lower limits of the bandpass frequency are ω_m1 and ω_m2, respectively, both of which are converted to m₁ and m₂ according to m=NΔtω_m/(2π). The sequence {W_m} is generated as follows:

When an inverse Fourier transform is performed on the sequence {W_m}, the real part of the resulting sequence is the signal {x_n^(k)} of the filtered bandpass.

In Eq.7, take m₂=N/2 and m₁=m₂–1, m₂–2, …, and calculate the filtered signal {x_n⁽¹⁾} for different values of m₁. In addition, every calculation result is examined for being an intrinsic mode function (IMF). The procedure is repeated until an m₁ value is taken that makes the signal filtered through the bandpass (ω_m1, ω_m2) an IMF while the signal filtered through the bandpass (ω_m1–1, ω_m2) is no longer an IMF. In this way, the first IMF {x_n⁽¹⁾} is decomposed from the wave height sequence {x_n}.

Similarly, take m₂=m₁–1, and repeat the previous procedure, and the second IMF {x_n⁽²⁾} can be decomposed.

The same procedure is repeated until the k^th IMF is decomposed from the wave height sequence {x_n}, in which the last signal through the last bandpass (0, ω_k–1) is a monotonic function {M_n}. In this way, the signal {x_n} is decomposed into multiple IMFs and one monotonic function, as shown below:

where are all IMFs, and {M_n} is the trend function that represents the overall trend of the signal.

2.2 An improved MF-DFA method

The improvement focuses on replacing the segmentation procedure of the MF-DFA method for fitting polynomial trends with mode decomposition of trends, while the other steps remain essentially the same.

For a wave height sequence {x_n} of a given length N, the major steps of the improved MF-DFA method are as follows:

1. Calculate the cumulative deviations {y_n} of the wave height sequence:

2. Calculate the mode trend function of the cumulative deviations and segment this cumulative deviation sequence as well as the trend function: the fast filtering method of signal decomposition is used to calculate the mode trend {T(n)} of the cumulative sequence {T(n)}. In order to maintain the remaining part, segmentation is conducted from opposite ends, dividing the sequences {y(n)} and {T(n)} into 2Ns segments of equal length s, in which Ns=[N/s] (rounding up the result).

3. Calculate the q^th order fluctuation function: first, calculate the mean of the squared residuals, i.e.,

Therefore, the q^th order fluctuation function is

4. Calculate the generalized Hurst exponent: the relationship between the fluctuation function F_q(s) and the length s of each segment is as follows:

In the expression, H(q) is called the generalized Hurst exponent. For each value of q, the corresponding H(q) can be obtained from the linear fitting of the lnF_q(s) vs lns plot.

When the value of the generalized Hurst exponent H(q) is independent of q, the wave height sequence is a monofractal series; when the value of H(q) changes with q, the sequence is a multifractal series.

3 EMPIRICAL ANALYSES

In this section, we analyze the observational data from the Chaolian Island Hydrological Station for the wave height sequence (significant wave heights) using the partition function method, the traditional MF-DFA method, and the improved MF-DFA method.

3.1 Data information

The data in this study are the wave height sequence (significant wave heights), measured four times a day in the Chaolian Island Hydrological Station from 1963 to 1989. Data are missing or insufficient in some years. Figure 1 shows the scatter plot of the original wave heights, illustrating that starting from 1963, there are some data anomalies in the first few years. In the subsequent discussion, the self-similarity and scale invariance of wave height through all the measured data are discussed. Based on a shape made of parts similar to the whole in some way, when extrapolating the design wave height, the single extremal data over a certain threshold will be selected to estimate the parameters of the distribution function.

To prevent interference from the abnormal data, the data of the first three years were excluded from the original data. For the sake of the simplification of discussions, this paper selects the wave height over the given threshold. If other thresholds are selected, it does not affect the self-similarity of the wave height data. As the scatter plot illustrates, the observed data in the first three years were belong to anomalies. Therefore, the abnormal data were excluded from the subsequent discussion and the data of the later 23 years were selected. In addition, since the research focuses on larger wave heights, a threshold selection was performed. A total of 13 546 effective wave height data whose threshold values were above average were used for the analysis, as shown in the scatter plot in Fig. 2. As the figure illustrates, the overall wave height sequence is statistically similar to its local parts, and therefore it can be further analyzed using fractal theory.

3.2 The fractal characteristic analysis method

Using a mass exponent τ(q) and multifractal spectrum f(α), we describe the characteristics of the normalized wave height sequence in intervals of different scales, to demonstrate the fractal properties of the sequence (Fig. 3).

First, the lnX_q(s) vs ln(s) relationship was examined. Figure 3a shows the double logarithmic plots of the sequence for different values of q, in which a total of 21 integer values ranging from -10 to 10 are taken for q. Figure 3a shows only the patterns for even values of q. The data points produce a nearly straight line regardless of the value of q, suggesting that for any given q, the wave height sequence has fractal scaling properties. The figure also shows that for different values of q, slopes of the lines in which the data points lie are different from each other, suggesting a multiscale property of the sequence; i.e., the wave height sequence of Chaolian Island has a multifractal characteristic.

Table 1 shows the differences between the value of τ(q) and that of the corresponding straight line for different values of q. Combined to Fig. 3b, it is found that the τ(q) vs q curve deviates slightly from a straight line, with the total sum of squared deviations reaching 1.650 7, implying insignificant multifractality of the wave height sequence.

When Fig. 3c and Table 2 are combined, the multifractal spectrum f(α) mainly appears as a hookshaped quadratic curve to the left. The Holder singularity exponent α is featured with a scale range of (0.948, 1.093)—the entire scale range is larger than 0.5, while Δα is only 0.145, and the maximum value of f(α) is 1.004, suggesting a relatively narrow distribution of the multifractal spectrum curve. This result not only demonstrates the presence of longterm memory behavior in the fractal pattern of this wave height sequence but also depicts its multiscaling property implied by the variability in the scaling exponent versus volatility, both of which further indicate the relatively weak multifractal characteristic of the wave height sequence at Chaolian Island.

3.3 A comparison between the analytical results obtained by the conventional MF-DFA method and those by the improved MF-DFA method

A comparison is made between the conventional MF-DFA method and the improved MF-DFA method on the basis of two aspects: trend-fitting effect and final pattern.

The trend-fitting effect obtained by mode decomposition is the overall trend-fitting result for the cumulative deviation sequence. For comparison, third-order polynomials fitted to the cumulative deviation sequence were added. In the conventional MF-DFA method, trend fitting for cumulative deviations is conducted using piecewise polynomials, normally with a first-, second-, or third-order polynomial trend that depends on the volume of data. In this case, the cumulative deviation data were divided into segments, each with 50 data points, and a second-order polynomial was fitted to each segment (Fig. 4).

As Fig. 4 shows, the third-order polynomial fitting is too general compared to that obtained by mode decomposition and to that by piecewise polynomial fitting—it cannot accurately represent the localfl uctuations of the cumulative deviations of wave heights, and therefore, it provides the worst fit. To further compare the results obtained by mode decomposition with those by piecewise polynomial fitting of the MF-DFA method, the region containing data numbers 4 960–5 160 in Fig. 4 is enlarged and shown in Fig. 5.

Figure 5 can be analyzed from two perspectives. First, within a complete segment—data points 5 001–5 050, for example—the trend obtained by mode decomposition is smoother than the MF-DFA piecewise polynomial trend, and it works better to depict the projections and declines of the fluctuations, implying a better fitting effect; second, at the connection points (i.e., data points 5 001 and 5 050), there is a large jump present in the curve of the MFDFA method, which means that the fitting is discontinuous at connection points of adjacent segments, which generates new pseudo-fluctuation errors that interfere with the subsequent calculations. However, such a discontinuity does not occur in the mode decomposition method, and therefore no pseudo-fluctuation error is generated.

Next, we examine the detrending results for the wave height cumulative deviations (Fig. 6). Figure 6a shows the normalized cumulative deviations, and Figure 6b shows the cumulative deviations after the overall third-order polynomial function is detrended. It reveals that the goal of detrending is not achieved. Figure 6c and d shows the cumulative deviation results obtained by detrending the piecewise polynomial as well as mode decomposition, respectively, both of which reach the goal of detrending. In addition, as the figure illustrates, because the trend obtained by mode decomposition implies a better fitting effect at the connection points, when the long-term measured data and the short-term measured data are nested each other, the advantage of the fitting, based on the signal mode decomposition method, will more evident with the increase of shortterm nested data nested. The means and variances of the cumulative deviations in the latter two groups are calculated for further analysis (Table 3).

As shown in Table 3, both means are close to zero. Since both methods need to solve the fluctuation function after detrending the cumulative deviations, the mean value does not have a significant effect on the succeeding calculations; the major impact is the volatility of the cumulative deviations. Variance reflects data dispersion, i.e., the magnitude of volatility. The variance of the first method is five times that of the latter, showing that the dispersion of the detrended cumulative deviations obtained by mode decomposition is smaller, with smaller fluctuations, and therefore that this trend tends to stabilize more easily.

In the previous discussions, the advantages and disadvantages of the MF-DFA method and its improved version in terms of fitting effect have been discussed; it is demonstrated that the latter can achieve a better fitting effect while avoiding some drawbacks of the original method.

In the following discussion, we further analyze the advantages and disadvantages of these two methods by comparing the final graphs. For both methods, a total of 21 integer values from -10 to 10 were set for the parameter q.

Figures 7 and 8, respectively, show the lnF_q(s) vs lns double logarithmic plots obtained with the MFDFA method and with the improved MF-DFA method. Figures 7a and 8a show the ln F_q(s) vs lns double logarithmic plots, and Figs. 7b and 8b show the corresponding linear curve fitting results. When comparing the two figures, we see that the double logarithmic plots obtained with both methods are of good linearity, indicating that the wave height sequence is long-range power-law correlated and thus that it has a multifractal characteristic, which is consistent with the results obtained by the partition function method. In addition, both of the wave height fluctuation functions obtained tend to converge at some level, the difference being that the former tends to converge slowly whereas the latter has a faster convergence speed; its fluctuation becomes smaller and smaller and therefore tends to stabilize more easily. This is inconsistent with the results shown in Table 3, which further corroborates the good performance of the improved method.

Figures 9 and 10 show the curves of the generalized Hurst exponent H(q) calculated using the two methods. When comparing the two figures, we see that the generalized Hurst exponent H(q) decreases as the partition order q increases. According to the theory, this suggests that the value of H(q) actually changes as q changes and that the wave height sequence has a multifractal characteristic.

Meanwhile, it is roughly shown that the H(q) curve of the improved method has a more marked inflection point compared with the original method—the inflection point occurs at q=-3 and decays much faster.

As Table 4 shows, for the same set of q values, the fluctuation ΔH(q) of H(q) calculated with the original method is only 0.284 3, while that of the improved method reaches 0.694 3, which is almost 2.5 times that of the original method. This further confirms that the improved MF-DFA method has a faster decay rate as well as better stability.

4 CONCLUSION

The multifractal spectrum f(α) of Chaolian Island appears as a hook-shaped quadratic curve to the left and its Holder singularity exponents α obtained are all larger than 0.5. However, the entire scale range is relatively small: Δα is only 0.145, and the maximum value of f(α) is 1.004. The multifractal spectrum has a relatively narrow distribution, showing weak multifractality in the wave fluctuations of Chaolian Island.

In this paper, we have applied fractal theory to the analysis of the fluctuation characteristics of waves. A partition function-based multifractal analysis of the observational wave height data of Chaolian Island was conducted, and it was observed that the multifractal spectrum f(α) appears as a hook-shaped quadratic curve to the left. In addition, the Holder singularity exponents α obtained are all larger than 0.5. However, the entire scale range is relatively small: Δα is only 0.145, and the maximum value of f(α) is 1.004. The multifractal spectrum has a relatively narrow distribution, showing weak multifractality in the wave fluctuations of Chaolian Island.

In addition, in this paper, the signal mode decomposition method is proposed to solve the trend function against the existing problems of the piecewise polynomial method used during the detrended fluctuation analysis of the MF-DFA method. An improved MF-DFA method is established based on mode decomposition, which is applied to the analysis of the wave height observational data from Chaolian Island. A comparison with the conventional MF-DFA method shows that:

(1) Compared to the original method, the improved method achieves finer trend fitting, enabling it to better depict the projections and declines of the fluctuations, thus achieving a better fitting effect. In addition, the use of the improved method avoids the problems of the original one wherein the fitted trend is discontinuous at connection points of adjacent intervals; thus, it does not generate a new pseudofluctuation error, which is beneficial to the succeeding calculations of fluctuation functions.

(2) The cumulative deviations of the wave height data have a smaller dispersion after the detrending procedure by the improved method and thus have smaller fluctuations—the variance obtained is only 1/5 that of the original method, suggesting that the improved method can better meet the requirements of detrending fluctuations.

(3) A comparison between the double logarithmic plots lnF_q(s) vs lns obtained with the two methods shows that both methods can achieve good linearity, indicating that the wave height sequence is longrange power-law correlated and thus has a multifractal characteristic, which is consistent with the results obtained by the partition function method. However, the wave height fluctuation function calculated with the improved method has a more obvious convergence tendency, with a faster convergence speed, and thus tends to stabilize more easily.

(4) The generalized Hurst exponent H(q) curve calculated with the improved method has a more marked inflection point, which occurs at q=-3. Under the same set of q values, the fluctuation range ΔH(q) of H(q) calculated with the original method is only 0.284 3, while that of the improved method reaches 0.694 3, nearly 2.5 times that of the original method, implying that the generalized Hurst exponent has a faster decay rate, thus indicating better stability of the improved method.

In this paper, we introduced fractal to explore the fluctuation characteristics of sea waves. The introduction of the method can better analyses the internal statistical characteristics of waves. By discussing the statistical self-similaritities of the multifractal and the probability distribution models known (Pearson-Ⅲ, Gumbel Distribution), we can not only calculate the design wave height and the design water level, but also get more details about the internal statistical characteristics of waves.The research work in this paper facilitates to accurately classify the disaster grades of storm surges. By discussing the statistical characteristics of dense discrete data, it provides feasible theoretical support for extrapolating the design wave height with discrete extremum data.

5 DATA AVAILABILITY STATEMENT

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.