Relocating Mining Microseismic Earthquakes in a 3-D Velocity Model Using a Windowed Cross-Correlation Technique

Microseismic source location is a fundamental research interest in real-time microseismic monitoring and hazard risk assessment, and it provides the basis for determining the fracture zones and calculating seismic source parameters of the microseismicity (e.g., the event magnitude, the focal mechanism). In the present work, we carefully relocate some microseismic earthquakes that occurred in the Yongshaba mine (China) using the shooting 3-D ray-tracing (3D-RT) method based on a high-resolution 3-D velocity model, which is determined by a large amount of seismic data. Furthermore, a semiautomatic waveform cut method based on the cross-correlation (CC) technique (WCC) is developed for a quick, robust and precise determination of the direct P-phase relative delay times. Compared with the full waveform cross-correlation (FWCC)-based method, the WCC-based method is free from the influence of complex later-phase waveform spectra. To verify the proposed method, we artificially generate the locations of 892 synthetic microseismic events, the P-phase travel times of which are calculated based on the given 3-D velocity model and the 3-D ray-tracing technique. The relocated hypocentres of these synthetic events obtained by the developed method are improved compared with those obtained by a homogeneous velocity model and the straight line-ray tracing based L1 norm time difference (TD) method (HV-SL-TD1), as well as those obtained by the 3-D velocity model and the straight line-ray tracing based L1 norm and L2 norm TD methods (3D-SL-TD1 and 3D-SL-TD2). Finally, we apply the proposed method to eight experimental blasts with pre-measured locations. The synthetic and application tests show that, despite some multi-path phenomena existing in the ray-tracing methods, the WCC-based method performs as well as the experienced manual picking method, and the results obtained by the 3D-RT location have smaller location errors (~30 m) than those of the homogeneous velocity model-based methods (>40 m). The TD-based method is slightly better than the trigger time (TT)-based method when using the same norm, and the location precision of the L1 norm-based methods have a better resilience to larger picking errors than that of the L2 norm-based methods. Further improvements can be obtained by improving the resolution of the 3-D velocity model and including spatial voids in the model (e.g., tunnels and cavities) for the inversion.


I. INTRODUCTION
Microseismic monitoring, as an advanced and effective tool for monitoring rock bursts and other dynamic disasters, has been widely applied in mining engineering [1]. Microseismic The associate editor coordinating the review of this manuscript and approving it for publication was X. Huang . source location is a fundamental research interest in microseismic monitoring and hazard assessment and provides the basis for determining the fracture zones and calculating source parameters of the microseismicity (e.g., event magnitude, focal mechanism). Generally, a problem in source location is minimizing the objective function, which is based on the time difference between the theoretical and picked arrival times. In other words, the accuracy of the source location is influenced by different source location methods and velocity models, and the picked arrival time/relative delay time errors can affect the location precision. To handle this issue, many researchers have focused on iterative techniques to minimize the objective function, such as the grid search [2], Newton's method [3], Geiger's method [4], simplex algorithm [5], OcTree approach [6], differential evolution approach [7], simulated annealing (SA) and particle swarm optimization (PSO) algorithm [8]. A comprehensive event location review can be found in [9]. For the above studies, the L1 and L2 norms are mostly selected to qualify the objective functions, where the L1 method is usually less sensitive to large outliers than the L2 method [3], [9]. Therefore, in this work we apply the L1 norm time difference method (method TD1 in Eq. 5) and a combined global grid search technique to locate the microseismic source.

A. THE VELOCITY MODEL
The accuracy of a microseismic source location can be strongly dependent on the velocity model of the rock mass. It is generally assumed that the rock mass is homogeneous in small regions, where a homogeneous velocity model is usually adopted for microseismic source location [3], [10], [11]. However, Dong et al. [12], [13] noted that the errors of a pre-measured velocity can influence the location result. To address this problem, they treated the homogeneous velocity model as an unknown parameter in the minimization of the objective function, and their results outperformed the pre-measured velocity locations. The above methods treat the velocity model as a homogeneous velocity model for one event calculation, which may lead to source location errors since there are usually velocity variations in both the horizontal and vertical directions in the field. Feng et al. [8] horizontally divided a tunnel into several homogeneous sections (sectional velocity model), and then applied the PSO technique to search for both the microseismic source location and the sectional velocity. A ray tracing location based on a 1-D layered velocity model was used in [9], [14], [15]. The 1-D velocity model has a better location accuracy than that of the homogeneous velocity model, but remains insufficient for accurate source locations in regions of complex structural zones, such as mining engineering [16]. Thus, a 3-D velocity model that accommodates an in-situ status should be adopted. Many researchers have applied a high-resolution 3-D velocity model for earthquake locations [6], [17], [18] and microseismic event locations of reservoirs [19], [20], petroleum exploration [21] and rock slopes [22]. These results have demonstrated the 3-D model's advantages over a simple 1-D velocity model. To the best of our knowledge, only a few studies have applied simple 3-D velocity models for microseismic source location in mines. Collins et al. [23] tried to use a 1-D layered velocity model with a low velocity zone for microseismic source location in a mine, and Peng and Wang [24] treated the P wave velocity in the mined-out regions and tunnels as the wave velocity of air and other zones as a homogeneous velocity model. Gharti et al. [25] considered velocities of air, rock and ore in their mining applications. However, in mining zones 3-D high-resolution velocity models are rarely considered for microseismic source location.

