Prediction and Optimization of Blasting-Induced Ground Vibration in Open-Pit Mines Using Intelligent Algorithms

: Prediction and parameter optimization are effective methods for mine personnel to control blast-induced ground vibration. However, the challenge of effective prediction and optimization lies in the multi-factor and multi-effect nature of open-pit blasting. This study proposes a hybrid intelligent model to predict ground vibrations using a least-squares support vector machine (LSSVM) optimized by a particle swarm algorithm (PSO). Meanwhile, multi-objective particle swarm optimization (MOPSO) was used to optimize the blast design parameters by considering the vibration of particular areas and the bulk rate of blast fragmentation. To compare the prediction performance of PSO-LSSVM, a genetic-algorithm-optimized BP neural network (GA-BP), unoptimized LSSVM, and BP were used, by applying the same database. In addition, the root-mean-squared error (RMSE), the mean absolute error (MAE), and the correlation coefﬁcient ( r ) were regarded as the evaluation indicators. Furthermore, the optimization results of the blasting parameters were obtained by quoting the established vibration prediction model and bulk rate proxy model in MOPSO and veriﬁed by ﬁeld tests. The results indicated that the PSO-LSSVM model provided the highest efﬁciency in predicting vibrations with an RMSE of 1.954, MAE of 1.717, and r of 0.965. Furthermore, the blasting vibration can be controlled by using the two-objective optimization model to obtain the best blasting parameters. Consequently, this study can provide more speciﬁc recommendations for vibration hazard control.


Introduction
Blasting is the preferred technique for open-pit mining activities, primarily aimed at achieving rock fragmentation.However, the release of blast energy will have adverse effects, e.g., seismic waves, flying rocks, and dust [1].Hence, pre-blast prediction and optimization are essential tools to control blasting effects.In addition, research on blasting effects is shifting from traditional empirical formulas to intelligent algorithmic predictions.
To predict and optimize blasting effectiveness, researchers have conducted many studies.The traditional prediction methods are based on empirical models.Zhu et al. [2] modified the Kuz-Ram empirical model to improve the accuracy of blasting block prediction.Matidza et al. [3] compared the performance of the Sadovsky formulation with other empirical models in predicting blast vibration.Recently, intelligent algorithms have gained widespread use in predicting blast-induced outcomes, including the predicting of adverse effects, e.g., backbreak [4], dust emissions [5], and vibration [6], the predicting of direct blasting results, e.g., block degree [7] and throwing distance [8].Moreover, blasting impacts the operations of various production segments in open-pit mines [9].For example, bulk ore produced by blasting implies additional loading, transportation, and crushing costs [10].In addition, blast cracking and bulkiness affect crushing and grinding [11].Therefore, various indirect effects can be predicted based on blasting, such as the prediction of shovel loading time, truck payload [12], production energy consumption [13], etc. Above all, the objects of prediction and optimization for blasting effectiveness are multifaceted [14].As the research of intelligent algorithms in mining and natural resources is intensifying, the prediction and optimization of multiple objects by intelligent algorithms are receiving more attention [15][16][17].
The different contributions of explosive energy are essentially responsible for the diversity of blasting results, i.e., a smaller fraction of energy (about 20-30%) is used for rock fragmentation and removing, whereas most of the energy is wasted in various adverse effects [18].Among the various adverse effects, blasting-induced ground vibrations are considered the most harmful [19], accounting for approximately 40% of the explosive energy [20].For example, blasting vibrations affect crushers and explosives magazines in open-pit mines [21].Therefore, when there are important buildings and unstable terrain around the blasting area, blasting-induced ground vibration should be strictly controlled [22].Peak particle velocity (PPV) is generally considered a crucial indicator of blast-induced ground vibration intensity [23].Over time, numerous empirical models have been developed and proven effective in predicting PPV [3,24].However, a literature review found that empirical models are less accurate [25,26].Therefore, adopting intelligent algorithms and incorporating more factors affecting PPV into prediction models improves prediction accuracy [27][28][29].
Scholars have paid significant attention in recent years to enhancing intelligent algorithms for accurate PPV prediction.For example, Faradonbeh and Monjezi [20] used the cuckoo optimization algorithm (COA) to optimize gene expression programming (GEP), and Zhang et al. [30] optimized XGBoost using the particle swarm algorithm (PSO).Therefore, the hybrid models predict PPV more effectively, e.g., by optimizing the parameter configuration of prediction models.In addition, preprocessing the input data using algorithms (e.g., clustering techniques and feature selection) can enable predictive models to exhibit a dominant performance in a comprehensive comparison [19,31,32].Comparing the hybrid model with a single model, some scholars found a significant improvement in PPV prediction accuracy with the hybrid approach [33].However, these situations are common in that the prediction results based on blasting parameters cannot directly give reference values, especially since the adjustment of blasting design parameters still needs to rely on the engineers' experience.
To improve the prediction accuracy, metaheuristic algorithms are utilized in hybrid models for predicting blast vibration [20,30,33].Additionally, metaheuristic algorithms are also widely applied in mine blasting optimization, e.g., the grasshopper optimization algorithm (GOA) is employed to identify the blasting solution that minimizes dust generation [34].Meanwhile, multi-objective particle swarm optimization (MOPSO), developed from PSO, is capable of handling multiple metrics in mining applications [35].
Therefore, establishing a highly reliable prediction model and giving reference values of blasting design parameters is particularly important for engineers in controlling blasting vibration hazards and guaranteeing blasting effectiveness.This study proposes a hybrid model of PSO and least-squares support vector machine (PSO-LSSVM) to predict the PPV of a sand and gravel mine.Furthermore, the LSSVM, GA-BP, and BP were considered for comparison purposes.In addition, the optimization strategy of the blasting design parameters is proposed based on multi-objective particle swarm optimization, and a blast vibration control scheme for sensitive areas is given.

