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

Authors: Tengfei Jiang; Annan Jiang; Shuai Zheng; Mengfei Xu
Journal: Shock and Vibration
 (2021)

---

## Abstract

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.


---

## Body

## 1. Introduction

At present, the widespread application of blasting technology has penetrated into many areas of the national economy. Because of its high efficiency, economy, and speed, it has long been favored in the field of engineering construction [1–3]. However, blasting also produces a series of harmful effects, mainly including blast-induced ground vibration (BIGV), blasting flying rocks, noise, shock wave, etc. [4, 5].BIGV adversely affects surrounding rock masses and nearby structures and even causes damage [6]. Most countries in the world use peak particle velocity (PPV) as the measurement index of BIGV [7, 8]. However, the propagation medium of blasting vibration waves is heterogeneous rock mass, and there are many influence factors of PPV [9, 10]. How to accurately predict PPV has become the primary issue in the field of blasting construction safety technology and scientific research [11].In the past blasting engineering, the most commonly used method for PPV prediction is the empirical formulas. The empirical formulas are used to perform regression analysis on the measured data of blasting vibration. It mainly depends on distance from the blast-face and maximum charge per delay [12, 13]. However, the prediction accuracy of the empirical formulas is generally low, and the regression calculation effect for some blasting projects is not ideal. As a result, the blasting vibration attenuation formulas regressed by the empirical formulas is prone to excessive deviation in the guidance of site construction [14, 15].Because the prediction of PPV is a complex nonlinear problem, researchers generally use intelligent algorithm models with powerful processing capability to predict PPV, such as artificial neural network (ANN), genetic algorithm (GA), particle swarm optimization (PSO), and support vector machine (SVM)[1, 16–27], as shown in Table 1.Table 1
Prediction of PPV by intelligent algorithms.
ReferencesModelsInput parametersNumber of samplesCoefficient of determination (R2)Armaghani et al. [1]ANFIS; ANNQm, D1090.973 (ANFIS); 0.949 (ANN)Hajihassani et al. [16]ANN; ICA-ANNB/S, SL, Qm, E, Vp, D950.911 (ANN);0.976 (ICA-ANN)Hajihassani et al. [17]SVMQm, D800.957Amiri et al. [18]ANN; ANN-KNNQm, D750.82 (ANN); 0.88 (ANN-KNN)Ghasemi et al. [19]ANFIS-PSO; SVMB, S, SL, Qm, D1200.950 (ANFIS-PSO); 0.924 (SVM)Ebrahimi et al. [20]ANNB, S, SL, H, PF340.78Taheri et al. [21]ANN; ABC-ANNQm, D990.91 (ANN); 0.95 (ABC-ANN)Hasanipanah et al. [22]PSOQm, D800.9019Faradonbeh et al. [23]GEPB, S, SL, d, H, PF, Qm, D1150.876Armaghani et al. [24]ICAQm, D730.9334Zhang et al. [25]PSO-XGBoostB, S, PF, Qm, D1750.968Tian et al. [26]GA-linear; GA-powerQm, D850.961 (GA-linear); 0.967 (GA-power)Rezaeineshat et al. [27]ICA-ANNB, S, Qm, D, RQD1120.9002Distance 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).It can be seen from Table1 that intelligent algorithms have been widely used in the prediction of blasting vibration velocity. Research on the combination of various artificial intelligence methods has gradually increased, and the hybrid artificial intelligence algorithms have high accuracy. Most of the above research works are based on a hybrid model composed of ANN algorithm. However, the ANN algorithm has some disadvantages such as being easy to fall into a local optimum, requiring a large number of parameters during calculation, and slow training speed. At present, optimization methods such as balancing composite motion optimization (BCMO) [28] and Gaussian process (GP) have begun to be gradually applied. GP has the advantages of good adaptability, strong generalization ability, simple implementation process, and adaptive acquisition of hyperparameters for processing nonlinear regression problems. At the same time, the differential evolution algorithm (DE) has strong global optimization capabilities and can quickly approach the global optimum. Therefore, it can achieve rapid acquisition of GP hyperparameters by introducing the DE into the GP model and form a new model of PPV prediction based on DE-GP, which will realize the effective control of BIGV.Aiming at the problem of unpredictable PPV in the bustling downtown area, this study uses DE to optimize GP, constructs a model of mixed algorithm DE-GP predicting PPV, and analyzes the accuracy of the model. Finally, the 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 this article is divided as follows: The second part introduces the methods used in blasting prediction, including empirical formulas and the intelligent algorithms studied in this article. The third part introduces the engineering situation and the acquisition of blasting vibration velocity data. The fourth part uses conventional empirical formulas and the proposed intelligent algorithm to predict PPV. The fifth part analyzes the influence of the hyperparameters of the GP model, the selection of DE algorithm parameters, the comparison of the predictive capabilities of different algorithms, and the sensitivity between blasting influencing factors. Finally, the sixth part draws a conclusion.


## 2. Methods Used

### 2.1. Conventional Prediction Methods of Blasting

PPV is the most important indicator to measure the impact of blasting 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 establish a PPV prediction formula that considers all influence factors. A large number of blasting engineering monitors show that the maximum charge per delay and the distance from the blast-face are the two main factors affecting PPV. Therefore, many domestic and foreign scholars have proposed different PPV formulas for predicting BIGV [29–34]. In this study, several commonly used empirical formulas to characterize the attenuation of blasting vibration were selected, as shown in Table 2. In the formula, K, α, β, and n are parameters related to rock properties, blasting operation parameters, and terrain conditions. We often use the field blasting vibration measured data obtained in the actual engineering as the basis and perform regression analysis according to the empirical formulas in Table 2 to obtain the values of the unknown parameters K, α, β, and n.Table 2
Conventional empirical formulas for blasting prediction.
The name of the empirical formulasFormulasUnited States Bureau of Mines (USBM) [29]PPV=KR/Qmax−αLangefors–Kihlstrom (LK) [30]PPV=KQmax/R2/3αGeneral predictor [31]PPV=KR−αQmaxβAmbraseys–Hendron (AH) [32]PPV=KR/Qmax3−αBureau of Indian Standard (BIS) [33]PPV=KQmax/R23αCMRI predictor (CMRI) [34]PPV=n+KR/Qmax−1PPV 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.According to the blasting prediction empirical formulas in Table2, the site constants suitable for this project are fitted by using the actual measurement data of the construction site, and the blasting vibration velocity prediction formulas considering the maximum charge per delay and the distance from the blast-face is constructed. After that, the corresponding PPV can be obtained according to the maximum charge per delay and the distance from the blast-face of a certain blasting construction.


### 2.2. Blasting Vibration Velocity Prediction Model Based on DE-GP

This section selects the blasting parameters according to the influencing factors of blasting vibration. Through a brief description of the GP, a nonlinear prediction model of blasting vibration velocity based on the GP is formed. Limited by the difficulty of selecting hyperparameters for the GP, the DE algorithm is introduced to optimize the GP and applied to the blasting prediction of rock foundation pits. Then, the principle, frame structure, and solution process are described.


#### 2.2.1. Selection of Blasting Prediction Parameters

In the actual blasting of rock foundation pits, there are many factors that affect the blasting effect, and various different factors also cause different degrees of impact. In summary, there are four main aspects: explosive performance, rock properties, explosive and rock-related factors, and blasting technology. Combining the actual engineering situation and related construction experience of the Labor Park Station of Dalian Metro Line 5, the distance from the blast-face, the maximum charge per delay, the total charge, the hole depth, the spacing, the burden, the stemming length, and the powder factor are used as the input parameters of the blasting vibration velocity prediction model, and the peak particle velocity at the monitoring point is used as the output parameter. The above indicators can completely describe the blasting effect of the rock foundation pit of the Suoyuwan South Station of Dalian Metro Line 5 and fully reflect the blasting effect, so as to realize the blasting prediction of PPV.


#### 2.2.2. Nonlinear Prediction Model of PPV Based on GP

