*Article* **Prediction of Blast-Induced Ground Vibration at a Limestone Quarry: An Artificial Intelligence Approach**

**Clement Kweku Arthur 1, Ramesh Murlidhar Bhatawdekar 2,\*, Edy Tonnizam Mohamad 2, Mohanad Muayad Sabri Sabri 3, Manish Bohra 4, Manoj Khandelwal <sup>5</sup> and Sangki Kwon 6,\***


**Abstract:** Ground vibration is one of the most unfavourable environmental effects of blasting activities, which can cause serious damage to neighboring homes and structures. As a result, effective forecasting of their severity is critical to controlling and reducing their recurrence. There are several conventional vibration predictor equations available proposed by different researchers but most of them are based on only two parameters, i.e., explosive charge used per delay and distance between blast face to the monitoring point. It is a well-known fact that blasting results are influenced by a number of blast design parameters, such as burden, spacing, powder factor, etc. but these are not being considered in any of the available conventional predictors and due to that they show a high error in predicting blast vibrations. Nowadays, artificial intelligence has been widely used in blast engineering. Thus, three artificial intelligence approaches, namely Gaussian process regression (GPR), extreme learning machine (ELM) and backpropagation neural network (BPNN) were used in this study to estimate ground vibration caused by blasting in Shree Cement Ras Limestone Mine in India. To achieve that aim, 101 blasting datasets with powder factor, average depth, distance, spacing, burden, charge weight, and stemming length as input parameters were collected from the mine site. For comparison purposes, a simple multivariate regression analysis (MVRA) model as well as, a nonparametric regression-based technique known as multivariate adaptive regression splines (MARS) was also constructed using the same datasets. This study serves as a foundational study for the comparison of GPR, BPNN, ELM, MARS and MVRA to ascertain their respective predictive performances. Eighty-one (81) datasets representing 80% of the total blasting datasets were used to construct and train the various predictive models while 20 data samples (20%) were utilized for evaluating the predictive capabilities of the developed predictive models. Using the testing datasets, major indicators of performance, namely mean squared error (MSE), variance accounted for (VAF), correlation coefficient (*R*) and coefficient of determination (*R*2) were compared as statistical evaluators of model performance. This study revealed that the GPR model exhibited superior predictive capability in comparison to the MARS, BPNN, ELM and MVRA. The GPR model showed the highest VAF, *R* and *R*<sup>2</sup> values of 99.1728%, 0.9985 and 0.9971 respectively and the lowest MSE of 0.0903. As a result, the blast engineer can employ GPR as an effective and appropriate method for forecasting blast-induced ground vibration.

**Keywords:** artificial intelligence; backpropagation neural network; blast-induced ground vibration; Gaussian process regression

**Citation:** Arthur, C.K.; Bhatawdekar, R.M.; Mohamad, E.T.; Sabri, M.M.S.; Bohra, M.; Khandelwal, M.; Kwon, S. Prediction of Blast-Induced Ground Vibration at a Limestone Quarry: An Artificial Intelligence Approach. *Appl. Sci.* **2022**, *12*, 9189. https://doi.org/ 10.3390/app12189189

Academic Editor: José António Correia

Received: 28 July 2022 Accepted: 26 August 2022 Published: 14 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 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/).

#### **1. Introduction**

Ground vibration is one of the main adverse blasting outcomes that has received significant attention in the mining and civil industries [1,2]. Ground vibration is known to have a lot of adverse impacts on the environment (cracks on building structures) and the stability of pit walls. It is worth mentioning that several factors contribute to the occurrence of these blast-induced ground vibrations. These factors can be categorized into controllable factors and uncontrollable factors [3,4]. The controllable factors are those that the blast engineer has control over and can change. These include the blast design parameters of stemming length, hole depth, spacing, burden, hole inclination and explosive parameters of delay timings, a maximum charge per delay, and total charge. The uncontrollable factors are those the blast engineer has no control over, and they include both geotechnical and geomechanical parameters such as rock strength, faults, and folds [5–9]. The peak particle velocity (PPV) is the index for assessing ground vibration induced by blasting [10]. When detonation of explosives takes place, high energy is released in the blast hole which fractures the rock surrounding the blasthole [11]. Some of the energy released is used to fragment and displace the rock mass. The rest of the energymove through the ground as ground vibration andimpacts surrounding structures.

Due to the adverse impact of blast-induced ground vibration, it has always been in the interest of the blast engineer to model and predicts its occurrence to minimize vibration level as much as possible. In that regard, a lot of research has been conducted since the 1950s [12] to develop models for predicting ground vibration arising out of blasting operations. These models have been developed using empirical techniques through to the use of artificial intelligence (AI) techniques [13]. These AI techniques have been found to produce more accurate results than the empirical techniques and hence have received worldwide attention due to their unique capabilities [14]. AI techniques that have been developed and used in the prediction of blasting outcomes (ground vibration, air overpressure, and flyrock) are outlined in Table 1. It is worth noting that all abbreviations used in this work are presented in the Abbreviations Section.


