An improved framework to predict river flow time series data

Decomposition stage

EMD: To decompose non-linear and non-stationary data, the EMD method introduced by Huang et al. (1998) which decomposed data into IMF’s by satisfying two conditions as follows: (a) The number of zero crossings and extreme, from complete data, must be equal or differ to one at most; (b) At all levels, mean value of envelope must be zero.

The complete EMD procedure is defined as:

1. Find all local maxima and minima from data $y\left(t\right)$, ( t = 1, 2, .., N). Use cubic splines interpolation to find an upper envelope from maxima e_max(t) and lower envelope from minima e_min(t).

2. Calculate the average of upper and lower envelope $m\left(t\right)=\left(e_{max\left(t\right)}+e_{min\left(t\right)}\right)∕2$. Take the deviations between original time series data and calculated envelope mean as:

3. (1)${\displaystyle g\left(t\right)=y\left(t\right)-m\left(t\right)}$

4. Match the requirements of $g\left(t\right)$which is defined in (a) and (b) as an IMF, if the conditions are satisfied then mark this g(t) as ith IMF.

5. In the next step replace original time series $y\left(t\right)$ by $r\left(t\right)=y\left(t\right)-g\left(t\right)$, if g(t) does not meet the (a) and (b) then just replace $y\left(t\right)$ with $g\left(t\right)$.

6. The process of (1–4) is repeated until no IMF is being extracted from the residue $r\left(t\right).$

In the end, original time series data will be written as the summation of all extracted IMF’s and residue as: (2)${\displaystyle y\left(t\right)=\sum _i^l{g}_i\left(t\right)+r\left(t\right)}$

where $g_i\left(t\right)$ is the ith IMF and r(t) is the trend of the signal. However, EMD techniques have two drawbacks. The one is the endpoint problem that is the extreme values of two ends of the series can’t be determined properly which distort the IMF’s and the other is mode mixing aliasing in which same IMF contains even more than two frequencies.

EEMD: To resolve the mode mixing problem of EMD, (Wu & Huang, 2004) proposed EEMD. The EEMD decompose the non-linear signals into IMF’s as follows:

1. Add a white Gaussian noise series to the original data set as follows:

2. (3)${\displaystyle y_2\left(t\right)=y\left(t\right)+n\left(t\right)}$

3. Decompose the new $y_2\left(t\right)$ with EMD and obtains the IMFs.

4. Repeat step (a) and (b) mth time i.e., (j = 1,2,…,m) with different white noises to get IMF’s from new series.

5. Find the ensemble means of all IMFs obtained as mthensemble time as follows: where k = 1, 2, .., K is kth IMF.

6. (4)${\displaystyle {\overline{imf}}_k=\sum _{j=1}^{m}IM{F}_{jk}^m∕m}$

where k = 1, 2, .., K is kth IMF.

CEEMDAN: Although the mode mixing problem is alleviated with EEMD, taking the simple averages of IMF’s without considering the independent addition of white noises could not completely remove noises from IMF’s. To tackle the simple averaging problem of EEMD, the CEEMDAN function is introduced by Torres et al. (2011). Here in our study, CEEMDAN is used to decompose the river flow data which is briefly described as follows:

1. The CEEMDAN add Gaussian noises like EEMD in signals as follows:

2. (5)${\displaystyle y_3\left(t\right)=y\left(t\right)+w_0{n}^j\left(t\right)}$

3. where w₀ is the amplitude of the added white noises and (j = 1,2,…,m). Find the first IMF using simple EEMD defined as:

4. (6)${\displaystyle I\widetilde{M{F}_1}=\sum _{j=1}^{m}IM{F}_{j1}^m∕m}$

5. Compute the deviation of original signals from the first IMF as:

6. (7)${\displaystyle r_1\left(t\right)=y\left(t\right)-{\overline{imf}}_1}$