GP has good adaptability to the processing of complex classification and regression problems with high dimensionality, small samples, and nonlinearity. Therefore, it is used to establish the blasting PPV prediction model of the rock foundation pit.For a set of known blasting influencing factorsXx1,x2,x3,x4 and corresponding blasting vibration velocity y, the blasting PPV learning samples set D=X,y can be established, and the prior distribution of blasting PPV y can be established through the GP as follows:(1)y∼N0,K+σn2I.In the formula,K=KX,X is a symmetric positive definite covariance matrix of order n×n; and any item Kij in the K=KX,X matrix measures the correlation between xi and xj; 0 is the mean; σn2 is the variance; and I is the identity matrix.From equation (1), the joint GP prior distribution composed of n training sample output y and one test sample output y∗ is derived as follows:(2)yy∗∼N0,KX,X+σn2IKX,x∗KX,x∗kx∗,x∗.In the formula,KX,x∗ is the n×1-order covariance matrix of the blasting test input point x∗ and all input points X in the training set; and kx,x∗ is the covariance of the blasting test input point x∗ itself.Further solving to obtain the mean and variance ofy∗:(3)y^x∗=kTx∗K+σn2I−1y,(4)σ^2x∗=kx∗,x∗−kTx∗K+σn2I−1kx∗.In the formula,kx∗ is the covariance matrix of order n×1, kx∗=KX,x∗, namely kx1,x∗,…,kxn,x∗ΤThe nonlinear prediction model of the blasting PPV is described by the square exponential covariance function:(5)kxp,xq=σf2exp−12l2xp,xq2+σn2δpq.In the formula,σf, l, and σn are hyperparameters. σf2 is the signal variance of the kernel function, which is used to control the degree of local correlation; l is the hyperparameter for correlation measurement, the larger the value, the smaller the correlation between input and output; σn2 is the variance of noise; and δpq is the symbol of Kronecker, when p=q, δpq=0, otherwise δpq=1Hyperparameters can be obtained by maximizing the log-likelihood of the training samples (maximum posterior distribution). The logarithmic form of the likelihood function is as follows:(6)L=lgpy|X=−12yTK+σn2I−1y−12lgK+σn2I−n2lg2π.Obtain the optimal hyperparameters by formula (6), and substitute formula (3) to obtain the blasting PPV prediction value y∗ corresponding to the blasting test sample x∗. According to the establishment steps of the above model, the corresponding calculation program is compiled using MATLAB language. So far, the modeling and regression process of the nonlinear prediction model for blasting PPV of rock foundation pit based on the GP has been completed.


#### 2.2.3. Hybrid DE-Optimized GP and Its Architecture Design

DE is a genetic algorithm based on real number coding. Compared with the selection operation of GA algorithm, DE updates the population according to a one-to-one elimination mechanism and uses directional perturbation on individuals to achieve a reduction in the value of individual functions. In terms of applicability, the GA algorithm and PSO algorithm are mainly used to solve continuous optimization problems, and the ACO algorithm is mainly used to solve discrete optimization and combinatorial optimization problems. However, DE not only does not use the gradient information of the objective function but also has no requirements on the derivability and continuity of the function. So, DE has strong applicability. In terms of parameter setting, DE has two main parameters to be adjusted, and the parameter setting has less obvious influence on the result, so it is easier to use. Compared with the GA algorithm and the PSO algorithm, there are many parameters, and different parameter settings have a greater impact on the final result. Therefore, in actual use, it needs to be adjusted continuously, which increases the difficulty of using the algorithm. When solving multivariable function optimization problems, the GA algorithm converges slowly or even hard for high-dimensional problems, while the PSO algorithm and DE algorithm have the advantages of fast convergence speed and high accuracy. But the PSO algorithm and ACO algorithm are easy to fall into a local optimal solution. In addition, DE considers the correlation between multiple variables to a certain extent and has great advantages in the problem of variable coupling.The study found that the regression results of the nonlinear prediction model for blasting PPV of rock foundation pit based on GP are significantly affected by the hyperparametersσf, l, and σn. In this study, DE is introduced into GP. It can quickly approach the global optimization with the help of DE’s strong global optimization ability. This can realize the optimal selection of three hyperparameters and form a new model for controlling BIGV. RMSE is used as the evaluation index of the prediction results of the GP model, which is used to describe the optimization effect of the GP hyperparameters. Compared with the MSE evaluation index, the calculation result of RMSE can not only keep the error evaluation result consistent with the problem described by GP in the order of magnitude but also intuitively describe the data. Compared with the MAE evaluation index, RMSE is the square root of the sum of squares of errors, which means that the single error has a quadratic increase. It enlarges the gap between errors and makes the error evaluation result more obvious. Therefore, RMSE is used as an evaluation index for the prediction effect of the GP model. The realization process is shown in Figure 1.Figure 1
Flow chart of DE-GP blasting PPV regression model.For the aforementioned GP blasting PPV prediction model, the specific steps of using DE optimization are as follows:Step 1.
Initialize the DE algorithm, and set the DE algorithm’s mutation factorF, cross factor CR, population size NP, and difference strategy.Step 2.
Chooseσf, l, and σn in GP as the optimization target, and randomly generate the initial population in the optimization interval.Step 3.
Use the parameters expressed by each individual in the initial population to perform GP training separately.Step 4.
Evaluate the fitness of each individual in the population. The evaluation index is the root mean square error of the prediction result of the test sample of the trained GP model:(7)F0=1Np∑i=1Npvi−vi′2.
In the formula,vi and vi′ are the actual vibration velocity and the predicted vibration velocity, Np is the number of test samples, and F0 is the individual fitness value.Step 5.
For the minimum fitness valueFmin and expected response value Fbest of the current population, if the current individual satisfies Fbest>Fmin, then enter Step 10; otherwise, enter Step 6.Step 6.
Perform a mutation operation on the current population, and combine the difference values between several individuals in the population with the individuals to be mutated according to certain rules. Then, obtain the mutation vector of each individual in the population.Step 7.
Perform crossover operation on the current population; crossover the mutation vector and its corresponding target vector to obtain the test vector.Step 8.
Select the current population, use the fitness function of the optimized problem to compare each test vector with its corresponding target vector, and put the better one into the next-generation population.Step 9.
Judge whether the maximum number of iterations is reached, if it reaches the maximum number of iterations, enter Step10; otherwise, return to Step 4.Step 10.
Take theσf, l, and σn represented by the minimum fitness individual as the initial parameters of the GP model, so complete the GP optimization process.


## 2.1. Conventional Prediction Methods of Blasting

PPV is the most important indicator to measure the impact of blasting 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 establish a PPV prediction formula that considers all influence factors. A large number of blasting engineering monitors show that the maximum charge per delay and the distance from the blast-face are the two main factors affecting PPV. Therefore, many domestic and foreign scholars have proposed different PPV formulas for predicting BIGV [29–34]. In this study, several commonly used empirical formulas to characterize the attenuation of blasting vibration were selected, as shown in Table 2. In the formula, K, α, β, and n are parameters related to rock properties, blasting operation parameters, and terrain conditions. We often use the field blasting vibration measured data obtained in the actual engineering as the basis and perform regression analysis according to the empirical formulas in Table 2 to obtain the values of the unknown parameters K, α, β, and n.Table 2
Conventional empirical formulas for blasting prediction.
The name of the empirical formulasFormulasUnited States Bureau of Mines (USBM) [29]PPV=KR/Qmax−αLangefors–Kihlstrom (LK) [30]PPV=KQmax/R2/3αGeneral predictor [31]PPV=KR−αQmaxβAmbraseys–Hendron (AH) [32]PPV=KR/Qmax3−αBureau of Indian Standard (BIS) [33]PPV=KQmax/R23αCMRI predictor (CMRI) [34]PPV=n+KR/Qmax−1PPV 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.According to the blasting prediction empirical formulas in Table2, the site constants suitable for this project are fitted by using the actual measurement data of the construction site, and the blasting vibration velocity prediction formulas considering the maximum charge per delay and the distance from the blast-face is constructed. After that, the corresponding PPV can be obtained according to the maximum charge per delay and the distance from the blast-face of a certain blasting construction.


## 2.2. Blasting Vibration Velocity Prediction Model Based on DE-GP

This section selects the blasting parameters according to the influencing factors of blasting vibration. Through a brief description of the GP, a nonlinear prediction model of blasting vibration velocity based on the GP is formed. Limited by the difficulty of selecting hyperparameters for the GP, the DE algorithm is introduced to optimize the GP and applied to the blasting prediction of rock foundation pits. Then, the principle, frame structure, and solution process are described.


### 2.2.1. Selection of Blasting Prediction Parameters

In the actual blasting of rock foundation pits, there are many factors that affect the blasting effect, and various different factors also cause different degrees of impact. In summary, there are four main aspects: explosive performance, rock properties, explosive and rock-related factors, and blasting technology. Combining the actual engineering situation and related construction experience of the Labor Park Station of Dalian Metro Line 5, the distance from the blast-face, the maximum charge per delay, the total charge, the hole depth, the spacing, the burden, the stemming length, and the powder factor are used as the input parameters of the blasting vibration velocity prediction model, and the peak particle velocity at the monitoring point is used as the output parameter. The above indicators can completely describe the blasting effect of the rock foundation pit of the Suoyuwan South Station of Dalian Metro Line 5 and fully reflect the blasting effect, so as to realize the blasting prediction of PPV.


### 2.2.2. Nonlinear Prediction Model of PPV Based on GP