B. THE ARRIVAL TIME
The source location precision can also be influenced by the picked absolute arrival time/relative delay time. The short term average to long term average (STA/LTA) method by Allen [26] is the most famous absolute arrival time picking method, whose characteristic function (CF) considers both the amplitude and frequency changes when the P-phase arrives. The Akaike information criterion (AIC) algorithm has also been widely adopted. The algorithm divides a signal into two parts [1, k] and [k + 1, N ], and then calculates its AIC value at the dividing point k. This algorithm assumes that the windowed waveforms before and after the P-phase arrival time are formed by two different processes, and thus have the smallest AIC value at arrival onset [27]. Another commonly applied method is the higher-order statistics-based Kurtosis/Skewness method [28]. In addition to the above methods, the damped predominant period (T pd ) method [29], multiband frequency analysis [30], and wavelet-based picking method [31] have also been developed for the absolute arrival time picking.
The absolute arrival time picking usually utilizes the abrupt change of phase arrival amplitude of a single waveform. However, the background, stationary noises and wavefront healing effect can make the P-phase rather unclear, which results in a problematic dataset for the manually and automatically picked absolute arrival times. Wavefront healing is a phenomenon where waves bend around an enclosed low-velocity/high-velocity anomaly, so that the plane wave outside of the low-velocity/high-velocity anomaly transmits quicker/slower than that in the low-velocity/high-velocity anomaly [32]. However, the wave fronts undergo healing with an increasing distance from the anomaly, which makes the P phase arrival amplitude smaller and unclear and the picked P phase arrival time becomes different from the theoretical arrival time. The relative delay time picking approach takes advantage of the time difference of the main peak between the two windowed waveforms, which is less affected by the arrival P-phase amplitude (better noise resilience). Luo and Schuster [33] proposed a cross-correlation technique to measure the time shift between the observed and synthetic waveforms. Since then, this technique has attracted widespread attention. Studies have shown that cross-correlation is a good technique for calculating the relative delay time between two earthquake waveforms [34], [35]. Trojanowski et al. [36] proposed a multi-channel convolution filter (MCCF) for correlated noise, while VanDecar and Crosson [37] used a multi-channel cross-correlation (MCCC) to pick the relative phase arrival times for teleseismic events. De Meersman et al. [38] proposed an iterative algorithm-based crosscorrelation to refine the initially picked arrival times for VOLUME 8, 2020 microseismic data, which was accomplished using Python software [39]. Leeuwen and Mulder [40] used a weighted norm of the cross-correlation to implicitly measure the time shift, which is insensitive to the differences in the amplitude spectra and can be executed automatically. Instead of using the time series cross-correlation, a phase cross-correlation (PCC) based on the similarity of instantaneous phases was proposed by Schimmel [41], where the results showed that the PCC technique is able to detect a weak phase. Recently, Huang et al. [10], [11] took advantage of the cross wavelet transform to determine the relative delay time, and their results showed a better performance over many time series cross-correlation approaches.
In the cross-correlation technique, the size of the correlation window is a critical parameter. Akram and Eaton [42] recommended a window size ≤10 times the dominant period of the signal for large deviations from the true times, otherwise, a window size ≤3 times the dominant period of the signal was suggested. VanDecar and Crosson [37] chose a 2 to 4 second window starting from the preliminary arrival times for the teleseismic events. A combined window (first applied to a large window, then changed to a small window)based cross-correlation was proposed in Raymer et al. [43]. Huang et al. [10] selected a 1.0 s window in the time interval [t − 0.1 t + 0.9] s for a mining microseismic signal cross-correlation, where t is the roughly estimated preliminary P phase arrival time. In this paper, a windowed crosscorrelation (WCC) technique that only utilizes the first half cycle after the P phase arrival is proposed for determining the direct P-phase relative delay time by considering the high frequency concentration property of the first-arriving phase efficiently.

