Optimizing Extreme Learning Machine for Drought Forecasting: Water Cycle vs. Bacterial Foraging

Abstract

Machine learning (ML) methods have shown noteworthy skill in recognizing environmental patterns. However, presence of weather noise associated with the chaotic characteristics of water cycle components restricts the capability of standalone ML models in the modeling of extreme climate events such as droughts. To tackle the problem, this article suggests two novel hybrid ML models based on combination of extreme learning machine (ELM) with water cycle algorithm (WCA) and bacterial foraging optimization (BFO). The new models, respectively called ELM-WCA and ELM-BFO, were applied to forecast standardized precipitation evapotranspiration index (SPEI) at Beypazari and Nallihan meteorological stations in Ankara province (Turkey). The performance of the proposed models was compared with those the standalone ELM considering root mean square error (RMSE), Nash-Sutcliffe efficiency (NSE), and graphical plots. The forecasting results for three- and six-month accumulation periods showed that the ELM-WCA is superior to its counterparts. The NSE results of the SPEI-3 forecasting in the testing period proved that the ELM-WCA improved drought modeling accuracy of the standalone ELM up to 72% and 85% at Beypazari and Nallihan stations, respectively. Regarding the SPEI-6 forecasting results, the ELM-WCA achieved the highest RMSE reduction percentage about 63% and 56% at Beypazari and Nallihan stations, respectively.

1. Introduction

Attributable to climate change, the extent and frequency of extreme weather events are changing at both local and global scales so that an increase in the severity and intensity of drought in arid and semi-arid regions has been projected [1]. Considering the duration of rainfall deficit, drought is often categorized into four classes: meteorological, agricultural, hydrologic, and socioeconomic [2]. When a natural ecosystem is injured by water deficit, this extreme condition is called ecological drought [3]. Each class has own detecting variables such as precipitation and temperature for meteorological drought [4], evaporation stress and soil moisture for agricultural drought [5], surface and subsurface water shortage for hydrological one [6], resilience of inflow-demand and water storage for socioeconomic drought [7], and difference between ecological water requirement and consumption for ecological drought [8]. As these parameters, particularly water cycle components, are highly nonlinear and non-stationary in nature, the associated drought indices have chaotic characteristics, and hence are challenging to model and forecast [9,10].

One of the most frequently asked questions in sustainable water resources management and operation is: “How often does a drought occur?” The answer to these questions is important since it determines the amount of time that we need to bridge with current water resources. Surprisingly, several approaches exist to model and predict a drought occurrence. The relevant literature shows that drought prediction models are broadly categorized into three groups: (a) physically-based (dynamical) models [11,12]; (b) statistical models including data-driven and machine learning (ML) models [13,14,15,16]; and (c) hybrid statistical–dynamical techniques [17]. Status and prospects of each type were elaborately discussed in [18] highlighting a growing tendency to use multiple models and combine their forecasts.

Fundamentals of several ML methods and their potential applications in drought forecasting have been reviewed in a recent study [10] demonstrating that an efficient drought mitigation plan and adaptation strategy needs an efficient drought forecasting model. Among stochastic, probabilistic, and ML approaches, the study has concluded that the ML approaches provide more accurate and efficient models. Given that ML methods are, by definition, data-driven, they simply capture short-term evolution of the underlying manifold. Such models, however, do not capture and explain extreme phenomena including both wet and dry extremes [19]. Regarding meteorological drought prediction in Iran, a study showed that artificial neural networks (ANNs) could be used for Standardized Precipitation Index (SPI) and Effective Drought Index prediction in Iran [20]. Satisfactory implementation of the group method of data handling, adaptive neuro-fuzzy inference system (ANFIS), generalized regression neural network, least square-support vector machine models, and ANFIS coupled with three nature-inspired optimization algorithm models for multivariate standardized precipitation index prediction in Iran were also reported by [13]. Another article demonstrated that extreme learning machine (ELM) can be used for Standardized Precipitation Evapotranspiration Index (SPEI) prediction in eastern Australia [16]. The study revealed that ELM outperforms ANN in their study area. The ANN and XGBoost methods were used for SPEI prediction using hydroclimatological predictors and showed that XGBoost outperforms ANN [21]. More recently, [22] compared the efficiency of ELM, random forest, and support vector machine models for SPEI prediction in different basins and showed that a precise forecasting model can be developed via the use of sea surface temperature as an SPEI predictor.