GP has good adaptability to the processing of complex classification and regression problems with high dimensionality, small samples, and nonlinearity. Therefore, it is used to establish the blasting PPV prediction model of the rock foundation pit.For a set of known blasting influencing factorsXx1,x2,x3,x4 and corresponding blasting vibration velocity y, the blasting PPV learning samples set D=X,y can be established, and the prior distribution of blasting PPV y can be established through the GP as follows:(1)y∼N0,K+σn2I.In the formula,K=KX,X is a symmetric positive definite covariance matrix of order n×n; and any item Kij in the K=KX,X matrix measures the correlation between xi and xj; 0 is the mean; σn2 is the variance; and I is the identity matrix.From equation (1), the joint GP prior distribution composed of n training sample output y and one test sample output y∗ is derived as follows:(2)yy∗∼N0,KX,X+σn2IKX,x∗KX,x∗kx∗,x∗.In the formula,KX,x∗ is the n×1-order covariance matrix of the blasting test input point x∗ and all input points X in the training set; and kx,x∗ is the covariance of the blasting test input point x∗ itself.Further solving to obtain the mean and variance ofy∗:(3)y^x∗=kTx∗K+σn2I−1y,(4)σ^2x∗=kx∗,x∗−kTx∗K+σn2I−1kx∗.In the formula,kx∗ is the covariance matrix of order n×1, kx∗=KX,x∗, namely kx1,x∗,…,kxn,x∗ΤThe nonlinear prediction model of the blasting PPV is described by the square exponential covariance function:(5)kxp,xq=σf2exp−12l2xp,xq2+σn2δpq.In the formula,σf, l, and σn are hyperparameters. σf2 is the signal variance of the kernel function, which is used to control the degree of local correlation; l is the hyperparameter for correlation measurement, the larger the value, the smaller the correlation between input and output; σn2 is the variance of noise; and δpq is the symbol of Kronecker, when p=q, δpq=0, otherwise δpq=1Hyperparameters can be obtained by maximizing the log-likelihood of the training samples (maximum posterior distribution). The logarithmic form of the likelihood function is as follows:(6)L=lgpy|X=−12yTK+σn2I−1y−12lgK+σn2I−n2lg2π.Obtain the optimal hyperparameters by formula (6), and substitute formula (3) to obtain the blasting PPV prediction value y∗ corresponding to the blasting test sample x∗. According to the establishment steps of the above model, the corresponding calculation program is compiled using MATLAB language. So far, the modeling and regression process of the nonlinear prediction model for blasting PPV of rock foundation pit based on the GP has been completed.


### 2.2.3. Hybrid DE-Optimized GP and Its Architecture Design

DE is a genetic algorithm based on real number coding. Compared with the selection operation of GA algorithm, DE updates the population according to a one-to-one elimination mechanism and uses directional perturbation on individuals to achieve a reduction in the value of individual functions. In terms of applicability, the GA algorithm and PSO algorithm are mainly used to solve continuous optimization problems, and the ACO algorithm is mainly used to solve discrete optimization and combinatorial optimization problems. However, DE not only does not use the gradient information of the objective function but also has no requirements on the derivability and continuity of the function. So, DE has strong applicability. In terms of parameter setting, DE has two main parameters to be adjusted, and the parameter setting has less obvious influence on the result, so it is easier to use. Compared with the GA algorithm and the PSO algorithm, there are many parameters, and different parameter settings have a greater impact on the final result. Therefore, in actual use, it needs to be adjusted continuously, which increases the difficulty of using the algorithm. When solving multivariable function optimization problems, the GA algorithm converges slowly or even hard for high-dimensional problems, while the PSO algorithm and DE algorithm have the advantages of fast convergence speed and high accuracy. But the PSO algorithm and ACO algorithm are easy to fall into a local optimal solution. In addition, DE considers the correlation between multiple variables to a certain extent and has great advantages in the problem of variable coupling.The study found that the regression results of the nonlinear prediction model for blasting PPV of rock foundation pit based on GP are significantly affected by the hyperparametersσf, l, and σn. In this study, DE is introduced into GP. It can quickly approach the global optimization with the help of DE’s strong global optimization ability. This can realize the optimal selection of three hyperparameters and form a new model for controlling BIGV. RMSE is used as the evaluation index of the prediction results of the GP model, which is used to describe the optimization effect of the GP hyperparameters. Compared with the MSE evaluation index, the calculation result of RMSE can not only keep the error evaluation result consistent with the problem described by GP in the order of magnitude but also intuitively describe the data. Compared with the MAE evaluation index, RMSE is the square root of the sum of squares of errors, which means that the single error has a quadratic increase. It enlarges the gap between errors and makes the error evaluation result more obvious. Therefore, RMSE is used as an evaluation index for the prediction effect of the GP model. The realization process is shown in Figure 1.Figure 1
Flow chart of DE-GP blasting PPV regression model.For the aforementioned GP blasting PPV prediction model, the specific steps of using DE optimization are as follows:Step 1.
Initialize the DE algorithm, and set the DE algorithm’s mutation factorF, cross factor CR, population size NP, and difference strategy.Step 2.
Chooseσf, l, and σn in GP as the optimization target, and randomly generate the initial population in the optimization interval.Step 3.
Use the parameters expressed by each individual in the initial population to perform GP training separately.Step 4.
Evaluate the fitness of each individual in the population. The evaluation index is the root mean square error of the prediction result of the test sample of the trained GP model:(7)F0=1Np∑i=1Npvi−vi′2.
In the formula,vi and vi′ are the actual vibration velocity and the predicted vibration velocity, Np is the number of test samples, and F0 is the individual fitness value.Step 5.
For the minimum fitness valueFmin and expected response value Fbest of the current population, if the current individual satisfies Fbest>Fmin, then enter Step 10; otherwise, enter Step 6.Step 6.
Perform a mutation operation on the current population, and combine the difference values between several individuals in the population with the individuals to be mutated according to certain rules. Then, obtain the mutation vector of each individual in the population.Step 7.
Perform crossover operation on the current population; crossover the mutation vector and its corresponding target vector to obtain the test vector.Step 8.
Select the current population, use the fitness function of the optimized problem to compare each test vector with its corresponding target vector, and put the better one into the next-generation population.Step 9.
Judge whether the maximum number of iterations is reached, if it reaches the maximum number of iterations, enter Step10; otherwise, return to Step 4.Step 10.
Take theσf, l, and σn represented by the minimum fitness individual as the initial parameters of the GP model, so complete the GP optimization process.


## 2.2.1. Selection of Blasting Prediction Parameters

In the actual blasting of rock foundation pits, there are many factors that affect the blasting effect, and various different factors also cause different degrees of impact. In summary, there are four main aspects: explosive performance, rock properties, explosive and rock-related factors, and blasting technology. Combining the actual engineering situation and related construction experience of the Labor Park Station of Dalian Metro Line 5, the distance from the blast-face, the maximum charge per delay, the total charge, the hole depth, the spacing, the burden, the stemming length, and the powder factor are used as the input parameters of the blasting vibration velocity prediction model, and the peak particle velocity at the monitoring point is used as the output parameter. The above indicators can completely describe the blasting effect of the rock foundation pit of the Suoyuwan South Station of Dalian Metro Line 5 and fully reflect the blasting effect, so as to realize the blasting prediction of PPV.


## 2.2.2. Nonlinear Prediction Model of PPV Based on GP

GP has good adaptability to the processing of complex classification and regression problems with high dimensionality, small samples, and nonlinearity. Therefore, it is used to establish the blasting PPV prediction model of the rock foundation pit.For a set of known blasting influencing factorsXx1,x2,x3,x4 and corresponding blasting vibration velocity y, the blasting PPV learning samples set D=X,y can be established, and the prior distribution of blasting PPV y can be established through the GP as follows:(1)y∼N0,K+σn2I.In the formula,K=KX,X is a symmetric positive definite covariance matrix of order n×n; and any item Kij in the K=KX,X matrix measures the correlation between xi and xj; 0 is the mean; σn2 is the variance; and I is the identity matrix.From equation (1), the joint GP prior distribution composed of n training sample output y and one test sample output y∗ is derived as follows:(2)yy∗∼N0,KX,X+σn2IKX,x∗KX,x∗kx∗,x∗.In the formula,KX,x∗ is the n×1-order covariance matrix of the blasting test input point x∗ and all input points X in the training set; and kx,x∗ is the covariance of the blasting test input point x∗ itself.Further solving to obtain the mean and variance ofy∗:(3)y^x∗=kTx∗K+σn2I−1y,(4)σ^2x∗=kx∗,x∗−kTx∗K+σn2I−1kx∗.In the formula,kx∗ is the covariance matrix of order n×1, kx∗=KX,x∗, namely kx1,x∗,…,kxn,x∗ΤThe nonlinear prediction model of the blasting PPV is described by the square exponential covariance function:(5)kxp,xq=σf2exp−12l2xp,xq2+σn2δpq.In the formula,σf, l, and σn are hyperparameters. σf2 is the signal variance of the kernel function, which is used to control the degree of local correlation; l is the hyperparameter for correlation measurement, the larger the value, the smaller the correlation between input and output; σn2 is the variance of noise; and δpq is the symbol of Kronecker, when p=q, δpq=0, otherwise δpq=1Hyperparameters can be obtained by maximizing the log-likelihood of the training samples (maximum posterior distribution). The logarithmic form of the likelihood function is as follows:(6)L=lgpy|X=−12yTK+σn2I−1y−12lgK+σn2I−n2lg2π.Obtain the optimal hyperparameters by formula (6), and substitute formula (3) to obtain the blasting PPV prediction value y∗ corresponding to the blasting test sample x∗. According to the establishment steps of the above model, the corresponding calculation program is compiled using MATLAB language. So far, the modeling and regression process of the nonlinear prediction model for blasting PPV of rock foundation pit based on the GP has been completed.


