Next Article in Journal
A Discriminative Multi-Output Gaussian Processes Scheme for Brain Electrical Activity Analysis
Next Article in Special Issue
Improvement in Accuracy of a Multi-Joint Robotic Ultrasonic Inspection System for the Integrity of Composite Structures
Previous Article in Journal
A New Technology for Smooth Blasting without Detonating Cord for Rock Tunnel Excavation
Previous Article in Special Issue
Microseismic Signal Denoising and Separation Based on Fully Convolutional Encoder–Decoder Network
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Fast Ray-tracing Method for Locating Mining-Induced Seismicity by Considering Underground Voids

1
School of Resources and Safety Engineering, Central South University, Changsha 410083, China
2
Digital Mine Research Center, Central South University, Changsha 410083, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2020, 10(19), 6763; https://doi.org/10.3390/app10196763
Submission received: 24 August 2020 / Revised: 25 September 2020 / Accepted: 25 September 2020 / Published: 27 September 2020

Abstract

:

Featured Application

For underground mines in multi-voids, the proposed method can effectively improve the accuracy and efficiency of ray tracing, thereby improving the accuracy of the microseismic event location.

Abstract

The accurate localization of mining-induced seismicity is crucial to underground mines. However, the constant velocity model is used by traditional location methods without considering the great difference in wave velocity between rock mass and underground voids. In this paper, to improve the microseismicity location accuracy in mines, we present a fast ray-tracing method to calculate the ray path and travel time from source to receiver considering underground voids. First, we divide the microseismic monitoring area into two categories of mediums—voids and non-voids—using a flexible triangular patch to model the surface model of voids, which can accurately describe any complicated three-dimensional (3D) shape. Second, the nodes are divided into two categories. The first category of the nodes is the vertex of the model, and the second category of the nodes is arranged at a certain step length on each edge of the 3D surface model to improve the accuracy of ray tracing. Finally, the set of adjacent nodes of each node is calculated, and then we obtain the shortest travel time from the source to the receiver based on the Dijkstra algorithm. The performance of the proposed method is tested by numerical simulation. Results show that the proposed method is faster and more accurate than the traditional ray-tracing methods. Besides, the proposed ray-tracing method is applied to the microseismic source localization in the Huangtupo Copper and Zinc Mine. The location accuracy is significantly improved compared with the traditional method using the constant velocity model and the FMM-based location method.

1. Introduction

The dramatic increase of the global population and extensive urbanization has similarly increased the requirement for minerals from the mining industry. With the depletion of near-surface ore deposits, underground mining is being increasingly employed all around the world. Microseismic monitoring is an emerging technology used in underground mines to identify seismic sources, to achieve an understanding of the seismic hazard and, consequently, to ensure safety inside the mine [1,2,3,4].
Accurate location of mining-induced seismicity is crucial to the successful application of microseismic monitoring, which is usually performed using a constant velocity model [5,6]. However, the excavation of ore deposits using underground mining inevitably results in lots of underground voids (also known as “goaves” or “stopes”). A constant velocity model without considering the great difference between rock mass and underground voids will result in large location error. To solve this problem, Feng [7] proposed a sectional velocity model in the tunnel, but the voids of the mine are many and complicated, and the velocity from each sensor to the source is different. Furthermore, Feng [8] proposed an anisotropic velocity model for location. This method requires association with rockburst events to solve the underdetermined equation, and it is not suitable for situations where there is no rockburst information. From another perspective, Virieux and Zhao [9,10] use a ray-tracing method to greatly improve the accuracy of earthquake location. For small-scale underground mines with a large number of known voids, it is appropriate to use ray tracing to improve the accuracy of the microseismic location.
Ray tracing is a method to find the propagation path of the wave from the source to the receiver. Traditional ray-tracing methods include shooting methods [11,12,13] and bending methods [14,15]. However, the accuracy and efficiency of such methods in complicated models are notably low, and they cannot satisfy the computational requirements of complicated rock masses in engineering. Newer methods include wavefront reconstruction methods [16], slowness matching methods [17], shortest path ray-tracing (SPR) methods [18,19], simulated annealing methods [20], etc. Among them, the finite-difference (FD) approximation to the eikonal equation is very popular to calculate the ray path. In the field of earthquakes, the eikonal solver is used for many seismological problems with very large data sets. The eikonal solvers based on the fast-marching method (FMM) [21] are also common methods in the field of microseismic. SM Hassouna developed multi-stencils fast marching methods (MSFMM) [22] based on FMM to improve accuracy. Moreover, the FMM is also applied to the microseismic location by many scholars [23,24,25]. The FMM used to solve the eikonal equation implies a plane wave approximation, so it causes larger errors in places close to the source [26]. To overcome this singularity at the source, Zhang [27] and Fomel [28] use a simpler algorithm based on the factored eikonal equation, but this algorithm causes the error far away from the source to become larger. To fulfill all these objectives, Noble [29] develops a hybrid algorithm that combines a spherical approximation by solving the factored eikonal equation when close to the source and a plane wave approximation by solving the standard eikonal equation when far away from the source. However, one disadvantage of the above method is that it needs to be meshed to calculate the ray path, and a regular cuboid grid is generally used. The complicated geometry of the geological media requires fine-scale grids that may make computationally intensive, particularly for 3D models [30]. The more grid nodes, the lower the computational efficiency.
Therefore, the above methods are not suitable for underground mines with many voids. To solve this problem, for the special situation of underground mines, we propose a ray-tracing method that does not need to be divided into grids, which not only improves efficiency but also ensures accuracy. Our proposed method has the following advantages: (1) Any geological model can be accurately described by the use of triangular patches so that a velocity model can be established accurately. (2) There is no need to divide a regular cuboid grid, which reduces a lot of calculations and improves efficiency. (3) Based on the idea of the Dijkstra algorithm, while improving efficiency, the accuracy of ray tracing is higher than that of traditional methods.

