## Abstract

Due to the nonlinear characteristics of runoff data and the poor performance of the single prediction model, a weighted integrated modified complementary ensemble empirical mode decomposition (MCEEMD)-based model was proposed to predict the monthly runoff of three hydrological stations in the lower reaches of the Yellow River. In this model, particle swarm optimization (PSO) was used to optimize the parameters of support vector regression (SVR), back propagation neural network (BP), long short-term memory neural network (LSTM) that constitute it. The weight coefficients and frequency terms decomposed by MCEEMD were used to obtain the final prediction results. Results indicated that this model performs better than other models, with the Nash–Sutcliffe efficiency (NSE) reaching above 0.92, qualification rate (QR) reaching above 75% and all error indicators being minimal. In addition, considering the influence of extreme weather and climate anomalies, the integrated model combined the atmospheric circulation anomalies factors and the results can still be improved. It can be verified that this weighted integrated model can be used for the stable and accurate predication of medium- and long-term runoff.

## HIGHLIGHTS

Decomposition of runoff into smooth sequences using MCEEMD reduces the accumulation of errors associated with the CEEMD.

A weighted integrated method was used to develop predictive models based on the MCEEMD decomposition.

In view of the influence of extreme weather and abnormal climate on runoff prediction accuracy, the atmospheric circulation anomaly factors were used as the input data of the model to predict.

### Graphical Abstract

## INTRODUCTION

Runoff is an important condition of regional industrial and agricultural water supply, and also a restricting factor of regional socio-economic development. The measurement, calculation, and forecast of runoff are important tasks of water conservancy construction. For example, they are of great significance for the operation, planning, and dispatching of hydropower stations to predict monthly runoff accurately (Huang *et al.* 2015a, 2015b). Runoff prediction has received a lot of attention in recent years.

The traditional model usually uses mathematical statistics. Typically, Tesfaye *et al.* (2006) used a periodic autoregressive moving average model to produce a realistic simulation for monthly runoff data of the Fraser River in British Columbia. Ouyang *et al.* (2010) used the dynamic time warping distance method to perform a similarity search of floods at Shaliguilanke station in the Tarim River basin. With the development of data science, there are more and more studies on runoff prediction based on machine learning and neural networks. Typically, Cui (2013a, 2013b) proposed to improve the Elman neural network prediction model and the hidden layer back propagation neural network (BP) prediction model. Lu and Zhou (2014) screened input factors of a forecast model by using the mutual information method. In the BP neural network model, mean square error and mutual information are used as objective functions to measure the correlation between factors, optimize the final prediction factors, and apply the prediction factors to the runoff prediction of Biliu River in flood season; the model can identify multiple complex correlations between forecast factors and forecast quantity. Xiang *et al.* (2020) used long short-term memory neural network (LSTM) and the seq2seq structure to estimate hourly rainfall-runoff.

However, the generation process of runoff tends to be uncertain, highly nonlinear, and time varying, especially when extreme weather appears, thus the monthly streamflow series contains different frequency components (Huang *et al.* 2015a, 2015b). Similarly, the nonlinear character of runoff contradicts the data requirements of many models. Moreover, the single prediction model (i.e., a model using only one algorithm) has limitations, such as the over-fitting phenomenon and local optimality in neural network, and there is no uniform standard in kernel function selection and parameter calibration in the support vector machine method (Han *et al.* 2017).

In order to solve the above defects, this study has developed a weighted integrated model based on the modified complementary ensemble empirical mode decomposition (MCEEMD) to predict monthly runoff. In this study, the particle swarm optimization optimize support vector regression (PSO-SVR), particle swarm optimization optimize back propagation neural network (PSO-BP), and particle swarm optimization optimize long short-term memory neural network (PSO-LSTM) methods with excellent runoff prediction performance were selected as the compositions of the integrated model, and the runoff data decomposed by MCEEMD were used as the training data. Then, the weights of each model's frequency terms were calculated according to the errors of the single algorithms, and the runoff prediction results of the weighted integrated model were finally obtained by adding each single model's weighted frequency terms. In addition, the optimal factors of atmospheric circulation anomaly factors were selected as the additional item of training data of the integrated model to predict runoff, which improved the prediction accuracy and stability of the model in extreme weather periods. Compared with other models, the weighted integrated model had a better prediction performance, especially during periods of extreme weather and weather anomalies, which can provide a reference for regional water resources allocation and regional water resources optimization. It also provided a new idea for the development of hydrological prediction.

