High-Accuracy Location of Microseismic Events in a Strong Inhomogeneous Mining Environment by Optimized Global Full Waveform Inversion

Abstract

High-accuracy determination of a microseismic (MS) location is the core task in MS monitoring. In this study, a 3D multi-scale grid Green’s function database, depending on recording wavefield frequency band for the target mining area, is pre-generated based on the reciprocity theorem and 3D spectral element method (SEM). Then, a multi-scale global grid search strategy is performed based on this pre-stored Green’s function database, which can be effectively and hierarchically processed by searching for the spatial location. Numerical wavefield modeling by SEM effectively overcomes difficulties in traditional and simplified ray tracing modeling, such as difficult wavefield amplitude and multi-path modeling in 3D focusing and defusing velocity regions. In addition, as a key step for broadband waveform simulation, the source-time function estimated from a new data-driven singular value decomposition averaged fractional derivative based wavelet function (DD-SVD-FD wavelet) was proposed to generate high-precision synthetic waveforms for better fitting observed broadband waveform than those by simple and traditional source-time function. Combining these sophisticated processing procedures, a new robust grid search and waveform inversion-based location (GSWI location) approach is integrated. In the synthetic test, we discuss and demonstrate the importance of 3D velocity model accuracy to waveform inversion-based location results for a practical MS monitoring configuration. Furthermore, the average location error of the 3D GSWI location for eight real blasting events is only 15.0 m, which is smaller than error from 3D ray tracing-based location (26.2 m) under the same velocity model. These synthetic and field application investigations prove the crucial role of 3D velocity model, finite-frequency travel-time sensitivity kernel characteristics and accurate numerical 3D broadband wavefield modeling for successful MS location in a strong heterogeneous velocity model that are induced by the presence of ore body, host rocks, complex tunnels, and large excavations.

1. Introduction

Mine microseismic (MS) monitoring can effectively evaluate the risks of a mining area and optimize production plan, in which MS location is a fundamental indicator parameter to study the fracturing location of rock mass, understand crack extension law, focal mechanism inversion, and analyze characteristics of mining activities [1,2,3]. At present, ray tracing-based location methods based on travel time are generally used for MS location, in which the objective function is built by the difference between theoretical and picked arrival time of specified seismic phases. The ray tracing-based location methods have been widely used in the literature due to their simple theory and operation. Detailed reviews of related studies can refer to [4,5]. However, ray tracing-based location methods are set up on the basis of geometrical ray approximation for high-frequency wave propagating, which have inherent difficulties in modeling multi-path effects and focusing and defocusing effects of realistic broadband wavefield propagation in complex velocity structure. Moreover, they could be easily affected by large arrival time picking errors. The above shortcomings increase difficulties in an accurate location with ray tracing-based location technique [6]. In the past few years, with the development of computer technologies and improvement of numerical modeling methods, wave equation-based location methods have been gradually introduced and applied to a MS location in the fields of mining and oil-gas exploration. Wave equation-based modeling can overcome the above shortcomings of simplified ray tracing-based location methods. In addition, the wave equation-based location can directly utilize the recording broadband waveform data, that is, it does not require the high-accuracy arrival time picking of specified seismic phase. Wave equation-based location methods usually fall into two categories, i.e., wavefield reverse-time migration-based location [7] and waveform inversion-based location [8]. The former takes the maximum energy spatial point of stacking backpropagated waveform from all sensors as the location result, and it only requires a small amount of simultaneous wavefield simulations. For example, Witten and Shragge [9,10] developed a full wavefield-based reverse time imaging method for locating microseismic events. The wavefield reverse-time migration-based location methods are conceptually simple and easy to perform, but the critical imaging conditions for location can be easily affected by noises [11]. As the ultimate means of geophysical inversion, waveform inversion-based location methods takes advantage of the fitting degree between synthetic waveforms and observed waveforms, which fully considers waveform complexity and can comprehensively constrain the source location.

1.1. Calculation of the Green’s Function

A high accuracy Green’s function is the basis for a successful waveform inversion (see Section 2.1 for more details). Due to limited computing resources in the early years, the frequency wavenumber (FK) method that adopts a one-dimensional (1D) horizontal layered velocity model was usually used to model Green’s function. Moreover, Tkalčić et al. [12] utilized different 1D velocity models from the source to each station to reduce the influence of velocity models errors. However, in complex structural areas including mines, there is strong three-dimensional (3D) heterogeneity, thus, a 1D velocity model usually cannot accurately model broadband seismic waveforms in such regions. In particular, it performs unsatisfactorily in characterizing the high frequency waveform information.

In recent years, with the improvement of parallel computing capabilities, the Green’s function numerical computation based on 3D velocity model (3D Green’s function) has been efficiently developed. Hsu et al. [13], Zeng and Song [14], Eisner and Clayton [15], Bignardi et al. [16] and Zhao et al. [17] separately modeled 3D Green’s functions by utilizing finite element method (FEM), finite difference method (FDM), boundary element method (BEM) and pseudo spectral method (PSM). However, the FEM is usually based on low-order polynomials that have poor performance in modeling high-frequency signals. The FDM is easy for implementation and has high calculation efficiency, but the FDM is hard to adapt to complex geological structures and irregular boundary conditions for preset numerical accuracy. The BEM can effectively deal with complex boundary problems, but it is often limited to simulate linear wavefield response in relatively simple media. The PSM has a very high calculation precision, but it is also difficult to handle complex regions and interfaces [18]. The spectral element method (SEM), combining the high accuracy of spectral methods and flexibility in grid generation of FEM, is becoming more and more popular in simulating broadband seismic wavefields in 3D complex structures. Komatitsch and Vilotte [19] and Komatitsch et al. [20] developed open source parallel programs (1D~3D SEMs) for modeling seismic waveforms. On this basis, Fichtner et al. [21] and Tape et al. [22] conducted the 3D SEM for imaging large-scale 3D velocity structures. While Liu et al. [23], Hejrani et al. [24] and Hejrani and Tkalčić [25] obtained focal mechanisms and locations of natural earthquakes in large-scale regions through the 3D SEM and full waveform inversion (FWI). In conclusion, 3D SEM-based waveform inversions have been successfully applied in imaging large-scale velocity structures, and inversions of focal mechanisms and source location for natural earthquakes. Inspired by this, we attempted to introduce the 3D SEM into broadband waveform modeling and source location in a mine, which has a smaller spatial scale but with more complex structures.