The purpose of any supervised ML algorithm is to get low bias and low variance. Linear ML algorithms, which assume a linear relationship between the inputs and the output variable, generally have a high bias but a low variance. In contrast, nonlinear ML algorithms, such as ANNSs and ELMs, usually have a low bias but a high variance. Thus, the parameterization of nonlinear ML algorithms often needs a trade-off to balance out bias and variance. Recent studies have indicated that different factors such as bias-variance tradeoff and number of training data may alter the performance of nonlinear supervised ML methods [10]. In addition, they might get trapped in local optimum when classic gradient-based algorithms are used to train them [23,24]. To cope with these problems, recent studies have attempted to develop hybrid ML models in which either an optimization algorithm is utilized for structural improvement of the desired model [25,26] or an effective data pre-processing approach is used to denoise or transform raw predictor variables (inputs) to more effective ones [27,28]. In the present study, the former approach is taken to develop two novel hybrid ML models based on a combination of ELM with two emerging optimization algorithms, namely bacterial foraging optimization (BFO [29]) and water cycle algorithm (WCA [30]) to investigate how much these bio-inspired techniques can improve the forecasting accuracy of ELM-based SPEI prediction model. To this end, we used 46 years (1971–2016) near surface air temperature and precipitation datasets from two synoptic stations in Ankara province, Turkey. Then, the SPEI time series in two accumulation times (3-month and 6-month) were calculated and modeled using the attained SPEI time series. Finally, a comparative analysis among these hybrid models (i.e., ELM-BFO and ELM-WCA) and the standalone ELM was performed. It is worth mentioning that the evolved hybrid algorithms are, in fact, nonlinear approximations of the short-term evolution of the system; however, we call them models throughout this study.

Although WCA and BFO algorithms have been effectively used to optimize a variety of ML models in different problems [31,32,33], current literature has indicated that SPEI forecasting using hybrid ELM-BFO and ELM-WCA models have not been explored yet. Therefore, the main contributions of this study are twofold. First, for the first time, the study demonstrates predictive capabilities of these hybrid models for SPEI forecasting in two accumulation times. Second, a comprehensive comparison among these hybrid models and the standalone ELM is presented for the case study area. The rest of this paper is organized as below: Section 2 provides the information required for SPEI calculation at each region. Section 3 the details of the proposed hybrid methods. In Section 4 and Section 5, the forecasting results/comparisons and the conclusions are presented, respectively.

2. Materials and Methods

2.1. SPEI Calculation Using the Ground Truth Data

The SPEI is a meteorological drought indicator that takes both precipitation and temperature data into account to characterize spatial and temporal drought condition [34]. The index is suitable for monitoring and predicting drought under different type of climates [35]. To calculate SPEI, climatic water balance, also known as deficit, is first calculated via subtracting potential evapotranspiration from precipitation at each month and then, a desired distribution model is fitted to monthly moisture deficit. The months with negative SPEI values are considered as dry months. When log-logistic function is selected as the desired distribution model, SPEI values in the range [−1.1 to −1.42], [−1.43 to −1.23], and less than −1.83 respectively represent the moderate, severe, and extreme drought events [36].

To attain SPEI series, long-term precipitation, and temperature data (from 1971 to 2016) at Beypazari and Nallihan synoptic stations (Figure 1) were used. As shown in Figure 2, we considered three- and six-months accumulation periods so that all the training and verification tasks will be done for SPEI-3 and SPEI-6, separately.

2.2. The Benchmark ELM Model

Despite various effective implementations of regular neural networks [37] such as feedforwards, they are criticized because parameter adjustment requires so much time [38]. In case of complicated systems, use of parameter optimization and data pre-processing techniques to obtain the best ANN structure/inputs could be necessary. Huang et al. [38] recommended a novel neural training method, called ELM, to solve such drawbacks. It applies an innovative feed-forward ANN algorithm having greater generalization capabilities and more speed than ordinary feed-forward neural networks. In fact, ELM has a single hidden layer in which the weights W and bias b of its input layer (see Figure 3) are set randomly and remain unchanged during its training process. As the input weights remain constant, the output weights (denoted by β in the figure) are independent of them and are directly determined without iteration. This is the reason responsible for ELM’s accelerated learning process [39]. As illustrated in Figure 3, the relation between the inputs ($X_{i^{\mathrm{i}}}$) and outputs ($t_i$) is expressed by Equation (1) in which f is the activation function (sigmoid function in the present study).

