Effective Assessment of Blast-Induced Ground Vibration Using an Optimized Random Forest Model Based on a Harris Hawks Optimization Algorithm

Abstract

Most mines choose the drilling and blasting method which has the characteristics of being a cheap and efficient method to fragment rock mass, but blast-induced ground vibration damages the surrounding rock mass and structure and is a drawback. To predict, analyze and control the blast-induced ground vibration, the random forest (RF) model, Harris hawks optimization (HHO) algorithm and Monte Carlo simulation approach were utilized. A database consisting of 137 datasets was collected at different locations around the Tonglvshan open-cast mine, China. Seven variables were selected and collected as the input variables, and peak particle velocity was chosen as the output variable. At first, an RF model and a hybrid model, namely a HHO-RF model, were developed, and the prediction results checked by 3 performance indices to show that the proposed HHO-RF model can provide higher prediction performance. Then blast-induced ground vibration was simulated by using the Monte Carlo simulation approach and the developed HHO-RF model. After analyzing, the mean peak particle velocity value was 0.98 cm/s, and the peak particle velocity value did not exceed 1.95 cm/s with a probability of 90%. The research results of this study provided a simple, accurate method and basis for predicting, evaluating blast-induced ground vibration and optimizing the blast design before blast operation.

1. Introduction

In blasting, only 25% to 30% of the explosive energy is spent in rock mass fragmentation, and more than 70% is wasted and causes side effects such as back break, fly rock, and blast-induced ground vibration [1,2,3]. Among these blast-induced side effects, blast-induced ground vibration is considered as the most common, important and dangerous side effects for nearby structures, human life and the environment [4,5,6]. Therefore, establishing the prediction methods to predict and analyze the blast-induced ground vibration and then to optimize the blasting design and building structure by using that model is of great significance to decrease the damage caused by blast-induced ground vibration.

Blast-induced ground vibration can be described by using peak particle velocity (PPV) and frequency. According to previous literature [7], the peak particle velocity is normally considered to be more important than frequency and always selected to describe the blast-induced ground vibration.

After reviewing the previous literature, some empirical models [8,9,10,11] were proposed and most of them predict PPV by using the distance between blast blocks to monitor location and max charge per delay. Although some studies show these empirical models can predict PPV in some mines with high precision, there are also some research works [5,12,13] that show these models are not sufficient for all field sites.

In recent years, artificial intelligence (AI) technology has been widely developed and applied in many areas for solving special problems such as the prediction of blast-induced rock movement [14], flyrock distance [15,16], overbreak [17,18], air overpressure [19,20], rockburst [21,22,23], mechanical properties of material [24,25,26,27,28,29], tunnel convergence [30] and tunnel boring machine penetration rate [31]. In practice, AI technology is also applied to predict blast-induced ground vibration and yield high prediction precision [32,33,34]. For example, Zhang et al. [2] combined the particle swarm optimization (PSO) algorithm and extreme gradient boosting machine (XGBoost) for peak particle velocity prediction. In his paper, the PSO was utilized to search the optimal hyperparameters of XGBoost and a promising result was obtained. Shang et al. [35] proposed a hybrid model that is combined with the firefly algorithm (FFA) and an artificial neural network (ANN). The prediction performance of FFA-ANN was compared with classification and regression tree (CART), support vector machine (SVM), and k-nearest neighbor (KNN) models, and the FFA-ANN model provided the best performance. Azimi [36] utilized the genetic algorithm (GA) to optimize the ANN, and have developed an efficient model for PPV prediction after comparing the performance of the neuro-fuzzy inference system by using 70 blast vibration datasets. In 2009, Khandelwal and Singh [37] developed an ANN model for blast-induced ground vibration prediction, and the prediction results of the ANN model is highly encouraging. Bui et al. [38] tested the feasibility of the FCM-QRNN (fuzzy C-means clustering and quantile regression neural network) model and found that the model has higher accuracy than other soft computing models such as ANN and QRNN. Nguyen et al. [39] combined the k-means clustering algorithm (HKM) and ANN and utilized that model to predict blast-induced ground vibration. In this model, the input datasets should be clustered first by the HKM algorithm, and the result obtained will be used to train the ANN model. In addition to these studies, the blast-induced ground vibration in some famous mines such as the Nui Beo open-cast mine [7], Miduk copper mine [40] and Gol-E-Gohar iron ore mine [41] were studied by using the artificial intelligence technology in previous literature. It can be found that the artificial intelligence technology shows strong prediction ability in blast-induced ground vibration prediction, and research into new AI technology and the new prediction model is necessary and meaningful for hazard control.

