Annual precipitation analysis and forecasting – taking Zhengzhou as an example

China's total water resources are less than 3 trillion m³, while the per capita water resources are only 0.22 million m³, only 0.25 times of the global average (Li 2014). However, Zhengzhou has a large population base, and the total amount of water per capita is only 0.2 times of the average level of China, belonging to the area of extreme shortage of water resources. Precipitation is an important source of regional water resources and precipitation anomaly is the direct driver of regional flood and drought disasters. Analysis of precipitation periodic characteristics and precipitation forecasting can provide data support for rational utilization of regional water resources, flood control and drought relief, and have important guiding significance for social and economic development.

Many eastern and western scholars have carried out annual variation analysis of precipitation series in various regions, using the wavelet analysis method (Duan et al. 2018; Partal 2018; Zhao et al. 2019). The wavelet analysis method can not only clearly show the intensity distribution and precipitation change trend of various time-scales, but also analyze the main precipitation cycle of the region. In this paper, the wavelet transform is used to analyze the annual precipitation in Zhengzhou, and thus the periodic variation characteristics of precipitation in Zhengzhou were obtained.

Wavelet analysis mainly reflects the change of precipitation period on time-scales, but it cannot reflect the specific changes in the next few years. In order to further reflect the precipitation in the next few years, it is necessary to make a precise forecast of precipitation. Precipitation forecasting at home and abroad mainly focuses on remote sensing observation and mathematical models. In terms of remote sensing observation, the studies of Ganguly & Bras (2003), Diomede et al. (2008) and Pombo et al. (2015) are the most representative. In terms of mathematical models, many scholars mainly focus on the study of precipitation in the aspects of artificial neural networks (Teng et al. 2016; Boonyuen et al. 2019), regression analysis (Nayagam et al. 2008; Nicholson 2015) and combined models (Spiridonov et al. 2017; Xu & Wang 2019). However, according to the characteristics of the precipitation sequence itself, it is rare to build a coupled forecasting model by reducing the non-stationarity of the sequence. Precipitation evolution is a complex multidimensional system driven by many uncertain factors (Song et al. 2018; Li et al. 2019). The development trend and direction of regional precipitation forecasting are combined with modern nonlinear theory, which can effectively improve the forecasting accuracy. Complementary Empirical Mode Decomposition (CEEMD) can decompose non-stationary time-series into Intrinsic Modal Function (IMF) components and trend, and the amplitude, fluctuation and trend of the IMF components are smaller than the original sequence. Thus, the non-stationarity of the sequence can be reduced. Combined with the inherent advantages of the Elman neural network in dealing with nonlinear problems (Liu & Zhu 2019), the Autoregressive Integrated Moving Average model (ARIMA) has a strong approximation ability to stationary sequences (Liu et al. 2017). The CEEMD − Elman + ARIMA model is constructed, which provides a new method for precipitation forecasting.

Morlet in 1984 proposed wavelet analysis with time-frequency multi-resolution function, which provided the possibility of better research on time-series problems (Goupillaud et al. 1984). It can clearly reveal various change cycles hidden in time-series, fully reflect the change trend of the system under different time-scales, and make qualitative estimation of the future development trend of the system.

In accordance with Equations (3) and (4), the basic principle of wavelet analysis is to obtain signal information at various frequencies by adjusting the extension scale a. Then the local information of the signal is analyzed so as to realize the analysis of time-scales and local spatial characteristics of different signals.

The change process of the wavelet squared with scale a is called the wavelet variance diagram. According to Equation (6), it can reflect the energy distribution of signal fluctuation with scale a. Therefore, the wavelet variance diagram can be used to determine the relative intensity of different scale disturbances in the signal and the main time-scale, namely the main period.

The wavelet is a kind of continuous complex wavelet whose real part reflects the time distribution and phase information of the signal at different scales. In order to simplify the calculation, the precipitation time-series was processed centrally. The centralization process (⁠⁠; ⁠) and Morlet wavelet were substituted into Equation (4), and the value of a and b to calculate the wavelet transformation coefficient ⁠.

It can be seen that the signal after CEEMD decomposition is complete, and the accuracy of the reconstructed signal is also guaranteed.

