Hybrid Wavelet Neural Network Approach for Daily Inflow Forecasting Using Tropical Rainfall Measuring Mission Data

To get a representative precipitation data set, the 169 TRMM rainfall time series were analyzed according to a hierarchical cluster analysis (Cattell 1943). Hierarchical cluster analysis groups data over several scales by creating a cluster tree, which is known as dendrogram (Fig. 2). Here the agglomerative strategy for hierarchical clustering was used, which is a bottom up approach, i.e., each rainfall time series starts in its own cluster, and pairs of clusters are merged as one moves up the hierarchy. In Fig. 2, the numbers along the horizontal axis correspond to the indices of the clusters in the original data set. The upside-down U-shaped lines represent the links between objects, and the height of such a U-shaped line indicates the distance between the objects. Thus, the link representing, for example, the cluster containing Clusters 6 and 10 has a height of 0.403; whereas, the link representing the cluster that groups Object 1 together with Objects 2, 3, 4, and 9 has a height of 0.508. The height represents the distance linkage computed between clusters, which is calculated here as a correlation distance $d_{st}$ (Székely and Rizzo 2014) between individuals $s$ and $t$, constructed from values of $n$ variables, translated by the vectors $x_s$ and $x_t$, i.e., one minus the sample correlation between points (treated as sequences of values)where ${\overline{x}}_s=\frac{1}{n}\sum _j{x}_{sj}$ and ${\overline{x}}_t=\frac{1}{n}\sum _j{x}_{tj}$, i.e., ${\overline{x}}_s$ = average of each vector $x_{sj}$; and ${\overline{x}}_t$ = average of each vector $x_{tj}$.

This is the task of grouping the time series in such a manner that ones in the same cluster (group) are more similar (in the sense of correlation distance) to each other than to those in other clusters (groups). It is a common technique for statistical data analysis, which is used in many fields, and can be considered as a main task in exploratory data mining. The tree is not a single set of groups but rather a multilevel hierarchy whereby groups at one level are joined as groups at the next level.

Such a procedure allows one to decide the level or scale of grouping that is most suitable for the application. Based on the size of the basin and on the mean annual rainfall distribution, the 169 rainfall time series were grouped into 10 clusters [Fig. 3(a)]. Finally, an average daily rainfall time series over each cluster was calculated based on the Thiessen method (Boots 1999), i.e., the TRMM grid point weights based on the relative areas of each measurement grid in the Thiessen polygon network were calculated [Fig. 3(b)]. The individual weights were multiplied by the grid TRMM rainfall estimates, and the values were summed to obtain the areal average precipitation time series for each cluster, which will be referred to as raw rainfall data.

An artificial neural network (ANN) could be defined as a set of simple processing units (neurons), which work as a parallel distributed processor. Such neurons are responsible for storing experimental knowledge for further disposal. The biological nervous system inspired the ANNs, and they learn throughout examples, as the brains, and store the acquired knowledge in the connection weights among the neurons.

In the present study, a feed-forward network with one hidden layer was selected, in which the input data ($x_1,x_2,\dots ,x_n$) are included in the first layer, and the network progressively processes those data throughout subsequent layers to produce the results ($y_1,y_2,\dots ,y_k$) in the output layer. The input neurons are linked to those in the intermediate layer by $w_{ji}$ weights (weight connecting the $i$th neuron in the input layer and the $j$th neuron in the hidden layer), and the neurons in the intermediate layer are linked to those in the output layer by $w_{kj}$ weights (weight connecting the $k$th neuron in the output layer and the $j$th neuron in the hidden layer). The ANNs, based on the nonlinear activation functions, map the relationship between the inputs and the output. Thus, the explicit correlation for the output values is expressed bywhere $f_h$ = activation function of the nodes in the hidden layer; $f_o$ = activation function of the nodes in the output layer; $s$ and ${\mathrm{s}}^{\prime }$ = number of nodes in the input and hidden layers, respectively; $b_j$ = bias for the $j$th hidden neuron; and $b_k$ = bias for the $k$th output neuron.

The ANN weights were obtained by calibration, when several sets of input–output are used in the ANN, and the weights are iteratively modified until the output values differ little from the target outputs. As the Levenberg–Marquardt (LM) algorithm is assumed to be one of the fastest methods for training ANNs, it was chosen here. The LM algorithm takes advantage of internal recurrence to dynamically incorporate past experience in the training process (Krishna et al. 2012).

A discrete wavelet transform is a signal filtering technique that could also be used in raw hydrologic data. Usually, it is more reasonable to calculate only a subset of scales and positions to reduce the quantity of data if the wavelet coefficients are calculated at every possible scale. It turns out that if scales and positions based on powers of two are chosen then the analysis will be much more efficient and accurate. Thus, practical filters are used to fast decompose the signal into high-frequency and low-frequency components.