2. Methodology

As the propagation medium of the wave in the voids is air, the wave velocity in the voids is 340 m/s. Besides, the energy attenuation is more significant than that in the non-voids. The propagation medium of the wave in the non-voids is the rock mass, and the propagation velocity of seismic waves in rocks is much higher than 340 m/s [31,32]. Obviously, for small-scale multi-voids of an underground mine, the wave propagation is mainly affected by the voids. Therefore, we divide the monitoring area into two media. The wave velocity in these two media is very different and the shortest path of the wave from the source to the receiver generally bypasses the void. Therefore, we can assume that the shortest path of the wave does not pass through voids.
The steps of the proposed method are as follows: (1) To establish the surface model of the underground voids, this paper uses the 3D modeling software based on the triangular patch modeling. (2) Interpolate nodes on the edge of the surface model according to a certain step length, and the step length controls the accuracy of ray tracing. (3) The key step of the proposed method is to calculate the adjacent nodes set of each node, and the definition of adjacent points is in the method section of the paper. (4) Based on Fermat’s principle, the Dijkstra algorithm is used to calculate the minimum value of the arrival time of each node, and then the shortest path from the source to the sensor is obtained. In summary, we propose a new ray-tracing method without meshing. This new method avoids dividing a large number of grids and provides a fast and accurate method for microseismic source location in underground mining.

2.1. Fermat Principle

The Fermat principle [33] is one of the basic principles in the ray theory of seismic waves, which reflects the propagation law of seismic waves. The popular expression of the Fermat principle is that the travel time of the seismic wave along a ray is minimal compared to the travel time along other paths. As shown in Figure 1, in the uniform velocity model, the propagation path from point A to point B is a straight line (Figure 1a); in the layered velocity model, the propagation path from point A to point B is a polyline. Moreover, the inflection point is on the interface (Figure 1b). In the continuously changing velocity model, the propagation path from points A to B is a curve (Figure 1c). The propagation path of the wave varies with the velocity model. However, no matter how it changes, the propagation path always satisfies the Fermat principle. For underground mines, we only need to find the shortest path from the source to the receiver in the non-voids.

2.2. A Fast Ray-tracing Method

Since the monitoring area is divided into two parts: voids and non-voids, we only need to build the voids model, and the rest are the non-voids. Finally, we describe the specific process of the proposed method in detail.

2.2.1. Modeling

Minghui Zhang uses triangular grids to describe the irregular surface [34] and obtains good results. The shape of the voids in the actual engineering is notably complicated, so the triangular patch is used to describe the complicated surface of the 3D voids. As shown in Figure 2, a simple void is created. The black circle is the vertex of the void, and there are 4 vertices, 6 edges and 4 faces in this void.

2.2.2. Implementation

After modeling, we perform ray tracing. First, to improve the accuracy of ray tracing, we arrange nodes on the edge according to the step length s. The nodes are divided into two categories. The first category of nodes is the vertex of the void, and the second category of nodes is arranged on the edge. The step size s controls the accuracy of ray tracing. The smaller the s, the more the number of the second category of nodes. Generally, we set the step size s according to the actual situation.
  • Arrange the second category nodes