Project Overview
This study was conducted at a gravel mine on a sea island in Zhoushan, Zhejiang Province.The mine has an original high-steep slope (Area e), as shown in Figure 1.Being mined from the top of the slope means the steps are negative relative to the elevation of the detonation point, which proves that the propagation of blasting seismic waves has a negative elevation effect under such geological conditions [36].The mine mines the ore to produce sand and gravel aggregate, with a production scale of 4.74 million m 3 /year.Notably, due to the topography, the ore-crushing plant (Area c) is adjacent to the slope and is a typical blasting-vibration-sensitive area.
In this mine, the ore is mainly crystalline bosom tuff, rhyolite porphyry, and sunken tuff, with a dense structure and hard rock type.Besides, the fault and joints are not developed in the mine area, and the rock body is generally intact.Therefore, the area with a single lithology (rhyolite porphyry) was used as the test area (Area b).The main physical and mechanical parameters of the rock body are shown in Table 1.Rock blasting was carried out by deep-hole step blasting with a continuous charge structure and a step height of 15 m.The blasting parameters for the mine rock are shown in Figure 2. The arrow (d) in Figure 2 shows the direction of blast throwing, i.e., back from the existing slope.

Project Overview
This study was conducted at a gravel mine on a sea island in Zhoushan, Zhejiang Province.The mine has an original high-steep slope (Area e), as shown in Figure 1.Being mined from the top of the slope means the steps are negative relative to the elevation of the detonation point, which proves that the propagation of blasting seismic waves has a negative elevation effect under such geological conditions [36].The mine mines the ore to produce sand and gravel aggregate, with a production scale of 4.74 million m 3 /year.Notably, due to the topography, the ore-crushing plant (Area c) is adjacent to the slope and is a typical blasting-vibration-sensitive area.
In this mine, the ore is mainly crystalline bosom tuff, rhyolite porphyry, and sunken tuff, with a dense structure and hard rock type.Besides, the fault and joints are not developed in the mine area, and the rock body is generally intact.Therefore, the area with a single lithology (rhyolite porphyry) was used as the test area (Area b).The main physical and mechanical parameters of the rock body are shown in Table 1.Rock blasting was carried out by deep-hole step blasting with a continuous charge structure and a step height of 15 m.The blasting parameters for the mine rock are shown in Figure 2. The arrow (d) in Figure 2 shows the direction of blast throwing, i.e., back from the existing slope.

