#### Abstract

The fruit fly optimization algorithm-general regression neural network (FOA-GRNN) coupled model and the Finite Element Method-Smoothed Particle Hydrodynamics (FEM-SPH) numerical calculation method are comprehensively used. The control problem of blasting vibration in the process of mining hidden resources under complex environmental conditions has been studied. Taking a lead-zinc mine as the engineering background, the development of hidden resources in the collapse area due to unreasonable mining was studied. Based on the establishment of the first mining stope and its mining method in this area, biosimulation and generalized neural networks were introduced to solve this problem, the coupling of blasting parameters was analyzed, and the 3D nonlinear dynamic coupling model was constructed for numerical simulation. The results show that the blasting parameters of deep-hole mining were optimized, including the values of six output quantities: hole distance, row spacing, side hole distance, explosive unit consumption, minimum resistance line, and interval ratio (the Root Mean Squared Error value is only 0.051). The error between the network optimization parameters and the empirically obtained values was controlled to within 0.05; five possible edge-hole charge structures were designed (the interval ratio is 0.696), and the vibration velocity peak and pressure peak variations with time after detonation are reflected by the simulation results. The dynamic evolution law of the rock mass velocity vector and the damage of the rock damage was revealed. According to the analysis in this paper, the smallest and optimal edge-hole charge structure of the surrounding rock was obtained.

#### 1. Introduction

The underground goaf is common in mining projects [1–3]. They are temporary or permanent underground spaces formed for ore mining [4]. The roof and the two gangs of the goaf are susceptible to pressure changes and redistribution caused by the sudden collapse of the overburden. Unstable underground goaf may damage the stability of the surrounding diverticulum and cause instability, which not only causes huge economic losses but also causes serious harm to personnel [5–7]. Compared with the instability damage caused by excavating intact rock masses, the goaf of the gangs which are filled with complex materials is more likely to collapse. In fact, although the engineer will have the mining sequence in the early stage of the excavation, in the mining environment, if the pillar is not enough to balance the external force and prevent the loosening and deformation of the two-package, the roof is likely to occur and collapse. In addition, the control of the stability of the two-pack backfill plays an important role in the mining; especially for the larger underground goaf, improper mining design will directly lead to the collapse of the roof. Generally, the influence of blasting vibration on the stability of the surrounding rock in the collapsed area is very important. The study on the stability of surrounding rock under blasting impact in the collapsed area is a popular and difficult problem in ore body mining.

When mining in the collapse zone, the surrounding rock mechanical properties are poor, the strength is low, and the disturbance caused by the excavation will be very obvious. The blasting dynamic impact load may generate more pressure and more settlement problems. Even a slightly inappropriate excavation scheme can lead to a larger area of collapse. This will seriously affect the safety of mining. Therefore, the careful comparison of mining methods and proper selection of blasting schemes is particularly important for exploiting hidden resources under such conditions. There are still many unanswered questions about how to efficiently and safely exploit hidden resources under unfavorable mining conditions. The key to dealing with such problems is to select reasonable blasting parameters to reduce the impact damage to surrounding rocks. At this point, it is necessary to conduct an in-depth study of the blasting damping scheme.

There are many methods for analyzing the stability of surrounding rock under blasting impact under complex environmental conditions. In summary, there are three classic methods: (1) the physical model experimental method [8], (2) the mechanical model analysis method [9], and (3) the numerical simulation method [10, 11]. These methods have some limitations according to their respective characteristics. For example, using physical model experimental methods, due to a large number of different rock mass materials in the collapsed area, it is difficult to construct physical models consistent with the actual conditions, and this work consumes many manhours. Additionally, the material resources and the impact of the loading and measurement of blasting power are difficult to determine. In addition, mechanical analysis methods are solutions that usually rely on certain assumptions. The surrounding rock faces complex hydrogeological and engineering geological conditions in the collapsed area, so this analysis method is also very limited [12–15].

