Suspended sediment load prediction in consecutive stations of river based on ensemble pre-post-processing kernel based approaches Open Access

Sediment transport is one of the most important issues in river engineering. So far, various and complex relationships have been proposed to predict the amount of suspended sediment transport rate such as velocity and critical shear stress based equations. However, the complex nature of sediment transport and lack of validated models make it difficult to model the suspended sediment concentration and suspended sediment discharge carried by rivers. Bhattacharya et al. (2004) stated that it is difficult to express the transport process through a deterministic mathematical framework. Based on laboratory experiments, Vongvisessomjai et al. (2010) studied the non-cohesive sediment transport in uniform flow at a non-deposition state. Harrington & Harrington (2013) evaluated the efficiency of the Sediment Rating Curve (SRC) method in modeling the suspended sediment load of the Bandon and Owenabue rivers in Ireland. Rajaee et al. (2009) and Chen & Chau (2016) indicated that the sediment rating curve and the auto-regressive integrated moving average model are inadequate to predict suspended sediment concentration under extreme hyperconcentrated flow conditions. Although the mentioned models led to promising results in sediment transport prediction, however, due to the importance of sediment transport and its impact on hydraulic structures, it is necessary to use other methods with higher efficiency (Morianou et al. 2017).

In recent years, the intelligence techniques such as Artificial Neural Networks (ANNs), Neuro-Fuzzy (NF), Genetic Programming (GP), Multivariate Adaptive Regression Splines (MARS), Kernel Extreme Learning Machine (KELM), and Gaussian Process Regression (GPR) have been used in assessing the complex hydraulic and hydrological phenomena (Roushangar & Ghasempour 2018) such as estimation of reference evapotranspiration (Yin et al. 2017), daily suspended sediment concentration modeling (Kaveh et al. 2017), bedload discharge modeling (Roushangar & Ghasempour 2018), side weir discharge coefficient modeling (Azamathulla et al. 2017), relative energy dissipation prediction (Saghebian 2019), roughness coefficient prediction in sever pipes (Roushangar et al. 2020), and modeling form resistance coefficient of movable bed channels (Saghebian et al. 2020). In general, the task of a machine learning algorithm can be described as follows: given a set of input variables and the associated output variable(s), the objective is learning a functional relationship for the input-output variables set. It should be noted that artificial intelligence models typically do not really represent the physics of a modeled process; they are just devices used to capture relationships between the relevant input and output variables. However, when the interrelationships among the relevant variables are poorly understood, finding the size and shape of the ultimate solution is difficult, and conventional mathematical analysis methods do not (or cannot) provide analytical solutions; these methods can predict the variable of interest with more accuracy.

On the other hand, hybrid models based on signal decomposition can be effective in increasing the time series prediction methods efficiency (Pachori et al. 2015). Wavelet analysis is one of the most common approaches used by researchers for this aim. Additionally, the Empirical Mode Decomposition (EMD) method, which is suitable for nonlinear and non-stationary time series (Huang et al. 1998), has been used recently. Unlike wavelet decomposition, empirical mode decomposition extracts the data oscillatory mode components without a priori determining the basis functions or level of decomposition (Labate et al. 2013). Also, according to Partalas et al. (2008) and Cloke & Pappenberger (2009) there is no unique model that is superior to others in all cases and the performances of different models may be different according to the condition of each intended parameter. Therefore, it is verified that the combination of outputs (from different models) through an ensemble method may lead to more accurate results. Integrating different models using an ensemble model as a post-processing method will represent different aspects of the underlying patterns more accurately (Zhang 2003).