## 2.2.3. Hybrid DE-Optimized GP and Its Architecture Design

DE is a genetic algorithm based on real number coding. Compared with the selection operation of GA algorithm, DE updates the population according to a one-to-one elimination mechanism and uses directional perturbation on individuals to achieve a reduction in the value of individual functions. In terms of applicability, the GA algorithm and PSO algorithm are mainly used to solve continuous optimization problems, and the ACO algorithm is mainly used to solve discrete optimization and combinatorial optimization problems. However, DE not only does not use the gradient information of the objective function but also has no requirements on the derivability and continuity of the function. So, DE has strong applicability. In terms of parameter setting, DE has two main parameters to be adjusted, and the parameter setting has less obvious influence on the result, so it is easier to use. Compared with the GA algorithm and the PSO algorithm, there are many parameters, and different parameter settings have a greater impact on the final result. Therefore, in actual use, it needs to be adjusted continuously, which increases the difficulty of using the algorithm. When solving multivariable function optimization problems, the GA algorithm converges slowly or even hard for high-dimensional problems, while the PSO algorithm and DE algorithm have the advantages of fast convergence speed and high accuracy. But the PSO algorithm and ACO algorithm are easy to fall into a local optimal solution. In addition, DE considers the correlation between multiple variables to a certain extent and has great advantages in the problem of variable coupling.The study found that the regression results of the nonlinear prediction model for blasting PPV of rock foundation pit based on GP are significantly affected by the hyperparametersσf, l, and σn. In this study, DE is introduced into GP. It can quickly approach the global optimization with the help of DE’s strong global optimization ability. This can realize the optimal selection of three hyperparameters and form a new model for controlling BIGV. RMSE is used as the evaluation index of the prediction results of the GP model, which is used to describe the optimization effect of the GP hyperparameters. Compared with the MSE evaluation index, the calculation result of RMSE can not only keep the error evaluation result consistent with the problem described by GP in the order of magnitude but also intuitively describe the data. Compared with the MAE evaluation index, RMSE is the square root of the sum of squares of errors, which means that the single error has a quadratic increase. It enlarges the gap between errors and makes the error evaluation result more obvious. Therefore, RMSE is used as an evaluation index for the prediction effect of the GP model. The realization process is shown in Figure 1.Figure 1
Flow chart of DE-GP blasting PPV regression model.For the aforementioned GP blasting PPV prediction model, the specific steps of using DE optimization are as follows:Step 1.
Initialize the DE algorithm, and set the DE algorithm’s mutation factorF, cross factor CR, population size NP, and difference strategy.Step 2.
Chooseσf, l, and σn in GP as the optimization target, and randomly generate the initial population in the optimization interval.Step 3.
Use the parameters expressed by each individual in the initial population to perform GP training separately.Step 4.
Evaluate the fitness of each individual in the population. The evaluation index is the root mean square error of the prediction result of the test sample of the trained GP model:(7)F0=1Np∑i=1Npvi−vi′2.
In the formula,vi and vi′ are the actual vibration velocity and the predicted vibration velocity, Np is the number of test samples, and F0 is the individual fitness value.Step 5.
For the minimum fitness valueFmin and expected response value Fbest of the current population, if the current individual satisfies Fbest>Fmin, then enter Step 10; otherwise, enter Step 6.Step 6.
Perform a mutation operation on the current population, and combine the difference values between several individuals in the population with the individuals to be mutated according to certain rules. Then, obtain the mutation vector of each individual in the population.Step 7.
Perform crossover operation on the current population; crossover the mutation vector and its corresponding target vector to obtain the test vector.Step 8.
Select the current population, use the fitness function of the optimized problem to compare each test vector with its corresponding target vector, and put the better one into the next-generation population.Step 9.
Judge whether the maximum number of iterations is reached, if it reaches the maximum number of iterations, enter Step10; otherwise, return to Step 4.Step 10.
Take theσf, l, and σn represented by the minimum fitness individual as the initial parameters of the GP model, so complete the GP optimization process.


## 3. Data Used

### 3.1. Description of Study Area

Dalian Metro Line 5 has a total length of 23.8 km and has 18 stations. The entire line was laid by underground lines with an average distance of 1.38 km. Labor Park Station is located at the intersection of Jiefang Road and Ziwei Street in Zhongshan District. The main body was arranged along the north and south of Jiefang Road. The starting mileage of the station was ZK7 + 532.919 (YK7 + 531.658), the ending mileage was K7 + 722.859, the length was 191.2 m, and the width of standard section was 23.3 m. The main body of the station was an underground two-story island station with a platform width of 14 m, and the station roof was covered with soil about 18.7–24.0 m. The 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 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 the west side of the site, and the schematic diagram of the project is shown in Figure2. The minimum distance between the main structure of the station and the library of the 24th Middle School was 5.9 m, and the minimum distance from 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 a dense area of pedestrians and vehicles.Figure 2
Sketch map of project overview.The main part of the station was constructed by the PBA construction method, and two construction foundation pits were set up. The foundation pits were strictly in accordance with the principle of layered, segmented, symmetrical, balanced, and timely excavation, and supports were erected in time. The foundation pit adopted a combined support 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 16 mm. The initial support system of the concealed excavation section was advanced small pipe grouting to reinforce the ground or deep hole grouting, grid steel frame, and mesh shotcrete support. Because most of the surrounding rock during construction was in strongly weathered slate and moderately weathered slate, in order to ensure the safety of ground buildings, underground structures, and the surrounding environment during drilling and blasting, we should adopt microseismic controlled blasting construction 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, and the location of excavation, different blasting parameters were adopted and adjusted according to the actual blasting effect during the blasting construction.The foundation pit of Labor Park adopted a triangular hole pattern and vertical drilling. The detonation adopted a millisecond-delay nonelectrically conductive detonation system, and the method of detonation outside the hole with slight difference, cluster connection method, five holes with one sound. Used ms-8 in the hole, ms-4 for the connection between the surface holes, and ms-5 for the connection between the surface rows. During construction, the blast hole depth of each layer was 1.4 m∼2.8 m, and the blast hole diameter was 42 mm. The blast hole spacing was 0.6 m∼1.5 m, and the blast hole burden was 0.5 m∼1.3 m. The rock drill was used for drilling, and the explosive was 2# rock emulsion explosive, with a specification of Ф32 mm ∗ 200 cm, 0.2 kg or 0.3 kg each, and the maximum charge per delay was 1.4 kg. During the blasting process, it is necessary to control the maximum charge per delay, reduce the number of holes in the same section, and determine reasonable blasting parameters. In addition, measures such as excavation of damping trenches, covering, and protection were taken in special areas to reduce the impact of blasting on surrounding buildings.


### 3.2. Monitoring and Samples Construction

Combined with the actual situation on site, this monitoring adopted the TC-4850N intelligent blasting vibration monitor, which is a portable remote sensing vibration monitor with four parallel channels specially designed for engineering blasting. It has the advantages of lightweight, waterproof, dustproof, high pressure resistance, high precision, and wide application. Figure3 shows the working diagram of the TC-4850N wireless blasting monitor.Figure 3
Schematic diagram of TC-4850N wireless blasting monitor.Several blasting vibration test monitoring points were arranged near the blasting area of the foundation pit of Suoyuwan South Station. Then, begin to monitor and record the blasting vibration velocity at different locations and different blasting parameters during the initial stage of excavation, as shown in Figure4.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.
(a)(b)(c)(d)By on-site monitoring, the monitoring data of blasting vibration at different positions and different blasting parameters during the excavation of the foundation pit were obtained. There are 88 sets of test data, as shown in Table3. This paper takes the distance from the blast-face, the 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 the output variable to build the samples. The first 72 groups are training samples, and the last 16 groups are test samples.Table 3
Model statistical parameters of training samples and test samples.
No.ParameterUnitSymbolCategoryMinMaxMeanSDMedianTraining data (72 data)1Distance from the blast-facemDInput6.551.927.9612.6828.452Maximum charge per delaykgQmInput0.31.350.730.290.683Total chargekgQtInput11.477.242.718.342.44Hole depthmHInput1.42.82.10.42.15SpacingmSInput0.61.51.10.21.16BurdenmBInput0.51.30.90.21.07Stemming lengthmSLInput0.81.71.20.31.28Powder factorkg/m3PFInput0.310.820.580.140.599Peak particle velocitycm/sPPVOutput0.173.471.280.900.96Test data (16 data)1DistancemDInput7.347.127.211.326.62Maximum charge per delaykgQmInput0.301.200.770.270.803Total chargekgQtInput12.173.841.818.035.94Hole depthmHInput1.52.62.00.42.15SpacingmSInput0.81.41.10.21.16BurdenmBInput0.61.31.00.21.07Stemming lengthmSLInput0.81.61.10.31.28Powder factorkg/m3PFInput0.330.790.560.140.549Peak particle velocitycm/sPPVOutput0.732.081.370.421.30


