**About the Editors**

### **Mahdi Hasanipanah**

Dr. Mahdi Hasanipanah obtained his BSc degree in Mining Engineering from the University of Kashan, Iran, in 2008; his MSc degree in Mining Engineering from the Islamic Azad University—South Tehran Branch, Iran, in 2010; and his PhD in Mining Engineering from University of Kashan, Iran, in 2020. His research interests include rock blasting, rock mechanics, machine learning methods, and optimization algorithms. He has cooperated in a large number of projects in Iran and Malaysia. He is also the chief Guest Editor of *Sustainability* (MDPI)'s Special Issue "Advances in Rock Mechanics and Geotechnical Engineering". His citation and H-index are 5600 and 49, respectively, which indicate the significant impact and influence of his research in the academic and scientific communities. He was among the top 2% of scientists for two consecutive years, including 2022 and 2023, according to Stanford University.

#### **Danial Jahed Armaghani**

Dr. Danial Jahed Armaghani is a prominent researcher in the field of civil and geotechnical engineering. His work has made significant contributions to the understanding and mitigation of geotechnical and geological hazards, earning him a reputation as an excellent researcher in his field. Dr. Danial's research focuses on a wide range of topics, including slope stability analysis, rock mechanics, tunnel construction, surface and deep excavations, and applying machine learning models and optimization algorithms to solve various geotechnical problems. Dr. Danial obtained his MEng (2012) and PhD (2015) in the area of Civil-Geotechnics, from Universiti of Teknologi Malaysia (UTM), Malaysia. He is currently working as a Postdoctoral Research Fellow at the School of Civil and Environmental Engineering, University of Technology Sydney (UTS), Australia. He published more than 300 articles in well-established ISI and Scopus journals and at national and international conferences. His published works have received more than 19,000 citations, which indicate the significant impact and influence of his research in the academic and scientific communities. He was among the top 2% of scientists for four consecutive years from 2020 to 2023, according to Stanford University.

### **Jian Zhou**

Assoc. Prof. Jian Zhou obtained his BSc degree (2008) and PhD (2015) from Central South University (CSU), China, and was a visiting scholar with Mine Design Laboratory at McGill University from 2013 to 2014. Currently, he is an associate professor in the School of Resources and Safety Engineering at CSU, China. His current research interests include geological and geotechnical hazards prediction and mitigation, applying predictive models in rock mechanics, and mining engineering. Dr. Zhou is the Highly Cited Researcher in the field of Cross-Field (Clarivate); a Highly Cited Chinese Researcher in the field of Mining Engineering (Elsevier); one of the World's Top 2% Scientists according to Stanford University; and a recipient of the Distinguished Young Scholars Fund of Hunan Province, China. He has published more than 180 papers in international journals on mining and geotechnical issues and has received China's 100 Most Influential International Academic Papers Award, the *Journal of Rock Mechanics and Geotechnical Engineering's* Best Paper Award, and the *Journal of Central South University's* Best Paper Award. His citation and H-index are 9600 and 55, respectively.

**Xiaohua Ding 1,2, Mehdi Jamei 3,4, Mahdi Hasanipanah 5,6,\*, Rini Asnida Abdullah <sup>6</sup> and Binh Nguyen Le 5,7**


**Abstract:** Using explosive material to fragment rock masses is a common and economical method in surface mines. Nevertheless, this method can lead to some environmental problems in the surrounding regions. Flyrock is one of the most dangerous effects induced by blasting which needs to be estimated to reduce the potential risk of damage. In other words, the minimization of flyrock can lead to sustainability of surroundings environment in blasting sites. To this aim, the present study develops several new hybrid models for predicting flyrock. The proposed models were based on a cascaded forward neural network (CFNN) trained by the Levenberg–Marquardt algorithm (LMA), and also the combination of least squares support vector machine (LSSVM) and three optimization algorithms, i.e., gravitational search algorithm (GSA), whale optimization algorithm (WOA), and artificial bee colony (ABC). To construct the models, a database collected from three granite quarry sites, located in Malaysia, was applied. The prediction values were then checked and evaluated using some statistical criteria. The results revealed that all proposed models were acceptable in predicting the flyrock. Among them, the LSSVM-WOA was a more robust model than the others and predicted the flyrock values with a high degree of accuracy.

**Keywords:** blast-induced flyrock; LSSVM; optimization; prediction models

### **1. Introduction**

Drilling and blasting is an indispensable technique for breakage and displacement of rock masses in open-pit mines. Nevertheless, some undesirable phenomena, such as ground vibration, airblast, flyrock (FR), and backbreak are produced by blasting operations [1–5]. Any blasting event produces a sudden ejection of rock pieces, which are referred to as "FR". This phenomenon is one of the most hazardous environmental issues induced by blasting which may lead to various problems for humans, including fatalities [6–9]. As mentioned in previous studies, some blast design factors, such as burden (*B*), spacing (*S*), stemming (*ST*), weight charge (*WC*), and powder factor (*PF*), are the effective factors in the intensity of *FR* [8–10]. Aside from the aforementioned factors, the properties of the rock mass, such as rock density and uniaxial compressive strength, are considered the effective factors on the *FR*, called uncontrollable factors [9,10]. The *FR* can occur based on three different mechanisms, i.e., rifling, face burst, and cratering [11,12]. The poor *ST* material, the small ratio of *ST* to blast-hole diameter, and inadequate *B* are the most important causes for the rifling, cratering, and face burst mechanisms [9,13].

**Citation:** Ding, X.; Jamei, M.; Hasanipanah, M.; Abdullah, R.A.; Le, B.N. Optimized Data-Driven Models for Prediction of Flyrock due to Blasting in Surface Mines. *Sustainability* **2023**, *15*, 8424. https://doi.org/10.3390/su15108424

Academic Editor: Rajesh Kumar Jyothi