We first make the following statement: Nc, Np represent the number of vertices, and the number of triangular faces of the model, respectively. v represents the propagation velocity of the wave in the non-voids. The number of edges in the model is 3Np/2. The number of nodes we arrange on the i-th edge is:
L N i = [ D i / s ]
where LNi represents the number of nodes arranged on the i-th edge. Di is the length of the i-th edge and [X] represents the largest integer not greater than X. Then the total number of nodes arranged on the edges is:
N a = i = 1 n L N i , n 3 N p / 2
where Na is the total number of the second category nodes. The number of nodes of the first category plus that of the second category is the total number of nodes.
N s = N c + N a
where Ns is the total number of nodes.
Let us take the model in Figure 3 as an example. There are 4 nodes of the first category (black circle) and 13 nodes of the second category (red circle). If we want to improve the accuracy of ray tracing, we only need to reduce s. The number of nodes of the second category increases, but the number of nodes of the first category remains unchanged.
2.
Adjacent nodes of each node
In addition to the two categories of nodes on the above model, we must have a source and a receiver. Therefore, we get a collection of all data, denoted as Datas.
D a t a s = { N , S , R }
where N is the collection of nodes on the model. S and R are the source and receiver, respectively.
The definition of adjacent nodes is as follows: If the straight line from node A to node B does not pass through the voids, then A and B are adjacent nodes to each other; otherwise, they are not adjacent nodes.
The method of judging adjacent nodes is as follows: We turn this question into a question of whether a certain point is in a void. See reference [23] for the algorithm to determine whether the point is in the void. As shown in Figure 4, we first judge whether the midpoint C of node A and node B is in the void. If C is in the void, A and B are not adjacent nodes; otherwise, continue to judge whether D and E are in the void until the allowable tolerance c is reached. Please see Figure 5 for the detailed process.
3.
Ray tracing
Finally, we perform the shortest path from the source to the receiver based on the Dijkstra algorithm. We divide the Datas obtained from Equation (4) into three sets. P represents the set of nodes that have not got the arrival time. Q represents the set of nodes that have got the arrival time but have not been the source. K represents the set of nodes that have been the source.
Step 1: Initialization
P = { N , R } , Q = , K = S
Step 2: Calculate the arrival time from the source S to each of its adjacent nodes according to Equation (6).
t i = ( x s x i ) 2 + ( y s y i ) 2 + ( z s z i ) 2 / v + t s , i J ( s )
where t s and t i are the arrival time of S and the i-th node, respectively. J(s) is the set of adjacent nodes of S. Then move J(s) to the set Q and record the source that minimizes the arrival time of the i-th node.
Step 3: In the set Q, the node with the smallest arrival time is selected as the source in the next iteration.
S c = min ( t j ) , j Q
where S c is the source at the c-th iteration.
Calculate the arrival time from S c to all its adjacent nodes. The arrival time of the i-th node is the minimum of the arrival time t i c obtained in the c-th iteration and its arrival time t i before this iteration.
t i = min ( t i , t i c )
Step 4: Continue to perform step 3 until the set Q is empty or all receivers R are in the set K, exit the loop.
Step 5: Starting from the receiver, search for the node with the smallest arrival time in the record until the source S, connect all its points, and get the shortest path.
Below we use Figure 6 to illustrate all the steps of ray tracing. Step 1—initialization, P = { A , B , C , 1 , 2 , 3 , R } , Q = , K = { S } . Step 2—use Equation (6) to calculate the arrival time from source S to nodes A, B, and 1. That is P = { C , 2 , 3 , R } , Q = { A , B , 1 } , K = { S } . Step 3—select the node B with the smallest arrival time from the set Q as the source of the second iteration. That is P = { 2 , R } , Q = { A , 1 , 3 , C } , K = { S , B } . Continue to step 3 until the 8th iteration. That is P = , Q = , K = { S , B , 1 , 3 , A , C , 2 , R } . Step 5—the shortest path from S to R is S-A-R.

3. Synthetic Tests

In this part, we verify the accuracy and efficiency of the proposed method. First, we use two cases to test whether the proposed method is feasible. Then, the proposed method is compared with the classic FMM and SPR. Finally, we analyze the factors that affect the ray tracing accuracy of the proposed method.

3.1. Experiments

3.1.1. Case 1

We first established a void model, as shown in Figure 7. This model has 180 triangular faces and 92 vertices (the first category of nodes). The vertices are represented by green circles in the figure. Then according to Equations (1)–(3), the number of nodes of the second category is calculated as 3304 (step size is 1). The red circle in the figure represents the second category of nodes. The total number of nodes is 3396.
The scope of the monitoring area is x [ 0 , 100 ] , y [ 0 , 100 ] , z [ 0 , 100 ] , as shown in Figure 8. Eighteen receivers (green inverted triangle) and one source (red circle) are arranged in the monitoring area. The propagation velocity of the wave in the non-void is set to 5000 m/s. Then, we use the proposed method to calculate the shortest paths from the source to the receivers, as shown by the red line in Figure 8. We can see that the ray paths from the source to the receivers bypass the void and finds the shortest paths that propagate in the non-void along the surface of the void.

3.1.2. Case 2