$$t_i+{\mathsf{\epsilon }}_i={{\displaystyle \sum }}_{j=1}^L{\beta }_{j}f\left(W_j{X}_i+b_j\right)$$

(1)

where f is a user-defined activation function (the sigmoid function was used in the present study), L is the number of hidden neurons and i indicates an individual input in the training set.

In this study, ELM was simply coded in MATLB and was trained using the training data set of each meteorology station. As the SPEI signals vary in the limited range [–3, 3], the hidden weights and bias were selected using the empirical randomization function in the range [−1, 1], and the number of hidden neurons was determined via a trial-error procedure. Such ELMs are sub-optimal and might produce a high variance error as we previously mentioned. Thus, optimizing the matrix of weights and biases of the ELM trough WCA and BFO algorithms was investigated in this study.

2.3. Water Cycle Algorithm

The WCA is a relatively new algorithm that uses the conceptual formulation of the hydrological water cycle to determine the optimal solution for parameter optimization problems [30]. The algorithm is based on the interconnectedness of streams, rivers, and a sea [40,41]. The initial assumption is that a population of raindrops is generated at random by a rainfall process, and the best raindrop is the one that produces the best fitness (i.e., minimum cost function) [41,42]. Streams, rivers, and the sea are labeled as “good or the best raindrops” based on their performance compared to the other raindrops. The best raindrop is selected as a sea and a river is considered as a good raindrop, and the rest of drops are called streams [41]. Like nature, streams flow towards the river and rivers flow to the sea. The algorithm can be used to solve any N_var-dimensional optimization problem [30].

To begin the optimization, the initial population of randomly distributed raindrops (i.e., candidate solution matrix X) is made with the dimensions N_pop × N_var.

$$\mathrm{Population}\mathrm{of}\mathrm{Raindrops}=\left[\begin{array}{c}{\mathrm{Raindrop}}_1\\ {\mathrm{Raindrop}}_2\\ {\mathrm{Raindrop}}_3\\ ⋮\\ {\mathrm{Raindrop}}_{N_{\mathrm{pop}}}\end{array}\right]=\left[\begin{array}{cccccc}{x}_1^1& x_2^1& x_3^1& & \cdots & x_{N_{\mathrm{var}}}^1\\ x_1^2& x_2^2& x_3^2& & \cdots & x_{N_{\mathrm{var}}}^2\\ ⋮& ⋮& ⋮& & ⋮& ⋮\\ x_1^{N_{\mathrm{pop}}}& x_2^{N_{\mathrm{pop}}}& x_3^{N_{\mathrm{pop}}}& & \cdots & x_{N_{\mathrm{var}}}^{N_{\mathrm{pop}}}\end{array}\right]$$

(2)

where N_var is the number of design variables that must be optimized, and N_pop is the population of raindrops (potential solutions).

As illustrated in Equation (2), each raindrop—a row in the matrix—is designated by a vector of position also known as raindrop layer. Then, the fitness value of each raindrop layer is calculated by any given cost function (here root mean square error (RMSE)) was chosen). The raindrop layer with the smallest cost function is considered the “sea”, and the raindrop layers with cost function values close to the sea are considered “rivers.” The remaining raindrop layers are considered “streams”.

The position of the streams and rivers is changed at each iteration, and then, their cost value is updated. If the cost function of a stream is less than the cost function of a river, the positions of the stream and river are swapped. Similarly, if the cost function of a river is less than the cost function of the sea, the positions of the river and sea are swapped. This process continues until the maximum number of iterations is reached or the minimum error is achieved. The final position of the sea is the optimal weight for the ELM.

To start changing the positions, the summation of the number of rivers and a single sea is calculated.

N_sr = Number of rivers + 1