With the rapid development of computer technology, numerical simulation has been used in recent years to study the stability of surrounding rock in mining. Several models, such as finite element models (FEMs) [16], boundary element models (BEMs) [17], and discrete element models (DEMs) [18], have been proposed for numerical simulation, and all three methods have advantages and disadvantages. The Finite Element Method is the most popular; it has been widely used to analyze practical problems in geotechnical engineering.

In general, numerical simulation methods for studying the stability of surrounding rock under blasting impact loads in a collapsed area are very popular. Considering the fact that the exploitation of hidden resources in a special environment requires higher control of the blasting vibration [19], the coupling method Smoothed Particle Hydrodynamics- (SPH-) DAM-Finite Element Method (FEM) reproduces the overall blasting responses of presplitting blasting and step blasting. The accuracy and feasibility of this coupling technique are verified by on-site monitoring. The results show that this method has an obvious control effect on blasting vibration. Based on this method, SPH-FEM is coupled in AUTODYNA numerical simulation software [20] and can be used to simulate and explore mining blasting vibration control under complex environmental conditions. Due to a large number of mining sites in previous mining efforts, a range of blasting parameters can be selected. When mining the hidden resources, it is necessary to learn from the previous blasting parameters and then optimize the mining sites for the collapsed area. Reference [21] proposed an ISVR model for mapping the parameters of smooth blasting parameters. By coupling a genetic algorithm to search the optimal parameters of the ISVR model, a control optimization method for smooth blasting parameters is proposed, and all the blasting parameters are determined. Preferably, this process takes a long time, which is not practical. Therefore, the appropriate preferred method is selected to optimize the preliminary blasting parameters, and then the numerical simulation is carried out so that the blasting parameters can be accurately designed.

In the research presented here, a nonlinear three-dimensional numerical model is established to simulate the dynamic process of stope excavation under complex environmental conditions. The collapsed area of a lead-zinc mine was selected as the survey object. First, the first mining stope in the collapsed area is established, and then the blasting parameters in the collapsed area are optimized. Then, the Finite Element Method is used to analyze the deformation and stress distribution characteristics of the surrounding rock during the mining process in the collapsed area. The stability of surrounding rock under blasting impact load in the collapsed area is studied. The research results reveal the mechanism of the influence of the surrounding rock under the blasting impact load on the stability of the collapsed area. In the present study, the above five blasting damping schemes compare the dynamic evolutions of the stress and deformation, the displacement field of the surrounding rock mass, and the distribution characteristics of the damaged area. Finally, a suitable blasting charge structure scheme is proposed for the recovery of hidden resources in the collapsed area. The research results can provide necessary guidance for the safe and efficient recovery of hidden resources with surrounding rock control in the collapsed area. In addition, the results can provide a reference for the optimization of the blasting parameters of stopes in similar mine collapse areas.

#### 2. Hidden Resources Mining Status

The bottom of the No. 1 stope in a metal mine is approximately 22.5 m high and collapses after mining stops (Figure 1), which affects the safety of the adjacent six stopes. The No. 7 stope is a relatively safe stope. Safe mining of the stope and adequate and high-quality filling work is the key to recovering more mining resources. Because the No. 7 stope is very high (55 m), the mining sequence is from bottom to top, and for the bottom 22 m high ore body (out of the virtual coil in Figure 1), the mining method is bottomless deep-hole back mining. After filling, the blasting vibration should be reduced during mining to protect the surrounding rocks as much as possible.

#### 3. Optimization of Blasting Parameters Based on FOA-GRNN Coupling

Determining reasonable blasting parameters is a prerequisite for effectively controlling the blasting vibration of the abovementioned stope. The traditional empirical formula has a large range of blasting parameters and is subjective, while the field test is expensive and the experimental location is difficult to select. By fully considering the nonlinear relationship between the blasting parameters and their influencing factors [22], the neural network was introduced, and its high degree of nonlinear mapping ability can avoid the establishment of mathematical equations in the optimization process of blasting parameters, allowing the adaptation of the nonlinear relationship between the blasting parameters and their mapping parameters. A fruit fly optimization algorithm (FOA) that simplifies the optimization process and improves engineering applicability was introduced. With the fruit fly optimization algorithm (FOA) and optimization of the neural network [23], better results could be achieved. The generalized regression neural network (GRNN) is a neural network model based on nonlinear regression. It has advantages over traditional neural networks in terms of data fitting and learning speed, especially when the training samples are small, and it has achieved good regression effects. It has been successfully applied in research in fields of study such as rock burst.