In this study, AI technology was tested to predict the blast-induced ground vibration of Tonglvshan open-cast mine. According to the no free lunch theorem, if the performance of algorithm A is better than algorithm B in this problem, there must be other problems that algorithm B has better performance at resolving. Hence, although a series of prediction models such as HKM-ANN [39], MLPNN (multilayer perceptron neural network) [42] and CART [43] were utilized in some mines, the proposed algorithms are far from solving all possible problems and being suited for all field sites. As far as the authors know, the combination of the Harris hawks optimization algorithm (HHO) and random forest (RF) has not been tried in predicting blast-induced ground vibration. In this paper, a new hybrid model namely HHO-RF is proposed to predict and control the PPV value of Tonglvshan open-cast mine. Considering the uncertainty of blasting parameters in field sites, a Monte Carlo process was considered as a stochastic approach to analyze the probability distribution with various combinations of blasting parameters for hazard control. Eventually, a sensitivity investigation was carried out to measure the impact of various parameters on peak particle velocity.

After analyzing, engineers can use the accurate prediction model proposed in this study to predict and evaluate the PPV value caused by blast operation before a blast, which is helpful for optimizing blast design and reducing disputes between mines and surrounding residents. In addition, the security department can stipulate measures for protecting the buildings and industrial facilities based on the probability distribution of PPV because all the probable combinations of variables will be examined. Last, the analysis method and research results in this paper can provide a reference for controlling blast-induced ground vibration in similar mines.

2. Materials and Methods

2.1. Data Set

To achieve the objectives of this study, blast-induced ground vibration of Tonglvshan open-cast mine was monitored by using an accurate blast-induced ground vibration monitoring machine (EXP3850-3) in Hubei province, China (see Figure 1). That mine is situated 3 km west of Daye city, and the drilling and blasting technique is utilized for rock fragmentation. In addition, the geological structure of the Tonglvshan open-cast mine is a skarn deposit, and the ore deposit is mainly composed of granodiorite porphyry, sodium porphyry, skarn, marble, copper-bearing skarn, copper-iron ore, etc. Meanwhile, the value of Protodyakonov impact strength coefficient of different rock type varies greatly, which leads to the significant difficulty of accurate peak particle velocity prediction.