Based on Case 1, we created a new void within the monitoring range, as shown in Figure 9. The void 2 is a rectangular parallelepiped model with a range of x [ 40 , 70 ] ,   y [ 40 , 70 ] ,   z [ 40 , 70 ] . The model consists of 8 vertices and 12 triangular faces. The number of second category nodes on void 2 is 792 (step length is 1). One hundred receivers are arranged at 10-meter intervals in the x and z directions (the green inverted triangle in Figure 9) and one source (the red circle). The ray paths (the red line) from the source to the receivers are calculated by the proposed method. It can be seen that the ray path bypasses the voids 1 and 2, and finds the shortest path from the source to the receiver. The effect of ray tracing is very good, reaching the expected result.

3.2. Comparison and Analysis

We use an example to compare the proposed method with the classic FMM and SPR. As shown in Figure 10, this void model consists of 8 vertices and 12 faces, with a range of x [ 40 , 70 ] ,   y [ 40 , 70 ] ,   z [ 40 , 70 ] . The source coordinates are (0,50,50) and the coordinates of the 25 receivers R1–R25 are shown in Table 1.
We arranged the second category of nodes on the edges of the model with a step length of 1. According to Equation (4), the total number of nodes is 634. Then, we used the proposed method for ray tracing. The travel time of the wave from the source to the receivers is shown in Table 1. Similarly, we use classic FMM and SPR for ray tracing. The grid size is 1 × 1 × 1, and there is a total of 100 × 100 × 100 grid nodes. The results obtained with FMM and SPR are shown in Table 1. The difference between the arrival time calculated by the proposed method and the theoretical value is less than 10−5 s, so the theoretical arrival time is not shown in Table 1.
In Figure 10, we take the receiver R20 as an example to show the ray paths calculated by three methods. The ray path calculated by the proposed method coincides with the theoretical ray path. It also can be seen that the ray path calculated by the proposed method is similar to that calculated by FMM, but is quite different from that calculated by SPR. We also compared the error between the arrival time of R1–R20 and the corresponding theoretical value, as shown in Figure 11. The error of arrival time calculated by our proposed method is the smallest and the maximum error is less than 0.0093%. According to Fermat’s principle, the wave propagates along the minimum arrival path, so the ray path calculated by the proposed method is more accurate. SPR and FMM have the same grid size, but the error of the ray path calculated by SPR is greater.
Above we conclude that the ray-tracing accuracy of the proposed method is higher than that of the classic FMM and SPR. Next, we compare the efficiency of the proposed method with the classic FMM and SPR. In this example, the number of nodes in the proposed method is 634, and the number of nodes in FMM and SPR is 1 million. We greatly reduce the number of nodes. The proposed method takes 2 s, FMM takes 4 s, and SPR takes 34 s. In summary, the accuracy and efficiency of the proposed method are higher than those of classical FMM and SPR.
How is the accuracy of the proposed method affected? Let us look at another example. Figure 12 is a void model with a length of 50 m, a width of 30 m, and a height of 40 m. The source and receiver are arranged on the surface, and the coordinates are shown in Table 2. The propagation velocity of the wave is 5000 m/s. We set the step size s to four different values of 1, 5, 10, and 20 to arrange the second category of nodes. The results of ray tracing using these step lengths s are shown in Table 2. Since the void model is a regular cuboid, we calculate the theoretical arrival time from the source to the receiver based on the cuboid expansion principle, and the results obtained are shown in Table 2.
According to the analysis of Table 2, when step length s = 20 m, the error of arrival time is between 0.3630 and 1.892 ms. When step size s = 10 m, the error of arrival time is between 0.0180 and 0.3154 ms. When the step length s = 5 m, the error of arrival time is between 0.007 and 0.1007 ms. When step size s = 1 m, the error of arrival time is between 0.0006 and 0.0056 ms. It also can be seen from Figure 13 that as the step size s decreases, the arrival time error is smaller. Moreover, when s is less than 5, the change of arrival time error is no longer obvious. The step length s = 5 m, the error of the maximum arrival time is about 0.1 ms, and the corresponding distance error is 0.5 m, which fully meets the needs of microseismic location in underground mines. In short, the accuracy of our proposed method is only related to the step size s.
In conclusion, our proposed method is superior to the classic FMM and SPR in terms of accuracy and efficiency and is suitable for microseismic location in complex and multi-voids in underground mines. However, our proposed method also has certain limitations. It cannot be applied to diversified velocity models. Table 3 lists the comparison of the proposed method with FMM and SPR.

4. Field Application