Received: 6 April 2023 Revised: 13 May 2023 Accepted: 18 May 2023 Published: 22 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In the literature, numerical simulations are considered to be the common methods to study the blasting mechanism in rock masses. According to Kutter and Fairhurst [14], three main zones, i.e., the crushed, cracked, and elastic vibration zones, can be formed after each blasting event. Additionally, the in-situ stress has an important effect on the propagation of cracks produced by blasting. To study the failure responses of rocks, the distinct element method (DEM), finite difference method (FDM), and finite element method (FEM) have been extended in recent years by many scholars [15–20]. As an example, a twodimensional discrete element method was used to numerically simulate the mechanism of rock fragmentation produced by blasting in the study conducted by Hajibagherpour et al. [19]. They showed that the proposed numerical model can be effectively employed to simulate the crack propagation process around a blast-hole. Aside from numerical modelling, several empirical models have been employed to predict flyrock [20]. These empirical models have been formulated based on considering only one or two of the effective factors of flyrock. For this reason, the accuracy of the mentioned empirical models is not good enough. Therefore, the use of artificial intelligence methods can be a good solution to predict *FR* with a high degree of performance. Additionally, the use of artificial intelligence methods in different fields of mingling and civil engineering indicates the effectiveness of these methods for predicting and optimizing aims [21–31].

An artificial neural network (ANN) model was employed to predict *FR* in the study conducted by Monjezi et al. [32], and its performance was compared with statistical models. Their results indicated the performance of ANN was better than statistical models in predicting *FR*. For the same purpose, Ghasemi et al. [33] employed ANN and fuzzy system (FS) and showed better prediction capability of FS over ANN. Moreover, the ANN model was compared with the adaptive neuro-fuzzy inference system (ANFIS) for the prediction of *FR* by Trivedi et al. [34]. They revealed higher performance in respect to accuracy of ANFIS compared with the ANN model. In another study, a genetic programming (GP) model was employed by Faradonbeh et al. [35] to predict *FR*. In their study, the non-linear regression models were also used for comparison aims. They concluded that the GP predicted *FR* with a higher performance in comparison to non-linear regression models. The GP model was also employed by Ye et al. [36] and its results was compared with a random forest model. According to their results, the performance of GP was better than the random forest model.

The present study attempts to propose several efficient hybrid models through the cascaded forward neural network (CFNN) and also the least squares support vector machine (LSSVM) in combination with three optimization algorithms, including artificial bee colony (ABC), gravitational search algorithm (GSA), and whale optimization algorithm (WOA), for the prediction of *FR*.

The rest of this article includes the following sections. More details about the source of the database and the developed models are explained in the second section. Then, the setting parameters in the modelling processes are explained in the third section. The results and discussions are provided in the fourth section; finally, the last section presents the conclusions of the study.

#### **2. Research Significance**

FR is considered as an environmental and hazardous problem in mine blasting, which may result in human injuries, fatalities, property damage, and instability of slopes. Hence, a valid and reliable prediction of *FR* has critical implications in mitigating and controlling the adverse effects along with sustainable development and responsible mining. In other words, the control and minimization of *FR* can lead to sustainability of surroundings environment in blasting sites. For the aforementioned aims, the present study attempts to propose several efficient hybrid models through the CFNN and also LSSVM in combination with three optimization algorithms. To the best of our knowledge, this is the first work that predicts the *FR* by using the proposed models.

#### **3. Materials and Methods**

#### *3.1. Materials*

The database used in this study was collected from three granite quarry sites located in Malaysia, including the Ulu Tiram, Pengerang, and the Masai quarry sites. The values of the rock quality designation (RQD) of the aforementioned quarry sites ranged from 45 to 80, 50 to 70, and 40 to 75, respectively. Additionally, the values of the rock strength ranged from 30 to 110 MPa, respectively. In total, 80 datasets including some effective parameters on the *FR* were used in constructing the predictive models. In this regard, the *S*, *B*, *ST*, *PF*, and density were used as the input parameters, and the *FR* was used as the output parameter. More details about the statistical properties of datasets will be provided in Section 5.

### *3.2. Methods*

In this study, the LSSVM is combined with the ABC, GSA, and WOA to predict *FR*. Additionally, the CFNN model is also used for comparison aims. In this section, the mentioned models are briefly explained.

#### 3.2.1. LSSVM Model

LSSVM is a robust machine learning technique. This method was proposed as an upgraded form of the SVM, which suffered from some drawbacks in its learning stage, namely the demand in calculability and the limitation in dealing with inequality constraints. Therefore, LSSVM becomes an efficient ML technique after fixing the aforesaid issues [37,38]. For a regression task which aims at finding a suitable correlation that emulates the behaviour of a system defined by a set of data having inputs *x* = {*x*1, *x*2, . . . , *xN*} that *x<sup>j</sup>* ∈ *R <sup>D</sup>* and *N* is the number of samples in the set, and targets *t* defined on *R* as *y* = {*y*1, *y*2, . . . , *yN*}, the first step in the LSSVM method consists of formulating the following minimization problem [39]:

$$\begin{aligned} \text{minimize } & \frac{1}{2} w^T w + \frac{1}{2} \gamma \sum\_{j=1}^N \left( e\_j^2 \right) \\ \text{s.t. } & y\_j = w^T \boldsymbol{\varrho}(\boldsymbol{x}\_j) + \boldsymbol{b} + e\_{j\prime} \ \boldsymbol{i} = 1, 2, \dots, N \end{aligned} \tag{1}$$

where *e<sup>j</sup>* denotes the regression error, *γ* represents the regularization parameter, *<sup>T</sup>* points out the transpose operator, *w* and *b* are the weight and bias parameters, respectively, and *ϕ* is a nonlinear mapping function.

The learning phase of LSSVM passes through finding the proper values of *w* and *b*. To this end, the formulated minimization problem is transformed into a Lagrangian function using the formula shown below [40]:

$$L(w, b, a, e) = \frac{1}{2}w^T w + \frac{1}{2}\gamma \sum\_{j=1}^{N} \left(e\_j^2\right) - \sum\_{j=1}^{N} a\_j \left(w^T \boldsymbol{\varphi}(\mathbf{x}\_j) + b + e\_j - y\_j\right) \tag{2}$$

where the coefficients *α<sup>i</sup>* are called Lagrangian multipliers. The solution of L is obtained by solving the following system of equations:

$$\begin{cases} \frac{\partial L(w, b, a, \varepsilon)}{\partial w} = 0 \Rightarrow w = \sum\_{j=1}^{N} \mathfrak{a}\_{j} \eta\_{j}(\mathfrak{x}\_{j})\\ \frac{\partial L(w, b, a, \varepsilon)}{\partial b} = 0 \Rightarrow \sum\_{j=1}^{N} \mathfrak{a}\_{j} = 0\\ \frac{\partial L(w, b, a, \varepsilon)}{\partial \varepsilon\_{j}} = 0 \Rightarrow \mathfrak{a}\_{j} = \gamma e\_{j}, \quad j = 1, 2, \dots, N\\ \frac{\partial L(w, b, a, \varepsilon)}{\partial \mathfrak{a}\_{j}} = 0 \Rightarrow w^{T} \varrho(\mathfrak{x}\_{j}) + b + e\_{j} - y\_{j} = 0, \quad j = 1, 2, \dots, N \end{cases} \tag{3}$$

The above system which defines the vanishing of the partial derivatives of *L* with regard to *w*, *b*,*e*, and *α* can be arranged in the following matrix scheme:

$$
\begin{bmatrix}
\Omega & \mathbf{1}\_N^T \\
\mathbf{1}\_N & \Omega + \gamma^{-1} I\_N
\end{bmatrix}
\begin{bmatrix}
b \\
\mathbf{a}
\end{bmatrix} = \begin{bmatrix}
\mathbf{0} \\
\mathbf{y}
\end{bmatrix} \tag{4}
$$

In the above equation, *I<sup>N</sup>* points out *N* × *N* size identity matrix, *y* = [*y*1, *y*2, . . . , *yN*] *T* , *α* = [*α*1, *α*2, . . . , *αN*] *T* , 1*<sup>N</sup>* = [1, 1, . . . , 1] *T* , and Ω is the kernel matrix. The elements of this latter term are expressed as follows:

$$
\Omega\_{j,l} = \varrho(\mathfrak{x}\_{\mathfrak{l}})\,\varrho(\mathfrak{x}\_{\mathfrak{l}}) = \mathcal{K}(\mathfrak{x}\_{\mathfrak{l}}, \mathfrak{x}\_{\mathfrak{l}}) \tag{5}
$$

where *K* is the kernel function. Gaussian radial basis function (RBF) is among the frequently considered kernel functions in LSSVM.

In the last step, the gained LSSVM paradigm can predict the investigated target using the following expression:

$$f(\mathbf{x}) = \sum\_{j=1}^{N} \alpha\_j \mathbb{K}(\mathbf{x}\_{\mathbf{j}}, \mathbf{x}\_l) + b \tag{6}$$

where (*α<sup>j</sup>* , *b*) are determined from Equation (4).

It is worth mentioning that the robustness of LSSVM is related to the proper selection of its hyper-parameters, viz.,

*σ* <sup>2</sup> and *γ*. To do so, in the present work, three rigorous metaheuristic algorithms were suggested to tune these control parameters.

In this study, three optimization algorithms, including the GSA, WOA, and ABC, are used to improve the LSSVM performance. The aforementioned algorithms are briefly explained in this part.

#### (A) Gravitational search algorithm (GSA)

The GSA is a metaheuristic algorithm developed by Rashedi et al. [41] based on Newton's law of gravity [37]. The GSA is a population-based algorithm, and this means that a population of possible solutions is considered during the optimization process. The particles of the population are subjected to positions updating using the main governing equations of the algorithms. In this regard, the position of each particle is denoted by a vector *x*, while the force between two elements *i* and *j* at iteration *g*, is expressed as follows [41]:

$$F\_{ij}^{\mathcal{S}} = G^{\mathcal{S}} \frac{\mathcal{M}\_i^{\mathcal{S}} \, \mathcal{M}\_j^{\mathcal{S}}}{R\_{ij}^{\mathcal{S}} + \varepsilon} \left(\mathfrak{x}\_i^{\mathcal{S}} - \mathfrak{x}\_j^{\mathcal{S}}\right) \tag{7}$$

where *ε* points out a constant with a small value, *R* denotes the Euclidian distance between the two particles, while *G* represents the gravitational constant defined as:

$$G^{\mathcal{S}} = G^{\mathcal{S}0} \frac{g\_0^{\mathcal{X}}}{\mathcal{S}} \quad \mathcal{X} < 1 \tag{8}$$

where *G <sup>g</sup>*<sup>0</sup> is the initial value of the gravitational constant. The overall force resulted from the particles of the population on each particle *i* is determined using the following formula [41]:

$$F\_i^{\mathcal{S}} = \sum\_{\substack{j \in I\_{\text{best}} \ j \neq i}}^N r\_{1j} F\_{i\mathbf{j}}^{\mathcal{S}} \tag{9}$$

where *Jbest* represents a set of best particles in the population and *r*1*<sup>j</sup>* is a random determined uniformly over iterations from [0, 1].

In another step, the motion law is considered as per following formulas [41]:

$$a\_{\bar{l}} = \frac{F\_{\bar{l}}^{\mathcal{S}}}{M\_{\bar{i}}^{\mathcal{S}}} \tag{10}$$

where *a<sup>i</sup>* points out the acceleration of mass and *M<sup>i</sup>* is the inertia mass which is determined using the equation below:

$$M\_i^{\mathcal{S}} = \frac{m\_i^{\mathcal{S}}}{\sum\_{j=1}^{N} m\_j^{\mathcal{S}}} \tag{11}$$

and