A total number of 137 blasting operations were analyzed, and 8 variables including peak particle velocity (PPV), Protodyakonov impact strength coefficient (f), max charge per delay (Q_max), total charge (Q_total), delay time of detonator (T), burden of the first row (B), horizontal (D_H) and vertical (D_V) distance between blast block and monitoring station were collected. After analyzing by using the box plot, 20 outliers were moved from the database to improve the model performance [44] and the remaining 117 datasets were utilized in this study (see

Table 1. Some of the datasets used for peak particle velocity (PPV) prediction.

Number	Qmax	Qtotal	DH	DV	B	T	f	PPV
1	160	1440	125.30	52.30	5	50	5	0.343
2	312	3120	311.90	42.00	7	25	8	0.753
3	312	3120	389.40	108.00	7	25	8	0.572
4	312	3120	362.30	86.00	7	25	8	1.214
5	312	3120	199.20	30.00	7	25	8	2.148
6	350	1050	114.70	16.80	5	100	5	0.417
7	350	1050	143.50	52.30	5	100	6	0.609
8	350	1050	133.30	62.00	5	100	5	0.194
9	370	2150	70.30	42.00	7	25	8	4.754
10	370	2150	52.40	30.00	7	25	8	4.585
11	370	2150	101.20	54.00	7	25	6	1.554
12	380	1550	162.40	73.00	4	75	5	0.102
13	380	1550	147.80	62.00	4	75	5	0.201
14	380	1550	104.90	28.60	4	75	5	0.463
15	380	1550	336.20	58.90	5	50	6	0.101
16	390	2730	120.90	46.90	4	100	5	0.143
17	390	2730	69.90	53.90	6	75	8	0.394
18	390	2730	50.40	30.00	6	75	8	0.126
19	390	2730	137.80	63.10	5	50	7	0.318
20	390	2730	63.30	27.60	5	50	7	0.657
21	390	2730	128.90	73.00	5	50	7	0.435
22	400	3600	110.60	16.80	6	50	6	0.898
23	400	3600	47.10	16.90	6	50	6	5.371
24	400	3600	137.90	52.30	6	50	6	0.815
25	400	3600	31.50	16.80	6	50	6	5.436
26	456	1860	284.10	85.30	6	50	6	0.368
27	456	1860	300.50	102.40	6	50	6	0.270
28	456	1860	330.40	58.90	6	50	6	0.144
29	460	4600	443.60	108.00	7	25	8	0.510
30	460	4600	363.30	42.00	7	25	8	0.476
31	460	2760	323.60	30.00	6	50	6	0.203
32	460	2760	313.20	86.00	6	50	6	0.391
33	468	936	85.30	54.00	6	50	6	2.486
34	468	936	75.50	30.00	6	50	6	3.350
35	494	3952	122.60	62.10	5	100	5	0.608
36	494	3952	66.70	28.60	5	100	6	3.029
37	494	3952	183.20	28.90	5	100	5	0.105
38	494	1482	170.40	46.90	5	50	6	0.195
39	494	1482	132.50	63.10	5	50	6	0.413
40	494	5434	149.20	73.30	6	25	7	0.923
41	494	5434	96.40	27.60	6	25	7	1.413
42	494	1482	142.90	73.30	5	50	6	0.231
43	494	1482	81.50	51.90	5	50	6	0.806
44	494	1482	80.20	27.60	5	50	6	0.598
45	507	7650	193.50	53.00	5	100	6	0.106
46	507	7650	235.20	86.00	5	100	6	0.196
47	507	7650	332.30	86.00	5	100	8	0.292
48	532	6084	73.00	40.00	6	25	7	2.372
49	532	6084	124.20	73.00	6	25	7	1.442
50	550	4400	123.40	53.90	5	75	8	1.562
51	550	4400	89.50	42.00	5	75	8	1.279

).

2.2. Random Forest (RF) Model

The RF model was invented by Breiman [45]. In this model, randomly selected input variables and datasets were utilized to train the classification or regression model, and a large number of decision trees were developed and combined to make a final decision. Limited by the article length, the detailed introduction was not shown here but readers can find it in the published papers of Breiman [45] and Liaw and Wiener [46].

2.3. Harris Hawks Optimization (HHO) Algorithm

As an important branch of artificial intelligence, the random forest model has been widely used in many areas such as civil areas [47], energy areas [48] and medical areas [49]. Meanwhile, previous literature [50] shows that the using of meta-heuristic optimization technique can significantly improve prediction performance by search the optimal combination of prediction models. In fact, many meta-heuristic optimization techniques are inspired by natural phenomena and have widely applied in mining [51], healthcare [52,53,54,55], finance [56] and query optimization problems [57] with good optimization performance. Hence, a nature-inspired algorithm (the Harris hawks optimization (HHO) algorithm) was tested in this paper because of its superior performance compared with the other 11 methods through 29 unconstrained benchmark problems and 6 constrained engineering design tasks carried out by Heidari et al. [58].

The Harris hawk (also called Parabuteo unicinctus), a well-known bird in the southern half of Arizona, USA, often be considered as a truly cooperative predator because these birds usually discover, attack and share prey with other family members [58,59] (Figure 2). Inspired by the social behavior of Harris hawks, the HHO algorithm was proposed for solving optimization problems with the proper formulation. The HHO algorithm includes three phases: the exploration phase, the transition from exploration to exploitation and exploitation phase. After reviewing the study of Heidari et al. [58], a brief introduction of the HHO algorithm was performed here.