The Elman neural network was proposed by Elman in the 1990s (Elman 1990), the network is a kind of feedback neural network composed of input layer, hidden layer, acceptor layer and output layer. Different from the BP neural network, the Elman neural network adds a connective layer in the hidden layer as a delay operator to realize the dynamic memory of the system. Because of the good memory and stability characteristics of the Elman neural network, it has been widely used in various fields (Jiang et al. 2016; Yates & Keedwell 2018; Jia et al. 2019; Wu & Xu 2019). Aiming at the advantages of the Elman neural network, it can also be applied to precipitation forecasting. The network structure is shown in Figure 1.

The Elman neural network learning algorithm adopts a momentum gradient descent back-propagation algorithm. The weights and thresholds are optimized by using the difference between the network output values and the output samples, to minimize the sum of squared errors of the network output layer.

Box et al. (1997) proposed the ARIMA model in the 1970s, which has been widely used in time-series analysis. The data can be processed smoothly and normally by studying the probability distribution of noise, thus solving the problem of sequence random disturbance.

The modeling ideas of ARIMA: the forecast object is regarded as a set of random series, which are approximately described by a certain mathematical model. Once the model is identified, the past and present values of the series can be used to forecast the future values. In the ARIMA (p, d, q) model, AR is the autoregressive component (Zhao & Chen 2015); I is difference; MA is the moving average component; p is the order of the autoregressive component; d is the differential times; q is the order of the moving average component.

The sequence of precipitation can be decomposed into low-frequency and high-frequency components by CEEMD decomposition, thus reducing the non-stationarity of the precipitation sequence. For high-frequency components, the non-stationarity of the sequence is relatively large. If the linear model is used, the forecasting effect will be insufficient. Therefore, the Elman neural network suitable for dealing with nonlinear problems was selected for forecasting. The low-frequency components show good regularity (periodicity and linearity). In the case of maximum decomposition scale, the low-frequency component is close to the original sequence in numerical value and shows good stationarity. Therefore, the ARIMA model suitable for linear problems can be selected for forecasting. The precipitation forecasting technical route is shown in Figure 2.

Zhengzhou is located in the south of the north China plain, the lower Yellow River, and central Henan province in the north. Zhengzhou is affected by the temperate continental monsoon climate instability and the variability of the weather system, resulting in great variation in annual precipitation. The data comes from the national meteorological information center of China. The location as shown in Figure 3.

As can be seen from Figure 4, the annual precipitation series in Zhengzhou shows an irregular fluctuation trend with rising and falling, and the series shows a greater randomness. Precipitation variation is affected by regional topography, land and sea position, underlying surface, monsoon climate, etc, so the sequence presents non-stationarity. However, in some years, precipitation may be less affected by external factors, and the sequence shows certain linear characteristics. Therefore, CEEMD suitable for nonlinear processing is selected as the decomposition tool, the Elman network suitable for nonlinear modeling and ARIMA suitable for linear modeling as the forecasting tool.

The scale change of precipitation is mainly reflected by the wavelet real-part contour map of precipitation, where the contour value represents the real part of the wavelet coefficient. The wavelet coefficient is positive, indicates that the precipitation is heavy; the wavelet coefficient is negative, means the precipitation is low. The value of wavelet variance reflects the proportion of a certain time-scale in the whole time-series. The larger the variance value, the stronger the oscillation, and the more significant the periodic change of the corresponding time-scale. The smaller the variance value, the less significant the periodic change of the corresponding time-scale. The wavelet series numerical modulus reflects the energy density corresponding to different time-scales. The larger the numerical modulus of the wavelet series, the stronger the periodicity of the time-scale, and vice versa.

Wavelet coefficients are calculated by MATLAB. After obtaining the wavelet coefficient, Surfer 12.0 software was used to make the contour map. The results are shown in Figure 5.

As can be seen from Figure 5, there are various time-scale structures in the precipitation series from 1951 to 2017, which are roughly 6a–20a, 30a–47a and 48a–64a. Among them, the corresponding wavelet coefficients of 48a–64a are large, the period is obvious, and the oscillation center is 55a. Precipitation changes on the central scale of 55a show the characteristics of less-more-less-more-less. The oscillation center is 40a at the scale of 30a–47a, and the variation characteristics of precipitation are more-less-more-less-more-less. The oscillation center is 15a on the scale of 6a–20a, and the precipitation changes have undergone several alternation cycles of less and more. In addition, from the perspective of the 55a scale, the contour line of the real part of the wavelet of precipitation in 2017 is not closed. It can be seen that the precipitation after 2017 will show a decreasing trend, but it is unknown how much the value is reduced. Therefore, the CEEMD − Elman + ARIMA model is used to forecast precipitation after 2017.

