Improving Accuracy of River Flow Forecasting Using LSSVR with Gravitational Search Algorithm

River flow prediction is essential in many applications of water resources planning and management. In this paper, the accuracy of multivariate adaptive regression splines (MARS), model 5 regression tree (M5RT), and conventional multiple linear regression (CMLR) is compared with a hybrid least square support vector regression-gravitational search algorithm (HLGSA) in predicting monthly river flows. In the first part of the study, all three regression methods were compared with each other in predicting river flows of each basin. It was found that the HLGSA method performed better than the MARS, M5RT, and CMLR in river flow prediction.The effect of log transformation on prediction accuracy of the regressionmethods was also examined in the second part of the study. Log transformation of the river flow data significantly increased the prediction accuracy of all regression methods. It was also found that log HLGSA (LHLSGA) performed better than the other regression methods. In the third part of the study, the accuracy of the LHLGSA and HLGSAmethods was examined in river flow estimation using nearby river flow data. On the basis of results of all applications, it was found that LHLGSA and HLGSA could be successfully used in prediction and estimation of river flow.


Introduction
River flow forecasting plays a vital role in planning of water projects, irrigation systems, hydropower system, and optimized utilization of water resources [1].Due to continuous increase of population growth, industrial uses, and irrigation needs, the river flow forecasting has received great attentions of researchers for operational river management [2].Forecasting of river flow provides alerts of approaching floods and also assists in controlling the outflows of reservoir during low flows days of river.Floods affect countless lives, infrastructure, and property and cause limitless damage more than any other natural disaster.Due to no assessment of flood magnitude, a flood resulted in a loss of thousand lives and damage of agriculture land of million dollars in Pakistan in 2010 [3].It is not possible to provide complete safety from flood, but high amounts of money and many lives can be saved by providing accurate flood predictions, flood magnitude, and flood duration [4].The importance of water measurement compelled researchers to apply various types of forecasting methods to estimate and forecast river flows.
From the last three decades of the previous century, the statistical methods were applied successfully in the field of hydrology including the river flow forecasting.Statistical methods try to find inherent relationships within the actual data.The autoregressive integrated moving average (ARIMA) and seasonal autoregressive integrated moving average (SARIMA) methods are the most popular in the statistical methods category and have been extensively used to model different variables in the field of hydrology [5][6][7][8][9][10][11][12][13][14][15].Ahlert and Mehta [5] and Kurunc ¸et al. [7] used ARIMA statistical models for modeling river flows data.Ahmad et al. [6] and Mirzavand and Ghazavi [8] applied ARIMA statistical methods to analyse water quality and groundwater data, 2 Advances in Meteorology respectively.Otok and Suhartono [11], Rabenja et al. [13], and Valipour [14] forecasted runoff data in Indonesia and USA, respectively, by applying SARIMA model and compared with ARIMA.Psilovikos and Elhag [12] and Yang et al. [15] applied ARIMA models successfully in modeling different processes of evaporation data.Mishra and Desai [9] and Modarres [10] used SARIMA method efficiently for drought forecasting in India and Iran, respectively.In the previous two decades, the artificial neural networks (ANN) have been replaced with the statistical methods in solving different problems due to their flexible nature and capturing the nonlinearity in the data.In the literature, many researchers compared the ANN with statistical methods in solving many problems of hydrology and reported that ANN outperformed the statistical methods [16][17][18][19][20][21].Huang et al. [18] used the ANN method to forecast the river flows of Apalachicola River, USA, using the previous rainfall and river flow data.They compared the quarterly and yearly river flow forecasts results with the ARIMA method's results and found that the ANN performed better than the ARIMA method in prediction of river flow.The detailed discussion of all ANN applications in comparison with statistical methods to forecast different variables in hydrology is not possible in this paper.However, ANN also have some major weakness, that is, overfitting, falling into local minimum, slowing convergence speed, and requiring large number of training data.Thus, in the last decade, the support vector regression (SVR) took priority over ANN due to its parallel distributed processing, self-learning features, avoiding the overfitting problems, and providing globally optimal solutions [22][23][24][25][26][27][28].Ahmad et al. [22] applied SVR model to forecast runoff of Bakhtiyari Basin, Iran, and results explored that SVR showed better accuracy than the ANN methods for daily runoff forecasting especially in case of prediction of higher values of river flows.However, SVM faces computationally difficulties in determining optimal solution due to use of quadratic programming with nonlinear equation.This procedure is time consuming.
Recently least square support vector regression (LSSVR), the improved version of SVR, received much attention in the field of prediction methods due to use of linear squares principle for the loss function instead of the quadratic programming in the SVR method and fast computational speed [29][30][31][32][33][34][35][36][37][38].Shabri and Suhartono [37] and Kisi [34] applied LSSVR successfully to forecast river flows.Shabri and Suhartono [37] compared the prediction accuracy of LSSVR with ANN and multivariate linear regression methods whereas Kisi [34] compared it with adaptive neurofuzzy embedded fuzzy -means clustering (ANFIS-FCM) method and they both found that LSSVR performed better than the other methods.Kisi [33] and Goyal et al. [30] forecasted the reference evapotranspiration and pan evaporation by using LSSVR method.Kisi [33] compared its results with feed forward ANN whereas Goyal et al. [30] compared it with ANN and ANFIS methods and they found that LSSVR performed better than the other methods.Okkan and Serbes [36] and Bhagwat and Maity [29] successfully used LSSVR method to forecast runoff data by using the previous meteorological and river flows data.Kisi [32] estimated the suspended sediment by using the river flow data through LSSVR method and reported that the LSSVR gave better estimates in comparison with ANN and sediment rating curve (SRC) methods.Hwang et al. [31] predicted the daily water demand of the Seoul City, Korea, and daily mean inflow of Cheng-ju Dam by using LSSVR model.He compared the predicted results of LSSVR with the conventional multiple linear regression (CMLR) and back propagation neural network methods in both cases and found that LSSVR showed superiority in prediction accuracies.Wu et al. [38] and Mellit et al. [35] applied LSSVR method to predict the different meteorological variables and found that LSSVR performed better than ANN method.Motivated by these successful applications of the LSSVR, the LSSVR method was selected as a forecasting method in this research.LSSVR method has parameters which play a vital role in determining the prediction accuracy of the method.Determining suitable value of these parameters will produce better river flow prediction results.Still, there is no specific way to determine optimal parameters for LSSVR method in the literature of river flow forecasting.Thus, the novelty of this study is to generate a hybrid LSSVR-gravitational search algorithm (HLGSA) river flow forecasting method.In this study, GSA is used to find the optimal values of LSSVR method to increase the prediction accuracy of the method.GSA was preferred in this study over other heuristic algorithms such as simulated annealing algorithm, genetic algorithm, memetic algorithm, differential evolution, and particle swarm optimization due to their premature convergence, parameter sensitivity, and consuming too much time to obtain global optimal solution.Instead of these heuristic algorithms, GSA improves the global search ability and optimization speed by using the principle of gravity and motion.To the best knowledge of the authors, there is not any published work in the literature that predicts the river flow using hybrid LSSVR-gravitational search algorithm (HLGSA) method.Recently, researchers preferred hybrid methods for solving different problems in the field of hydrology.
In addition to the HLGSA method, multivariate adaptive regression splines (MARS) is another popular regression method used to model the complex nonlinear relationships among the variables.MARS is a nonparametric regression method and it has been applied extensively nowadays in the field of hydrology to predict different variables [39][40][41][42][43][44][45].To determine the benefits of using MARS over other conventional regression methods, MARS method was compared with CMLR and model 5 regression tree (M5RT) methods in this study.Cross validation (CV) technique was used to better see the prediction accuracy of all applied methods.Log transform function was also utilized in this study to see its effect on the prediction accuracy of these methods.