## METHODS

### Modified complementary ensemble empirical mode decomposition (MCEEMD)

Huang *et al.* (1998) of the National Aeronautics and Space Administration proposed a method of data decomposition called empirical mode decomposition (EMD), which essentially identifies all vibration modes contained in a signal through characteristic time scales. The basic principle of EMD is to decompose the complex sequence of input into a finite number of intrinsic mode functions (IMF) from high to low frequencies and a residual term. The decomposed IMF contains local characteristic signals of the original signal at different time scales, and reduces the mutual interference between different trend information. IMF must meet with two characteristics:

- 1.
the number of extremum points and zero crossings must be the same or at most one different throughout the data segment;

- 2.
in any case, the mean value of the envelope defined by the local maxima and that of the local minima is zero.

However, EMD will lead to the frequent occurrence of mode aliasing in the decomposition process, which will undermine the physical significance of IMF (Han *et al.* 2017). Huang & Wu (2008) proposed the ensemble empirical mode decomposition (EEMD) to solve the defects of EMD. Although EEMD effectively solves the modal aliasing phenomenon of EMD by adding white noise, it is affected by residual noise in signal reconstruction. Yeh *et al.* (2010) improved the EEMD method and proposed the complementary ensemble empirical mode decomposition (CEEMD), which not only solved the problem of mode aliasing that EMD is prone to, but also avoided the influence of residual noise of EEMD in signal reconstruction. In recent years, CEEMD has been applied to many models for runoff prediction (Zhang *et al.* 2019).

However, the components obtained after the decomposition of CEEMD are too much, especially the nonlinear and non-stationary features of IMF1, which will bring a large number of calculations and error accumulation to the prediction model. MCEEMD reconstructed the IMFs decomposed by CEEMD through the fluctuation frequency and amplitude to obtain the high frequency (HF) term, medium frequency (MF) term, low frequency (LF) term, and residual term (Res) to address the above deficiencies. The process of decomposition for MCEEMD is as follows:

- 1.Adding a random sequence of positive and negative white noise (i.e., the original signal plus the white noise plus the original signal minus the white noise) with the same amplitude and phase angle difference of to the original signal . Adding a new noise with the same amplitude every time to get a new signal.
- 2.
- 3.
Repeating step 1 and step 2, adding different positive and negative white noise sequences each time to get M (i.e., number of times white noise was added) group of IMF components and residual terms.

- 4.
- 5.
- 6.

### Support vector regression (SVR)

*et al.*2015a, 2015b; Chu

*et al.*2016). The calculation formula of SVR is as follows:where

*l*is the number of the data set; and are lagrange multiplier; is the kernel function;

*b*is the deviation between the true value and the predicted value.

*w*is the weight vector; and are relaxation factor;

*C*is the penalty coefficient, an appropriate

*C*can make the model have a better generalization ability; is the mapping from the input space to the feature space; is the constant deviation;

*n*is the sample size. The optimized process is subjected to the constraints under it.

### Back propagation neural network (BP)

BP is a multi-layer feed-forward network trained by error back propagation. BP uses the gradient search technology to minimize the mean square error between the actual output value and the expected output value of the network. BP is usually combined with other algorithms for runoff prediction (Lu & Zhou 2014). The structure of BP is shown in Figure 1.

The input signal X through the hidden layer node affects the output node, each nerve input includes input vector and the desired output vector. If the deviation vector is not as expected, the weight W and threshold value are adjusted in the opposite direction to reduce the error along the gradient direction. After repeated transmission, reduce the error to expectations and get the optimal solution Y.

### Long short-term memory network (LSTM)

LSTM is an improved recurrent neural network (RNN). LSTM uses gates to control the memory process, making up for the loss caused by RNN gradient explosion and gradient disappearance to a large extent, and solving the problem that RNN cannot handle long-distance dependence. LSTM adds input gate, forgetting gate, output gate, and internal memory unit on the basis of RNN (Sepp & Jürgen 1997). The unit structure diagram of the LSTM model is shown in Figure 2.

*W*and

*U*are weight parameters;

*b*is the bias; is the cell state. The input gate is obtained through the activation function after the transformation between input layer and output of the previous hidden layer; the values of weight parameters

*W*and

*U*, bias