(3)

where N_sr denotes number of rivers plus a sea. Therefore, the rest of the population of initial raindrops form streams, denoted by N_Raindrops, that flow into rivers or directly into the sea.

N_Raindrops = N_pop − N_sr

(4)

After that, the following formula are be used to calculate the number of streams that flows to rivers and the sea, where NS_n represents the total number of streams:

$${NS}_n=round\left\{\left|\frac{Cos{t}_n}{{{\displaystyle \sum }}_{i=1}^{N_{sr}}Cos{t}_i}\right|\times N_{Raindrops}\right\},i=1,2,3,\cdots ,N_{sr}$$

(5)

The following equations determine the movement of streams and rivers:

$${\overrightarrow{X}}_{Stream}\left(t+1\right)={\overrightarrow{X}}_{Stream}\left(t\right)+rand\times C\times \left({\overrightarrow{X}}_{Sea}\left(t\right)-{\overrightarrow{X}}_{Stream}\left(t\right)\right)$$

(6)

$${\overrightarrow{X}}_{Stream}\left(t+1\right)={\overrightarrow{X}}_{Stream}\left(t\right)+rand\times C\times \left({\overrightarrow{X}}_{River}\left(t\right)-{\overrightarrow{X}}_{Stream}\left(t\right)\right),$$

(7)

$${\overrightarrow{X}}_{River}\left(t+1\right)={\overrightarrow{X}}_{River}\left(t\right)+rand\times C\times \left({\overrightarrow{X}}_{Sea}\left(t\right)-{\overrightarrow{X}}_{River}\left(t\right)\right),$$

(8)

$$X\in \left(0,C\times d\right),\hspace{1em}C>1$$

(9)

where the constants rand and C are very important. The rand is a random number between 0 and 1 that is expanded out the same way everywhere. The C is a constant value that can be between 1 and 2. (near to 2). The d is how far away the stream is from the river right now.

The evaporation process causes seawater to be lost as freshwater rivers and streams make their way to the ocean. This assumption is made so that one does not become trapped in a situation where a local optimum exists. There must be a progressive reduction in the d_max during the optimization process. The following requirement is utilized to determine whether the river empties into the sea:

$$if\Vert {\overrightarrow{X}}_{Sea}^t-{\overrightarrow{X}}_{Rive{r}_j}^t\Vert <d_{\max }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}or\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}rand<0.1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,3,\cdots ,N_{sr}-1$$

(10)

$$d_{\max }\left(t+1\right)=d_{\max }\left(t\right)-\frac{d_{\max }\left(t\right)}{I{t}_{\max }},t=1,2,3,\cdots ,I{t}_{\max }$$

(11)

where d_max is a very small value used for controlling the “intensification level” near the sea. It_max is the maximum number of iterations. A new position of the streams is created as follows:

$${\overrightarrow{X}}_{Stream}^{New}\left(t+1\right)=\overrightarrow{LB}+rand\times \left(\overrightarrow{UB}-\overrightarrow{LB}\right)$$

(12)

where $\overrightarrow{LB}$ and $\overrightarrow{UB}$ are the lower and upper bounds of the problem, respectively.

In the training phase, weight and bias values of the ELM were optimized to hybridize. The RMSE was designated as the fitness function for minimizing the difference between measured and predicted values in the WCA algorithm.

In WCA, N_Raindrops and N_sr are the parameters that may control the number of raindrops and the number of rivers that the simulated rain system is run, respectively. The significance of these parameters is that they can affect the accuracy and efficiency of the optimization algorithm. For example, a larger value of N_Raindrops may lead to a more accurate simulation of the water cycle, but it may also increase the computational cost of the algorithm. Similarly, a larger value of N_sr may improve the convergence of the optimization, but it may also increase the run time of the algorithm. The optimal values for these parameters will depend on the specific application and the desired trade-off between accuracy and computational efficiency. The variable d_max is a small number close to zero. If the distance between a river and the sea is less than d_max, it means the river has reached the sea. In this case, the evaporation process is applied, and after sufficient evaporation, rainfall will begin. A large value for d_max reduces the search, while a small value increases search intensity near the sea. Therefore, d_max controls the search intensity near the sea to find the optimal solution. The value of d_max adaptively decreases as the search progresses.