**Table 1.** AI Models developed and applied to predict ground vibration, air overpressure and flyrock.

More recently in ground vibration studies, other researchers have applied evolutionary and metaheuristic optimization algorithms to optimize simple AI techniques. Some of these works are presented in Table 2.

**Table 2.** Hybrid Models developed and applied to predict ground vibrations.


Table 3 provides a detailed summary of some research on ground vibration prediction.


**Table 3.** Input parameters, size of data and AI techniques for prediction of ground vibration.

Nevertheless, the application of single AI techniques is still of interest in this evergrowing technological world. ANN has been developed by [69] to predict the earth surface deformation. Thus, the predictive capacities of three artificial intelligence algorithms, backpropagation neural network (BPNN), ELM, and GPR, are investigated in this study using blasting data from a quarry (Ras Limestone Mine of Shree Cement) in India to estimate PPV values. A multivariate adaptive regression spline (MARS) approach, as well as a multivariate regression analysis (MVRA) model, was developed and used for comparison purposes. Studies have been made to compare the GPR and BPNN [22], MARS and BPNN [67], ELM and BPNN [70], GP and MARS [63], GPR and MVRA [71] and BPNN and MVRA [55]. However, little has been done in the literature to compare the predictive performance of GPR, MARS, BPNN, ELM and MVRA in ground vibration prediction studies. In that regard, this study is exploratory. It is worth mentioning that the empirical models developed for predicting blast-induced ground vibration were not considered

in this study. The reason is that studies done by [17,40,53,57,59,72,73] have proved that these empirical models do not produce accurate results. The models used in this study consider seven effective parameters, namely the average depth, a maximum charge per delay, powder factor, spacing, burden, distance and stemming length, because, as shown in [5–7], they significantly affect the intensity of ground vibration.

#### **2. Study Site and Data Description**

The Ras Limestone Mine of Shree Cement is located 30 km from Beawar City, Ajmer District, Rajasthan, India. The mining concession of 750.0 ha lies between longitude E 74◦10 5.96 to E 74◦11 9.62 and latitudes N 26◦16 57.13 to N 26◦15 36.23, on toposheet No. 45 J/3 & 45 J/4 of the survey of India.

The projected production capacity of the mine is 25.3 million tons of limestone per year. The mining area is generally rocky with no overburden. A general strike of limestone at Ras Mine is North-South direction and dips in the eastern direction. Limestone has four major folds and one reverse fault. Limestone strata are massive, blocky and fractured in different portions of the deposit. HRB 150 (INDUS Make) drills are used for drilling hole diameter of 165 mm. ANFO with cast booster/slurry explosives and nonel detonators are used as explosives for blasting limestone. Figure 1 shows a blasting round view with Figure 2 showing the close-up view of blasted limestone at Shree Cement Ras Limestone Mine in India.

**Figure 1.** Blasting Round View.

**Figure 2.** Close-up View of Blasted Limestone at Shree Cement Ras Limestone Mine in India.

As a part of this study, for the establishment of the various models described therein, a total of 101 sets of data were collected from the Ras limestone mine. The data collected consisted of parameters such as average depth (m), spacing (m), burden (m), powder factor (t/kg), the distance between the blasting point and the monitoring station (m), stemming length (m), a maximum charge per delay (kg) and PPV (mm/s). In the creation of the various models, the input parameters were average depth (m), spacing (m), burden (m), powder factor (t/kg), the distance between the blasting site and the monitoring station (m), stemming length (m), and maximum charge per delay (kg), while the output parameter was PPV. Table 4 shows the statistical description of the dataset collected.


**Table 4.** Description of dataset parameters.

The values for the maximum charge per delay, stemming length, powder factor, spacing, burden, and average depth as statistically described in Table 2 were obtained from the daily blast plans of the mine. The distance values were calculated using the coordinates of the blasting face and monitoring locations obtained using a Global Positioning System (GPS). As shown in Figure 3, the PPV values were monitored using an Instantel Micromate ISEE Std/XM seismograph [74].

**Figure 3.** Instantel Micromate ISEE Std/XM seismograph.

It is worth mentioning that the mine has no permanent monitoring location due to different blasting positions. Thus, in monitoring the ground vibration due to blasting, the seismograph is positioned using pegs with an arrow on the geophone pointing towards the blast site. Figure 4 shows the portable monitoring station used by the mine. It is worth noting that the terrain of the Ras Limestone Mine is generally hilly.

The correlation coefficient matrix shows how strong the interaction between the input parameters (average depth, burden, spacing, distance, powder factor, stemming length, and maximum charge per delay) and the measured PPV is, as shown in Table 5.

**Figure 4.** Portable ground vibration monitoring station in realistic conditions (at Shree Cement Ras Limestone Mine in India).

