Hindawi

Shock and Vibration

Volume 2021, Article ID 2847180, 18 pages
[Link]

Research Article
Prediction of BlastInduced Ground Vibration (BIGV) of Metro
Construction Using Difference Evolution AlgorithmOptimized
Gaussian Process (DE-GP)

Tengfei Jiang , Annan Jiang , Shuai Zheng , and Mengfei Xu

College of Transportation Engineering, Dalian Maritime University, Dalian 116026, China

Correspondence should be addressed to Annan Jiang; jiangannan@[Link]

Received 2 June 2021; Revised 21 August 2021; Accepted 20 September 2021; Published 4 October 2021

Academic Editor: Hung Nguyen-Xuan

Copyright © 2021 Tengfei Jiang et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Rock blasting often has an irreversible impact on the surrounding environment and threatens the safety of life and property.
Therefore, accurate prediction of blast-induced ground vibration (BIGV) is a prerequisite for safe construction. In view of the fact
that traditional blasting peak particle velocity (PPV) empirical formulas cannot be accurately predicted, this study selected 88 sets
of blasting monitoring data, based on distance from the blast-face, maximum charge per delay, total charge, hole depth, spacing,
burden, stemming length, and powder factor being used as input variables and PPV being used as output variable to characterize
BIGV. First, a nonlinear mapping relationship between input variables and output variable is established through the Gaussian
process (GP). The differential evolution algorithm (DE) is used to optimize the hyperparameters σ f, σ n, and l of the GP, and a
blasting PPV model based on the DE-GP is constructed. The proposed model is compared with the empirical formulas, least
square support vector machine (LSSVM), artificial neural network (ANN), and GP model, and its prediction performance is
evaluated by statistical indicators such as root mean square error (RMSE). Finally, the cosine amplitude method (CAM) is used to
analyze the sensitivity of blasting parameters. The results show that the DE-GP algorithm for blasting vibration velocity prediction
has higher precision and accuracy, which is significantly better than other models, and is the closest to the measured PPV.
Distance from the blast-face, total charge, and maximum charge per delay have a greater impact on the prediction of PPV, while
stemming length and powder factor have a smaller impact on the prediction of PPV. The DE-GP model proposed by this research
has certain reference value for the prediction and control of PPV in blasting construction.

1. Introduction rock mass, and there are many influence factors of PPV
[9, 10]. How to accurately predict PPV has become the
At present, the widespread application of blasting tech- primary issue in the field of blasting construction safety
nology has penetrated into many areas of the national technology and scientific research [11].
economy. Because of its high efficiency, economy, and speed, In the past blasting engineering, the most commonly
it has long been favored in the field of engineering con- used method for PPV prediction is the empirical formulas.
struction [1–3]. However, blasting also produces a series of The empirical formulas are used to perform regression
harmful effects, mainly including blast-induced ground analysis on the measured data of blasting vibration. It mainly
vibration (BIGV), blasting flying rocks, noise, shock wave, depends on distance from the blast-face and maximum
etc. [4, 5]. charge per delay [12, 13]. However, the prediction accuracy
BIGV adversely affects surrounding rock masses and of the empirical formulas is generally low, and the regression
nearby structures and even causes damage [6]. Most calculation effect for some blasting projects is not ideal. As a
countries in the world use peak particle velocity (PPV) as the result, the blasting vibration attenuation formulas regressed
measurement index of BIGV [7, 8]. However, the propa- by the empirical formulas is prone to excessive deviation in
gation medium of blasting vibration waves is heterogeneous the guidance of site construction [14, 15].
2 Shock and Vibration

Because the prediction of PPV is a complex nonlinear distance from the blast-face are the two main factors af-
problem, researchers generally use intelligent algorithm fecting PPV. Therefore, many domestic and foreign scholars
models with powerful processing capability to predict PPV, have proposed different PPV formulas for predicting BIGV
such as artificial neural network (ANN), genetic algorithm [29–34]. In this study, several commonly used empirical
(GA), particle swarm optimization (PSO), and support formulas to characterize the attenuation of blasting vibration
vector machine (SVM)[1, 16–27], as shown in Table 1. were selected, as shown in Table 2. In the formula, K, α, β,
It can be seen from Table 1 that intelligent algorithms and n are parameters related to rock properties, blasting
have been widely used in the prediction of blasting vibration operation parameters, and terrain conditions. We often use
velocity. Research on the combination of various artificial the field blasting vibration measured data obtained in the
intelligence methods has gradually increased, and the hybrid actual engineering as the basis and perform regression
artificial intelligence algorithms have high accuracy. Most of analysis according to the empirical formulas in Table 2 to
the above research works are based on a hybrid model obtain the values of the unknown parameters K, α, β, and n.
composed of ANN algorithm. However, the ANN algorithm According to the blasting prediction empirical formulas
has some disadvantages such as being easy to fall into a local in Table 2, the site constants suitable for this project are fitted
optimum, requiring a large number of parameters during by using the actual measurement data of the construction
calculation, and slow training speed. At present, optimiza- site, and the blasting vibration velocity prediction formulas
tion methods such as balancing composite motion opti- considering the maximum charge per delay and the distance
mization (BCMO) [28] and Gaussian process (GP) have from the blast-face is constructed. After that, the corre-
begun to be gradually applied. GP has the advantages of sponding PPV can be obtained according to the maximum
good adaptability, strong generalization ability, simple charge per delay and the distance from the blast-face of a
implementation process, and adaptive acquisition of certain blasting construction.
hyperparameters for processing nonlinear regression
problems. At the same time, the differential evolution al-
gorithm (DE) has strong global optimization capabilities and 2.2. Blasting Vibration Velocity Prediction Model Based on
can quickly approach the global optimum. Therefore, it can DE-GP. This section selects the blasting parameters
achieve rapid acquisition of GP hyperparameters by intro- according to the influencing factors of blasting vibration.
ducing the DE into the GP model and form a new model of Through a brief description of the GP, a nonlinear prediction
PPV prediction based on DE-GP, which will realize the model of blasting vibration velocity based on the GP is
effective control of BIGV. formed. Limited by the difficulty of selecting hyper-
Aiming at the problem of unpredictable PPV in the parameters for the GP, the DE algorithm is introduced to
bustling downtown area, this study uses DE to optimize GP, optimize the GP and applied to the blasting prediction of
constructs a model of mixed algorithm DE-GP predicting rock foundation pits. Then, the principle, frame structure,
PPV, and analyzes the accuracy of the model. Finally, the and solution process are described.
sensitivity analysis method is used to discuss the sensitivity
of blasting influencing factors, which provides a certain basis
for the optimization of blasting parameters. The main part of 2.2.1. Selection of Blasting Prediction Parameters. In the
this article is divided as follows: The second part introduces actual blasting of rock foundation pits, there are many
the methods used in blasting prediction, including empirical factors that affect the blasting effect, and various different
formulas and the intelligent algorithms studied in this ar- factors also cause different degrees of impact. In summary,
ticle. The third part introduces the engineering situation and there are four main aspects: explosive performance, rock
the acquisition of blasting vibration velocity data. The fourth properties, explosive and rock-related factors, and blasting
part uses conventional empirical formulas and the proposed technology. Combining the actual engineering situation and
intelligent algorithm to predict PPV. The fifth part analyzes related construction experience of the Labor Park Station of
the influence of the hyperparameters of the GP model, the Dalian Metro Line 5, the distance from the blast-face, the
selection of DE algorithm parameters, the comparison of the maximum charge per delay, the total charge, the hole depth,
predictive capabilities of different algorithms, and the sen- the spacing, the burden, the stemming length, and the
sitivity between blasting influencing factors. Finally, the powder factor are used as the input parameters of the
sixth part draws a conclusion. blasting vibration velocity prediction model, and the peak
particle velocity at the monitoring point is used as the output
2. Methods Used parameter. The above indicators can completely describe the
blasting effect of the rock foundation pit of the Suoyuwan
2.1. Conventional Prediction Methods of Blasting. PPV is the South Station of Dalian Metro Line 5 and fully reflect the
most important indicator to measure the impact of blasting blasting effect, so as to realize the blasting prediction of PPV.
on the surrounding environment. There are many factors
that affect blasting vibration, and the relationship between
the factors is also intricate. It is practically impossible to 2.2.2. Nonlinear Prediction Model of PPV Based on GP.
establish a PPV prediction formula that considers all in- GP has good adaptability to the processing of complex
fluence factors. A large number of blasting engineering classification and regression problems with high dimen-
monitors show that the maximum charge per delay and the sionality, small samples, and nonlinearity. Therefore, it is
Shock and Vibration 3