C. HOW TO ADDRESS THE ABOVE PROBLEMS
The main objective of this work is to introduce a high resolution 3-D velocity model into the source location of mining microseismicity, whilst using ray-tracing via the shooting method. Another objective is to propose a WCC technique for a quick, robust and precise determination of the P-phase relative delay time. The remainder of this paper is organized as follows. In section II, we introduce the basic location algorithms, the seismic ray tracing and the WCC technique. In section III, we artificially generate 892 microseismic events to test the importance of the 3-D velocity model in locating mining microseismicity. Its application to eight location pre-measured blasting experiments, comparisons with previous studies and discussions are given in section IV. Finally, in section V, we present brief conclusions of this work and discuss future work.

A. BASIC LOCATION METHODS
The Cartesian coordinates and original time of a microseismic event are denoted as (x 0 , y 0 , z 0 , t 0 ), and the Cartesian coordinates and arrival time measured of the i th station are denoted as (x i , y i , z i , t i ). The theoretical travel time from the source to the i th and j th stations are then written as equations (1) and (2), respectively.
where l i is the ray propagation path from the source to the i th station, and v i is the velocity along the ray propagation path. The ray propagation path is discretized into M cells, where d im is the propagation distance in the m th cell, and s im = 1/v im is the velocity slowness for the m th cell. Then, the difference in the theoretical travel time between the i th and j th sensors can be calculated by subtracting equation (2) from equation (1), i.e., The time difference between the measured and theoretical travel times of the i th and j th stations can be written as: By minimizing the summation of t ij for all pairs of stations, we can obtain the source location by avoiding the estimation of original time of the event, which is called the time difference (TD) location method. Its L1 and L2 norms are usually adopted, which are given by: where K is the number of P-phase arrivals used, (x 0 , y 0 , z 0 ) denote the unknown parameters to be solved, and the symbol ''minimize'' means to find the (x 0 , y 0 , z 0 ) that minimizes the objective function. The trigger time (TT) method is also shown in the text for a better illustration of the discussion part. The trigger time method is slightly different from the time difference method by directly minimizing the residuals between the measured and theoretical arrival times of the station and source. Its L1 and L2 norm formulas are given by: (8) where (x 0 , y 0 , z 0 , t 0 ) means the unknown parameters to be solved. Equations (5)∼(8) can be solved by the minimization methods mentioned in the first paragraph of the Introduction. In this paper, a combined global grid search is adopted to solve these equations. The L2 norm-based location methods can obtain a good location result when the picking precision of the P phase arrival follows an approximate Gaussian distribution. If there is a large picking error, its square misfit will dominate the minimization and may result in a large location error. Thus, the L1-based methods usually have a better noise resilience than that of the L2 norm-based method [44]. Compared with the TT method, the TD method can eliminate the effect of the event's original time, which could further reduce the location error.

B. RAY TRACING
The ray tracing method can be treated as solving the eikonal equation (ray kinematics equation). The eikonal equation represents the travel time τ (x, y, z) for a ray passing through the point (x, y, z) in a medium with velocity model v(x, y, z), and is given by [45]: Specifically, τ (x, y, z) = T 0 represents the wave-front (the spatial points of the wave propagating through a medium at the same instant) at instant time T 0 . The wave is propagated from one wave-front to the next by the ray path, which is perpendicular to the wave-front. In the infinite high frequency approximation, the eikonal equation provides a simplification from wave optics to ray optics. In other words, we can solve the wave seismology using geometric seismology.
Snell's law (sin θ 1 / sin θ 2 = v 1 /v 2 ) in geometric seismology states that the sine ratio of the incidence angle and the refraction angle is equivalent to the ratio of the phase velocities in the two adjacent media, and an extended form of Snell's law can be written as: where θ i is the angle measured from the normal of the boundary, v i is the seismic velocity in the i th media along the ray path between the two medium, and p is the ray parameter. A ray propagation example using the shooting method of ray tracing is shown in Fig. 1, where the 3-D velocity model was updated from Wang et al. [46]. To implement the shooting method of ray tracing, we discretized the 3-D space onto a grid with 10 m spacing, then an initial ray set with a 1 • interval from the source was adopted for the first time search of ray tracing, which followed Snell's law between grids. The two closest rays to the geophone were selected, then we narrowed the ray path region by the bisection method until the distance between the ray path and the geophone was smaller than 0.1 m, which was set as the ray path from the source to the geophone.
It is easy to see that the ray propagation path is not a simple straight line as commonly used for a homogeneous velocity model case, but it is always perpendicular to the wave-front, as the eikonal equation illustrates. There may be a focusing and defocusing phenomena in the wave-front evolution (blue and yellow circles in Fig. 1) and multi ray paths to one point at a time, due to the existence of strong low/high velocity zones. The ray propagation example in the 3-D velocity model demonstrates that the homogeneous velocity model-based ray tracing methods may have large picking errors due to the real complex ray path distribution and strong velocity contrasts in the study region.

