Multi-scale fluctuation analysis of precipitation in Beijing by Extreme-point Symmetric Mode Decomposition

With the aggravation of the global climate change, the shortage of water resources in China is becoming more and more serious. Using reasonable methods to study changes in precipitation is very important for planning and management of water resources. Based on the time series of precipitation in Beijing from 1951 to 2015, the multi-scale features of precipitation are analyzed by the Extreme-point Symmetric Mode Decomposition (ESMD) method to forecast the precipitation shift. The results show that the precipitation series have periodic changes of 2.6, 4.3, 14 and 21.7 years, and the variance contribution rate of each modal component shows that the inter-annual variation dominates the precipitation in Beijing. It is predicted that precipitation in Beijing will continue to decrease in the near future.

1 Introduction

Beijing is one of the cities with the most serious water shortage in China. In the context of global warming, the annual precipitation is significantly reduced, which aggravates the water shortage crisis in Beijing. Therefore, studying the changing law of precipitation and the trend of its future development in Beijing has become a central issue of widespread concern. At present, many scholars have done a lot of researches on the variation and evolution of precipitation in Beijing. Zhu et al. (2012) used Empirical Mode Decomposition method to analyze the multi-scale oscillations of the precipitation time series from 1951 to 2009 in Beijing, the results show that the annual precipitation in Beijing will decrease continuously in the short term (Zhu et al., 2012). Through statistical calculations and frequency analysis of precipitation data from 16 rainfall stations in Beijing, Sun et al. (2007) explored the spatial and temporal distribution characteristics of precipitation in Beijing and its changing trend. It is found that continuous drought in Beijing after the 1990s and annual precipitation after 2008 may enter a relatively abundant period (Sun et al., 2007). Wang et al. (2009) analyzed the changes of temperature and precipitation over the past 48 years in Beijing and found that the annual precipitation showed a decreasing trend in both suburban and urban areas. The above research results have analyzed the cycle and variation characteristics of temperature and precipitation in Beijing from different perspectives, which is not only helpful for understanding the process and laws of climate change, but also great significance for the social and economic development of Beijing.

The climate system is a nonlinear, non-stationary, and hierarchical system (Xue et al., 2013). In the traditional trend analysis, large-scale cycles and trend changes are likely to mix together, can not tell the trend of changes or periodic oscillations. Extreme-point Symmetric Mode Decomposition (ESMD) is a locally adaptive time series analysis technique developed in recent years. Compared with the wavelet analysis method, it is free from the constraints of Fourier transform and is very suitable for analyzing non-stationary and nonlinear time series (Xia and Liu, 2017). This paper uses ESMD method to conduct multi-scale analysis of the annual precipitation series from 1951 to 2015 in Beijing, and compares the results with wavelet analysis to provide reference for future precipitation forecast and water resources planning and management in Beijing.

Figure 1 Annual precipitation time series from 1951 to 2015 at Beijing station.

2 Research methods

2.1 The basic principle of Empirical Mode Decomposition

The Empirical Mode Decomposition (EMD) method proposed by Huang et al. (1998), which is very suitable for dealing with nonlinear and non-stationary signals. The EMD method separates the fluctuations of different periods from the original signal and finally obtains the trend component. Fluctuations at different scales are defined as Intrinsic Mode Function (IMF). The modal decomposition of the EMD method is as follows:

1. Find all maxima and minima of the sequence x(t), through the cubic spline interpolation, the upper envelope (u₁(t)) and the lower envelope (u₂(t)) of x(t) are calculated.