7. Decompose $r_1\left(t\right)+w_0{n}^j\left(t\right)$ to get first IMF and find the second IMF as: (8)${\displaystyle I\widetilde{M{F}_2}=\sum _{j=1}^{m}IM{F}_{j1}^m∕m}$

8. Repeat the (2–3) until stoppage criteria are met and the residual contains not more than two extremes. Finally, the residual is defined as:

9. (9)${\displaystyle R\left(t\right)=y\left(t\right)-\sum _{k=1}^{K}I\widetilde{M{F}_k}}$

However, the selection of a number of ensemble and amplitude of white noise is still an open challenge but here in our study we used a number of ensemble members as 100 and standard deviation of white noise is settled as 0.2 according to Di, Yang & Wang (2014).

Identification of noisiest IMFs: After deriving the IMF’s, next step is to screen out the noise only IMF’s and noise free IMFs (Wei et al., 2016). Two types of IMF’s are derived from CEEMDAN/EEMD; the IMF’s contains high frequency which is corrupted with sparsity and noises (Wei et al., 2016) and the second part of IMF’s comprised on low frequencies which are free from noises (Wei et al., 2016). To get numerical validation, cross-correlation is calculated between all extracted IMF’s and original river flow data. The low cross-correlation implies that the high-frequency IMF’s are overwhelmed with noises. Then, the noise only IMF’s are further denoised through the appropriate threshold to get the important features from them and to get rid of noises.

Denoising stage

After decomposing the non-linear and non-stationary data into IMF components, an appropriate estimator is chosen to remove noises and sparsities from extracted IMFs. The reason for selecting an appropriate estimator is to find an optimal value of threshold as the highest threshold value would lead to biases, whilst the lowest threshold value would increase the noise variance. The IMFs which are extracted are mostly empirical Bayesian threshold estimator is adapted to denoised the noisiest IMFs. Later, the existing soft and hard threshold and improved threshold functions are used for comparison purposes. Details of all estimators are given below:

CEEMDAN/EEMD-based Empirical Bayesian Threshold: To estimate the sparseness and noises from decomposed IMF’s (Wei et al., 2016), an Empirical Bayes Threshold (EBT) inspired from wavelet denoising (Johnstone & Silverman, 2005; Jansen & Bultheel, 2014) is used. For the successful application of EBT, first, all the data is scaled transformed to efficiently select the prior distribution for noises and sparsities. After the scaled transformation data follows $N\left({\theta }_i,1\right)$. Then a mixture of priors for θ_i are considered as follows: (10)${\displaystyle f_{prior}\left(\theta \right)=\left(1-w\right){\delta }_0\left(\theta \right)+w\gamma \left(\theta \right)}$

where ${\delta }_0\left(\theta \right)$ is a zero part of scaled data and $\gamma \left(\theta \right)$ is a density of non-zero part. The density of prior should be chosen in such a way that it must belong to a family of distributions whose tails decays at polynomial rates. The parameters and weights of a mixture of prior distributions are estimated through maximum likelihood approach. The reason for using this method to estimate unknowns is that it estimates weights and parameters to be proportional to the likelihood function evaluated at the estimators based on data (Hossain, Kozubowski & Podgórski, 2018).

Finally, the posterior median ${\tilde{\theta }}_i\left(imf,w\right)$ is calculated from a mixture of prior distribution is given as follows: (11)${\displaystyle \widetilde{F_1}\left(\mu |imf\right)={\int }_{\mu }^{\infty }{f}_1\left(\mu |imf\right)d\mu }$

which is used as a threshold rule for $\tilde{\mu }$ given data (Johnstone & Silverman, 2005). In general, an estimation rule comprised on η (imf,t) defined for all t > 0, is a thresholding rule if and only if for all t > 0, $\eta \left(imf,t\right)$ is an antisymmetric and increasing function of data and η (imf,t)=0 if and only if $\left|imf\right|\le t$ where t is defined as a median value which is calculated through the Eq. (11).