1.2. Waveform Inversion-Based Location Methods

Waveform inversion-based location methods constantly update synthetic waveforms by changing input source parameter models, and the source location is determined when the synthetic and observed waveforms have the best fitting. As early as 1982, McMechan [26] numerically modeled wavefield propagation in a two-dimensional (2D) uniform-velocity model using the FDM, and he inverted the source location and origin time of MS events. Then, Wu and McMechan [27] inverted source location and magnitude based on a more efficient PSM and approximate linear iterative algorithm. Sen et al. [28] developed a 1D velocity model and waveform fitting-based full moment tensor inversion and source location method, which has been applied to microseismic event location in the fields of mining engineering [29] and gas injection [30]. Michel and Tsvankin [31] simultaneously obtained the source location, origin time and point source moment tensor in a 2D inhomogeneous velocity model by the utilizing waveform inversion method. In addition, Michel and Tsvankin [32] and Igonin and Innanen [33] used waveform inversion to jointly inverse the microseismic event location and velocity model. Wang and Alkhalifah [34] put forward a source function independent waveform inversion to obtain a microseismic source location, origin time, and velocity model. Kim et al. [8], Rodriguez et al. [35], and Hejrani et al. [24,25] performed joint inversion on point source moment tensor and location for natural earthquakes in regional scales under high-resolution 3D velocity models. Multi-parameter joint inversion greatly increases the search dimensions and the nonlinearity of the waveform inversion problem, and the multi-parameters are essentially highly coupled in the inversion process, which make the nonlinear waveform inversion difficult to converge. Therefore, more search times or iteration numbers are needed to find the global optimal location.

Therefore, some studies have simplified the source parameter model into an inversion that only considers source locations (the simplification process and description are shown in Section 2.1), which improves inversion speed and reduces the nonlinearity of convergence. Shekar and Sethi [36] formulated an FWI scheme for determining hydraulic fracturing and mining microseismic event locations, and synthetic test on the 2-D SEG/EAGE overthrust model demonstrates its feasibility. Tong et al. [37] derived a waveform inversion-based location method using the time-domain wave equation and cross-correlation travel time measurement between simulated and observed waveform as the misfit. They tested the proposed location algorithm through 2D and 2.5D (a 2D velocity model is adopted in the vertical plane determined by the source location and sensor location) field and regional velocity models. While Huang et al. [38,39] derived a waveform inversion-based location method that is independent from source-time characteristics, and they successfully applied the fast convergent truncated Newton method for microseismic source location based on 2D and 3D velocity models. Iteration-based optimization technique is often used in the above nonlinear waveform inversion-based location problem, which may be strongly affected by the initial model choice and cannot fully search the whole model space to obtain global convergence. In view of the above difficulties, we adopted the SEM for accurate wavefield modeling and global grid search idea from Hejrani et al. [24] for global source location. The difference is that a multi-scale grid generation and search strategy depending on recording waveform frequency content (the fine grid wavefield only stored near the active mining zones) was used in this study, which greatly reduces the computational time and storage for global grid search.

1.3. Green’s Function Database Generated Based on the Reciprocity Theorem

The general idea of global grid search and waveform inversion-based location (GSWI location) is to calculate the Green’s function from each candidate grid point to each sensor, and then the synthetic waveform can be obtained through the product of Green’s function, source-time function, and point source moment tensor. However, such a large number of wavefield modeling in a 3D space is impossible for existing parallel computing. An accurate and wise choice is to use the reciprocity theorem in elastic dynamics [40] for reducing the computation of Green’s functions. According to the time-space reciprocity of point source Green’s function between source point and sensor point, it only needs to set the sensor locations as virtual sources for wavefield modeling. In this way, the Green’s function of each point in the whole space is generated and stored as a database. The reciprocity theorem can significantly reduce the computation of waveform modeling with tens of thousands of candidate source grid points, and it will be easy to recalculate the Green’s function database for any updated velocity model.

Taking advantage of the reciprocity theorem, many researchers have carried out a series of studies on GSWI location. Bouchon [41] investigated the radiation mode of teleseismic body waves by wavefield modeling combining the reciprocity theorem and flat layer velocity model. Graves and Wald [42] generated a Green’s function database using the reciprocity theorem and found that the 3D velocity model can more accurately invert finite fault sources than the 1D velocity model. Eisner and Clayton [43,44] used the reciprocity theorem to simulate the response of 300 earthquakes to Southern California faults along with site and complicated wavefield propagation path effects in complex heterogeneous media. In some other studies, focal mechanisms, source location, and velocity model are jointly studied by utilizing the reciprocity theorem. Lee et al. [45] and Lu and Zhou [46] carried out lowpass filtered waveform inversion to focal mechanism solutions of Southern California and Lushan earthquakes, respectively, in which the 3D Green’s function database is generated by the reciprocity theorem. By utilizing the reciprocity theorem and grid search algorithm, Hejrani et al. [24] and Hejrani and Tkalčić [25] conducted high-resolution waveform inversion on focal mechanism solutions and locations of earthquakes in Papua New Guinea, the Solomon Islands, and Australia. The results show that the waveform inversion results based on 3D velocity models are improved compared with those obtained by 1D velocity model. The above source parameter studies almost focus on natural earthquakes at regional scales, and there are few applications in broadband waveform modeling and source parameter inversion in a local scale (e.g., in the mine region) with more complex tectonic settings. With the improvement of computation capability and quality of broadband seismic records, it has become an inevitable trend to comprehensively carry out the 3D GSWI location study in a mine.

1.4. What Will Be Objectives of This Work

In this study, we introduced 3D waveform inversion for mine microseismic event location and applied a multi-scale grid search algorithm based on global optimization into source location problem in a complex 3D local region (abbreviated as MS-GSWI location). Our main contributions and innovations are shown as follows:

(1)

The SEM was used to effectively generate high-accurate 3D Green’s function in a mine region with complex structures, which overcomes the difficulties of modeling multi-path effects and focusing and defocusing effects in ray tracing methods.

(2)