2. Calculate the median curve (m₁(t)).

   $$\begin{array}{}\text{(1)}& m_{\mathrm{1}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(u_{\mathrm{1}}\left(t\right)+u_{\mathrm{2}}\left(t\right)\right)\end{array}$$

3. Subtracting the median curve from the original signal yields the remaining signal.

   $$\begin{array}{}\text{(2)}& h_{\mathrm{1}}\left(t\right)=x\left(t\right)-m_{\mathrm{1}}\left(t\right)\end{array}$$

4. Using h₁(t) as a new signal sequence, repeat the above operation until h₁(t) becomes a zero-mean process.

5. After the zero-mean process h₁(t) is obtained, it is taken as the first IMF component c₁(t), It represents the highest frequency component of the original signal.

6. Subtracting c₁(t) from the original signal x(t) results in a new signal process r₁(t).

   $$\begin{array}{}\text{(3)}& r_{\mathrm{1}}\left(t\right)=x\left(t\right)-c_{\mathrm{1}}\left(t\right)\end{array}$$

7. The sequence r₁(t) is taken as the original data, repeat steps (1)–(5) and obtain the sequence c₂(t), c₃(t), c₄(t) and so on. The original signal is finally reconstructed.

   $$\begin{array}{}\text{(4)}& x\left(t\right)=\sum _{j=\mathrm{1}}^{n-\mathrm{1}}{c}_j\left(t\right)+r_n\left(t\right)\end{array}$$

   Where: The sequence r_n(t) is the trend term, representing the average trend of the signal.

Figure 2 Beijing annual precipitation model ratio coefficient curve.

2.2 The basic principle of Extreme-point Symmetric Mode Decomposition

The Extreme-point Symmetric Mode Decomposition (ESMD) method is a new development of the EMD method. This method solves the problem of “modal aliasing” in EMD and it has a big advantage in climate data analysis (Wang and Li, 2014). This method is good at finding trends and can separate inter-annual trends and the general trend of climate change in the observed sequence, helping to probe the issue of global warming. And this method can find the anomalous period and frequency band from the decomposition mode, which is good for the research of climate anomaly. The modal decomposition of the ESMD method is as follows:

1. Find all maxima and minima of the sequence X(t), denoted by $E_i(\mathrm{1}\le i\le n)$.

2. Connect all the adjacent poles with line segments and mark the mid-points as $F_i(\mathrm{1}\le i\le n-\mathrm{1})$, and add the boundary mid-points F₀ and F_n at the left and right ends.

3. The interpolation curve L₁, …, L_p(p≥1) is constructed by using n+1 midpoints, and calculated the mean curve $L^{*}=(L_{\mathrm{1}}+\mathrm{\dots }+L_p)/p$.

4. The above steps are repeated for the $X\left(t\right)-L^{*}$ sequence until $\left|L^{*}\right|\le \mathit{\epsilon }$ (ε is the allowable error) or the number of screened times reaches the preset maximum value K, resulting in the first empirical mode M₁(t).

5. The above four steps are repeated for the remaining sequences X(t)−M₁(t), and the empirical modes M₂(t)−M₃(t) … are obtained respectively until the remaining sequence R(t) remains only for a certain number of poles.

6. Change the maximum number of screening times K within the defined interval [$K_{min},K_{max}$] and repeat the above five steps. Then calculate the variance ratio σ∕σ₀, and draw it with K changes in the map to find the minimum σ∕σ₀ corresponding K₀. Repeat the above five steps again with K₀ as the limiting condition. The last remaining term R(t) is the adaptive global averaging of the sequence X(t).

After decomposition, the original time series X(t) is reconstructed as:

$$\begin{array}{}\text{(5)}& {\displaystyle }X\left(t\right)=\sum M_i\left(t\right)+R\left(t\right)\end{array}$$

3 Case study

3.1 Study area profile

Beijing is located in the northwestern end of the North China Plain, the central Haihe River Basin, the land area of 16 410 km², is a typical northern temperate semi-humid monsoon climate. The complex and varied topography has created a diversity of weather conditions in Beijing. Summer heavy rainfall and strong convective weather are the main causes of frequent droughts and floods in Beijing. Shown in Fig. 1 for the Beijing 1951–2015 nearly 65 years of precipitation. Among them, the largest precipitation was 1407 mm in 1957, the smallest precipitation was 267 mm in 1965, the average annual rainfall was 589.75 mm, and the rate of linear precipitation decrease is 3.47 mm a⁻¹.

In this study, the annual rainfall data are provided by the National Climate Center of China Meteorological Administration. The data are released after the reorganization and review of the institution, so the reliability of the data can be basically guaranteed. As shown in Fig. 2, plot annual precipitation model ratio coefficient curve in Beijing, and analyze the representativeness of the data. It can be seen from the figure that with the increase of time, the amplitude of the curve is getting smaller and smaller, and when the annual rainfall series is over 50 years, the cumulative mean of the modulus ratio tends to 1, indicating that this series of data has certain stability.

3.2 Multi-scale analysis of precipitation changes

The ESMD method was used to decompose the annual mean precipitation series from 1951 to 2015 in Beijing. The best screening frequency was 26 times when the ratio of variance was the smallest, and 4 modes (Mode 1–4) and a trend margin R, as shown in Fig. 3.

Figure 3 Modal components and trend of precipitation time series by ESMD

In order to verify the reliability of the decomposition results, this paper synthesize a reconstructed sequence from the decomposed Mode 1–4 and the trend series R, and find that the reconstructed sequence is completely consistent with the original precipitation sequence, indicating that the decomposition result of ESMD method is credible of. In Fig. 3, the four mode components of ESMD decomposition each reflect the oscillations of different characteristic scales inherent in the original sequence, and the influence of signal fluctuation frequency of each scale on the overall characteristics of the original data is represented by the variance contribution rate. The Periodic Diagram Method is used to estimate the average period of each decomposition mode. Table 1 shows the mean period and the variance contribution rate of each modal component. It can be seen that annual precipitation in Beijing has annual variations of 2.6 and 4.3 years, and inter-decadal changes of 14 and 21.7 years. In addition, the main periods of the precipitation time series does not change over time during the study period.

Combined with Table 1 and Fig. 3 quantitative analysis: The precipitation series of 65 years in Beijing contain many time-scale features, and the precipitation changes are mainly determined by the three higher-frequency oscillations of Mode 1–3. Among them, the Mode 1 (Fig. 3a) with a period of 2.6 years has the largest variance contribution rate, reaching 32.47 %. The signal oscillation of Mode 1 is very obvious and indicates the precipitation has a decreasing-increasing cycle change, which basically reflects the alternation of precipitation in Beijing before the 1970s. Mode 2 (Fig. 3b) indicates that there is a 4.3-year periodic change in the precipitation series in Beijing. Its variance contribution rate is about 16.74 %. Throughout the study period, Mode 2 shows that the amplitude of rainfall in Beijing fluctuated greatly before the 1970s and gradually decreased after the 1970s. Mode 3 (Fig. 3c) with a period of 14 years has only a smaller variance contribution rate of Mode 1, about 27.5 %. As can be seen from Fig. 3c, the amplitude changes steadily and fluctuates relatively little in the late 1960s and mid-1990s, but it significantly increased from the early 1950s to the late 1960s and after the 1990s. This shows that the precipitation variation in these two periods is great compared with other periods. Mode 4 (Fig. 3d) with a period of 21.7 years has the smallest variance contribution rate, only 1.37 %. Its amplitude is relatively stable over the entire time span. Comparing Mode 3 with Mode 4 can find that Mode 3 has a wide range of amplitude variations and the amplitude range of Mode 4 is significantly smaller. This indicates that the anomalous precipitation mainly occurs on the time scale of 14 years.The variance contribution rate of trend item R (Fig. 3e) is about 16.74 %. It reflects the overall trend of precipitation in Beijing over the past 65 years. As shown in Fig. 3e, trend item show a non-linear trend of gradual decline over the entire time scale.

3.3 Wavelet analysis and verification

Wavelet analysis is a hot frontier field in recent years because of its special advantages for signal processing, it is widely used in the field of time-frequency structure analysis of weather and climate sequences (Yang and Shu, 2017). This paper uses the Morlet wavelet transform to analyzes the annual precipitation anomaly sequence. Figure 4 shows the real part time frequency distribution of wavelet coefficients of the precipitation anomaly sequence in Beijing in recent 65 years. When the wavelet coefficient is positive, it indicates that the precipitation is relatively large and the greater the value, the more precipitation. When the wavelet coefficient is negative, the precipitation is less than normal.

Table 1 Contribution rates of ESMD decomposition for mean precipitation.

Table 2 Comparison of the cycle result by ESMD and Wavelet Analysis.

Figure 4 The real part distribution of wavelet coefficients.

It can be seen from Fig. 4 that precipitation in Beijing has multi-scale changes of 2, 5–8, 10–15 and 25 years. Throughout the study period, the wavelet coefficients alternated between positive and negative, indicating that the change of precipitation with the alternation of abundance and dryness over time. Table 2 shows the comparison between the wavelet and the ESMD analysis cycle results, the results of the two methods have some differences, and the longer the cycle shows the greater the difference.

As a time-frequency analysis method, wavelet transform is the same as other time-frequency analysis methods, such as short-time Fourier transform and classical spectral estimation, which is based on Fourier transform theory as a basis. There are some limitations, such as the selection of wavelet basis functions, constant multi-resolution and other issues. Compared with the wavelet analysis, the ESMD method based entirely on the characteristics of data changes decomposition and got rid of the shackles of the Fourier transform theory. This method not only has the advantages of multi-resolution wavelet transform, but also breaks through restrictions of the wavelet basis functions, it has strong flexibility and adaptability. The decomposition process of this method is easier than wavelet analysis, and the components can clearly characterize the variation characteristics of signals on different time scales. It is more helpful to explore the inherent laws of things.

4 Conclusion

ESMD method is a suitable method for analyzing non-linear and non-stationary sequence signals. Applying the ESMD method to climate element time series not only separates the annual and inter-decadal trends, but also separates the general trend of climate change. This study shows that the 65-year precipitation series in Beijing are mainly composed of Mode 1, Mode 2 and Mode 3, indicating that the 2.6-, 4.3-, and 9-year scale oscillations play a major role in the change of the whole sequence. Combining several mode components and the trend term R predicts that the annual precipitation in Beijing will continue to decrease in the next few years.

Under the background of global warming, the influence of natural factors and human activities are superimposed on each other, and the variation of hydrological time series is more complicated. Based on the study of precipitation series in Beijing, this paper shows that ESMD method has better analysis effect on time series of hydrology. However, the current research is only in the initial stage. How to determine the influence factors of each mode component, calculate its weight and finally apply it to the prediction of time series will be the key work in future.