In this section, we first verify the performance of proposed method through synthetic tests and then apply our proposed method to real data.
3.1. Synthetic Tests
In this section, we first build a simple tunnel model and then test the accuracy and efficiency of the HVM-based approach with the synthetic microseismic data. Finally, the noise immunity of the HVM-based approach is analyzed by adding the arrival time error, and the SVM-based approach is compared to the HVM-based approach.
3.1.1. Establishment of the Tunnel Engineering Model
As shown in
Figure 2, the microseismic monitoring area range was a cube, along x coordinates from 0 m to 200 m, y coordinates from −30 m to 30 m, and z coordinates from −30 m to 30 m. A rectangular section of the tunnel, perpendicular to the x direction, was 5 m × 5 m. A total of six receivers were arranged on both sides of the tunnel, which are represented by the green triangles and denoted as R1, R2,..., R6. The three designed seismic sources, represented using red circles, are denoted S1, S2, and S3. The specific parameters of the receivers and seismic sources are shown in
Table 1. The wave velocity in the monitoring area was 5000 m/s. However, since the medium in the tunnel after excavation was air, and the energy of the waves in air is greatly attenuated, the propagation velocity in the air was much lower than the propagation velocity in the rock. Therefore, the velocity of the empty area after the tunnel excavation was set to 340 m/s.
The propagation of waves follow the principle of minimum travel time. As the wave travels slowly in the tunnel after excavation, it bypasses the empty area and travels in the rock. In this example, the rock is homogeneous, and the wave velocity is the same throughout the rock; therefore, the minimum path from the source to the receiver is the shortest path. According to the knowledge of spatial analytic geometry, we can obtain the minimum travel path, namely, the theoretical path, which is represented by a solid black line in
Figure 3.
The specific calculation method of the theoretical path is as follows.
Figure 4a is the tunnel model. The faces ADHE and BFGC are the two sides of the tunnel. R is the receiver that is very close to the face ADHE, and S is the source near the face BFGC. R′ and S′ are the vertical projections of R and S on the faces ADHE and BFGC, respectively. We first find the shortest distance from R′ to S′. We expand the tunnel model vertically to obtain
Figure 4b. According to the principle that the line segment between two points is the shortest, we obtain the shortest path from S′ to R′, which intersects the line segments AD and BC at points M and N. We then utilize plane DNMR to obtain the shortest path from S to R (green line in
Figure 4a). According to the theoretical path, the theoretical arrival time of each source to each receiver can be calculated, as shown in
Table 2.
3.1.2. MEL Based on an HVM
Based on the abovementioned tunnel engineering model, we used an HVM to locate the three sources. First, we meshed the monitoring area with a grid size of 0.5 m × 0.5 m × 0.5 m, thus obtaining 400 × 120 × 120 unit blocks. The tunnel geometry model was used as a constraint to assign velocity values to the unit block. The velocity of the unit block inside the tunnel was 340 m/s, and that of the unit block outside the tunnel was 5000 m/s.
BL was carried out using the theoretical arrival time in
Table 2 as input parameters. Each block was a potential source, and the residual arrival time of each block was calculated according to Equation (5). The theoretical travel time for the waves between the sources and receivers was calculated using the FMM, as shown in
Figure 4. The red lines in the figure represent the calculated path from the sources to the receivers. We used the grid search method to assign the block with the minimum value of Equation (5) as the target block, and the center coordinates of the target block represent the BL results, as shown in
Table 3.
AL was then performed in the target block according to Equations (6)–(9). The minimum value of Equation (9) was found using PSO, and the AL results are shown in
Table 3.
Figure 5 shows that the AL results of the three sources were very close to the theoretical position. The PSO iteration parameters were as follows: The maximum number of iterations was 2000, and the number of seeds was 80. The acceleration parameters of the algorithm were 2 and 2, which affected the local and global optimal values, respectively. The weighted values for the initial and convergence moments were 0.9 and 0.4, respectively. The threshold of the termination algorithm was 1 × 10
−25, and when the minimum value of the target function was less than this value, the algorithm stopped. The change in the value of the objective function in the PSO iteration is shown in
Figure 6, and the three sources converged after 25 iterations.
A description of the computational efficiency is as follows. First, we used the FMM to calculate the travel time of the waves from each receiver to all grids and saved these travel times in a database. As long as the receiver position and velocity model were not changed, each subsequent microseismic event was located via the travel time database, such that the FMM solution needed to be calculated only once. In this paper, the FMM code was based on C++ programming, and the other code was based on MATLAB programming. In this case, it took 51 s to construct the travel time database using the FMM. The BL and AL computation times for the three events are shown in
Table 3. After obtaining the travel time database, locating an event took approximately 17 s. All the above programs ran on a 3.6 GHz Intel Core i9-9900k CPU.
As shown in
Figure 7, we compared the BL results of the three sources with the AL results. The teal columns in
Figure 7 represent the errors of the BL results, and the pink columns represent the errors of the AL results. The location results include the X, Y, and Z coordinate errors, spatial error, origin time error, and minimum value of the target function. The AL errors of S1 and S2 were larger than the BL errors in the Y direction, and the AL errors of S3 were larger than the BL errors in the
Z direction. We calculated the spatial location error using Equation (10):
where
represents the spatial location error. , , and represent the location errors of the X, Y, and Z coordinates, respectively. , , and represent the location results of the X, Y, and Z coordinates, respectively. , , and represent the theoretical values of the X, Y, and Z coordinates, respectively.
However, the AL results of the three sources were clearly smaller than the BL results in terms of the spatial error, which indicates that the AL results were closer to the theoretical position. In addition, in terms of the error of the origin time and the minimum value of the objective function, the AL results were also better than the BL results.
3.1.3. Comparison and Analysis
Due to the influence of noise in the tunnel, there was a certain error in the arrival picking in actual projects. To further test the practicality of the HVM-based method, we added a certain noise in the theoretical arrival time (see
Table 2). We designed three contrasting experiments. First, the theoretical arrival time was used as the location input parameter, and an SVM was adopted to locate the three sources, which was expressed using the SVM. Second, the noisy arrival time was used as the location input parameter, and an HVM was adopted to locate the three sources, which was represented using the HVM(N), where N is stand for noisy arrival time. Third, the noisy arrival time is used as the location input parameter, and an SVM was adopted to locate the three sources, which is represented using the SVM(N). The location results of the three experiments are shown in
Table 4, which clearly shows that the origin time obtained by these three experiments was very close to the theoretical origin time. However, the location results of the SVM and SVM(N) were significantly different from those of the HVM(N) in the Y direction.
We used these three experiments to compare the errors of the location results of the HVM and SVM in detail, as shown in
Figure 8. The origin time error of these four experiments was within 0.6 ms, and the accuracy was very high, which was not used as a criterion. In the case where no noise was added, the location error of the SVM of the three sources in the X direction was smaller than that of the HVM, but the location error of the HVM in the Y direction was smaller than that in the SVM. The location errors of the SVM in the Z direction of S1 and S3 were smaller than those of the HVM, but the opposite was true for S2. Considering the location errors in the X, Y, and Z directions, it was impossible to determine which method had a higher location accuracy. The spatial location error was a comprehensive error in the integrated X, Y, and Z directions. Therefore, the spatial location error calculated using Equation (10) was used as the criterion.
Figure 8 clearly shows that the spatial location error of the HVM was smaller than the spatial location error of the SVM. For the average spatial location errors of the three sources, HVM (2.06 m) < SVM (5.51 m). In addition, from the minimum value of the objective function, the value of the HVM-based MEL was smaller. In the case with the added noise, the spatial location error and the residual of the HVM-based method were smaller than those of the SVM-based method. Regarding the average spatial location error of the three sources, HVM (4.95 m) < SVM (6.23 m). In summary, it can be concluded that the location accuracy of the HVM was higher than that of the SVM. Below, we use the proposed method to analyze real data.
3.2. Application to Real Data
The Qinling No. 4 tunnel of the Yinhanjiwei project is located south of the Qin Mountains in southern Shaanxi Province, China. The length of this shaft is 5820 m, with a section 6.5 m high and 6.7 m wide. The maximum slope is 11.96%, and the maximum depth is 1600 m. Drilling and blasting are used in this project, resulting in frequent rockbursts that may both damage infrastructure and injure people. To provide safety guidance, a microseismic system was used for 24 h of continuous monitoring. Four accelerometers with a sensitivity of 10 V/g were embedded along two sides of the tunnel. The spatial arrangement of the receiver is shown in
Figure 9a. The sampling frequency was set to 10 kHz.
The working face of the actual project is shown in
Figure 9d. Due to the poor lighting inside the tunnel, to distinguish the microseismic hole, the sensor installation position was marked with red paint for visibility, as shown in
Figure 9c. Some ejected fragments were found at the top of the working face on 15 June 2017, as shown in
Figure 9b. According to engineering experience, these fragments formed due to the destruction of the roof rock mass caused by the excavation of the tunnel. Below, we verify this through microseismic monitoring.
During the tunneling process, the microseismic monitoring system monitored a large number of events. The microseismic event and the blasting event could be clearly distinguished by the waveform, as shown in
Figure 10. We selected 51 events with better waveforms during the period from 6 June 2017 to 13 June 2017, including 7 blast events and 44 microseismic events. We verified the location effect in two ways: (1) due to the known blasting position of the working face, we verified the location effect based on 7 blasting events; and (2) we verified the location effect based on the spatial relationship between the location results of 44 microseismic events and the working surface.
We used the HVM to locate these 51 events. The monitoring area covered the x coordinates from 3,727,271 m to 3,727,516 m, y coordinates from 502,564 m to 502,761 m, and z coordinates from 558 m to 598 m. The block size was 0.5 m × 0.5 m × 0.5 m, so there were a total of 490 × 394 × 80 unit blocks. There are many ways to calculate the propagation velocity of waves in rock masses. For example, Wang et al. optimized the seismic wave velocity in the deep mining area of a coal mine by using a combination method, residual error optimization method, location error optimization method, location residual optimization method, and combined inversion method [
33]. The acquisition of the velocity model was not the focus of this paper. In this application, the wave propagated at a velocity of approximately 6000 m/s in the rock mass. Therefore, we set the velocity of the unexcavated rock mass to 6000 m/s and the velocity of the empty area after excavation to 340 m/s, as shown in
Figure 11. The coordinates of the receivers and the arrival time of the seven blasting events are shown in
Table 5. The arrival time of the 44 microseismic events is shown in
Appendix B.
The location results are shown in
Figure 12. In this figure, red indicates the location result of the HVM, and blue indicates the location result of the SVM; five-pointed stars indicate the blast events, and circles indicate the microseismic events. The results of the seven blasting events are shown in
Table 6. The results of the 44 microseismic events are shown in
Appendix C.
First, we analyzed the results of the HVM. From the location results of the blasting events, it can be clearly seen that the seven blasting events occurred in the vicinity of the working face, which was consistent with the actual engineering excavation. The location results of the microseismic events showed that 44 microseismic events were concentrated around the blasting events. The blasting of the tunnel face caused damage to the surrounding rock mass. The microseismic signal from the rock mass damage could be received by the sensor; therefore, in theory, most of the events occurred near the working face. The location results of the 44 microseismic events were consistent with this theoretical derivation.
Then, we analyzed the location results of the SVM. As shown in
Figure 12, the location results of the blasting events suggested that the SVM results were very scattered. Four of the blasting events were located behind the propulsion direction of the working face, and two blasting events were located approximately 85 m in front of the propulsion direction of the working face. From the location results of the microseismic events, the spatial distribution of the 44 microseismic events was consistent with that of the seven blasting events. The spatial distribution was very scattered, irregular, and not concentrated near the working surface.
In summary, we can conclude that the HVM had a high location accuracy and good effect, with a location accuracy that was much higher than that of the SVM.
What caused the location results of the SVM-based method to be so poor? We compared these results with the HVM-based MEL and found that the velocity model error decreased the location accuracy. After excavation, the tunnel was filled with air, and the propagation of the microseismic signal was greatly affected. The wave did not pass directly through the empty zone to the receiver, but rather bypassed the empty zone and reached the receiver by travelling through the unexcavated rock mass. Therefore, in the MEL in the tunnel, we should have fully considered the impact of the empty zone and used the HVM for location.
The HVM-based method proposed in this paper had a high precision and was suitable for MEL during tunnel excavation. The location accuracy of the HVM varied with the size of the mesh. The larger the mesh size was, the lower the accuracy. Conversely, the smaller the mesh size was, the higher the accuracy. However, the smaller the grid was, the lower the computational efficiency. Therefore, in practical engineering applications, we should consider both accuracy and efficiency when determining the appropriate grid size.