Improving the prediction accuracy of river inflow using two data pre-processing techniques coupled with data-driven model

SSA for denoising

For time series analysis, the SSA method is known as a powerful non-parametric method (Golyandina, Nekrutkin & Zhigljavsky, 2001). SSA combines the principals of time series analysis, multivariate statistics, dynamical and signal processing (Suhartono et al., 2018). The reason for using SSA as it is a model-free technique (Romero et al., 2015), which can be applied on any type of data without any assumption. The main function of SSA is to decompose the time series data into a trend, seasonal, oscillations, and aperiodic noises and then reconstruct it after removing aperiodic noises from time series data (Traore et al., 2017). Unlike other methods of time series analysis, SSA assumes no statistical assumption about noises while performing analysis and investigating its properties (Traore et al., 2017).

VMD as decomposition

The VMD is a non-recursive signal decomposition estimation method introduced by (Dragomiretskiy & Zosso, 2014; Rezaie-Balf et al., 2019a; Rezaie-Balf et al., 2019b). The VMD adaptively decomposes complicated original non-linear, non-stationary and multi-scale signals into band-limited IMFs i.e., u_k with a specific bandwidth in the spectral domain. To achieve a bandwidth of each IMF, the constrained variational optimization problem is solved as follows (Dragomiretskiy & Zosso, 2014): (5)${\displaystyle \underset{\left\{u_k\right\},w_k}{min}\begin{array}{c}\hfill \left\{\sum _k{∥{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-jwkt}∥}_2^2\right\}\hfill \\ \hfill s.t\sum _k{u}_k=f\hfill \end{array}}$ where w_k is the center frequency of kth IMF, $\delta \left(t\right)$ is the Dirac function, t is the time script and k is the number of modes. Moreover, $\left(\delta \left(t\right)+\frac{j}{\pi t}\right)$ is the Hilbert transformation function which transform u_k into an analytical signal to form a one-side frequency. The spectrum of each mode can be shifted to a base-mode with the e^−jwkt term. The above constrained problem is converted into unconstrained by making use of a quadratic penalty term i.e., α and Lagrangian multipliers λ, which is easier to address described as follows: (6)${\displaystyle L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \right)=\alpha \sum _k{∥{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-jwkt}∥}_2^2+{∥f\left(t\right)-\sum _{k=1}^K{u}_k\left(t\right)∥}_2^2+〈\lambda \left(t\right),f\left(t\right)-\sum _{k=1}^K{u}_k\left(t\right)〉}$ where α denotes balancing parameter. The Eq. (5) can also be solved by an alternate direction method of multipliers. It is implied that updating u_k ω_k and λ_k in two directions is conducive for realizing the analysis process of VMD, and the solutions of u_k ω_k and λ_k can be calculated as follows: (7)${\displaystyle u_k^{n+1}\left(\omega \right)=\frac{\hat{f}\left(\omega \right)-\sum _{i\ne k}{\hat{u}}^n\left(\omega \right)+\left(\frac{\hat{\lambda }\left(\omega \right)}{2}\right)}{1+2\alpha {\left(\omega -{\omega }_k\right)}^2}}$ (8)${\displaystyle {\omega }_k^{n+1}=\frac{{\int }_0^{\infty }\omega {\left|u_k^{n+1}\left(\omega \right)\right|}^{2}d\omega }{{\oint }_0^{\infty }{\left|u_k^{n+1}\left(\omega \right)\right|}^{2}d\omega }}$ and (9)${\displaystyle {\hat{\lambda }}_k^{n+1}\left(\omega \right)={\hat{\lambda }}_k^n\left(\omega \right)+\tau \left(\hat{f}\left(\omega \right)-\sum _{i\ne k}{\hat{u}}^n\left(\omega \right)\right)}$ where $\hat{f}\left(\omega \right)$, $u_k^{n+1}\left(\omega \right),u_k^n\left(\omega \right)$ and $\hat{\lambda }\left(\omega \right)$ are the Fourier transforms of f and n denotes the number of iterations. The termination condition of VMD is defined as follows: (10)${\displaystyle \frac{\sum _k{∥u_k^{n+1}-u_k^n∥}_2^2}{{∥u_k^n∥}_2^2}<ϵ}$ where ϵ is the tolerance level of the convergence criterion.

From VMD, the IMF u_k is obtained from the entire decomposition process according to the following steps:

1. Set iteration number n = 1, and initialize parameters for VMD including $u_k^1$ , ω¹ and λ¹.

2. Using the Eqs. (8) and (9), calculate $u_k^{n+1}\left(\omega \right)$ and ${\omega }_k^{n+1}$.

3. After calculating $u_k^{n+1}\left(\omega \right)$ and ${\omega }_k^{n+1}$, update Lagrangian multiplier using the Eq. (9).

4. If the convergence condition of Eq. (10) is met, the iteration will be stopped, otherwise n moves to n + 1, and again return to step 2. Finally, the IMF are obtained.

Support vector machine (SVM) as a prediction method