The decomposition process can be performed with successive approximations. Consequently, the original time series can be broken down into many lower-resolution components. Hence, in the filtering process, at the most basic level, the raw signal passes throughout two complementary filters (high-pass and low-pass) and turns into two signals, which are called details and approximations, whereby the details ($D$) are the low-scale/high-frequency components, and the approximations ($A$) are the high-scale/low-frequency components of the signal. One may note that the raw signals ($Q$ and $P$) are formed by those $A$ and $D$ components. Fig. 4 shows exactly this aspect, up to five levels of decomposition, that the raw inflow and TRMM rainfall signals ($Q$ and $P$) are formed by the sum of approximation (left side of Fig. 4) and details (right side of Fig. 4), which could be considered as the signal noise in some cases. Then, at the first level, the raw signal is decomposed in $A_1$ and $D_1$, e.g., $Q=A_1+D_1$. However, at the second level, $Q$ is equal to the sum of $A_2$, $D_2$, and $D_1$. Thus, as $A_1=A_2+D_2$, $A_2=A_3+D_3,\dots ,A_{n-1}=A_n+D_n$, one may note that $Q=A_1+D_1$, $Q=A_2+D_2+D_1$, $Q=A_3+D_3+D_2+D_1,\dots ,Q=A_n+D_n+D_{n-1}+\cdots +D_1$. This approach produces the wavelet transform of the input data at all dyadic scales. However, rather than relying on an upsampling procedure, this relies more on downsampling, which is an excellent technique for denoising. Krishna et al. (2012) stated that, in general, DWT could be used to smooth and explore denoise time series to improve forecasting, and Quiroz et al. (2011) also used such a technique to reconstruct daily rainfall from rain-gauge data and the normalized difference vegetation index based on the fact that both signals are proportional and periodic.

There are several mother-wavelets that could be used in such a procedure of transformation. In this work, 22 mother-wavelets were tested: (1) Daubechies (db1, db2, db3, db4, db5, db6, db7, db8, db9, db10); (2) Symlets (sym2, sym3, sym4, sym5, sym6, sym7, sym8); and (3) Coiflets (coif1, coif2, coif3, coif4, coif5). Although all tested mother-wavelets improved the performance of the ANN models, the coif5 ($N=5$) presented the best results, presumably because the Coiflets are discrete wavelets with scaling functions and vanishing moments, they are near symmetric, and the wavelet functions have $N/3$ vanishing moments and the scaling functions have $N/3–1$, in which $N$ is the wavelet index set here to 5. Thus, only the results of the hybrid models with coif5 were used to be compared with the regular ANN.

Two models, called henceforth as regular ANN and WA-ANN, with three variants each, were built to forecast inflows seven days ahead because ONS has been working with such a lead time. Regular ANN model stands for the models that use as input the raw rainfall ($P_t$) of each cluster and/or raw inflow data recorded in the current day ($Q_t$) and in the four previous days ($P_{t-1}$, $P_{t-2}$, $P_{t-3}$, $P_{t-4}$, $Q_{t-1}$, $Q_{t-2}$, $Q_{t-3}$, $Q_{t-4}$) while WA-ANN model stands for the models that use as input data the subseries of approximations ($A$) of those raw data obtained through a multiresolution analysis (considered also as filtered signal). Several quantities of antecedent days were tested, but it is noted that when increasing them to more than four days, there was a loss in computer time performance and no improvement was observed in the results. The decomposition was performed in 10 levels ($A_1,D_1,\dots ,A_{10},D_{10}$), for each situation, and the correlation between the raw data (inflow and precipitation) and their approximations and details were respectively calculated (Table 2), as well as the correlation between the observed inflows and the calculated ones using the inflow approximations ($A_1,\dots ,A_{10}$) as input data for the WA-ANN model (Table 3). For example, in Table 2, the first value of the first row (0.995) is the correlation between raw inflow data and its approximation $A_1$ time series, and the fifth value of the fourth row (0.471) is the correlation between the raw rainfall data of Cluster #2 and its approximation $A_5$, and the last value (0.007) is the correlation between raw rainfall data of Cluster #10 and its details $D_{10}$. Thus, the approximation at level three was selected ($A_3$) because it shows simultaneously a correlation greater than using level four ($A_4$) of both inflow and rainfall data (Table 2) and the best results when using the transformed inflows as inputs of the WA-ANN (Table 3). Note also that the high-frequency components (details) present low correlation with the original signal (Table 2). Therefore, the proposed WA-ANN will use as input values the $A_3$ of inflows and TRMM rainfall data for each cluster in the four previous days and in the current day ($A_{3t},A_{3t-1},A_{3t-2},A_{3t-3},A_{3t-4}$). For illustration, Fig. 4 shows the approximations and details of the inflow and TRMM data up to level five. In both models (regular ANN and WA-ANN), the output layer will have only one neuron, which is the forecasted inflow seven days ahead ($Q_{t+7}$).