Hybrid LSSVR-Gravitation Search
Algorithm (HLGSA) Method for River Flow Prediction  by using linear equation instead of quadratic equations [47].Figure 1 demonstrates the process of LSSVR algorithm.By using time series inputs   (lagged river flows) and output   (predicted river flow), the function of nonlinear LSSVR is given as where ⟨⋅, ⋅⟩ represents dot product, () is a nonlinear function that employs regression, and  and  are the weight vector and bias term, respectively [48].The cost function () of LSSVR can be minimized as Subject to   = ⟨,  ()⟩ +  +   ( = 1, 2, . . ., ) , where ,   represent the regularization constant and the training error for   , respectively.To solve (2), the Lagrange function is used to find the solutions of  and .The Lagrange function can be calculated as [49] where   ∈  is the Lagrange multipliers.The solution of above equation can be achieved by determining the partial differential of Lagrange function and applying the kernel function (KF) to satisfy Mercer's condition.To solve regression problems, there are many types of KF including polynomial, radial basis, Gaussian, sigmoid, Mexican hat, Meyer, and Morlet.The KF type plays a vital role in constructing high accurate LSSVR model [50].This study used the radial basis KF (RBKF) due to its effectiveness for the nonlinear regression problems [51].The performance of the RBKF with other KFs is shown in Section 7. RBKF can be expressed as After selecting the RBKF for the LSSVR method, finding proper values for penalty factor parameter that is  and RBKF parameter that is  2 is necessary.There is no specific way to obtain the optimal values of parameters.Due to these reasons, GSA is adopted in the study to calculate the suitable parameter values.