Project Overview
This study was conducted at a gravel mine on a sea island in Zhoushan, Zhejiang Province.The mine has an original high-steep slope (Area e), as shown in Figure 1.Being mined from the top of the slope means the steps are negative relative to the elevation of the detonation point, which proves that the propagation of blasting seismic waves has a negative elevation effect under such geological conditions [36].The mine mines the ore to produce sand and gravel aggregate, with a production scale of 4.74 million m 3 /year.Notably, due to the topography, the ore-crushing plant (Area c) is adjacent to the slope and is a typical blasting-vibration-sensitive area.
In this mine, the ore is mainly crystalline bosom tuff, rhyolite porphyry, and sunken tuff, with a dense structure and hard rock type.Besides, the fault and joints are not developed in the mine area, and the rock body is generally intact.Therefore, the area with a single lithology (rhyolite porphyry) was used as the test area (Area b).The main physical and mechanical parameters of the rock body are shown in Table 1.Rock blasting was carried out by deep-hole step blasting with a continuous charge structure and a step height of 15 m.The blasting parameters for the mine rock are shown in Figure 2. The arrow (d) in Figure 2 shows the direction of blast throwing, i.e., back from the existing slope.

Data Collection
In order to establish the prediction model of PPV and optimize the blasting parameters, a statistical scheme of the PPV data, including the blast design parameters and negative elevation (H), was developed.The values of the maximum explosive charge capacity (W), monitoring distance (R), negative elevation (H), burden (B), spacing (S), subdrilling (H 0 ), powder factor (PF), and PPV were collected and recorded.In this statistical scheme, W, R, and B have been used by other scholars as controllable parameters to predict PPV [37], while S and PF were also shown by Nguyen and Bui [38] to be related to blasting vibration.In addition, negative elevation steps can have a decreasing effect on vibration velocity [36].Therefore, H and H 0 (meaning the maximum depth at which explosives are placed) were counted and used as input variables.The monitoring equipment for PPV caused by blasting vibration was the TC-4850 blasting [39] vibration tester, and the monitoring points were arranged as shown in Figure 3. R and H were determined by considering laser ranging and calibration on the construction drawings.A total of 30 blasting events were investigated.Due to topographic constraints, PPV in high-steep slope and blasting steps were monitored separately for different periods of blasting.Five measurement points were arranged in front of each blasting event.In addition, less-effective recovery data were obtained about the high-steep slope because of the large number of falling rocks.Therefore, a total of 50 sets of PPV data were recorded.The characteristics of the recorded input and output data are shown in Table 2.

Data Collection
In order to establish the prediction model of PPV and optimize the blasting parameters, a statistical scheme of the PPV data, including the blast design parameters and negative elevation (H), was developed.The values of the maximum explosive charge capacity (W), monitoring distance (R), negative elevation (H), burden (B), spacing (S), subdrilling (H0), powder factor (PF), and PPV were collected and recorded.In this statistical scheme, W, R, and B have been used by other scholars as controllable parameters to predict PPV [37], while S and PF were also shown by Nguyen and Bui [38] to be related to blasting vibration.In addition, negative elevation steps can have a decreasing effect on vibration velocity [36].Therefore, H and H0 (meaning the maximum depth at which explosives are placed) were counted and used as input variables.The monitoring equipment for PPV caused by blasting vibration was the TC-4850 blasting [39] vibration tester, and the monitoring points were arranged as shown in Figure 3. R and H were determined by considering laser ranging and calibration on the construction drawings.A total of 30 blasting events were investigated.Due to topographic constraints, PPV in high-steep slope and blasting steps were monitored separately for different periods of blasting.Five measurement points were arranged in front of each blasting event.In addition, less-effective recovery data were obtained about the high-steep slope because of the large number of falling rocks.Therefore, a total of 50 sets of PPV data were recorded.The characteristics of the recorded input and output data are shown in Table 2.

Particle Swarm Optimization Algorithm
PSO belongs to the intelligent swarm algorithm (SI), which is derived from the observation and simulation of insect populations in nature.The basic principle is that the particle dynamically adjusts the movement speed and position according to its optimal position and the globally known optimal position, and the selection of the optimal position is determined by the fitness function [40].In dealing with the optimization of indicators with multiple parameters, PSO has a powerful function.The dynamic adjustment of the particles is achieved by two core formulas, i.e., the velocity update formula and the position update formula, as follows: Among them, ε is the inertia weight coefficient, i.e., the weight of the updated particle velocity vector; c 1 is the individual learning acceleration coefficient; c 2 is the global learning acceleration coefficient; rand(1) and rand(2) are random functions, ranging from 0 to 1; LP ik and GP ik are the local optimal position and the global optimal position of the kth iteration of the ith particle, respectively.