C. WINDOWED CROSS-CORRELATION
The cross-correlation technique takes advantage of the similarity between the two selected waveforms, and studies have demonstrated its effectiveness in determining the P phase relative delay time T between two waveforms excited by the same earthquake. T between two earthquake waveforms x 1 (t) and x 2 (t) can be obtained by searching for the maximum of the cross-correlation function defined by [47]: To implement the cross-correlation method, a semiautomatic window selection technique is proposed in the next paragraph. The WCC between the windowed specific waveforms u(x 1 , t) and u(x 2 , t) is defined by: An example of the WCC picking is shown in Fig. 2. First, we manually picked a rough time T 1 before the P-phase arrival and a rough time T 2 after the first peak amplitude (the brown window in Fig. 2a). Then, T 2 is automatically extended to the next zero amplitude time T 3 (the green window). Finally, the windowed waveform u(x, t) corresponds to the time series in the time interval [T 1 , T 3 ]. The WCC lag time series between the first windowed waveform and each windowed waveform are shown in Fig. 2b. There are clear maximum values (the marked points) that correspond to the relative delay time between two selected windowed waveforms. The picking time difference of the WCC method and the manual picking method is only 0.0023 s. This very small difference demonstrates the effectiveness of the WCC method.

III. SYNTHETIC TESTS A. SYNTHETIC DATASET
A synthetic test is a good strategy to verify the performance of a proposed method. We used the 3-D P-wave velocity model (Fig. 3) updated from Wang et al. [46] and the same geophone locations as in section IV. There are remarkably high-V and low-V zones in the mine, and the P-wave velocity can vary from 3,500 m/s to 5,500 m/s, which is a large velocity difference. Thus, the homogeneous velocity model and the 1-D velocity model may result in a large location error. In the synthetic tests, we aimed to generate events distributed across the entire 2D plane of Fig. 4a. Using a spacing of 40 m between events in each direction on the 2D plane, would generate approximately 900 events. We artificially generated more synthetic events in the low velocity zones, to account for possible roof collapses and rock bursts that are more likely to occur within or close to the spatial voids that are represented by at least some of the low velocity zones. This created a total of 892 events (red circles in Fig. 4a). The theoretical travel time is obtained by the ray tracing method mentioned in the 4 th paragraph of part B in section II.
The homogeneous velocity model and straight line ray tracing-based TD1 location method (HV-SL-TD1) and 3-D velocity and straight line ray tracing TD1-based location method (3D-SL-TD1) are selected to validate the 3D-RT-TD1 method. The 3D-SL-TD1 method is designed to take into account the effect of travel time perturbations along a ray by due to 3-D velocity variations, whilst using a straightline ray-path as valid for a homogeneous velocity model. The straight line (SL) travel time is calculated by the following procedure: First, we set the mining zone to 380,800 m ≤ X ≤ 381,700 m, 2,995,500 m ≤ Y ≤ 2,998,500 m, 850 m ≤ Z ≤ 1,300 m. Then, this zone is discretized by 3-D grids with 10 m spacing. For simplicity, we connected the source and geophone using a straight line. The travel distance along the straight line in cube i is marked as l i , then the travel time t i in cube i is l i /V i P , where V i P is the P wave velocity in cube i. Supposing that the straight line passed through M cubes, then the straight line travel time can be calculated