*b*depend on the training results of the model; the result of the input gate is a vector, which is responsible for calculating the current input state based on the last output and this input; the calculation process of forgetting gate and output gate is similar to that of input gate . Such a control process enables the LSTM model to quickly and accurately learn the long-term dependence between sequences, which has great advantages in the processing of time series runoff data.

### Particle swarm optimization (PSO)

*i*in the

*j*dimension at the

*t*iteration; is the position of

*i*; is the weight of inertia; and are learning factors; is the individual extremum point at the

*t*iteration of the particle swarm; is the global extremum; and are random numbers between 0 and 1.

The optimization process of PSO is as follows:

- 1.
Initialize the parameters of the particle swarm, including population size, number of iterations, learning factor, and range of speed and location.

- 2.
Choose the fitness function of PSO.

- 3.
Create a particle randomly, including the various parameters of the algorithm.

- 4.
Get the local and global optimal positions of the particle according to the fitness function.

- 5.
Update the particle speed and position according to Equations (14) and (15). Update the local and global optimal values according to the new fitness function.

- 6.
After reaching the maximum number of iterations, the optimal particle is taken as the parameter of the algorithm.

### Weighted integrated model

*n*is the number of the single model; , , , and are the errors between the predicted and true values of the frequency terms for different single models. The average relative error of the predicted values was used as the error in this study.

In order to improve the accuracy of runoff prediction, PSO optimize SVR based on MCEEMD (MCEEMD-PSO-SVR), PSO optimize BP based on MCEEMD (MCEEMD-PSO-BP), and PSO optimize LSTM based on MCEEMD (MCEEMD-PSO-LSTM) were used to develop a weighted integrated model in this study. The weighted integrated model is a process of ‘Decomposition – Reconstruction – Optimization – Prediction – Weighted Integration – Reconstruction’. Figure 3 shows the framework of the model, including:

- 1.
Decompose the original runoff data

*R*into three frequency terms (, , ) and one residual item by using MCEEMD. - 2.
Use each frequency term as the training data of each single model.

- 3.
Optimize the different single models by using PSO.

- 4.
Calculate the weights (, , , ) of each frequency term from different single models by using the average relative error.

- 5.
Develop the weighted integrated model based on the weights and reconstruct the predicted results of the frequency terms to obtain the final prediction result .

### Model evaluation

*et al.*2016). The calculation formulas are as follows:where is the predicted value; is the measured value; is the mean value of the measured value;

*n*is the total number of predicted samples; and

*m*is the total number of qualified predicted samples.

## STUDY AREA AND DATA

Lijin hydrological station is located in Dongying, Shandong province, 104 km from the mouth of the Yellow River, which is located in the lowest reaches of the Yellow River. Gaocun hydrological station is an important control station for the Yellow River flowing into Shandong province, with a section 579.1 km from the estuary and a catchment area of 734,146 km^{2}. Aishan hydrological station is located in Liaocheng, Shandong province. The above three hydrological stations are responsible for providing water conditions for flood control and water resources dispatching in the lower reaches of the Yellow River, studying and exploring the changing rules of hydrological factors, and collecting hydrological data for river regulation and water and sand resource utilization in the lower reaches of the Yellow River. The data from these stations are complete and consistent with the law of hydrology. Figure 4 shows the location of these hydrological stations.

This study used monthly runoff data for a total of 780 months from January 1950 to December 2014. The training set used runoff data from January 1950 to December 2001, a total of 624 months. The validation set runoff data from January 2002 to December 2014, a total of 156 months, was used to verify results of the model. Table 1 shows the characteristics of monthly runoff data of the three hydrological stations. Coefficient of variation (CV) is a relative index to measure the dispersion degree of runoff data, skewness represents the degree of asymmetry in the distribution of runoff data, and kurtosis represents the peak shape characteristic of probability density distribution curve. Table 1 shows the following:

- 1.
The CVs of each hydrological station in the lower reaches of the Yellow River are 0.86–1.07, which indicates that the change of monthly runoff is dramatic and the distribution is uneven in the year.

- 2.
All skewness are greater than 0, and all are positive deviations, indicating that a small number of monthly runoff data are large.

- 3.
All kurtosis are greater than the kurtosis of normal distribution and uniform distribution, indicating that the monthly runoff data differ greatly from the mean value and have many extreme values.

Station . | Average monthly runoff . | CV . | Skewness . | Kurtosis . |
---|---|---|---|---|