Gravitational Search Algorithm (GSA).
GSA is one of the effective optimization algorithms compared with other evolutionary algorithms.It is based on the law of gravity and motion and first proposed by Rashedi et al. [52].In GSA, each agent has four parameters: position, velocity, inertial mass, and gravitational mass.The location of the agent corresponds to a solution of the problem whereas its gravitational and inertia masses are obtained utilizing a fitness function [53].
The location of particle can be expressed as where    represents the location of the th agent in the kth dimension.The mass of each agent is computed after calculating the fitness of current population as [52,54] where The total gravitational acceleration of the th agent can be calculated using the law of motion as follows: where    () represents the gravitational acceleration of the agent  in the kth dimension and rand indicates a random variable with uniform distribution in the interval [0, 1].With the help of (9), the total gravitational force exerted on the agent  in the kth dimension can be calculated as Then the speed and location of the agent are updated as follows: It is clearly seen from the brief description of the GSA that it utilizes the gravitational force as the direct form to communicate the agents' cooperation.The heavy agents in GSA are processed, infer good solutions, and move more gradually than lighter ones, which guarantee the algorithm's exploitation step.In other words, the GSA searches for the ideal solution by appropriately calibrating the inertia and gravitational masses of agents where every agent provides a solution.As time progresses, the heaviest agent will exhibit an ideal solution in the search space [55].

HLGSA (Hybrid LSSVR-GSA).
The process of constructing the river flow prediction model HLGSA by using the hybrid of LSSVR and GSA methods is described in this section and shown in Figure 2. The process is as follows: (i) Firstly, divide all river flow data sets into training and test parts.
(ii) Select the RBF kernel function and initial parameters for the HLGSA method to build the initial LSSVR model.The initial value of the parameters is set as follows: the range of penalty factor  is 0.1 to 2000, the range of RBF parameter  2 is 0.001 to 20, number of iterations is 15, the number of particles can be set up to 40, and constant alpha is found to be better in range of 16 to 20, whereas initial gravitational constant  0 is found to be better in range from 105 to 115.
(iii) Compute the particle fitness value of each agent.In this paper, RF RMSE is selected as the fitness function.
The fitness function for this method can be defined as (iv) Choose the best parameters combination through GSA to obtain the optimal values of the LSSVR parameters.
(v) If it does not meet the stopping criterion, then utilize the new combination of parameters to reconstruct the LSSVR.Compute the fitness until it suits the stopping criterion.
(vi) The ideal parameter values are achieved to build the optimal LSSVR model for forecasting river flow.Now, the testing values are used for the optimal LSSVR model to get river flow prediction results.

Regression Methods
In this study, the performance accuracy of a hybrid nonlinear optimized regression method (HLGSA) was compared with a nonlinear, nonparametric regression method (MARS), with a piecewise linear regression method (M5RT), and with a conventional linear regression method (CMLR) in forecasting monthly river flow.

Multivariate Adaptive Regression
Splines.MARS is a flexible method which finds relationships that are nearly additive or involve interactions with fewer parameters.The general MARS method is introduced by Friedman [56] and is expressed by the following equation: where RF  is the forecasted river flow by the MARS that is dependent variable,  0 is a constant,   are the model coefficients calibrated to provide the best fit to the used data,  is the quantity of basis functions (BFs),  is the "splits" quantity that generates the mth BFs, and   gets values of 1 or −1 and represents the (right/left) sense of the associated step function. (,) is the independent variable's label [57].Two-step MARS provides optimal MARS model.MARS develops a huge number of BFs chosen to overfit the data at first step, where variables are permitted to enter-as continuous, categorical, or ordinal-the formal system by

Obtain the optimal parameters
Build the optimal LSSVR model

Yes
No which variable ranges are characterized, and they can interact with each other or be restricted to enter only as additive components.In the second step, BFs are erased in the order of minimum effect utilizing the generalized cross validation criterion (GCV).A measure of variable significance can then be evaluated by watching the decrement in the computed GCV when a parameter is excluded from the model.This procedure proceeds until the rest of the BFs all satisfy the predecided necessities.The GCV can be computed as [58]: where RF  and RF  are the actual and predicted river flow values and () is a complexity penalty function.
After building the MARS model, the relative importance of a variable in terms of its contribution to the fit of the model can be estimated.MARS is capable of tracking very complex data structures, so selected in this study for modeling river flow time series.