##### 3.1. GRNN and Its Standard of Error Analysis

The above General Regression Neural Network (GRNN) model is composed of four layers: input layer (input layer), mode layer (pattern layer), summation layer (summation layer), and output layer (output layer) (Figure 2). The number of neurons in the input layer is equal to the dimension of the input vector in the sample. The input vector was directly submitted to the pattern layer with the number of *n* elements. Each element of the layer corresponds to a training sample, and the Gaussian function was generally selected as the kernel function. *X*_{i} is the center vector of each unit kernel function. The summation layer contains molecular units and denominators; the output value of each unit is in the denominator unit summation mode layer, and the molecular unit weights each unit in the pattern layer, where the weight of each training sample is the value of each training sample; the output layer divides the output of the summation numerator and the denominator to obtain an estimate of .

In the joint density function *f*(*x*, *y*) of two random vectors, let *x*_{0} denote the observation of *x*. The value of the function *y* based on *x*_{0} is *y*_{0}, and its mathematical description is

Assume *f* (*x*, *y*), with a subjective and normal distribution, and then estimate the density function.

Here, and ( = 1, 2, 3, …, ) are the input and output values of the *i* sample, *d* is the dimension of the random vector *x*, and is the smoothing factor, which is also known as the spread parameter value.

Substituting with in (1),

The General Regression Neural Network (GRNN) model is simpler than the backpropagation(BP) neural network model, and only one spread parameter value is needed, which reduces the influence of human factors and reduces the randomness of the network structure design. The fruit fly optimization algorithm (FOA) was used to optimize the GRNN model, and the optimal smoothing factor was selected to minimize the mean square error between the output value and the actual value of the GRNN model [24]. The Root Mean Squared Error (RMSE) was used as the analytical standard, that is, the difference between the parameter analysis value and the true value of the parameter [25]. RMSE could be used to measure the average error while analyzing the degree of change in the data. The smaller the RMSE value was, the higher the concentration of preserved fruit flies [26], and the more accurate the analytical model.

Here, is the original value and is the analytical value.

For the GRNN model that was based on the FOA optimization algorithm, the mean square error of the output value of the training sample and the actual value was used as a moderate function to find the optimal smoothing factor. The specific steps were as follows [27].

Set the size of the population to N, which is the population of randomly formed initial fruit flies [28]:

Here, is an evolutionary individual encoded by real numbers:(i)Based on the FOA optimization algorithm, the optimization process is shown in Figure 3. It was programmed by MATLAB, and the parameter values in GRNN were minimized so that the individual flies and the fruit fly population were optimal and consistent.(ii)Substituting the optimal solution into the parameters of the GRNN, the MATLAB function newgrnn() was used to input the training data and obtain the network output value.

##### 3.2. Establishment of an FOA-GRNN Blasting Parameter Optimization Model

The input layer included ore bulk density, elastic modulus, tensile strength, rock solidity, internal friction angle, and cohesion. The output layer includes hole pitch, row spacing, edge-hole distance, explosive unit consumption, minimum resistance line, and interval ratio. In the early stage, the measured model of the stope detected by 3D Cavity Monitoring System (CMS) was sorted out, and 105 sets of blasting parameter data with good blasting effect were obtained, which were divided into two parts:(i)One hundred parameter datasets were used as training samples.(ii)20 datasets were used as testing samples. The MATLAB software was used to optimize the training of the fruit fly optimization algorithm-general regression neural network (FOA-GRNN) model sample. For the relationship between the lengths, Table 1 lists 100 sets of training sample data. In this table, *ρ*_{r} is the density of rock (g/cm^{3}), *E* is the elastic modulus (GPa), *R*_{m} is the compressive strength (MPa), *F* is the rock solidity, *φ* is the internal friction angle (°), *c* is the cohesion (MPa), *H*_{d} is the hole distance (m), *d*_{b} is the distance between two rows (m), *S*_{h} is the side hole spacing (m), *E*_{u} is the explosive unit consumption (kg/m^{3}), *W* is the minimum resistance line (m), and *I*_{r} is the interval ratio.