Least-Squares Support Vector Machine
SVM is an important branch of machine learning algorithms, minimizing structural risk as a fundamental principle, i.e., achieving a balance between empirical error and confidence intervals to improve the generalization ability of the algorithm.LSSVM was improved by SVM [41], which converts the quadratic inequality constraint of the optimization function into the equation constraints of the linear equation system, and the training time of the algorithm is shortened.At the same time, LSSVM inherits the advantages of SVM in dealing with problems such as small samples, nonlinearity, and local minima.The main principles of LSSVM are as follows: Given a training sample dataset D = {(x 1 , y 1 ), (x 2 , y 2 ), . . ., (x N , y N )}, i = 1 to N, where x i ∈ R d , i.e., a total of N d dimensions for input vectors x i , y i , is the output indicator.A nonlinear function maps the low-dimensional samples to the high-dimensional feature space.The nonlinear function estimation modeling is given as where ω is the weight vector and b is the bias term.
In order to evaluate and classify the sample set, the following equality constraint was set as an optimization function in LSSVM.
Among them, γ is the regularization parameter, which is used to determine the tradeoff between model complexity and accuracy; e i is the error vector, which represents the difference between the actual value of the output indicator and the predicted value.
Solve the above optimization function by constructing the Lagrangian function and finding partial derivatives: where ξ i is the Lagrange multiplier.Take the partial derivative for the parameters ω, b, e i , and ξ i , and set the partial derivative to zero [42]; we have: By eliminating ω and e i , the four linear equations of Equation ( 6) are simplified as: and Ω is an N × N symmetry matrix of kernel function, as follows: Equations ( 7) and ( 8) lead to the decision function of the LSSVM model: where the radial basis function (RBF) is widely used as the kernel function.Furthermore, the RBF kernel has also been widely used and shown to achieve the best prediction performance in LSSVM prediction model studies of open-pit blasting vibration [43,44].
where σ is the key parameter of the kernel function.

Optimization of LSSVM by Particle Swarm Algorithm
It can be seen from Equations ( 7) and ( 8) that γ and σ are the critical parameters in determining the weight of the error vector in the LSSVM model [45].Thus, improving the ability to learn and generalize is vital when applying the LSSVM model to predict nonlinear data with a small sample.In order to improve the prediction accuracy, this paper introduced the PSO algorithm into the LSSVM.Figure 4 shows the optimization flow of the hyperparameters (γ and σ) in the LSSVM model using PSO.The root-mean-squared error (RMSE) was chosen as the fitness function for the optimization of the two hyperparameters γ and σ: Among them, n is the number of samples; y i is the actual value of the sample; ŷi is the predicted value of the sample.
Among them,  is the number of samples;  is the actual value of the sample;  is the predicted value of the sample.In the model's training, the particle velocity was updated by Equation (1), and the particle position was updated by Equation (2).The individual and global fitness were updated by calculating and finding the minimum RMSE value in each step.The minimum RMSE was recorded in the iterative process to find the corresponding optimal hyperparameters.The time cost of the algorithm was appropriately improved, while the prediction accuracy of the algorithm was enhanced.

Multi-Objective PSO Optimization
The MOPSO is recognized as an efficient tool for optimizing the proxy model of study objects, which was proposed by Coello and Lechuga [46] based on the regular PSO.Furthermore, the particles are updated according to Equations ( 1) and ( 2).The conventional PSO can only set one objective function (fitness function) to obtain a unique optimal solution.In contrast, the MOPSO can handle multiple functions.Therefore, MPSO can obtain a set of optimal solutions, i.e., Pareto solutions [35,47].MOPSO can quickly converge to the Pareto front for its practical searchability.Then, the optimal values (function value and variables) are obtained by extracting the Pareto solution [48].
PPV and the bulk rate (Br) were taken as the two optimized objects to control the ground vibration in a sensitive area and ensure the blasting fragmentation effect.Therefore, the developed PSO-LSSVM (PPV) is called in MOPSO as one of the proxy models.Furthermore, an empirical model, which combined the Kuznetsov model with the Rosin-Rammler model (Kuz-Ram model), was considered as the proxy model for Br.The Kuz- In the model's training, the particle velocity was updated by Equation (1), and the particle position was updated by Equation (2).The individual and global fitness were updated by calculating and finding the minimum RMSE value in each step.The minimum RMSE was recorded in the iterative process to find the corresponding optimal hyperparameters.The time cost of the algorithm was appropriately improved, while the prediction accuracy of the algorithm was enhanced.