B. SYNTHETIC RESULTS
The location results of the synthetic tests are shown in Fig. 4b in boxplots, where the X-axis displays the P-phase average travel time from the source to geophones. Overall, the smaller the average travel time (usually a better geophone coverage), the better the location accuracy. However, because of the 3-D velocity model, the smaller average travel time may have a larger location error. Although the HV-SL-TD1 method has the same ray propagation path to the 3D-SL-TD1 method, the 3D-SLTD1 method has a much better location result than that of the HV-SL-TD1 method, with average location errors of 26.93 m and 49.92 m, respectively. In addition, in events using the 3-D velocity model, the straight-line ray tracingbased location still has moderate errors. The synthetic   Table 1. The microseismic data is acquired from the Institute of Mine Seismology (IMS) system. There are 28 geophones in this system with a sampling frequency at 6,000 Hz: 12 geophones located at the 930 m level, 12 geophones located at the 1,080 m level, and 4 geophones located at the north of the 1,120 m level. The orebody has a thickness of approximately 6 m with a dip angle of 30 • . Its upper wall is mainly dolomite, while its lower wall is mainly shale and sandstone. There are several main faults and many small faults. The mine has a phosphate production capacity that exceeds 2,000,000 tons/year, and its open stope mining method has left many mined-out regions. Therefore, the geological conditions and mining activities result in a complicated velocity model (Fig. 3), which provides a good dataset to verify the 3D-RT-TD1 method for locating mining microseismicity.

B. WCC FOR THE RELATIVE DELAY TIME
Our former studies have shown that there may be obvious high frequencies (>200 Hz) and stable power frequency noises (50 Hz) in the microseismic waveforms [46], [48], which could affect the manual picking of the P phase arrival time. Therefore, a two passes, four-poles, <200 Hz lowpass Butterworth filter and a cosine function filter of 50 Hz are applied to the selected waveforms, with examples of the filtered waveforms shown in Fig. 5a. Then, the WCC is VOLUME 8, 2020  conducted on the filtered waveforms following the steps shown in part C of section II. The lag times in Fig. 5b show clear maximum values, which again demonstrates the feasibility of the WCC technique in determining the relative delay time for a real dataset.
To further illustrate the picking performance of the WCC, the P-phase arrival times/relative delay times from the manual picking and full waveform cross-correlation (FWCC) picking are also included for comparison. The manual picking method is often chosen as a reference for comparing the performance of alternative picking techniques, since it is believed to be the most stable picking method [48]. For the full waveform selection, we take advantage of the length of full window defined in Huang et al. [10]. However, we only used the first wave packet of the microseismic signal. The occurrence of either single or multi-wave packets is due to the complicated effect of refraction, reflection, attenuation and dispersion during the signal propagation through the 3-D velocity field. The signal envelope, based on the Hilbert transform, was calculated for improved identification of the first wave packet. An approximate start and end time were then manually picked for the first wave packet, which were marked as T 1 and T 2 , respectively. Finally, a time closest to T 2 with an amplitude larger than 10% bounds of the maximum amplitude in the first wave packet is marked as T 3 , and the time interval in [T 1 , T 3 ] is selected as our full waveform window.
A comparison of the cross-correlation derived arrival times (both the WCC and FWCC methods) with the manually picked first arrivals is shown in Fig. 6. It is clear that there is only a small time difference (mean = 3.3 ms and standard deviation = 2.8 ms) between the manually picked first arrivals and the WCC picked relative delay times, which may be due to the unclear P-phase arrival amplitude and the personal subjectivity of the manual picking method. While there is a large time difference (mean = 21.5 ms and standard deviation = 23.2 ms) between the manually picked first arrivals and FWCC picked relative delay times, due to that the FWCC method can be affected by the strong waveform coda (Fig. 5a), which can result in a relatively large error in determining P-phase delay times. The picking examples show the effectiveness of the WCC method over the manual picking method and the FWCC method.