Table 1: Prediction of PPV by intelligent algorithms.

References Models Input parameters Number of samples Coefficient of determination (R2)
Armaghani et al. [1] ANFIS; ANN Qm, D 109 0.973 (ANFIS); 0.949 (ANN)
0.911 (ANN);
Hajihassani et al. [16] ANN; ICA-ANN B/S, SL, Qm, E, Vp, D 95
0.976 (ICA-ANN)
Hajihassani et al. [17] SVM Qm, D 80 0.957
Amiri et al. [18] ANN; ANN-KNN Qm, D 75 0.82 (ANN); 0.88 (ANN-KNN)
Ghasemi et al. [19] ANFIS-PSO; SVM B, S, SL, Qm, D 120 0.950 (ANFIS-PSO); 0.924 (SVM)
Ebrahimi et al. [20] ANN B, S, SL, H, PF 34 0.78
Taheri et al. [21] ANN; ABC-ANN Qm, D 99 0.91 (ANN); 0.95 (ABC-ANN)
Hasanipanah et al. [22] PSO Qm, D 80 0.9019
Faradonbeh et al. [23] GEP B, S, SL, d, H, PF, Qm, D 115 0.876
Armaghani et al. [24] ICA Qm, D 73 0.9334
Zhang et al. [25] PSO-XGBoost B, S, PF, Qm, D 175 0.968
Tian et al. [26] GA-linear; GA-power Qm, D 85 0.961 (GA-linear); 0.967 (GA-power)
Rezaeineshat et al. [27] ICA-ANN B, S, Qm, D, RQD 112 0.9002
Distance from the blast-face (D); maximum charge per delay (Qm); total charge (Qt); hole depth (H); spacing (S); burden (B); Stemming length (SL); powder
factor (PF); burden-to-spacing ratio (B/S); Young’s modulus (E); p-wave velocity (Vp ); rock quality designation (RQD); hole diameter (d); coefficient of
determination (R2); imperialist competitive algorithm (ICA); adaptive neuro-fuzzy inference system (ANFIS); K-nearest neighbors (KNN); gene expression
programming (GEP); artificial bee colony (ABC); extreme gradient boosting machine (XGBoost).

Table 2: Conventional empirical formulas for blasting prediction.