A 3D Green’s function database in the mine was efficiently generated by the reciprocity theorem and served as multi-scale grid-search strategy, which largely improves the computational efficiency of 3D GSWI location.

(3)

Both the data-driven SVD-averaged fractional derivative based wavelet function (DD-SVD-FD wavelet) obtained by windowed specified seismic phases and traditional fractional-order Gaussian wavelet can be used for estimating the source-time function to generate synthetic waveforms for waveform inversion well. By comparing them, waveforms synthesized by the DD-SVD-FD wavelet have a higher accuracy and more reasonable physical explanation.

(4)

Synthetic test results show that the accuracy of velocity model is important for waveform inversion-based location, and the average location error of eight blasting events is only 15.0 m. The location error is smaller than the 3D ray tracing-based location method (average location error 26.2 m) under the same 3D velocity model.

The rest of this research is arranged as follows: Section 2 introduces main theories including synthetic waveform modeling based on the reciprocity theorem, estimation of source-time function with DD-SVD-FD wavelet, waveform modeling considerations, and multi-scale grid search strategy. In Section 3, dependence of GSWI location results to the accuracy of velocity model is tested by using velocity models with different smoothing degrees. Section 4 further tests effectiveness of the 3D proposed MS-GSWI location method by eight blasting events with known locations and compares it with previous published location results. In Section 5, the critical issues of 3D Green’s function database modeling, estimation of source-time function, and effects of 3D velocity model on location accuracy are discussed. Section 6 makes summaries and puts forward prospects for future work.

2. Methodology

2.1. Synthetic Waveforms Based on the Reciprocity Theorem

The core mission of waveform inversion is to find synthetic waveforms that fit the recording waveforms as much as possible. For a mine MS location, it is usually reasonable to treat the event as a point-source, and the synthetic waveform from earthquake location $x_s$ to sensor location $x_r$ is set as $u(x_r;x_s)$. Solving the elastic wave equation is a precise way to simulate three-component synthetic waveform, however, it is quite hard to obtain a high-resolution elastic parameter model for our study zone (limited by sparse single-component sensor network and low signal-to-noise ratio records). Therefore, in this study, results are all based on the acoustic wave equation. The synthetic waveform data can be simplified and calculated by the following equation:

$$u(x_r;x_s)=\left[G(x_r;x_s)*S^{\prime }(t-t_0)\right]M(x_s)$$

(1)

where, the Green’s function $G(x_r;x_s)$ refers to the time-space point pulse displacement response excited at the earthquake location $x_s$ on the sensor location $x_r$, which is determined by $x_s$, $x_r$ and propagation media distribution; $S^{\prime }\left(t\right)$ and $t_0$ represent the normalized source-time function and origin time of a MS event, respectively; $M\left(x_s\right)$ characterizes the influences of source amplitude and focal mechanism.

Complicated focal mechanism (such as double-couple, explosion) could significantly affect the waveform through source-time function form. The double-couple component, in contrast to the non-double-couple components usually reflecting collapse or implosive events, is an important indicator of shear fracture for a microseismic event. Its impact on waveform records is complicated. However, the influence of the double-couple component is generally better reflected and quantified in the waveform recorded by multi-component monitoring network with a higher signal-to-noise ratio through the complex polarity distribution character. Considering that it is usually difficult to constrain the double-couple component under a mining monitoring network with sparse single-component sensors and low signal-to-noise ratio records, blasting events were selected for synthetic and application tests of the GSWI location method in the study, and so the scalar acoustic wave approximation is acceptable. Therefore, the source amplitude $M\left(x_s\right)$ can be naturally incorporated into the source-time function S(t), thus, the source-time function characterizes moment release rate, strength property, and rupture evolution at the source location. The origin time $t_0$ of an MS event can be regarded as an overall time-shift of all the recording waveforms. It can be eliminated in the data preprocessing step [11], then the $M\left(x_s\right)S^{\prime }(t-t_0)$ is shortened as $S\left(t\right)$ for convenience and Equation (1) can be rewritten as Equation (2). The source-time function estimation strategy will be defined in Section 2.2. After determining the source-time function, our focus moves to an accurate and efficient modeling of the Green’s function $G(x_r;x_s)$.

$$u(x_r;x_s)=G(x_r;x_s)*S(t)$$

(2)

Assuming that the wavefield modeling is carried out at a preset highest frequency of 100 Hz, and the average P wave velocity is set as 4500 m/s, then the shortest wavelength of modeled wavefield will be 45 m. According to the sampling conditions of accurate wavefield simulation described in Section 2.3: for an accurate SEM numerical simulation, each wavelength should contain at least 4.5 grids [47], that is, the grid length should be less than 10 m. By referring to the mining zone of the Yongshaba mine (details in Section 4.1), the dimension is about 3000 × 1000 × 500 m³, then the number of grids reaches 300 × 100 × 50 = 1,500,000. If the recording waveforms are fitted through forward modeling waveforms, the Green’s function at each candidate source grid point will be modeled in the GSWI with relatively large computation.

To solve this computational problem, it is necessary to store the Green’s function in wavefield modeling as the database. Only a small number of wavefield-forward modeling is carried out based on the reciprocity theorem, that is, the number of wavefield modeling is directly proportional to the number of sensors instead of all searching grid points as candidate source. There are usually 10~50 sensors in a mine MS monitoring system, thus, it is easy to obtain the Green’s function from each sensor as virtual source to each grid point. In accordance to the reciprocity theorem in scalar acoustic wave representation, Equation (2) can be rewritten as:

$$u(x_r;x_s)=G(x_s;x_r)*S(t)$$

(3)

where, $G(x_s;x_r)$ denotes the Green’s function from the sensor location $x_r$ to source location $x_s$. In a practical application, coarse grids are used in the regions far away from the mining area, and fine grids are used in the regions near the active mining area (where MS and blasting events occur frequently in practice). The combination of coarse and fine grids greatly reduce modeling computation and storage space of Green’s function database, which is called the multi-scale strategy. It is also conducive to improve the grid-search speed in the GSWI process.

