Hybrid wavelet transform – MLR and ANN models for river flow prediction: Case study of Brahmaputra river (Pancharatna station)

Wavelet theory is discussed thoroughly in works by Mallat (1998), and Labat, Ababou and Mangin (2000). Wavelet transform of a raw signal has the ability of providing time-frequency, using a range of window sizes. It divides the input signal into wavelets (small waves), which are scaled and shifted versions of the original wavelet, called mother wavelet. There are two forms of wavelet transformation: the continuous wavelet transform (CWT) and the discrete wavelet transform (DWT). The CWT of a signal x(t) is described: (1) CWT(a,b)=1a∫−∞∞x(t)⋅ψ*t−ba⋅dt. CWT(a,b) = {1 \over {\sqrt a}}\mathop \smallint \nolimits_{- \infty}^\infty x(t) \cdot {\psi ^*}\left({{{t - b} \over a}} \right) \cdot {\rm{d}}t.

In Eq. (1), the transformed signal is a function of two variables, a and b, which represent the scale and translation factor of the function ψ(t), respectively; * corresponds to complex conjugate (Mallat, 1998).

The ψ(t) is the transforming function, called as mother wavelet, is mathematically defined as: (2) ∫−∞∞ψtdt=0. \mathop \smallint \nolimits_{- \infty}^\infty \psi \left(t \right){\rm{d}}t = 0.

From a data set of length (L), CWT produces L² coefficients. Therefore, duplicate data is trapped within the coefficients, which may or may not be a desired quality (Nourani et al., 2009; Rajaee, Nourani, Zounemat-Kermani & Kisi, 2011). The analysis will be significantly more precise and efficient, resulting in N transform coefficients, if scales and positions are chosen based on the powers of two. The resulting transform is called DWT, and has the form: (3) ψm,nt=1aomψt−nboaomaom, {\psi _{m,n}}\left(t \right) = {1 \over {\sqrt {a_o^m}}}\psi \left({{{t - n{b_o}a_o^m} \over {a_o^m}}} \right), where the wavelet translation and dilation are each controlled by an integer, m, and the other by n; b_o which must be larger than 0, is the location parameter; a_o is a specified fixed dilation step greater than 1. The most common and simplest choice for parameters a_o and b_o are 2 and 1 (time steps), respectively. This power of two logarithmic scaling of the translations and dilations is known as the dyadic grid arrangement (Mallat, 1989). For parameters a_o and b_o, 2 and 1 (time steps), respectively, are the most typical options. According to Mallat (1989), this power-of-two logarithmic scaling of the translations and dilations is referred to as the dyadic grid layout.

Two sets of functions, known as high-pass (wavelet function) and low-pass (scaling function) filters, are operated by DWT. The original time series undergo processing via high-pass and low-pass filters (as shown in Fig. 3), followed by a down-sampled process that discards every second data point (Deka & Prahlada, 2012). The detailed (D₁, D₂, ..., D_n) and approximation (A₁, A₂, ..., A_n) coefficients, which represent the low frequency and high frequency components of the original signal, respectively, are derived after the signal is passed through high pass and low pass filters. There will be a total of n + 1 coefficients at every given n^th decomposition level, where there will be one series of approximation coefficients at the n^th level (A_n) and n series of detailed coefficients (D₁, D₂, ..., D_n). The sum of A_n + D₁ + D₂ +... + D_n is equal to the original signal.

FIGURE 3.

Wavelet decomposition tree

Source: Khandekar (2014).

Neural networks are interconnected groups of artificial neurons that may be utilized as a computational model for information processing based on the connectionist approach to computing. An ANN may be thought of mathematically as a universal approximator with the capacity to learn from instances without the need of explicit physics.

Three-layer-feedforward artificial neural networks are most frequently utilized in hydrologic time series modeling (Jain & Chalisgaonkar, 2000; Raghuwanshi, Singh & Reddy, 2006; Tayfur & Singh, 2006; Bajirao, Kumar, Kumar, Elbeltagi & Kuriqi, 2021; Katipoǧlu, 2023a). In the current study, the Lavenberg–Marquardt (LM) learning function and tangent sigmoid as the transfer function were employed in the feedforward ANN. The LM technique was used to train the ANN since it is more effective and faster than the traditional gradient descent technique (Kisi, 2009b; Mohseni & Muskula, 2023). In three-layer-feedforward ANN, information flow takes place from the input to the output side. Weights regulate the strength of a signal that passes across connections between neurons, altering the information that is transferred between them. While not connected to nodes in the same layer, nodes in one layer are connected to those in the next one. As a result, a node's output in a layer depends exclusively on the inputs and weights it receives from the previous layers. Every node multiplies each input by its weight, adds up the result, and then runs the total through a transfer function to get the desired outcome. Typically, this transfer function is represented as a sigmoid function, which is an S-shaped curve that increases continuously. The “S” upper and lower limb attenuation keeps the raw sums smoothly within predetermined bounds. The transfer function additionally adds a nonlinearity that improves the network's capacity to represent complicated functions (Jain & Chalisgaonkar, 2000). The sigmoid function is monotonically growing, continuous, and differentiable everywhere.