Among them, the FOA program is needed to dynamically adjust the smoothing factor in GRNN and construct an analytical model. To test the difference in model classification detection capabilities, the value needs to be normalized to be in the range [0∼1]. Therefore, the input parameters are normalized according to the following formula to prevent input variables with different physical meanings and dimensions from causing unequal use. At the same time, this normalization can speed up the convergence of the network.

In the formula, is the actual input of the *i* input indicator at the *i* sample point, and *X*_{max}(*i*) and *X*_{min}(*i*) are the actual maximum and minimum values of the *i* input indicator.

Because the network output is a preferred result after normalization, the network output needs to be denormalized.

##### 3.3. FOA-GRNN of the Optimization Blasting Parameter

We assign P1 to P6 in the FOA-GRNN model as input neurons and E1 to E6 as output neurons. The position of the fruit fly population is initialized to [0, 1]; the iterative fruit fly searched for the random flight direction of the food, and the range of the distance was [-10, 10]; the number of iterations is 100 (shown in Figure 4). After searching for the best P1 to P6 through 100 iterations, the analysis results approached the RMSE convergence trend graph between the analytical value and the actual target value. Figure 4 shows the trajectory of the Drosophila group used in the optimization; the group of fruit flies was concentrated and flew to the location near the highest concentration of flavor [4.4, 3.8] (shown in Figure 5). The results suggest that the optimization process was orderly and converged in a uniform direction, so the results obtained were accurate and reliable. The RMSE value was only 0.051, and it is considered that the parameter preference model was stable and effective.

Five sets of test samples were used to determine the blasting parameters using the FOA-GRNN model, and the test sample data were input (Table 2). The comparison between the FOA-GRNN preferred results and the actual engineering blasting parameters is shown in Table 3 (E∗1 to E∗6 were the output neurons of the test samples).

From the comparison between the FOA-GRNN parameter optimization results and the actual blasting parameters of the mine, it was found that the error of the network optimization parameters and the empirically obtained values were controlled to within 0.05, and the preferred results were highly reliable. On this basis, the mining section parameters were input into the trained FOA-GRNN model, and the following optimized blasting parameters were obtained, as shown in Table 4.

#### 4. Numerical Model of Three-Dimensional Coupled Fluid-Solid Nonlinear Dynamics

##### 4.1. Principle of the Coupled Fluid-Solid Model

Practice proved that optimizing the edge-hole charge structure could reduce the blasting vibration, so in addition to the preferred blasting parameters, the charge structure of the side hole should be optimized. A numerical model to reduce the occurrence of earthquakes in blasting was constructed by using AUTODYN-3D, in which the explosive and stope units used fewer smooth-mesh Smoothed Particle Hydrodynamics (SPH) particles. The coupled Finite Element Method (FEM) and SPH algorithm was used to constrain the boundary. In addition, the SPH particle and the LAGRANGE unit directly adopt point-point contraction, and the absorption boundary constraint was selected for the side and bottom of the LAGRANGE unit. Figure 6 is a coupling technique flow chart [29].

The JWL equation of state could be used to indicate the relationship between the specific volume (*υ*) and pressure (*P*_{0}) produced by explosives throughout the detonation process [30]:

Here, *P*_{0}, *V*_{x}, and *E*_{c} are the pressure, relative volume, and initial specific energy, and *A*_{c}, *B*_{c}, *R*_{1}, *R*_{2,} and *ω* are the material constants. Determining these correlation coefficients yields a state equation for the explosive. From the CJ condition, the detonation pressure *P*_{b}, detonation rate *P*_{b} and chemical energy *E*_{b} of a given explosive could be determined.

The ideal gas equation of state is used for the air