$$m\_i^{\mathcal{S}} = \frac{f\_i^{\mathcal{S}} - w^{\mathcal{S}}}{b^{\mathcal{S}} - w^{\mathcal{S}}} \tag{12}$$

where *f* is the fitness value of the element *i*, and w and b represent the worst and best fitness values in the population, respectively. Finally, the velocity and the position of the elements are:

$$
v\_{i}^{g+1} = r\_{2i}v\_{i}^{g} + a\_{i}^{g} \tag{13}$$

$$\mathbf{x}\_{i}^{g+1} = \mathbf{v}\_{i}^{g+1} + \mathbf{x}\_{i}^{g} \tag{14}$$

where *r*2*<sup>i</sup>* is a random generated uniformly from [0, 1], and *v* and *x* point out the velocity and position of elements, respectively.

The steps of GSA based on the stated equations are repeated until a stopping criterion is fulfilled.

#### (B) Whale optimization algorithm (WOA)

The WOA is another population-based algorithm introduced by Mirjalili and Lewis [42]. The main steps and the governing equations of WOA mimic the hunting process of back whales [37]. Initially, an initial population of whales is created randomly. The positions of the whales represent possible solutions of the optimization problem. In order to evaluate the quality of these positions, a fitness function that emulates the objective function to optimize is applied. Based on the evaluation step, the whales are subjected to update their positions at a given generation (*g* + 1). To do so, the associated shapes of the position are changed to spiral or circular forms with respect to a probability *p* using the following equation [42]:

$$X^{g+1} = \begin{cases} D'e^{bl}\cos(2\pi t) + X^{\*g} & \text{if } \begin{array}{l} p \ge 0.5\\ \end{array} \\\ X^{\*g} - AD & \text{if } \begin{array}{l} p < 0.5 \end{array} \end{cases} \tag{15}$$

In the above equation, *X* ∗ points out the best whale which is located nearby the prey. *D*<sup>0</sup> = |*X<sup>i</sup>* − *X* ∗ | represents the distance between the whale *i* and the prey, *b* is a constant for specifying the spiral shape, and *l* is a random number from [−1, 1]. The other terms of the above equation are defined as per following equations [42]:

$$A = 2ar - a \tag{16}$$

$$\mathbb{C} = \mathbf{2}r \tag{17}$$

$$D = \left| \mathbb{C}X^\* - X\_i \right| \tag{18}$$

where *a* is a number decreasing linearly from 2 to 0 over the distance, and *r* is a random from [0, 1]. According to the value of *A*, if it is mainly outside the interval [−1, 1], a circular form is considered for upgrading the positions based on a randomly selected whale *X g r* . The following equation shows this process:

$$X^{g+1} = X\_r^g - AD \tag{19}$$

The new positions obtained after carrying out the update of the shapes are assessed using the fitness function. Lastly, if the best whale shows an enhancement in its fitness values, it will change its position to this newest one, otherwise, the best existing position will be conserved.

The described steps of this optimization algorithm are repeated until a stopping criterion is achieved.

(C) Artificial Bee Colony (ABC)

The ABC is an intelligent swarm-based metaheuristic algorithm proposed by Karaboga [43]. As indicated in its appellation, ABC emulates the steps performed by honeybees as they search for nectar sources. In this regard, three groups, including employed, onlookers, and scout bees are considered in this process. The role of employed group is to amass the information and expose it to the onlooker group. The scout group consists of changing the positions once no improvement is noticed in some sources of bees.

The main steps of ABC are given below [44]:


$$\mathbf{x}\_{i}^{g+1} = \mathbf{x}\_{i}^{g} + \theta\_{i} \left(\mathbf{x}\_{i}^{g} - \mathbf{x}\_{\omega}^{g}\right) \tag{20}$$

where *ϑ<sup>i</sup>* is a random from [0, 1] and *<sup>ω</sup>* <sup>∈</sup>{1, 2, . . . , *colony size*}*ω*6=*<sup>i</sup>* .

Afterward, the quality of each new position is examined using the fitness function, and if an improvement is obtained, this new position is conserved, otherwise, it is abandoned.


$$P\_i = \frac{ft\_i}{\sum\_{i=1}^{E} ft\_i} \tag{21}$$

where *f t* points out the fitness value and *E* represents the number of employed bees. As in the case of employed bees, if an improvement is obtained, this new position is conserved, otherwise, it is abandoned.


The optimum solution of the problem is represented by the fittest bee. The abovedescribed steps are recurring until a stopping condition is reached.

### 3.2.2. CFNN Model

The CFNN is a type of ANN which is recognized by its flexible-based structure [45]. This advantage allows CFNN to generate accurate predictive models for many systems with different degrees of complexity. The structure of CFNN neurons is distributed into three kinds of layers, including input, hidden, and output layers [46]. The input layer receives the data, then this latter is transformed and processed in one or more hidden layers using the so-called activation functions (such as *tansig* and *logsig*), while the results of the paradigm are obtained from the output layer. The number of hidden layers and their involved neurons depend on the complexity of the system as one hidden layer is generally enough for low to medium complicated cases, while more than one hidden layer is requested for highly complex cases. CFNN is characterized by its specific cascaded scheme for linking the neurons to the others [46]. This scheme is ensured by linking each neuron from a preceding layer to the nodes of the subsequent layers [46].

As the other kinds of ANN, the learning phase of CFNN aims at achieving the suitable weight and bias values of its architecture. Backpropagation-based algorithms, such as the Levenberg–Marquardt algorithm (LMA), are known to be highly efficient for this kind of optimization. In this investigation, LMA algorithms were considered in the optimization of the bias and weights of CFNN. More information about LMA can be found in published literature [47,48].

It is worth mentioning that the considered soft computing approaches in this study, namely CFNN and LSSVM, differ from each other mainly on the learning strategy where in LSSVM, the learning process is done after the formulation of the minimization problem (shown in Equation (1)) and then the problem is resolved by finding the control parameters of the model, while in CFNN, the learning approach is gained by finding the suitable topology and the appropriate weights linking between the neurons of different layers until reaching the low function error (such as root mean square error) value.