The MLR is one of the statistical methods which represents a mathematical equation expressing the relation between two or more explanatory variable and a response variable. The MLR attempts to model the link between two or more explanatory variables and a response variable by fitting a linear equation to the observed data (Sharma, Isik, Srivastava & Kalin, 2013). The relationship between the dependent variable y and the independent variable x is described by the MLR equation as follows: (4) y=A0+A1x1+A2x2+⋯+Anxn y = {A_0} + {A_1}{x_1} + {A_2}{x_2} + \cdots + {A_n}{x_n} where: x_n is the value of n^th predictor, A₀ is the regression constant, A_i is the coefficient of n^th predictor.

For deciding the number of input parameters, no fixed rule is available. The values of the future time step (Q_t+n) in any time-series forecasting are inevitably reliant on the antecedent values Q_t, Q_t−1, Q_t−2, …, Q_t−j (Table 2). But because the value of j (lag) is not known beforehand, it is difficult to determine how many lags in the past will lead to higher efficiency. In hydrologic time series forecasting, determining j is crucial since it can help to reduce information loss and the omission of critical input variables that could interfere with training.

Numerous researchers (Sudheer, Gosain & Ramsastri, 2002; Budu, 2013; Nayak et al., 2013; Khazaee Poul et. al., 2019; Sun, Niu & Sivakumar, 2019; Katipoǧlu, 2023a, 2023b, 2023c) have used the technique based on the statistical features such as cross-correlation, autocorrelation and partial correlation of the data series in order to identify a distinct input vector. In the current research, based on the autocorrelation coefficient between the relevant variables (shown in Table 2), the input vectors to the models are chosen. To forecast the value of discharge, based on the following various input combinations were taken into account (Khandekar, 2014):

• Q_t.

• Q_t, Q_(t−1).

• Q_t, Q_(t−1), Q_(t−2).

• Q_t, Q_(t−1), Q_(t−2), Q_(t−3).

• Q_t, Q_(t−1), Q_(t−2), Q_(t−3), Q_(t−4).

• Q_t, Q_(t−1), Q_(t−2), Q_(t−3), Q_(t−4), Q_(t−5).

At the first stage, to forecast discharge, ANN and MLR models without data preprocessing were created. In the present study it was decided to predict discharge for multiple lead times such as two days, four days, seven days and 14 days because the Brahmaputra river is one of the major rivers in India which carries heavy flood during monsoon season, so it is necessary to design flood warning system well in advance. First 70% data was used to train the model and last 30% data was used for testing the model. Initially, a three-layer feedforward backpropagation ANN were applied in the study. For each lead time, an optimal input combination is obtained by providing each of one to six input combination obtained through autocorrelation as input to ANN and by varying number of neurons in the hidden layer from two to 10 through trial and error (Budu, 2013; Moosavi, Vafakhah, Shirmohammadi & Behnia, 2013). The optimal input combination that gave the minimum root mean square error (RMSE) during testing period for each lead time is shown in Table 3. The output layer has only one neuron which is the discharge value for the given lead time. From viewpoint of comparing all models, same input combination is employed for the MLR and all hybrid models. The time series data normalized between zero and one (Nourani et al., 2009) by dividing the discharge value by the maximum one.

After developing the ANN and MLR models, hybrid models were developed. First of all, using DWT, for developing hybrid models, the normalized input data were decomposed into approximation and detail coefficients. Since all hydrological data are observed at discrete time intervals, all hybrid models used DWT to process time series data in the form of approximations and details at various levels so that gross and small features of a signal can be separated (Deka & Prahlada, 2012). Then, to obtain the output at a predetermined lead time, the approximation and detail coefficient were fed as input to ANN and MLR. Without decomposition, the output signals were preserved as normalized original series. The complete flow chart of model development is shown in Figure 4.

FIGURE 4.