Assessment of suspended sediment discharge is not a trivial matter, due to the multitude of factors influencing this parameter (e.g. flow and sediment conditions, bed material, bedforms, prophase sediment concentration, etc.). Due to the complex nature of sediment transport, the main aim of this study was temporal and special investigation of the SSD for successive hydrometric stations. Therefore, the Kernel Extreme Learning Machine (KELM) and Gaussian Process Regression (GPR) were used, as kernel based approaches. Kernel based approaches, based on quadratic optimization of convex function, can easily switch from linear to nonlinear separation. This is realized by nonlinear mapping using so-called kernel functions. Kernel based approaches such as KELM and GPR are a relatively new important method based on the different kernel types. Such models are based on statistical learning theory and are capable of adapting themselves to predict any variable of interest via sufficient inputs. Their training is fast, has high accuracy, and the probability of occurrence of data overtraining in this method is less. Also, Discrete Wavelet Transform (DWT) and Ensemble Empirical Mode Decomposition (EEMD) were used as pre-processing methods to improve the models' efficiency. In this regard, temporal features of the time series were decomposed via DWT and further breaking down was done via EEMD to obtain features with higher stationary properties. The time series decomposition to various periodicity scales with WT and further decomposition of sub-series using EEMD can lead to more stationary sub-series. Then, the sub-series energy values were calculated and the sub-series with higher energy were used as inputs in kernel based approaches to produce a multi-scale model for modeling the SSD time series. In the next step, the impacts of linear and nonlinear ensemble methods as post-processing approaches for improving the models' efficiency were assessed. For this aim, Simple Linear Averaging (SLAM) model was used as a linear ensemble method and Nonlinear Kernel Extreme Learning Machine Ensemble (NKELME) model was used as a nonlinear ensemble method. Between the GPR and KELM models, the KELM method was used for nonlinear ensembling. For evaluating the applied methods, daily sediment and flow data of the Mississippi River in the period of 2005–2008 were used and under two states various models were developed, in which at first, the intended parameter of each station was estimated using the station's own data, and then, the SSD was estimated using the previous stations' data.

The Mississippi River is the second longest river and chief river of the second-largest drainage system on the North American continent. From its traditional source of Lake Itasca in northern Minnesota, it flows generally south for 3,730 km to the Mississippi River Delta in the Gulf of Mexico. In the current study, daily data of the Mississippi River during the period of 2005–2008 were used. The Mississippi River basin has been repeatedly impacted by glacial ice during the last 2.5–3.0 million years. The continental glaciers eroded and deposited great masses of sediment and in many instances they were responsible for major adjustments in the direction and patterns of river drainage. Data on daily sediment load were available for areas in the USGS site where sufficient information on the sediment concentration and water discharge were obtained. The database on sediment records includes daily streamflow, daily mean concentration of suspended sediment, and daily suspended-sediment discharge. Three consecutive stations, namely station 1 (7,010,000), station 2 (7,020,500), and station 3 (7,022,000) were selected and suspended sediment discharge was investigated. During modeling, data from previous days at a determined point were used to predict the current suspended sediment discharge values and compare them with observed data. Although information about any changes in the study area (i.e. land use changes, erosion or deposition of the river bed, sediment conditions, etc.) was not available; however, the impact of these factors was reflected on the datasets and they had impact on modeling process. Table 1 shows the statistical characteristics of the selected stations. In this table, parameters Cs, Qs, and Wd are suspended sediment concentration, suspended sediment discharge, and flow discharge, respectively. Also, Figure 1 shows the location of the selected stations.

One of the most popular approaches in time series processing is the Wavelet Transform (WT) (Farajzadeh & Alizadeh 2017). The WT uses the flexible window function (mother wavelet) in signal processing. The flexible window function can be changed over time according to the signal shape and compactness (Mehr et al. 2013). After using WT, the signal will decompose into two approximation (large-scale or low-frequency component) and detailed (small-scale component) components. An illustration of a 3-levels WT is shown in in Figure 2. In the first level, the original signal (x) is decomposed to two components of approximation (A1) and detailed (D1). In the second level, A1 is again decomposed to approximation (A2) and detailed (D2) components. Finally, in the third level, A2 is decomposed to A3 approximation and D3 detailed components. The sum of all detailed sub-series and approximation series obtained from the third level will be the original signal (i.e. X = D1 + D2 + D3 + A3). The other approach for time series processing is Empirical Mode Decomposition (EMD). The EMD method is an effective self-adaptive dyadic filter bank which is applied to the white noise (a random signal which has equal intensity at different frequencies). By applying this method, each signal can be decomposed into a number of Inherent Mode Functions (IMFs), which can be used to process nonlinear and non-stationary signals. One of the advantages of this method is the ability to determine the instantaneous frequency of the signal. At each step of the signal decomposition into its frequency components, the high frequency components are separated first and this process must continue until the component with the lowest frequency remains (see Lei et al. 2009 for more details). EEMD is developed based on EMD. The main benefit of EEMD is solving the mode mixing problem of EMD, which determines the true IMF as the mean of an ensemble of trials (Wu & Huang 2009). The EEMD algorithm can be described as: (1) for a given signal x(t), random white noise is added to the signal, (2) the noise-added signal is decomposed using EMD for obtaining IMF series, (3) steps 1 and 2 are repeated until the number of added white noises is greater than or equal to the number of trials, (4) for obtaining the ensemble IMF, the average of the sum of all IMFs is computed (I_j(t)), and (5) the original signal is formed as ⁠. For selecting the most effective IMFs and using them as inputs in the modeling process, their energy values can be calculated and the IMFs with higher energy can be used as inputs.