C. LOCATION RESULTS
The 3D-RT-TD1 locations of the first and third experimental blasts are shown in Fig. 7. Source locations were obtained by applying a global grid search strategy (light blue stars in Fig. 7). A first approximate location was obtained using a coarse grid (with 50 m spacing) over the entire search volume. A finer grid (with 10 m spacing) covering ±25 m in each direction around the initial location was used to refine the location estimate. The shooting method follows the instruction of ray tracing in part B of section II. Fig. 7 shows that the ray path tends to propagate along high-speed zones and travel away from low speed zones, which can ensure that it propagates from the source to the geophone in the shortest time.   Fig. 7).
The location results for the eight experimental blasts based on manual picking and WCC picking are listed in Table 2, where the M-TD1 method corresponds to manual picking along with 3-D ray tracing-based TD1 algorithm, and the WCC-TD1, WCC-TD2, WCC-TT1 and WCC-TT2 denote the WCC picking and 3-D ray tracing-based TD1, TD2, TT1 and TT2 algorithms, respectively. All these methods obtained good location results, which demonstrate the quality of the picked data and the effectiveness of the 3-D velocity model in locating mining microseismicity. It is interesting to see that the WCC-TD1 and WCC-TT1 methods had a slightly better location accuracy than that of the M-TD1 method, which demonstrates that the WCC can have a better noise resilience than that of the manual picking. There are relatively large location errors for the 8 th event when using the WCC-TD2 and WCC-TT2 methods, which may be caused by large picking errors and the L1 norm-based methods have a better noise resilience to larger picking errors [3]. The location results of the TD-based result are slightly better than those of the TT-based location when using the same norm, which may be due to the fact that the TT-based methods should estimate the influence of the event's original time. Overall, the WCC-TD1 method has the best location performance with a maximum location error of 32.6 m and an average location error of 26.2 m, which are acceptable errors in mining engineering.
Several researchers have studied these eight blasts. Shang et al. [48] applied the TT1 method to these eight events, where the third event has large location errors for different AIC-based pickings. Li et al. [3] showed that the TD2 (53.6 m), TT2 (53.5 m) and TT1 (52.7 m) methods have an average location error of approximately 50 m for a good P-phase arrival dataset, while the TT1 (40.6 m) and virtual field optimization method (VFOM, 42.0 m) methods have location errors of approximately 40 m. Dong et al. [13] proposed an analytical method to remove the large picking errors, and their results showed that the average location error of the TD2 method and CLMAI method (filtered large errors + TD2 method) are 67.3 m and 44.6 m, respectively. Dong et al. [12] proposed an analytical solution based on the logistic probability density function, which can deduce the effect of the large picking errors. They then applied it to the third, 4 th and 6 th blasts, and their average location error was 47.4 m. VOLUME 8, 2020 FIGURE 7. The 3D-RT-TD1 location for the first and third experimental blasts using the WCC pickings. The coloured background shows the 3-D velocity model projected onto a specific 2-D slice plane, depending on the smallest least square distance to all geophones (just for illustration purposes). The velocity scales are shown on the right. The black triangle represents the geophone projection position, the red star corresponds to the best solution position, the blue star is the pre-measured location, and the small light blue stars show the grids of the search method. The red line is the ray propagation path from the located source to the geophone. The left corner map shows the misfit of the objective function with the coordinates corresponding to the pre-measured locations.
All these studies have adopted a homogeneous velocity model or treated it as an unknown parameter in the location problem. There may be large picking errors for the AIC-based pickings, while the WCC picking can obtain a reliable picking. Although researchers have attempted different location approaches, the average location errors are still larger than 40 m. However, for a good 3-D velocity model, the average location errors are approximately 30 m. The above results demonstrate that the WCC and 3-D velocity model combined method is a good way to improve the source location accuracy.

A. COMPARISONS AMONG PICKING METHODS
The ray tracing theory is only a good approximation for wavefield modelling at a very high-frequency limit. However, the frequency band of a real microseismic signal is practically limited. The high frequency band of a wave field is more likely to be absorbed and more strongly scattered than the low frequency band when the wave propagates through a void. In addition, the P-phase amplitude will become smaller by the effect of wave-front healing [49]. Therefore, the absolute P-phase arrival time can be different from the theoretical arrival time of ray tracing for wave propagation in complicated 3-D media. Compared with the manual picking, the cross-correlation technique, taking advantage of the main peak of the windowed waveforms, is less affected by the P-phase arrival amplitude. Moreover, the semiautomatic picking potential of cross-correlation measurement makes it quicker than the manual picking.
The applicability of cross-correlation is usually limited to scenarios where single phases can be distinctly identified and where compared waveform segments have similar waveforms [33]. Practically speaking, the direct cross-correlation delay time measurements for band-limited signals only provide good estimates around the dominant frequency of the waveforms, and do not capture the differences in the details of the waveform shapes very well. In fact, the differences in the details of the direct P waveform reflect the high-frequency travel-time information of the involved phases. The more complicated situation is that the cross-correlation technique may not extract the correct delay time information when two different phases interfere, which is a common case in complex inhomogeneous media where waveforms disperse due to the presence of strong heterogeneities.
The FWCC is very convenient in semiautomatic picking. However, it will inevitably contain various complex wave field propagation effects, which leads to the travel-time obtained by the FWCC lacking a clear physical meaning. The FWCC will be sensitive to the differences in the amplitude spectra and cannot solve multi-arrival problems reliably [40]. For a mining microseismic signal, there will be strong interfered signal and dispersion in the first wave packet of the FWCC [33]. To improve the picking quality and adaptability for a complicated interfered signal, our WCC technique only uses the first half cycle after the direct P phase arrival. Compared with the FWCC, the WCC has evident high-frequency concentration characteristics, which can ensure a better and more reasonable picking result.