Flow chart of model development

Source: own elaboration.

The mother wavelet that will be used depends on the data to be analyzed (Nayak et al., 2013). As per authors' knowledge, there is no comparative study on prediction of hydrological parameter using group of Daubechies wavelets. So, in this research, an irregular wavelet, Daubechies (db) of orders 1 (db1), 2 (db2), 3 (db3), 8 (db8), and 10 (db10), illustrated in Figure 5, was chosen to deal with very irregular signal form. Daubechies wavelets of order N (dbN) are all asymmetric, orthogonal, and biorthogonal. They are compactly supported wavelets with extremal phase and highest number of vanishing moments for a given support width (Misiti, Misiti, Oppenheim & Poggi, 2010). Daubechies wavelet has a support width of 2N – 1. Except for db1 (haar), no Daubechies wavelet has an explicit expression. Wavelets with compact support or a narrow window function are appropriate for a local analysis of the signal.

FIGURE 5.

Daubechies wavelets

Source: Khandekar (2014).

The optimal decomposition level was obtained using the formula l = int [log (L)], (Nourani et al., 2014; Khazaee et al., 2019; Tarate et al., 2021; Katipoǧlu, 2023a, 2023b, 2023c), where l is the decomposition level, L represents the number of time series data, int represents the integer part function, and log represents logarithm with base 10. In the present study, L is 3,650 so l is approximately four. However, in our study, we decomposed raw signal up to fifth level. At any l^th decomposition level, DWT produces l detailed coefficients and one approximation coefficient at l^th level.

The following evaluation measures were used to compare model performance. (5) RMSE=∑i=1L(Qobs−Qcom)2L, RMSE = \sqrt {{{\sum\nolimits_{i = 1}^L {{{({Q_{{\rm{obs}}}} - {Q_{{\rm{com}}}})}^2}}} \over L}}, (6) R2=1−∑i=1L(Qobs−Qcom)2∑i=1L(Qobs−Qobs¯)2, {R^2} = 1 - {{\sum\nolimits_{i = 1}^L {{{({Q_{{\rm{obs}}}} - {Q_{{\rm{com}}}})}^2}}} \over {\sum\nolimits_{i = 1}^L {{{({Q_{{\rm{obs}}}} - Q\overline {_{{\rm{obs}}}})}^2}}}}, (7) MAE=1n∑i=1LQobs−Qcom, MAE = {1 \over n}\sum\nolimits_{i = 1}^L {\left| {{Q_{{\rm{obs}}}} - {Q_{{\rm{com}}}}} \right|}, (8) B=∑i=1LQcomp∑i=1LQobs, B = {{\sum\nolimits_{i = 1}^L {{Q_{{\rm{comp}}}}}} \over {\sum\nolimits_{i = 1}^L {{Q_{{\rm{obs}}}}}}}, (9) SI=RMSE∑i=1LQobs¯, SI = {{RMSE} \over {\sum\nolimits_{i = 1}^L {Q\overline {_{{\rm{obs}}}}}}}, where: RMSE, R², MAE, B, SI, L, Q_obs, Q_com, and Qobs¯ Q\overline {_{{\rm{obs}}}} are root-mean squared error (RMSE), determination coefficient (R²), mean absolute error (MAE), bias (B), scatter index (SI), number of observations (L), observed data (Q_obs), computed values (Q_com), mean of observed data ( Qobs¯ Q\overline {_{{\rm{obs}}}} ), respectively. A R² of 0.9 or more is considered very satisfactory, 0.8 to 0.9 represents a fairly good model, and less than 0.8 is considered unsatisfactory (Dawson & Wilby, 2001).

The results of all models are presented in Tables 4–7. It can be seen from Tables 4–7 that in comparison with MLR models, the ANN model has performed better for all lead times, except for the fourteen-day one. It is observed that the effectiveness of both models decreases while increasing lead time. This could be due to considerable fluctuations in the data around mean values, such as a high standard deviation. Tables 4–7 also indicates the optimal ANN structure (e.g. for a two-day lead time, 2-8-1 means two neurons in the input layer, eight neurons in the hidden layer, and one neuron in the output layer).

The normalized observed data was decomposed using Daubechies wavelets of order 1 (db1), 2 (db2), 3(db3), 8 (db8) and 10 (db10) up to fifth level decomposition, which were fed as input to ANN and MLR, making the models as WANN(dbNli) and WMLR(dbNli), respectively, where N is the order of Daubechies wavelet and i is decomposition level. The performances of these hybrid models only for best decomposition level (low RMSE) are presented in Tables 4–7. Results of Tables 4–7 show that wavelet-based hybrid models perform significantly better than the standalone model.