Multi-Objective PSO Optimization
The MOPSO is recognized as an efficient tool for optimizing the proxy model of study objects, which was proposed by Coello and Lechuga [46] based on the regular PSO.Furthermore, the particles are updated according to Equations ( 1) and ( 2).The conventional PSO can only set one objective function (fitness function) to obtain a unique optimal solution.In contrast, the MOPSO can handle multiple functions.Therefore, MPSO can obtain a set of optimal solutions, i.e., Pareto solutions [35,47].MOPSO can quickly converge to the Pareto front for its practical searchability.Then, the optimal values (function value and variables) are obtained by extracting the Pareto solution [48].
PPV and the bulk rate (Br) were taken as the two optimized objects to control the ground vibration in a sensitive area and ensure the blasting fragmentation effect.Therefore, the developed PSO-LSSVM (PPV) is called in MOPSO as one of the proxy models.Furthermore, an empirical model, which combined the Kuznetsov model with the Rosin-Rammler model (Kuz-Ram model), was considered as the proxy model for Br.The Kuz-Ram model of the Br is described as Equation ( 12) [49].Moreover, the process of two-objective optimization by MOPSO is shown in Figure 5.After initializing the population, PSO-LSSVM (PPV) and Kuz-Ram (Br) are invoked through function calls.Then, the non-dominated solution is determined by comparing the fitness values and updated to the repository set.Meanwhile, the local optimum and the global optimum are selected in this update.The particles move according to Equations ( 1) and ( 2) until the maximum number of iterations.Br = e −( x 0 xe ) where x is the average block size, β is the uniformity factor, x e is the characteristic block size, and Br is the ratio of ores with a block size greater than x 0 .
repository set.Meanwhile, the local optimum and the global optimum are selected in this update.The particles move according to Equations ( 1) and ( 2) until the maximum number of iterations.
where x is the average block size, β is the uniformity factor, xe is the characteristic block size, and Br is the ratio of ores with a block size greater than x0.

PSO-LSSVM Model Building
In this study, intelligent algorithms were built based on the MATLAB software.Moreover, the RMSE was used as a performance index for modeling.The performance of the prediction and RMSE were negatively correlated.The smaller the RMSE value, the better the prediction effect is.In order to obtain the best prediction performance of the

PSO-LSSVM Model Building
In this study, intelligent algorithms were built based on the MATLAB software.Moreover, the RMSE was used as a performance index for modeling.The performance of the prediction and RMSE were negatively correlated.The smaller the RMSE value, the better the prediction effect is.In order to obtain the best prediction performance of the PSO-LSSVM model, the maximum number of iterations (I max ) and the number of particles in the swarm (n pop ) of PSO were determined according to the experimental method, and the other parameters were taken empirically.In addition, the inertia weight coefficient (ε) was iterated according to formula ε − (ε − 0.4)T/k.As shown in Figure 6, after 100 iterations, there was no significant change in the fitness curve.Among them, the fitness curve with n pop = 20 had the smallest RMSE at the final steady state and reached stability at approximately 70 iterations.Therefore, I max was taken as 100, and n pop was taken as 20.The initial parameters of the model are shown in Table 3.
in the swarm (npop) of PSO were determined according to the experimental method, and the other parameters were taken empirically.In addition, the inertia weight coefficient (ε) was iterated according to formula . As shown in Figure 6, after 100 iterations, there was no significant change in the fitness curve.Among them, the fitness curve with npop = 20 had the smallest RMSE at the final steady state and reached stability at approximately 70 iterations.Therefore, Imax was taken as 100, and npop was taken as 20.The initial parameters of the model are shown in Table 3.Therefore, the optimal values of  and  were 17.315 and 0.917, respectively, when PPV was used as the output indicator.Based on the optimized hyperparameters and training set, the regression fit of the PSO-LSSVM model was obtained, as shown in Figure 7.It can be seen that the PSO-LSSVM model fit the PPV better.Therefore, the optimal values of γ and σ 2 were 17.315 and 0.917, respectively, when PPV was used as the output indicator.Based on the optimized hyperparameters and training set, the regression fit of the PSO-LSSVM model was obtained, as shown in Figure 7.It can be seen that the PSO-LSSVM model fit the PPV better.