Table 1. The initial parameters of the WCA evolutionary algorithm.

Parameter	Value
Sum of rivers and sea (NSr)	4
Evaporation condition constant (dmax)	1 × 10−5
Maximum iteration (Itmax)	100
Population (Npop)	5
Number of design variables (Nvar)	2
Search range	[−1, 1]

listed the parameters used for WCA setting in this study. The detailed pseudo-code of the complete WCA algorithm is given in Appendix A (Algorithm A1).

2.4. BFO Algorithm

BFO is another powerful metaheuristic algorithm inspired by the social foraging and chemotactic behaviors of the Escherichia coli (E. coli) bacteria [29]. Like WCA algorithm, several improvements were introduced for BFO [28]. In this study, we implemented the original version of BFO algorithm [29] that considers lifetime movements (i.e., chemotaxis step $j$, reproduction step $k$, and elimination-dispersal step $l$), and group behavior (i.e., swarming) of the E. coli bacteria to optimize the matrix of weights and biases of the ELM.

Mimicking the chemotactic behavior of E. coli, the virtual bacteria search for nutrients in a manner to maximize energy. The searching process simulates the movement of bacteria through swimming or tumbling. In swimming, they can move for a long period in a same direction, and in tumbling, a unit wake is done in a random direction. Equation (13) expresses the mathematical equation for the movement of the bacteria.

$$x_i^{j+1}\left(k,l\right)=x_i^j\left(k,l\right)+C_i\frac{\mathsf{\Delta }\left(i\right)}{\sqrt{{\mathsf{\Delta }}^T\left(i\right)\mathsf{\Delta }\left(i\right)}}$$

(13)

where $x_i^{j+1}\left(k,l\right)$ represents the new position of i-th bacterium at k-th reproductive and l-th elimination-dispersal step. The C_i represents the constant step size taken in the random direction vector $\mathrm{\Delta }$ indicates as specified by tumble. Like WCA, for the randomness, we used normally distributed random values within the domain −1 and 1.

Mimicking the group behavior of E. coli, a group of virtual E. coli cells arrange themselves in a traveling ring and each bacterium communicates in the group by sending signals. The cell-to-cell signals of the bacteria is the objective function value. The detailed pseudo-code of the complete BFO algorithm is given in Appendix A (Algorithm A2). Accordingly, the main steps of the algorithm include:

1.

Initializing parameters of elimination probability (P), swarm (population) size (S), number of chemotaxis steps (N_s; includes swimming and tumbling movements), number of reproduction step (N_re), number of elimination-dispersal step (N_ed). The N_ed denotes maximum number of iterations.

2.

Performing elimination-dispersal step.

3.

Performing reproduction step.

4.

Performing chemotaxis step.

5.

Computing health status (J_h) for each bacterium i and sort bacteria.

$$J_h={\displaystyle \sum }_{j=1}^{{\mathrm{N}}_s}J\left(i,j,k,l\right)$$

(14)

6.

Selecting the healthiest E. coli group, splitting, and repeating steps II to V until k = N_re and l = N_s

In this study, we sort the bacteria in reverse order according to theire health status. The more the bacteria are healthier, the lower the fitness value. Then, the first half of the population was survived to net iteration, and the least healthy bacteria will eventually die.

Table 2. The initial parameters of the BFO algorithm.

Parameter	Value
Population size (S)	20
Step size (Ci)	0.01
Number of chemotaxis steps (Ns)	2
Number of reproduction step (Nre)	2
Number of elimination-dispersal step (Ned)	1000
Elimination Probability (Ped)	0.9

listed the other parameters used for BFO setting in this study.

2.5. Hybrid ELM-BFO and ELM-WCA Models

As illustrated in Figure 4, two new hybrid models, called ELM-WCA and ELM-BFO, are developed for modeling and predicting SPEI-3 and SPEI-6 with the lead time of one month. To this end, the most effective SPEI lags are selected regarding mutual information criterion (MI). The input vectors and associated target SPEI at each station are then separated into training (the first 70%) and testing (the last 30%) subsets. The classic ELM is run to generate the benchmark model. In this phase, the weight and bias values of the model are stored in a matrix. Finally, BFO and WCA algorithms are performed to improve the predictive accuracy standalone ELM’s through optimizing the matrix of weights and biases of the ELM (See Equation (1)). To run BFO (WCA) optimization, the population size of 20 E. coli (four streams) and maximum number of 1000 generations (100 iteration) were used.