Lijin | 25.65 (10^{8}m^{3}) | 1.07 | 2.05 | 5.11 |

Gaocun | 29.01 (10^{8}m^{3}) | 0.86 | 2.10 | 5.02 |

Aishan | 27.93 (10^{8}m^{3}) | 0.96 | 2.14 | 5.51 |

Station . | Average monthly runoff . | CV . | Skewness . | Kurtosis . |
---|---|---|---|---|

Lijin | 25.65 (10^{8}m^{3}) | 1.07 | 2.05 | 5.11 |

Gaocun | 29.01 (10^{8}m^{3}) | 0.86 | 2.10 | 5.02 |

Aishan | 27.93 (10^{8}m^{3}) | 0.96 | 2.14 | 5.51 |

The variation process of monthly runoff is shown in Figure 5. It can be seen from Table 1 and Figure 5 that the monthly runoff had a large amplitude of variation and a weak periodicity, presenting a nonlinear and non-stationary characteristic state.

## RESULTS

### The original series decomposed by MCEEMD

MCEEMD was used to decompose the monthly runoff data from January 1950 to December 2014. The decomposed results included three frequency terms and one residual term. The components were reduced to four after the reconstruction, and the characteristics and objective laws of the original time series data were still retained. In addition, the high frequency term obtained by combining IMF1 and IMF2 eliminated the nonlinear characteristics of IMF1. The monthly runoff data after decomposition by MCEEMD are shown in Figure 6. As can be seen from Figure 6, although the reconstructed high frequency term retains some features of IMF1 and IMF2, it has reduced the nonlinear feature, and the features of the intermediate frequency term and low frequency term are similar to those of the IMF. It can be seen that HF has the highest frequency, the largest fluctuation, and the shortest wavelength. The frequency gradually decreases in MF and LF, the fluctuations gradually weaken, the wavelength gradually increases and the periodicity becomes stronger. The residual term represents the trend of monthly runoff at stations over time from 1950 to 2014. As can be seen in Figure 6, the monthly runoff shows a declining trend year to year, with a low decline rate during the early months and a higher decline rate in the later months.

### Data processing

After prediction, the data can be re-scaled following the contrary procedure of Equation (25).

### Parameter selection by PSO

Proper parameter selection can directly determine the performance of the training process. The parameters of the predictive algorithm are usually selected manually. However, manual selection is too subjective. In addition, the numerous parameters of the algorithm can make a lot of unnecessary work. Therefore, algorithms in this study were optimized using PSO with global optimization capabilities.

After optimization, the radial basis function (RBF) kernel function of SVR was selected as the kernel function, the optimal penalty factor *C* of SVR was 8.3598, and the parameter *g* of the optimal kernel function was 0.031413 for SVR. The input layer of BP was set to 7 and the output layer to 1. The activation function of LSTM was selected as rectified linear unit (ReLU), the units of the LSTM were selected as 20, the loss function of LSTM was selected as MAE, and the optimizer of LSTM was selected as adaptive moment estimation (Adam).

The activation function ReLU of LSTM is a non-saturated activation function. Compared with other activation functions, ReLU has the advantages of fast convergence and fast calculation speed (Ramachandran *et al.* 2017). The optimizer Adam, which is computationally efficient and has little memory requirements, combines first moment estimation (the mean of the gradients) and second moment estimation (the uncentralized variance of the gradients) (Kingma & Ba 2014).

### Predicted results

According to Equation (17), the weights of each frequency term for different single models in the weighted integrated model were obtained by calculating average relative errors of the predicted results for the same frequency term of different models, and the results are shown in Table 2. The greater the weight, the greater the contribution of the single model to the integrated model.

Stations . | Terms . | MCEEMD-PSO-SVR . | MCEEMD-PSO-BP . | MCEEMD-PSO-LSTM . |
---|---|---|---|---|

Lijin | HF | 0.32 | 0.26 | 0.43 |

MF | 0.36 | 0.30 | 0.33 | |

LF | 0.31 | 0.32 | 0.37 | |

Res | 0.29 | 0.34 | 0.37 | |

Gaocun | HF | 0.36 | 0.24 | 0.40 |

MF | 0.40 | 0.22 | 0.38 | |

LF | 0.34 | 0.28 | 0.37 | |

Res | 0.30 | 0.35 | 0.35 | |

Aishan | HF | 0.23 | 0.22 | 0.55 |

MF | 0.35 | 0.29 | 0.36 | |