For the broadband waveform inversion-based location method, it needs to define a reasonable objective function to determine the fitting degree of the target observed waveform and theoretical modeling waveform. After that, inversion algorithm searches for the minimum value corresponding point as the source location. Due to the complex geological structures and velocity distribution in a mining area, the inversion problem will be highly nonlinear if the L2 norm misfit between the synthetic and observed waveforms is directly used as the objective function. This research utilized the L2 norm difference of cross-correlation travel time as the objective function of GSWI location study, which is more robust and linear for inversion problems. The cross-correlation function of two discrete waveform time series is defined as follows:

$$f(x_r,u,d)={\displaystyle \sum _{t=-\infty }^{\infty }u(x_r,t)\cdot d(x_r,t+\tau )}$$

(4)

The cross-correlation travel time difference is determined by the maximum value of Equation (4) based on the first-order approximation, that is, $\Delta \tau (x_r,u,d)=\frac{1}{N_r}{\displaystyle {\int }_0^{T}w(t)}\frac{\partial u(x_r,t)}{\partial t}[d(x_r,t)-u(x_r,t)]dt$ where, $N_r={\displaystyle {\int }_0^{T}w(t)}u(x_r,t)\frac{\partial u(x_r,t)}{\partial t}dt$, $w(t)$ is the time window function of specified seismic phases in the interval [0 T]. Through the summation of cross-correlation travel time difference $\Delta \tau (x_r,u,d)$ of all windowed waveforms, the target misfit function for waveform inversion location is shown as follows:

$$\chi (u,d)=\frac{1}{2}{{\displaystyle \sum _{i=1}^M\left[\Delta \tau (x_r,u,d)\right]}}^2$$

(5)

where, M represents the number of windowed waveforms for location.

Compared with the ray tracing methods where one can generate a distribution of locations assuming some uncertainty on the arrival picks, the waveform-based location method uses wavefield information with limited frequency band distribution. Many factors will affect the uncertainty of the final location result in a coupled way, such as the recording clock deviation between different sensors, error of velocity model, and spatial distribution of sensor array with respect to microseismic events. In this study, the spatial distribution characteristics of the misfit function around the GSWI location point are used to qualitatively evaluate the uncertainty of location parameters. More sophisticated quantitative uncertainty analysis for waveform-based location is based on second-order Hessian matrix with respect to a target misfit function, which will be an independent future work.

2.2. Source-Time Function Estimation

The source-time function $S\left(t\right)$ absorbing origin time and event magnitude is usually estimated by the Gaussian and Ricker wavelet [8,37], which is calculated as follows:

$$S(t)=-\sqrt{\pi }{\omega }_0\exp \left(-\frac{{\omega }_0^2}{4}{(t-{\tau }_0)}^2\right)$$

(6)

where, ${\tau }_0$ is the time value corresponding to wavelet central time; ${\omega }_0$ denotes the controlled dominant frequency, which is related with specified sensor frequency response.

While Wang [48] proposed a kind of Gaussian family wavelet based on fractional-order derivatives, whose frequency domain form is expressed as follows:

$${\mathsf{\Phi }}^{(u)}(\omega )={\left(\frac{u}{2}\right)}^{-u/2}\frac{{\omega }^u}{{\omega }_0^u}\exp \left(-\frac{{\omega }^2}{{\omega }_0^2}+\frac{u}{2}\right)\times \exp \left(-i\omega {\tau }_0+i\pi \left(1+\frac{u}{2}\right)\right)$$

(7)

where, $u$ stands for the fractional-order of time differential. The fractional-order Gaussian function can change fractional order $u$ and $w_0$ at the same time, thus, it has more freedom than the pure Gaussian or Ricker wavelet. However, for quite dissimilar waveform data, there are still some limitations in wavefield modeling using the two parameters ($u$, $w_0$) fractional-order Gaussian function as source-time function. In order to further improve fitting performance for complex waveforms, we put forward a new approach to obtain the source-time function based on observed waveforms themselves, namely the DD-SVD-FD wavelet estimation. The detailed steps are shown as follows: firstly, each recording waveform containing the first half period of P-phase arrival is cut by the semiautomatic windowing method proposed by Wang et al. [6]: (1) A two passes, four-poles, 200 Hz lowpass Butterworth filter was conducted to the recorded signal. (2) We manually picked a rough time T₁ before the P-phase arrival and a rough time T₂ after the first peak amplitude; (3) The T₂ is automatically extended to the next zero amplitude time T₃; (4) The windowed waveform corresponds to the filtered time series in the time interval [T₁, T₃]. The waveforms are aligned based on the cross correlation time difference between different windowed waveforms. Then, the SVD average of windowed and aligned signals is taken as a basis function, where SVD is a classical signal processing approach based on coherence between different data traces. It can enhance the signal-to-noise ratio and suppress incoherent background noise [49]. The aligned and windowed data traces can be represented by a data matrix $S{(x_r,t_i)}_{M\times N}$ (dimension is M × N), with M traces and N time samples per trace (usually M < N). Assuming that the data matrix $S{(x_r,t_i)}_{M\times N}$ has a rank R < M, the SVD of $S{(x_r,t_i)}_{M\times N}$ is a orthogonal expansion of the data space:

$$S{(x_r,t_i)}_{M\times N}={\displaystyle \sum _{i=1}^R{\sigma }_i{u}_i{v}_i{}^{\mathrm{T}}}$$

(8)

where, ${\sigma }_i$ is the ith nonzero singular value, and $u_i$ and $v_i$ are the corresponding left and right singular eigenvectors with dimensions M and N, respectively. The product ${\sigma }_i{u}_i{v}_i{}^{\mathrm{T}}$ is a (M × N) matrix with rank one. If the windowed dataset has a high degree of trace-to-trace coherence, then the most coherent or average trace can be constructed from the largest singular eigenvector with the largest eigenvalue. The SVD acts as a data-driven, lowpass filter by keeping the main character of the dataset and rejecting highly incoherent noises.

Finally, time scaling factor and fractional-order are introduced to estimate the frequency-time domain characteristics and obtain the source-time function in the frequency-time domain as (Equation (9)).