The Huangtupo Copper and Zinc Mine is located in the southwest of Hami City, Xinjiang Uygur Autonomous Region, China. As the mine uses non-pillar sublevel caving, three large voids have been formed so far, as shown in Figure 14. The volumes of the first, second, and third voids are 120,068.60, 42,633.25, and 183,483.19 m3, respectively. We have arranged a total of eight sensors in the middle of 210 and 260 m (Figure 14), with a sensitivity of 10 V/g, a sampling frequency of 10 V/g, and a microseismic system for continuous 24-h monitoring. The propagation velocity of the wave in the non-voids is about 5500 m/s.
As the existence of the three voids, the location of the microseismicity that propagates through voids using a constant velocity model will result in a large error. Thus, to show the improvement by integrating our proposed ray-tracing method into the location procedure, controlled blast experiments with small amounts of explosives were carried out in three different locations in the monitoring area. The coordinates of these blasts were pre-measured and listed in Table 4. The location method used in this study is proposed by Peng [25]. The only difference is that the ray-tracing method mentioned in this article is used instead of FMM. The calculated blasting source location and location error are shown in Table 4. Then, we compared the location results based on the proposed method with those based on FMM and those using constant velocity. As shown in Table 4, the mean location error is 5.88 m by integrating our proposed ray-tracing method, while the mean location error is 7.22 m using the FMM-based location method, and the mean location error of the method with constant velocity is 22.14 m. The accuracy is greatly improved, which indicates that our method is effective in the location of mining-induced seismicity.

5. Conclusions

In this paper, we propose a new algorithm of ray tracing for the location of mining-induced seismicity considering underground voids. The triangular patch is used to describe the complicated void model. The accuracy of ray tracing is controlled by arranging the step length of nodes on the edges of the model. A fast and accurate ray-tracing method is proposed based on Fermat’s principle and Dijkstra algorithm. We tested the effectiveness of the new method with models of single-void and multiple-voids and compared the new method with classic FMM and SPR. The results show that the accuracy and efficiency of the new method are higher than the classic ray-tracing method. The biggest advantage of our proposed method is that there is no need to divide the grid, and only a small number of nodes can be used to obtain high-precision ray paths. Finally, our method was applied to locate the pre-measured blasts in Huangtupo Copper and Zinc Mine. By integrating our proposed ray-tracing method, the mean location error is improved from 22.14 to 5.88 m compared with that using a constant velocity model, which indicates that our method is effective in the location of mining-induced seismicity.

Author Contributions