The SVM is a supervised machine learning method that comprised of statistical learning principles for nonlinear classification, function estimation, and pattern recognition applications (Vapnik, 1998). After introducing loss function, SVM can be used as a time series forecasting as well (Yaseen, Kisi & Demir, 2016;Sanghani, Bhatt & Chauhan, 2018). The concept behind the SVM is that it maps the complex non-linear high dimensional data into a high feature space through a nonlinear mapping. After mapping data into a high feature space, linear regression is performed by SVM in that feature space. Let we have a training set consists of N sample points, ${\left\{x_i,y_i\right\}}_i^N$, where x_i is lagged input vector and y_i is the estimated value of a time series data, then the SVM is formulated as follows: (13)${\displaystyle Y=f\left(x\right)=w^T\varphi \left(x\right)+b}$ where $\varphi \left(x\right)$ is a non-linear transfer function projecting the input data into high dimensional space, w_i are the weight vectors and b_i is a bias. Estimating the sampled values with the range of allowed precision is considered as the problem of finding the minimum value for $∥w∥$. This can be summarized as convex programming: (14)${\displaystyle min\left(\frac{∥w^2∥}{2}+C\sum _i^N\left(\xi +{\xi }^{\ast }\right)\right)}$ subject to (15)${\displaystyle \left\{\begin{array}{c}f\left(x_i\right)-y_i\le \epsilon +{\xi }^{\ast }\hfill \\ y_i-f\left(x_i\right)\le ϵ+\xi \hfill \\ \xi ,{\xi }^{\ast }\ge 0\hfill \end{array}\right.}$ where C is the user-defined penalty coefficient which represents the dispersion between the weights and objective function. The ξ and ξ^∗ termed as the slack variables which describes how much data exceeds from tolerance. The Lagrangian function is further applied that uses regression function to replace weight vector and $\varphi \left(x\right)$ given in Eq. (13) as follows: (16)${\displaystyle f\left(x\right)=\sum _i\left({\alpha }_i-{\alpha }_i^{\ast }\right)K\left(x,x_i\right)+b}$ where α_i and ${\alpha }_i^{\ast }$ are the Lagrangian multipliers and K is called the kernel function. The possible tested kernels includes linear, polynomial, Gaussian and sigmoid kernels which are defined respectively as follows: ${\displaystyle K\left(x,x_i\right)=x.x_i}$ ${\displaystyle K\left(x,x_i\right)={\left(\gamma \left(x.x_i\right)+r\right)}^d}$ ${\displaystyle K\left(x,x_i\right)=exp\left(-\gamma {\left|x-x_i\right|}^2\right)}$ (17)${\displaystyle K\left(x,x_i\right)=tanh\left(\gamma \left(x.x_i\right)+r\right)}$ where γ is the structural parameter, d is a polynomial degree and r represents the residuals of the system. Different values of , γ and penalty parameter C is used in this study. The quadratic structure of the Eq. (16) is defined as: (18)${\displaystyle W\left({\alpha }_i,{\alpha }_i^{\ast }\right)=\sum _i{y}_i\left({\alpha }_i-{\alpha }_i^{\ast }\right)-\epsilon \sum _i\left({\alpha }_i+{\alpha }_i^{\ast }\right)-\frac{1}{2}\sum _i\sum _j\left({\alpha }_i-{\alpha }_i^{\ast }\right)\left({\alpha }_j-{\alpha }_j^{\ast }\right)K\left(x_i,x_j\right)}$ With the following constraints: ${\displaystyle \sum _i^N\left({\alpha }_i-{\alpha }_i^{\ast }\right)=0}$ ${\displaystyle 0\le {\alpha }_i\le C,i=1,2,\dots ,N}$ (19)${\displaystyle 0\le {\alpha }_i^{\ast }\le C,i=1,2,\dots ,N.}$

Evaluation assessment methods

We assessed and compared the predction performance of our proposed hybrid model SSA-VMD-EBT-SVM with other existing models (EMD-SVM, EEMD-SVM, VMD-SVM, SSA-EMD-SVM, SSA-EEMD-SVM, SSA-VMD-SVM) as a benchmark using following four measures: the Nash-Sutcliffe Efficiency (NSE), Mean Square Error (MSE), Root Mean Square Error (RMSE) (Ghorbani et al., 2018) and Mean Absolute Error (MAE) (Yaseen et al., 2018) with following equations respectively; (20)${\displaystyle NSE=1-\left[\frac{\sum _{t=1}^N{\left(y_{ot}-y_{pt}\right)}^2}{\sum _{t=1}^N{\left({\bar{y}}_{ot}-{\bar{y}}_{ot}\right)}^2}\right]}$ (21)${\displaystyle MSE=\frac{\sum _{t=1}^N{\left(y_{ot}-y_{pt}\right)}^2}{N}}$ (22)${\displaystyle RMSE=\sqrt{\frac{\sum _{t=1}^N{\left(y_{ot}-y_{pt}\right)}^2}{N}}}$ (23)${\displaystyle MAE=\frac{\sum _{t=1}^N|y_{ot}-y_{pt}|}{N}}$ where y_ot is the observed values, ${\bar{y}}_{ot}$ is the mean of observed values and y_pt is predicted value of model. Moreover, Taylor diagram is used to prepare a visual comprehension with the help of polar plot for the evaluation of modeling results. The Taylor diagram represents the normalization statndard deviation between simulated and observed values with normalized origin and R² are represented as directional angles (Darbandi & Pourhosseini, 2018). The interpretation of Taylor diagram is that an observed point is shown on graph and the closer the simulated performance measures to the observe point, the better the model performance (Al-Sudani, Salih & Yaseen, 2019).

Benchmark models for the evaluation of the proposed hybrid model

The proposed hybrid model i.e., SSA-VMD-EBT-SVM is compared with six benchmark models described as follows:

1. Without denoising: this type of existing models comprised on decomposition and prediction stages only in which VMD and two different data decomposition methods i.e., EMD, and EEMD are chosen which decompose non-linear, non-stationary and multi-scale data into multiple IMFs with the different sound of time-frequency components. For prediction, the extracted IMFs through EMD, EEMD, and VMD are predicted with the same prediction method i.e., SVM as used in our proposed hybrid model. Then, the performance of proposed hybrid model i.e., SSA-VMD-EBT-SVM is compared with existing benchmark models i.e., EMD-SVM (Yu et al., 2017) and EEMD-SVM (Rezaie-Balf et al., 2019b) and VMD-SVM (Wu & Lin, 2019).

2. With denoising: these models use denoising-decomposition and prediction stages to predict river inflow data. For denoising, SSA is selected with the same decomposition and prediction stages as described in (a). Then, the performance of the proposed hybrid model i.e., SSA-VMD-EBT-SVM is compared with existing benchmark models i.e., SSA-EMD-SVM, SSA-EEMD-SVM, and SSA-VMD-SVM.

Data

The daily river inflow dataset used in this study is comprised on 1st January to 31st March for the period of 2015–2019. To the application of proposed objective, the daily inflow of Indus River at Tarbela with its two principal left and one right bank tributaries: Jhelum River at Mangla, Chenab River at Marala and Kabul River at Nowshera respectively are selected. The daily inflow data is measured in 1,000 ft/s which was acquired from the site of Pakistan Water and Power Development Authority (WAPDA).

Discussion

In this article, we proposed a novel hybrid model to efficiently predict the river inflow time series data. Our proposed method comprised of SSA as denoising, VMD as a data decomposition with EBT threshold, and SVM as a prediction method.

In order to understand the applicability of our proposed model i.e., SSA-VMD-EBT-SVM, two different model strategies are adapted (see Tables 3 and 4). First, without denoising model strategy is implemented on which the same decomposition method i.e., VMD and two different decomposition methods i.e., EMD and EEMD are used. Moreover, for prediction purpose, SVM which is also adapted in proposed methodology is used here to compare the performance of proposed denoised-decomposed strategy with only ‘without denoising models i.e., VMD-SVM (Wu & Lin, 2019), EMD-SVM (Yu et al., 2017) and EEMD-SVM (Rezaie-Balf et al., 2019b). From Tables 3 and 4, it can be seen that the overall performance of without denoising models are poor for all river inflow data with high NSE, MSE, RMSE, MAE and MAPE values. Specifically, EMD-SVM performed worst among all without denoising models due to the fact that EMD suffers from the mode mixing problem and fail to produce noiseless IMFs (Di, Yang & Wang, 2014). The proposed model performed well with the lowest NSE, MSE, RMSE, MAE and MAPE values as compared to without denoising models.

The second strategy used the concept of denoising and decomposition, here the same method of denoising which is used in proposed methodology is employed with different decomposition methods i.e., SSA-EMD-SVM, SSA-EEMD-SVM and with same decomposition method but without thresholding is adapted i.e., SSA-VMD-SVM. From Tables 3 and 4, it is observed that the proposed model i.e., SSA-VMD-EBT-SVM performed well for Indus, Jhelum, and Chenab, but for Kabul river inflow the results of SSA-VMD-EBT-SVM and SSA-VMD-SVM are same as thresholding the IMF did not enhance the prediction performance of SVM. Moreover, SSA-EEMD-SVM also performs well among existing with and without denoising methods for all case studies. Overall, the proposed SSA-VMD-EBT-SVM model showed a much better agreement between predicted and observed river inflow data which demonstrates the suitability of SSA, VMD, and EBT in pre-processing inputs/output data over other decomposition methods i.e., EMD and EEMD. Thus, it is concluded that the appropriate way of denoising, decomposition and thresholding can effectively enhance the performance of non-linear, non-stationary and multi-scale time series data.

By applying proposed simulation models in IRB, it is expected that this will provide new tools for improving inflow prediction over what is possible with the current generation of statistical models as well as help with other land and water management questions. It may also be helpful in setting policies regarding what appropriate methods should be chosen for ‘denoising-decomposition and prediction’ and in assessing the effects of climate warming. These modeling efforts are therefore significant both for the scientific issues involved as well as for the practical relevance of the results.