(1) Exploration phase:

In this part, the Harris hawks try to discover, search and find the prey in the desert site by using eyes. The position of a Harris hawk in the mathematical model means the closest position to the prey and can be obtained by using the following equations [58]:

$$X\left(t+1\right)=\{\begin{array}{ll}{X}_{rand}\left(t\right)-r_1\left|X_{rand}\left(t\right)-2r_{2}X\left(t\right)\right|& q\ge 0.5\\ \left(X_{rabbit}\left(t\right)-X_m\left(t\right)\right)-r_3\left(LB+r_4\left(UB-LB\right)\right)& q<0.5\end{array}$$

(1)

where X(t + 1) and X(t) are the positions of Harris hawk in iteration t + 1 and t, respectively. X_rabbit(t) and X_rand(t) are the position of prey in the t iteration and randomly selected Harris hawk in the current population. r₁, r₂, r₃, r₄, q are the random values, UB and LB are the upper and lower bounds, respectively.

Besides these parameters, X_m(t) is the mean position of Harris hawks.

$$X_{\mathrm{m}}\left(t\right)=\frac{1}{N}{\displaystyle \sum _{i=1}^N{X}_i}\left(t\right)$$

(2)

where N is the number of Harris hawks; X_i(t) is the position of each hawk in iteration t.

(2) The transition from exploration to exploitation:

According to the escaping energy of prey, the behavior of Harris hawks in the HHO algorithm will be changed. For example, the exploitation behavior will be performed when the escaping energy is less than 1, and the Harris hawks continue to search for prey when escaping energy is greater than or equal to 1. In that model, the escaping energy (E) can be obtained as [58]:

$$E=E_0\left(1-\frac{t}{T}\right)$$

(3)

where E₀ is the initial escaping energy, and E is the escaping energy in iteration t, and T is the maximum number of iterations.

(3) Exploitation phase

After searching for prey in the desert site, a Harris hawk will attempt to attack the prey and the prey will correspondingly attempt to escape the surprise pounce. According to the different escaping behavior of the prey, several possible strategies may be performed during the exploitation phase. To evaluate the possibility of successfully escaping, r, a random value range from 0 to 1, was utilized in this algorithm, and four situations were divided with different range of r value and E value.

(a)

$$r\ge 0.5\&\left|E\right|\ge 0.5$$

In this situation, the prey still has enough energy, so the Harris hawks hope to exhaust it and then perform the attack. The following equations can be used to calculate the position of the Harris hawk in iteration t + 1 [58].

$$X\left(t+1\right)=\left(X_{rabbit}\left(t\right)-X\left(t\right)\right)-E\left|2X_{rabbit}\left(t\right)\left(1-r_5\right)-X\left(t\right)\right|$$

(4)

where r₅ is a random value range from 0 to 1.

(b)

$$r\ge 0.5\&\left|E\right|<0.5$$

In this situation, the prey is tired and the Harris hawks can perform the surprise pounce. The position can be updated using the following equation [58]:

$$X\left(t+1\right)=X_{\mathrm{ra}bbit}\left(t\right)-E\left|\mathsf{\Delta }X\left(t\right)\right|$$

(5)

(c)

$$r<0.5\&\left|E\right|\ge 0.5$$

In this situation, the prey still has enough energy and has a low chance of successfully escaping. The position of Harris hawks can be updated by the following rules [58]:

$$Y=\mathsf{\Delta }X\left(t\right)-E\left|J{X}_{rabbit}\left(t\right)-X\left(t\right)\right|$$

(6)

The hawks can dive according to the following rule:

$$Z=Y+S\times 0.01\times \frac{u\times \sigma }{{\left|v\right|}^{\frac{1}{\beta }}}$$

(7)

Subjected to

$$\sigma ={\left(\frac{\mathsf{\Gamma }\left(1+\beta \right)\times \sin \left(\frac{\pi \beta }{2}\right)}{\mathsf{\Gamma }\left(\frac{1+\beta }{2}\right)\times \beta \times 2^{\left(\frac{\beta -1}{2}\right)}}\right)}^{\frac{1}{\beta }}$$

here, $\mathsf{\Gamma }$ is Gamma function; β is 1.5; S is a random vector; u and v are the random value range from 0 to 1.

After calculation, the final position can be shown as follows:

$$X\left(t+1\right)=\{\begin{array}{l}\begin{array}{cc}Y& ifF\left(Y\right)<F\left(X\left(t\right)\right)\end{array}\\ \begin{array}{cc}Z& ifF\left(Z\right)<F\left(X\left(t\right)\right)\end{array}\end{array}$$

(8)

(d)

$$r<0.5\&\left|E\right|<0.5$$

In this situation, the prey is exhausted and the escaping chance is low. Therefore, the Harris hawks can perform the hard siege behavior and can be modeled as follows [58]:

$$X\left(t+1\right)=\{\begin{array}{l}\begin{array}{cc}Y& ifF\left(Y\right)<F\left(X\left(t\right)\right)\end{array}\\ \begin{array}{cc}Z& ifF\left(Z\right)<F\left(X\left(t\right)\right)\end{array}\end{array}$$

(9)

where

$$\begin{gathered} Y=X_{rabbit}\left(t\right)-E\left|J{X}_{rabbit}\left(t\right)-X_m\left(t\right)\right| \\ Z=Y+S\times LF\left(D\right) \end{gathered}$$

here, LF is levy flight function.

2.4. HHO-RF Model

The procedure of developing an HHO-RF model is shown in Figure 3. Firstly, the collected database should be divided into training datasets and testing datasets with a ratio of 80% and 20% [60,61,62,63]. Then the Harris hawks optimization algorithm will be used to search the best combination of the number of variables used to grow each try (m_try) and the number of trees (n_tree) in the RF model, and the final RF model will be established. Finally, the testing datasets can be used to evaluate the RF model by using some performance metrics such as correlation coefficient (R²), mean absolute error (MAE) and root mean square error (RMSE).

2.5. Performance Metrics

In this section, an RF model and a HHO-RF model were developed to establish a stable and accurate relationship between the input variables (f, B, T, D_V, D_H, Q_total, Q_max) and output variable (PPV). To evaluate the predictive performance, three performance indices including RMSE, R² and MAE, were introduced and utilized here [64,65,66,67,68].

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(f\left(x_i\right)-x_i\right)}^2}}$$