## 3.1. Description of Study Area

Dalian Metro Line 5 has a total length of 23.8 km and has 18 stations. The entire line was laid by underground lines with an average distance of 1.38 km. Labor Park Station is located at the intersection of Jiefang Road and Ziwei Street in Zhongshan District. The main body was arranged along the north and south of Jiefang Road. The starting mileage of the station was ZK7 + 532.919 (YK7 + 531.658), the ending mileage was K7 + 722.859, the length was 191.2 m, and the width of standard section was 23.3 m. The main body of the station was an underground two-story island station with a platform width of 14 m, and the station roof was covered with soil about 18.7–24.0 m. The 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 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 the west side of the site, and the schematic diagram of the project is shown in Figure2. The minimum distance between the main structure of the station and the library of the 24th Middle School was 5.9 m, and the minimum distance from 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 a dense area of pedestrians and vehicles.Figure 2
Sketch map of project overview.The main part of the station was constructed by the PBA construction method, and two construction foundation pits were set up. The foundation pits were strictly in accordance with the principle of layered, segmented, symmetrical, balanced, and timely excavation, and supports were erected in time. The foundation pit adopted a combined support 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 16 mm. The initial support system of the concealed excavation section was advanced small pipe grouting to reinforce the ground or deep hole grouting, grid steel frame, and mesh shotcrete support. Because most of the surrounding rock during construction was in strongly weathered slate and moderately weathered slate, in order to ensure the safety of ground buildings, underground structures, and the surrounding environment during drilling and blasting, we should adopt microseismic controlled blasting construction 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, and the location of excavation, different blasting parameters were adopted and adjusted according to the actual blasting effect during the blasting construction.The foundation pit of Labor Park adopted a triangular hole pattern and vertical drilling. The detonation adopted a millisecond-delay nonelectrically conductive detonation system, and the method of detonation outside the hole with slight difference, cluster connection method, five holes with one sound. Used ms-8 in the hole, ms-4 for the connection between the surface holes, and ms-5 for the connection between the surface rows. During construction, the blast hole depth of each layer was 1.4 m∼2.8 m, and the blast hole diameter was 42 mm. The blast hole spacing was 0.6 m∼1.5 m, and the blast hole burden was 0.5 m∼1.3 m. The rock drill was used for drilling, and the explosive was 2# rock emulsion explosive, with a specification of Ф32 mm ∗ 200 cm, 0.2 kg or 0.3 kg each, and the maximum charge per delay was 1.4 kg. During the blasting process, it is necessary to control the maximum charge per delay, reduce the number of holes in the same section, and determine reasonable blasting parameters. In addition, measures such as excavation of damping trenches, covering, and protection were taken in special areas to reduce the impact of blasting on surrounding buildings.


## 3.2. Monitoring and Samples Construction

Combined with the actual situation on site, this monitoring adopted the TC-4850N intelligent blasting vibration monitor, which is a portable remote sensing vibration monitor with four parallel channels specially designed for engineering blasting. It has the advantages of lightweight, waterproof, dustproof, high pressure resistance, high precision, and wide application. Figure3 shows the working diagram of the TC-4850N wireless blasting monitor.Figure 3
Schematic diagram of TC-4850N wireless blasting monitor.Several blasting vibration test monitoring points were arranged near the blasting area of the foundation pit of Suoyuwan South Station. Then, begin to monitor and record the blasting vibration velocity at different locations and different blasting parameters during the initial stage of excavation, as shown in Figure4.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.
(a)(b)(c)(d)By on-site monitoring, the monitoring data of blasting vibration at different positions and different blasting parameters during the excavation of the foundation pit were obtained. There are 88 sets of test data, as shown in Table3. This paper takes the distance from the blast-face, the 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 the output variable to build the samples. The first 72 groups are training samples, and the last 16 groups are test samples.Table 3
Model statistical parameters of training samples and test samples.
No.ParameterUnitSymbolCategoryMinMaxMeanSDMedianTraining data (72 data)1Distance from the blast-facemDInput6.551.927.9612.6828.452Maximum charge per delaykgQmInput0.31.350.730.290.683Total chargekgQtInput11.477.242.718.342.44Hole depthmHInput1.42.82.10.42.15SpacingmSInput0.61.51.10.21.16BurdenmBInput0.51.30.90.21.07Stemming lengthmSLInput0.81.71.20.31.28Powder factorkg/m3PFInput0.310.820.580.140.599Peak particle velocitycm/sPPVOutput0.173.471.280.900.96Test data (16 data)1DistancemDInput7.347.127.211.326.62Maximum charge per delaykgQmInput0.301.200.770.270.803Total chargekgQtInput12.173.841.818.035.94Hole depthmHInput1.52.62.00.42.15SpacingmSInput0.81.41.10.21.16BurdenmBInput0.61.31.00.21.07Stemming lengthmSLInput0.81.61.10.31.28Powder factorkg/m3PFInput0.330.790.560.140.549Peak particle velocitycm/sPPVOutput0.732.081.370.421.30


## 4. Results

### 4.1. Prediction Result Analysis of PPV Using Conventional Predictors

According to the previous discussion, we can use the empirical formulas provided in Table2 to predict the blasting PPV of the monitoring points around the rock foundation pit. These empirical formulas are mainly determined by the distance from the blast-face and the maximum charge per delay. They belong to the nonlinear blasting PPV prediction model of rock foundation pits. In this study, 72 sets of blasting training samples were used for regression analysis to determine the rock properties and terrain conditions in the commonly used blasting PPV prediction formula, as shown in Figure 5. The site constants in the various blasting PPV prediction formulas are shown in Table 4, where PPV is the peak particle velocity (cm/s) at the monitoring point, and Qmax is the maximum charge per delay (kg), R is the distance from the blast-face (m), and K, α, β, and n are site constants related to rock properties and geological conditions.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.
(a)(b)(c)(d)(e)(f)Table 4
The conventional empirical formulas and their site constants.
NameEquationSite constantsKαβnR2USBMPPV=KR/Qmax−α7.7070.5410.319LKPPV=KQmax/R2/3α5.4711.1790.261General predictorPPV=KR−αQmaxβ4.1850.3020.6430.213AHPPV=KR/Qmax3−α7.6570.5490.302BISPPV=KQmax/R23α5.4710.5890.26CMRIPPV=n+KR/Qmax−115.5810.6490.257After determining the site constants of the conventional empirical formulas, we use 16 sets of blasting test samples to measure the obtained empirical formula. The relationship between the predicted PPV and the measured PPV of different empirical formulas on the equality line with a slope of 1 is shown in Figure6.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) Relationship 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.
(a)(b)(c)(d)(e)(f)By comparing the coefficient of determination (R2), it is found that the R2 of the above empirical formulas varies between 0.219 and 0.617. Among them, the empirical formula of USBM has higher prediction accuracy for PPV, with the largest coefficient of determination, and the predicted PPV is the closest to the actual measured PPV. The second is the empirical formulas of AH and CMRI, and the prediction accuracy is 0.614 and 0.565, respectively. However, the general predict empirical formula has the lowest prediction accuracy, only 0.219, and the coefficient of determination is the smallest.


### 4.2. Prediction Result Analysis of PPV Using DE-GP

As mentioned earlier, about 80% of the blasting data samples (72 sets of data) are regarded as training samples, which are used to establish a nonlinear mapping relationship for predicting PPV; the remaining 20% of the blasting data samples (16 sets of data) are regarded as test samples.During the blasting construction of rock foundation pits, the DE-GP model is used to predict PPV. By learning the training samples, taking the minimum fitness value of the population as the goal, the DE algorithm is used for global optimization, and the optimal hyperparameters in the GP are quickly searched. Among them, set the population sizeNP = 60 and the variation factor F, crossover factor CR, and maximum number of iterations are 0.7, 0.8, and 150, respectively. According to the application experience of machine 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 target value between iteration steps is less than 1.0 × 10−6 [35–37]. Through the learning and prediction of the DE-GP model, the relationship between the actual measured PPV and the predicted PPV in the training phase and the test phase is shown in Figure 7.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.
(a)(b)It can be seen from Figure7 that the DE-GP model is used to predict the training samples and test samples, and R2 is 0.937 and 0.955, respectively. It can be seen from the measured PPV and the predicted PPV that the prediction effect is good, and the predicted PPV is in good agreement with the actual PPV.