$$\{\begin{array}{l}\tilde{S}(w)={\displaystyle {\int }_{-\infty }^{+\infty }S}(t)e^{-iwt}\mathrm{d}t\\ \tilde{S}(at)=\tilde{S}(w/a)/a\end{array}$$

(9)

where, $S(x_r,t_i)$ indicates the denoised time sequence of the rth windowed and aligned waveform at the time $t_i$ (r = 1, 2, ..., M, and i = k, k + 1, ..., k + N − 1), k represents the start point of the windowed waveform. $\tilde{S}(w)$ and $\tilde{S}(at)$ separately stand for the Fourier transforms of S(t) and $S\left(at\right)$, in which a is a scaling factor in time domain. The optimal source-time function can be estimated by updating the two parameters [u, a] searching for the maximum value of normalized cross-correlation coefficient between recording waveform and synthetic waveform modeled from corresponding source-time function.

2.3. Wavefield Modeling Based on Acoustic Wave Equation

Tong et al. [37] and Huang et al. [38] derived wavefield modeling and waveform inversion-based location methods, by taking advantage of numerically solving the acoustic wave equations in time domain and frequency domain, respectively, and this study will carry out wavefield modeling under the same framework. The scalar acoustic wave equation in frequency domain can be expressed as follows:

$$-\frac{{\omega }^2}{c^2\left(\mathbf{x}\right)}u(\mathbf{x},\omega )={\nabla }^{2}u(\mathbf{x},\omega )+S(\omega )\delta (\mathbf{x}-{\mathbf{x}}_s)$$

(10)

where, $c\left(\mathbf{x}\right)$ indicates the compressional P wave velocity at point $\mathbf{x}$ in the medium; $S\left(\omega \right)$ represents the frequency spectrum of the source-time function $S\left(t\right)$, usually showing complex spectrum characteristics.

Through a series of derivations (see [37] for more details), sensitivity kernels with respect to velocity model parameters for wavefield disturbance $\delta u$ at a sensor caused by the source spatial location disturbance are obtained by using the popular adjoint method. On this basis, in accordance with the objective function defined by waveform cross-correlation travel time difference and sensitivity kernels, $\mathsf{\Delta }{T}_i\left({\stackrel{\rightharpoonup }{X}}_r\right)$ and disturbance of model parameters in the wavefield propagation space (especially in first Fresnel zone surrounding the main energy ray path) can be obtained.

$$\mathsf{\Delta }{T}_i(x_r)={\displaystyle {\int }_{\mathsf{\Omega }}K(x_i,x_r,x_s)}\frac{\delta c\left(x\right)}{c\left(x\right)}\mathrm{d}x$$

(11)

where, K is the cross-correlation travel time sensitivity kernel and $\delta c\left(x\right)$ represents the P wave velocity perturbation at point $\mathbf{x}$.

$$K(x_i,x_r,x_s)=-{\displaystyle {\int }_0^T[2c^2(x)\nabla u^{+}(x,t^2;x_r)\cdot \nabla u(x,t^2;x_r)]}\mathrm{d}t$$

(12)

where, $u^{+}(x,t^2;x_r)$ is the adjoint wavefield in the media excited by a virtual adjoint source defined at the sensor $x_r$ (the detailed derivation refers to [37]). The sensitivity kernel of travel time has very complex intensity and spatial distribution, which reflects the frequency-dependent weighted influence of velocity perturbation in the Fresnel zone on the final travel time anomaly $\mathsf{\Delta }{T}_i\left(x_r\right)$ recorded by sensor. It reveals the nature of finite frequency effects in actual wavefield propagation.

2.4. Numerical Wavefield Modeling Conditions

The accuracy of wavefield modeling generally needs to meet two key numerical conditions: the Courant–Friedrichs–Lewy (CFL) condition and the least sampling number of grid points for shortest wavelength in wavefield modeling. The CFL condition satisfies the following inequality:

$$\frac{\mathsf{\Delta }t}{\underset{{\mathsf{\Omega }}_e}{\min }\left(h_e/v\right)}\le C$$

(13)

where, $v$, $h$ and $\mathsf{\Delta }t$ indicate the velocity model, grid size and time step for numerical modeling, respectively; $C$ is the Courant number, generally ranging from 0.3 to 0.5.

Seriani and Priolo [47] pointed out that for a high accurate modeling by SEM, the number of grid points $N_e$ in the shortest simulated wavelength should meet the following condition:

$$N_e=p_{se}\cdot \underset{{\mathsf{\Omega }}_e}{\min }\left(\frac{v}{{\overline{h}}_e}\right)\approx N_{emp}$$

(14)

where, $p_{se}$ is the shortest period of frequency band based on target modeled wavefield, ${\overline{h}}_e$ is the average grid size. $N_{emp}$ represents the number of grid points in a single wavelength that can maintain the numerical simulation accuracy, and it is fixed at 5 in the SPECFEM3D program used for this study. As shown in Equation (13), the smaller the shortest period of modeled wavefield, the smaller the required grid dimension.

Based on the above conditions, we know that after determining the highest frequency of modeled wavefield, the average grid size can be determined through Equation (14) and time step for modeling is then limited by Equation (13). Equation (14) shows that when the highest frequency of modeled wavefield is doubled, the grid size in each dimension will be reduced to half of the original one, and the number of 3D grid points will become eight times of the original number. Equation (13) indicates that as the grid dimension is reduced by half, the time step for modeling that maintains the stability of the time scheme will also be decreased by half. Therefore, when the highest frequency of modeled wavefield is doubled, the wavefield calculation will be 16 times the original one. This is why we propose the multi-scale strategy (also corresponding to the multi-grid operation depending on wavefield frequency content) to reduce the storage requirements and the computational cost.

3. Synthetic Test

To test the proposed 3D GSWI location method, the 3D velocity model inverted through ray tracing tomography based on P wave travel-time data by Wang et al. [50] and shown in Section 4.1 was selected as the target velocity model. Moreover, velocity models with different smoothing scales were used to verify the importance of velocity model accuracy to the 3D GSWI location method. In this method, a 3D spatial Gaussian function has been used for smoothing the original high-resolution velocity model, and the standard deviation is used to define the length of smoothing scale. This different smoothing length approximately reflects the different resolution limit of the ray-based tomographic model obtained from transmission wave data with various signal frequency content [51], while high and low velocity anomalies that are incorrectly positioned may result in a large location error are not considered in this paper. That is, we only conducted a specific qualitative assessment to consider the influence of insufficient resolution of inverted velocity model on GSWI location errors. Velocity model 2D slices under different smoothing scales are shown in Figure 1a. The synthetic test was carried out using the premeasured locations (locations measured before the blasting events) of eight practical blasting events and their corresponding sensors (Table 1 in [6]), and the waveform dataset was generated by the unsmoothed target 3D velocity model with 120 Hz dominant frequency referring to the data spectrum analysis.

The location of blasting event 8 is taken as the example and its waveform inversion-based location results are demonstrated in Figure 1. The figure shows that with the increase of smoothing lengths, the smoothed velocity models gradually tend to be homogeneously distributed, which are characterized by the ray paths from the source to each sensor that gradually become straight lines. The 17.5 m smoothing obtained a 0 m location error, as when a relatively accurate velocity model is used for location, the waveform information related to source location could still be reasonably explored and inverted by the grid search and waveform inversion-based location method. When the model error increases (the smoothing length continues to increase), the source information contained in the data cannot be inverted or stacked correctly, and the location error introduced by velocity error appears. Location errors caused by reduction of velocity model accuracy do not follow a monotonous change law, due to location result is also being controlled by some other factors, such as specific triggered sensor geometric distribution and inhomogeneity of velocity structures near the source location. Furthermore, the overall location error increases with the decrease of velocity model accuracy (Figure 2), which confirms the necessity of an accurate 3D velocity model to waveform inversion-based location in regions with complex geological structural property. In addition, the errors of GSWI location results in this synthetic test with a small number of test events are all within 40 m. It is complicated to perform a completely quantitative analysis explaining the coupling effect between the location error and the current velocity model inaccuracy (such as by smoothing). The systematic location error analysis will be left to the future.

4. Application

4.1. Engineering Background and Multi-Scale Grid Generation

In this study, eight blasting events obtained by the MS monitoring system (Figure 3) in the Yongshaba mine (China) were selected for application test. In this system, there are 28 sensors: sensors T1 and T2 are triaxial sensors, sensors 1–26 are uniaxial sensors, and they are set at the 930 m level (12 sensors), 1080 m level (12 sensors), and 1120 m level (4 sensors) along the mining tunnels. The main working principles of the Institute of Mine Seismology (IMS) are shown as follows: firstly, a sensor can monitor the signal all the time. Then, the continuous signal is transported to its administering NetADC (network analogue-to-digital converter) through 4 core signal cable (uniaxial sensor) and 8 core signal cables (triaxial sensor), and the NetADC samples the signal 6000 times per second (it means that the sampling frequency of a microseismic signal that we used is 6000 Hz). After that, the NetSP (network seismological processor) does some pre-processing to the discrete signal, which comes from its administrated NetADC, including the detection of MS event signals based on threshold value. Finally, the detected discrete signal is transferred to the data exchange center, and then the data processing center. Based on microseismic signals, we can calculate microseismic event location, magnitude, focal mechanism, etc.

The dip angle of ore body in the mine was about 30°. The hanging wall of ore body is dolomite with a high stability, while crushed shale and sandstone were presented in the footwall, the rock mass of hanging wall and footwall have large differences in P-wave velocity [54,55]. After years of mining, a large number of low-velocity regions and high-velocity regions have been formed. Therefore, MS events including fault slip and caving collapse occurred frequently. It can be seen that there are high- and low-velocity anomalies, and the velocity difference can reach as much as 2000 m/s (Figure 4a). Therefore, it is very significant and necessary to carry out 3D waveform inversion-based MS location. Considering that a mine MS event is more likely to occur in the regions with large velocity contrast, the modeling space for mining area was divided into coarse-grid and fine-grid modeling regions according to the degree of velocity model variation and field mining data (dotted lines in Figure 4b show the boundary of the fine-grid range). The coarse grids with a dimension of 20 × 20 × 20 m³ are outside the dashed lines, while the fine grids with a dimension of 10 × 10 × 10 m³ were limited within the dashed lines, so as to reduce storage and communication time of the Green’s function database. If we want to further model the wavefield database with higher frequency bands to improve location precision, finer grids can be refined in the small region near the coarse grid-based location point for new-stage MS-GSWI location. The multi-scale grid generation introduced above combining with the multi-frequency-band grid search strategy, as described in Section 2.4, can gradually optimize the location results and keep a reasonable amount of Green’s function storage and computational cost.

4.2. Case of Waveform Modeling and Source-Time Function Estimating

Taking the black star in Figure 4 as the source location and setting the source as an explosion source, and considering that the main frequency of a blasting event is usually within 250 Hz and the dominant frequency is usually within 120 Hz, a Ricker wavelet with dominant frequency 120 Hz is taken as the source-time function for SEM modeling. The time interval between adjacent wavefronts shown in Figure 4 is 0.004 s. It can be seen that the SEM method can precisely model the focusing and defocusing effects in complicated velocity regions compared with the ray tracing technique (Figure 1 in [6]): The low-velocity regions have focusing effects on the transmitted waves and produce diffraction and multi-path effects. There are defocusing effects on waves propagating through the high-velocity regions, which are usually difficult to rapidly and appropriately calculate travel-time field of waves passing through the strong high-velocity regions by using the ray-shooting method. In addition, the wavefield attenuates rapidly with the increase of propagation distance, indicating that the SEM modeling can correctly consider the variation of waveform amplitude so as to make full use of information contained in the observed waveforms.

To test location performance, eight blasting events with pre-measured locations were selected for application test, and the event location information is listed in Table 1 of [6]. The blasting event two was taken as the example to illustrate the source-time function estimation performance when using the DD-SVD-FD wavelet. After estimating the source-time function by the DD-SVD-FD wavelet and fractional-order Gaussian function, the corresponding modeled waveforms are shown as red and blue curves in Figure 5a, respectively. When estimating source-time function with the fractional-order Gaussian function, the optimal fractional-order $u$ and dominant frequency $w_0$ are obtained by a grid search procedure for maximum cross-correlation coefficient between data waveform and modeled waveform (Equation (7)). Source-time function estimation using the DD-SVD-FD wavelet is achieved by grid search of fractional order $u$ and time scaling factor a. The normalized amplitude spectra corresponding to the colored waveforms are shown in Figure 5b, in which the vertical dotted line indicates the peak frequency, and the shaded region represents the half-bandwidth frequency distribution range. Moreover, the intersection point between the right side of shaded range and curve of amplitude spectrum corresponds to the highest cut-off frequency w_h used for wavefield modeling and GSWI location. It can be seen that both synthetic waveforms obtained by the DD-SVD-FD source-time wavelet and fractional-order Gaussian source-time function can reasonably fit the observed waveforms. After carefully comparing the two results, the synthetic waveforms obtained by the DD-SVD-FD wavelet obtain a better fitting degree to the observed waveforms in the time domain and frequency domain through the normalized cross-correlation coefficient improvement. Generally, some small cycle skipping phenomena may occur in the synthetic waveform simulation with inaccurate source-time function for certain frequency bands, which in turn affected the location results. Therefore, a more precise source-time function estimation leads to a better matching between synthetic waveforms and observation waveforms, which can be used to more effectively constrain the source location and is beneficial to the study of microseismic activity.

4.3. MS Location Examples

Taking the blasting events 1 and 3 as the cases to illustrate 3D MS-GSWI location method, their location results are shown in Figure 6 and Figure 7, respectively. Figure 6b–d and Figure 7b–d also demonstrate that velocity parameter sensitivity kernels of cross-correlation travel time anomaly $\mathsf{\Delta }{T}_i\left(x_r\right)$ have broader and more reasonable illumination range compared with the simple infinite-frequency assumption-based geometric ray tracing methods (Figure 6b and Figure 7b). The sensitivity kernels of cross-correlation travel time anomaly are closely related to the frequency content of the wavefield, and here only the distribution of sensitivity kernels calculated by the Gaussian source-time function with main frequency of 120 Hz is displayed. If the frequency bands tend to infinite high frequency, the sensitivity kernels of travel time anomaly will shrink to a very narrow tube just like the geometric ray. The sensitivity kernels of travel time based on numerical waveform modeling can better reveal multi-path effects around strong low-velocity anomalies (Figure 6c,d and Figure 7c) and evident wavefield interference influence (Figure 7d). In strongly heterogeneous media, the changes of sensitivity kernels of travel time corresponding to the changes of source location will also present very sophisticated characteristics.

The whole velocity field perturbation in the Fresnel zone is weighted by amplitude of sensitivity kernels of travel time anomaly (the red to blue color scale shows the intensity value of positive and negative strength of the sensitive kernel), which finally affects waveforms at sensors and thus leads to complex frequency-dependent perturbation of arrival time and waveform interference. Both the SEM modeled waveforms based on MS-GSWI location result and premeasured true source location can fit the observed waveforms well, with an overall correlation coefficient about 0.8. It should be noted that an appropriate band-pass filtering has been conducted to the waveforms modeled for different sensors to better fit the corresponding observed waveforms. In addition, the total cross-correlation travel time misfit of the modeled waveforms excited from MS-GSWI source location is smaller than that of the premeasured source location, which may be caused by currently inaccurate velocity models. According to misfit distribution of the L2 norm difference of cross-correlation travel time (inserted map in the upper left corner of Figure 6b and Figure 7b), the location determined by the MS-GSWI location is very close to the premeasured source location, and the location errors of blasting events 1 and 3 are 25.0 m and 14.0 m, respectively. The direction-dependent gradient and curvature characteristics of misfit function near MS-GSWI location point could indicate whether the location result along corresponding direction is well constrained.

4.4. MS Location Results

MS-GSWI location results of the eight blasting events are shown in

Table 1. Location errors of the eight premeasured blasting events using different location methods.

Event id	Event Location (m)			Location Error (m)
	X	Y	Z	Shang et al. [56]	Dong et al. [57]	Li et al. [59]	Wang et al. [6]	Wang et al. [7]	GSWI Location
1	381,683	2,997,760	1107	84.8	61.1	31.1	32.6	23.0	25.0
2	381,653	2,997,405	1099				22.4	21.1	14.8
3	381,194	2,996,224	1014	27.6	28.7	49.4	14.3	5.1	14.0
4	381,684	2,997,777	1107	3878.2		30.8	29.8	24.2	15.6
5	381,503	2,997,036	1028				28.5	7.2	7.8
6	381,590	2,997,278	1053	70.1	52.4	49.0	30.8	23.8	21.5
7	381,526	2,997,584	1044	55.0		48.8	24.5	18.2	6.7
8	381,442	2,998,029	1017	143.5		42.8	26.6	13.8	14.6
Average				709.9	47.4	42.0	26.2	17.0	15.0

. The table demonstrates that the maximum location error of the waveform inversion-based location method is 25.0 m and the average error is 15.0 m. There are some uniform velocity model-based location results on the eight blasting events [56,57,58,59], and their average location error is larger than 39.5 m. Particularly, Shang et al. [56] show that location error of the blasting event three in their paper is very large when using automatically picked data. Dong et al. [57] found that the average location error reaches 47.4 m for their selected three events. Li et al. [59] developed a new location method called VFOM to reduce the influence of large picking error, but still obtained an average error larger than 39.5 m. The average error obtained by Wang et al. [6] using 3D ray tracing-based location method is reduced to 26.2 m, which confirms the importance of velocity model accuracy to a mine MS location. Furthermore, Wang et al. [7] utilized a 3D Gaussian beam-based reverse-time migration location method, which has an average location error of only 17.0 m.

The research [56,57,58,59] used homogeneous velocity models and straight ray paths to locate an event, while Wang et al. [6] used a 3D ray tracing-based shooting method. For the shooting method, we discretized the 3D 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, and the two closest rays to the sensor were selected, then we narrowed the ray path region by the bisection method until the distance between the ray path and the sensor was smaller than 0.1 m, which was set as the ray path from the source to the sensor. Furthermore, there are some other ray tracing methods and codes/software packages, such as Lee and Lahr [60] published a 1D velocity model-based location method-HYP071 (software package is available through http://pubs.usgs.gov/of/1999/ofr-99-0023/), and its updated vision that considers different phase arrival time is also available (HYPOSAT, ftp://ftp.norsar.no/pub/outgoing/johannes/hyposat/); Rawlinson and Sambridge [61] developed a multi-stage fast marching method (FMM) for 3D travel time calculation (FMTOMO, http://rses.anu.edu.au/~nick/fmtomo.html), and it can be used for source location and tomography; another widely used ray tracing location method is the double difference approach proposed by Waldhauser and Ellsworth [62] (hypoDD, http://geophysics.eas.gatech.edu/people/cwu/teaching/hypoDD/hypoDD.html#part1_2), they assumed that the propagation paths of wavefields generated by two adjacent seismic events are similar, which effectively reduces the influences of structural anomalies due to the similar ray path from receiver to adjacent earthquakes.

The location results in this study are similar to those obtained by the 3D Gaussian beam-based method, proving that the location method considering finite-frequency effects of wave propagation is significantly superior to the geometric ray tracing-based location. According to Wang et al. [50] and Virieux and Operto [51], the velocity model used in this paper has a resolution of about 50 m. In addition, the Gaussian beam-based method approximately solves the wave equations in a local region around the ray path, thus, partial ray tracing operation is still needed. The modeling accuracy of the Gaussian beam method is lower than that of numerically solving the complete wave equation. If finite-frequency effects are considered by using challenging full waveform inversion, a 3D velocity model with higher resolution could be obtained, which will result in better location results using the proposed MS-GSWI location method.

5. Discussions

5.1. 3D Green’s Function Database Modeling

In waveform inversion, the 3D Green’s function database modeling for a single source point with the finest grid (10 × 10 × 10 m³) is about 120 cpu-hours. Especially when the full-space grid search strategy is adopted, the workload of modeling the Green’s functions from all candidate source grid points to sensors is directly proportional to the total number of spatial grid points. Therefore, this study submitted parallel tasks on local clusters with 2 eight-cores Intel E5-2680 2.7 GHz processors and introduced the reciprocity theorem into wavefield modeling in a mine, which can rapidly compute the Green’s functions using fewer global simulations. The current supercomputing resources can finish the construction of the 3D Green’s function database with fine grid for a whole mining region in two or three days. To further improve the computational efficiency, the workload of waveform inversion and storage requirements of the Green’s function database are reduced by the multi-scale mesh modeling strategy that is similar to Fichtner et al. [63]. The next steps can be realized through a laptop at a quick calculation. If there is more prior knowledge about the regions where MS events occur frequently, the combination of reciprocity theorem and multi-scale grid search could further reduce the wavefield modeling workload and storage requirement of the Green’s function database effectively, which can be generalized for simultaneous multi-source locating in the future.

In addition, under the current limitations of single-component sparse monitoring network and dataset with insufficient signal-to-noise ratio, we have not been able to obtain an appropriate elastic parameter model better than the currently used isotropic Vp tomographic model. Therefore, pure scalar acoustic wave modeling with a rough or simple velocity model is used (corresponding to acoustic pressure field records) to carry out studies in many current exploration-scale microseismic applications. However, the simplified velocity model is not accurate enough in the estimation of the strength of velocity anomalies and geometrical characteristic of the strong structure heterogeneities, such that the next stage task is to further improve the resolution of the structure parameter model. According to the complicated relationship between the spatial resolution of the target model and various control factors, the waveform inversion aiming at the broadband waveform dataset (such as joint tomography or full waveform inversion) can be performed based on accurate wavefield modeling considering that the practical waveform records presents a limited frequency band distribution. High-resolution structure parameter waveform inversion approach considers the wavefield finite frequency effect to improve the model resolution. In the future, when the mining region deploys multi-component monitoring network to obtain waveform data with higher signal-to-noise ratio, a more complete three-component Green’s function database based on the elastic wave equation modeling will be generated for further implementing full centroid inversion based on multi-component waveform data.

5.2. Estimation of the Source-Time Function

In many studies on a MS location and natural seismic source parameter investigation, the second-order Gaussian wavelet (Ricker wavelet) or simple Gaussian function is often used as the source-time function to generate synthetic waveforms. However, the time domain Gaussian or Ricker function is simply distributed in a symmetrical shape with respect to the main peak, and it has limited fitting performance to asymmetric or more complex waveforms. For the fractional-order Gaussian function proposed by Wang [48], better waveform fitting performance is obtained through searching for the best fractional-order derivate u and dominant frequency w₀ through normalized cross-correlation coefficients between the modeled waveforms and target data waveforms. However, its application is still focused on the monitoring operation with densely distributed stations and highly similar waveforms. The stations in a mine MS monitoring network are usually sparsely distributed and have large waveform variation with wave propagating distances, which result in more complex and dissimilar waveforms. In this paper, we introduced the fractional-order derivative and time scaling factor into a source-time function estimation procedure based on DD-SVD-FD wavelet. This technique fully considers variability and complexity of target-observed waveforms and can better fit the recording waveforms in a mine. The waveform-based grid search location method proposed here is a kind of source scanning algorithm by performing grid search strategy to evaluate the matching between the recording and synthetic waveforms, and the misfit minimization corresponding event location. Such that, an improved source-time function estimation and corresponding forward modeling is important for precisely computing misfit function value for each candidate source grid.

6. Conclusions

This study proposed a 3D multi-scale grid search and waveform inversion-based location method (GSWI location) into mine MS location, and the main conclusions are made as follows: (1) The combination of reciprocity theorem and multi-scale grid generation and search strategy effectively reduces modeling computational cost and storage requirement of the 3D Green’s function database. The establishment of the Green’s function database is favorable to the multi-scale grid search, which is helpful to find the globally optimal source location that fits the observed waveforms best. (2) The SEM method can effectively realize precise waveform modeling in complex structures, and the frequency-dependent sensitivity kernels of travel time anomaly have more accurate spatial distribution and intensity characteristics than the ray tracing methods. (3) Both the proposed DD-SVD-FD wavelet and fractional-order Gaussian function can well estimate the source-time function to generate synthetic waveforms fitting the recording waveforms. The synthetic waveforms obtained by the DD-SVD-FD source wavelet taking advantage of recording waveforms have better fitting performance, and have better data adaptability and versatility. (4) The average location error of the 3D GSWI for eight blasting events is only 15.0 m, which is smaller than that of the 3D ray tracing location method (26.2 m) with the same 3D velocity model. This proves that the waveform inversion-based location method can better constrain source location in complex regions than the simplified ray tracing method only based on arrival time pickings. In the future, taking advantage of resolution improved 3D velocity model through full waveform inversion approach relying on numerical wave-equation-solver can further raise location accuracy.