B. IMPROVING THE 3-D VELOCITY MODEL RESOLUTION
The 3-D velocity model adopted in this paper is based on geometric ray tracing tomography under an infinite high frequency approximation assumption, and the inversion uses a manually picked dataset. Although the generalization of the asymptotic ray-based method (e.g., paraxial approximation and corresponding dynamic ray tracing) could, to some degree, account for the finite frequency effects on wave propagation, it still suffers from important shortcomings. First, ray theory cannot be used to describe a complex diffracted wavefield produced by a sharply changed interface or strong scatters. This considerably limits its ability to image structures with a very high velocity contrast, which is possible in mining field applications. Second, the asymptotic approximations of geometric ray theory makes it hard to fully characterize the complexity of various wave propagation phenomena (e.g., P-to-S conversions, reflection on complex interface and interference of different types of seismic phases) in a 3-D medium. Third, the computation of wave amplitudes with asymptotic ray theory is rather difficult in heterogeneous media. This problem is important because an inaccurate estimation of the wave amplitudes may translate into errors in travel time sensitivity kernels for high-resolution tomography, and thus into distortion and artefacts in the inverted models. These will result in significant differences between the ray-based time travel tomography model and the real model, as it lacks a sufficient resolution for lowspeed anomalous bodies [50], which will further affect the location accuracy. Therefore, to improve the resolution of the 3-D model and source location, it is important to develop efficient high precision numerical methods that can model the complete (and broadband) wavefield within complex media [51]- [53]. Such methods can describe all phases that can potentially interfere and contribute to the finitefrequency travel time measurement for the time window of interest [54].

C. INCLUSION OF LOW VELOCITY ZONES
There are various extremely low velocity zones (e.g., tunnels, caves and mined out zones) in a mine that are filled with air and backfill material. Some researchers have considered these voids as a simple 3-D velocity model with regular geological shapes. Collins et al. [23] showed an improvement in the location accuracy when accounting for mining stopes at a hard rock mine. Gharti et al. [25] considered velocities of air, rock and ore in their mining applications. Especially, a more significant improvement can be obtained when seismic events occur near the boundaries of the voids [24]. These studies used a 1-D velocity model or a homogeneous velocity model plus the voids, which are still quite different from a real 3-D velocity model in the ray tracing.
The seismic wave actually propagates around the void instead of through part of the void, and a detailed description is given as follows. For a mining microseismic signal, its main frequency of the P phase arrival is approximately 200 Hz, and the average P wave velocity in a mining zone is approximately 4,000 m/s. Therefore, the wavelength corresponding to the main frequency is approximately 20 m, which is usually similar to or larger than the size of a void. Thus, there will be a strong and complicated wave-front healing effect in these low velocity regions. Therefore, it is not very suitable to use geometric ray tracing based on an infinite high frequency approximation for wave propagation through these extremely low velocity zones. In our study, we treated the void as a low velocity zone, which is a compromise between its surrounding relative high velocity and the void's extremely low velocity (sound speed in the air). The travel time of the ray tracing in an average low velocity zone causes a lower arrival time error than that of the travel time of the ray directly passing through the voids by considering the diffraction of the actual wave propagation through the void zone. Thus, it is reasonable to treat the voids as a low velocity zone for the ray tracing modelling.
To accurately consider the void shapes and low velocity zones (e.g., backfill material), an accurate wave field modelling could be adopted, which can account for the finite frequency nature of the wave propagation and improve the location accuracy. Currently, wavefield modelling still has quite a high computational cost and can involve complex mesh generation: within the air-filled void V P equals the air sound speed (∼340 m/s), which is an extreme value for the numerical modelling that is hardly achievable due to the restriction of numerical dispersion [55], [56]. If the above problems can be solved, the wave field modelling could be a promising way to improve the 3-D model resolution and the accuracy of locating mining microseismicity.