The initial objective, when an ANN is designed, is in finding an optimal architecture that allows for capturing the relation among the input and output variables. However, there is not a general rule to indicate the quantity of neurons for the input and hidden layers. Usually, the number is achieved through a trial-and-error procedure. Therefore, the quantity of neurons in the hidden layer was increased from 2 to 30, with an increment of two unities, and based on the root-mean-square deviation of the training set, the quantity of neurons in the hidden layer was finally set in 20 neurons to keep the same architecture for the six analyzed models. Adamowski and Sun (2010) had used 22 neurons in their study to develop a hybrid wavelet transform and artificial neural network model for flow forecasting of nonperennial rivers in semi-arid basins.

Yonaba et al. (2010) tested the Elliott sigmoid, the bipolar sigmoid, and the tan-sigmoid transfer functions to propose multilayer perceptron artificial neural networks for multistep ahead inflow forecasting over five different river basins and lead times, and the results showed that the tan-sigmoid was the most pertinent transfer function for inflow forecasting. Also, the results confirmed that the universal approximation theorem, i.e., a linear transfer function, is suitable for the output layer. Thus, the tan-sigmoid function was selected as the activation function for the hidden neurons ($f_h$), and a linear activation function ($f_o$) for the output layer neuron. The model performance during the validation was evaluated by using the correlation coefficientwhere $Q_{\mathrm{o}}$ = observed inflow; ${\overline{Q}}_o$ = mean of the observed inflow; and $Q_{\mathrm{c}}$ = calculated inflow.

To forecast the inflows seven days ahead into the Três Marias reservoir, the re-naturalised daily inflows and rainfall data of 15 years (January 1, 1998 to December 31, 2012) were used, a total of 5,468 records for each time series (i.e., 5,479 days minus 7 days for forecasting and minus 4 days for input), of which 70% was used for training, 15% for validation, and 15% for testing. Thus, six models were built with different input data set, namely ${\mathrm{M}}_1$ using raw TRMM rainfall data as input, ${\mathrm{M}}_2$ using raw inflow, ${\mathrm{M}}_3$ using the approximation at level three ($A_3$) of rainfall data, ${\mathrm{M}}_4$ using $A_3$ of inflow, ${\mathrm{M}}_5$ using raw TRMM rainfall and raw inflow data, and ${\mathrm{M}}_6$ using $A_3$ of TRMM rainfall and $A_3$ of inflow data (see the two first columns in Table 4).

Several statistics that describe the degree of similarity among the data forecasted by the model and those observed can be used to assess the model efficiency (Santos and Silva 2014). In this work, besides the correlation coefficient, five indices were used, the Nash–Sutcliffe model efficiency coefficient (NASH), the percent bias (PBIAS), the standard deviation ($\sigma$), the root-mean-square deviation (RMSD), as defined and discussed by Moriasi et al. (2007) and the mean absolute percentage error (MAPE) as defined by Renno et al. (2015)where $n$ = number of records in the time series; and $\overline{Q_c}$ = mean of the calculated inflow time series. The Nash–Sutcliffe efficiency is a normalized statistic that is used to determine the relative magnitude of the residual variance compared to the observed data variance. It indicates how the scatter plot for the observed and simulated data well fits the equality line. It ranges between $-\infty$ and 1.0, and the optimal value is $NASH=1$. Values between 0.0 and 1.0 are usually acceptable levels of performance, whereas unacceptable levels of performance are for values less than 0.0, which means that the mean observed value is a better predictor than the simulated value. Pearson’s correlation coefficient ($r$) indicates the degree of collinearity between predicted and measured values. This coefficient ranges between $-1$ and 1, and if $r$ is equal to 1 or $-1$, a perfect positive or negative linear relationship exists, whereas if $r$ is equal to 0, no linear relationship exists. Percent bias measures the average tendency of the simulated values to be larger or smaller than the respective observed values. Then, the optimal value of PBIAS is 0.0%; therefore, the low-magnitude values would indicate an accurate forecasting. Thus, positive values would indicate model underestimation bias, whereas the negative values would indicate model overestimation bias. The root-mean-square deviation computes the standard deviation of the model prediction error. The smaller the RMSD value, the better the model performance; however, it has limited ability to clearly indicate poor model performance, and the mean absolute percentage error usually expresses forecast error as a percentage, and the smallest values show the best forecasting.