The name of the empirical formulas Formulas
􏽰���� − α
United States Bureau of Mines (USBM) [29] PPV � K(R/ Qmax )
􏽰��������
Langefors–Kihlstrom (LK) [30] PPV � K( Qmax /R2/3 )α
− α β
General predictor [31] PPV � KR􏽰 Qmax
���� − α
Ambraseys–Hendron (AH) [32] PPV � K(R/ 3 Q√ ���) α
max
3
Bureau of Indian Standard (BIS) [33] R2 )
PPV � K(Qmax􏽰/ ����
CMRI predictor (CMRI) [34] PPV � n + K(R/ Qmax )− 1
PPV is peak particle velocity (cm/s); R is the distance from blast-face to monitoring point (m); Qmax is the maximum charge per delay (kg); and K, α, β, and n
are the site constants.

used to establish the blasting PPV prediction model of the Further solving to obtain the mean and variance of y∗ :
rock foundation pit.
− 1
For a set of known blasting influencing factors y􏽢 x∗ 􏼁 � kT x∗ 􏼁􏼐K + σ 2n I􏼑 y, (3)
X(x1 , x2 , x3 , x4 ) and corresponding blasting vibration ve-
locity y, the blasting PPV learning samples set D � (X, y) − 1
σ􏽢2 x∗ 􏼁 � k x∗ , x∗ 􏼁 − kT x∗ 􏼁􏼐K + σ 2n I􏼑 k x∗ 􏼁. (4)
can be established, and the prior distribution of blasting PPV
y can be established through the GP as follows: In the formula, k(x∗ ) is the covariance matrix of order
Τ
y ∼ N􏼐0, K + σ 2n I􏼑. (1) n × 1, k(x∗ ) � K(X, x∗ ), namely 􏼈k(x1 , x∗ ), . . . , k(xn , x∗ )􏼉
The nonlinear prediction model of the blasting PPV is
In the formula, K � K(X, X) is a symmetric positive described by the square exponential covariance function:
definite covariance matrix of order n × n; and any item Kij in
the K � K(X, X) matrix measures the correlation between xi 1 2
k􏼐xp , xq 􏼑 � σ 2f exp􏼠− 2
2 􏼐xp , xq 􏼑 􏼡 + σ n δpq . (5)
and xj ; 0 is the mean; σ 2n is the variance; and I is the identity 2l
matrix.
From equation (1), the joint GP prior distribution In the formula, σ f , l, and σ n are hyperparameters. σ 2f is
composed of n training sample output y and one test sample the signal variance of the kernel function, which is used to
output y∗ is derived as follows: control the degree of local correlation; l is the hyper-
parameter for correlation measurement, the larger the value,
y
⎝0, ⎡ K(X, X) + σ 2n I K X, x∗ 􏼁 ⎠ the smaller the correlation between input and output; σ 2n is
􏼢 ∗ 􏼣 ∼ N⎛ ⎣ ⎤⎦⎞. (2)
y K X, x∗ 􏼁 k x∗ , x∗ 􏼁 the variance of noise; and δpq is the symbol of Kronecker,
when p � q, δpq � 0, otherwise δpq � 1
In the formula, K(X, x∗ ) is the n × 1-order covariance Hyperparameters can be obtained by maximizing the
matrix of the blasting test input point x∗ and all input points log-likelihood of the training samples (maximum posterior
X in the training set; and k(x, x∗ ) is the covariance of the distribution). The logarithmic form of the likelihood
blasting test input point x∗ itself. function is as follows:
4 Shock and Vibration

L � lgp(y | X) magnitude but also intuitively describe the data. Compared

with the MAE evaluation index, RMSE is the square root of
(6)
1 − 1 1 􏼌􏼌 􏼌􏼌 n the sum of squares of errors, which means that the single
� − yT 􏼐K + σ 2n I􏼑 y − lg􏼌􏼌K + σ 2n I􏼌􏼌 − lg(2π). error has a quadratic increase. It enlarges the gap between
2 2 2
errors and makes the error evaluation result more obvious.
Obtain the optimal hyperparameters by formula (6), and Therefore, RMSE is used as an evaluation index for the
substitute formula (3) to obtain the blasting PPV prediction prediction effect of the GP model. The realization process is
value y∗ corresponding to the blasting test sample x∗ . shown in Figure 1.
According to the establishment steps of the above model, the For the aforementioned GP blasting PPV prediction
corresponding calculation program is compiled using model, the specific steps of using DE optimization are as
MATLAB language. So far, the modeling and regression follows:
process of the nonlinear prediction model for blasting PPV
of rock foundation pit based on the GP has been completed. Step 1. Initialize the DE algorithm, and set the DE algo-
rithm’s mutation factor F, cross factor CR, population size
NP, and difference strategy.
2.2.3. Hybrid DE-Optimized GP and Its Architecture Design.
DE is a genetic algorithm based on real number coding. Step 2. Choose σ f , l, and σ n in GP as the optimization target,
Compared with the selection operation of GA algorithm, DE and randomly generate the initial population in the opti-
updates the population according to a one-to-one elimi- mization interval.
nation mechanism and uses directional perturbation on
individuals to achieve a reduction in the value of individual Step 3. Use the parameters expressed by each individual in
functions. In terms of applicability, the GA algorithm and the initial population to perform GP training separately.
PSO algorithm are mainly used to solve continuous opti-
mization problems, and the ACO algorithm is mainly used Step 4. Evaluate the fitness of each individual in the pop-
to solve discrete optimization and combinatorial optimi- ulation. The evaluation index is the root mean square error
zation problems. However, DE not only does not use the of the prediction result of the test sample of the trained GP
gradient information of the objective function but also has model:
no requirements on the derivability and continuity of the 􏽶�������������
􏽴
function. So, DE has strong applicability. In terms of pa- 1 p
N
2
rameter setting, DE has two main parameters to be adjusted, F0 � 􏽘 v − v′i􏼁 . (7)
and the parameter setting has less obvious influence on the Np i�1 i
result, so it is easier to use. Compared with the GA algorithm
and the PSO algorithm, there are many parameters, and In the formula, vi and v′i are the actual vibration velocity
different parameter settings have a greater impact on the and the predicted vibration velocity, Np is the number of test
final result. Therefore, in actual use, it needs to be adjusted samples, and F0 is the individual fitness value.
continuously, which increases the difficulty of using the
algorithm. When solving multivariable function optimiza- Step 5. For the minimum fitness value Fmin and expected
tion problems, the GA algorithm converges slowly or even response value Fbest of the current population, if the current
hard for high-dimensional problems, while the PSO algo- individual satisfies Fbest > Fmin , then enter Step 10; otherwise,
rithm and DE algorithm have the advantages of fast con- enter Step 6.
vergence speed and high accuracy. But the PSO algorithm
and ACO algorithm are easy to fall into a local optimal Step 6. Perform a mutation operation on the current
solution. In addition, DE considers the correlation between population, and combine the difference values between
multiple variables to a certain extent and has great advan- several individuals in the population with the individuals to
tages in the problem of variable coupling. be mutated according to certain rules. Then, obtain the
The study found that the regression results of the mutation vector of each individual in the population.
nonlinear prediction model for blasting PPV of rock
foundation pit based on GP are significantly affected by the Step 7. Perform crossover operation on the current pop-
hyperparameters σ f , l, and σ n . In this study, DE is intro- ulation; crossover the mutation vector and its corresponding
duced into GP. It can quickly approach the global optimi- target vector to obtain the test vector.
zation with the help of DE’s strong global optimization
ability. This can realize the optimal selection of three Step 8. Select the current population, use the fitness function
hyperparameters and form a new model for controlling of the optimized problem to compare each test vector with
BIGV. RMSE is used as the evaluation index of the pre- its corresponding target vector, and put the better one into
diction results of the GP model, which is used to describe the the next-generation population.
optimization effect of the GP hyperparameters. Compared
with the MSE evaluation index, the calculation result of Step 9. Judge whether the maximum number of iterations is
RMSE can not only keep the error evaluation result con- reached, if it reaches the maximum number of iterations,
sistent with the problem described by GP in the order of enter Step 10; otherwise, return to Step 4.
Shock and Vibration 5

Begin
northeast side of the site is the Unicom Building, the
northwest side is the stadium under construction and Labor
Park, the east and southwest sides of the site are the Dalian
DE algorithm initialization No. 9 Middle School and the Dalian No. 24 Middle School,
the southeast side is the residential area, and the south side is
the Shengli East Road viaduct. There is a 2 m × 4 m culvert on
Generate an initial population
representing GP hyperparameters
the west side of the site, and the schematic diagram of the
project is shown in Figure 2. The minimum distance between
the main structure of the station and the library of the 24th
Learning GP model Middle School was 5.9 m, and the minimum distance from
samples training the No. 26 building in the residential area was 7.2 m. Jiefang
Road is the main traffic arterial road in Dalian. This road has
a large traffic volume, and many pipelines are laid down. It is
Fitness Selection a dense area of pedestrians and vehicles.
evaluation

DE search for optimal

The main part of the station was constructed by the PBA

hyperparameters
construction method, and two construction foundation pits
Judgment of minimum were set up. The foundation pits were strictly in accordance
GP model training

Crossover
and prediction

fitness of population with the principle of layered, segmented, symmetrical,

balanced, and timely excavation, and supports were erected
in time. The foundation pit adopted a combined support
Fbest>Fmin
Mutation
system with retaining piles with a diameter of 600 mm and
an interval of 1300 mm∼1500 mm; the outer diameter of the
steel supports was 609 mm, and the pipe wall thickness was
No 16 mm. The initial support system of the concealed exca-
Whether to reach the maximum
evolutionary algebra vation section was advanced small pipe grouting to reinforce
the ground or deep hole grouting, grid steel frame, and mesh
Yes shotcrete support. Because most of the surrounding rock
during construction was in strongly weathered slate and
Obtain the GP hyperparameters moderately weathered slate, in order to ensure the safety of
corresponding to the optimal fitness value ground buildings, underground structures, and the sur-
rounding environment during drilling and blasting, we
Establishment of DE-GP regression model should adopt microseismic controlled blasting construction
of blasting vibration velocity technology, low-detonation explosives, and millisecond
detonators. According to the type of surrounding rock, the
integrity of the rock, the development of joints, the bedding,
End and the location of excavation, different blasting parameters
were adopted and adjusted according to the actual blasting
Figure 1: Flow chart of DE-GP blasting PPV regression model.
effect during the blasting construction.
The foundation pit of Labor Park adopted a triangular
Step 10. Take the σ f , l, and σ n represented by the minimum hole pattern and vertical drilling. The detonation adopted a
fitness individual as the initial parameters of the GP model, millisecond-delay nonelectrically conductive detonation
so complete the GP optimization process. system, and the method of detonation outside the hole with
slight difference, cluster connection method, five holes with
3. Data Used one sound. Used ms-8 in the hole, ms-4 for the connection
between the surface holes, and ms-5 for the connection
3.1. Description of Study Area. Dalian Metro Line 5 has a between the surface rows. During construction, the blast
total length of 23.8 km and has 18 stations. The entire line hole depth of each layer was 1.4 m∼2.8 m, and the blast hole
was laid by underground lines with an average distance of diameter was 42 mm. The blast hole spacing was
1.38 km. Labor Park Station is located at the intersection of 0.6 m∼1.5 m, and the blast hole burden was 0.5 m∼1.3 m. The
Jiefang Road and Ziwei Street in Zhongshan District. The rock drill was used for drilling, and the explosive was 2# rock
main body was arranged along the north and south of Jiefang emulsion explosive, with a specification of
Road. The starting mileage of the station was ZK7 + 532.919 V32 mm ∗ 200 cm, 0.2 kg or 0.3 kg each, and the maximum
(YK7 + 531.658), the ending mileage was K7 + 722.859, the charge per delay was 1.4 kg. During the blasting process, it is
length was 191.2 m, and the width of standard section was necessary to control the maximum charge per delay, reduce
23.3 m. The main body of the station was an underground the number of holes in the same section, and determine
two-story island station with a platform width of 14 m, and reasonable blasting parameters. In addition, measures such
the station roof was covered with soil about 18.7–24.0 m. The as excavation of damping trenches, covering, and protection
6 Shock and Vibration

Labor Park
Ganjingzi Subway Station
District

Shahekou
District Zhongshan
Xigang District
District

Dalian City

Other
School
projects

Labor Park
Subway Station Jiefang Road

Residential
area
Stadium
Arrange the blasting area
Commercial

Figure 2: Sketch map of project overview.

were taken in special areas to reduce the impact of blasting This paper takes the distance from the blast-face, the
on surrounding buildings. maximum charge per delay, total charge, hole depth,
spacing, burden, stemming length, and powder factor as
input variables. The PPV of the measuring point is used as
3.2. Monitoring and Samples Construction. Combined with the output variable to build the samples. The first 72 groups
the actual situation on site, this monitoring adopted the TC- are training samples, and the last 16 groups are test samples.
4850N intelligent blasting vibration monitor, which is a
portable remote sensing vibration monitor with four parallel 4. Results
channels specially designed for engineering blasting. It has
the advantages of lightweight, waterproof, dustproof, high 4.1. Prediction Result Analysis of PPV Using Conventional
pressure resistance, high precision, and wide application. Predictors. According to the previous discussion, we can use
Figure 3 shows the working diagram of the TC-4850N the empirical formulas provided in Table 2 to predict the
wireless blasting monitor. blasting PPV of the monitoring points around the rock
Several blasting vibration test monitoring points were foundation pit. These empirical formulas are mainly de-
arranged near the blasting area of the foundation pit of termined by the distance from the blast-face and the
Suoyuwan South Station. Then, begin to monitor and record maximum charge per delay. They belong to the nonlinear
the blasting vibration velocity at different locations and blasting PPV prediction model of rock foundation pits. In
different blasting parameters during the initial stage of this study, 72 sets of blasting training samples were used for
excavation, as shown in Figure 4. regression analysis to determine the rock properties and
By on-site monitoring, the monitoring data of blasting terrain conditions in the commonly used blasting PPV
vibration at different positions and different blasting pa- prediction formula, as shown in Figure 5. The site constants
rameters during the excavation of the foundation pit were in the various blasting PPV prediction formulas are shown
obtained. There are 88 sets of test data, as shown in Table 3. in Table 4, where PPV is the peak particle velocity (cm/s) at
Shock and Vibration 7

Th
and e
tra sens
ig nal nt ns or
s e
the rum Explosion source m re
its ce
e ives e inst it t ive
rec to th o t s th
or he e
s ens its it ins sign
e
Th ransm tru al
m an
t en d
t
Sensor Sensor Sensor

TC-4850N blasting vibration TC-4850N blasting vibration TC-4850N blasting vibration

measuring instrument measuring instrument measuring instrument
Figure 3: Schematic diagram of TC-4850N wireless blasting monitor.

Signal antenna WIFI antenna

TC-4850N wireless Three vector

blasting vibrometer speed sensor

(a) (b)

(c) (d)

Figure 4: Blasting monitoring chart. (a) TC-4850N blasting monitor. (b) Vibration speed monitoring near high-rise residential buildings.
(c) Vibration speed monitoring of adjacent engineering slopes. (d) Vibration speed monitoring of ground around the foundation pit.

the monitoring point, and Qmax is the maximum charge per different empirical formulas on the equality line with a slope
delay (kg), R is the distance from the blast-face (m), and K, α, of 1 is shown in Figure 6.
β, and n are site constants related to rock properties and By comparing the coefficient of determination (R2), it is
geological conditions. found that the R2 of the above empirical formulas varies
After determining the site constants of the conventional between 0.219 and 0.617. Among them, the empirical for-
empirical formulas, we use 16 sets of blasting test samples to mula of USBM has higher prediction accuracy for PPV, with
measure the obtained empirical formula. The relationship the largest coefficient of determination, and the predicted
between the predicted PPV and the measured PPV of PPV is the closest to the actual measured PPV. The second is
8 Shock and Vibration

Table 3: Model statistical parameters of training samples and test samples.

No. Parameter Unit Symbol Category Min Max Mean SD Median
Training data (72 data)
1 Distance from the blast-face m D Input 6.5 51.9 27.96 12.68 28.45
2 Maximum charge per delay kg Qm Input 0.3 1.35 0.73 0.29 0.68
3 Total charge kg Qt Input 11.4 77.2 42.7 18.3 42.4
4 Hole depth m H Input 1.4 2.8 2.1 0.4 2.1
5 Spacing m S Input 0.6 1.5 1.1 0.2 1.1
6 Burden m B Input 0.5 1.3 0.9 0.2 1.0
7 Stemming length m SL Input 0.8 1.7 1.2 0.3 1.2
8 Powder factor kg/m3 PF Input 0.31 0.82 0.58 0.14 0.59
9 Peak particle velocity cm/s PPV Output 0.17 3.47 1.28 0.90 0.96
Test data (16 data)
1 Distance m D Input 7.3 47.1 27.2 11.3 26.6
2 Maximum charge per delay kg Qm Input 0.30 1.20 0.77 0.27 0.80
3 Total charge kg Qt Input 12.1 73.8 41.8 18.0 35.9
4 Hole depth m H Input 1.5 2.6 2.0 0.4 2.1
5 Spacing m S Input 0.8 1.4 1.1 0.2 1.1
6 Burden m B Input 0.6 1.3 1.0 0.2 1.0
7 Stemming length m SL Input 0.8 1.6 1.1 0.3 1.2
8 Powder factor kg/m3 PF Input 0.33 0.79 0.56 0.14 0.54
9 Peak particle velocity cm/s PPV Output 0.73 2.08 1.37 0.42 1.30

3.5 3.5

3.0 3.0

2.5 2.5
PPV (cm/s)

PPV (cm/s)

2.0 2.0

1.5 PPV = 7.707(R/√Qmax )–0.541 1.5

R2 = 0.319 PPV = 5.471(√Qmax/R2/3)1.179
1.0 1.0
R2 = 0.261
0.5 0.5

0.0 0.0
0 10 20 30 40 50 60 70 80 90 0.1 0.2 0.3 0.4 0.5 0.6
USBM scaled distace, R/(Qmax )1/2 (m/kg1/2) LK scaled distace, (Qmax /R2/3)1/2 ((kg/m2/3)1/2)

Measure PPV Measure PPV

Fitted curve Fitted curve
(a) (b)

3.5 3.5

3.0 3.0

2.5 2.5
PPV (cm/s)

PPV (cm/s)

2.0 2.0 PPV = 7.657(R/3√Qmax)–0.549

1.5 PPV = 4.185R–0.302Qmax
0.643 1.5 R2 = 0.302

1.0 R2 = 0.213
1.0
0.5 0.5
0.0 0.0
0.0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 60 70
General predictor scaled AH scaled distace, R/(Qmax)1/3 (m/kg1/3)
distace, R–0.302Qmax0.643 (m–0.302kg0.643)
Measure PPV
Measure PPV Fitted curve
Fitted curve
(c) (d)
Figure 5: Continued.
Shock and Vibration 9

3.5 3.5

3.0 3.0

2.5 2.5
PPV = 0.649+15.581(R/√Qmax)–1
PPV (cm/s)

PPV (cm/s)
2.0 2.0
R2 = 0.257
1.5 1.5
PPV = 5.471(Qmax/3√R2)0.589
1.0 1.0
R2 = 0.260
0.5 0.5

0.0 0.0
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0 10 20 30 40 50 60 70 80 90
BIS scaled distace, Qmax/(R)2/3 (kg/m2/3) CMRI scaled distace, R/(Qmax)1/2 (m/kg1/2)

Measure PPV Measure PPV

Fitted curve Fitted curve
(e) (f )

Figure 5: Relationship between PPV and scaled distance of the conventional empirical formulas. (a) Relationship between PPV and scaled
distance of the USBM empirical formula. (b) Relationship between PPV and scaled distance of the LK empirical formula. (c) Relationship
between PPV and scaled distance of the general predict empirical formula. (d) Relationship between PPV and scaled distance of the AH
empirical formula. (e) Relationship between PPV and scaled distance of the BIS empirical formula. (f ) Relationship between PPV and scaled
distance of the CMRI empirical formula.

Table 4: The conventional empirical formulas and their site constants.

Site constants
Name Equation
K α β n R2
􏽰���� − α
USBM PPV � K(R/ Qmax )
􏽰�������� 7.707 0.541 0.319
LK PPV � K( Qmax /R2/3 )α 5.471 1.179 0.261
− α β
General predictor PPV � KR􏽰 Qmax
���� 4.185 0.302 0.643 0.213
− α
AH PPV � K(R/ 3 Q√ ���) α
max 7.657 0.549 0.302
3
BIS R2 )
PPV � K(Qmax􏽰/ ���� 5.471 0.589 0.26
CMRI PPV � n + K(R/ Qmax )− 1 15.581 0.649 0.257

2.2 2.2

2.0 2.0
Predicted PPV by USBM (cm/s)

Equality line Equality line

Predicted PPV by LK (cm/s)

1.8 1.8

1.6 1.6

1.4 1.4

1.2 1.2

1.0 y = 0.704x + 0.335 1.0 y = 0.491x + 0.652

R2 = 0.617 R2 = 0.372
0.8 0.8

0.6 0.6
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2
Measured PPV (cm/s) Measured PPV (cm/s)
(a) (b)
Figure 6: Continued.
10 Shock and Vibration

Predicted PPV by general predictor (cm/s) 2.2 2.2

2.0 2.0

Predicted PPV by AH (cm/s)

Equality line Equality line
1.8 1.8

1.6 1.6

1.4 1.4

1.2 1.2

1.0 y = 0.362x + 0.832 1.0 y = 0.717x + 0.314

R2 = 0.219 R2 = 0.614
0.8 0.8

0.6 0.6
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2
Measured PPV (cm/s) Measured PPV (cm/s)
(c) (d)

2.2 2.2

2.0 2.0

Predicted PPV by CMRI (cm/s)

Equality line Equality line
Predicted PPV by BIS (cm/s)

1.8 1.8

1.6 1.6

1.4 1.4

1.2 1.2

1.0 1.0
y = 0.491x + 0.653 y = 0.630x + 0.403
0.8 R2 = 0.372 0.8 R2 = 0.565

0.6 0.6
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2
Measured PPV (cm/s) Measured PPV (cm/s)
(e) (f )

Figure 6: Relationship between measured and predicted PPV by USBM (a), LK (b), general predict (c), AH (d), BIS (e), and CMRI (f ). (a)
Relationship between measured and predicted PPV by USBM. (b) Relationship between measured and predicted PPV by LK. (c) Rela-
tionship between measured and predicted PPV by general predict. (d) Relationship between measured and predicted PPV by AH. (e)
Relationship between measured and predicted PPV by BIS. (f ) Relationship between measured and predicted PPV by CMRI.

the empirical formulas of AH and CMRI, and the prediction are quickly searched. Among them, set the population size
accuracy is 0.614 and 0.565, respectively. However, the NP � 60 and the variation factor F, crossover factor CR, and
general predict empirical formula has the lowest prediction maximum number of iterations are 0.7, 0.8, and 150, re-
accuracy, only 0.219, and the coefficient of determination is spectively. According to the application experience of ma-
the smallest. chine learning models in the field of geotechnical
engineering, the convergence criterion is set to reach the
maximum number of iterations or the difference of the
4.2. Prediction Result Analysis of PPV Using DE-GP. As target value between iteration steps is less than 1.0 × 10− 6
mentioned earlier, about 80% of the blasting data samples [35–37]. Through the learning and prediction of the DE-GP
(72 sets of data) are regarded as training samples, which are model, the relationship between the actual measured PPV
used to establish a nonlinear mapping relationship for and the predicted PPV in the training phase and the test
predicting PPV; the remaining 20% of the blasting data phase is shown in Figure 7.
samples (16 sets of data) are regarded as test samples. It can be seen from Figure 7 that the DE-GP model is
During the blasting construction of rock foundation pits, used to predict the training samples and test samples, and R2
the DE-GP model is used to predict PPV. By learning the is 0.937 and 0.955, respectively. It can be seen from the
training samples, taking the minimum fitness value of the measured PPV and the predicted PPV that the prediction
population as the goal, the DE algorithm is used for global effect is good, and the predicted PPV is in good agreement
optimization, and the optimal hyperparameters in the GP with the actual PPV.
Shock and Vibration 11

4.0 4.0

3.5 3.5
Equality line
3.0 3.0
Predicted PPV (cm/s)

2.5 2.5

PPV (cm/s)
2.0 2.0

1.5 1.5
y = 0.832x + 0.226
1.0 R2 = 0.937 1.0

0.5 0.5

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 1 8 16 24 32 40 48 56 64 72
Measured PPV (cm/s) Dataset Number

Measured
DE-GP
(a)

2.2 2.2
Predicted PPV by DE-GP (cm/s)

2.0 2.0
Equality line
1.8 1.8

1.6 1.6
PPV (cm/s)

1.4 1.4
1.2
1.2
y = 0.981x – 0.005 1.0
1.0
R2 = 0.955 0.8
0.8
0.6
0.6
0.4
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
Measured PPV (cm/s) Dataset Number

Measured
DE-GP
(b)

Figure 7: Analysis and comparison of the measured PPV and the predicted PPV by the DE-GP model in the training phase and the test
phase. (a) Training phase. (b) Test phase.

5. Discussion seen from Figure 8(b) that when l � 0.3 and σ f � − 1.9, the
RMSE of the test samples reaches the minimum. After that, as
5.1. Influence of GP Parameters on Regression Accuracy. l and σ f change, the RMSE of the test samples becomes larger
The 72 groups of blasting monitoring data mentioned above and larger. It can be seen from Figure 8(c) that when σ n is
are used as training samples for GP model learning. Among greater than 0, the RMSE of the test samples reaches above 0.5;
them, the noise variance σ n, the signal variance σ f of the as σ n approaches − 1.8 and l approaches − 0.5, the RMSE of the
kernel function, and the correlation determination hyper- test samples gradually decreases and finally reaches the
parameter l are important hyperparameters in the GP model. minimum. The value is 0.06. Changes in the values of σ n, σ f,
Under the influence of different hyperparameters, the root and l will have a greater impact on the root mean square error
mean square error analysis of 16 sets of test samples is of the test samples. Therefore, it is necessary to select ap-
performed. It can obtain the effects of σ n and σ f, σ f and l, and propriate parameters to ensure the regression accuracy of the
σ n and l on the prediction accuracy of the GP model, as GP model. Because the GP algorithm cannot search for the
shown in Figure 8. optimal hyperparameters quickly and accurately, the DE
Figure 8(a) shows that when σ f approaches − 3 and σ n algorithm is introduced to optimize the hyperparameters,
approaches 0, the greater the RMSE of the GP model test thereby avoiding the blindness of the hyperparameters se-
samples; and when σ f � − 1.6 and σ n � − 2, the RMSE of the test lection in the GP and achieving the purpose of rapid selection
samples takes the minimum value, which is 0.114. It can be and reasonable value selection.
12 Shock and Vibration

Value of RMSE Value of RMSE

0.610 0.411

0.560 0.376

0.511 0.340
0.64 0.50
0.461 0.45 0.305

rror
re error

0.56 0.40

square e
0.48 0.412 0.35 0.270
0.30
a

0.40 0.362 0.235

ean squ

0.25
0.32

an
0.20 0.199
0.312

root me
0.24 0.15
root m

0.10 -2.5
0.16 0.0 0.263 0.164
0.05 -2.0
0.08 -0.6

n he
-1.0 -1.5 0.129

tio f t
-1.2

sn
-3.0 0.213

nc e o
-0.5

ce

sf
-2.5 -1.8 Cor -1.0
ian

fu c
Sig rela 0.0 0.5

el ian
nal -2.0-1.5 -2.4 r 0.164 -0.5 0.093
hyp tion m
va

rn ar
ker varia -1.0 -3.0 1.0

ke al v
se

erp e 0.0
nel n
fun ce of -0.5 0.0 -3.6 ara asure 1.5 2.0 0.058
oi

gn
0.114
N

me m

Si
ctio the ters ent
ns l
f

(a) (b)
Value of RMSE
0.616
0.560
0.505
0.7
0.449
rror

0.6
0.394
e

0.5
square

0.4 0.338
0.3
ean

0.282
0.2
root m

-1.2 0.227
t

0.1
s l en

-0.6
et rem

0.80.0 0.0 0.171

am su

-0.8 0.6
er
ar ea

No -1.6 0.116
rp m

ise
var -2.4-3.2 1.2
pe on

ian 1.8
ce s -4.0-4.8
hy lati

0.060
e

n
rr
Co

(c)

Figure 8: The influence of GP hyperparameters on regression accuracy. (a) The influence of σ n and σ f on regression accuracy. (b) The
influence of σ f and l on regression accuracy. (c) The influence of σ n and l on regression accuracy.

5.2. Influence of DE Parameters on Optimization Effect. development ability of the DE and improve the convergence
When using the DE-GP model, there are many factors that speed of the DE, and also it is easy to make the population
affect the optimization effect of the DE, and they are directly converge prematurely. The cross factor CR affects the weight
related to the selection of hyperparameters in the GP model. of evolutionary information. Increasing CR can improve the
The control parameters in the DE mainly include mu- diversity of the population and speed up the convergence
tation factor F, crossover factor CR, and population size NP. speed of the DE; reducing CR is beneficial to the analysis of
Usually, NP is 5∼10D (D stands for dimensionality), F is individual separable problems. Select different variation
0.5∼0.9, and CR is 0.5∼0.9. It can not only ensure a higher factor F and cross factor CR to perform optimization search
success rate of optimization but also increase the speed of trial calculation, and the corresponding iterative conver-
convergence. The mutation factor F and the cross factor CR gence curve is shown in Figure 9.
have a greater impact on the speed and stability of the search. It can be seen from Figure 9 that when CR � 0.7 and F ∈
Improper selection of parameters may lead to evolutionary [0.5, 0.9], DE iterative search can converge, but the iteration
stagnation and premature maturity. The mutation factor F speed is different. When F � 0.8, the convergence speed is the
mainly affects the search steps. Increasing F can increase the fastest, the number of iteration steps is less, and the iteration
search range of the DE and increase the diversity of the steps reach 54 steps to complete convergence. When F � 0.6,
population but, at the same time, weaken the convergence the convergence speed is the slowest, and the iterative steps
speed of the algorithm; reducing F can increase the reach 85 steps before it converges completely. Therefore,
Shock and Vibration 13

0.25 0.25

0.20 0.20
Convergence value

Convergence value
0.15 0.15

0.10 0.10

0.05 0.05

0.00 0.00
0 40 80 120 160 200 0 40 80 120 160 200
Number of iterations Number of iterations

F=0.5 F=0.8 CR=0.5 CR=0.8

F=0.6 F=0.9 CR=0.6 CR=0.9
F=0.7 CR=0.7
(a) (b)

Figure 9: Iterative curves of different variation factors and crossover factors. (a) When CR � 0.7, iterative curves of different variation factors
F. (b) When F � 0.8, iterative curves of different crossover factors CR.

choosing appropriate search parameters can greatly save 0.25

time and improve the search speed. When the variation
factor F � 0.8 is fixed and CR ∈ [0.5, 0.9], the convergence
curve is relatively stable. When CR � 0.9, the search result 0.20
reaches the fastest convergence, and it is completely con-
Convergence value

verged in 39 iterations. When CR � 0.6, the convergence 0.15

speed is the slowest, and the iterative steps reach 80 steps
before it starts to converge. CR and F should be dynamically
adjusted according to actual problems, and further com- 0.10
parisons should be made within the selection range of
possible convergence parameters in order to obtain the 0.05
optimal convergence result.
When the variation factor F � 0.8 and the crossover
factor CR � 0.9 are fixed, increasing NP can improve the 0.00
diversity of the population but, at the same time, reduce the 0 40 80 120 160 200
convergence speed of the population. Reducing NP can Number of iterations
increase the convergence speed, but it is easy to lead to
NP=40 NP=100
premature convergence. When the population number NP is
NP=60 NP=120
equal to 40, 60, 80, 100, and 120, respectively, the iterative
NP=80
convergence curves are shown in Figure 10.
It can be seen from Figure 10 that when NP is 40 and 60, Figure 10: Iterative convergence curves of different NP.
the optimal value is searched for 54 iterations and 48 it-
erations, respectively. But when NP is 80, 100, and 120, it
takes 57, 73, and 85 iterations to search for the optimal value. the search to stagnate. It can also be seen from Figure 10 that
NP mainly reflects the amount of population information in the fitness convergence speed does not decrease mono-
DE. The larger the population size, the richer the population tonically with the increase of NP. When NP is set to 60,
information contained, it is more likely to find the global convergence can be achieved at the fastest speed. So, NP is set
solution, and it is less likely for individual optimal solutions to 60.
to dominate the overall solution. However, the calculation Choosing appropriate initial parameters can not only
time is relatively long, and the calculation convergence is save time and improve convergence speed but also avoid the
relatively slow, which lead to an increase in the calculation algorithm from falling into local optimality. If the selection
cost of a single population. On the contrary, the diversity of of parameters is not appropriate, it may cause the algorithm
the population is restricted, which is not conducive to DE search to stagnate due to excessive emphasis on the search
obtaining the global optimal solution and may even cause ability or cause the algorithm to converge prematurely due
14 Shock and Vibration

to excessive emphasis on the development ability. In this PPV is also closest to the measured PPV. Plot the measured
research, the final DE algorithm control parameters settings: and predicted PPV of 16 sets of test samples obtained by
NP � 60, F � 0.8, and CR � 0.9. various prediction methods as shown in Figure 12. The
specific index evaluation effects of the forecast results are
shown in Table 7.
5.3. Computing Ability Evaluation of Different Prediction Generally, we take δ � (y′i − yi /yi ) × 100% as the judg-
Methods. Since the constructed blasting PPV model will ment of the relative error of the test samples, which can
directly affect the accuracy of the prediction results, it is better reflect the credibility of the prediction. The relative
necessary to establish the evaluation index of the model to error and calculation time of the test samples using the
test the prediction effect. From a statistical point of view, it is LSSVM, ANN, GP, and DE-GP models and USBM empirical
more one sided to comprehensively evaluate a predictive formula are shown in Figure 13.
model with only one performance index. Therefore, in this It can be seen from Figure 13 that the relative error range
research, statistical indicators such as coefficient of deter- of the GP model prediction results is − 15.81%∼24.29%, and
mination (R2), root mean square error (RMSE), mean square the average relative error is 2.62%, which is significantly
error (MSE), mean absolute error (MAE), mean absolute smaller than the prediction results of the LSSVM and ANN
percentage error (MAPE), variance account for (VAF), and models and USBM empirical formula. The regression ac-
standard deviation (SD) are used to evaluate the perfor- curacy of the GP model optimized by the DE is significantly
mance of the model, as shown in Table 5. improved, and the relative error range of the prediction
Among them, yi is the measured PPV, yi is the predicted results is only − 17.78%∼11.05%. It shows that the DE-GP
PPV, y is the average value of the measured PPV, and m is hybrid prediction model established in this research can be
the total number of samples. It should be pointed out that effectively applied to the prediction of PPV in the con-
the closer R2 is to 1, the closer the values of RMSE, MSE, struction field of rock foundation pit.
MAE, MAPE, and SD are to 0, and the closer the VAF is to In addition, the computation cost of different models in
100, indicating that the performance of the prediction this research refers to the time spent under the conditions of
models is better. six-thread CPU processor, 8 GB running memory, and
According to the site constants in the traditional em- RTX2060 graphics card. It can be seen from Figure 13 that
pirical formulas obtained above, 16 sets of test samples are the computation time cost of different models is different.
used to test the predictive ability of the above several Compared with the GP model, the computation time cost of
commonly used empirical formulas. Each index is shown in the DE-GP model has increased by 1′34″. However, from the
Table 6. perspective of the improvement of calculation accuracy, this
It can be seen from Table 6 that when using the USBM part of the increase in computation time cost is completely
empirical formula, R2 � 0.62, RMSE � 0.26, MSE � 0.07, acceptable.
MAE � 0.21, MAPE � 15.27, VAF � 60.49, and SD � 0.29,
and all evaluation indicators are better than other empirical
formulas. Therefore, the empirical formula of USBM is 5.4. Sensitivity Analysis. In order to use the PPV to measure
selected among the six empirical formulas for comparison BIGV, the most important step is to determine the sensitivity
with other blasting intelligent algorithms. of the eight blasting influence factors. The cosine amplitude
In this research, the performance of the hybrid DE-GP method (CAM) is a commonly used sensitivity test method
prediction model is compared with LSSVM, ANN, and GP [38]. Among them, the influence of each blasting influence
models and USBM empirical formula. In the ANN model, factor on PPV can be calculated by the correlation coefficient
the number of hidden layer nodes Nf of the algorithm is 5, rij , that is, the degree of influence of each blasting influence
the number of iterations is Ep � 100, and the learning rate factor on PPV can be obtained from the sensitivity analysis:
is Lr � 0.1. In the GP model, the variance of noise is
􏽐m n�1 Xin Xjn
σ 2n � 3.72, the signal variance of the kernel function is rij � 􏽱���������������, (8)
σ f � − 1.7, and the hyperparameter for correlation mea- 􏽐n�1 X2in 􏽐m
m 2
n�1 Xjn
surement is l � 0.3. In the core parameters of the LSSVM
model, the regular parameter is c � 9.54, and square where rij is the degree of each blasting influence factor; Xin is
bandwidth is σ 2 � 2.17. In the DE-GP model, the pop- the i-th blasting influence factor of the n-th sample data; Xjn
ulation number NP is 60, mutation factor F is 0.8, is the blasting PPV of the n-th sample data; and m is the
crossover factor CR is 0.9, and the upper limit of evo- number of samples.
lutionary algebra is 100. The relationship between mea- The value of the correlation coefficient rij is between 0
sured PPV and predicted PPV using LSSVM, ANN, GP, and 1. The closer the correlation coefficient is to 1, the
and DE-GP models is shown in Figure 11. greater the influence of input parameters on PPV. Similarly,
It can be seen from Figure 11 that using the ANN model the smaller the influence of input parameters on output
to predict the PPV of rock foundation pits has the lowest parameter, the closer the correlation coefficient is to 0.
accuracy, R2 is 0.724. While using the LSSVM model and the The sensitivity of blasting input parameters are analyzed
GP model, R2 is 0.832 and 0.885, respectively. After finally based on the blasting sample set data and the established
using the GP model optimized by the DE to predict, the LSSVM, ANN, GP, and DE-GP models, as shown in
accuracy is the highest, reaching 0.955, and the predicted Figure 14.
Shock and Vibration 15

Table 5: Statistical indicators of blasting PPV prediction effect.

Statistical indicators Formulas Desired value
Coefficient of determination R2 � 1 − 􏽐m (y ′2 m
i − yi) /􏽐i�1 (yi − yi )
􏽱�����������������
i�1
′2 1
m 2
Root mean square error RMSE � (1/m) 􏽐i�1 (yi − y′i) 0
Mean square error MSE � (1/m)􏽐m i�1 (yi − yi)
′2 0
Mean absolute error MAE � (1/m) 􏽐m i�1 |y i − y ′i| 0
Mean absolute percentage error MAPE � (1/m) 􏽐m i�1 |y i − y′i i × 100%
/y | 0
Variance account for VAF � [1 − var(yi − y′i)/var(yi )] × 100% 100
􏽱��������������������������������
Standard deviation SD � (1/m)􏽐m i�1 [yi − (1/m)(y1 + · · · + ym )]
2
0

Table 6: Predictive effects of several empirical formulas.

2
Name R RMSE MSE MAE MAPE VAF SD Rank
USBM 0.62 0.26 0.07 0.21 15.27 60.49 0.53 1
LK 0.39 0.35 0.12 0.26 22.15 33.32 0.69 4
General predict 0.22 0.41 0.16 0.3 26.24 12.61 0.83 6
AH 0.61 0.28 0.08 0.21 15.85 59.72 0.57 2
BIS 0.37 0.35 0.12 0.27 24.62 27.51 0.75 5
CMRI 0.57 0.3 0.09 0.23 17.18 55.79 0.64 3

2.2 2.2

2.0
Predicted PPV by LSSVM (cm/s)

y = 0.984x + 0.063 2.0

Predicted PPV by ANN (cm/s)

Equality line
1.8 R2 = 0.832 1.8

1.6 1.6

1.4 1.4

1.2 1.2

1.0 Equality line 1.0

y = 0.771x + 0.279
0.8 0.8 R2 = 0.724

0.6 0.6
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2
Measured PPV (cm/s) Measured PPV (cm/s)
(a) (b)

2.2 2.2

2.0 2.0
Predicted PPV by DE-GP (cm/s)

y = 1.037x – 0.016 Equality line

Predicted PPV by GP (cm/s)

1.8 R2 = 0.885 1.8

1.6 1.6

1.4 1.4

1.2 1.2

1.0 1.0 y = 0.981x – 0.005

Equality line R2 = 0.955
0.8 0.8

0.6 0.6
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2
Measured PPV (cm/s) Measured PPV (cm/s)
(c) (d)

Figure 11: Relationship between measured and predicted PPV by LSSVM (a), ANN (b), GP (c), and DE-GP (d). (a) Relationship between
measured and predicted PPV by LSSVM. (b) Relationship between measured and predicted PPV by ANN. (c) Relationship between
measured and predicted PPV by GP. (d) Relationship between measured and predicted PPV by DE-GP.
16 Shock and Vibration

2.4

2.1

1.8
PPV (cm/s) 1.5

1.2

0.9

0.6
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
Test sample number

Actural PPV GP
ANN LSSVM
USBM DE-GP
Figure 12: Relationship between measured and predicted PPV by various prediction methods.

Table 7: Evaluation effect of various prediction methods.

2
Name R RMSE MSE MAE MAPE VAF SD Rank
USBM 0.62 0.26 0.07 0.21 15.27 60.49 0.53 5
LSSVM 0.83 0.19 0.04 0.16 11.37 80.52 0.46 3
ANN 0.72 0.23 0.05 0.17 13.45 72.19 0.48 4
GP 0.89 0.16 0.02 0.11 8.07 86.93 0.46 2
DE-GP 0.96 0.1 0.01 0.08 6.45 95.44 0.42 1

60

45
Relative error/(%)

30

15

-15

-30
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
Test sample number

LSSVM (1′18″) DE–GP (2′37″)

ANN (1′25″) USBM
GP (1′03″)
Figure 13: The relative error and calculation time of different PPV prediction methods.

It can be seen from Figure 14 that the blasting in- The less influential blasting factors are stemming length
fluence factors of the LSSVM, ANN, GP, and DE-GP and powder factor, and the correlation coefficient rij is
models have different effects on the blasting PPV. between 0.62 and 0.78. Therefore, the sensitivity analysis
However, the three most influential factors are distance result can be used to optimize the blasting design pa-
from the blast-face, maximum charge per delay, and total rameters, which can be used to guide project
charge, all of which have a degree of influence above 0.85. construction.
Shock and Vibration 17

distance from the blast-face Data Availability

1.0
The data used to support the findings of this study are
0.9 available from the corresponding author upon request.
stemming length total charge
0.8
0.7 Conflicts of Interest
0.6 The authors declare that there are no conflicts of interest
burden hole depth regarding the publication of this paper.

Acknowledgments
This work was supported by the National Natural Science
spacing powder factor Foundation of China (Grant nos. 52078093 and 51678101),
Central University Basic Research Fund Special Funds
(Grant no. 3132014326), and Liaoning Revitalization Talents
maximum charge per delay Program (Grant no. XLYC1905015).
DE-GP LSSVM
GP ANN References
Figure 14: PPV sensitivity analysis of different prediction methods. [1] D. J. Armaghani, E. Momeni, and S. V. A. N. K. Abad,
“Feasibility of ANFIS model for prediction of ground vi-
brations resulting from quarry blasting,” Environmental Earth
6. Conclusions Sciences, vol. 74, no. 4, pp. 2845–2860, 2015.
[2] L. Ma, K. Li, and X. Ding, “Development and application of
The present research established a PPV prediction model blast casting technique in large-scale surface mines: a case
based on the DE-GP intelligent algorithm. It was compared study of Heidaigou surface coal mine in China,” Shock and
with empirical formulas and other intelligent algorithms, Vibration, vol. 2016, Article ID 8496742, 11 pages, 2016.
taking the blasting of the foundation pit at the Suoyuwan [3] S. Lu, C. Zhou, and Z. Zhang, “PPV criterion of a rock slope
South Station as an example. The following conclusions are imbedded with a fault subjected to blasting P-waves,” Shock
and Vibration, vol. 2020, Article ID 8865981, 7 pages, 2020.
obtained:
[4] Y. Gou, X. Shi, and X. Huo, “Motion parameter estimation
(1) The conventional empirical formulas are mainly and measured data correction derived from blast-induced
related to distance from the blast-face and the vibration: new insights,” Measurement, vol. 135, pp. 213–230,
maximum charge per delay. However, the calculated 2019.
PPV deviates greatly from the measured PPV. [5] M. Hasanipanah, A. Shahnazar, and H. B. Amnieh, “Pre-
diction of air-overpressure caused by mine blasting using a
(2) The blasting PPV prediction model based on the DE- new hybrid PSO-SVR model,” Engineering with Computers,
GP method was established for the first time, and it vol. 33, no. 1, pp. 23–31, 2017.
was found to have high accuracy by comparing the [6] D. Ainalis, L. Ducarne, and O. Kaufmann, “Improved analysis
measured PPV and the predicted PPV. of ground vibrations produced by man-made sources,” The
Science of the Total Environment, vol. 616, pp. 517–530, 2018.
(3) The DE-GP blasting prediction model used for the
[7] X. Wu, Q. Guo, and Y. Zhang, “A new mathematical model
rock blasting prediction. Compared with the for predicting the surface vibration velocity on the step to-
LSSVM, ANN, and GP prediction models and em- pography,” Shock and Vibration, vol. 2018, Article ID
pirical formulas, the prediction accuracy was im- 8356254, 8 pages, 2018.
proved, and the relative error was significantly [8] H. Agrawal and A. K. Mishra, “Probabilistic analysis on
reduced. scattering effect of initiation systems and concept of modified
(4) The CAM was used to analyze the sensitivity of charge per delay for prediction of blast induced ground vi-
blasting influence factors for different prediction brations,” Measurement, vol. 130, pp. 306–317, 2018.
[9] N. Fouladgar, M. Hasanipanah, and H. B. Amnieh, “Appli-
models. The impact of the distance from the blast-
cation of cuckoo search algorithm to estimate peak particle
face, maximum charge per delay, and total charge are velocity in mine blasting,” Engineering with Computers,
the largest; stemming length and powder factor are vol. 33, no. 2, pp. 181–189, 2017.
relatively small. [10] D. Ainalis, O. Kaufmann, and J.-P. Tshibangu, “Modelling the
The PPV prediction model based on DE-GP can be source of blasting for the numerical simulation of blast-in-
duced ground vibrations: a review,” Rock Mechanics and Rock
effectively applied to the blasting prediction in geotech-
Engineering, vol. 50, no. 1, pp. 171–193, 2017.
nical engineering. It can control the adverse effects of [11] X. Qiu, X. Shi, and Y. Gou, “Short-delay blasting with single
blasting and reduce the harm of BIGV to the surrounding free surface: results of experimental tests,” Tunnelling and
environment. This method can provide strong technical Underground Space Technology, vol. 74, pp. 119–130, 2018.
support for the blasting construction of similar [12] M. Khandelwal and M. Saadat, “A dimensional analysis ap-
engineering. proach to study blast-induced ground vibration,” Rock
18 Shock and Vibration

Mechanics and Rock Engineering, vol. 48, no. 2, pp. 727–735, the aim of reducing ground vibration using robust tech-
2015. niques,” Geomechanics and Geophysics for Geo-Energy and
[13] X. Shi, X. Qiu, and J. Zhou, “A comparative study of ground Geo-Resources, vol. 6, no. 2, p. 40, 2020.
and underground vibrations induced by bench blasting,” [28] L.-D. Thang, S. Quoc-Hung Nguyen, and H. Nguyen-Xuan,
Shock and Vibration, vol. 2016, Article ID 2143575, 9 pages, “Balancing composite motion optimization,” Information
2016. Sciences, vol. 520, pp. 250–270, 2020.
[14] N. Hoang, B. Xuan-Nam, and T. Quang-Hieu, “Evaluating [29] W. I. Duvall and B. Petkof, Spherical Propagation of Explosion-
and predicting blast-induced ground vibration in open-cast Generated Strain Pulses in Rock, U.S. Department of the
mine using ANN: a case study in Vietnam,” Sn Applied Interior, Bureau of Mines, Washington, DC, USA, 1958.
Sciences, vol. 1, no. 1, 2019. [30] B. K. Langefors U and K. S. Kenky� ujo, The Modern Technique
[15] E. Ghasemi, M. Ataei, and H. Hashemolhosseini, “Develop- of Rock Blasting, Wiley, New York, NY, USA, 1963.
ment of a fuzzy model for predicting ground vibration caused [31] F. I. W. Davies and P. B. Attewell, “Ground vibrations from
by rock blasting in surface mining,” Journal of Vibration and shallow sub-surface blasts,” Engineering, vol. 217, pp. 553–
Control, vol. 19, no. 5, pp. 755–770, 2013. 559, 1964.
[16] M. Hajihassani, D. J. Armaghani, and A. MartoEdy, “Ground [32] N. Ambraseys, Rock Mechanics In Engineering Practice, Wiley,
vibration prediction in quarry blasting through an artificial New York, NY, USA, 1968.
neural network optimized by imperialist competitive algo- [33] BOI Standard, Criteria for Safety and Design of Structures
rithm,” Bulletin of Engineering Geology and the Environment, Subjected to Underground Blast, ISI-6922, New Delhi, India,
vol. 74, no. 3, pp. 873–886, 2015. 1973.
[17] M. Hasanipanah, M. Monjezi, and A. Shahnazar, “Feasibility [34] S. K. Biyosi and B. K. Pal, “Putting ground vibration pre-
of indirect determination of blast induced ground vibration dictions into practice,” International Journal of Rock Me-
based on support vector machine,” Measurement, vol. 75, chanics and Mining Science & Geomechanics Abstracts, vol. 30,
pp. 289–297, 2015. no. 5, pp. 302–310, 1993.
[18] M. Amiri, H. B. Amnieh, and M. Hasanipanah, “A new [35] Y. Zhang, G. Su, and B. Liu, “A novel displacement back
combination of artificial neural network and K-nearest analysis method considering the displacement loss for un-
neighbors models to predict blast-induced ground vibration derground rock mass engineering,” Tunnelling and Under-
and air-overpressure,” Engineering with Computers, vol. 32, ground Space Technology, vol. 95, Article ID 103141.1, 2020.
no. 4, pp. 631–644, 2016. [36] G. Su, J. Jiang, and B. Yu, “A Gaussian process-based response
[19] E. Ghasemi, H. Kalhori, and R. Bagherpour, “A new hybrid surface method for structural reliability analysis,” Structural
ANFIS-PSO model for prediction of peak particle velocity due Engineering & Mechanics, vol. 56, no. 4, pp. 549–567, 2015.
to bench blasting,” Engineering with Computers, vol. 32, no. 4, [37] S. Y. Su Guoshao and L. Yan, “A new method for forecasting
pp. 607–614, 2016. of blasting effect in rock mass,” Chinese Journal of Rock
[20] E. Ebrahimi, M. Monjezi, and M. R. Khalesi, “Prediction and Mechanics and Engineering, vol. 26, no. 1, pp. 3509–3514,
optimization of back-break and rock fragmentation using an 2007.
artificial neural network and a bee colony algorithm,” Bulletin [38] P. F. Asl, M. Monjezi, and J. K. Hamidi, “Optimization of
of Engineering Geology and the Environment, vol. 75, no. 1, flyrock and rock fragmentation in the Tajareh limestone mine
pp. 27–36, 2016. using metaheuristics method of firefly algorithm,” Engi-
[21] K. Taheri, M. Hasanipanah, and S. B. Golzar, “A hybrid ar- neering with Computers, vol. 34, no. 2, pp. 241–251, 2018.
tificial bee colony algorithm-artificial neural network for
forecasting the blast-produced ground vibration,” Engineering
with Computers, vol. 33, no. 3, pp. 689–700, 2017.
[22] M. Hasanipanah, R. Naderi, and J. Kashir, “Prediction of
blast-produced ground vibration using particle swarm opti-
mization,” Engineering with Computers, vol. 33, no. 2,
pp. 173–179, 2017.
[23] R. S. Faradonbeh and M. Monjezi, “Prediction and mini-
mization of blast-induced ground vibration using two robust
meta-heuristic algorithms,” Engineering with Computers,
vol. 33, no. 4, pp. 835–851, 2017.
[24] D. J. Armaghani, M. Hasanipanah, and H. B. Amnieh,
“Feasibility of ICA in approximating ground vibration
resulting from mine blasting,” Neural Computing & Appli-
cations, vol. 29, no. 9, pp. 457–465, 2018.
[25] X. Zhang, H. Nguyen, and X.-N. Bui, “Novel soft computing
model for predicting blast-induced ground vibration in open-
pit mines based on particle swarm optimization and
XGBoost,” Natural Resources Research, vol. 29, no. 2,
pp. 711–721, 2020.
[26] E. Tian, J. Zhang, M. S. Tehrani, and A. Surendar, “Develop-
ment of GA-based models for simulating the ground vibration
in mine blasting,” Engineering with Computers, vol. 35, no. 3,
pp. 849–855, 2019.
[27] A. Rezaeineshat, M. Monjezi, and A. Mehrdanesh, “Opti-
mization of blasting design in open pit limestone mines with