D. COMBINING WITH OTHER EFFECTIVE LOCATION ALGORITHMS
In this paper, a combined global grid search technique is proposed to minimize the location objective function, which can usually obtain a global optimization. However, the grid VOLUME 8, 2020 search technique is very time consuming. Therefore, iterative techniques, as shown in the first paragraph of the introduction, as well as some more advanced methods (e.g., Markov Chain Monte Carlo (MCMC), Transdimensional Bayesian inversion [57], Hamiltonian Monte Carlo (HMC)), could be combined with the 3-D velocity model for a more efficient source location. However, the location results can still be affected by large picking errors. Therefore, a location method that accounts for these picking errors (e.g., the VFOM method in [3]) could be combined with the 3-D velocity model. Alternatively, a technique that wipes out the large picking errors before the location procedure (e.g., Dong et al. [13]) could be used in combination with a 3-D velocity model. The removal of large picking errors can also be achieved by using subsets of the pick data. Such a method could involve, for example, generating data subsets that each have 4 P-phase arrivals. For a dataset with a total of n P-phase arrivals, this would generate c 4 n data subsets. The analytical solutions of equation (8) based on a homogeneous velocity model for data subsets without large picking errors should be similar, while the solutions of equation (8) for data subsets with large picking errors should be very different. We can then use the highly similar solutions (e.g. their average) to decrease the effect of large picking errors and obtain an approximate event location and origin time. The arrivals with large picking errors can then be removed by using a threshold value of the time difference between the theoretically calculated travel time and the observed travel time based on the approximate event location and origin time. After their removal, a more accurate event location can be obtained by locating with the 3-D velocity model. Moreover, more effective and precise waveequation-based location approaches could be also adopted for a large number of source location applications.

VI. CONCLUSION
In this paper, a WCC technique and a high resolution 3-D velocity model were used for locating mining microseismicity. The picking results of the WCC method were similar to those of conventional manual picking, while the WCC required lower picking precision than those of the conventional manual picking method and had a better noise resilience. In addition, the WCC method could effectively remove the effect of the complex waveform spectra for relative delay time determination compared with the FWCC method with a much better picking performance. Then, a high resolution 3-D velocity model was applied to source location in a mining microseismic zone. The synthetic and application results demonstrated that the location accuracy of the 3D-RT-TD1 method was better than that of the HV-SL-TD1 method. The TD-based method was slightly better than the TT-based method when using the same norm, and the L1 norm-based method had a better noise resilience to larger picking errors than that of the L2 norm-based method. The average location errors of traditional homogeneous velocity-model based methods were all larger than 40 m, while the location errors were within 30 m for a simple 3D-RT-TD1 method. Further developments can consider a higher resolution 3-D velocity model in combination with other effective location methods, and include spatial voids (e.g., tunnels and cavities) in the inversion.
YI WANG received the Ph.D. degree from the School of Universe, Université de Toulouse III Paul Sabatier, Toulouse, France, in 2017. He holds a postdoctoral position at Sun Yat-sen University. His current research interests are seismic tomography, mechanical multiparameter inversion, and high-precision mathematical physics modeling techniques. He has served as a Reviewer of IEEE ACCESS and Geophysics. VOLUME 8, 2020 XUEYI SHANG (Member, IEEE) received the Ph.D. degree from the School of Resource and Safety Engineering, Central South University, Changsha, China, in 2019. From 2017 to 2019, he was supported by the CSC to study in the Research School of Earth Sciences, Australia National University, Canberra, Australia. He is currently a Lecturer with Chongqing University. His current research interests are microseismic monitoring and rock mechanics. He has published more than ten peer-reviewed articles and nine invention patents as a Young Researcher, among which one article received the prize of the best article at the University of Seville, Spain. He has served as a Reviewer of some journals, such as IEEE ACCESS, Computers and Geosciences, Soil Dynamics and Earthquake Engineering, and Computers and Electrical Engineering.
KANG PENG received the Ph.D. degree from the School of Resource and Safety Engineering, Central South University, Changsha, China, in 2014. He is currently an Associate Professor with Chongqing University. He has published more than 20 peer-reviewed articles and 20 patents as a Young Researcher. His current research interests are rock mechanics and mining science. He received three projects of National Science and Technology Awards and more than ten provincial and ministerial awards. He has served as a Reviewer of some journals, such as IEEE ACCESS, Transactions of Nonferrous Metals Society of China, and Rock Mechanics and Rock Engineering.