To differentiate and discuss the model results, several performance metrics could be used. Among potential metrices [43], RMSE and Nash-Sutcliffe efficiency (NSE) goodness of fit measures were used in the study.

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{(SPE{I}_o-SPE{I}_p)}^2}{n}}$$

(15)

$$NSE=1-\left[\frac{{{\displaystyle \sum }}_{i=1}^n{\left(SPE{I}_o-SPE{I}_p\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(SPE{I}_o-\overline{SPE{I}_o}\right)}^2}\right]$$

(16)

where n denotes the number of months applied to the train or test the models. $SPE{I}_o$, $SPE{I}_p$, and $\overline{SPE{I}_o}$ indicate the observed, predicted, and average $SPE{I}_o$.

2.6. Identification of Optimum Predictors (Lagged SPEI Vectors)

To identify the most effective predictors, both linear and nonlinear correlation metrices such as autocorrelation function MI were suggested in the literature. For the same SPEI data sets, [15] showed superiority of MI over the autocorrelation function as the former may detect high-order relations between the given time series and its preceding values. Thus, in the present study, the MI graphs developed by [15] were used to select the optimum predictors. While [15] confined the number of predictors to the optional MI threshold 0.2, we let the models use all the local maxima back to 13 lags. Therefore, the proposed models can be trained via longer history of observed drought events. Equations (17) and (18) express the functional forms of prediction models at Beypazari station. Similarly, Equations (19) and (20) express the functional forms of prediction models at Nallihan station.

$$SPEI-3_t=f \left(SPEI-3_{t-1}, SPEI-3_{t-2},SPEI-3_{t-4},SPEI-3_{t-7},SPEI-3_{t-10}\right)$$

(17)

$$SPEI-6_t=f \left(SPEI-6_{t-1}, SPEI-6_{t-2},SPEI-6_{t-3},SPEI-6_{t-4},SPEI-6_{t-8},SPEI-6_{t-12}\right)$$

(18)

$$SPEI-3_t=f \left(SPEI-3_{t-1}, SPEI-3_{t-2},SPEI-3_{t-4},SPEI-3_{t-7},SPEI-3_{t-10}\right)$$

(19)

$$SPEI-6_t=f \left(SPEI-6_{t-1}, SPEI-6_{t-2},SPEI-6_{t-4},SPEI-6_{t-7},SPEI-6_{t-13}\right)$$

(20)

where $SPE{I}_t \mathrm{and} SPE{I}_{t-i}$ are the predictand and predictors, respectively. The letter $i$ denotes lag time.

3. Results

Using the most effective lags, an ELM model as developed for each station. As previously mentioned, the best structure for ELM was determined through a trade-off analysis so that the number of hidden neurons was increased from 1 to 10 (one-by-one) to figure out the model having the highest accuracy in the testing stage. The ELM models were trained up to 200 epochs with a learning rate of 0.01 at each trial. Then, BFO and WCA were implemented to improve ELM accuracy via optimizing the matrix of weights and biases of the best ELM.

Table 3. The new models’ performance for drought prediction at two study sites.

Model (Structure)	Beypazari		Nallihan
	RMSE	NSE	RMSE	NSE
Training stage of SPEI-3
ELM	0.650	0.540	0.597	0.590
ELM-WCA	0.405	0.821	0.363	0.848
ELM-BFO	0.503	0.695	0.526	0.681
Testing stage of SPEI-3
ELM	0.825	0.481	0.801	0.438
ELM-WCA	0.436	0.829	0.463	0.812
ELM-BFO	0.618	0.652	0.663	0.615
Training stage of SPEI-6
ELM	0.494	0.736	0.434	0.780
ELM-WCA	0.219	0.948	0.205	0.951
ELM-BFO	0.0.38	0.854	0.306	0.891
Testing stage of SPEI-6
ELM	0.641	0.605	0.564	0.712
ELM-WCA	0.235	0.947	0.249	0.944
ELM-BFO	0.343	0.887	0.446	0.820

listed the most successful models’ accuracy in both training (from January 1972 to June 2003) and testing (from July 2003 to December 2016) periods. The results proved that the integrated models are significantly superior to their standalone structure. Since a hybrid method serves two or more models within a unique framework, it can benefit the advantages of the sole models while ignoring the disadvantages.

Considering the results attained for SPEI-3, the ELM-WCA is superior to ELM-BFO at both Beypazari and Nallihan stations. The Table 3 also indicates that the algorithms modeled SPEI-6 with higher accuracy than SPEI-3. This may be owing to the greater time accumulation in SPEI-6, which results in smoother time series with a smaller standard deviation.

It is necessary to investigate the significance of performance improvement achieved by the proposed hybrid models. To this end, we compared the percentage of RMSE reduction and NSE increase attained by each model in the testing phase in

Table 4. Percentage of improvement in forecasting accuracy in the testing period attained at each station.

Station	Models	SPEI-3		SPEI-6
		RMSE (%)	NSE (%)	RMSE (%)	NSE (%)
Beypazari	ELM-WCA vs. ELM	47.1	72.0	63.3	56.2
	ELM-BFO vs. ELM	25.0	35.0	46.5	46.6
Nallihan	ELM-WCA vs. ELM	42.2	85.0	55.8	32.6
	ELM-BFO vs. ELM	17.2	40.0	21.0	15.2

. Considering the SPEI-3 forecasting at Beypazari (Nallihan) station, the ELM-WCA reported a boosted drought forecasting by accuracy of 47% (42%), and 72% (85.4%), over the ELM model in terms of RMSE and NSE, respectively. Likewise, the ELM-BFO showed an augmentation of accuracy by 25% (17%), and 35% (40%) over the ELM. Regarding the SPEI-6 forecasting results, the ELM-WCA achieved the highest RMSE reduction percentage about 63% and 56% at Beypazari and Nallihan stations, respectively. The ELM-BFO also showed 46% and 21% decrease RMSE at Beypazari and Nallihan stations, respectively. The Table 4 indicates a noticeable difference in the hybrid models’ performance, particularly in terms of NSE percentage improvement. Inasmuch as the SPEI series at higher accumulation time has lower frequency, it can be concluded that standalone ELM had a considerable predictive performance for the low frequency SPEI series. Hence, one might not need a hybrid model for SPEI series modeling in higher accumulations such as SPEI-9 or SPEI-12.

Figure 5 and Figure 6 respectively show the time series and scatter plots of the SPEI predictions at two meteorological stations over the testing period. Figure 5 implies that all the ML models can detect the behavior of both SPEI-3 and SPEI-6 time series according to the actual occurrence of drought. However, for the case of SPEI-3, the models failed to detect extreme droughts events (SPEI-3 < −2.0), albeit the hybrid models are more successful than ELM. Higher density of black circles above the 1:1 line in Figure 6 implies that ELM overestimates SPEI values. Thus, the models tend to represent the study region more wet. Therefore, the use of standalone ELM for droughts forecasting, in which accurate prediction of troughs are more important than peaks is not recommended. Higher density of green triangles along 1:1 line proves the superiority of ELM-WCA over the other models implemented in this study. While the observed SPEI > 1.0, the hybrid models slightly tend to underrate the wet events.

To assess the ML models’ efficiency regarding potential phase error, which is seen in the reported forecasts in Figure 5, the RMSE values of the ML methods were compared to those of a naïve forecast—namely the forecast whereby the future SPEI value is made equal to the last observed value—model at each station (not given in this article). Considering SPEI-3, while the ELM error at Beypazari (Nallihan) is 0.825 (0.801), the corresponding naïve forecasting model error was 0.925 (0.882). Similarly regarding SPEI-6, the ELM error at Beypazari (Nallihan) is 0.641 (0.564), the corresponding naïve forecasting model error was 0.669 (0.607). It is clear that the ML models are superior to naïve forecasting model.

The Taylor diagram was used for analyzing the result of drought forecasting at both stations (Figure 7). The diagram provides a concise statistical summary of how well models’ prediction and the observed data match each other in terms of their correlation, their RMSE and the ratio of their variances. Generally, all the models provided better capability in SPEI-6 forecasting at Nallihan station, as their statistical metrics indicated between the correlation of 0.85 to 0.95 by minimum RMSE (red curves). Figure 7 shows that the results of the ELM and ELM-BFO techniques are far away from calculated SPEI values (red point). Result of SPEI-3 at Beypazari station showed that the weakest drought forecasting model has the correlation about 0.65. Also, all applied ML models showed dependable capability in drought forecasting, but the ELM-WCA (green point) reported an excellent simulated result in comparison with the counterparts at both stations. In the Taylor diagram, the distance between points of the hybrid models and the ELM indicates the ability of WCA and BFO algorithms as an enhancing mechanism for boosting the ability of ELM for drought forecasting.

4. Discussion

Over the past few years, ELM and hybrid ELM models have been successfully utilized for drought modeling and forecasting [9,15,41]. Despite the greater efficiency of ELM over ANN [9], there is a growing consensus that the standalone ELM has failed to precisely model and forecast SPEI series [15,44]. Like the present study, [44] implemented two optimized ELM techniques: namely online sequential ELM and self-adaptive evolutionary ELM for SPEI prediction in Vietnam. The authors demonstrated that the hybrid models may improve forecasting accuracy; however, the augmentation of the accuracy of the new methods are in the range of 10% to 20% in different locations. Likewise, [15] showed that optimizing ELM via Bat algorithm can improve the forecasting accuracy by 20%. Comparing to these studies, it is concluded that WCA showed the best performance to optimize ELM. From a model accuracy perspective, wise selection of predictors is also an important factor. Although the main mechanism for the outperformance of the hybrid models is optimizing the matrix of weights and biases of the ELM, the higher performance of ELM-WCA over that of previous studies, might partly be due to the way of selected input parameters. Thus, identifying the most effective inputs needs to be considered in developing hybrid models.

5. Conclusions

Precise prediction of droughts is known as a challenging issue in hydrology and sustainable watershed management. To address the problem, the coupled ML approaches with optimization algorithms were recommended in the relevant literature as the potential boosted tools for drought forecasting. However, efficiency of these techniques is highly depending on climate of study areas (i.e., observed datasets). In this paper, two new coupled ML approaches (namely ELM-WCA and ELM-BFO) were suggested for forecasting meteorological drought events. The proposed models were implemented for one month ahead forecast of SPEI-3 and SPEI-6 at two meteorological stations in Turkey. They were trained and tested using inputs selected via MI approach. Our findings showed that the proposed hybrid models are more accurate than the benchmark. Regarding the SPEI-3 forecasting, the proposed ELM-WCA and ELM-BFO models showed approximately 25% and 17% RMSE reduction over the ELM at Beypazari and Nallihan stations, respectively. Correspondingly for the SPEI-6, the hybrid models showed approximately 46% and 21% RMSE reduction over the ELM at Beypazari and Nallihan stations, respectively. However, the accuracy of drought forecasting indicated that there were some limitations in the cases of SPEI-3, then the current study recommends that application of pre-processing approaches on the dataset before modeling process maybe increase accuracy of forecasts.

As previously discussed, to improve ELM accuracy, one needs to identify optimum values of ELM parameters including total number of hidden neurons and maximum number of iterations. In this study, the WCA and BFO algorithms merely used to optimize the matrix of weights and biases of hidden-to-output-layer of the ELM. Future studies may apply similar algorithms for fine tuning of the other ELM parameters. Regarding the possible applications of the suggested models, it is worth mentioning that the provided technique may be used for comparable forecasting problems involving various stochastic hydrological phenomena. Future research may be necessary to explore the suggested model’s capabilities to forecast agricultural and hydrological drought, since the present work was limited to meteorological drought forecasting. Before ELM training process, one may use other optimization approaches or a decomposition strategy to improve drought prediction accuracy. All the models implemented in this study are from shallow learning neural networks family. For longer lead time forecasting, use of deep learning models could be recommended as a topic for future studies. Investigation on the effect of modified versions of WCA [45] or BFO algorithms on ELM optimizing could also be informative.