Y.J. and P.P. participated in data analysis, participated in the design of the study and drafted the manuscript, carried out the statistical analyses; Z.H. and S.T. collected field data; L.W. conceived of the study, designed the study, coordinated the study, and helped draft the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by The National Key R&D Program of China (2017YFC0602905) and the Fundamental Research Funds for the Central Universities of Central South University (2020zzts194).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wang, H.-L.; Ge, M.-C. Acoustic emission/microseismic source location analysis for a limestone mine exhibiting high horizontal stresses. Int. J. Rock Mech. Min. Sci. 2008, 45, 720–728. [Google Scholar] [CrossRef]
  2. Le Gonidec, Y.; Sarout, J.; Wassermann, J.; Nussbaum, C. Damage initiation and propagation assessed from stress-induced microseismic events during a mine-by test in the Opalinus Clay. Geophys. J. Int. 2014, 198, 126–139. [Google Scholar] [CrossRef] [Green Version]
  3. Cheng, J.; Song, G.; Sun, X.; Wen, L.; Li, F.; Wena, L.; Li, F.; Wen, L.; Li, F. Research developments and prospects on microseismic source location in mines. Engineering 2018, 4, 653–660. [Google Scholar] [CrossRef]
  4. Kinscher, J.; Bernard, P.; Contrucci, I.; Mangeney, A.; Piguet, J.P.; Bigarre, P. Location of microseismic swarms induced by salt solution mining. Geophys. J. Int. 2015, 200, 337–362. [Google Scholar] [CrossRef] [Green Version]
  5. Feng, G.-L.; Feng, X.-T.; Chen, B.-R.; Xiao, Y.-X. Performance and feasibility analysis of two microseismic location methods used in tunnel engineering. Tunn. Undergr. Space Technol. 2017, 63, 183–193. [Google Scholar] [CrossRef]
  6. Xu, N.W.; Li, T.B.; Dai, F.; Zhang, R.; Tang, C.A.; Tang, L.X. Microseismic monitoring of strainburst activities in deep tunnels at the Jinping II hydropower station, China. Rock Mech. Rock Eng. 2016, 49, 981–1000. [Google Scholar] [CrossRef]
  7. Feng, G.-L.; Feng, X.-T.; Chen, B.-R.; Xiao, Y.-X.; Jiang, Q. Sectional velocity model for microseismic source location in tunnels. Tunn. Undergr. Space Technol. 2014, 45, 73–83. [Google Scholar] [CrossRef]
  8. Feng, G.L.; Feng, X.T.; Chen, B.R.; Xiao, Y.X. A highly accurate method of locating microseismic events associated with rockburst development processes in tunnels. IEEE Access 2017, 5, 27722–27731. [Google Scholar] [CrossRef]
  9. Virieux, J.; Farra, V.; Madariaga, R. Ray tracing for earthquake location in laterally heterogeneous media. J. Geophys. Res. 1988, 93, 6585. [Google Scholar] [CrossRef]
  10. Zhao, A.-H.; Ding, Z.-F.; Bai, Z.-M. Improvement of the ray-tracing based method for calculating hypocentral Loci for earthquake location. Chin. J. Geophys. 2015, 58, 701–717. [Google Scholar]
  11. Langan, R.T.; Lerche, I.; Cutler, R.T. Tracing of rays through heterogeneous media: An accurate and efficient procedure. Geophysics 1985, 50, 1456–1465. [Google Scholar] [CrossRef]
  12. Virieux, J.; Farra, V. Ray tracing in 3-D complex isotropic media: An analysis of the problem. Geophysics 1991, 56, 2057–2069. [Google Scholar] [CrossRef] [Green Version]
  13. Xu, T.; Zhang, Z.; Zhao, A.; Zhang, A.; Zhang, X.; Zhang, H. Subtriangle shooting ray tracing in complicated 3D VTI media. J. Seism. Explor. 2008, 17, 133–146. [Google Scholar]
  14. Xu, T.; Xu, G.; Gao, E.; Li, Y.; Jiang, X.; Luo, K. Block modeling and segmentally iterative ray tracing in complex 3D media. Geophysics 2006, 71, T41–T51. [Google Scholar] [CrossRef]
  15. Um, J.; Thurber, C. A fast algorithm for two-point seismic ray tracing. Bull. Seismol. Soc. Am. 1987, 77, 972–986. [Google Scholar]
  16. Vinje, V.; Iversen, E.; Gjoystdal, H. Traveltime and amplitude estimation using wavefront construction. Geophysics 1993, 58, 1157–1166. [Google Scholar] [CrossRef]
  17. Symes, W.W.; Qian, J. A Slowness matching eulerian method for multivalued solutions of eikonal equations. J. Sci. Comput. 2003, 19, 501–526. [Google Scholar] [CrossRef]
  18. Moser, T.J. Shortest path calculation of seismic rays. Geophysics 1991, 56, 59–67. [Google Scholar] [CrossRef]
  19. Zhou, B.; Greenhalgh, S.A. “Shortest path” ray tracing for most general 2D/3D anisotropic media. J. Geophys. Eng. 2005, 2, 54–63. [Google Scholar] [CrossRef]
  20. Bóna, A.; Slawinski, M.A.; Smith, P. Ray tracing by simulated annealing: Bending method. Geophysics 2009, 74, 25–32. [Google Scholar] [CrossRef] [Green Version]
  21. Sethian, J.A. A fast marching level set method for monotonically advancing fronts. Proc. Natl. Acad. Sci. USA 1996, 93, 1591–1595. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Hassouna, M.S.; Farag, A.A. Multistencils fast marching methods: A highly accurate solution to the Eikonal equation on Cartesian domains. IEEE Trans. Pattern Anal. Mach. Intell. 2007, 29, 1563–1574. [Google Scholar] [CrossRef] [PubMed]
  23. Peng, P.; Jiang, Y.; Wang, L.; He, Z. Microseismic event location by considering the influence of the empty area in an excavated tunnel. Sensors 2020, 20, 574. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Peng, P.; Wang, L. 3DMRT: A computer package for 3D model-based seismic wave propagation. Seismol. Res. Lett. 2019, 90, 2039–2045. [Google Scholar] [CrossRef]
  25. Peng, P.; Wang, L. Targeted location of microseismic events based on a 3D heterogeneous velocity model in underground mining. PLoS ONE 2019, 14, e0212881. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Qian, J.; Symes, W.W. An adaptive finite-difference method for traveltimes and amplitudes. Geophysics 2002, 67, 167–176. [Google Scholar] [CrossRef] [Green Version]
  27. Zhang, L.; Rector, J.W.; Hoversten, G.M. Eikonal solver in the celerity domain. Geophys. J. Int. 2005, 162, 1–8. [Google Scholar] [CrossRef] [Green Version]
  28. Fomel, S.; Luo, S.; Zhao, H. Fast sweeping method for the factored eikonal equation. J. Comput. Phys. 2009, 228, 6440–6455. [Google Scholar] [CrossRef]
  29. Noble, M.; Gesret, A.; Belayouni, N. Accurate 3-D finite difference computation of traveltimes in strongly heterogeneous media. Geophys. J. Int. 2014, 199, 1572–1585. [Google Scholar] [CrossRef]
  30. Zhao, D.; Lei, J. Seismic ray path variations in a 3D global velocity model. Phys. Earth Planet. Inter. 2004, 141, 153–166. [Google Scholar] [CrossRef]
  31. Jiang, F.; Ye, G.; Wang, C.; Zhang, D.; Guan, Y. Application of high-precision microseismic monitoring technique to water inrush monitoring in coal mine. Chin. J. Rock Mech. Eng. 2008, 27, 1932–1938. [Google Scholar]
  32. Wang, J.; Jiang, F.; Lu, W.; Wang, C. Microseismic wave propagation velocity insitu experiment and calculation. J. China Coal Soc. 2010, 35, 2059–2063. [Google Scholar]
  33. Vesnaver, A.L. Ray tracing based on Fermat’s principle in irregular grids. Geophys. Prospect. 1996, 44, 741–760. [Google Scholar] [CrossRef]
  34. Zhang, M.; Xu, T.; Bai, Z.; Liu, Y.; Hou, J.; Yu, G. Ray tracing of turning wave in elliptically anisotropic media with an irregular surface. Earthq. Sci. 2017, 30, 219–228. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Ray path from point A to B in different velocity models, (a) Uniform velocity model; (b) Layered velocity model; (c) Continuously changing velocity model.