#### 3.2.3. Model Performance Evaluation

As stated in the previous section, the modelling task of *FR* using CFNN and LSSVM was extended in this study by investigating the suitable input parameters that can give the most accurate predictions of this vital factor. Before showing the main finding of the modelling tasks using the aforesaid ML models, it is worth mentioning that during the performance evaluation of the models, the following statistical indexes were calculated using the equations shown below [1,4,5,49–52]:

Average Absolute Relative Error (*AARE*):

$$AARE\% = \frac{1}{N} \sum\_{i=1}^{N} \left| \frac{FR\_{imea} - FR\_{ipred}}{FR\_{imea}} \right| \times 100\tag{22}$$

Coefficient of Determination (*R* 2 ):

$$R^2 = 1 - \frac{\sum\_{i=1}^{N} \left( FR\_{i \text{mean}} - FR\_{i \text{pred}} \right)^2}{\sum\_{i=1}^{N} \left( FR\_{i \text{pred}} - \overline{FR} \right)^2} \tag{23}$$

Root Mean Square Error (*RMSE*):

$$RMSE = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( FR\_{ima} - FR\_{i\,pred} \right)^2} \tag{24}$$

In the above equations, the subscripts *mea* and *pred* denote the measured and estimated *FR*, respectively, *FR* mean of *FR* values, and *N* points out the number of samples.

#### **4. Development of Predictive Models**

The main purpose of the suggested ML learning methods in this study was to deliver robust models that can estimate the *FR* under different circumstances. For a better investigation, several input parameters were considered in each of the proposed paradigms. Accordingly, six different schemes were involved in the development of these ML-based predictive models. These schemes are summarized in the following equations:

$$f(M\_1) \qquad \quad \quad \quad FR = f(\mathcal{S}\_\prime \mathcal{B}\_\prime \mathcal{S} \mathcal{T}\_\prime \mathcal{P} \mathcal{F}\_\prime \mathcal{D} \textit{ensity}) \tag{25}$$

$$f\_1(M\_2) \qquad \qquad FR = f(S, B, ST, PF) \tag{26}$$

$$f(M\_3) \qquad \qquad FR = f(\mathcal{S}, \mathcal{B}, ST, \; Density) \tag{27}$$

$$f(M\_4) \qquad \qquad FR = f(S, B, PF, Density) \tag{28}$$

$$f(M\_5) \qquad \qquad FR = f(\mathcal{S}\_\prime \mathcal{S} T\_\prime \mathcal{P} \mathcal{F}\_\prime \mathcal{D}ensity) \tag{29}$$

$$f\left(M\_{\text{f}}\right) \qquad \qquad FR = f\left(B, ST, PF, Density\right) \tag{30}$$

In order to properly implement CFNN-LMA and the proposed hybridization LSSVMmetaheuristic algorithms, including the LSSVM-WOA, LSSVM-GSA, and LSSVM-ABC models, some necessary steps were carried out as shown in the two flowcharts of Figures 1 and 2, respectively. According to these flowcharts, the first step consisted of normalizing the database. This step significantly improves the performance of the considered ML techniques. The normalization procedure of the database is given in the following equation: *Sustainability* **2023**, *15*, x FOR PEER REVIEW 8 of 20 considered ML techniques. The normalization procedure of the database is given in the fol‐ lowing equation:

$$\chi\_n = \frac{2(\mathbf{x}\_i - \mathbf{x}\_{\min})}{(\mathbf{x}\_{\max} - \mathbf{x}\_{\min})} - 1 \tag{31}$$

where *x* and *x<sup>n</sup>* point out the variable and the normalized value, respectively, while *xmax* and *xmin* represent the maximum and minimum values of the variable, respectively. where *x* and point out the variable and the normalized value, respectively, while ௫ and represent the maximum and minimum values of the variable, respectively.

**Figure 1. Figure 1.** The workflow of the suggested CFNN paradigm. The workflow of the suggested CFNN paradigm.

**Figure 2.** The workflow of the suggested LSSVM paradigms. **Figure 2.** The workflow of the suggested LSSVM paradigms.

After performing the normalization of the data, the latter was divided into train and test sets. The aim of these two groups was to train the models (train set) and certify their robustness of unseen measurements (test set). These two sets covered 80% and 20% of the whole measurements, respectively. During the learning phase of the CFNN‐ and LSSVM‐based models, their control parameters were investigated using the above‐discussed algorithms and some other techniques. In this regard, the trial and error method was considered for selecting the best topology of CFNN, while LMA was ap‐ plied for optimizing the weights and bias values of the network. For the LSSVM model, three metaheuristic algorithms including ABC, GSA, and WOA were implemented in the optimization of the two impacting LSSVM control parameters, namely and ଶ. It is necessary to add that it was proven in some previous works that it is more suitable to consider some specific trust region algorithms such as Levenberg–Marquardt algorithms rather than applying metaheuristic algorithms in the training phase of some feedforward networks such as MLP and CFNN [53]. However, for many other soft‐computing ap‐ proaches such as LSSVM and SVM, metaheuristic algorithms are much more appropriate for finding the control parameters of these techniques [54]. Table 1 reports the main After performing the normalization of the data, the latter was divided into train and test sets. The aim of these two groups was to train the models (train set) and certify their robustness of unseen measurements (test set). These two sets covered 80% and 20% of the whole measurements, respectively. During the learning phase of the CFNNand LSSVM-based models, their control parameters were investigated using the abovediscussed algorithms and some other techniques. In this regard, the trial and error method was considered for selecting the best topology of CFNN, while LMA was applied for optimizing the weights and bias values of the network. For the LSSVM model, three metaheuristic algorithms including ABC, GSA, and WOA were implemented in the optimization of the two impacting LSSVM control parameters, namely *γ* and *σ* 2 . It is necessary to add that it was proven in some previous works that it is more suitable to consider some specific trust region algorithms such as Levenberg–Marquardt algorithms rather than applying metaheuristic algorithms in the training phase of some feedforward networks such as MLP and CFNN [53]. However, for many other soft-computing approaches such as LSSVM and SVM, metaheuristic algorithms are much more appropriate for finding the control parameters of these techniques [54]. Table 1 reports the main setting of these metaheuristic algorithms. In our suggested workflows (Figures 1 and 2), the constraints and/or evalua-