Kernel based (KB) approaches are new methods that are used for classification and regression purposes (Roushangar et al. 2019). Two important kernel based approaches are Gaussian Process Regression (GPR) and Kernel Extreme Learning Machine (KELM), which work based on the different kernel types such as Linear, Polynomial, Radial Basis Function (RBF), and Sigmoid functions in SVM and Exponential, Squared Exponential, Rational Quadratic, ARD Rational Quadratic, and so on in GPR. Kernel based approaches are based on statistical learning theory initiated and can be used for modeling the complex and nonlinear phenomenon. The kernel type affects the training and classification precision; therefore, the most important step in kernel based approaches is the appropriate selection of kernel type. These methods are memory intensive, and trickier to tune due to the importance of picking the right kernel (Babovic 2009). The aim of the KB approaches is to determine a function that has the most deviation from the actual target vectors for all given training data (Roushangar et al. 2020).

The matrix T is identified as the target matrix of the hidden layers of the neural network. H is considered as the output matrix of the neural network. Huang et al. (2012) also introduced kernel functions in the design of ELM. Presently, a number of kernel functions are used in the design of ELM such as linear, radial basis, normalized polynomial, and polynomial kernel functions. Kernel function based ELM design is called Kernel Extreme Learning Machine (KELM). For more detail about KELM, readers and researchers are referred to Huang et al. (2012).

Also, deriving from the Bayesian framework, GPR can be considered as a random process to carry out nonparametric regression with the Gaussian process (Rasmussen & Nickisch 2010). In most prediction issues, GPR is preferred owing to its flexibility to provide the uncertainty representations.

It should be noted that each artificial intelligence method has its own parameters for achieving the desired results, and the optimized amount of these parameters should be determined. For example, in designing the KELM and GPR approaches, selection of the appropriate type of kernel function is needed. There are various kernel functions that can be used based on the nature of the studied phenomenon. In this research, in order to determine the best performance of the GPR and KELM models and selection of the best-kernel function, the model (III) was tested for the SSD prediction in station 1 via various kernel functions. Figure 4 indicates the results of statistical parameters of different kernel functions for this model. According to Figure 4, the RBF kernel function in the KELM model and squared exponential in the GPR model were found to be the best kernel function. Therefore, the RBF and squared exponential kernels were selected as the core tool of the KELM and GPR, which were applied for the rest of the models.