LF | 0.29 | 0.35 | 0.36 | |

Res | 0.33 | 0.36 | 0.31 |

Stations . | Terms . | MCEEMD-PSO-SVR . | MCEEMD-PSO-BP . | MCEEMD-PSO-LSTM . |
---|---|---|---|---|

Lijin | HF | 0.32 | 0.26 | 0.43 |

MF | 0.36 | 0.30 | 0.33 | |

LF | 0.31 | 0.32 | 0.37 | |

Res | 0.29 | 0.34 | 0.37 | |

Gaocun | HF | 0.36 | 0.24 | 0.40 |

MF | 0.40 | 0.22 | 0.38 | |

LF | 0.34 | 0.28 | 0.37 | |

Res | 0.30 | 0.35 | 0.35 | |

Aishan | HF | 0.23 | 0.22 | 0.55 |

MF | 0.35 | 0.29 | 0.36 | |

LF | 0.29 | 0.35 | 0.36 | |

Res | 0.33 | 0.36 | 0.31 |

The predicted values of the weighted integrated model were calculated according to the weights obtained in Table 2. Table 3 indicates the prediction performance of the weighted integrated model on each station in training set and verification set. It can be seen that the NSEs of each station during the training period and the verification period were all greater than 0.5, indicating that the predicted results were reliable. NSE was greater than 0.98 and QR was greater than 85% in each station during the training period. NSE was greater than 0.92 and QR was greater than 75% in each station during the verification period. The predicted results accord with the standard of hydrological forecast.

Indicators . | Training set . | Verification set . | ||
---|---|---|---|---|

Stations . | NSE . | QR/% . | NSE . | QR/% . |

Lijin | 0.98 | 87 | 0.93 | 78 |

Gaocun | 0.99 | 86 | 0.95 | 79 |

Aishan | 0.98 | 85 | 0.92 | 75 |

Indicators . | Training set . | Verification set . | ||
---|---|---|---|---|

Stations . | NSE . | QR/% . | NSE . | QR/% . |

Lijin | 0.98 | 87 | 0.93 | 78 |

Gaocun | 0.99 | 86 | 0.95 | 79 |

Aishan | 0.98 | 85 | 0.92 | 75 |

Since the data in this study were divided into training set and verification set, based on the following considerations, we did not pay much attention to the prediction results of the training set:

- 1.
when establishing the prediction model, the performance of the verification set can represent the real application effect;

- 2.
for the adaptive prediction model, we do not care about the prediction performance of the training set;

- 3.
the output of the training set is not a real prediction, it is only an input to the model of the verification set.

Therefore, in the following sections, we just show the prediction results of the verification set. The errors of the weighted integrated model were compared with those of the MCEEMD-PSO-SVR model, the MCEEMD-PSO-BP model, and the MCEEMD-PSO-LSTM model, and the results are shown in Table 4. It can be obviously seen from Table 4 that the weighted integrated model has the best performance and it outperforms the other models in terms of all the error coefficients, acquiring the best MAE, MAPE, and RMSE. The MAE, MAPE, and RMSE of the integrated model in Lijin, Gaocun, and Aishan hydrological stations were reduced by 24.36%–30.04%, 18.54%–24.82%, 8.63%–16.37%; 29.31%–34.66%, 10.91%–35.21%, 15.73%–26.52%; and 16.17%–24.23%, 12.44%–13.55%, 13.08%–16.22%, respectively, compared with the other three models. In addition, the values of MaxAE and MinAE for the integrated model and the difference between them were the smallest of all models.

Stations . | Indicators . | MCEEMD-PSO-SVR . | MCEEMD-PSO-BP . | MCEEMD-PSO-LSTM . | Integrated model . |
---|---|---|---|---|---|

MaxAE (10^{8}m^{3}) | 12.43 | 13.06 | 18.17 | 10.71 | |

MinAE (10^{8}m^{3}) | 0.05 | 0.12 | 0.09 | 0.01 | |

Lijin | MAE (10^{8}m^{3}) | 2.41 | 2.53 | 2.34 | 1.77 |

MAPE (%) | 30.86 | 29.43 | 28.48 | 23.20 | |

RMSE (10^{8}m^{3}) | 3.37 | 3.42 | 3.13 | 2.86 | |

MaxAE (10^{8}m^{3}) | 13.93 | 17.89 | 9.28 | 7.82 | |