Here, *P*_{0} is the pressure, *e*_{0} is the specific heat, *ζ* is the polytropic exponent, and *P*_{s} is the pressure offset. Here, *ζ* = 1.4, *ρ* is 1.23 mg/cm^{3}, and *e*_{0} is 2.07 × 105 mm/mg/ms.

The Riedel-Hiermaier-Thoma model (RHT) constitutive parameters of rock were determined by laboratory tests and solved according to the determination method of RHT material constitutive model parameters [31], as shown in Table 5 and Table 6. In Table 7 [32], *A*_{1}, *A*_{2}, *A*_{3}, *B*_{0}, *B*_{1}, *T*_{1}, *T*_{2}, and *T*_{0} are the material parameters; *ρ*s_{0} is the density of the material; *F*_{t} is the tensile strength; *F*_{s} is the shear strength; *F*_{c} is determined by *F*_{s}; *A*, *N*, and *Q* are the failure strength parameters; *α* and *δ* are the tensile strain rate index; uniaxial compression (TENSRAT), the ratio of the elastic limit to the uniaxial compressive strength, and COMPRAT, the ratio of the uniaxial tensile elastic limit to the uniaxial tensile strength, are the parameters on the elastic limit surface; *B* is the residual strength surface constant; *M* is the residual strength surface index; PREFACT is the ratio of the elastic shear modulus to the plastic hardening shear modulus; SHRATD is the ratio of the residual shear modulus to the shear modulus before the damage.

##### 4.2. Three-Dimensional Numerical Dynamic Calculation Model for Coupled Flow-Solid to Explore the Hidden Danger in the Stope

The mining process is shown in Figure 7. The specific steps are as follows. First, a track-free track roadway is formed. A rock drilling chamber approach and a top rock drilling chamber are also formed. The first step blasting forms a bottom space, and then blasting upward occurs until the roof is broken, at which point lateral blasting occurs. The ore is transported from the lower chamber. Figure 8 shows the five feasible edge-hole charging structures designed with the same charge (0.5 rolls or 1 roll, while the air gap was 0.5 m, 0.6 m, 1.0 m, 1.1 m, or 1.2 m). Each roll of explosive has a length of 0.67 m and a diameter of 90 mm, with an interval ratio of 0.696. The numerical analysis model was established by simplifying Figure 7, as shown in Figure 9. The distance from the side hole to the filling body was 1.2 m, the minimum resistance line was 1.5 m, the width was 6.6 m; the model length was 20 m, the width was 24 m, the height was 20 m, and the boundary was set as the transmission boundary. Figure 10 shows the arrangement of the Gaussian point monitoring.

#### 5. Optimization of Nonlinear Dynamic Disturbance

In this study, the original nonlinear dynamic disturbance is optimized by setting Gaussian point elements in the blasting process of five schemes on the original basis and then reading the peak values obtained. Two directions, horizontal and vertical, are set, respectively, and the most effective charging structure scheme in lateral state protection is analyzed according to the knot of blast vibration velocity and pressure peak.

##### 5.1. Law of Movement and Velocity Vector Variation in the Rock Mass

It can be seen from Figure 11 that after the explosives in the blast hole are detonated in turn, the cracks that formed in the rock mass first appeared at the detonation point of the first row of blast holes, and then, the propagation of the detonation wave appeared to be ellipsoidal, rapidly expanding to the periphery, with a peak value at the focus of the ellipsoid [19]. When the detonation wave decayed, the velocity vector appeared elliptical and expanded around the blast hole. Because the products of detonation in the holes were mutually squeezed, there are more significant speed vector concentration phenomena in the two rows of the first row of blast holes [33]. The second row of blast holes detonates and forms a new velocity field before the first row of blast hole velocity fields dissipate [34], superimposing with the previous velocity field. As the stress waves continued to diffuse and decay, the velocity field becomes more uniform. When the explosive was just detonated [35], the maximum deformation of the rock mass was mainly concentrated on the wall of the hole subjected to the detonation pressure. Under the combined action of the detonation wave and the explosion stress wave, the wall of the blast hole was damaged by the power disturbance [36]. When the magnitude of the tensile stress was greater than the ultimate strength of the rock, the ore near the free surface was thrown out. When the crack penetrated, the particles between the first row of blast holes and the free surface were destroyed and thrown out [37].