Accurate prediction of the suspended sediment discharge in rivers or streams is very crucial for sustainable water resources and environmental systems. In this study, the suspended sediment discharge in the selected stations was assessed via the KELM and GPR kernel based approaches. The previous parameters of flow and suspended sediment were used for model development. The results obtained from the developed GPR and KELM models are listed in Table 3 and shown in Figure 5. The results indicated that in the state of modeling based on each station's own data, the model (III) with input parameters of C_S(t-1), C_S(t-2), Q_s(t-1), and W_d(t-1) performed more successfully than the other models. Kisi et al. (2012) showed that in the suspended sediment modeling via genetic programming, the model whose inputs were the current water discharge and one previous water discharge and sediment load performed better. Also, the models whose inputs were the current and one immediate previously recorded water discharge and one and two previous sediment loads, as well as models whose input were the current water discharge and one previous sediment load, were ranked as the second and third best models, respectively. In this study, from the obtained results, it could be stated that the model (II) with two input parameters of Qs(t-1), Qs(t-2) yielded the desired accuracy. Therefore, the sediment discharge could be predicted using only the previous one and two days' sediment discharge variables. In the second state, it could be seen that between the selected stations, the SSD modeling of the third station based on data from both first and second stations performed more successfully. A comparison between the results of two considered states showed that modeling based on the each station's own data led to more desirable results. However, using the previous stations' data in modeling process yielded relatively accurate results; therefore, via the kernel based approaches, the previous stations' data could be used when the station's own data were unavailable. Artificial intelligence methods are very powerful tools and when the interrelationships among the relevant variables are difficult to understand and conventional mathematical analysis methods cannot provide analytical solutions, these methods can be used successfully. Choubin et al. (2018) evaluated the use of a Classification and Regression Tree (CART) model to estimate the SSD based on hydro-meteorological data. They indicated that the CART as an artificial intelligence model can be a helpful tool in basins where hydro-meteorological data are readily available. In this study, the applied methods had desirable efficiency in the SSD modeling under both considered states. It should be noted that, in the second state, an attempt was made to investigate the issue that the existing sub-basins between the consecutive stations may have noticeable impacts on the flow regime of the downstream station. For this state, if there be special conditions between the stations (such as diversion dams, intake structures, tributary, etc.) the relationship between the flow regime of the upstream and downstream may be affected. In this study, there was a tributary between stations 2 and 3. However, since artificial intelligence methods are black-box methods and do not show the physics of the phenomenon and, also, the effect of the tributary is reflected on all data; therefore, by comparing the obtained results for two states of the SSD modeling (i.e. the SSD modeling of station 2 using the station 1 data with the SSD modeling of station 3 using the station 2 data), it seems that the effect of the tributary on the modeling process cannot be accurately expressed.

In the next step, for evaluating the impact of the independent parameters, sensitivity analysis was performed for the best model of each state. In this regard, the input variables of the best models were omitted one by one and the KELM model, which had higher efficiency compared with the GPR model, was rerun. Based on the results from Figure 5(b), it could be induced that in the state 1 parameter Qs(t-1) and in state 2 the parameter Q_s1(t) for station 2 and the parameter Q_{s1, 2}(t-1) for station 3 had the most significant impact on the SSD prediction.

In this section, the effect of time series pre-processing on increasing the models' accuracy was investigated. Therefore, the time series were decomposed via the WT and EEMD methods. In time series decomposition by WT, a mother wavelet that is more similar to the signal should be selected. In this study, the Daubechies (db2, db4) and Symlets (sym2, sym4) mother wavelets were trained for decomposition of Qs(t-1) in station 1. According to Figure 6(a), it was found that the db4 mother wavelet led to better outcomes. Therefore, the db4 mother wavelet was used for time series decomposition. Also, in the second step, data was decomposed via EEMD. The principle of EEMD is decomposition of the signal to different IMFs and one residual signal. The sum of these signals will be the same original signal. The formation of IMFs is based on subtracting the basic function from the original signal. This process continues until the residual signal remains almost constant. In this study, time series were decomposed into 10 IMFs and one residual signal. In this study, details 1 and 2 sub-series were decomposed into nine IMFs and one residual signal. Further decomposition using EEMD leads to sub-series with more stationarity and with less noise.

Since the number of input data increased after time series decomposition, the energy values of sub-series were calculated and the sub-series with higher energy were selected as inputs. Then, the selected sub-series were used as inputs in kernel based methods to predict the SSD. In Figure 6(b), the Qs(t-1) sub-series decomposed by db4-EEMD is shown for station 1. According to this figure, the IMFs 2 and 3 from detail 1 and the IMFs 4 and 5 from detail 2 were selected as appropriate IMFs based on their higher energy values in comparison with the other IMFs.

