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

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.


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

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.

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.
Appl. Sci. 2020, 10, x 3 of 15 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.

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.

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.

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.

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.

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.

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: 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: 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.

Ns Nc Na
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.

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.
where N is the collection of nodes on the model. S and R are the source and receiver, respectively.

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.

1.
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: where LN i 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: 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.
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. Appl. Sci. 2020, 10, 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.

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.
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. 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: Step 2: Calculate the arrival time from the source S to each of its adjacent nodes according to Equation (6).
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.
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 c i obtained in the c-th iteration and its arrival time t i before this iteration.
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 2-use Equation (6)  Step 5-the shortest path from S to R is S-A-R.
Step 2-use Equation (6) to calculate the arrival time from source S to nodes A, B, and 1. That is Step 3-select the node B with the smallest arrival time from the set Q as the source of the second iteration. That is Step 5-the shortest path from S to R is S-A-R.

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.

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.   , 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. 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. , 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.  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.

Case 2
Based on Case 1, we created a new void within the monitoring range, as shown in Figure 9.
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.

Comparison and Analysis
We use an example to compare the proposed method with the classic FMM and SPR. As shown in Figure 10,

Comparison and Analysis
We use an example to compare the proposed method with the classic FMM and SPR. As shown in Figure 10 Table 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.  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. 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. 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 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. 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.
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.

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 m 3 , respectively. We have arranged a total of eight sensors in the middle of 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.

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 m 3 , 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. Appl. Sci. 2020, 10, x 13 of 15 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.  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.

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.