##### 5.2. Dynamic Evolution Law of Surrounding Rock Blasting Damage

Figure 12 shows that after the explosive was detonated, the high-temperature and high-pressure detonation products pushed the surrounding rock mass to expand outward, the violent shock wave pressure directly acted on the blast hole, and the wall of the blast hole in the blast zone first broke. When the pressure around the wall of the blast hole was greater than the compressive strength of the two rams, it was subjected to compression yielding until crushing, thereby producing a compression comminution zone. Over time, the shock wave propagated in the two-package body, causing radial cracks in the tangential tensile stress zone; the cracks expanded toward the periphery of the sample. When the shock wave reached the boundary of the two filling bodies, the wave impedance of the two filling bodies was different from that of the surrounding air so that the shock wave was reflected in the two filling bodies, thereby forming a tensile wave and accelerating the destruction of the two filling bodies. Because of the tensile stress caused by the lateral deformation of the two-package body, the part was damaged earlier. Under the dual action of tensile stress and cracking, the initial crack that formed expanded to the periphery, resulting in large-scale fracture of the two-layer filling body. The damaged area was also expanding. Subsequently, the second row of blast holes detonated and interacted with the ore blocks that had collapsed due to detonation from the previous row of blast holes. When the blast wave was terminated, the internal fluctuations of the two filling bodies would not stop under the action of inertia, and finally, the two gang filler bodies lost part of their bearing capacity.

##### 5.3. Vibration Velocity and Pressure-Time Curve Analysis of Surrounding Rock Blasting

The blasting speed-time and pressure-time curves in Figures 13 and 14 show that due to the reduction in the charging of the northern gang, the peak blasting vibration and the blasting pressure are smaller than those of the south gang for the two main blast holes. The value of the corresponding position indicates that the blasting vibration could be reduced to some extent by reducing the charge of the side hole.

By extracting the Gaussian point elements on the northern gang during the blasting process of the five schemes, the peaks obtained by reading are read. Through the comparative analysis of the blasting vibration velocity and the pressure peak in the horizontal direction and the vertical direction, it can be concluded that the charging structure scheme that is the most effective protection for the side states is exhibited in the following order: scheme 5, scheme 3, scheme 4, scheme 2, and scheme 1, shown in Table 8.

#### 6. Conclusions

In this study, the FOA-GRNN coupling model was used to optimize the blasting parameters of deep-hole mining. On this basis, a three-dimensional numerical calculation model of the nonlinear display dynamics of the first mining stope was constructed, the velocity vector of the ore rock and the dynamic change law of the rock damage were revealed, and the optimal edge-hole charge structure that creates the least damage to the surrounding rock was determined. The conclusions are as follows:(1)The error between the optimized parameters of the FOA-GRNN model network and the empirically obtained values is controlled within 0.05, and the preferred results are highly reliable.(2)The more uniform the explosive is distributed upward in the gun aperture, the smaller the blasting vibration effect is. The order of the effectiveness of the five kinds of side hole charging structures to control the blasting vibration effect is as follows: scheme 5 > scheme 3 > scheme 4 > scheme 2 > scheme 1, suggesting that the charging structure of the side hole adopts the half-volume medicine roll plus 0.6 m bamboo tube charge structure.

The FOA-GRNN model has certain limitations in the optimization and 3D numerical simulation of blasting parameters and does not consider the influence of objective factors such as deposit thickness, dip angle, geological structure, mining depth, and weak interlayer. And in this study, the fuzzy and nonlinear characteristics of the above influencing factors are fully considered in the one-step research, which expands the sample information base, improves the optimization model, and increases the optimization accuracy, which can give some scientific support to mining enterprises and related scientific research work.

#### Data Availability

All the data included in this study are available upon request to the corresponding author.

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Acknowledgments

This research was funded by the Hunan Province Science Foundation, under grant number “2021JJ30679,” and Hunan Provincial Department of Education General Project, under grant number “19C1744.”