Parameter Configuration of the Models Used for Comparison
In order to validate the predictive performance of the proposed joint model (PSO-LSSVM), other hybrid models (GA-BP) and single models (LSSVM and BP) were used for comparison purposes.BP neural networks were proposed by Rumelhar et al. [50].Their structure consists of an input layer, a hidden layer, and an output layer.In this study, a

Parameter Configuration of the Models Used for Comparison
In order to validate the predictive performance of the proposed joint model (PSO-LSSVM), other hybrid models (GA-BP) and single models (LSSVM and BP) were used for comparison purposes.BP neural networks were proposed by Rumelhar et al. [50].Their structure consists of an input layer, a hidden layer, and an output layer.In this study, a three-layer neural network was used.BP has the disadvantage of high workload and reliance on empirical values when assigning the connection weights and thresholds between neurons [51].GA belongs to the metaheuristic algorithm [33].Using the algorithm function of GA (e.g., selection, crossover, and mutation) to optimize BP is considered adequate [52].
According to the method in [53], the LSSVM was trained three times for the training datasets with three different sets of γ and σ 2 , i.e., (50, 1.9), (10, 0.6), and (4, 0.04).The one with the smallest RMSE was taken to participate in the comparison.The research experience in [52] was introduced to show that the BP and GA-BP methods were applied in the prediction of PPV.The accuracy of the training model was set to 0.0001; the learning rate was 0.1; the maximum number of iterations was 1000.The weight threshold was set in the range of [−2, 2] for the BP neural network without optimization.The multiple predictions (10 times) of BP were performed, and the group with the smallest RMSE was taken to participate in the discussion.All the RMSEs of LSSVM and BP in the modeling process are shown in Figure 8.

Parameter Configuration of the Models Used for Comparison
In order to validate the predictive performance of the proposed joint model (PSO-LSSVM), other hybrid models (GA-BP) and single models (LSSVM and BP) were used for comparison purposes.BP neural networks were proposed by Rumelhar et al. [50].Their structure consists of an input layer, a hidden layer, and an output layer.In this study, a three-layer neural network was used.BP has the disadvantage of high workload and reliance on empirical values when assigning the connection weights and thresholds between neurons [51].GA belongs to the metaheuristic algorithm [33].Using the algorithm function of GA (e.g., selection, crossover, and mutation) to optimize BP is considered adequate [52].
According to the method in [53], the LSSVM was trained three times for the training datasets with three different sets of γ and σ 2 , i.e., (50, 1.9), (10, 0.6), and (4, 0.04).The one with the smallest RMSE was taken to participate in the comparison.The research experience in [52] was introduced to show that the BP and GA-BP methods were applied in the prediction of PPV.The accuracy of the training model was set to 0.0001; the learning rate was 0.1; the maximum number of iterations was 1000.The weight threshold was set in the range of [−2, 2] for the BP neural network without optimization.The multiple predictions (10 times) of BP were performed, and the group with the smallest RMSE was taken to participate in the discussion.All the RMSEs of LSSVM and BP in the modeling process are shown in Figure 8.

Two-Objective Optimization Model
PSO-LSSVM (PPV) for the sensitive area (c) and Kuz-Ram (Br) in the chosen blasting area (a) were used as a proxy model in the MOPSO.After the site survey, the R of the sensitive area (c) and blasting area was 200 m, and H was 90 m.According to the actual ore, the bulk size x 0 was 80 cm.Therefore, the two-objective optimization strategy required considering five variables, i.e., W, PF, S, B, and H 0 .Therefore, the number of decision variables (n Var ) was five.The parameter set is shown in Table 4.
In the MOPSO, all the data were treated the same way as the PPV prediction model, i.e., the problem of different magnitudes was avoided by means of normalization and anti-normalization.The data were manipulated by Equation (13).For example, the random variable (X) assigned to the PSO-LSSVM model was within 0 and 1, while the predicted PPV was anti-normalized to the original magnitudes.Besides, the randomly generated variables were anti-normalized to the values with the range of [x min , x max ] of measured variables and quoted by Equation ( 13) to obtain the bulk rate.Therefore, the optimization was conducted on 30 blasting experiences.

X =
x − x min x max − x min (13)