Figure 1. Ray path from point A to B in different velocity models, (a) Uniform velocity model; (b) Layered velocity model; (c) Continuously changing velocity model.
Applsci 10 06763 g001
Figure 2. A simple void model.
Figure 2. A simple void model.
Applsci 10 06763 g002
Figure 3. Two categories of nodes on the model. The black circle and the red circle represent the first category node and the second category node, respectively.
Figure 3. Two categories of nodes on the model. The black circle and the red circle represent the first category node and the second category node, respectively.
Applsci 10 06763 g003
Figure 4. Schematic diagram of dividing the line AB into n equal parts. C and D, E are the points where the line segment AB is halved and quartered, respectively.
Figure 4. Schematic diagram of dividing the line AB into n equal parts. C and D, E are the points where the line segment AB is halved and quartered, respectively.
Applsci 10 06763 g004
Figure 5. The process of judging whether node A and node B are adjacent nodes.
Figure 5. The process of judging whether node A and node B are adjacent nodes.
Applsci 10 06763 g005
Figure 6. Schematic diagram of calculating the shortest path.
Figure 6. Schematic diagram of calculating the shortest path.
Applsci 10 06763 g006
Figure 7. A void model and nodes layout. The green circle and the red circle represent the first category of nodes (vertex) and the second category of nodes, respectively.
Figure 7. A void model and nodes layout. The green circle and the red circle represent the first category of nodes (vertex) and the second category of nodes, respectively.
Applsci 10 06763 g007
Figure 8. In the case of a single void, the ray path is calculated by the proposed method.
Figure 8. In the case of a single void, the ray path is calculated by the proposed method.
Applsci 10 06763 g008
Figure 9. In the case of two voids, the ray path is calculated by the proposed method.
Figure 9. In the case of two voids, the ray path is calculated by the proposed method.
Applsci 10 06763 g009
Figure 10. The ray path of the receiver R20 is calculated by the proposed method, FMM, and SPR. (a) The view on XY plane, (b) the view on YZ plane, and (c) the view on XZ plane.
Figure 10. The ray path of the receiver R20 is calculated by the proposed method, FMM, and SPR. (a) The view on XY plane, (b) the view on YZ plane, and (c) the view on XZ plane.
Applsci 10 06763 g010aApplsci 10 06763 g010b
Figure 11. The relative error of the arrival time of the three methods compared to the analysis solution.
Figure 11. The relative error of the arrival time of the three methods compared to the analysis solution.
Applsci 10 06763 g011
Figure 12. A cuboid void model and its source and receiver layout.
Figure 12. A cuboid void model and its source and receiver layout.
Applsci 10 06763 g012
Figure 13. The relative error of the arrival time of the different step size s compared to the analysis solution.
Figure 13. The relative error of the arrival time of the different step size s compared to the analysis solution.
Applsci 10 06763 g013
Figure 14. Distributions of voids and sensors in the Huangtupo Copper and Zinc Mine. P1, P2 are collection devices.
Figure 14. Distributions of voids and sensors in the Huangtupo Copper and Zinc Mine. P1, P2 are collection devices.
Applsci 10 06763 g014
Table 1. Receiver coordinates and travel time calculated using the proposed method, FMM, and SPR.
Table 1. Receiver coordinates and travel time calculated using the proposed method, FMM, and SPR.
R1R2R3R4R5R6R7R8R9
x (m)100100100100100100100100100
y (m)0000021212121
z (m)0214263840214263
Travel time (ms)Proposed method24.5023.1022.4222.5123.3723.1021.6220.8920.99
FMM24.5523.1622.4622.5623.4323.1621.6720.9321.03
SPR24.4923.6022.7022.9223.8123.6022.6121.7121.92
R10R11R12R13R14R15R16R17R18
x (m)100100100100100100100100100
y (m)214242424242636363
z (m)8402142638402142
Travel time (ms)Proposed method21.9122.4220.8920.3220.4321.3322.5120.9920.43
FMM21.9722.4620.9320.7120.6121.4922.5621.0320.61
SPR22.8422.7021.7121.0021.2221.9822.9221.9221.22
R19R20R21R22R23R24R25
x (m)100100100100100100100
y (m)63638484848484
z (m)6384021426384
Travel time (ms)Proposed method21.2721.4323.3721.9121.3321.4322.33
FMM22.8321.5923.4321.9721.4921.5922.68
SPR21.8622.1923.8122.8421.9822.1923.09
Table 2. The result of ray tracing corresponding to step length s.
Table 2. The result of ray tracing corresponding to step length s.
x (m)y (m)z (m)Arrival Time (ms)
Theorys = 1s = 5s = 10s = 20
S0253000000
R1507716.2716.3116.4116.6319.72
R25071216.0116.0716.0717.1119.48
R35071715.8215.8415.8416.8718.71
R45072215.6915.6915.8216.6617.76
R55072715.0415.1015.1015.6716.85
R65073214.0714.1414.1414.7216.02
R75012715.3115.3815.4216.5618.90
R850121215.0415.0915.1416.1818.89
R950121714.8314.8514.8515.8918.08
R1050122214.6914.7014.8215.8617.33
R1150122714.6114.6514.6515.6416.68
R1250123213.7913.8914.0614.6616.18
R135017714.3614.4414.4415.9118.00
R1450171214.0714.1114.2515.2918.41
R1550171713.8513.8713.8714.9017.52
R1650172213.7013.7213.8314.8716.68
R1750172713.6213.6713.6715.1915.91
R1850173213.4113.4813.5115.1315.29
R195022713.4113.4813.4815.1317.17
R2050221213.1113.1513.4214.5017.93
R2150221712.8712.9112.9113.9417.11
R2250222212.7012.7612.8513.8916.18
R2350222712.6212.7012.7114.3615.29
R2450223212.6112.6512.6514.5014.50
Table 3. Comparison of our proposed method, FMM, and SPR.
Table 3. Comparison of our proposed method, FMM, and SPR.
AccuracyEffectivenessGridScope of Application
Proposed methodHighFastNoUnderground mine with a lot of voids and the wave velocity in the non-voids is not much different.
FMMHigherFasterYesAll fields
SPRlowslowYesAll fields
Table 4. Location results of the controlled blast experiments.
Table 4. Location results of the controlled blast experiments.
BlastsB1B2B3Mean
Actual positionx (m)31,412,506.5031,412,516.9231,412,514.32
y (m)4,719,746.394,719,743.594,719,815.21
z (m)146.68119.19134.94
Location method with constant velocityx (m)31,412,512.1931,412,531.6231,412,533.47
y (m)4,719,758.604,719,726.154,719,802.70
z (m)137.57127.94146.69
Location error (m)16.2624.4425.7222.14
Location method based on FMMx (m)31,412,509.4431,412,511.3531,412,521.42
y (m)4,719,749.604,719,749.674,719,813.89
z (m)143.57122.49136.74
Location error (m)5.358.877.447.22
Location method based on the proposed ray tracingx (m)31,412,508.9831,412,520.5231,412,517.76
y (m)4,719,743.054,719,746.214,719,818.79
z (m)149.03114.29131.16
Location error (m)4.786.626.245.88

Share and Cite

MDPI and ACS Style

Peng, P.; Jiang, Y.; Wang, L.; He, Z.; Tu, S. A Fast Ray-tracing Method for Locating Mining-Induced Seismicity by Considering Underground Voids. Appl. Sci. 2020, 10, 6763. https://doi.org/10.3390/app10196763

AMA Style

Peng P, Jiang Y, Wang L, He Z, Tu S. A Fast Ray-tracing Method for Locating Mining-Induced Seismicity by Considering Underground Voids. Applied Sciences. 2020; 10(19):6763. https://doi.org/10.3390/app10196763

Chicago/Turabian Style

Peng, Pingan, Yuanjian Jiang, Liguan Wang, Zhengxiang He, and Siyu Tu. 2020. "A Fast Ray-tracing Method for Locating Mining-Induced Seismicity by Considering Underground Voids" Applied Sciences 10, no. 19: 6763. https://doi.org/10.3390/app10196763

APA Style

Peng, P., Jiang, Y., Wang, L., He, Z., & Tu, S. (2020). A Fast Ray-tracing Method for Locating Mining-Induced Seismicity by Considering Underground Voids. Applied Sciences, 10(19), 6763. https://doi.org/10.3390/app10196763

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop