Abstract
Microseismic (MS) source location can help us obtain the fracture characteristics of a rock mass under a thermal-hydraulic-mechanical (THM) coupling environment. However, the commonly used ray-tracing-based location methods are easily affected by the large picking errors, multipath effects of travel time, and focusing and defocusing effects of rays in wavefield propagation, which are caused by the strong three-dimensional (3D) heterogeneity in mining areas. In this paper, we will introduce the rapidly developed waveform inversion-based location method into a mine MS field study. On this basis, the wavefields were modeled by utilizing a high-resolution 3D velocity model, a fractional-order Gaussian wavelet source-time function, and spectral element method (SEM) wavefield modeling. In order to reduce the computation cost of wavefield modeling, the 3D ray-shooting method based on a coarse grid was adopted to obtain an approximate MS event location. Based on this initial location, the multiscale waveform inversion strategy (from coarse to fine grids) and the L-BFGS iteration optimization algorithm were separately jointly selected to improve wavefield modeling speed efficiency and iterative convergence rate. Then, the IMS MS monitoring system set in the Yongshaba mine (China) and its tomographic 3D velocity model were used to conduct the synthetic test and application study. Results show that the source-time function based on the fractional-order Gaussian function wavelet can better fit complex recording waveforms compared with the conventional Ricker wavelet-based source-time function, and the waveform misfit during the L-BFGS iteration decreased rapidly. Furthermore, the multiscale waveform inversion method can obtain a similar location accuracy compared with the waveform inversion based on a single fine grid, and it can significantly decrease the iteration times and wavefield modeling computational cost. The average location error of the eight premeasured blasting events by the proposed approach is only 17.6 m, which can provide a good data research basis for analyzing MS event location and rock mass fracture characteristics in a mine.
1. Introduction
Microseismic (MS) source location, taking advantage of waveforms recorded by highly sensitive, broadband, and wide dynamic response sensors, plays an important role in the acquisition of fracturing location, displacement characteristics of a rock mass, inversion of a MS source focal mechanism, and magnitude along with velocity structures [1–6]. Furthermore, based on the information of MS locations, more accurate geological information can be extracted and further engineering plans can be guided. Motivated by this, researchers have conducted lots of studies on improving the precision and accuracy of MS source location methods, and the commonly used ones can be divided into ray-tracing-based location methods based on picking the travel time of a specified phase and wave equation-based location methods using waveform information directly. For the former method, the objective functions take advantage of residuals between theoretical and observed arrival time of a specified seismic phase. It is easy to realize and quick to compute. However, the location results are easily influenced by the arrival data quality and ray-tracing simplification in representing the finite-frequency nature of wavefield propagation. The ray-tracing-based location methods have been comprehensively reviewed in researches [7, 8]; thus, we will not describe it in more detail. Compared with the ray-tracing-based location methods, the wave equation-based location methods consider inherent finite-frequency effects of wavefield modeling in complex media. If the problem pertaining to computational efficiency is solved, the wave equation-based location method can have broader prospects for application.
1.1. Wave Equation-Based Location Methods
The wave equation-based location methods generally fall into two categories: wavefield migration- and waveform inversion-based location techniques.
1.1.1. Wavefield Migration-Based Location Methods
The wavefield migration-based location methods take the maximum spatial energy point of stacking back-propagation waveform from all sensors as the location result. In this way, the effects of large picking errors can be effectively reduced. Furthermore, it is not necessary to pick specific phase arrival times, which is conducive to automatic location of a huge number of MS events. Rentsch et al. [9] introduced the wavefield migration technique based on the Gaussian beam wavefield modeling into a hydraulic fracturing MS location. The Gaussian beam modeling technique can locally solve the wave equations in heterogeneous media, in which the Gaussian ray beams have a frequency-dependent width and can effectively model the wavefield focusing/defocusing and multipath effects. Li et al. [10] proposed a simplified fast migration-based location method with beam-forming-based Gaussian beams, which further increases the calculation computational efficiency. The method modeling complicated wavefields with the Gaussian beam approach presents a relatively high calculation computational efficiency and accuracy. Nevertheless, it can only solve wave equations around the central ray, and the modeling accuracy is still limited. Thus, the Gaussian beam migration-based location method is a fast reverse-time migration solution but without the highest location accuracy.
In contrast, McMechan [11] had utilized the high-accuracy finite difference method to model wave equations, and they solved the wave equations for reverse-time migration imaging to determine a MS event location. After that, some researchers modified the above wave equation-based location method by improving imaging methods [12, 13] and conditions [14, 15]. These wave equation-based reverse-time migration location methods have been widely used in various fields, such as natural earthquakes [16, 17], oil- and gas-exploitation-induced earthquakes [18], and volcanic earthquakes [19]. In the above applications, a regularly distributed and dense sensor array is often used to obtain high-resolution images. However, monitoring sensors are distributed rather irregularly and are quite limited in number (generally less than 50) for mining engineering. In addition, as a key factor for reverse-time migration location results, the selection of imaging conditions usually depends on experience and field condition, which limits the generality of the migration-based location methods.
1.1.2. Waveform Inversion-Based Location Methods
The second family of location techniques based on wave equations is the waveform-based iterative inversion of source parameters, which requires a numerical wave equation solver to model broadband wavefields in 3D media. The reverse-time migration location technique only needs a small number of full-waveform simulations, while the waveform-based iterative inversion, also called full-waveform inversion (FWI), requires more wavefield modeling for iterative inversion and a sophisticated data preprocessing. It considers the time-dependent frequency changes of different seismic phases as well as the complexity of wavefield propagation in 3D heterogeneous media completely; thus, it constrains source parameters by extracting as much waveform information as possible. As the ultimate technique in seismic inversion, it should be mentioned that FWI needs to deal with the problem of low signal-to-noise ratio (SNR) waveform recordings, which is more challenging than the wavefield migration-based location methods. Therefore, some advanced data signal processing techniques are introduced into the waveform inversion-based location methods, such as signal denoising and adaptive waveform windowing.
For an 3D velocity model, the highly accurate wavefield modeling required by FWI is closely related to source location, origin time of events, source moment tensor, and source-time function. Before the 21st century, due to the limited computing resources, FWI was rarely applied in a location study. Since then, with the development of numerical modeling and computer technology, a series of breakthroughs has been made in source location based on FWI. Kim et al. [16] and Rodriguez et al. [20] fixed the horizontal location of a source, and then they jointly inverted the source moment tensor and source depth in a high-resolution 3D velocity model, through iterative fitting of theoretical numerically modeling waveforms and observed waveforms. However, the source moment tensor and source depth are strongly coupled in the modeled waveforms, which increases the nonlinearity of the inversion problems and requires more iteration times to ensure the inversion convergence. Tong et al. [21] further simplified the number of source parameters in the inversion problem. They treated the moment tensor as a known variable and introduced it into the source-time function modeling. On this basis, considering that the cross-correlation travel time of a specific phase is not very sensitive to source-time function, they derived a waveform inversion-based location method using 2D acoustic wave equation solving and cross-correlation travel time measurement. Tong et al. [21] obtained good location results in regional seismic location by applying a simplified 2D sliced velocity model. In addition, Kim et al. [16] and Rodriguez et al. [20] adopted three-component waveform recordings in FWI, while Tong et al. [21] used single-component waveform recordings for inversion. The latter is more commonly used in the fields of oil and gas exploration and mining engineering.
In the above waveform inversion-based location methods, the elegant adjoint method is used to calculate the source-parameter-gradient vectors of an objective function. To further reduce the computational cost, Huang et al. [22, 23] utilized the scattering-integral method [24, 25] to halve the computation of wavefield modeling. They used the truncated Newton iterative algorithm for inversion, which has a faster convergence rate than the conjugate gradient method used by Kim et al. [16] and Tong et al. [21]. Huang et al. [22] effectively separated the unknown origin time of MS events from source-time function, and their study mainly focuses on the strategy of further improving waveform inversion efficiency and reasonably separating source parameters. In their separation strategy of source-time function, Huang et al. [22, 23] used a reference channel for cross-correlation operation and defined the corresponding objective function, thus decreasing the nonlinearity of inversion problems. Finally, they discussed influences of the velocity model accuracy and signal-to-noise ratio (SNR) on the FWI location results through 2D and 3D synthetic tests. They successfully performed the FWI-based location for an actual MS event recorded by a very dense single-component sensor array and used a relatively simple 1D velocity model.
1.2. Wavefield Modeling Based on Numerical Wave Equation Solver
The resolution and accuracy of a waveform inversion result heavily depends on the accuracy of wavefield modeling. For a waveform inversion-based location application in complicated heterogeneous media, the wave equations will be numerically solved in a 3D velocity model. The generally used numerical methods for wave equation solution include the finite element method (FEM) [26, 27], the finite difference method (FDM) [28], the boundary element method (BEM) [29], the pseudospectral method (PSM) [30], and the spectral element method (SEM) [31–37]. The finite element method is flexible in grid generation and can accurately simulate wavefields in various complex structure models. However, the conventional finite element method is usually used for modeling waves through a low-order polynomial interpolation, which results in poor performance in modeling high-frequency signals [38]. The finite difference method is simple and has high calculation efficiency, while the commonly used FDM grid subdivision is difficult to use when handling complex geological structures and interfaces for preset numerical accuracy. The boundary element method can effectively deal with complex boundary problems, but it is often limited to model linear wavefield response in very simple media. The pseudospectral method shows high calculation accuracy, but it is difficult to use when processing complex regions and boundaries.
For the spectral element method, the field parameter and variable in a wave equation are interpolated in subdivided grids with a high-order polynomial. It combines the high accuracy of the pseudospectral method and the flexible grid subdivision property of the finite element method, thus having unique advantages in modeling broadband wavefields in a 3D complex structure. Komatitsch and Vilotte [31] and Komatitsch et al. [32] programmed SEM parallel computing and open programs for modeling acoustic and elastic waves. Fichtner et al. [33] and Tape et al. [34] modeled wide-frequency-band wavefields using the 3D SEM, and they developed the most popular adjoint tomography and FWI approach for imaging a large-scale 3D velocity structure. Furthermore, Kim et al. [16] and Liu et al. [34] carried out the SEM-based FWI for natural earthquake focal mechanisms. The waveform misfit and focal mechanisms are obviously improved in comparison with the inversion results obtained through simplified 1D velocity models. In recent years, SEM has been developed to successfully model and image 3D subsurface voids, which further improves its scope of application in more complex shallow areas [39].
Utilizing SEM, theoretical seismic broadband waveforms in 3D media were modeled. SEM has been widely used in the inversion of large-scale velocity structures and focal mechanisms of natural earthquakes. However, this study applied 3D SEM into a complicated waveform modeling and source location in a small-scale complex mine region for the first time.
1.3. What Will Be Done in This Paper?
Inspired by ideas of Tong et al. [21] and Huang et al. [22, 23], this study will conduct waveform inversion-based source location in a mining region. Compared with natural earthquake and oil exploration signals, mine MS signals attenuate more seriously with propagation distance. In addition, the complex velocity structures in a mining area can lead to very complicated wavefield interference and strong scattering. Furthermore, due to the limited number of available sensors in a monitoring network, the source-time function cannot be simply estimated by the traditional average strategy of stacking waveforms. Moreover, the nonuniform distributed monitoring sensors and large propagation distance differences from the source to different sensors will increase the nonlinearity property of the FWI problem searching for the best source location parameters, which minimize the fitting difference between theoretical and observed waveforms. To handle this, an objective function based on cross-correlation travel time for broadband target phase waveforms was established to enhance the convergence of waveform inversion location. Furthermore, this study introduced a high-resolution tomographic 3D velocity model and the high-accuracy SEM method for modeling broadband wavefields in a complicated mine region. When the highest cut-off frequency for the SEM wavefield modeling is doubled, the required grid size in three directions and the time step should be reduced by half, in which the computational cost is 16 times that of the original lower frequency band-based modeling. The peak frequency of the actual waveforms is quite different due to influences of propagation distance difference and strong attenuation. In order to reduce the computation of wavefield modeling and storage requirement for scattering-integral calculation, a set of fractional-order Gaussian function- [40, 41] based multiband source-time function estimate strategy and multiscale waveform inversion strategy [42, 43] are developed.
The contributions and innovations of this study are shown as follows: (1) The waveform inversion-based location method was introduced to a complex mine MS location problem, and a tomographic 3D velocity model and a high-accuracy SEM method were utilized for modeling broadband wavefields. (2) A coarse-grid 3D ray-shooting method and multiscale waveform inversion strategy were combined to reduce the computation and ensure location accuracy of waveform inversion-based locations. In addition, the L-BFGS algorithm was used to reduce storage of the Hessian matrix required by the second-order Newton method. (3) The fractional-order Gaussian wavelet can obtain a better source-time function compared with the conventional and simple Ricker wavelet, thus improving the fitting performance between theoretical waveforms and complicated recorded waveforms. (4) With the above integrated technology, the average location error of eight experimental blasting events was only 17.6 m, which illustrates that the proposed method can provide a good basis for analyzing fracturing zones of a rock mass and further inverting the source mechanism. The rest of the sections are arranged as follows: Section 2 introduces the waveform inversion-based location method, the L-BFGS algorithm, the fractional-order Gaussian wavelet for the source-time function estimate, and the multiscale waveform inversion strategy. In Section 3, the wavefield modeling performance of the SEM method and the location accuracy and validity of the multiscale waveform inversion are tested using a 3D velocity model. In Section 4, the waveform inversion-based location method is applied to locate eight blasting events with known locations, and the location results are compared with previous studies. Section 5 discusses the significance of the waveform inversion-based location method, the multiscale waveform inversion strategy, and the fractional-order Gaussian wavelet source-time function estimate in a MS location. Finally, Section 6 summarizes the research and presents future prospects.
2. Theory
2.1. Waveform Inversion-Based Location Method
2.1.1. Source Model Parameterization and Waveform Inversion Misfit
Waveform inversion for source location is solved by iteratively fitting observed waveform data by modeled theoretical waveform data . If a point-source approximation is adopted, the modeling of theoretical waveform data at the sensor location is closely related to the source location , origin time of source events, point-source moment tensor , and source-time function . However, it is very difficult to simultaneously invert so many coupled source parameters; thus, proper simplification is needed. For a simple blasting or MS event, the moment tensor can be approximated by some data preprocessing procedure or incorporated into the source-time function . The origin time of a MS event can be considered as an overall translation shift parameter included in the source-time function , which is jointly determined by source magnitude, geometric configuration between the source and receiver locations, and the neighboring medium information. In conclusion, the point-source model parameters for simulating seismic waveforms can be simplified, and this mainly depends on the estimation accuracy of the source-time function and event spatial location. When the source-time function is accurately and correctly simulated, only the inversion of source location is needed [22], that is, the source model parameters are finally simplified as .
The waveform inversion technique needs to define a stable and reasonable criterion for measuring waveform misfit, namely, an objective misfit function needs to be established at first. In general, the model parameters in waveform inversion are directly evaluated and inverted by the waveform difference between observed and synthetic seismograms by using the least square norm, where this kind of waveform misfit requires relatively high-quality waveform data. In view of the inversion theory, there is a strong nonlinearity in the least square waveform inversion. Compared with the direct waveform difference measurement for inversion, there will be a higher robustness and quasilinearity convergence when using the cross-correlation between multiband data and synthetic waveforms to measure the travel time differences of a specified phase for the corresponding waveform inversion [36]. Therefore, the L2 norm of cross-correlation travel time difference was adopted as a band-pass-filtered waveform misfit function to invert the source location in this research. The detailed derivation for this kind of waveform-based cross-correlation travel time measurement and relevant waveform inversion procedure will be introduced below.
The cross-correlation travel time difference of the theoretical and observed waveforms for the target seismic phase is denoted as the objective function , and waveform inversion needs to iteratively solve model parameters that minimizes the . In general, the following equation is used to update the model parameters based on the gradient of an objective misfit function (i.e., , ), where indicates the total iteration number and represents the updating step length of the th iteration, which can be determined by the classical Wolfe conditions [44]. Moreover, denotes the direction vector, which is correlated with the gradient of misfit function with respect to source location parameter . It is clear to see that accurately solving is the key to a successful inversion of source location. Whether the first-order gradient-based descent algorithm or the second-order Newton iterative algorithm is used, the direction vector is the result of certain mathematical operations on the gradient for the objective function. Therefore, the critical step is to efficiently and accurately calculate the gradient vector based on the chain rule: where represents the Frechet derivative (gradient) of wavefields with respect to source location , and the derivation process will be shown in Section 2.1.2.
2.1.2. Frechet Derivatives Based on Acoustic Wave Equation
To guarantee accurate location result and avoid the high computational cost of solving the three-component elastic wave equation but still taking finite-frequency effects into consideration, we assume that the propagation of mine microseismic waves satisfy the acoustic wave equation, and its frequency domain can be expressed as follows: where denotes the wave velocity value at an arbitrary point in the medium model; indicates the frequency spectrum of the source-time function , which usually shows a complex spectral characteristic.
Perturbation at the source location can lead to a wavefield disturbance . Supposing that wavefields excited by the source at and are and , respectively, the following equation is obtained by defining a perturbation wavefield :
If the source-time function remains unchanged when the source location changes, the acoustic wave equation that satisfies can be obtained through equation (2) and expressed as follows:
Furthermore, due to the spatial translation invariance of the differential operator in a wave equation, the source wavefield updated from to is precisely equivalent to the wavefield received at in the medium model but excited by as the source location and with a shifted velocity model of [44]. By referring to equation (2), it is assumed that wavefield of source location and velocity model satisfies the following:
According to the above analysis, is obtained, thus further reaching the following equation through the first-order Taylor series expansion:
By combining equations (6) and (3), we obtain
On the other hand, the velocity model in equation (5) can be regarded as a result after the velocity model in equation (2) experiences a disturbance , namely, . Substitute this first-order approximation into equation (5), we get
By performing the first-order Taylor series expansion again on equation (8), we can obtain
By substituting equation (7) into equation (9) and subtracting equation (2) from equation (9), we can obtain
Consider that in the first-order Born approximation, the wavefield excited after disturbance of source location in equation (10) can be approximated with the background wavefield solved in the original velocity model before the disturbance, which gives us
By comparing equation (11) with equation (2), it can be inferred that the perturbation wavefield is excited by the following spatial distributed source:
Taking advantage of the corresponding wave equation and representation and reciprocity theorems of elastic mechanics, through convolution between a spatial distributed source and the point-source Green’s function, the Frechet derivative (also called sensitivity kernel) of the wavefield at the receiver with respect to source location disturbance is finally obtained, which is shown as follows: where is the point-source Green’s function at the receiver corresponding to the frequency domain in the background velocity model .
Furthermore, equation (13) can be divided into two terms: where and are respectively given as follows:
According to equations (1), (13), and (14), the gradient vector of the objective function required in the iterative waveform inversion can be calculated as follows:
Then, the updating direction can be obtained through an appropriate iterative optimization algorithm, and the source model will be successively updated until the waveform misfit satisfies the preset threshold value. can be considered as a wavefield perturbation caused by geometrical disturbance of the source spatial location, which mainly reflects the effects of the spatial difference of the source location and the first-order effect of location error. However, is more complex, and it is not only the wavefield perturbation induced by the corresponding source region velocity change resulting from the inaccurate source location but is also considered the overall travel time and waveform perturbation of the global velocity model variation on the main energy propagation path caused by the source location change. If the velocity model is very heterogeneous, the influences of on waveforms can be further enhanced [22]. For example, when there is a spatial error in source location, complicated focusing and defocusing effects on the corresponding wavefield propagation path can be “recorded and carried,” and finally reflected by records received by relevant sensors. Hence, reflects the finite-frequency nature of wavefield propagation and complex interference of different seismic phases, thus affecting the cross-correlation travel time location inversion and tomography using broadband waveforms of specified seismic phases.
In conclusion, the FWI and cross-correlation travel time-based waveform inversion are both based on computation of the finite-frequency sensitivity kernel (equation (13)). Back to equation (1), the objective function defined in this study is the cross-correlation travel time difference of windowed waveforms of selected seismic phases. The cross-correlation function of discrete waveform time series is defined as follows:
The cross-correlation travel time difference is determined by the maximum value of equation (18), which is shown as where and is the time window function of specified seismic phases in the interval [21]. Through the summation of cross-correlation travel time difference of all windowed waveforms, the objective misfit function for location waveform inversion is shown as follows: where represents the number of windowed waveforms for location inversion. By utilizing equations (1), (13), and (21), as well as the iteration with the appropriate optimization algorithm, source location based on waveform inversion is performed.
2.1.3. FWI Using the L-BFGS Algorithm
The second-order quasi-Newton algorithm-L-BFGS was selected for our iterative location waveform inversion, which involves the calculation of the approximate Hessian matrix. The approximate Hessian inverse (equation (22)) of the BFGS algorithm is used for computing the model updating direction: where indicates the product of vectors; represents the difference between the current model and the previous model in the iterative process; denotes the adjacent gradient change of the objective function.
Based on the L-BFGS algorithm, the updating direction of the th iteration is shown as follows:
Each Hessian approximation of the commonly BFGS algorithm needs to use all the historical vectors and , which will take up a lot of memory. Therefore, the more efficient and accurate L-BFGS [44] was utilized here to approximate the Hessian matrix. The L-BFGS only saves the latest iterations to calculate the . By this way, the iteration computation and required memory will be obviously reduced, which brings a unique advantage to the calculation efficiency.
2.2. Fractional-Order Gaussian Function for Source Time Function Estimation
Before calculating gradient of the objective misfit function (equation (1)) and modeling the corresponding waveforms, it is necessary to build a source-time function as accurately as possible. The accuracy of source-time function modeling will greatly affect the convergence and results of nonlinear waveform inversion. In most natural earthquake and oil exploration MS event locations, the popular Ricker wavelet (the second-order derivative of the Gaussian function) is usually selected as the source-time function for wavefield modeling only by tuning its dominant frequency. The definition of the Gaussian function is shown as follows: where and indicate the time series and moment corresponding to the time of wavelet center; represents the control dominant frequency, which is correlated with sensor frequency response.
By transforming equation (24) from time domain to frequency domain, we obtain:
The symmetric Ricker wavelet has a good analytical characteristic. The precise source-time function modeling by the Ricker wavelet is based on high-similarity waveform data and the relatively consistent frequency response of recorded waveforms. However, the more general asymmetric distribution of the time domain waveform can reflect the practically frequency-dependent scattering and attenuation behavior of seismic wave propagation in realistic media. Moreover, if the recorded waveforms for the wavefields excited by the same source are very different and the corresponding dominant frequencies are rather inconsistent (which is quite the usual observation in a mine field), it is not appropriate to use such a uniform and simple Ricker wavelet to model the source-time function. In view of this problem, Wang [40] proposed a fractional-order Gaussian function to model the source-time function for observed data waveform modeling and spectral analysis. The spectra of the fractional-order Gaussian function relating to time can be obtained through multiplying by frequency factor [40], shown as follows: where indicates the fractional order of a time domain derivative. Changing and at the same time can increase the freedom degrees for estimating the source-time function to better model the real complex signal waveform through wave equation solving. Then, equation (26) is normalized by multiplying the coefficient , and the normalized fractional-order Gaussian function in the frequency domain is expressed as follows:
Wang [40] found that the time domain fractional-order Gaussian functions with different orders have very different shapes: , the corresponding standard Ricker wavelet is symmetrically distributed on both sides of the main peak; as gradually decreases from 2, the first time domain main peak gradually rises, while the last main peak gradually reduces; an opposite law is obtained for and . The introduced fractional-order Gaussian function can better fit the observed broadband waveforms, which are controlled by strong attenuation and scattering in the complex media. Furthermore, he proposed a set of rapid analytical algorithms to determine the fractional order and optimum dominant frequency of the corresponding source-time function; then, the rapid spectrum of the data waveform was analyzed. Owing to the simple form of fractional-order Gaussian function, this study will take advantage of the scaling transform properties and cross-correlation coefficient between a target and the synthesized waveform to determine the fractional order and optimum domain frequency of each waveform, by using a two-parameter grid search algorithm. The application examples show that the computational cost of the two-parameter grid search can be completed instantaneously.
2.3. Multiscale FWI Strategy
The FWI for MS event location requires modeling a broadband wavefield in 3D complex media. Generally, the number of grids in the shortest wave length of a modeled wavefield and the Courant-Friedrichs-Lewy (CFL) condition are used to evaluate and control the accuracy and stability of the numerical modeling method. The number of grids per shortest wave length in numerical wavefield modeling is usually governed by the following empirical equation [45]: where is the number of grids within a single wave length in the grid subdivision space , is the shortest period of the frequency band of the modeling wavefield, denotes the average grid size, represents the empirical number of grids per wave length. This value is usually set to about 4.5 for the elastic-dynamic equation modeling. Thus, equation (28) can be rewritten as follows:
When the shortest period of modeling frequency bands is smaller than the , the modeled wavefield using the average grid size becomes increasingly inaccurate, especially for the high-frequency content. From equation (28), we know that in order to accurately calculate the theoretical high-frequency (short period) waveforms, the grid size should be reduced proportionally when is fixed.
The CFL condition is the critical constraint to control the numerical modeling stability of wave equations, which is given as follows: where indicates the spacing between the adjacent grids, which will be changed in space for the adaptive variation grid technique. and represent the wave velocity of the model and the time step of finite difference for the time scheme, respectively. stands for the Courant number and is usually set to .
In accordance with equation (30), we can know that as the grid size decreases, the time step will reduce correspondingly. Thus, equations (28) and (30) show that for a 3D grid subdivision, when the frequency band of wavefield modeling is doubled, the grid sizes in the , , and directions and time steps will be reduced by half, which results in the computation workload increased 16 times that of the original coarse grid.
Waveform inversion is an optimization problem requiring nonlinear iteration. It is difficult to avoid a large amount of wavefield forward modeling in 3D complex media. This is also the core limiting factor when using waveform inversion as the ultimate technology for seismic inversion. Therefore, this study introduced the multiscale grid FWI strategy, which has been used in resource exploration inversion [43] and natural earthquake adjoin tomography [42]. The idea of the multiscale grid FWI strategy is as follows: firstly, we perform waveform inversion on the low-pass filtered (corresponding to a coarse grid subdivision) waveform record to obtain a good initial location after inversion convergence. According to the above discussion, waveform inversion for low-frequency content only requires a coarse grid and a longer time step, and the computational cost is much lower than that of waveform inversion for high-frequency content. After the low-frequency waveform inversion has converged, the spatial grid is finely divided and the velocity model is interpolated. Then, the next stage waveform inversion moves to the high-frequency content of waveform recording, which obtains a higher precision source location. The multiscale waveform inversion strategy can effectively reduce the nonlinear property of the inversion problem and save a lot of computation resources. For a mine MS inversion in this study, the waveform received by the long-distance sensor is affected by absorption attenuation and strong scattering on the wavefield propagation path, resulting in the domination of the corresponding recording waveform by low-frequency information. Therefore, we perform this multiscale inversion strategy from low-frequency to high-frequency content of the waveform step by step.
We applied a low-pass filtering- (LPF-) coarse grid to the long travel distance waveform with a low SNR, and which can achieve an overall low-frequency waveform inversion of the source location. Then, the waveform inversion was performed on high-frequency data (fine grid) to further improve the accuracy of a source location.
3. Synthetic Test of Multiscale Waveform Inversion for Source Location
Many synthetic tests on waveform inversion for source location have been done to verify the effectiveness and robustness of location results [21, 22]. The complex 3D tomographic velocity model of the Yongshaba mine [46] shown in Section 4 was used for this synthetic test. The synthetic test was conducted on the multiscale waveform inversion algorithm. The 3D spatial distribution of MS monitoring sensors for the target mining area is shown in Figure 1, and the cyan triangles indicate the vertical projection locations of the sensors on an inclined 2D plane, which is determined by the least square value of summation of the vertical distances from the 28 sensors to certain spatial 2D plane. All the sensors have a sampling frequency of 6,000 Hz, and they are located in three mining sublevels ( side view on the upper left corner of Figure 1 and coordinate represents the altitude). The sensors are distributed unevenly, and there are many complex structures (e.g., cavities and tunnels) in the mine region, which make it very difficult to locate earthquakes with high accuracy.
Eight test events (black stars in Figure 2) were selected from strong high- and low-velocity anomaly regions to verify the effectiveness of multiscale waveform inversion-based location strategy; the triggered sensors were, respectively, the same as the sensor distribution of the eight blasting events described in Section 4. The Ricker wavelet with a dominant frequency of 120 Hz was selected for the SEM wavefield modeling, and the schematic diagram of the wavefront evolution is shown in Figure 3 with a time interval of 0.04 s. The wavefield amplitude is normalized by the maximum amplitude of full space wavefields at the entire time series. It can be seen that the wavefield energy attenuates rapidly with the increase of propagation distance. The background of Figure 3 is the 2D slice heterogeneous velocity model shown in Figure 1, and it can be found that the low-velocity anomaly has remarkable focusing effects on transmitted waves, while the high-velocity anomaly exerts defocusing effects. The wavefield focusing results in multipath effects for the simplified ray-tracing modeling, and it is difficult to use the ray-tracing method for modeling the defocusing region effectively [47]. Compared with the ray-tracing method, the SEM numerical method can effectively simulate wavefields in focusing and defocusing regions, which lays a good foundation for waveform inversion-based location in complex structures by completely considering the complicated practical wave propagation effects.
Test event 7 was selected to illustrate the multiscale waveform inversion-based location strategy. The focal mechanism was set as an explosion source, and the Ricker wavelet with a dominant frequency of 120 Hz was selected to generate the theoretical waveforms. Firstly, the initial location was obtained by the 3D ray-tracing method in our recently published paper [47] (yellow star in Figure 4). In this study, a two-stage waveform inversion was used for the multiscale iterative inversion because the complex influences of 3D viscoelastic absorption in the media were not considered. A more complex, multistage and multiscale wavefield modeling and inversion strategy will be developed and carried out in the future. The coarse grid size in the first-stage FWI for the waveform with low-frequency content was twice that of the fine grid in the high-frequency waveform modeling, and the computational cost of the corresponding waveform modeling was only one-sixteenth that of modeling in the fine grid. Figure 4 shows that the first-stage waveform inversion of this test event based on the coarse grid converges after six iterations, and the waveform misfit decreases rapidly. Then, the location results based on the coarse grid inversion were set as the initial location in the subsequent fine-grid FWI iteration, with the location results stabilizing after just three more iterations. The L2 norm waveform misfit between data and synthetic waveforms is defined in equation (31), and the misfits are normalized through division by those from initial locations at the beginning of each stage of FWI. In order to illustrate the efficiency of multiscale waveform inversion, the broadband waveform inversion using a uniform fine grid was selected as a comparison. Waveform inversion based on a single uniform fine grid stabilizes after 11 iterations, and the wavefield modeling computational cost and iteration times increase largely for this uniform fine grid inversion. However, the location errors of multiscale waveform inversion are nearly equal to those of waveform inversion based on a uniform fine grid, and the location errors are 1.5 m and 1.1 m, respectively. Location results of the eight test events based on different waveform inversion strategies are listed in Table 1. It can be seen that location errors of all the test events based on multiscale FWI strategy are smaller than 2 m. However, the computational cost of multiscale waveform inversion is only one-quarter to one-third of the waveform inversion based on a single fine grid, that is, the calculation efficiency rises obviously. where indicates the number of windowed waveforms involved in the MS location inversion.
(a)
(b)
(c)
4. Application
In this section, eight blasting events with known locations recorded by the IMS monitoring system of the Yongshaba mine, Guizhou province (China) were taken as application test data, and the blasting event locations are shown in the paper [47]. The hanging wall of ore body in the Yongshaba mine is dolomite with a good stability [48], while the crushed shale and sandstone were presented in the footwall, which showed an obvious difference in the -wave velocity model. Furthermore, a large number of cavities were left in the early stage of mining by using the open-stope mining method, and the stress was adjusted accordingly. For these reasons, obvious high-velocity and low-velocity regions existed [48–52] in the mining area (Figure 3), which increased the difficulty of an accurate MS event location. Therefore, a proper 3D velocity model should be considered in the location investigation.
4.1. Source Time Function Estimate by Fractional-Order Gaussian Function Fitting
Blasting event 2 was taken as the illustration example in this section and in Subsection 4.2. To reduce nonlinearity of waveform inversion, the waveforms of specified seismic phases were extracted using the fast windowing method proposed by Wang et al. [47] (black curve in Figure 5(a)). The windowed waveforms can effectively truncate the subsequent complex multiple scattered wave or interference signal. It is known that with the increase of propagation distance, the high-frequency band attenuation is quicker than that of the low-frequency band in the waveform, thus reducing the overall effective frequency bands. Therefore, the main frequency and spectral distribution characteristics of the windowed waveforms of the initial phase vary greatly for different sensor positions (Figure 5(a)). The maximum cut-off frequency of wavefield modeling is defined as the highest frequency value corresponding to half of the peak value of the waveform amplitude spectrum (Figure 5(b)). The largest of the recorded waveform determines the grid size and time step of numerical modeling. We generally need to find a proper waveform fitting evaluation criterion when using the fractional-order Gaussian function to estimate a global optimal source-time function. In this study, the fitting performance is evaluated based on the normalized cross-correlation coefficient between modeled and observed waveforms, which is also frequency dependent.
(a)
(b)
For a windowed waveform, the observed waveforms can be fitted by adjusting the fractional order and the dominant frequency of a fractional-order Gaussian function. In a two-parameter grid search, the waveforms recorded by all receivers were simultaneously fitted with the fractional-order Gaussian function to obtain an overall fixed optimal fractional order . All observed waveforms are fitted with this same value, while the dominant frequency changes for each waveform fitting, thus obtaining the red curves in Figure 5. Meanwhile, we also search for the optimal fractional order for fitting each individual waveform by using an independent single fractional order for comparison (blue curve in Figure 5). It should be noted that each waveform and its amplitude spectrum are normalized, and the -wave polarities are corrected. Figure 5 shows that for waveform records with complex spectrum characteristics, the fitting performance of waveforms using the optimal fractional-order Gaussian function by individually searching for the value for each sensor is better than that of using a unified fractional-order value. However, we selected the unified fractional-order Gaussian function to fit the recorded waveforms in this study, and the reason is shown in Discussions. In general, the red curve in Figure 5 shows that a grid search with a unified value can still yield good waveform fitting performance and hence can effectively estimate the source-time function. However, it shows moderate fitting performance for a few complex windowed waveforms, probably because the complex wavefield propagation effects and low SNR caused by background noise increase the recorded waveform complexity (e.g., sensor 3 in Figure 5).
To analyze the fitting performance of the fractional-order Gaussian source-time function under a uniform before and after filtering, nine windowed MS waveforms were randomly selected from MS event 7. Figure 6(a) shows that the minimum correlation coefficient between the synthetic waveform and the recorded waveform is 0.8296, where the synthetic waveform is generated by the fractional-order Gaussian function with a unified value of 0.35 and a dominant frequency of about 100 Hz as the source-time function. This waveform comparison indicates that a good fitting performance is obtained on the whole. For each recorded waveform, the corresponding LPF could be conducted on the synthetic waveforms with the highest cut-off frequency depending on its own optimal dominant frequency. Figure 6(b) demonstrates that the correlation coefficient of a waveform between a synthetic waveform and data under a roughly unified LPF (cut-off frequency of 120 Hz) is obviously improved and filtered waveform fitting performance becomes better. Therefore, these comparisons prove the robustness of the fractional-order Gaussian function for modeling source-time function under various types of frequency band filtering. After confirming the highest cut-off frequency of each record using the fitted fractional-order Gaussian source-time function and the relevant optimal dominant frequency, the wavefield modeling frequency bands for waveform inversion can be determined. Furthermore, the computational cost can be largely reduced by using the multiscale waveform inversion strategy, and the misfit of multiscale waveform inversion converges faster with a reduction in the inversion nonlinearity.
(a)
(b)
4.2. Location Example
The blasting event 2 was taken as the illustration example again; the initial location result obtained by the 3D ray-tracing method based on a 30 m grid spacing is shown by the yellow star in Figure 7(b). The SEM-modeled waveforms based on this initial source location are indicated by the blue dotted lines in Figure 8. The first-stage iterative LPF waveform inversion process based on the coarse grid is demonstrated by the red star in Figure 7. The convergence model of waveform inversion based on the coarse grid at a low frequency (LPF range: ) is utilized as the initial model of the second-stage broadband waveform inversion based on the fine grid. From this, iterative inversion is further performed by using broadband waveforms based on the fine grid (the blue star in Figure 7). The final convergence source location and corresponding modeling waveforms are shown by the green star in Figure 7 and the red waveform curve in Figure 8, respectively. It can be seen that the red curves fit the recorded waveforms better than the blue curves on the whole, with a few waveforms having equivalent fitting performance (e.g., traces 2, 11, and 13). Moreover, the following phenomenon may appear: the waveform fitting performance generated from the initial source location is even superior to that generated from the final FWI convergence source location. The reason may be that the tomographic velocity model obtained by the 3D ray-tracing method still has insufficient accuracy for numerical broadband wavefield modeling that requires a higher-resolution velocity model considering finite-frequency effects.
(a)
(b)
(c)
The source locations and waveform misfits during the two-stage multiscale waveform inversion for this event are shown in Figure 7; it shows that the waveform misfit decreases on the whole. Particularly, the waveform misfit reduces by about 50% during the low-frequency coarse-grid waveform inversion. However, the waveform misfit decreases slightly (only about 7%) based on the fine grid in the second-stage FWI. This is because the waveform data with a dominantly low-frequency band spectrum has a better azimuth coverage, and there are more sensors far away from the source already providing a sufficient and efficient constrain on source location to ensure a better location convergence for coarse-grid FWI. However, waveform inversion based on the fine grid is mainly a refinement process involving short-distance broadband waveform information explorations in the inversion, and the waveform located far away lacks information at high frequencies due to strong scatter and attenuation. The FWI converging locations approach the premeasured blasting location as a whole, and the final location error is only 15.4 m.
4.3. Location Results
Location errors of the eight blasting events using the multiscale waveform inversion are shown in Table 2. It can be seen that blasting event 6 has the largest location error of 31.6 m, while the average location error is only 17.6 m, which well satisfies the requirements for location accuracy of an engineering MS monitoring. These eight blasting events were also located in some other studies as a methodology validity test. The average location error using ray-tracing methods in a uniform velocity model is larger than 39.5 m [53–56], while that in the same 3D velocity model is 26.2 m [47]. It is obvious that location results of the waveform inversion-based location method proposed in this study are superior to those of the traditional ray-tracing-based location methods even using the same 3D tomographic velocity model.
5. Discussions
5.1. Waveform Inversion-Based Location
Currently, the travel time-based ray-tracing method is the most commonly used method in a mine MS event location, where the wavefield propagation is infinite high-frequency approximated and treated as a simple geometric ray of a specified seismic phase. However, the frequency bands of the practical spectrum of recorded seismic wavefields are finitely distributed, thus reflecting significant finite-frequency effects in wavefield propagation. The underground structure through which the seismic waves propagate is not a perfect elastic body, and the high-frequency energy of wavefields is severely absorbed and attenuated in the propagation process. The seismic signal spectra received by the sensors are actually affected and modulated by various complex factors, usually resulting in a complex band-limited distribution. Therefore, the simplified infinite high-frequency approximation ray theory fails to take into account travel time delay or advance caused by medium disturbance with the scale corresponding to the wave length of the dominant frequency content of the wavefield on the ray path. This is the famous wavefront healing effects resulting in the amplitude of first-arrival waves being compensated and interfered by the surrounding wavefield wavefronts, which reduce the accuracy of location results based on arrival time picking. Compared with ray-tracing-based location methods, the waveform inversion method precisely considering finite-frequency fluctuation effects can express the effects of various scale parameter disturbances on multifrequency-band travel time. These effects are quantified based on the distribution characteristics of sensitivity kernels of the cross-correlation travel time difference of frequency-dependent waveforms with respect to structure parameters or source locations [57]. In other words, the waveform inversion method does not have an inversion error caused by wavefront healing effects. It can more completely explore waveform information compared with the single travel time measurement and corresponding ray-tracing-based location methods, thus better inverting location results. The application examples in this study show that by combining with the high-accuracy 3D tomographic velocity model, waveform inversion can accurately locate MS events in a very complex mine. The wavefield modeling with the SEM method is not only a highly accurate and parallel efficient modeling method, but it is also a method that has unique advantages in future modeling of tunnels and cavities in a mining area. The location error in this application results from insufficient estimation of fine structures and anomaly amplitude of the 3D velocity model, so the FWI can be carried out for further structural refinement in the future.
5.2. MultiScale Waveform Inversion Strategy
Precise wavefield forward modeling is the most important aspect in iterative waveform inversion-based location. In the mine field, it needs to model broadband wavefields, that is, it should be capable of performing numerical modeling with a fine grid. According to the sampling condition of the accurate simulation of wavefields and CFL condition, with the increase of wavefield modeling frequency band, the grid size and time step reduce correspondingly and the computation of forward modeling of waveforms dramatically rises in a quartic manner. Therefore, an appropriate technique should be adopted to reduce the computational cost of waveform inversion, which should balance the high modeling precision and high calculation cost. Considering the strong dependence of iterative waveform inversion on the initial model, a fast 3D ray-tracing method based on a coarse grid was firstly implemented to initially locate the source. This could decrease the locating object range in comparison with a conventional global iteration algorithm by FWI, thus effectively reducing the times of wavefield modeling required by waveform inversion. On this basis, the strategy of multiscale waveform inversion based on variable grids was utilized. Firstly, the wavefields were modeled with a coarse grid corresponding to low-frequency filtered waveforms, and a favorable converged initial location was obtained through waveform inversion. Then, wavefield modeling based on a fine grid and the second-stage waveform inversion were performed to realize a refined location of the MS source. The synthetic test results show that equivalent location results are obtained by the multiscale waveform inversion and single fine-grid-based waveform inversion. However, the average computational cost of the two-stage multiscale waveform inversion is only about 30% of the single fine-grid-based waveform inversion, indicating the high efficiency of the multiscale waveform inversion. In addition, the second-order quasi-Newton iteration algorithm (i.e., the L-BFGS method), with higher efficiency than the first-order gradient-based iteration method, was used to reduce the storage of the huge Hessian matrix and improve the convergence rate during iteration. In subsequent research, multiscale waveform inversion with more stages and grid levels can be developed to further accelerate convergence and decrease nonlinearity of the inversion problem.
5.3. Source-Time Function Estimation by the Fractional-Order Gaussian Function
In a seismic waveform inversion, the symmetrically distributed Gaussian function with a second-order derivative (the Ricker wavelet) or the Gaussian function itself is usually treated as the source-time function modeling. However, the recorded waveforms in a mine region are very complex due to the effects of strong waveform attenuation and scattering at a long propagation distance. Wang [37] put forward the concept of the fractional-order Gaussian functions, and he discussed the influences of two control parameters, namely, the time fractional order derivative and dominant frequency, on shapes of the Gaussian functions and on fitting performance with recorded waveforms. Different from the optimal control parameters calculated by analytical expression from Wang [37], this study utilized a two-parameter grid search to fit more complicated waveforms. The results demonstrate that fractional-order Gaussian functions with a unified value or an independent for each sensor recording both can fit observed windowed waveforms very well. The fitting performance of the fractional-order Gaussian function with an independent for each sensor is slightly higher than the one using a unified derivative for all recordings. The time derivative of the fractional order is mainly related to absorption and attenuation effects of waveform propagation. If the physical meaning of the value is fully considered, it will greatly increase the complexity of the structural parameter modeling (need for accurate modeling and for inverting the value). A purely elastic model was used in this study for wavefield modeling, which mainly considered the relationship between wavefield attenuation and dominant frequency of recording as the signal propagation distance varies. Therefore, a unified fractional order was selected, but the dominant frequency for each sensor recording was adjusted to realize an optimal overall waveform fitting. In the future, for more complex waveforms, we will introduce a more elaborate attenuation model for modeling viscoelastic wavefields, as well as consider varying the fractional order and the dominant frequency for each recording, so as to better simulate the source-time function for synthetic waveforms.
6. Conclusions
This study successfully introduced the multiscale waveform inversion-based location method into the mine MS monitoring field with 3D complicated velocity structures. Furthermore, wavefields were modeled using the highly accurate spectral element method and the high-resolution 3D tomographic velocity model obtained by a 3D ray-tracing method. The SEM wavefield modeling overcomes the difficulties of traditional ray-tracing-based modeling methods in complicated broadband waveform simulation, namely, multipath, focusing, and defocusing finite-frequency effects. Compared with reverse-time migration location methods, the waveform inversion method greatly reduces the requirements for the number and distribution of sensors to iteratively improve the location accuracy. The optimal dominant frequency and fractional order derivative of a generalized Gaussian function are determined through a two-parameter grid search, and the synthetic waveforms based on this generalized Gaussian source-time function can fit the observed waveforms very well. This estimation method provides a new idea for precisely calculating the source-time function and a key to model a complicated wavefield. A more elaborate wavefield modeling approach needs to consider the physical significance and relationship between viscoelasticity and the fractional order derivative. In the L-BFGS iteration of the FWI, the fitting misfit decreases rapidly. By utilizing waveform data with multiple frequency bands step by step, the multiscale waveform inversion strategy can obtain location results equivalent to those obtained by the FWI with a single fine grid. Therefore, this multiscale strategy can significantly reduce iteration times and greatly shorten times for wavefield forward modeling. The average location errors of the eight premeasured blasting events in the Yongshaba mine from the proposed multiscale FWI method is only 17.6 m, which is much smaller than the location errors based on uniform velocity models in previous studies. Furthermore, the location results obtained here are also better than the ray-tracing-based location results when using the same 3D velocity model. In conclusion, the highly accurate location inverted from the multiscale broadband FWI can provide an important basis for further geological analysis, such as focal mechanisms of MS events and displacement characteristics.
Data Availability
The blasting event data used to support the findings of this study have not been made available because we signed a confidentiality agreement.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (NSFC) under Grant 41904083, Grant 52004041, and Grant 41804045; the China Postdoctoral Science Foundation under Grant 2020T130757, and the Basic Scientific Research Operating Expenses of Central Universities under Grant 2020CDJQY-A048 and Grant 2020CDJQY-A047.