Then, the selected sub-series were used as inputs in the KELM and GPR models to predict the SSD. The results of the integrated pre-processing models are listed in Table 4. According to the results presented in Tables 3 and 4, it could be induced that data pre-processing significantly improved the results' accuracy and integrated models were more accurate than single methods. In fact, the use of WT and further decomposition of the detailed series led to an improvement in the outcomes. For example, for station 3 in state 1, the RMSE values of the KELM and GPR models for testing sets of the model (III) were 0.048 and 0.05, respectively; the values of RMSE for the WT-EEMD-KELM and WT-EEMD-GPR models decreased to 0.028 and 0.033, respectively. In general, the data pre-processing increased the models' accuracy by between 10 to 15% in calibration sets and 9 to 20% in verification sets. Hazarika et al. (2020) used machine learning methods with wavelet conjunction and indicated that the wavelet transformations model helps to analyze the time and frequency information to estimate sediment load by decomposing data over several phases. Their study showed that the hybrid models based on the Coiflet wavelet offered good performance. Also, Nourani et al. (2019) used the hybrid Wavelet-M5 model to model suspended sediment discharge of two different rivers. The obtained results showed the better performance of Wavelet-M5 model in comparison with individual ANN and M5 models. In this study, from the results, it could be seen that the integrated methods resulted in desirable accuracy in both states and the use of previous stations' data could be utilized as a reliable scenario to predict the SSD values for the stations with a lack of observational data.

In the next part, an attempt was made to combine the results of pre-processing-kernel based methods and evaluate the simultaneous impact of pre-post-processing on the accuracy of the outputs. The obtained results are listed in Table 5 and shown in Figure 7. According to the results, it could be indicated that the combined methods yielded more accurate results regarding the single and pre-processing-kernel based approaches. It could be seen that both simple and nonlinear ensemble methods caused an increment in the applied models' efficiency, respectively by about 5% to 10% for verification sets and 5% to 15% for testing sets in comparison with pre-processing models. Based on the results listed in Tables 4 and 5, it seems that in enhancing the predictions' accuracy, the pre-processing methods performed more successfully than the post-processing methods. It should be noted that determining the extreme values of the SSD is the most important issue in river engineering. Therefore, in the use of different models, their performance in estimating minimum and maximum values of time series should be taken into account. According to Figure 7, it can be seen that the extreme values in the SSD series calculated more correctly by the pre-post-processing methods. The results showed that in the SSD series prediction, the nonlinear ensemble method performed better than the linear method. This issue stated that the efficiency of the linear method relates to the performance of individual models. Therefore, for artificial intelligence models with weak performance, the ensemble model results would not be desirable. However, using a nonlinear pattern for simulation would lead to more appropriate results.

The RMSE error criterion was used to graphically compare the performance of the single and integrated methods. The results are shown in Figure 8. As can be seen, in the SSD modeling process, the RMSE values were smaller for the integrated methods and the WT-EEMD- NKELME model led to more accurate results.

The accurate prediction of the SSD values of rivers is an important factor in improving water management. This study assessed the capability of time series pre-processing methods for SSD modeling. In this regard, in the first step, the time series without any processing were imposed to single kernel based methods. Then, time series were decomposed to several sub-series using WT and further decomposition was performed via EEMD. Finally, the obtained results from pre-processing methods were combined via two linear and nonlinear post-processing methods in order to enhance the forecasting efficiency. According to the results, it was found that using integrated methods increased the models' accuracy. The applied pre-processing method enhanced the KELM and GPR models' performance approximately between 10 and 20%. It showed that the SSD modeling based on each station's own data led to more desirable results. In this state, the model with inputs C_S(t-1), C_S(t-2), Q_s(t-1), and W_d(t-1) was the superior model. However, using the integrated approaches, the previous station's data can be used when the stations' own data were unavailable. The obtained results revealed that the pre-processing methods performed more successfully than the SLAM and NKELME post-processing methods. The simple and nonlinear averaging techniques improved the efficiency of the applied models by about 5% to 15% compared with pre-processing models and the performance of the NKELME was higher than the SLAM. Also, it was found that the maximum and minimum values of the SSD variables were well predicted using the integrated models. Therefore, the integration of hybrid pre-post-processing techniques could be useful for the daily SSD modeling.

All relevant data are available from an online repository or repositories at https://waterdata.usgs.gov/nwis/sw.