Comparison and Evaluation of Prediction Models
To verify the performance of the developed prediction model (PSO-LSSVM), the same training datasets and testing datasets were utilized.Furthermore, the root-mean-squared error (RMSE), mean absolute error (MAE), and correlation coefficient (r) were selected as the assessment indicators.The RMSE and MAE were negatively correlated with the performance of the models, while on the contrary for r.The efficiencies of PSO-LSSVM, GA-BP, the unoptimized LSSVM model, and BP were compared by calculating and evaluating the performance indicators.As shown in Table 5, it can be seen that the hybrid models (PSO-LSSVM and GA-BP) outperformed the single prediction models (LSSVM and BP).In the single model, LSSVM seemed to exhibit a slightly better predictive performance than BP, with smaller error metrics (RMSE and MAE) and larger fit coefficients (r).Importantly, Table 5 shows that PSO-LSSVM had a more reliable performance compared with the remaining models.Additionally, a PPV profile between the predicted and the measured for training and testing is shown in Figures 9 and 10.Based on the graphical representation of the regression fit in Figures 9 and 10, the predicted values of PSO-LSSVM and GA-BP had a tendency to be close to the fit line.For instance, the predicted values of the single model had more discrete points, with a lower r.Consequently, the PSO-LSSVM model was closest to the measured values with an r of 0.983 and 0.965 for the training and test sets.

Optimization of Blast Design Parameters
In this section, a strategy of two-objective optimization is presented, i.e., using MOPSO to format the Pareto front by assigning and solving PSO-LSSVM (PPV) and Kur-Ren (Br).Therefore, a stable Pareto front was obtained after 100 iterations.As shown in Figure 11, the Pareto solution meant a set of optimized function values, which were limited to each other.Spontaneously, the corresponding independent variable could be extracted from the repository set.Furthermore, Points A and C were the two endpoints of the Pareto-optimal solution set, which corresponded to the two cases of the minimum PPV and Br, respectively.Point B was the middle point, which had a more balanced consideration between the two objectives.
sion fit in Figures 9 and 10, the predicted values of PSO-LSSVM and GA-BP had a tendency to be close to the fit line.For instance, the predicted values of the single model had more discrete points, with a lower r.Consequently, the PSO-LSSVM model was closest to the measured values with an r of 0.983 and 0.965 for the training and test sets.

Optimization of Blast Design Parameters
In this section, a strategy of two-objective optimization is presented, i.e., using MOPSO to format the Pareto front by assigning and solving PSO-LSSVM (PPV) and Kur-Ren (Br).Therefore, a stable Pareto front was obtained after 100 iterations.As shown in

Optimization of Blast Design Parameters
In this section, a strategy of two-objective optimization is presented, i.e., using MOPSO to format the Pareto front by assigning and solving PSO-LSSVM (PPV) and Kur-Ren (Br).Therefore, a stable Pareto front was obtained after 100 iterations.As shown in Thus, the three classical solutions (A, B, and C) are shown in Table 6.It is worth noting that the transformation of PPV in Sensitive Region c was weaker than the Br of the blasting fragmentation when the parameters changed.In order to validate the optimization results, a blasting scheme was designed using the blasting parameters extracted from Point B. Then, limited blasting field tests were carried out, and a total of three set dates were recorded, finally.Therefore, considering the B solution in the blasting scheme, the measured PPV values for Sensitive Area c were 0.752 cm/s, 1.00 cm/s, and 1.01 cm/s, and the Brs of the blasting fragmentation bulk were 11.02%, 9.06%, and 8.42%.Consequently, the actual measured values were similar to the optimized values, indicating that the optimization results were convincing.
Figure 11, the Pareto solution meant a set of optimized function values, which were limited to each other.Spontaneously, the corresponding independent variable could be extracted from the repository set.Furthermore, Points A and C were the two endpoints of the Pareto-optimal solution set, which corresponded to the two cases of the minimum PPV and Br, respectively.Point B was the middle point, which had a more balanced consideration between the two objectives.Thus, the three classical solutions (A, B, and C) are shown in Table 6.It is worth noting that the transformation of PPV in Sensitive Region c was weaker than the Br of the blasting fragmentation when the parameters changed.In order to validate the optimization results, a blasting scheme was designed using the blasting parameters extracted from Point B. Then, limited blasting field tests were carried out, and a total of three set dates were recorded, finally.Therefore, considering the B solution in the blasting scheme, the measured PPV values for Sensitive Area c were 0.752 cm/s, 1.00 cm/s, and 1.01 cm/s, and the Brs of the blasting fragmentation bulk were 11.02%, 9.06%, and 8.42%.Consequently, the actual measured values were similar to the optimized values, indicating that the optimization results were convincing.