MinAE (10^{8}m^{3}) | 0.08 | 0.09 | 0.07 | 0.01 | |

Gaocun | MAE (10^{8}m^{3}) | 2.51 | 2.42 | 2.32 | 1.64 |

MAPE (%) | 15.58 | 18.66 | 13.57 | 12.09 | |

RMSE (10^{8}m^{3}) | 3.21 | 3.28 | 2.86 | 2.41 | |

MaxAE (10^{8}m^{3}) | 14.98 | 23.00 | 11.05 | 8.44 | |

MinAE (10^{8}m^{3}) | 0.03 | 0.05 | 0.02 | 0.01 | |

Aishan | MAE (10^{8}m^{3}) | 2.35 | 2.49 | 2.60 | 1.97 |

MAPE (%) | 16.42 | 16.32 | 16.53 | 14.29 | |

RMSE (10^{8}m^{3}) | 3.21 | 3.33 | 3.22 | 2.79 |

Stations . | Indicators . | MCEEMD-PSO-SVR . | MCEEMD-PSO-BP . | MCEEMD-PSO-LSTM . | Integrated model . |
---|---|---|---|---|---|

MaxAE (10^{8}m^{3}) | 12.43 | 13.06 | 18.17 | 10.71 | |

MinAE (10^{8}m^{3}) | 0.05 | 0.12 | 0.09 | 0.01 | |

Lijin | MAE (10^{8}m^{3}) | 2.41 | 2.53 | 2.34 | 1.77 |

MAPE (%) | 30.86 | 29.43 | 28.48 | 23.20 | |

RMSE (10^{8}m^{3}) | 3.37 | 3.42 | 3.13 | 2.86 | |

MaxAE (10^{8}m^{3}) | 13.93 | 17.89 | 9.28 | 7.82 | |

MinAE (10^{8}m^{3}) | 0.08 | 0.09 | 0.07 | 0.01 | |

Gaocun | MAE (10^{8}m^{3}) | 2.51 | 2.42 | 2.32 | 1.64 |

MAPE (%) | 15.58 | 18.66 | 13.57 | 12.09 | |

RMSE (10^{8}m^{3}) | 3.21 | 3.28 | 2.86 | 2.41 | |

MaxAE (10^{8}m^{3}) | 14.98 | 23.00 | 11.05 | 8.44 | |

MinAE (10^{8}m^{3}) | 0.03 | 0.05 | 0.02 | 0.01 | |

Aishan | MAE (10^{8}m^{3}) | 2.35 | 2.49 | 2.60 | 1.97 |

MAPE (%) | 16.42 | 16.32 | 16.53 | 14.29 | |

RMSE (10^{8}m^{3}) | 3.21 | 3.33 | 3.22 | 2.79 |

Figure 7 illustrates the runoff prediction of Lijin, Gaocun, and Aishan stations by the MCEEMD-PSO-SVR, MCEEMD-PSO-BP, MCEEMD-PSO-LSTM and the weighted integrated model in the verification set. Figure 8 shows the scatter plots of predicted value and observed value of the three stations using all models, respectively.

Some studies have found that atmospheric circulation anomaly factors have a strong correlation with extreme weather and climate anomalies (Wang 2009; Huang *et al.* 2015a, 2015b, 2017; Liu *et al.* 2017; Yang *et al.* 2016). Unfortunately, the present study mainly focused on the atmospheric circulation anomaly factors to the overall effect of the runoff, without considering the decomposition technique abnormal atmospheric circulation factors for size (Meng *et al.* 2017). Therefore, monthly scale data of Nino3.4, Arctic Oscillation (AO), Pacific Decadal Oscillation (PDO) and Atlantic Multidecadal Oscillation (AMO) were used as additions to each frequency term of the input data of the weighted integrated model to improve the accuracy and stability of runoff forecasting during some periods. In order to select optimal factors, the correlation coefficient method was used to analyze the runoff frequency terms decomposed by MCEEMD and the atmospheric circulation anomaly factors with different time delays. Results are shown in Figure 9, where the darker the shade represents the higher the correlation. The factors whose correlation coefficient is greater than 0.3 were selected as the additional term of input data. It can be seen from Figure 9(a) that in Lijin station, Nino3.4 has significant influence on the high frequency term and the residual term, while AMO has significant influence on the residual term. It can be seen from Figure 9(b) that in Gaocun station, Nino3.4 has significant influence on the high frequency term and medium frequency term, while PDO has significant influence on the medium frequency term and lower frequency term. It can be seen from Figure 9(c) that in Aishan station, Nino3.4 has significant influence on the high frequency term, lower frequency term, and residual term.