Other traditional threshold estimators: To compare the EBT estimator with other non-linear estimators to suppress the noises and sparsities from noisiest IMFs, Soft Threshold (ST), Hard Threshold (HT) and Improved Threshold Function (ITF) are used which are most widely used in literature (Jeng et al., 2007; Candes, Sing-Long & Trzasko, 2012; Jansen & Bultheel, 2014) listed as follows, respectively;

(12)${\displaystyle IM{F}_{t,k}^{\prime }=\left\{\begin{array}{c}IM{F}_{t,k}\left|IM{F}_{t,k}\right|\ge T_k\hfill \\ 0\left|IM{F}_{t,k}\right|<T_k\hfill \end{array}\right.}$

(13)${\displaystyle IM{F}_{t,k}^{\prime }=\left\{\begin{array}{c}sgn\left(IM{F}_{t,k}\right)\left(\left|IM{F}_{t,k}\right|-T_k\right)\left|IM{F}_{t,k}\right|\ge T_k\hfill \\ 0\left|IM{F}_{t,k}\right|<T_k\hfill \end{array}\right.}$

and

(14)${\displaystyle IM{F}_{t,k}^{\prime }=\left\{\begin{array}{c}sgn\left(IM{F}_{t,k}\right)\left(\frac{\left|IM{F}_{t,k}\right|-T_k}{exp\left\{3\mathrm{\alpha }\left(\frac{IM{F}_{t,k}-T_k}{T_k}\right)\right\}}\right)\left|IM{F}_{t,k}\right|\ge T_k\hfill \\ 0\left|IM{F}_{t,k}\right|<T_k\hfill \end{array}\right.}$

where, T_k is the threshold value calculated as $T_k=a\sqrt{2E_{k}ln\left(N\right)}$, where k = 1,2,…..,K and a is constant takes the values between 0.4 to 1.4 with a step of 0.1 and E_k = median(|IMF_t,k|, t = 1, 2, .., N)∕0.6745 is median deviation of kth IMF.

Prediction and aggregation

In the prediction stage, the decomposed-denoised IMFs of river flow data are predicted through some data-driven and statistical models. Specifically, the denoised IMFs are predicted through a data-driven model whenever noise-free IMFs and the residual are predicted through a simple stochastic model. To train the model, 70%, 80%, and 90% data are used, and the remaining 30%, 20% and 10% data are used to test the accuracy of the models. The selected models are briefly described as follows:

The denoised-IMF prediction with the neural network: A data-driven technique has been proved a powerful tool to model complex non-linear data (Campisi-Pinto, Adamowski & Oron, 2012). The Multi-Layer Perceptron (MLP) which is the most popular sub-model of NN (Talaee, 2014; Ali et al., 2017) consists of three layers of nodes used here for the prediction of denoised IMFs of river flow data. The complete layout of MLP is given in Talaee (2014). The structure of MLP comprised of NN structure. Three nodes are included in MLP as an input layer, hidden layer, and an output layer. First, the single output is calculated by using linear combinations of inputs which are further transferred to some non-linear activation functions mathematically defined as by (15)${\displaystyle y=\phi \left(\sum _i^n{w}_i{x}_i+b\right)}$

where w_i are the weights of inputs, x_i are the inputs, b is the biased value for each layer and φ is a non-linear activation function which supplies output to the next layer. The mostly used activation function i.e., logistic function is used as activation function defined as follows: (16)${\displaystyle \frac{1}{\left(1+e^{-\sum _i^n{w}_i{x}_i+b}\right)}}$

For the optimization of neurons, the supervised learning algorithms called back-propagation and forward-propagation can be used.

The noise-free IMF prediction with ARIMA model: To predict the noise-free IMF’s and the residual, an autoregressive moving average model (ARMA) is selected which is described as follows: (17)${\displaystyle IM{F}_t^k={\alpha }_{1}IM{F}_{t-1}^k+\dots +{\alpha }_{p}IM{F}_{t-p}^k+{\epsilon }_t^k+{\beta }_1{\epsilon }_{t-1}^k+\dots +{\beta }_q{\epsilon }_{t-q}^k}$

where, $IM{F}_t^k$ is the kth IMF, $IM{F}_{t-p}^k$ is pth lag value of kth IMF ${\epsilon }_t^k$ is the residual of kth IMF, p and q are autoregressive and moving average lags. Moreover, in some cases time series data is not stationary. To make such stationary, differences at an appropriate degree are used (Box & Jenkins, 1970). If such a situation occurs, then the model is known as ARIMA (p, d, q) where d is the differenced value used to make the non-stationary data stationary.

Selection of study area

In this research, the largest water system of Pakistan is considered for the application of the proposed strategy. The Indus Basin System (IBS) is the largest river of Pakistan and it plays an important role specifically in power generation and irrigation department. The major tributaries of IBS i.e., River Jhelum, River Chenab, and River Kabul are selected for the present study. Pakistan is facing a large amount of frequent river flooding each year due to monsoon rain and melting snow or glaciers. As in Pakistan, glacier-covered 13,680 km² area which is estimated 13% of the mountainous areas of Upper Indus Basin (UIB). Melted water from these 13% areas adds the significant contribution of water in these rivers. which leads to complex characteristics in river flow data. Therefore, for sustainable economic development and efficient water resources planning, there is a need for analyzing such complex characteristics and predict the behaviors of river flow data at IBS and its tributaries.

Data

To investigate the improved framework, four daily river flow data comprised on (1st-January to 19th-June) for the period of 2015–2018 is used in this study. The main river flow of Indus at Tarbela is considered with its three principal tributaries: Jhelum at Mangla, Chenab at Marala and Kabul at Nowshera. Data is measured in 1,000 ft/s. The described river flow data was obtained from the website of Pakistan Water and Power Development Authority (WAPDA).

Evaluation criteria

Evaluation of noise reduction methods: The performances of denoised series needs comprehensive evaluation after using appropriate noise reduction threshold methods. In this research, to check the performances of CEEMDAN/EEMD-ST CEEMDAN/EEMD-HT, CEEMDAN/EEMD-ITF, and proposed framework i.e., CEEMDAN/EEMD-EBT, Signal-to-Noise Ratio (SNR), Mean Square Error (MSE) and Mean Absolute Error (MAE) (Nazir et al., 2019) are employed which are given as follows respectively; (18)${\displaystyle SNR=10log10\left(\frac{\sum _{t=1}^N{\left(y_{ot}\right)}^2}{\sum _{t=1}^N{\left(y_{pt}-y_{ot}\right)}^2}\right)}$ (19)${\displaystyle MSE=\sqrt{\frac{\sum _{t=1}^N{\left(y_{ot}-y_{pt}\right)}^2}{N}}}$

and (20)${\displaystyle MAE=\frac{\sum _{t=1}^N|y_{ot}-y_{pt}|}{N}}$

where y_ot is the tth observed value and y_pt is the tth predicted value. Moreover, the performance of proposed model (i.e., CEEMDAN-EBT-MM) and all other models including CEEMDAN/EEMD-ST-MM, CEEMDAN/EEMD-HT-MM, CEEMDAN/EEM-ITF-MM and EEMD-ITF-MM are compared using three popular statistical measures: MSE, MAE defined in Eqs. (19)–(20) and Mean Absolute Percentage Error (MAPE) given as follows; (21)${\displaystyle MAPE=\frac{|y_{mo}-y_{mp}|}{\left|y_{mo}\right|}\ast 100}$

where, y_ot and y_pt is defined above, y_mo is the mean value of observed values, y_mp is the mean value of predicted values and y_sp is the standard deviation of predicted values.