setting of these metaheuristic algorithms. In our suggested workflows (Figures 1 and 2),

tion of overfitting were considered by checking the accuracy of the paradigms during both training and testing phases. It is clear from the statistical evaluation of the models that if the prediction accuracy of the latter is very satisfactory during the training and testing phases, the overfitting issue is avoided.


**Table 1.** The considered control parameters of the three employed metaheuristic algorithms.

### **5. Results and Discussion**

### *5.1. Exploratory Analysis*

The optimal strategy for solving the high non-linear problems significantly depend on the behavior of the dataset used in the simulation process. Thus, the statistical data description is one of the most crucial tasks of the pre-processing stage in ML-based modelling in engineering problems. Table 2 lists the statistical properties of datasets implemented in the *FR* predicting procedure. The low values of skewness and kurtosis confirmed that all inputs and targets are categorized as a pseudo-normal distribution. Figure 3 demonstrates the Pearson correlation coefficients in form of a correlogram. It can be concluded that the "*S*" parameter with respect to the largest correlation coefficient (0.65) has the most significance in prediction of the *FR* value. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 11 of 20

**Figure 3.** Correlogram for assessing the Pearson correlation coefficient between input and target parameters. Green: Histogram, Red: Regression line, and Blue: Distribution points. **Figure 3.** Correlogram for assessing the Pearson correlation coefficient between input and target parameters. Green: Histogram, Red: Regression line, and Blue: Distribution points.

As stated before, the database considered in this study was divided into training and testing sets. The first set is applied in the models' development, while the test set is de‐

By performing the steps described in the previous sections, it was found that 3 hid‐ den layers with *tansig* as an activation function, and 12, 11, and 9 neurons in each of them, respectively, represented the proper CFNN topology in all of the six schemes. For LSSVM models, the achieved ଶand values using the ABC, GSA, and WOA ranged between

The statistical evaluation of the performance of the obtained ML‐based models with respect to the stated six schemes in the previous sections is shown in Table 3. In this table, statistical criteria, namely *AARE*, *R*2, and *RMSE* are reported for the training set, the test set, and the whole dataset. According to this table, and based on the schemes, it can be seen that *M*<sup>1</sup> is the best one, followed by *M*2. In the combination of *M*<sup>1</sup> including all inputs for the testing phase, the LSSVM‐WOA in terms of (*R*<sup>2</sup> = 0.999, *RMSE* = 3.4209 m, and *AARE* = 1.3017) was identified as the superior predictive model, followed by CFNN‐LMA (*R*<sup>2</sup> = 0.9347, *RMSE* = 16.5215 m, and *AARE* = 7.512), LSSVM‐GSA (*R*<sup>2</sup> = 0.904, *RMSE* = 16.1775 m, and *AARE* = 5.5193), and LSSVM‐ABC (*R*<sup>2</sup> = 0.9049, *RMSE* = 19.439 m, and *AARE* = 8.0032), respectively. In the *M*<sup>2</sup> (in testing phase), as the second‐best combination, the LSSVM‐WOA with respect to the highest *R*<sup>2</sup> (0.9896) and smallest *RMSE* (10.4268) outperformed the other models. The result assessment demonstrated that *M*<sup>6</sup> on account of poorest performance (*R*<sup>2</sup> = 0.3616 and *RMSE* = 58.5744 for the LSSVM‐WOA) was rec‐ ognized as the worst scheme regardless of the ML type. Typically, it can be understood from this remark that *S* is the most impacting input parameter on *FR* as its exclusion from the input variables (*M*6) caused the worst prediction performance regardless of the type of ML techniques, while density has a small effect on *FR* since its elimination from the input parameters (*M*2) did not significantly affect the degree of prediction accuracy. In

403.15 to 1847.43 and 35,479,174.56 to 67,688,321.04, respectively.

*5.2. Modelling Results*


**Table 2.** Descriptive statistics of all features used in modelling the *FR*.

As stated before, the database considered in this study was divided into training and testing sets. The first set is applied in the models' development, while the test set is devoted to the validation and investigation of the accuracy behavior of the established models when dealing with unseen measurements.

#### *5.2. Modelling Results*

By performing the steps described in the previous sections, it was found that 3 hidden layers with *tansig* as an activation function, and 12, 11, and 9 neurons in each of them, respectively, represented the proper CFNN topology in all of the six schemes. For LSSVM models, the achieved *σ* <sup>2</sup> and *γ* values using the ABC, GSA, and WOA ranged between 403.15 to 1847.43 and 35,479,174.56 to 67,688,321.04, respectively.