Atmospheric circulation anomaly factors that have significant influence on each frequency term of runoff in the above analysis were selected and used as additional eigenvectors to the weighted integrated model input data. In the weighted integrated model, the frequency terms and the atmospheric circulation anomaly factors were used as its input data, and the model first obtained the prediction of each frequency term and then reconstructed each frequency term according to the weights to obtain the final runoff prediction.

It can be seen from Table 5 that the weighted integrated model combined with the information of atmospheric circulation anomaly factors improved the runoff prediction accuracy of each station in different degrees. The NSE of each station increased by 1.05%–3.26%, the QR increased by −2.67% to 8.97%, the RMSE decreased by 2.10%–3.94%, and the difference between MaxAE and MinAE decreased by 8.78%–24.11%, among which the QR of Aishan station decreased by 2.67%, which may be due to the weak influence of the atmospheric circulation anomaly factor on Aishan station.

Stations . | NSE . | QR (%) . | RMSE (10^{8}m^{3})
. | Difference (10^{8}m^{3})
. |
---|---|---|---|---|

Lijin | 0.93 → 0.95 | 78 → 85 | 2.86 → 2.80 | 10.70 → 8.12 |

Gaocun | 0.95 → 0.96 | 79 → 82 | 2.41 → 2.34 | 7.81 → 6.26 |

Aishan | 0.92 → 0.95 | 75 → 73 | 2.79 → 2.68 | 8.43 → 7.69 |

Stations . | NSE . | QR (%) . | RMSE (10^{8}m^{3})
. | Difference (10^{8}m^{3})
. |
---|---|---|---|---|

Lijin | 0.93 → 0.95 | 78 → 85 | 2.86 → 2.80 | 10.70 → 8.12 |

Gaocun | 0.95 → 0.96 | 79 → 82 | 2.41 → 2.34 | 7.81 → 6.26 |

Aishan | 0.92 → 0.95 | 75 → 73 | 2.79 → 2.68 | 8.43 → 7.69 |

## DISCUSSION

According to the residual term decomposed by MCEEMD from the monthly runoff data of Lijin, Gaocun, and Aishan hydrological stations, the overall runoff monitored by stations showed a declining trend, especially in recent years. It indicates that the lower reaches of the Yellow River water shortage and other problems are increasingly serious, one important reason it is necessary to improve the accuracy of runoff prediction to promote the optimal allocation of water resources and maximize the benefit of water resources.

It appears from Table 4 that the MCEEMD-PSO-LSTM model is second only to the weighted integrated model based on MCEEMD in performance. In addition, it can be found from Table 2 that the MCEEMD-PSO-LSTM model has the largest weight in the weighted integrated model, thus it can be explained that the weight determines the importance of a single model in the weighted integrated model and also indicates the performance of a single model. It can be seen from Table 4 that, although the models’ errors are different, the differences among MCEEMD-PSO-SVR, MCEEMD-PSO-BP, and MCEEMD-PSO-LSTM are not too large, which indicates that the three models are similar in the performance of runoff prediction. The integrated model requires that the models that comprise it be close in performance, and the discussion above reveals that the three single models that comprise the integrated model meet this requirement. Moreover, the use of the weight coefficient further reduces the difference between single models, so that the single model with good performance takes up a larger proportion. In this way, the advantages of the model are magnified and the disadvantages of the model are reduced, which is conducive to improving the accuracy of runoff prediction.

It can be seen in Figure 7(a) that the cycle and the trend of the prediction results of the weighted integrated model were completely consistent with the original monthly runoff data in Lijin hydrological station. Moreover, the weighted integrated model had the best predictive performance among all models. Similarly, it can be observed from Figure 7(b) and 7(c) that the weighted integrated model had the better performance compared with the other models. In flood season, the weighted integrated model had higher prediction accuracy than other models.

Figure 8 indicates that the weighted integrated model had the best performance for forecasting the monthly runoff, as the linear trend line of it was closest to the that goes through the point (0,0) compared with the MCEEMD-PSO-SVR, the MCEEMD-PSO-BP, and the MCEEMD-PSO-LSTM. Similarly, the above conclusions can also be observed from Figure 8(b) and 8(c). represents the explanatory power of the equation variable *x* to *y*, and the closer to 1, the better the model fits the data. It can be seen that the of the weighted integrated model was the largest and closest to 1, which proves that the prediction performance of it was the best of all models.