(10)

$$R^2=\frac{{\left(n{\displaystyle \sum _{i=1}^{n}f\left(x_i\right)x_i}-{\displaystyle \sum _{i=1}^{n}f\left(x_i\right)}{\displaystyle \sum _{i=1}^n{x}_i}\right)}^2}{\left(n{\displaystyle \sum _{i=1}^{n}f{\left(x_i\right)}^2}-{\left({\displaystyle \sum _{i=1}^{n}f\left(x_i\right)}\right)}^2\right)\left(n{\displaystyle \sum _{i=1}^n{x}_i^2}-{\left({\displaystyle \sum _{i=1}^n{x}_i}\right)}^2\right)}$$

(11)

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|x_i-f\left(x_i\right)\right|}$$

(12)

here x_i, f(x_i) and n are the actual PPV value, predicted PPV value and number of input datasets, respectively.

2.6. Monte Carlo Simulation

The Monte Carlo (MC) simulation method is a typical sampling technique in stochastic modeling and has the ability to realize the risk and uncertainty in many problems such as financial cost and purchase decision [69]. Monte Carlo simulation can achieve two important functions; the first is the quantitative investigation of uncertainty and variation of problems regarding risk; and the second is the investigation of parameters affecting variation and uncertainty and their percentage [70]. In the calculation process, the repeated sampling method was utilized to create a large number of combinations of input variables, then these input variables will be calculated by using the HHO-RF model developed. Finally, statistical analysis will be carried out for obtaining the probability distribution of PPV values (see Figure 4).

In recent years, the MC simulation approach has been utilized to solve many engineering problems such as flyrock assessment [70], backbreak analysis [71], compression coefficient prediction [72], blasting fragmentation [73] and drilling rate evaluation [74]. In the following section, the MC simulation method will be utilized to analyze the probability distribution of blast-induced ground vibration caused by the blasting operation in Tonglvshan open-cast mine.