Conclusions
Vibration hazard control is the same as rock fragmentation and is the goal pursued in open-pit blasting.Meanwhile, reliable and accurate blasting vibration prediction models and referenceable blasting design parameters are important for the control of blastinginduced ground vibration.To accurately predict and control blasting vibrations, this paper proposed the PSO-LSSVM prediction model and the MOPSO optimization model.The following conclusions can be drawn: 1.The PSO-LSSVM model outperformed GA-BP, unoptimized LSSVM, and BP in terms of the prediction efficiency and accuracy of PPV, with better prediction error (RMSE, MAE) and goodness of fit (r).More specifically, the hyperparameter-optimized PSO-LSSVM reduced the prediction error (RMSE, MAE) by more than 55% and improved

Conclusions
Vibration hazard control is the same as rock fragmentation and is the goal pursued in open-pit blasting.Meanwhile, reliable and accurate blasting vibration prediction models and referenceable blasting design parameters are important for the control of blastinginduced ground vibration.To accurately predict and control blasting vibrations, this paper proposed the PSO-LSSVM prediction model and the MOPSO optimization model.The following conclusions can be drawn: 1.
The PSO-LSSVM model outperformed GA-BP, unoptimized LSSVM, and BP in terms of the prediction efficiency and accuracy of PPV, with better prediction error (RMSE, MAE) and goodness of fit (r).More specifically, the hyperparameter-optimized PSO-LSSVM reduced the prediction error (RMSE, MAE) by more than 55% and improved the goodness of fit (r) by more than 10%, compared with the least-efficient one (unoptimized BP).

2.
With the multi-objective optimization function of PSO, a relatively optimal blasting design was found based on existing blasting tests.When considering the PPV of the sensitive area (R = 100 m, H = 90 m) and the Br, the blasting design parameters for Zhoushan mine were taken as W = 173.68kg, PF = 0.38 kg/m 3 , S = 5.96 m, B = 3.19 m, and H 0 = 1.91 m, corresponding to a PPV of 0.893 cm/s and Br of 9.8%.Therefore, it can provide engineers with a more specific and reliable blasting vibration control reference.

3.
In order to predict and control the PPV more accurately, more experimental data and influencing factors should be considered.Furthermore, reliable agent models are a requirement for more accurate multi-objective optimization results.

Figure 3 .
Figure 3. Monitoring point arrangement (a is the area for the verification of the optimized blasting parameters scheme).

Figure 3 .
Figure 3. Monitoring point arrangement (a is the area for the verification of the optimized blasting parameters scheme).

Figure 8 .
Figure 8. RMSEs of LSSVM and BP in the modeling.

Figure 8 .
Figure 8. RMSEs of LSSVM and BP in the modeling.

Figure 9 .
Figure 9. Regression curve of predicted and measured PPV for training by different methods.

Figure 10 .
Figure 10.Regression curve of predicted and measured PPV for testing by different methods.

Figure 9 .
Figure 9. Regression curve of predicted and measured PPV for training by different methods.

Figure 9 .
Figure 9. Regression curve of predicted and measured PPV for training by different methods.

Figure 10 .
Figure 10.Regression curve of predicted and measured PPV for testing by different methods.

Figure 10 .
Figure 10.Regression curve of predicted and measured PPV for testing by different methods.

Table 1 .
Mechanical parameters of rocks in the test area.

Table 1 .
Mechanical parameters of rocks in the test area.

Table 1 .
Mechanical parameters of rocks in the test area.

Table 2 .
Input and output parameter characteristics.

Table 2 .
Input and output parameter characteristics.

n Var * c 1 c 2 ε 0 * n pop pop max * pop min * v max * v min * I max
* n Var : number of decision variables; ε 0 : initial inertia weight coefficient; pop max : optimization parameter maximum value; pop min : optimization parameter minimum value; v max : maximum update speed; v min : maximum update speed.

Table 4 .
Parameters setting of PSO.n rep : repository size; n Grid : the grid dimension of the decision variable; gamma: elimination factor of nondominated solutions; beta: selection factor for non-dominated solutions; mu: mutation rate. *

Table 5 .
Performance indicators of the prediction models.

Table 6 .
Three sets of classical solutions for the Pareto front.

Table 6 .
Three sets of classical solutions for the Pareto front.