After comparing the results of hybrid models with respect to wavelet order, it was observed that all models' efficiency is increasing with wavelet order, highest being at 10^th order. For example, for four-day lead time (during testing period), from the Table 4 it was found that for WANN model the value of R² increased from 0.956 (WANN-db1l5) to 0.986 (WANN-db10l5) and for WMLR model the increase is from 0.945 (WMLR-db1l5) to 0.991 (WMLR-db10l5). Similar trend was observed for all lead times. Table 8 shows percent improvement in RMSE values with increase in Daubechies wavelet order from db1 to db10. The average percent improvement in RMSE is found to be 57.19%. Figure 6 shows a sample plot of effect of Daubechies wavelet order on RMSE for four-day lead time.

FIGURE 6.

Effect of Daubechies wavelet order on RMSE (four-day lead time)

Source: own elaboration.

In general, for all lead times, from the analysis it was found that for lower order wavelets (db1, db2, db3), WANN model performance was better as compared to the WMLR model, while for higher order wavelets (db8, db10), WMLR was found to be superior as compared to WANN. A careful study of Tables 4–7 reveals that the forecasting ability of each hybrid model improves with increasing wavelet order. According to the observed time series flow data (Fig. 2), it appears to occur in a high frequency during the monsoon season (June to September, i.e. four months out of the year) and in a low frequency during the other non-monsoon eight months. Wider-support wavelets can capture low frequencies, while wavelets with a smaller support can capture high frequencies. The support width of a Daubechies wavelet of order N is equal to 2N – 1 (Misiti et al., 2010), hence the support widths of db1, db2, db3, db8, and db10, wavelets are one, three, five, 15, and 19, respectively. In brief, compared to db1, db2, db3, and db8 wavelet-based forecast models, the db10 wavelet has a reasonable support and good time-frequency localization properties, which together allow the model to capture both the underlying trend and the short-term variations in the time series. This conclusion is consistent with the findings of (Nourani et al., 2013; Katipoǧlu, 2023b, 2023e), who demonstrated that higher order mother wavelet (db4 and db10) offered substantially better results than lower order Haar (db1) wavelet. Scatter plots for all lead times for best model [WMLR(db10)] are shown in Figure 7. The scatter plots showed that the majority of the points were quite close to the 45° line, with only a few having greater magnitudes of observed flow on the lower side, indicating model underestimate (BIAS slightly below 1.0). WMLR(db10) models performance was very satisfactory (R² > 0.9) (Dawson & Wilby, 2001) for all lead times. Time series plots are shown in Figure 8 for all lead times.

FIGURE 7.

Scatter plot for lead time during different testing period

Source: own elaboration.

FIGURE 8.

Time series plot for lead time during different testing period influence of decomposition level on model performance

Source: own elaboration.

As mentioned earlier, the optimum decomposition level was determined using the formula l = int [log (L)]. In our study, we decomposed raw signal up to fifth level. For the best model WMLR(db10), the effect of decomposition level on R² is presented in Table 9. A careful examination of Table 9 demonstrates that model efficiency enhances with decomposition level. The initial resolution level of the original time series captures the high frequency components, and as the decomposition level (scale) is increased, the signal becomes smoother and more stationary. As a result, prediction errors were not affected by scale. Figure 9 shows effect of decomposition level on R² for WMLR(db10) model for all lead times.

FIGURE 9.

Effect of decomposition level on determination coefficient (R²) for WMLR(db10) model

Source: own elaboration.

Finally, because the Brahmaputra river carries substantial flood during the monsoon season (June to September), an attempt was made to assess the model's accuracy during the monsoon period (for three years in testing from 1997 to 1999) for the WMLR(db10) model. The statistical parameter values for the WMLR(db10) model during testing period of monsoon season are shown in Table 10. As shown in Table 10, the WMLR(db10) model demonstrated good performance during monsoon season producing extremely satisfactory results for all lead times, despite substantial non-stationarity. Figure 10 compares flow series between observed, WMLR(db10) and WANN(db10) modelled flow for monsoon season for two-day lead time. From Figure 10 it is clear that the WMLR(db10) model has captured well almost all peaks, except the highest.

FIGURE 10.

Flow series comparison between observed, WMLR(db10) and WANN(db10) modelled flow for monsoon season

Source: own elaboration.