It can be seen from Figure 5 that precipitation in Zhengzhou has a certain characteristic of periodic change. However, it is puzzling that Figure 5 reflects the range of periodic changes, and the specific periodic value of the dominant position, which we cannot know. Therefore, the water variance is calculated, as shown in Figure 6.

As shown in Figure 6, there are three large peak points of annual precipitation in Zhengzhou, which indicates that there are three main cycles of annual precipitation. Among them, the three main cycles are 55 years, 15 years and eight years respectively, and the fluctuation of these three cycles controls the variation characteristics of precipitation in the whole time domain.

Figure 7 shows the evolution of annual precipitation in Zhengzhou. The wavelet series corresponding to 48a–64a has the largest numerical modulus value, indicating that the cycle change is the most obvious. The mode value of the wavelet coefficient corresponding to 6a–20a is followed by its periodic change. The wavelet coefficients of 21a–47a are the smallest and the corresponding periodic changes are the weakest.

The above periodic analysis is to forecast the future precipitation trend from the perspective of qualitative analysis. Therefore, from the perspective of quantitative analysis, the paper also predicts the precipitation value. That is, the CEEMD − Elman +ARIMA model is used to forecast precipitation in 2018–2022.

The sequence of precipitation is decomposed into sub-signals of different frequencies, namely IMF components and trend. The complex precipitation forecasting is equal to the sum of the forecasting values of different frequency subcomponents. By calculating the relative error of the sub-signal, the contribution rate of each sub-signal to the precipitation sequence can be analyzed, and the relative error of the sub-signal can be explained, whether it affects the relative error of precipitation.

After repeated testing, when the maximum decomposition is 5, the noise logarithm is 100 and the noise amplitude is 0.2. CEEMD has the best decomposition effect on precipitation, as shown in Figure 8.

As can be seen from Figure 8, the precipitation series is decomposed into four IMF components and one trend. Among them, the IMF₁ component has the lowest stationary, which is manifested as larger volatility, higher frequency and shorter wavelength. The amplitude and frequency of other the IMF components gradually decreases and the wavelength gradually increases. In addition, the change amplitude of the IMF component is significantly lower than the original series. Therefore, CEEMD decomposition greatly reduces the volatility and tendency of the sequence.

The sub-components obtained by CEEMD decomposition were divided into two parts for forecasting. For high-frequency components, the non-stationarity of the sequence is relatively large, so the Elman neural network suitable for dealing with nonlinear problems is selected for forecasting. For the low-frequency component (trend), it is close to the original sequence in value and has good stationarity. Therefore, the ARIMA model suitable for linear problems can be selected for prediction. In order to verify the forecasting accuracy of the model, the precipitation from 2013 to 2017 is forecast and error analysis is carried out. When the model is feasible, the future precipitation from 2018 to 2022 is forecast, so as to further study the evolution trend of precipitation.

The Elman neural network was used to forecast the four nonlinear IMF components. The trend has an obvious linear characteristic, and the ARIMA was used to forecast the trend. After a lot of calculation, the ARIMA model corresponding to the trend is ARIMA(1,1,13). The forecasting effect of each component is shown in Table 1.

As described in Table 1, the forecasting error of the IMF₁ component is relatively high, while the forecasting error of the other components is relatively small, and the forecasting effect of the trend is the best. IMF₄ and the trend are the most stable, while IMF₁ and IMF₂ are relatively weak, although the forecasting error of IMF₁ and IMF₂ is relatively high. Due to the small proportion of IMF₁ and IMF₂ in the precipitation series from 1951 to 2017, this does not affect the forecasting error of precipitation.

The forecasting value of precipitation is equal to the sum of the four IMF component forecasting and the trend forecasting. The forecasting results of precipitation in Zhengzhou are shown in Table 2.

Table 2 shows that the CEEMD − Elman + ARIMA model has a good forecasting effect. The mean square error is only 87.86 and the average relative error is only 14.1%. The relative error of the forecasting in 2013 is relatively high, 26.19%, but the overall forecasting effect is feasible.

When the CEEMD − Elman + ARIMA model was verified to be feasible, the precipitation in 2018–2022 was forecast, and the forecasting results are shown in Table 3.