3. Results and Discussion

3.1. Results of RF Model

In order to predict the PPV by using the RF model, two typical model parameters including m_try and n_tree should be determined first. Hence a parameter investigation was carried out by using the grid search method with m_try of 1 to 7 and n_tree of 2000, and 7 was selected as the optimal m_try value (Figure 5a). In addition, the parameter investigation with m_try of 7 and n_tree of 1 to 2000 show that the prediction performance will vary within a range and that range will be stable when n_tree = 2000. So 2000 was selected as the optimal value of n_tree (Figure 5b). After optimal RF model development, the testing datasets were inputted into the established RF model for prediction, the prediction performance shows the prediction performance with using the training datasets with R² of 0.92, MAE of 0.19 and RMSE of 0.27 and using the testing datasets with R² of 0.90, MAE of 0.35 and RMSE of 0.40 (see Figure 6).

3.2. Results of HHO-RF Model

To obtain the best combination of m_try and n_tree, 8 HHO-RF models with the maximum number of iterations of 500 and variable numbers of Harris hawks ranging from 20 to 250 were developed. The fitness curve of these HHO-RF models checked by RMSE can be found in Figure 7. As shown in Figure 7, there is no significant change after 350 iterations. But the escaping energy of the rabbit has a relationship with the maximum number of iterations in the HHO algorithm according to Equation (3). The change of the maximum number of iterations may affect the optimization performance of the HHO algorithm. Hence, 500 and 150 were chosen as the best parameter combination of HHO-RF modeling in view of the excellent prediction performance.

As shown in Figure 8, the established HHO-RF model can predict and evaluate the PPV with R² of 0.94, MAE of 0.17 and RMSE of 0.23 for training datasets and R² of 0.94, MAE of 0.29 and RMSE of 0.34 for testing datasets. After comparing the prediction performance obtained by the RF model and HHO-RF model, the comparison results show that the HHO-RF model can be used to predict the PPV with a higher level of prediction accuracy.

3.3. Results of Monte Carlo Simulation

As mentioned in the previous sections, an RF model and a HHO-RF model were developed for peak particle velocity prediction, and the results obtained show that the HHO-RF model can provide more accurate prediction performance. Hence, the HHO-RF model was used here for Monte Carlo simulation, and four input variables including D_V, D_H, Q_total, and Q_max were considered as the continuous probability distributions (CPD) and two input variables including T and B were considered as the discrete probability distribution (DPD) (see

Table 2. Probability distribution functions of seven input variables.

Input Variables	Function
T	Discrete ({0, 25, 50, 75, 100}, {0.06, 0.17, 0.34, 0.24, 0.19})
B	Discrete ({4, 5, 6, 7}, {0.17, 0.41, 0.36, 0.06})
f	Discrete ({5, 6, 7, 8}, {0.15, 0.42, 0.18, 0.24})
DV	Weibull (A, B): A = 65.02, B = 2.62
DH	Logarithmic (mean, std. dev): mean = 5.06 std. dev = 0.54
Qtotal	Weibull (A, B): A = 4613.05 B = 2.15
Qmax	Logarithmic (mean, std. dev): mean = 6.64 std. dev = 0.58

).

After inputting the probability distribution functions into the MC simulation, a large number of combinations should be randomly selected by using the sampling method. In this regard, two sampling methods including Latin hypercube sampling and simple random sampling are recommended. According to the previous literature, the Latin hypercube sampling method is much better than the simple random sampling method because it can achieve the same calculation performance with fewer simulations. It should be noted that the relationship between input variables should be considered in the MC simulation to guarantee the consistency of input variable combinations in MC simulation and monitoring databases.

In order to achieve a high level of simulation and probability distribution, 20,000 simulation runs were performed by using the Latin hypercube sampling method, and the relationship of 20,000 randomly selected combinations can be found in Figure 9. It is found that the correlations between input values in MC simulation and blast-induced ground vibration database agree well, so the 20,000 simulation runs developed can be used to analyze the probability distribution of blast-induced ground vibration.