The statistical evaluation of the performance of the obtained ML-based models with respect to the stated six schemes in the previous sections is shown in Table 3. In this table, statistical criteria, namely *AARE*, *R* 2 , and *RMSE* are reported for the training set, the test set, and the whole dataset. According to this table, and based on the schemes, it can be seen that *M*<sup>1</sup> is the best one, followed by *M*2. In the combination of *M*<sup>1</sup> including all inputs for the testing phase, the LSSVM-WOA in terms of (*R* <sup>2</sup> = 0.999, *RMSE* = 3.4209 m, and *AARE* = 1.3017) was identified as the superior predictive model, followed by CFNN-LMA (*R* <sup>2</sup> = 0.9347, *RMSE* = 16.5215 m, and *AARE* = 7.512), LSSVM-GSA (*R* <sup>2</sup> = 0.904, *RMSE* = 16.1775 m, and *AARE* = 5.5193), and LSSVM-ABC (*R* <sup>2</sup> = 0.9049, *RMSE* = 19.439 m, and *AARE* = 8.0032), respectively. In the *M*<sup>2</sup> (in testing phase), as the second-best combination, the LSSVM-WOA with respect to the highest *R* 2 (0.9896) and smallest *RMSE* (10.4268) outperformed the other models. The result assessment demonstrated that *M*<sup>6</sup> on account of poorest performance (*R* <sup>2</sup> = 0.3616 and *RMSE* = 58.5744 for the LSSVM-WOA) was recognized as the worst scheme regardless of the ML type. Typically, it can be understood from this remark that *S* is the most impacting input parameter on *FR* as its exclusion from the input variables (*M*6) caused the worst prediction performance regardless of the type of ML techniques, while density has a small effect on *FR* since its elimination from the input parameters (*M*2) did not significantly affect the degree of prediction accuracy. In addition, it can also be deduced that for each of the six schemes, the LSSVM-WOA yielded more accurate predictions compared with the other LSSVM-metaheuristic algorithms and the CFNN-LMA. According to Table 3, it was found that *M*<sup>1</sup> outperformed other combinations followed by *M*2, *M*4, *M*3, *M*5, and *M*6, respectively. Additionally, it can be said that the LSSVM-WOA in all input combinations was the best predictive model developed in this study for prediction of *FR*. For better comparison between the predictive performances of the provided models in all combinations, the probability density function violin plots are exhibited in Figure 4. According to this figure, considering the best agreement between measured and predicted values of *FR*, it can be clearly implied that the LSSVM-WOA was the superior model for accurately estimating *FR*, the CFNN-LMA was identified as the second-best model, and LSSVM-ABC yielded the worst results in all combinations. Regarding the mentioned analysis, the combination of *M*<sup>1</sup> was kept for further performance investigation and validations. It is necessary to add that

in order to confirm that the ANN scheme did not suffer from the overfitting issue, a 4-fold cross-validation was performed on our best ANN paradigm (the case of *M*1) to assess the generalization of the model when dealing with new sets of data. To do so, the database was randomly divided into 4 folds, then, the modelling was done by considering a sole fold as the test sub-data and devoting the rest for the training phase. In order to swap between the folds involved in the training and testing phases, the aforesaid step was repeated 4 times. The results gained from the 4-fold cross-validation are reported in Table 4. As can be seen, the consistency of the model is confirmed for all the folds, thus, the overfitting issue is avoided.


**Table 3.** Performance of the suggested ML-based models with respect to the six schemes.


**Table 3.** *Cont.*

**Figure 4.** Probability density function of the predictive models of *M*<sup>1</sup> for prediction of *FR* for all datasets used in the simulation process. **Figure 4.** Probability density function of the predictive models of *M*<sup>1</sup> for prediction of *FR* for all datasets used in the simulation process.


In order to extend the examination of the accuracy of the best implemented models,

**Table 4.** Results of the 4‐fold cross‐validation for the CFNN‐LMA model (*M*<sup>1</sup> scheme). **Table 4.** Results of the 4-fold cross-validation for the CFNN-LMA model (*M*<sup>1</sup> scheme).

some graphical evaluation techniques were considered. Figure 5 shows the physical trend variation and cross plots related to the *M*<sup>1</sup> which illustrate a comparison of the measured and predicted *FR* values during both the training and testing phases. Based on the scatter plots, a tight cloud of points is located nearby the line *Y* = *X* for all datasets. This means that the LSSVM‐WOA can predict *FR* values with a great degree of accuracy as its predictions are very close to the perfect case shown by the unit‐slop line. The left side of Figure 6 shows the ability of the predictive models to capture the non‐linearity behaviour of the datasets. The LSSVM‐WOA yields the best agreement with the meas‐ ured *FR* compared with other LSSVM models and the CFNN‐LMA model. In order to extend the examination of the accuracy of the best implemented models, some graphical evaluation techniques were considered. Figure 5 shows the physical trend variation and cross plots related to the *M*<sup>1</sup> which illustrate a comparison of the measured and predicted *FR* values during both the training and testing phases. Based on the scatter plots, a tight cloud of points is located nearby the line *Y* = *X* for all datasets. This means that the LSSVM-WOA can predict *FR* values with a great degree of accuracy as its predictions are very close to the perfect case shown by the unit-slop line. The left side of Figure 6 shows the ability of the predictive models to capture the non-linearity behaviour of the datasets. The LSSVM-WOA yields the best agreement with the measured *FR* compared with other LSSVM models and the CFNN-LMA model.

Another visual tool for inspecting the reliability of the LSSVM-WOA paradigm is dealing with the relative deviation (RD%) distribution diagram that is exhibited in Figure 6. The quartiles of RD for 25% and 75% of datasets are listed in Table 5. The LSSVM-WOA model, owing to having the least values of Q25% = −0.5566, Q75% = 1.093, and IQR = 1.650, yielded more reliable outcomes in comparison with the CFNN-LMA (Q25% = −2.976, Q75% = 2.105, and IQR = 5.081), LSSVM-ABC (Q25% = −2.778, Q75% = 2.504, and IQR = 5.282), and LSSVM-GSA (Q25% = −2.712, Q75% = 2.716, and IQR = 5.428), respectively. Obviously, the above

diagnostic analysis reveals that the LSSVM-WOA model has very low prediction errors regardless of the considered conditions. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 15 of 20

**Figure 5.** Comparing the predicted and measured *FR* values in the form of the cross plot and trend variation plot for all datasets. **Figure 5.** Comparing the predicted and measured *FR* values in the form of the cross plot and trend variation plot for all datasets.