As can be seen from Table 3, the precipitation in Zhengzhou showed a downward trend from 2018 to 2020. Precipitation has increased slightly since 2021, but annual precipitation forecasts for 2018–2022 are still generally lower than in 2017. Due to the lack of real precipitation values after 2017 as a reference, the forecast values of 2018–2020 are highly reliable, while the forecast value error of 2021 and 2022 may be large. Therefore, it is reasonable to conclude from the wavelet cycle periodic analysis that the precipitation will decline in the next few years starting from 2017.

In the introduction, a variety of forecasting methods are mentioned. Among them, Empirical Mode Decomposition (EMD) and the Elman combined model (EMD–Elman), Markov and ARIMA combined model (Markov–ARIMA), Back Propagation neural network (BPNN) and Weather Research and Forecasting (WRF) and Cloud theory (WRF-Cloud) are representative. To further assess the capabilities of the proposed model, the simulation results of several models in precipitation forecasting are compared. Due to the limited grasp of various forecasting methods, only the simulation results of CEEM − Elman + ARIMA, EMD–Elman, Elman and Markov–ARIMA are compared. If there is time, the next step will be a systematic and comprehensive comparison of the advantages of various forecasting methods. The forecasting errors of the four models were made, as shown in Table 4.

As can be seen from Table 4, the relative error and average error of CEEMD − Elman + ARIMA are lower than that of Markov–ARIMA, EMD–Elman and Elman, indicating that CEEMD − Elman + ARIMA is more advantageous in precipitation forecasting. Compared with the other forecasting methods, the decomposition–reconstruction forecasting model can reduce the forecasting error. The reasons are as follows: the original signal is decomposed into different frequency components by the decomposition tools (CEEMD, EEMD, EMD), and the forecasting value of precipitation is equal to the sum of the forecasting values of the high-frequency component and low-frequency component. That is, the overall forecasting error of the model is determined by the forecasting error of these high- and low-frequency components. In this way, even if the forecasting effect of some components is not good, the overall forecasting effect will not be affected.

CEEMD is proposed on the basis of EMD, and the model constructed adopts different forecasting methods for high- and low-frequency components. Therefore, the forecasting effect of CEEMD − Elman + ARIMA is more advantageous than the single decomposition–reconstruction forecasting model, such as EMD–Elman. Compared with the forecasting model without decomposition, the forecasting effect of Elman is determined by the original data itself. If the forecasting effect of Elman is not good, it will directly affect the whole error. Markov–ARIMA averages the forecasting results of the ARIMA and Markov model, which reduces the forecasting error to some extent, but the forecasting accuracy is lower than that of the decomposition–reconstruction forecasting model. In conclusion, CEEMD − Elman + ARIMA is more suitable for precipitation forecasting, and has more advantages and stability.

Wavelet analysis was applied to the annual precipitation cycle analysis of Zhengzhou from 1951 to 2017, and it was concluded that the annual precipitation of Zhengzhou was mainly affected by scale fluctuations of 55a, 40a and 15a. Furthermore, the forecasting of precipitation in the future years of 2018–2022 was carried out, and from the qualitative and quantitative perspective, the precipitation of Zhengzhou will show a downward trend from 2017.

The non-stationarity of precipitation time-series is reduced by CEEMD decomposition, and Elman and ARIMA are used to forecast the data after decomposition. The CEEMD − Elman + ARIMA forecasting model was constructed and applied to Zhengzhou precipitation forecasting, and the average relative error of the forecasting in 2013–2017 was 14.1%.

This is a new attempt to apply the CEEMD − Elman + ARIMA coupling model to the forecasting of precipitation in Zhengzhou. Its extension and the improvement of model accuracy still need further study. The model does not consider the physical mechanism of precipitation evolution and the long-term forecasting. Furthermore, how to deal with the subcomponent properly and improve the forecasting accuracy of the model are the next research direction and emphasis.

This work is financially supported by Collaborative Innovation Center of Water Resources Efficient Utilization and Protection Engineering, Henan Province, Water Environment Governance and Ecological Restoration Academician Workstation of Henan Province, Program for Science & Technology Innovation Talents in Universities of Henan Province (No. 15HASTIT049). Our gratitude is also extended to reviewers for their efforts in reviewing the manuscript and their very encouraging, insightful and constructive comments.