## 4.1. Prediction Result Analysis of PPV Using Conventional Predictors

According to the previous discussion, we can use the empirical formulas provided in Table2 to predict the blasting PPV of the monitoring points around the rock foundation pit. These empirical formulas are mainly determined by the distance from the blast-face and the maximum charge per delay. They belong to the nonlinear blasting PPV prediction model of rock foundation pits. In this study, 72 sets of blasting training samples were used for regression analysis to determine the rock properties and terrain conditions in the commonly used blasting PPV prediction formula, as shown in Figure 5. The site constants in the various blasting PPV prediction formulas are shown in Table 4, where PPV is the peak particle velocity (cm/s) at the monitoring point, and Qmax is the maximum charge per delay (kg), R is the distance from the blast-face (m), and K, α, β, and n are site constants related to rock properties and geological conditions.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.
(a)(b)(c)(d)(e)(f)Table 4
The conventional empirical formulas and their site constants.
NameEquationSite constantsKαβnR2USBMPPV=KR/Qmax−α7.7070.5410.319LKPPV=KQmax/R2/3α5.4711.1790.261General predictorPPV=KR−αQmaxβ4.1850.3020.6430.213AHPPV=KR/Qmax3−α7.6570.5490.302BISPPV=KQmax/R23α5.4710.5890.26CMRIPPV=n+KR/Qmax−115.5810.6490.257After determining the site constants of the conventional empirical formulas, we use 16 sets of blasting test samples to measure the obtained empirical formula. The relationship between the predicted PPV and the measured PPV of different empirical formulas on the equality line with a slope of 1 is shown in Figure6.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) Relationship 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.
(a)(b)(c)(d)(e)(f)By comparing the coefficient of determination (R2), it is found that the R2 of the above empirical formulas varies between 0.219 and 0.617. Among them, the empirical formula of USBM has higher prediction accuracy for PPV, with the largest coefficient of determination, and the predicted PPV is the closest to the actual measured PPV. The second is the empirical formulas of AH and CMRI, and the prediction accuracy is 0.614 and 0.565, respectively. However, the general predict empirical formula has the lowest prediction accuracy, only 0.219, and the coefficient of determination is the smallest.


## 4.2. Prediction Result Analysis of PPV Using DE-GP

As mentioned earlier, about 80% of the blasting data samples (72 sets of data) are regarded as training samples, which are used to establish a nonlinear mapping relationship for predicting PPV; the remaining 20% of the blasting data samples (16 sets of data) are regarded as test samples.During the blasting construction of rock foundation pits, the DE-GP model is used to predict PPV. By learning the training samples, taking the minimum fitness value of the population as the goal, the DE algorithm is used for global optimization, and the optimal hyperparameters in the GP are quickly searched. Among them, set the population sizeNP = 60 and the variation factor F, crossover factor CR, and maximum number of iterations are 0.7, 0.8, and 150, respectively. According to the application experience of machine 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 target value between iteration steps is less than 1.0 × 10−6 [35–37]. Through the learning and prediction of the DE-GP model, the relationship between the actual measured PPV and the predicted PPV in the training phase and the test phase is shown in Figure 7.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.
(a)(b)It can be seen from Figure7 that the DE-GP model is used to predict the training samples and test samples, and R2 is 0.937 and 0.955, respectively. It can be seen from the measured PPV and the predicted PPV that the prediction effect is good, and the predicted PPV is in good agreement with the actual PPV.


## 5. Discussion

### 5.1. Influence of GP Parameters on Regression Accuracy

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


### 5.2. Influence of DE Parameters on Optimization Effect

When using the DE-GP model, there are many factors that affect the optimization effect of the DE, and they are directly related to the selection of hyperparameters in the GP model.The control parameters in the DE mainly include mutation factorF, crossover factor CR, and population size NP. Usually, NP is 5∼10D (D stands for dimensionality), F is 0.5∼0.9, and CR is 0.5∼0.9. It can not only ensure a higher success rate of optimization but also increase the speed of convergence. The mutation factor F and the cross factor CR have a greater impact on the speed and stability of the search. Improper selection of parameters may lead to evolutionary stagnation and premature maturity. The mutation factor F mainly affects the search steps. Increasing F can increase the search range of the DE and increase the diversity of the population but, at the same time, weaken the convergence speed of the algorithm; reducing F can increase the development ability of the DE and improve the convergence speed of the DE, and also it is easy to make the population converge prematurely. The cross factor CR affects the weight of evolutionary information. Increasing CR can improve the diversity of the population and speed up the convergence speed of the DE; reducing CR is beneficial to the analysis of individual separable problems. Select different variation factor F and cross factor CR to perform optimization search trial calculation, and the corresponding iterative convergence curve is shown in Figure 9.Figure 9
Iterative curves of different variation factors and crossover factors. (a) When CR = 0.7, iterative curves of different variation factorsF. (b) When F = 0.8, iterative curves of different crossover factors CR.
(a)(b)It can be seen from Figure9 that when CR = 0.7 and F ∈ [0.5, 0.9], DE iterative search can converge, but the iteration speed is different. When F = 0.8, the convergence speed is the fastest, the number of iteration steps is less, and the iteration steps reach 54 steps to complete convergence. When F = 0.6, the convergence speed is the slowest, and the iterative steps reach 85 steps before it converges completely. Therefore, choosing appropriate search parameters can greatly save 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 reaches the fastest convergence, and it is completely converged in 39 iterations. When CR = 0.6, the convergence 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 comparisons should be made within the selection range of possible convergence parameters in order to obtain the optimal convergence result.When the variation factorF = 0.8 and the crossover factor CR = 0.9 are fixed, increasing NP can improve the diversity of the population but, at the same time, reduce the convergence speed of the population. Reducing NP can increase the convergence speed, but it is easy to lead to premature convergence. When the population number NP is equal to 40, 60, 80, 100, and 120, respectively, the iterative convergence curves are shown in Figure 10.Figure 10
Iterative convergence curves of differentNP.It can be seen from Figure10 that when NP is 40 and 60, the optimal value is searched for 54 iterations and 48 iterations, respectively. But when NP is 80, 100, and 120, it takes 57, 73, and 85 iterations to search for the optimal value. NP mainly reflects the amount of population information in DE. The larger the population size, the richer the population information contained, it is more likely to find the global solution, and it is less likely for individual optimal solutions to dominate the overall solution. However, the calculation time is relatively long, and the calculation convergence is relatively slow, which lead to an increase in the calculation cost of a single population. On the contrary, the diversity of the population is restricted, which is not conducive to DE obtaining the global optimal solution and may even cause the search to stagnate. It can also be seen from Figure 10 that the fitness convergence speed does not decrease monotonically with the increase of NP. When NP is set to 60, convergence can be achieved at the fastest speed. So, NP is set to 60.Choosing appropriate initial parameters can not only save time and improve convergence speed but also avoid the algorithm from falling into local optimality. If the selection of parameters is not appropriate, it may cause the algorithm search to stagnate due to excessive emphasis on the search ability or cause the algorithm to converge prematurely due to excessive emphasis on the development ability. In this research, the final DE algorithm control parameters settings:NP = 60, F = 0.8, and CR = 0.9.


### 5.3. Computing Ability Evaluation of Different Prediction Methods