Another visual tool for inspecting the reliability of the LSSVM‐WOA paradigm is dealing with the relative deviation (RD%) distribution diagram that is exhibited in Fig‐ ure 6. The quartiles of RD for 25% and 75% of datasets are listed in Table 5. The LSSVM‐WOA model, owing to having the least values of Q25% = −0.5566, Q75% = 1.093, and IQR = 1.650, yielded more reliable outcomes in comparison with the CFNN‐LMA (Q25% = −2.976, Q75% = 2.105, and IQR = 5.081), LSSVM‐ABC (Q25% = −2.778, Q75% = tions.

tions.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 16 of 20

**Figure 6.** Relative errors of prediction models versus the measured *FR* values. **Figure 6.** Relative errors of prediction models versus the measured *FR* values. **Figure 6.** Relative errors of prediction models versus the measured *FR* values.



5.428), respectively. Obviously, the above diagnostic analysis reveals that the LSSVM‐WOA model has very low prediction errors regardless of the considered condi‐

5.428), respectively. Obviously, the above diagnostic analysis reveals that the LSSVM‐WOA model has very low prediction errors regardless of the considered condi‐

LSSVM‐WOA versus the *S* and *PF*, as the two most significant variables, to assess the re‐ liability seeking zone for estimating the *FR*. The RD variation reveals that the most relia‐ ble and accurate results are obtained in the ranges of 3 ≤ *S* ≤ 3.5 and 0.8 ≤ *PF* ≤1.1. Lastly, Figure 7 demonstrates the 3D surface of RD variation concerning the LSSVM-WOA versus the *S* and *PF*, as the two most significant variables, to assess the reliability seeking zone for estimating the *FR*. The RD variation reveals that the most reliable and accurate results are obtained in the ranges of 3 ≤ *S* ≤ 3.5 and 0.8 ≤ *PF* ≤1.1. Lastly, Figure 7 demonstrates the 3D surface of RD variation concerning the LSSVM‐WOA versus the *S* and *PF*, as the two most significant variables, to assess the re‐ liability seeking zone for estimating the *FR*. The RD variation reveals that the most relia‐ ble and accurate results are obtained in the ranges of 3 ≤ *S* ≤ 3.5 and 0.8 ≤ *PF* ≤1.1.

Lastly, Figure 7 demonstrates the 3D surface of RD variation concerning the

**Figure 7.** RD variation versus S and PF for LSSVM‐WOA model. (Red points: Samples in distribu‐ tion). **Figure 7.** RD variation versus S and PF for LSSVM‐WOA model. (Red points: Samples in distribu‐ tion). **Figure 7.** RD variation versus S and PF for LSSVM-WOA model. (Red points: Samples in distribution).

#### **6. Conclusions**

FR is one of the most adverse effects induced by blasting in surface mines. In this study, an attempt was made to predict blast-induced *FR* through hybridizing three optimization

algorithms including the ABC, GSA, and WOA with the LSSVM model. In addition, CFNN-LMA, as a powerful tool for the prediction aims, was developed. For developing the models, six different schemes based on different combinations of input parameters were employed and in total, twenty-four different models were constructed. After that, three statistical indexes, i.e., *AARE*, *R* 2 , and *RMSE*, were used to check the performance of the models and to compare their results.

Some conclusions are drawn as follows: the results indicated that among the total constructed models, the LSSVM-WOA model was the most accurate model in all six schemes compared with the CFNN-LMA, LSSVM-GSA, and LSSVM-ABC models. The most accurate results of the LSSVM-WOA (*AARE* = 1.3017, *R* <sup>2</sup> = 0.9991 and *RMSE* = 3.4209), LSSVM-GSA (*AARE* = 5.5193, *R* <sup>2</sup> = 0.904 and *RMSE* = 16.1775), and CFNN-LMA (*AARE* = 7.512, *R* <sup>2</sup> = 0.9347 and *RMSE* = 16.5215) were obtained from the first scheme, while that for the LSSVM-ABC model (*AARE* = 7.083, *R* <sup>2</sup> = 0.9235 and *RMSE* = 16.7483) was obtained from the second scheme. It is important to note that the above results were related to the testing phase. On the other hand, the sixth scheme had the worst performance. In this scheme, the *S* parameter was removed from the modelling. Therefore, it can be suggested that the *S* was an effective parameter in the modelling. Additionally, according to RD%, the LSSVM-WOA model had very low prediction errors regardless of the considered conditions. The presented results in this study cannot be compared with results of the previous studies because different fields investigation as well as different range of input parameters were used in the previous studies. Nevertheless, for a comparison with the literature, the LSSVM-WOA presented in this study predicted the *FR* with a very good *R* 2 , while Koopialipoor et al. [55], Faradonbeh et al. [56], Zhou et al. [57], Nguyen et al. [58], and Marto et al. [59] predicted the *FR* with an *R* <sup>2</sup> of 0.959, 0.924, 0.944, 0.986, and 0.981, respectively. The aforementioned results indicate the effectiveness of the LSSVM-WOA model in predicting the *FR*.

It is worth mentioning that the proposed models herein are specific to studied cases and the use of these models in other surface mines requires some modification based on blasting and mining conditions. However, our best implemented model can predict the *FR* with high accuracy and the paradigm can be applied for cases filling the applicability conditions considered in this work, namely by respecting the maximum and minimum values of the involved input parameters.

For future works, it can be recommended to use other optimization algorithms, such as seagull optimization algorithm (SOA), sparrow search algorithm (SSA), loin swarm optimization (LSO), and moth-flame optimization (MFO) algorithms, in combination with CFNN and LSSVM models.

**Author Contributions:** Conceptualization, M.H.; Methodology, X.D. and M.J.; Validation, X.D., M.J. and M.H.; Investigation, M.H. and B.N.L.; Data curation, M.H.; writing—original draft preparation, X.D., M.J., M.H., R.A.A. and B.N.L.; writing—review and editing, X.D., M.J., M.H., R.A.A. and B.N.L.; Supervision, M.H. and R.A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper is supported by the National Natural Science Foundation of China (Grant No. 52174131).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in this study may be available on request from the corresponding author.

**Acknowledgments:** The authors would like to thank Danial Jahed Armaghani for providing the information and facilities required for conducting this research.

**Conflicts of Interest:** The authors have no conflict of interest to declare that are relevant to the content of this article.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