It can be observed in scatter plots that the peak values of the weighted integrated model were mostly above the line of , indicating that the forecasted peaks are lower than the observed ones. This phenomenon may result from the limitation of our model. However, the values, except the peak value, in our model were evenly distributed around the line compared with other models, indicating that its fitting result had a better stability. In addition, the minimum difference between MaxAE and MinAE also indicated that.

From Figure 9 it appears that the influences of abnormal atmospheric circulation factors on the runoff of each station are different. This is because in this paper, the locations of the three hydrological stations are from west to east, and the influence of different locations by the Asian summer monsoon are different, so the correlation between the runoff of hydrological stations and abnormal atmospheric circulation factors is different.

Atmospheric circulation anomaly factors are greatly affected by extreme weather and climate anomalies, and the violent change of runoff in flood season is related to extreme weather. Therefore, the integrated model has a stronger ability to simulate the flood season runoff and its extreme value after adding atmospheric circulation anomaly factors, which improves the stability of the model and the prediction accuracy.

For complicated and non-stationary monthly runoff, the single prediction model has low accuracy and poor prediction effect and is not suitable for medium- and long-term prediction. CEEMD reduces the mutual interference between different trend information, retains the objective law of the original data, and provides a suitable data basis for the prediction algorithm. However, if modeling and forecasting are carried out for each of the eight components obtained after the decomposition of CEEMD, a large amount of work and errors will be incurred. Therefore, this study proposed MCEEMD, which reconstructed the IMFs with similar fluctuation frequency and amplitude to obtain the frequency terms, which is then used as the input of the prediction algorithm. The weighted integrated model combines the advantages of each algorithm and determines the importance of each algorithm in the integration model according to the weight coefficient. By combining the MCEEMD with the weighted integrated model, runoff can be reduced into a stationary sequence to increase the performance of prediction model, and the defects of the single prediction algorithm were also reduced. At the same time, the addition of abnormal atmospheric circulation factors improves the stability and prediction accuracy of the model under extreme weather conditions.

## CONCLUSIONS

In order to improve the accuracy of runoff prediction, a weighted integrated model based on MCEEMD was proposed, which was verified by using the monthly runoff data from three hydrological stations (Lijin, Gaocun, Aishan). The MCEEMD reduces the calculation problem caused by excessive components by establishing frequency terms, and also solves the problem that prediction algorithm has poor prediction effect on IMF1. The weighted integrated model uses MCEEMD to process training data, combines PSO-SVR, PSO-BP, and PSO-LSTM, and integrates the advantages of each single model, thus improving the prediction accuracy, stability, and fitting effect of the single model. Compared with single prediction models, the weighted integrated model has a better performance in runoff prediction and higher accuracy of prediction results. The prediction results of the weighted integrated model are in line with the hydrological prediction standard and have a high reliability. Compared with the single prediction algorithms, the MAE, MAPE, and RMSE of the weighted integrated model were reduced by more than 15.81%, 10.91%, and 13.42%, respectively. Therefore, the combination of the MCEEMD method and weighted integrated model is viable in the medium- and long-term runoff prediction. In addition, it can also provide guidance for flood control and drought relief.

Considering that runoff is affected by climate and extreme weather, atmospheric circulation anomaly factors with strong correlation with runoff were added and selected as additional items of input data of the integrated model in this study. The results show that the model combined with atmospheric circulation anomaly factors has higher prediction accuracy and stability, especially in extreme weather and flood seasons.

## ACKNOWLEDGEMENTS

Research of this article is granted by NSFC under grant no. U1604152 and the key research project of Ministry of Water Resources of the People's Republic of China. The authors of this paper would like to express their thanks to CMDC for providing the data. The authors declare no conflict of interest. (http://61.163.88.227:8006/hwsq.aspx, http://data.cma.cn/, https://www.ncdc.noaa.gov/).

## DATA AVAILABILITY STATEMENT

All relevant data are available from an online repository or repositories. http://61.163.88.227:8006/hwsq.aspx, http://data.cma.cn/, https://www.ncdc.noaa.gov/.

## REFERENCES

*arXiv:1710.05941*

*PhD Thesis*