Model 5 Regression Regression Tree (M5RT).
In decision tree (DT), each branch node indicates a choice between a number of alternatives and a decision is made in every leaf node [59].Regression trees (RT) are applied to solve those forecasting problems having numeric response variable.They are different from the DT only in that they involve a numeric value rather than a class label combined with the leaves [60].The M5RT method combines the features of DT and RT methods because the construction of the M5RT is similar to the DT but, instead of the class labels, it has linear regression functions at the leaves.The M5RT is a piecewise linear method that was introduced by Quinlan [61] and has many successful applications in the field of water resources [41,[62][63][64][65][66][67][68] that compelled the authors to use M5RT method in this paper for river flow prediction.
The division criteria for the M5RT method are based on reducing the standard deviation of the class values that reach a node as an error measure and computing the estimated reduction in this error as a consequence of testing each attribute at that node.The standard deviation reduction (SDR) is computed by where  stands for set of samples that enters the node,   indicates the subset of samples that have the th output of the potential set, and sd is the standard deviation [69].

Conventional Multiple Linear Regression (CMLR).
The multiple linear regression methods forecast values of a dependent variable  based on independent variables ( 1 ,  2 , . . .,   ).Two main advantages of the CMLR are that it has simple structure and it is included in lots of statistical packages [70].In this study, after determining the independent lagged river flow values for dependent river flows of both basins, the CMLR can be constructed as follows: where RF () is the dependent variable,   − −  are the equation parameters for the linear relation, and RF (−1) , RF () are the independent lagged river flow value used to forecast river flow.However, CMLRs have some disadvantages in predicting nonlinear situations because of their linear structure [71].

Study Sites and Data Preprocessing
The study used the river flows data from two catchments, Astore and Shyok, on the Upper Indus Basin of Pakistan.
Figure 3 shows the location map of the catchments.In this research, cross validation (CV) technique was applied to better see the prediction accuracy of the applied methods.In CV technique, the whole data is divided into  equal data sets (DS), then the −1 DS is used to train, and the other one DS is used to test the accuracy of the method.This process is repeated  times till every DS of the data is used to test the applied method.This CV technique is preferred over -fold cross validation techniques due to usage of every data set for testing and  >  which makes it closer to the real world problem [69,72].Similarly, in this research, the whole river flow data was divided into four equal DS.In all the applications, the three DS were used to train and remaining one DS was adopted to test the method.This procedure was repeated four times till every DS of data was used to test the method.The monthly river flow time series statistics of Astore and Shyok catchments are reported in Table 1.
Here, the DS1, DS2, DS3, and DS4 represent four equal data sets of whole data for CV analysis whereas RF mean , RF sd , RF sc , RF min , and RF max represent mean, standard deviation, skewness coefficient, minimum, and maximum river flows, respectively.The recorded monthly river flows data show similarly high positive skewed distribution for Astore and Shyok catchments (RF sc = 1.51 and 1.88).However, the range of the flow data of Shyok Basin (36.7-2080.7 m 3 /s) is much higher than that of the Astore Basin (19.3-654.9m 3 /s).The lagged values of river flows show low persistence (e.g., Lag 1 = 0.735, Lag 2 = 0.266, and Lag 3 = −0.152).However, the lagged values of Astore Basin give a little better persistence than those of the Shyok Basin.are generally taken as IC in many researches of river flow forecasting [37,73,74].In this research, the effect of lagged river flows values of both catchments was examined through autocorrelation function and was reported in Figure 4.

Input Combination Selection and Performance Evaluation Criteria
According to the analysis, the following three IC were selected as inputs on the basis of most significant lagged river flows values for both basins, that is, In this paper, two error indices were selected to evaluate the performance of the models in prediction of monthly river flows including the root mean square error (RF RMSE ) and mean absolute error (RF MAE ).The similarity between the observed value and the forecasted value of river flow is measured by using the determination coefficient (RF DC ) index.These three indexes have been extensively applied in many problems of water resources for evaluating the model performance [19,26,75].All the performance evaluation indexes can be calculated as where  RF is the total size of observations of river flow time series, RF  is observed river flow, RF  is forecasted river flow, RF  is average of river flows, and RF  is average forecasted river flow.

River Flow Prediction Using Soft Computing Methods
In the first part of the paper, the performance of the proposed method hybrid LSSVR-GSA (HLGSA) was compared with other regression methods in predicting monthly river flows   For the sake of simplicity, the performance accuracy of all methods was evaluated by comparing the overall mean errors representing the mean error of all data sets and input combinations.Figures 5(a)-5(b) show the overall average errors statistics of all methods.Mean errors statistics of RF RMSE and RF MAE clearly indicate that the HLGSA method performs better by having relatively less value of error indexes than the other methods in prediction of Astore catchment's river flows.M5RT and CMLR give almost same mean errors statistics and provide the worst accuracy prediction in comparison to HLGSA and MARS due to having higher values of both errors indexes.This indicates the nonlinearity of the investigated phenomenon because both M5RT and CMLR have linear structures.HLGSA decreases the overall mean RF RMSE of the MARS, M5RT, and CMLR by 8.22%, 23.15%, and 24.49%, respectively.The observed and forecasted monthly river flows of Astore Basin by all the methods using their best model structures are reported in Figures 6(a)-6(d).The figure clearly explores that the HLGSA method is in good fit with the original river flows data in comparison to the other methods.The HLGSA gives higher value of correlation index (RF DC ) than the other methods.From Figure 6, it can be clearly seen that the RF DC value of the HLGSA method is 0.916, which is higher than RF DC of MARS, M5RT, and CMLR methods (which are 0.900, 0.888, and 0.882).
Table 3 reports the results of three performance evaluation statistical indexes for the HLGSA, MARS, M5RT, and CMLR methods in forecasting river flows of the Shyok catchment.In case of input combinations, here also IC1 gives the worst forecast results compared to the IC1 and IC3.However, in contrast to Astore application, here IC3 shows better accuracy than the IC2 for all the applied methods.In the case of data sets, it is evident from the table that DS4 gives the worst prediction results for all the methods.Similar to Astore Basin, here also the maximum river flow value of the Shyok Basin's test set, DS4 (RF max = 2080.7 m 3 /s), is higher than that of the training value (see Table 1).Another reason of this may be the fact that the DS2 has low correlations with the preceding river flow input data in comparison to the other data sets.It is obvious from the table that all methods give good prediction results for the DS3 among all input combination scenarios.The best model structures for all four applied methods in case of Shyok catchment are found for the DS3 and IC3.Similar to the Astore Basin, here also Table 3 clearly explores that the HLGSA method outperforms the other methods from RF RMSE , RF MAE , and RF DC viewpoints for all data sets and input combination scenarios.Here, also the MARS gives better accuracy than the M5RT and CMLR methods for all data sets and input combinations.However in contrast to Astore Basin, here M5RT provides better prediction results than the CMLR method.The best RF RMSE for the proposed method is better as 109.20 m 3 /s compared

Effect of Log Transform on the Prediction Accuracy of the Soft Computing Methods
In this section of the paper, the effect of log transform on the prediction accuracy of the applied methods was investigated   The mean errors statistics of all log models is illustrated in Figures 9(a According to comparison of scatter plots of log and normal methods (Figures 6 and 10), it is evident that the log methods have better fits with the original data in comparison to normal methods.On the basis of fit line equation, the log methods are closer to the exact line than the normal methods (see  0 and  1 coefficients in Figures 6 and 10).LHLGSA decreases the overall mean RF RMSE of the LMARS, LM5RT, and LCMLR by 11.46%, 24.16%, and 35.65%, respectively.
Test statistics of the LHLGSA, LMARS, LM5RT, and LCMLR methods for the Shyok Basin are reported in Table 6.Here, in contrast to previous Shyok application, IC2 gives better forecast results for the LHLGSA and M5RT methods whereas IC3 performs better for the LMARS and LCMLR methods.However, IC3 performs slightly better for the LM5RT, whereas IC2 performs better for the LMARS in the case of DS3 data set.Similar to the previous Shyok application, here also DS3 performs the best whereas DS4 performs worst among all data sets.To check the best log method among all log methods, the mean errors indexes (RF RMSE and RF MAE ) of all log models are plotted in Figures 11(a The figure clearly shows that the LHLGSA method gives less scatter estimates with a higher value of RF DC in comparison to the LMARS, LM5RT, and LCMLR methods.The figure also reveals that all log methods are closer to exact line than the normal methods (compare Figures 8 and 12).On the basis of comparison, it is obvious that the normal methods give more scattered forecasts than the log methods.In case of Shyok catchment, LHLGSA decreases the overall mean RF RMSE of the LMARS, LM5RT, and LCMLR by 14.49%, 21.73%, and 43.84%, respectively.The log transform does not equally affect all input combinations and data sets but it can be observed that this effect is more prominent when more inputs are used in the case of CMLR methods (see Tables 5  and 6).In contrast to Astore Basin application, for the Shyok Basin, the best RF RMSE and RF MAE values for the HLGSA (109.20 m 3 /s and 55.75 m 3 /s) are a little better compared to LHLGSA (109.41 m 3 /s and 59.01 m 3 /s), respectively.However in case of similarity index (RF DC ), the best RF DC value for LHLGSA (0.959) is better than that of the HLGSA (0.947) method.In case of LMARS, LM5RT, and CMLR method, the best RF RMSE , RF MAE , and RF DC values are better in comparison with MARS, M5RT, and CMLR methods, similar to Astore Basin application.
To evaluate the overall effect of log transform function on all applied methods, the overall mean error indexes (RF RMSE and RF MAE ) for both basins with and without log methods are compared in Figures 13(a), 13(b), 14(a), and 14(b).Both graphs clearly explore that the log methods give better accuracy than the normal methods except the LCMLR which gives worse accuracy compared to the CMLR with respect to RF RMSE error index.However, the LCMLR gives better accuracy than the CMLR from the RF MAE index viewpoint for both basins.The reason of the LCMLR's bad performance is due to its worse results in case of IC1 that affects the overall RF RMSE value.The bar graphs also prove that the proposed HLGSA method shows better accuracy than the other models in both cases of normal and logarithm transformed time series data.According to the bar graphs, the LHLGSA reduces the overall RF RMSE of the HLGSA by 5.66% and 4.87% for the Astore and Shyok Basins, respectively.The LMARS reduces the overall RF RMSE of the MARS by 2.20% and 1.52% for the Astore and Shyok Basins, respectively.The LM5RT reduced the overall RF RMSE of the M5RT by 4.41% and 2.22% for the Astore and Shyok Basins, respectively.However, in the case of multiple linear regression method, CMLR reduced the overall RF RMSE of the LCMLR by 9.05% and 7.45% for the Astore and Shyok Basins, respectively, while, in the case of RF MAE error index, the LCMLR reduced the overall RF MAE of the CMLR by 3.60% and 10.80%, respectively.In this last section of the research, the performance of the HLGSA and LHLGSA is evaluated in river flow estimation of a basin using flow data from a nearby basin.The river flows estimation using nearby basin's flow data is a vital issue because Pakistan is a developing country and many basins have long duration of missing flows data due to financial problems in the maintenance of the hydraulic gauging stations at higher altitudes.Since river flows play a key role in planning and designing of hydropower projects and for the flood mitigation, it is necessary to find a suitable way to fill these missing flows data.In this paper, the river flows data of the Astore Basin is used to estimate the flow data of the Shyok Basin.In this application, also the CV technique is applied to better see the accuracy of both methods in estimating river flows.

Discussion
On the basis of above results, the key findings can be summarized as follows.
(1) The Astore Basin reported lower values of RF RMSE and RF MAE in comparison with the Shyok Basin.This is due to the mean river flow of the basins; Astore Basin Advances in Meteorology is characterized by mean river flow of 142 m 3 /s while the mean river flow of Shyok Basin is 457 m 3 /s.
(2) It was also observed that the prediction accuracy of all methods including proposed method was mostly improved with increasing in input numbers which indicated that all input combinations have positive effects on predicting river flow especially in case of CMLR method.
(3) It was also found that the higher value of testing data set's maximum river flow in comparison with other training data set's maximum river flow values caused the extrapolation difficulties and produced worst prediction results for that data set.
(4) Overall, the HLGSA and LHLGSA methods outperformed the MARS, M5RT, CMLR and LMARS, LM5RT, and LCMLR methods, respectively.Moreover, the comparison between Figures 13 and 14 indicates that the prediction results with log transform function are better on the mean basis than the prediction results without log transform function using all regression methods including the proposed method which means the log transform function is suitable for denoising the river flow data.
(5) In literature, many studies reported that MARS performed better or equally in comparison with the LSSVR methods [43,44,[76][77][78].However, in this study, the hybrid LSSVR method with gravitational search algorithm performed better than MARS method.The main reason behind this may be the LSSVR's strong generalization capability and nonlinear fitting ability and the second reason may be the selection of optimal LSSVR control parameters ( and  2 ) through GSA that directly affects and improves the accuracy of the method.The powerful global search ability of GSA helps to find the optimal and suitable values for the LSSVR control parameters in a shorter time in comparison to other algorithms.We can conclude that, in application of LSSVM, the control parameters should be adequately optimized by using global optimization techniques.This will decrease the uncertainties in obtaining optimal LSSVM models.
(6) In general, the benchmark regression methods can be ranked according to their prediction accuracy as MARS, M5RT, and CMLR.The reason behind the worst results of M5RT and CMLR methods can be the linear structure of these models.(7) The RF RMSE , RF MAE , and RF DC results validate that the HLGSA and LHLGSA methods can be effectively applied for the prediction and estimation of river flow.

Conclusions
In the current study, river flow data of Astore and Shyok rivers was used to determine the forecasting capability of HLGSA, MARS, M5RT, and CMLR methods by using the for the Astore and Shyok Basins, respectively.The third part of the study evaluated the prediction performance of the HLGSA and LHLGSA methods in river flow estimation using river flow data of the nearby basin.The test results revealed that the LHLGSA performed better than the HLGSA in estimating river flows of Shyok Basin by using Astore Basin data.
In this study we forecasted river flows with only previous river flow values as inputs.These prediction accuracies of the applied methods could be improved if more input variables were available.Further studies may be conducted by including more inputs such as rainfall, snowpack, and temperature or/and building prediction models using more advanced modeling methods at these study sites.The proposed data driven methods may be applied for other regions with similar or different climates.In this case, however, the methods should be properly calibrated by using high number of river flow data.

Figure 2 :
Figure 2: Flow chart of HLGSA model for river flow forecasting.

Figure 3 :
Figure 3: The location map of the Astore and Shyok catchments.

Figure 4 :
Figure 4: Autocorrelation function of river flow series of (a) Astore catchment and (b) Shyok catchment.

Figure 5 :
Figure 5: Overall errors of forecasted river flow using all normal models for Astore catchment.

Figure 6 :Figure 7 :
Figure 6: The observed and forecasted river flow by the HLGSA, MARS, M5RT, and CMLR models using DS4 data set: Astore catchment.

Figures 7 (Figure 8 :
Figures 7(a)-7(b) compare the overall average errors statistics of all methods in prediction river flows of Shyok Basin.Similar to Astore, here the figure clearly shows that RF RMSE and RF MAE indexes of the HLGSA method give lower values in comparison to the MARS, M5RT, and CMLR

Figure 9 :
Figure 9: Overall errors of forecasted river flow using all log models for Astore catchment.
)-9(b).The figure clearly explores the performance dominancy of LHLGSA over the other corresponding methods.According to Table 5, LCMLR generally performs better than LM5RT method in mostly cases of input combinations.However, on the mean errors values basis, the LCMLR performs worse than the LM5RT method due to inaccurate forecast results of IC1 (DS1 LCMLR = 123.24> DS1 LM5RT = 70.22,DS2 LCMLR = 135.75> DS2 LM5RT = 80.68, DS3 LCMLR = 129.02 > DS3 LM5RT = 64.28,DS4 LCMLR = 96.04> DS4 LM5RT = 61.91 according to RF RMSE viewpoint) which affects the mean errors values.Figures 10(a)-10(d) illustrate the observed and forecasted river flows of the Astore Basin by using the log transform methods.The figure clearly shows that the LHLGSA method is in good fit with the original data.
)-11(b) in the form of bar graphs.The bar graphs of mean error indexes clearly show that the LHLGSA performs better than the other log methods due to having lower values of error indexes.Scatter plots of the observed and predicted river flows of Shyok catchment for all log methods by using their best model structures are shown in Figures 12(a)-12(d).

Figure 14 :Figure 15 :Figure 16 :
Figure 14: Comparison of overall errors of all log and normal models for Shyok catchment.
fit  () and   () represent the fitness value and mass of the th agent at time , respectively, whereas best() and   () and  are the gravitational and small constant,    () and    () indicate position of kth dimension of agents  and  at the  generation, and ‖  (),   ()‖ is Euclidean distance between agents  and .
The geographical location of Astore Basin is approximately between longitudes 74 ∘ , 24  and 75 ∘ , 14  E and between latitudes 34 ∘ , 45  and 35 ∘ , 38  N. The river covers a catchment area of about 3750 km 2 .Water and power development authority (WAPDA), Pakistan, has one flow gauging station, that is, Doyian in this area for flow record under Surface Water Hydrology Project (SWHP).The elevation of this gauging station is 1583 masl and its geographical location in the basin is 35 ∘ , 33  N latitude and 74 ∘ , 42  E longitude.The Shyok Basin covers drainage area of 68,458 km 2 with average basin elevation of 4940 m.WAPDA also installed one flow gauging station hydrometric station in this area for flow record at Yogo with an elevation of 2469 m and its geographical location in the catchment is 35 ∘ , 11  N latitude and 76 ∘ , 06  E longitude.The recorded monthly data of river flows of both catchments were collected through WAPDA for the duration of 1975-2006 and the total time span of this duration is 384 months.

Table 1 :
The monthly statistical parameters of river flow data for Astore and Shyok catchments.

Table 2 .
It is clear from the table that all four applied methods provide different prediction results for different DS and input combinations.In case of input combinations, IC1 comprising the two consecutive previous months' river flow values provides the worst prediction results for the HLGSA, MARS, and CMLR methods.IC2 comprising the two consecutive antecedent months' river flow values including antecedent eleventh month's flow value gives the best prediction results for the HLGSA, MARS, and M5RT methods whereas, for the CMLR method, IC3 comprising the two consecutive antecedent months' river flow values and antecedent eleventh and twelfth months' river flow value provides better performance than the other two input combinations.In the case of data sets, it is clear from the table that DS2 gives the worst forecasts results for all the regressions models including proposed HLGSA method.The reason of this is the maximum river flow value of Astore Basin's test set; DS2 (RF max = 654.9m 3 /s) is higher than the corresponding extreme value of the training DS value (see Table1).This indicates that all trained methods encounter problems in constructing extrapolation in higher value of DS2.The higher values of RF sd and RF sc parameters for the DS2 data set in comparison with other data sets can also be another reason for the worst results.It is evident that all methods provide good forecasts for the DS4 under all input combination scenarios.The best model structures for HLGSA, MARS, and M5RT methods were found for the DS4 and IC2.However, in case of CMLR method, the best model structure was found for the DS4 and IC3.The best RF RMSE (40.09 m 3 /s) for the HLGSA is better as compared to MARS, M5RT, and CMLR methods(43.46,57.26,and58.08m 3 /s), respectively.This is also true for RF MAE values where the best RF MAE for the HLGSA is 25.05 m 3 /s compared to MARS, M5RT, and CMLR methods (28.36, 29.83, and 32.70 m 3 /s), respectively.Table2clearly shows that the HLGSA method provides better prediction results than the other regression methods under all data sets and input combinations scenarios.MARS is ranked as the second best and performs better than the M5RT and CMLR whereas M5RT performs better than the CMLR under all data sets.In case of IC3, however, CMLR gives better forecasts than the M5RT for all data sets by having lower values of error indexes (RF RMSE and RF MAE ) and higher value of correlation index (RF DC ).The best prediction results of HLGSA, MARS, and M5RT for IC2 indicate that the river flow of the two preceding months and eleventh month highly affects the current month river flow.The results of IC3 justify this statement by adding the river flow data of twelfth month that negatively affect the method performance.However, the CMLR showed positive dependence on the input combinations by showing the best prediction results for IC3.
The accuracy of the HLGSA method produced by different types of KF is different.In this research, the RBKF was utilized for determining the prediction results of river flows of both basins.However, the KF has many other types, and mostly common KF types are polynomial KF, sigmoid KF, Gaussian KF, Morlet KF, Mexican hat KF, and Meyer KF.To evaluate the performance of the applied KF in this study, six different kernel functions were compared in prediction of river flows of both catchments by using the best HLGSA model structures and were reported in Table4.Table4clearly proves the superiority of the RBKF over the other kernel functions due to having smaller values of both errors indexes (RF RMSE = 40.09,RF MAE = 25.05 for Astore and RF RMSE = 109.20,RF MAE = 55.75 for Shyok) and higher value of RF DC (RF DC = 0.916 for Astore and RF DC = 0.947 for Shyok) for both basins.
DC .From Figure8, it can be clearly seen that the RF DC value of the HLGSA method is 0.947, which is higher than RF DC of MARS, M5RT, and CMLR methods (which are 0.917, 0.891, and 0.876).Selection of the proper kernel function (KF) is very important in obtaining highly accurate HLGSA method.

Table 4 :
Comparison of HLGSA performance for different kernel functions using best data set of Astore and Shyok catchments.

Table 5 :
Comparison of the LHLGSA, LMARS, LM5RT, and LCMLR methods: Astore catchment.due to having lower values of RF RMSE and RF MAE and higher values of RF DC .Similar to previous Astore application, here also the DS2 gives the worst results whereas DS4 performs the best among all data sets.The reason of the worst results of DS2 was already mentioned before.In case of Astore Basin, the best RF RMSE , RF MAE , and RF DC values (35.33 m 3 /s, 13.08 m 3 /s, and 0.921) of the LHLGSA are better than those of the HLGSA (40.09 m 3 /s, 25.05 m 3 /s, and 0.916), respectively.This is also true for the LMARS, LM5RT, and LCMLR methods where the best RF RMSE , RF MAE , and RF DC values of the LMARS, LM5RT, and LCMLR methods, respectively, are 38.78 m 3 /s, 20.88 m 3 /s, and 0.906, 46.18 m 3 /s, 24.80 m 3 /s, and 0.899, and 39.02 m 3 /s, 21.08 m 3/s, and 0.901, which are better in comparison to those of the MARS, M5RT, and CMLR methods (43.46 m 3 /s, 28.36 m 3 /s, and 0.900,Advances in Meteorology
Table 7 shows the RF RMSE , RF MAE , and RF DC values of both methods in estimating monthly river flows of Shyok Basin.The best input combination for the LHLGSA and HLGSA methods is IC3 while IC2 generally gives worse results for both methods.DS4 gives better accuracy than the other DS whereas DS1 provides the worst estimates for both methods.Table 7 clearly shows the superiority of the LHLGSA over HLGSA in case of accuracy under all data sets

Table 7 :
Comparison of the LHLGSA and HLGSA methods in estimating river flow of the Shyok catchment by using the data of Astore catchment.The figure clearly shows that the LHLGSA has higher value of RF DC (0.882 > 0.871) representing better estimation than the HLGSA method.Figure16compares the mean error indexes of both methods in estimation river flows of Shyok Basin.It can be seen from the figure that the LHLGSA shows better accuracy than the HLGSA by having lower values of errors indexes.LHLGSA decreases the overall mean RF RMSE of the HLGSA method by 4.10% in the estimation of Shyok Basin flows using the river flows of Astore Basin.
antecedent river flow values as inputs.Two error indexes (RF RMSE and RF MAE ) and one similarity index (RF DC ) were used for comparing the prediction accuracy of these methods.CV technique was used in all the applications to better see the prediction accuracy of the data sets.In the first part of the study, among four regression methods, HLGSA provided better results than the other methods in prediction of the monthly river flow data of both catchments.HLGSA improved the prediction accuracy of the MARS, M5RT, and CMLR by 8.22%, 23.15%, and 24.49% in Astore Basin, respectively, whereas, for Shyok Basin, HLGSA improved the prediction accuracy of the MARS, M5RT, and CMLR methods by 11.48%, 19.55%, and 36.19%,respectively.In the current study, radial basis kernel function was selected for HLGSA model due to its better prediction accuracy.In the second part of the study, the effect of logarithm transform function on prediction performance of all regression methods was also investigated.Results reported that after applying logarithm function on river flow time series data, all the regression methods provided better prediction accuracy for both basins.Prediction results also exposed that the HLGSA method outperformed the other methods.LHLGSA decreased the overall mean RF RMSE of the LMARS, LM5RT, and LCMLR by 11.46%, 24.16%, and 35.65%, respectively, for the Astore Basin, whereas, for Shyok Basin, LHLGSA decreases the overall mean RF RMSE of the LMARS, LM5RT, and LCMLR by 14.49%, 21.73%, and 43.84%, respectively.On the comparison of log transformed methods and normal methods, the LHLGSA reduced the overall RF RMSE of the HLGSA by 5.66% and 4.87% for the Astore and Shyok Basins, respectively.LMARS reduced the overall RF RMSE of the MARS by 2.20% and 1.52% for the Astore and Shyok Basins, respectively.The LM5RT reduced the overall RF RMSE of the M5RT by 4.41% and 2.22%