After analyzing, the probability distribution of blast-induced ground vibration is shown in Figure 10. As can be seen, the mean value of PPV in MC simulation results is determined as 0.98 cm/s while the mean value of PPV is 1.0 cm/s in the database, and the PPV value does not exceed 1.95 cm/s with a probability of 90%. Moreover, the measured PPV values, the predicted PPV values by the HHO-RF model, and the Monte Carlo simulation results are also shown in Figure 10. It can be observed that these three cumulative frequency curves are closely related and the curve performed by the MC simulation is smoother. It is also found that these cumulative frequency curves still have some differences, and the differences may be due to the small amount of the monitoring database being too small and the error of the prediction model.

After obtaining the probability distribution of PPV values, the engineer can adjust the blasting parameters such as the max charge per delay and burden to control the peak particle velocity and minimize the damage upon the building structure. Moreover, for a case that needs special protection, the engineer is advised to assume the most dangerous situation with the probability of 0.001% to arrange safeguards.

3.4. Variable Importance

To analyze the impact of various blasting parameters on blast-induced ground vibration, a sensitivity analysis was performed to measure the impact of various parameters on peak particle velocity. In this section, with the help of SPSS (Statistical Product and Service Solutions), the relationships between two variables can be calculated with the range of −1 to +1 after sorting the database and comparing the rank of each value. Figure 11 shows the impact of these input variables on PPV values both in blast-induced ground vibration database and MC simulation results. As shown, the D_H with the correlation coefficient of −0.795 and −0.563 is the most sensitive variable, and the input variables can be divided into controllable parameters and uncontrollable parameters. The controllable parameters can be listed in descending order as D_H, D_V, Q_total, Q_max, T and B in the PPV monitoring database and D_H, D_V, Q_total, Q_max, B, T in MC simulation results. It is found that there is a difference in the ranking of T and B in these two databases. The reason may be that the PPV monitoring database is too small and the ranking order obtained from the MC simulation results is more reliable according to the big data principle. Meanwhile, the roughly same order of input variables demonstrates that the Monte Carlo simulation is reliable.

3.5. Limitations

Four limitations can be found in this paper and need to be considered in the future study. First, only 137 datasets were collected and utilized to develop the RF model and HHO-RF model, a larger database should be established for evaluating, improving model reliability and obtaining more accurate prediction performance. Then, only 8 variables including PPV, f, Q_max, Q_total, T, B, D_H, and D_V were collected, more qualitative parameters such as the detonation velocity of the explosive may have an impact on the PPV prediction and should be considered in future investigations. Third, although the HHO algorithm shows good results in improving the prediction performance of the RF model, there are many optimization algorithms such as the GA, grey wolf optimization algorithm (GWO), artificial bee colony algorithm (ABC) and cuckoo search algorithm (CS) that can also be tried to select the hyperparameters of the random forest algorithm and no prediction performance comparison was carried out in this paper for these optimization algorithms. Finally, as a typical character of the machine learning method, the prediction performance is dependent on the fitness of the collected database and the complex relationship between prediction results and input variables is difficult to describe due to its “black box” nature.

4. Conclusions

This study aims to propose a new hybrid model for accurate blast-induced ground vibration prediction and use the stochastic modeling approach to obtain the probability distribution. At first, a simple RF model and a new proposed HHO-RF model were developed and utilized to predict the peak particle velocity; the prediction performance showed by HHO-RF model is much better than the prediction performance yield by RF model. In addition, the Monte Carlo simulation showed that the mean PPV value of the MC simulation results is 0.98 cm/s and the PPV value does not exceed 1.95 cm/s with the probability of 90%. Finally, the results of the sensitivity analysis showed that the D_H is the most influential factor on blast-induced ground vibration. Based on the high-precision prediction model developed and the results of Monte Carlo simulation obtained, the engineers of Tonglvshan open-cast mine can accurately predict blast-induced ground vibration and ameliorate blast design before a blast operation, which is beneficial for controlling the blast’s impact on nearby structures, human life, and the environment.