**Table 5.** Matrix of Correlation Coefficients Between Input Parameters and PPV Measured.


#### **3. Methodology**

In this section, the mathematical description of the different methods applied in this study will be briefly outlined. Furthermore, the procedure followed to develop the various models as well as the models' performance indicators will be outlined.

#### *3.1. Study Steps*

A systematic methodology was utilized in this study. First, the data collected were prepared by removing all outliers and then were partitioned into two sets (training set and testing set) and normalized into the interval [–1,1]. Then the various models were built by selecting the model's hyperparameter. The models were then trained using the training dataset. Finally, the model's results were assessed based on the test dataset by some performance indicators. The performance results were then analyzed to either finetune the model's hyperparameter or select the model as optimum. Figure 5 shows the flowchart applied in this study.

**Figure 5.** A Systematic Flowchart for Prediction of Blast-Induced Ground Vibration.

*3.2. Mathematical Description of the Different Methods*

3.2.1. GPR

Gaussian Process (GP)

GP is a nonparametric Bayesian technique that is used in regression modelling [75]. This GP process can be described as a finite assemblage of a set of arbitrary parameters that follow a multivariate Gaussian (normal) distribution [76]. That is, for every given input point from a set of input vectors *r* = (*r*1,*r*2,*r*3,...,*rm*), the probability distribution over its function *h*(*r*) follows a Gaussian distribution. Thus, a GP *h*(*r*) is precisely shown in Equation (1) as:

$$h(r) \sim \text{GP}\left(b(r), \lg(r, r')\right) \tag{1}$$

From Equation (1) it can be deduced that a GP is fully characterized by a covariance function *g*(*r*,*r* ) and a mean function (MF) *b*(*r*) as expressed in Equation (2).

$$\begin{cases} \begin{array}{l} b(r) = E[h(r)] \\ g(r, r') = E[(b(r) - h(r))(b(r') - h(r'))] \end{array} \tag{2} \end{cases} \tag{2}$$

For the basic GPR, the MF is normally set as 0, however, there many other MFs which can be applied in building the GPR model [77]. The noted MFs in literature have been categorized into two kinds, namely: simple and composite. The simple MFs include zero, one, constant, linear, polynomial, nearest neighbor MFs etc. whereas the composite ones include: the scaled version, sum, product, power and warped MFs [77]. It is worth noting that this study adopted an MF with a constant, *b*.

The covariance function on the other hand is the main component in the development of the GPR model. The best covariance function is dependent on the data being modelled. Literature is replete with a number of these covariance functions [70]. However, the notable ones include: the rational quadratic, matérn class, squared exponential and the exponential

covariance functions. The most often used covariance function is the squared exponential covariance function [77,78].

#### Prediction Using GP

In the case of a regression modelling problem, an output variable *q* can be approximated, given function *h*(*r*) with an additive noise *ε<sup>i</sup>* component inherent in the dataset as shown in Equation (3).

$$
\eta\_i = h(r\_i) + \varepsilon\_i \tag{3}
$$

Assuming this noise component *ε<sup>i</sup>* has a zero mean and variance *σ*<sup>2</sup> *<sup>n</sup>*, the prior on the noisy data is expressed in Equation (4) as:

$$\text{cov}(q) = \text{g}\left(r, r'\right) + \sigma\_n^2 I\_n \tag{4}$$

where *In* is a matrix of the n-dimensional unit.

The GP *h*(*r*) (see Equation (1)) is then precisely considered in Equation (5) as:

$$h(r) \sim \text{GP}\left(b(r), \lg(r, r') + \sigma\_n^2 I\right) \tag{5}$$

It should be emphasized that the GP model training, seeks to ascertain the best possible hyperparameter set *Θ* = *β*, *χ*, *υ*<sup>2</sup> *<sup>s</sup>*, *σ*<sup>2</sup> *n* that best fits the data sets. This can be done by the use of a maximum possible method [69] in which the log-likelihood function is maximized (Equation (6)).

$$\log\left(p(q|r,\Theta)\right) = \frac{1}{2}\log\left(\det\left(\mathbf{g}\left(r,r'\right) + \sigma\_n^2 I\right)\right) - \frac{1}{2}q^T\left(\mathbf{g}\left(r,r'\right) + \sigma\_n^2 I\right)^{-1}q - \frac{n}{2}\log 2\pi \tag{6}$$

Of all the maximum likelihood functions available, the conjugate gradient method is the most widely used [79] and hence was used in this study. It finds the optimal hyperparameter sets by using the partial differential of the log-likelihood function (Equation (6)) in relation to the hyperparameter set, *Θ* as shown in Equation (7).

$$\begin{split} \frac{\partial}{\partial \Theta\_i} \log \left( p(q|r, \Theta) \right) &= \frac{1}{2} q^T G^{-1} \frac{\partial G}{\partial \Theta\_i} G^{-1} q - \frac{1}{2} tr \left( G^{-1} \frac{\partial G}{\partial \Theta\_i} \right) \\ &= \frac{1}{2} tr \left( \left( \beta \beta^T - G^{-1} \right) \frac{\partial G}{\partial \Theta\_i} \right) \end{split} \tag{7}$$

where *β* = *G*−1*q* and *G* = *g*(*r*,*r* ).

Given the joint prior distribution of the training output variable, *q* at point a and the value *q*<sup>∗</sup> to be predicted at the test point *r*<sup>∗</sup> expressed in Equation (8), the GPR model is able to predict *q*<sup>∗</sup> by calculating the posterior distribution *p*(*q*∗|*r*, *q*,*r*<sup>∗</sup> ) (Equation (9)).

$$
\begin{bmatrix} q \\ q\_\* \end{bmatrix} \sim GP \left( \begin{bmatrix} b(r) \\ b(r\_\*) \end{bmatrix}, \begin{bmatrix} \mathbf{g}(r,r) + \sigma\_n^2 I & \mathbf{g}(r,r\_\*) \\ \mathbf{g}(r\_\*,r) & \mathbf{g}(r\_\*,r\_\*) \end{bmatrix} \right), \tag{8}
$$

$$\not{p}(q\_\*|r,q,r\_\*) \sim GP(\overline{q}\_\*,cov(q\_\*)),\tag{9}$$

Here *<sup>q</sup>*<sup>∗</sup> (Equation (10)) is the mean value which is the estimation of *<sup>q</sup>*<sup>∗</sup> and *cov*(*q*∗) (Equation (11)) is the predictive variance matrix of the test data, which reveals the credibility of the prediction values [79].

$$\overline{q}\_{\*} = b(r\_{\*}) + \mathcal{g}(r\_{\*}, r) \left[ \mathcal{g}(r, r) + \sigma\_{n}^{2} I \right]^{-1} (q - b(r)) \tag{10}$$

$$\text{cov}(q\_{\ast}) = \text{g}(r\_{\ast}, r\_{\ast}) \left[ \text{g}(r, r) + \sigma\_{\text{n}}^{2} I \right]^{-1} \text{g}(r, r\_{\ast}) \tag{11}$$

#### 3.2.2. BPNN

BPNN is a widely used AI technique that was developed to mimic the human brain. In this, there is an input layer that takes impulses from the outside environment as inputs to the network. These inputs *xk* are weighted by connecting weights *wk* and relayed to the hidden layer. The hidden layer contains processing units called neurons which transform the weighted input by a transfer function, *t*. It is noteworthy that biases *b* are added to the transfer function before the transformation process. The hidden layer's output is subsequently conveyed to the output layer, which is transformed by a transfer function operating inside the hidden layer. The network's predicted values are then derived from the output, *y*ˆ from the output layer as shown in Equation (12).

$$\mathfrak{F} = t\left(\sum\_{k=1}^{m} w\_k x\_k + b\right) \tag{12}$$

In training the BPNN, a training algorithm is used in updating weights and biases based on the backpropagation error, *e* (divergence in true and predicted value) as shown in Equation (13) so as to produce a network with a minimum propagation error.

$$e = y - \hat{y}\tag{13}$$

Several training algorithms have been developed for such purposes. However, the Levenberg–Marquardt algorithm [80] is the widely used training function due to its high convergence speed and accuracy and thus was used in this study.

#### 3.2.3. MARS

The MARS algorithm is a non-parametric algorithm developed by [81] to estimate the complex nonlinear correlation between model inputs and output. This estimating process is achieved by automatically building a series of linear piecewise regression models through the use of basis functions, to fit the given data pair.

In the general, the MARS model is of the form precisely considered in Equation (14):

$$\hat{f}(z) = \beta\_0 + \sum\_{k=1}^{N} \beta\_k \lambda\_k(z) \tag{14}$$

where ˆ *f*(*z*) signifies the estimated output parameter value, *β*<sup>0</sup> is constant, *λk*(*z*) is the *kth* basis function, *β<sup>k</sup>* signifies the *k*th basis function's coefficient and *z* signifies the input variable. The basis function act as a hinge function to split the data into separate sections, which can be modelled individually. Each basis function can be precisely considered in Equation (15) as:

$$\lambda\_k(z) = \prod\_{i=1}^{I\_k} \left[ s\_{ik} \cdot \left( z\_{v(i,k)} - h\_{ik} \right) \right]\_+ \tag{15}$$

where *Ik* is the quantity of splits that formed *λk*(*z*), *sik* is the selected sign with value ±, *v*(*i*, *k*) labels the predictor variable and *hik* is the knot value on the corresponding input variables.

The MARS algorithm adopts two main steps namely: the forward selection process and the backward deletion process; to develop its model. In the forward selection process, the model is initially constructed with a constant basis function. New pairs of basis functions are thereafter iteratively included in the model to reduce the training residual sum-ofsquares error; to improve the model. However, as many basis functions are added in the forward process; the model built becomes overfit and cannot generalize well with unseen data. The backward deletion process is then introduced to remove all redundant basis functions. It employs the generalized cross-validation (GCV) Equation (16) to evaluate the performance of individually created models as it eliminates the unwanted basis functions. The individually created model with the least value of GCV is then chosen as the optimal MARS model.

$$\text{GCV}(Q) = \frac{\frac{1}{H} \sum\_{j=1}^{H} \left( y\_j - \hat{f}\_Q(z\_j) \right)^2}{\left( 1 - \frac{C(Q)}{H} \right)^2} \tag{16}$$

where *yj* and ˆ *fQ zj* denotes the actual output and predicted values of the training samples, and *H* represents the total number of training samples. As shown in Equation (17), *C*(*Q*) is a penalty for model complexity that is proportional to the model's number of basis functions.

$$C(Q) = (Q+1) + pQ \tag{17}$$

where *p* is the penalty cost for the optimization of every single basis function which works as a smoothing variable. The details of MARS as well as the selection of the *p* are in [76].

#### 3.2.4. MVRA

MVRA is a statistical tool applied to fit a model to establish a linear relation between a set of input parameters (independent variables) and an output parameter (dependent variable) [82]. This fitted model can then be used to make predictions on new data. MVRA works by studying the correlation between the various input parameters and output parameters to construct simultaneous equations so as to acquire the best-fit equation. It uses an ordinary least squares fit on the dataset to find the best-fit equation. It forms a regression matrix in the process of solving simultaneous equations. The regression matrix is then solved using the backslash operator to obtain the regression coefficient as well as the intercept [83]. Generally, the MVRA is mathematically expressed in Equation (18) as:

$$\mathcal{Y} = \beta\_0 + \beta\_2 X\_2 + \beta\_3 X\_3 + \dots + \beta\_k X\_k \tag{18}$$

where *β*1, ..., *β<sup>k</sup>* are the regression coefficients, *β*<sup>0</sup> is the intercept *X*1, *X*2, ..., *Xk* is the independent variable and *Y* is the dependent variable.

#### 3.2.5. ELM

In 2004 Huang introduces the mathematical model of ELM. The ELM's basic principle is based on a single hidden layer feed-forward neural network (SLFN) (Figure 6). Because of its improved generality, simplicity, and efficient forecasting nature, the ELM has been employed in a variety of application areas [84].

**Figure 6.** ELM Architecture.

The basic premise of ELM is as follows: Given *N* as the number of hidden units, *K* as the number of training samples, and the activation function *f*( ) in the hidden units, the output of the ELM *om* for the *m*th training sample is depicted in Equation (19) as:

$$\rho\_{\mathfrak{m}} = \sum\_{i=1}^{N} \beta\_i f(w\_{k'} b\_{i'} \mathfrak{x}\_{\mathfrak{m}}) \ m = 1, \ldots, K \tag{19}$$

where *bi* is the hidden neurons' bias factor, *xm* denotes the number of inputs, *β<sup>i</sup>* denotes the output weight vectors and *wk* denotes input weight vectors. The sigmoid function is used as an activation function. The sigmoid function's output is essentially in a range of 1 to 0. To determine the output weights, the linear equation (Equation (20)) is employed.

$$\mathcal{S} = H^\dagger Y \tag{20}$$

where *H* denotes the output matrix of the hidden layer, *H*† the Moore–Penrose generalized inverse [85] of *H*, and *Y* denotes the ELM output targets. In Equation (21), Equation (20) is written as:

$$H\beta = \text{Y} \tag{21}$$

Equation (22) can be used to define *H*, *β*, *Y* as follows:

$$H = \begin{bmatrix} p(\mathbf{x}\_1) \\ \vdots \\ p(\mathbf{x}\_{N}) \end{bmatrix} = \begin{bmatrix} f(w\_1, b\_1, \mathbf{x}\_1) & \cdots & f(w\_N, b\_N, \mathbf{x}\_1) \\ \vdots & \cdots & \vdots \\ f(w\_1, b\_1, \mathbf{x}\_j) & \cdots & f(w\_N, b\_N, \mathbf{x}\_j) \end{bmatrix}\_{N \times K}, \\ \boldsymbol{\beta} = \begin{bmatrix} \beta\_1^Y \\ \vdots \\ \beta\_N^Y \end{bmatrix} \text{and } \mathbf{Y} = \begin{bmatrix} y\_1^Y \\ \vdots \\ y\_K^Y \end{bmatrix} \tag{22}$$

In this case, the hidden layer's feature mapping is *p*(*x*). *H* is the ELM's output.

#### *3.3. Procedures for Model Construction*

#### 3.3.1. Data Selection and Division

In modelling the various approaches presented in this study, the hold-out crossvalidation technique was employed to partition the entire 101 datasets. The datasets were split into the 80:20 ratio. The first 80% of the total datasets were used as the training set (representing 81 training datasets). The remaining 20% (representing 20 datasets) were used as the test set. This strategy was adopted because [86,87] have proved that a ratio of 80:20 or 70:30 will produce accurate prediction results and will not cause overfitting.

#### 3.3.2. Data Normalization

In the data preparation phase, it is expedient that the input parameters be normalized. This is because the input parameters have different input ranges order and those with the higher values have the potential to skew the prediction results to themselves. Thus, to avoid this predicament and give equal chances to each input parameter to influence the prediction outcome, the input parameters defined in Table 1 were normalized into the interval [–1,1] [88,89] utilizing Equation (23).

$$F\_l = F\_{\rm min} + \frac{(E\_l - E\_{\rm min}) \times (F\_{\rm max} - F\_{\rm min})}{(E\_{\rm max} - E\_{\rm min})} \tag{2.3}$$

where *Ei* signifies the actual data, *Emax* and *E*min refer to the maximum values and minimum of the actual data, *Fi* are the normalized data and *F*min and *Fmax* being the min-max values of −1 and 1 in that order.

#### 3.3.3. Model Development

For the development of the GPR model, five different models based on the squared exponential (Equation (24)), exponential (Equation (25)), rational quadratic (Equation (26)), matérn 3/2 (Equation (27)), and matérn 5/2 (Equation (28)), covariance function as well as the functions were developed. Each model had a constant MF.

$$\log\left(r, r'\right) = v\_s^2 \exp\left[\frac{-||r - r'||}{2\chi^2}\right] \tag{24}$$

$$\log\left(r, r'\right) = \upsilon\_s^2 \exp\left[\frac{-||r - r'||}{\chi}\right] \tag{25}$$

$$\mathcal{g}\left(r,r'\right) = \upsilon\_s^2 \exp\left[\frac{-||r-r'||}{2\beta\chi^2}\right]^{-\beta} \tag{26}$$

$$\log(r, r') = v\_s^2 \exp\left[1 + \frac{\sqrt{3}||r - r'||}{\chi}\right] \exp\left[-\frac{\sqrt{3}||r - r'||}{\chi}\right] \tag{27}$$

$$\log(r, r') = v\_s^2 \exp\left[1 + \frac{\sqrt{5}||r - r'||}{\chi} + \frac{5||r - r'||^2}{3\chi^2}\right] \exp\left[-\frac{\sqrt{5}||r - r'||}{\chi}\right] \tag{28}$$

where *β* is the rational quadratic covariance's shape parameter, *χ* is the length scale, and *υ*<sup>2</sup> *s* is the covariance function's signal variance.

The model with the lowest mean squared error and highest correlation coefficient on the test dataset was chosen as the optimum GPR model. For the BPNN model, a threelayered architecture was chosen—the first with the input layer, the second with a hidden layer and the thirdly with an output layer. A single hidden layer was used because it has been established to be a reliable predictor for any prediction problem [90]. Furthermore, in the case of hidden and output layers, hyperbolic and linear transfer functions were selected and used. The Levenberg–Marquardt algorithm was used to train this BPNN model. According to the suggested values by the previous researchers, a range of 1 to 40 for neurons was tried and the optimum number was the one that gives the lowest MSE on the test dataset [91,92]. The optimum number of neurons in the hidden layer that resulted in the lowest MSE on the test dataset was determined using a sequential experimental procedure in the construction of the ELM model. In that regard, 1 to 20 neurons were tried. It is worth stating that, the building of the MARS model, entails the choice of the highest number of basis functions to be used in the forward selection stage as well as the maximum degree of interaction. These serve as constraints in the development process. Based on their levels of interaction, three independent MARS models were built in this study–zero-degree, first degree and second-degree. Furthermore, a maximum of 20 basis functions were selected for the forward selection stage. The model with the highest correlation coefficient and lowest mean squared error (MSE) was chosen as the optimum MARS model. The MVRA model was developed using the same dataset for the development and testing of the GPR, BPNN, ELM and MARS models. The MVRA solves the multilinear regression equations established for the various input parameters and PPV using the least square technique in order to find the regression coefficient (Equation (18)) for each input parameter as well as the intercept.

#### 3.3.4. Performance Indicators

The performance of the various models constructed in this study was assessed using performance measures such as variance accounted for (VAF), correlation coefficient (*R*), coefficient of determination (*R*2) and mean squared error (MSE). These indicators are precisely shown in Equations (29)–(32) as:

$$\text{MSE} = \frac{1}{p} \left[ \sum\_{i=1}^{p} (s\_i - q\_i)^2 \right] \tag{29}$$

$$R = \frac{\sum\_{i=1}^{p} \left(s\_i - \overline{s}\right) \left(q\_i - \overline{q}\right)}{\sqrt{\sum\_{i=1}^{p} \left(s\_i - \overline{s}\right)^2} \times \sqrt{\sum\_{i=1}^{p} \left(q\_i - \overline{q}\right)^2}}\tag{30}$$

*R*<sup>2</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ *p* ∑ *i*=1 (*si* − *s*)(*qi* − *q*) *p* ∑ *i*=1 (*si* − *s*) 2 × *p* ∑ *i*=1 (*qi* − *q*) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ 2 (31)

$$VAF = \left(1 - \frac{var(s\_i - q\_i)}{var(s\_i)}\right) \tag{32}$$

where *q* represents the mean of the estimated values, *qi* represents the estimated values, *si* represents the measured values, *p* is the number of observations, while *s* denotes the average of the measured values.

#### **4. Results and Discussion**

*4.1. Developed Models*

4.1.1. Gaussian Process Regression

As shown in Table 6, the optimum GPR model that produced the MSE of 0.0903 and the highest *R*-value of 0.9986 for the testing dataset, had a matérn 3/2 covariance function with a noise variance of 0.06434, a length scale of 3.6019, and a signal variance of 7.0339. This indicates that the GPR-matérn 3/2 can generalize well with unseen datasets relative to the other GPR models. Hence, GPR-matérn 3/2 model was selected as the best GPR model in this study.

**Table 6.** Results of the Five different GPR Models.


#### 4.1.2. BPNN

As shown in Table 7, the optimal BPNN model has one neuron in the hidden layered network. Thus, having an architecture [7-1-1] which means seven input parameters and one neuron in the hidden layer, and an output layer. This is because it has the lowest MSE value on test datasets.

**Table 7.** Results of BPNN for Different Architectures.



**Table 7.** *Cont.*

#### 4.1.3. MARS

As shown in Table 8, the developed MARS model with the first order of interaction had the highest R values as well as the lowest MSE values on both the training and test datasets. Hence it was chosen as the optimum MARS model in this study.

**Table 8.** Results of Different MARS Models.


In the developmental process of the selected first order of interaction MARS model, only eight basis functions after the backward elimination stage were used out of the 20 basis functions employed in the forward selection stage. The eight basis functions of the selected MARS model and their respective equations are shown in Table 9.

**Table 9.** The Relationship Between Basis Functions and their Related Equations.


The developed optimum MARS model for predicting ground vibration as a result of blasting is provided in Equation (33).

$$\begin{array}{l} \text{PPV} = -2.85717 - (0.0211305 \times BF1) + (0.0270673 \times BF2) + (0.0190881 \times BF3) \\ + (0.033926 \times BF5) - (0.0570272 \times BF6) + (5.46015 \times 10^{-5} \times BF7) \\ + (3.56504 \times 10^{-5} \times BF10) + (2.79304 \times 10^{-5} \times BF10) \end{array} \tag{33}$$

#### 4.1.4. ELM

With respect to the experimental results shown in Table 10, the optimum ELM model developed had 12 neurons in the hidden layer with a sigmoid activation function. Thus, having a structure [7-12-1] that represents seven inputs with 12 neurons in the hidden layer and one output.


**Table 10.** Training and Testing *R* and MSE Results for ELM.

#### 4.1.5. MVRA

The developed MVRA model hsa an *R*-value of 0.7909 for the training dataset and 0.8310 for the test dataset. With respect to the MSE, the developed MVRA model had a value of 3.8341 for the training dataset and 3.2456 for the test dataset. Thus, the developed MVRA model using the training datasets for this study is shown in Equation (34).

$$\begin{array}{l} \text{PPV} = 7.237178 + 0.714419AD - 2.80436B + 3.443905S - 0.02705MC \\ \quad - 2.33861SL - 0.67419PF - 0.00284D \end{array} \tag{34}$$

#### *4.2. Assessment of Models Performance*

In evaluating the prediction capabilities of the five predictive models presented in the study, the statistical performance outcomes of the testing samples are outlined in Table 11.


**Table 11.** PPV Prediction Results of Various Models.

Notionally, a predictive model is said to be accurate if *R* and *R*<sup>2</sup> are 1, MSE is 0 and VAF is 100%. In that regard, it can the seen that the GPR with the MSE value of 0.0903 closest to 0, *R* values of 0.9985 closest to 1, *R*<sup>2</sup> values of 0.9971 closest to 1 and VAF value of 99.1728% closest to 100% outperformed all the techniques applied in this study. This shows the reliability of the GPR in predicting ground vibration. The MARS performed better than the ELM by having had MSE value of 0.1038 and a VAF value of 98.5469% with the ELM having an MSE value of 0.1381 and a VAF value of 98.2273%. The ELM also performed better than the BPNN with MSE and VAF values of 0.2178 and 98.1919%. It is worth mentioning that the GPR, MARS, ELM and BPNN were superior in predicting ground vibration to the simple MVRA model which had an MSE of 3.2456, *R*-value of 0.8310, the *R*<sup>2</sup> value of 0.6906 and VAF value of 66.0603%. Figure 7 depicts the interpretation of the obtained results.

**Figure 7.** Comparison of Predicted and Measured PPV for: (**a**) BPNN (**b**) GPR (**c**) ELM (**d**) MARS (**e**) MVRA.

As ground vibration is one of the most unfavorable environmental effects of blasting operations which can cause serious damage to neighboring residences and structures, a precise prediction of its severity is critical to managing and lessening its incidence. The *R*, *R*<sup>2</sup> and VAF values for the GPR, MARS, BPNN, and ELM may not vary significantly, but any predictive model that delivers the most accurate prediction is of paramount relevance to the blast engineer. Hence the need to develop different models. This study found that the GPR is more accurate in forecasting ground vibration than the MARS, BPNN, ELM and MVRA and that it can be used by blast engineers to predict blast-induced ground vibration.

#### *4.3. Sensitivity Analysis*

To determine the most and least effective parameters, sensitivity analysis is performed to examine how the model responds to changes in the input variables with respect to PPV. Hence, in this study, a sensitivity analysis approach implemented in [93] was adopted. Here, while keeping the ranges of all other parameters fixed, the mean value of one of the

input variables is increased (i.e., New mean = Old mean + 5% Old Mean) and subsequently the amount of changes in the predicted PPV using the GPR model is recorded. The obtained results are graphically illustrated in Figure 8.

**Figure 8.** The Input Parameters and PPV Relationship's Strength.

As can be seen in Figure 8, increasing the mean values of spacing and maximum charge per delay, increases PPV. Furthermore, increasing the mean values of distance and stemming length decreases PPV. Increasing burden slightly increased PPV. Nevertheless, increasing values of powder factor and average hole depth did not significantly impact the values of PPV. It can thus be said that the most influential parameters that can affect PPV greatly are spacing, a maximum charge per delay, distance and stemming length.

#### **5. Conclusions**

In this paper, three AI models of GPR, ELM and BPNN were developed and applied to predict blast-induced PPV. In showing the predictive capabilities of these AI techniques, a MARS and MVRA model were also developed. In developing and evaluating these models, 101 datasets obtained from Ras Limestone Mine of Shree Cement, India were utilized. Out of the 101 datasets, 81 were utilized to create the various models, while the remaining 20 were used as test sets for the models that were developed. The input parameters in the creation of the various models were average depth (m), burden (m), spacing (m), powder factor (t/kg), the distance between the monitoring station and the blasting site (m), stemming length (m), and maximum charge per delay (kg), while the output parameter was PPV. The various developed models were then evaluated using performance metrics of *R*, *R*2, MSE and VAF. The results obtained showed that the GPR model had the lowest MSE of 0.0903, and the highest *R*, *R*2, and VAF values of 0.9985, 0.9971 and 99.1728% respectively, indicating that it was superior to the other models in predicting blasting-induced ground vibration. This was followed by MARS which had MSE, *R*, *R*<sup>2</sup> and VAF values of 0.1038, 0.9953, 0.9906 and 98.8692% respectively. Then ELM had an MSE of 0.1381, *R*-value of 0.9957, *R*<sup>2</sup> value of 0.9915 and VAF value of 98.5469. Then the BPNN with an MSE, *R*, *R*<sup>2</sup> and VAF of 0.1714, 0.9924, 0.9848 and 98.2273% respectively. The MVRA performed very poorly as it had, with the highest MSE of 3.2456, and lowest *R*-value of 0.8310, the *R*<sup>2</sup> value of 0.6906 and the VAF value of 66.0603%. The results obtained show that the GPR model can be utilized to forecast blast-induced ground vibration in the mining industry. The sensitivity analysis of the dataset found that spacing, a maximum charge per delay, distance and stemming length had a great influence on PPV whereas burden, powder factor and average depth had slight to no influence on PPV.

**Author Contributions:** Conceptualization, E.T.M., R.M.B., C.K.A., M.K.; methodology, R.M.B., C.K.A.; software, R.M.B., C.K.A.; formal analysis, R.M.B., C.K.A.; resources, E.T.M., R.M.B., C.K.A.; data curation, R.M.B. writing—original draft, M.B., E.T.M., R.M.B., C.K.A., M.K., M.M.S.S., S.K.; writing—review and editing, M.B., E.T.M., R.M.B., C.K.A., M.K., M.M.S.S., S.K.; Supervision, E.T.M., M.K., S.K.; funding acquisition, M.M.S.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research is partially funded by the Ministry of Science and Higher Education of the Russian Federation under the strategic academic leadership program 'Priority 2030 (Agreement 075-15-2021-1333 dated 30 September 2021).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** Authors are thankful to Pankaj Agarwal, Assistant Vice President and Management of Shree Cement, Beawar, Rajashthan for providing data for the preparation of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:



#### **References**