Since the constructed blasting PPV model will directly affect the accuracy of the prediction results, it is necessary to establish the evaluation index of the model to test the prediction effect. From a statistical point of view, it is more one sided to comprehensively evaluate a predictive model with only one performance index. Therefore, in this research, statistical indicators such as coefficient of determination (R2), root mean square error (RMSE), mean square error (MSE), mean absolute error (MAE), mean absolute percentage error (MAPE), variance account for (VAF), and standard deviation (SD) are used to evaluate the performance of the model, as shown in Table 5.Table 5
Statistical indicators of blasting PPV prediction effect.
Statistical indicatorsFormulasDesired valueCoefficient of determinationR2=1−∑i=1myi−yi′2/∑i=1myi−yi′21Root mean square errorRMSE=1/m∑i=1myi−yi′20Mean square errorMSE=1/m∑i=1myi−yi′20Mean absolute errorMAE=1/m∑i=1myi−yi′0Mean absolute percentage errorMAPE=1/m∑i=1myi−yi′/yi×100%0Variance account forVAF=1−varyi−yi′/varyi×100%100Standard deviationSD=1/m∑i=1myi−1/my1+⋯+ym20Among them,yi is the measured PPV, yi is the predicted PPV, y¯ is the average value of the measured PPV, and m is the total number of samples. It should be pointed out that the closer R2 is to 1, the closer the values of RMSE, MSE, MAE, MAPE, and SD are to 0, and the closer the VAF is to 100, indicating that the performance of the prediction models is better.According to the site constants in the traditional empirical formulas obtained above, 16 sets of test samples are used to test the predictive ability of the above several commonly used empirical formulas. Each index is shown in Table6.Table 6
Predictive effects of several empirical formulas.
NameR2RMSEMSEMAEMAPEVAFSDRankUSBM0.620.260.070.2115.2760.490.531LK0.390.350.120.2622.1533.320.694General predict0.220.410.160.326.2412.610.836AH0.610.280.080.2115.8559.720.572BIS0.370.350.120.2724.6227.510.755CMRI0.570.30.090.2317.1855.790.643It can be seen from Table6 that when using the USBM empirical formula, R2 = 0.62, RMSE = 0.26, MSE = 0.07, 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 selected among the six empirical formulas for comparison with other blasting intelligent algorithms.In this research, the performance of the hybrid DE-GP prediction model is compared with LSSVM, ANN, and GP models and USBM empirical formula. In the ANN model, the number of hidden layer nodesNf of the algorithm is 5, the number of iterations is Ep = 100, and the learning rate is Lr = 0.1. In the GP model, the variance of noise is σn2 = 3.72, the signal variance of the kernel function is σf = −1.7, and the hyperparameter for correlation measurement is l = 0.3. In the core parameters of the LSSVM model, the regular parameter is γ = 9.54, and square bandwidth is σ2 = 2.17. In the DE-GP model, the population number NP is 60, mutation factor F is 0.8, crossover factor CR is 0.9, and the upper limit of evolutionary algebra is 100. The relationship between measured PPV and predicted PPV using LSSVM, ANN, GP, and DE-GP models is shown in Figure 11.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.
(a)(b)(c)(d)It can be seen from Figure11 that using the ANN model to predict the PPV of rock foundation pits has the lowest accuracy, R2 is 0.724. While using the LSSVM model and the GP model, R2 is 0.832 and 0.885, respectively. After finally using the GP model optimized by the DE to predict, the accuracy is the highest, reaching 0.955, and the predicted PPV is also closest to the measured PPV. Plot the measured and predicted PPV of 16 sets of test samples obtained by various prediction methods as shown in Figure 12. The specific index evaluation effects of the forecast results are shown in Table 7.Figure 12
Relationship between measured and predicted PPV by various prediction methods.Table 7
Evaluation effect of various prediction methods.
NameR2RMSEMSEMAEMAPEVAFSDRankUSBM0.620.260.070.2115.2760.490.535LSSVM0.830.190.040.1611.3780.520.463ANN0.720.230.050.1713.4572.190.484GP0.890.160.020.118.0786.930.462DE-GP0.960.10.010.086.4595.440.421Generally, we takeδ=yi′−yi/yi×100% as the judgment of the relative error of the test samples, which can better reflect the credibility of the prediction. The relative error and calculation time of the test samples using the LSSVM, ANN, GP, and DE-GP models and USBM empirical formula are shown in Figure 13.Figure 13
The relative error and calculation time of different PPV prediction methods.It can be seen from Figure13 that the relative error range of the GP model prediction results is −15.81%∼24.29%, and the average relative error is 2.62%, which is significantly smaller than the prediction results of the LSSVM and ANN models and USBM empirical formula. The regression accuracy of the GP model optimized by the DE is significantly improved, and the relative error range of the prediction results is only −17.78%∼11.05%. It shows that the DE-GP hybrid prediction model established in this research can be effectively applied to the prediction of PPV in the construction field of rock foundation pit.In addition, the computation cost of different models in this research refers to the time spent under the conditions of six-thread CPU processor, 8 GB running memory, and RTX2060 graphics card. It can be seen from Figure13 that the computation time cost of different models is different. Compared with the GP model, the computation time cost of the DE-GP model has increased by 1′34″. However, from the perspective of the improvement of calculation accuracy, this part of the increase in computation time cost is completely acceptable.


### 5.4. Sensitivity Analysis

In order to use the PPV to measure BIGV, the most important step is to determine the sensitivity of the eight blasting influence factors. The cosine amplitude method (CAM) is a commonly used sensitivity test method [38]. Among them, the influence of each blasting influence factor on PPV can be calculated by the correlation coefficient rij, that is, the degree of influence of each blasting influence factor on PPV can be obtained from the sensitivity analysis:(8)rij=∑n=1mXinXjn∑n=1mXin2∑n=1mXjn2,where rij is the degree of each blasting influence factor; Xin is the i-th blasting influence factor of the n-th sample data; Xjn is the blasting PPV of the n-th sample data; and m is the number of samples.The value of the correlation coefficientrij is between 0 and 1. The closer the correlation coefficient is to 1, the greater the influence of input parameters on PPV. Similarly, the smaller the influence of input parameters on output parameter, the closer the correlation coefficient is to 0.The sensitivity of blasting input parameters are analyzed based on the blasting sample set data and the established LSSVM, ANN, GP, and DE-GP models, as shown in Figure14.Figure 14
PPV sensitivity analysis of different prediction methods.It can be seen from Figure14 that the blasting influence factors of the LSSVM, ANN, GP, and DE-GP models have different effects on the blasting PPV. However, the three most influential factors are distance from the blast-face, maximum charge per delay, and total charge, all of which have a degree of influence above 0.85. The less influential blasting factors are stemming length and powder factor, and the correlation coefficient rij is between 0.62 and 0.78. Therefore, the sensitivity analysis result can be used to optimize the blasting design parameters, which can be used to guide project construction.


## 5.1. Influence of GP Parameters on Regression Accuracy

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


## 5.2. Influence of DE Parameters on Optimization Effect

When using the DE-GP model, there are many factors that affect the optimization effect of the DE, and they are directly related to the selection of hyperparameters in the GP model.The control parameters in the DE mainly include mutation factorF, crossover factor CR, and population size NP. Usually, NP is 5∼10D (D stands for dimensionality), F is 0.5∼0.9, and CR is 0.5∼0.9. It can not only ensure a higher success rate of optimization but also increase the speed of convergence. The mutation factor F and the cross factor CR have a greater impact on the speed and stability of the search. Improper selection of parameters may lead to evolutionary stagnation and premature maturity. The mutation factor F mainly affects the search steps. Increasing F can increase the search range of the DE and increase the diversity of the population but, at the same time, weaken the convergence speed of the algorithm; reducing F can increase the development ability of the DE and improve the convergence speed of the DE, and also it is easy to make the population converge prematurely. The cross factor CR affects the weight of evolutionary information. Increasing CR can improve the diversity of the population and speed up the convergence speed of the DE; reducing CR is beneficial to the analysis of individual separable problems. Select different variation factor F and cross factor CR to perform optimization search trial calculation, and the corresponding iterative convergence curve is shown in Figure 9.Figure 9
Iterative curves of different variation factors and crossover factors. (a) When CR = 0.7, iterative curves of different variation factorsF. (b) When F = 0.8, iterative curves of different crossover factors CR.
(a)(b)It can be seen from Figure9 that when CR = 0.7 and F ∈ [0.5, 0.9], DE iterative search can converge, but the iteration speed is different. When F = 0.8, the convergence speed is the fastest, the number of iteration steps is less, and the iteration steps reach 54 steps to complete convergence. When F = 0.6, the convergence speed is the slowest, and the iterative steps reach 85 steps before it converges completely. Therefore, choosing appropriate search parameters can greatly save 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 reaches the fastest convergence, and it is completely converged in 39 iterations. When CR = 0.6, the convergence 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 comparisons should be made within the selection range of possible convergence parameters in order to obtain the optimal convergence result.When the variation factorF = 0.8 and the crossover factor CR = 0.9 are fixed, increasing NP can improve the diversity of the population but, at the same time, reduce the convergence speed of the population. Reducing NP can increase the convergence speed, but it is easy to lead to premature convergence. When the population number NP is equal to 40, 60, 80, 100, and 120, respectively, the iterative convergence curves are shown in Figure 10.Figure 10
Iterative convergence curves of differentNP.It can be seen from Figure10 that when NP is 40 and 60, the optimal value is searched for 54 iterations and 48 iterations, respectively. But when NP is 80, 100, and 120, it takes 57, 73, and 85 iterations to search for the optimal value. NP mainly reflects the amount of population information in DE. The larger the population size, the richer the population information contained, it is more likely to find the global solution, and it is less likely for individual optimal solutions to dominate the overall solution. However, the calculation time is relatively long, and the calculation convergence is relatively slow, which lead to an increase in the calculation cost of a single population. On the contrary, the diversity of the population is restricted, which is not conducive to DE obtaining the global optimal solution and may even cause the search to stagnate. It can also be seen from Figure 10 that the fitness convergence speed does not decrease monotonically with the increase of NP. When NP is set to 60, convergence can be achieved at the fastest speed. So, NP is set to 60.Choosing appropriate initial parameters can not only save time and improve convergence speed but also avoid the algorithm from falling into local optimality. If the selection of parameters is not appropriate, it may cause the algorithm search to stagnate due to excessive emphasis on the search ability or cause the algorithm to converge prematurely due to excessive emphasis on the development ability. In this research, the final DE algorithm control parameters settings:NP = 60, F = 0.8, and CR = 0.9.


## 5.3. Computing Ability Evaluation of Different Prediction Methods

Since the constructed blasting PPV model will directly affect the accuracy of the prediction results, it is necessary to establish the evaluation index of the model to test the prediction effect. From a statistical point of view, it is more one sided to comprehensively evaluate a predictive model with only one performance index. Therefore, in this research, statistical indicators such as coefficient of determination (R2), root mean square error (RMSE), mean square error (MSE), mean absolute error (MAE), mean absolute percentage error (MAPE), variance account for (VAF), and standard deviation (SD) are used to evaluate the performance of the model, as shown in Table 5.Table 5
Statistical indicators of blasting PPV prediction effect.
Statistical indicatorsFormulasDesired valueCoefficient of determinationR2=1−∑i=1myi−yi′2/∑i=1myi−yi′21Root mean square errorRMSE=1/m∑i=1myi−yi′20Mean square errorMSE=1/m∑i=1myi−yi′20Mean absolute errorMAE=1/m∑i=1myi−yi′0Mean absolute percentage errorMAPE=1/m∑i=1myi−yi′/yi×100%0Variance account forVAF=1−varyi−yi′/varyi×100%100Standard deviationSD=1/m∑i=1myi−1/my1+⋯+ym20Among them,yi is the measured PPV, yi is the predicted PPV, y¯ is the average value of the measured PPV, and m is the total number of samples. It should be pointed out that the closer R2 is to 1, the closer the values of RMSE, MSE, MAE, MAPE, and SD are to 0, and the closer the VAF is to 100, indicating that the performance of the prediction models is better.According to the site constants in the traditional empirical formulas obtained above, 16 sets of test samples are used to test the predictive ability of the above several commonly used empirical formulas. Each index is shown in Table6.Table 6
Predictive effects of several empirical formulas.
NameR2RMSEMSEMAEMAPEVAFSDRankUSBM0.620.260.070.2115.2760.490.531LK0.390.350.120.2622.1533.320.694General predict0.220.410.160.326.2412.610.836AH0.610.280.080.2115.8559.720.572BIS0.370.350.120.2724.6227.510.755CMRI0.570.30.090.2317.1855.790.643It can be seen from Table6 that when using the USBM empirical formula, R2 = 0.62, RMSE = 0.26, MSE = 0.07, 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 selected among the six empirical formulas for comparison with other blasting intelligent algorithms.In this research, the performance of the hybrid DE-GP prediction model is compared with LSSVM, ANN, and GP models and USBM empirical formula. In the ANN model, the number of hidden layer nodesNf of the algorithm is 5, the number of iterations is Ep = 100, and the learning rate is Lr = 0.1. In the GP model, the variance of noise is σn2 = 3.72, the signal variance of the kernel function is σf = −1.7, and the hyperparameter for correlation measurement is l = 0.3. In the core parameters of the LSSVM model, the regular parameter is γ = 9.54, and square bandwidth is σ2 = 2.17. In the DE-GP model, the population number NP is 60, mutation factor F is 0.8, crossover factor CR is 0.9, and the upper limit of evolutionary algebra is 100. The relationship between measured PPV and predicted PPV using LSSVM, ANN, GP, and DE-GP models is shown in Figure 11.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.
(a)(b)(c)(d)It can be seen from Figure11 that using the ANN model to predict the PPV of rock foundation pits has the lowest accuracy, R2 is 0.724. While using the LSSVM model and the GP model, R2 is 0.832 and 0.885, respectively. After finally using the GP model optimized by the DE to predict, the accuracy is the highest, reaching 0.955, and the predicted PPV is also closest to the measured PPV. Plot the measured and predicted PPV of 16 sets of test samples obtained by various prediction methods as shown in Figure 12. The specific index evaluation effects of the forecast results are shown in Table 7.Figure 12
Relationship between measured and predicted PPV by various prediction methods.Table 7
Evaluation effect of various prediction methods.
NameR2RMSEMSEMAEMAPEVAFSDRankUSBM0.620.260.070.2115.2760.490.535LSSVM0.830.190.040.1611.3780.520.463ANN0.720.230.050.1713.4572.190.484GP0.890.160.020.118.0786.930.462DE-GP0.960.10.010.086.4595.440.421Generally, we takeδ=yi′−yi/yi×100% as the judgment of the relative error of the test samples, which can better reflect the credibility of the prediction. The relative error and calculation time of the test samples using the LSSVM, ANN, GP, and DE-GP models and USBM empirical formula are shown in Figure 13.Figure 13
The relative error and calculation time of different PPV prediction methods.It can be seen from Figure13 that the relative error range of the GP model prediction results is −15.81%∼24.29%, and the average relative error is 2.62%, which is significantly smaller than the prediction results of the LSSVM and ANN models and USBM empirical formula. The regression accuracy of the GP model optimized by the DE is significantly improved, and the relative error range of the prediction results is only −17.78%∼11.05%. It shows that the DE-GP hybrid prediction model established in this research can be effectively applied to the prediction of PPV in the construction field of rock foundation pit.In addition, the computation cost of different models in this research refers to the time spent under the conditions of six-thread CPU processor, 8 GB running memory, and RTX2060 graphics card. It can be seen from Figure13 that the computation time cost of different models is different. Compared with the GP model, the computation time cost of the DE-GP model has increased by 1′34″. However, from the perspective of the improvement of calculation accuracy, this part of the increase in computation time cost is completely acceptable.


## 5.4. Sensitivity Analysis

In order to use the PPV to measure BIGV, the most important step is to determine the sensitivity of the eight blasting influence factors. The cosine amplitude method (CAM) is a commonly used sensitivity test method [38]. Among them, the influence of each blasting influence factor on PPV can be calculated by the correlation coefficient rij, that is, the degree of influence of each blasting influence factor on PPV can be obtained from the sensitivity analysis:(8)rij=∑n=1mXinXjn∑n=1mXin2∑n=1mXjn2,where rij is the degree of each blasting influence factor; Xin is the i-th blasting influence factor of the n-th sample data; Xjn is the blasting PPV of the n-th sample data; and m is the number of samples.The value of the correlation coefficientrij is between 0 and 1. The closer the correlation coefficient is to 1, the greater the influence of input parameters on PPV. Similarly, the smaller the influence of input parameters on output parameter, the closer the correlation coefficient is to 0.The sensitivity of blasting input parameters are analyzed based on the blasting sample set data and the established LSSVM, ANN, GP, and DE-GP models, as shown in Figure14.Figure 14
PPV sensitivity analysis of different prediction methods.It can be seen from Figure14 that the blasting influence factors of the LSSVM, ANN, GP, and DE-GP models have different effects on the blasting PPV. However, the three most influential factors are distance from the blast-face, maximum charge per delay, and total charge, all of which have a degree of influence above 0.85. The less influential blasting factors are stemming length and powder factor, and the correlation coefficient rij is between 0.62 and 0.78. Therefore, the sensitivity analysis result can be used to optimize the blasting design parameters, which can be used to guide project construction.


## 6. Conclusions

The present research established a PPV prediction model based on the DE-GP intelligent algorithm. It was compared with empirical formulas and other intelligent algorithms, taking the blasting of the foundation pit at the Suoyuwan South Station as an example. The following conclusions are obtained:(1)
The conventional empirical formulas are mainly related to distance from the blast-face and the maximum charge per delay. However, the calculated PPV deviates greatly from the measured PPV.(2)
The blasting PPV prediction model based on the DE-GP method was established for the first time, and it was found to have high accuracy by comparing the measured PPV and the predicted PPV.(3)
The DE-GP blasting prediction model used for the rock blasting prediction. Compared with the LSSVM, ANN, and GP prediction models and empirical formulas, the prediction accuracy was improved, and the relative error was significantly reduced.(4)
The CAM was used to analyze the sensitivity of blasting influence factors for different prediction models. The impact of the distance from the blast-face, maximum charge per delay, and total charge are the largest; stemming length and powder factor are relatively small.The PPV prediction model based on DE-GP can be effectively applied to the blasting prediction in geotechnical engineering. It can control the adverse effects of blasting and reduce the harm of BIGV to the surrounding environment. This method can provide strong technical support for the blasting construction of similar engineering.



---

*Source: 2847180-2021-10-04.xml*