#### Abstract

The purpose of this review is to discuss the development of the methods in microseismic/acoustic emission (AE) source localization and detail the principles of the mainly developed methods. The important applications of the microseismic/AE in engineering practice and the history of the source localization (from the initial Geiger method of localization to the later double-difference method) are introduced briefly in the beginning of the review. The factors that influence the accuracy of source locating are discussed in Section 2 to provide references on how to improve the accuracy of the localization methods. Then, the main localization methods in the history of the source and AE research have been classified into four parts (i.e., the localization methods for single event, the joint inversion methods, the relative localization methods, and the localization methods in spatial domain); each part is extended and systematically reviewed, respectively, from Sections 3 to 6, where the merit and demerit of the methods are discussed. Finally, a brief summary and some aspects of prospective research are presented.

#### 1. Introduction

Seismic source localization is one of the most classic and fundamental problems in seismology, and it is of great significance for researching the basic problems of seismology and rock mechanics (e.g., the structure of seismic activity, the internal structure of the Earth, the geometry structure of the seismic source, the response of microseismic to the rock mass, rock support under seismic loads, and the evolution of energy in the rock mass) [1–6]. Furthermore, fast and accurate source locating is also crucial for rescue and relief work after the earthquake disaster and mine earthquake. The localization methods were also widely applied to the location of microseismic sources in mines, defects in materials, and AE sources of indoor experiments [7–9], which can effectively analyze the evolution law of the high stress area and explore the propagation of material cracks [10, 11]. Therefore, fruitful results have been achieved in the field of microseismic/AE source localization, and a series of effective localization methods have been developed.

Due to the limitation of experimental equipment and computational tools, the generally used method to locate the sources in the early days was the geometry graphing method which confined the dissemination of the microseismic/AE technique. The application and optimization of microseismic/AE source localization methods have been advanced greatly as a result of the rapid development of the computer technology (e.g., the Geiger localization method was proposed in 1910 [12, 13], but it was not widely used for source locating until the 1970s [14]). Meanwhile, new localization methods (e.g., HYPOINVERSE [15], HYPOCENTER [16], simplex method [17], genetic algorithm [18], joint inversion methods [19–21], and double-difference localization method [22]) have gradually emerged, of which the mainly used methods are detailed in this paper.

This review is arranged in seven sections: After the “Introduction,” the factors influencing the locating accuracy are discussed in Section 2 to provide references for improving the accuracy of the localization methods. Four groups of localization methods are critically assessed from Sections 3 to 6. The localization methods for a single event, including analytical methods, linear iterative methods (Geiger method and its optimized methods), and nonlinear methods, are reviewed in Section 3, followed by the joint inversion methods in Section 4 where the methods considering the correction of station and velocity structure are discussed to reduce the effects caused by the error of P-wave picking and the velocity structure model. Then, the relative localization methods for earthquake swarm are reviewed, and the principle of the methods is presented. The localization methods in spatial domain are briefly reviewed in Section 6. Finally, a brief summary and some perspectives are presented.

The objective of this review paper is to (1) reveal the factors influencing the locating accuracy; (2) detail the merit and demerit of different localization methods; (3) provide some proposal to improve the location accuracy; and (4) discuss the current challenges facing the field of microseismic/AE source localization.

#### 2. Factors Influencing the Locating Accuracy

All the microseismic/AE source localization methods follow the basic distance-time constraints and can be expressed aswhere *t* is the moment when the sensor receives the signal, *t*_{0} is the time when the event happens, *s* is the wave propagation path from the source to sensor, and is the wave velocity on the propagation path which is a path-dependent function. Each localization method is based on these constraint conditions using a variety of algorithms to obtain the source coordinate values. The error of calculated ∆t can be written as a function as follows:where , , and are the errors caused by the arrivals, velocity, and methods employed and is the total error in the locating system. It can be seen from the formula that the determination of the arrival time, the wave velocity, and the localization algorithm is crucial to the locating result. The impact of these factors is analyzed as follows.

##### 2.1. Noise

The microseismic signal is usually of small magnitude, which is different from the normal earthquake signal. In this case, the background noise signal can greatly affect the waveform of the microseismic signal and decrease the utility of the microseismic data. Therefore, it is extremely essential to preprocess the microseismic signal to separate the noise from the microseismic signal [23]. Forghani-Arani et al. [24] suppressed surface-wave noise in microseismic data by conducting signal-noise separation in the τ-p domain. Due to the distinct characteristics that microseismic signal and noise show in the τ-p domain, the noise can be separated from the microseismic signal. Chen [25] proposed an algorithm which uses the least squares inversion method to decompose the signal into a set of components; then, the noise can be filtered out based on the predictive property of the microseismic event along the time direction. Furthermore, the methods using mathematical morphology to detect the microseismic signal were also proposed by Huang et al. [26] and Li et al. [27, 28], which decompose data set into multiscale components and obtain more information by detecting more waves. Dong et al. [29] proposed two indicators using the origin time of blasts and established some discriminant models of blasts and seismic events in mine seismology, which has a reasonably good discriminating performance for classifying between events and blasts in mines.

##### 2.2. P-Wave Arrivals

The arrival data of the P-wave reaching the sensor are widely used to locate the microseismic/AE sources. Picking the arrivals of the P-wave triggered by events is the most fundamental and important step in source localization. Moreover, the impact of the error arrivals on localization results cannot be eliminated by algorithms. At present, the widely used identification method is STA/LTA algorithm [30–32]. The principle of this algorithm is to use the ratio of short-term signal average to long-term average to reflect the change of signal amplitude and energy. Gibbons and Ringdal [33] proposed the arrival identification method of P-wave based on the cross-correlation technique and used the RAM seismogram method to determine the P-wave arrivals. Based on DWT and STA/LTA methods, Li et al. [34] used the kurtosis (WS/LK) method to pick up P-wave arrivals precisely and automatically. This method can accurately identify P-wave signals in the presence of noise, with 97.5% events identification accuracy within 15 ms. Taken together, the accuracy of arrival-time picking up is the crucial basis for precise source localization, and the improvement of the arrival pick-up method greatly enhances the efficiency and accuracy of source localization. Furthermore, with the development of machine learning algorithms, the picking methods using machine learning algorithm were proposed [35, 36], which can pick the arrivals automatically even in the case of moderately strong background noise.

##### 2.3. P-Wave Velocity

The impact of velocity can be considered as two aspects, one is the simplification of the velocity structure and the other is the impact of P-wave velocity measurement. For most of the localization methods, the velocity structure is assumed to be a constant and the propagation path of the P-wave also becomes a straight line; then, equation (1) can be written aswhere () are the coordinates of the sensors and () are the coordinates of the source. It is obvious that the simplification of velocity structure can lead to a huge error while the velocity structure is complex.

The premeasured P-wave velocity is generally used in many source localization methods. In practical application, the P-wave velocity is usually measured by conducting indoor tests on the rock samples taken from the site. The disadvantage of this method is that the rock samples cannot represent the wave velocity characteristics of rock mass under actual geological conditions, so relatively big error still exists. Another method is to calculate the wave velocity of the rock mass by monitoring the data of the active blasting event, as shown in Figure 1, which is of engineering significance to a certain degree. However, the blasting event will create empty zones and plastic zones in the rock mass, which can change the wave velocity structure of the rock mass and result in a greater locating error.

The dynamic changes of the wave velocity in engineering practice have greatly confined the application of the premeasured P-wave velocity localization methods. Therefore, the localization method without the premeasured P-wave velocity was developed. By simplifying the nonlinear equations to linear equations, Dong et al. [37] developed an analytical localization method to solve the coordinates of the sources, in which the necessity of the premeasured P-wave velocity and the square root operations are avoided. To compare the accuracy of the locating methods with or without the premeasured velocity, Dong et al. [38] analyzed the locating accuracy using the virtual position test by considering the velocity with small error rates of 1%–5% floating. Results show that the floating of the velocity resulted in large errors of the source coordinates. The result reflects that localization methods without the premeasured velocity are likely to be the main direction of source localization.

##### 2.4. Localization Algorithm

Algorithm means an unambiguous specification of how to solve a class of problems. For the microseismic/AE source localization area, the locating accuracy is the most important criterion to evaluate the pros and cons of the localization algorithm. Inaccurate localization algorithm has no engineering application value. In addition to the requirements for localization accuracy and stability, the proposed real-time localization also puts forward certain requirements on the calculation time of the localization algorithm. The speed of source localization is also of great significance in engineering practice. Dong et al. [39] used the famous optimization algorithms including Levenberg–Marquardt (LM) algorithm, simplex method (SM), quasi-Newton method (QN), maximum inherit optimization (MIO), self-organizing migrating algorithms (SOMAs), and global optimization coupled with MIO, QN, and LM to locate the events of Dongguashan copper mine in China. The locating errors of algorithms are compared and listed in Figure 2. The result shows that the SOMA localization algorithm is the most stable, while MIO, QN, and LM are relatively easy to fall into the local optimal solution. It is noteworthy that source localization algorithms are formed in different engineering backgrounds with specific application conditions, so the appropriate localization algorithm should be selected based on the actual engineering environment rather than the accuracy.

#### 3. Localization Methods for Single Event

##### 3.1. Analytical Localization Methods

The coordinates and arrival time of each sensor are used in the analytical localization method, to derive the analytical solutions of AE with mathematical formulas. In 1928, Inglada [40] proposed an efficient and explicit method to solve the location equation, whose advantage lies in that fewer sensors (4) are needed to locate the AE. In addition, there is a unique analytic solution in the central region of the sensor array. Essentially, the analytical method is also used to obtain the localization result in the Inglada method. Although the unique solution can be determined, the wave velocity needs to be measured in advance. Christy [41] performed an iterative optimization of the Inglada method and achieved better localization results than noniterations. In 1968, Mogi [42] of the Japan Institute of Geological Survey carried out the transverse bending test of long strip samples with premeasured wave velocity, where two and four sensors were used for 1D and 2D localizations, respectively. It can be classified as a localization method of arrival-time difference due to the small number of deployed sensors.

In 1970, Leighton and Blake [43] of the US Bureau of Mines developed an analytical localization method with simple steps, namely, USBM. The arrival times and coordinates of 4 sensors are recorded to achieve better localization accuracy, which is widely applied in the localization of MS/AE. In 1974, Blake [44] proved the uniqueness of the USBM solution. Godson and Bridgs [45] designed a 32-channel localization system, where USBM were applied to locate AE with high accuracy. Strang [46] conducted an iterative study of USBM and obtained a better localization effect compared to the noniterative method.

Smith and Abel [47] developed a new analytical localization method using the sphere interpolation and least squares. Duraiswami et al. [48] proposed a set of analytical localization methods based on time difference localization. The localization solutions can be optimized using the L1 norm when the number of sensors is more than 4. Brandstein et al. [49] proposed a three-position method based on linear interpolation.

All the analytical localization methods mentioned above share the same disadvantage: the wave velocity needs to be determined in advance. In fact, in many practical engineering applications, the real-time wave velocities are unknown or difficult to measure, which will cause great temporal and spatial errors due to the dynamic velocities and complex environment. Focusing on the vital issue, many scholars have proposed a series of optimized analytical localization methods to eliminate the effects of premeasured velocity on localization accuracy. Dong et al. [50] proposed the 2D and 3D analytical localization methods without premeasured velocity. The main methodology of 3D analytical localization method is summarized as follows:

The governing equation for the coordinates of the AE source *P* (*x*, *y*, *z*) is shown in (4), where *t*_{i} (*i* = 1, 2, 3, 4, 5, 6) is the arrival time corresponding to the sensor *S*_{i} (*x*_{i}, *y*_{i}, *z*_{i}). *t*_{0} is the origin time of the AE source. The symbol represents the average P-wave velocity in the travel path.

By taking difference between the equation with *i* = 1 and the others (*i* = 2, 3, 4, 5, 6), we can obtain the following equation:where *m *= (2, 3, 4, 5, 6) and .

By substituting *V* and *S* for and *Vt*_{0}, respectively, (5) can be transformed into the following equation:

Equation (6) can also be rewritten aswhere

Thus, the source coordinates (*x*, *y*, *z*), average wave velocity , and origin time *t*_{0} can be obtained by solving the proposed linear equations with five unknowns.

The above method is under the condition that the number of sensors is 6. In fact, it is common that greater than 6 sensors are used to improve the locating accuracy in practical engineering. Therefore, it is a vital problem to take full advantage of the rest sensors to obtain a more accurate solution. Dong et al. [37, 50] proposed the multisensor analytical localization method based on the logistic probability density function. The proposed analytical localization method was applied to locate and verify the blast tests in the Dongguashan copper mine (China). Results show that the effect of premeasured wave velocity is eliminated and the location accuracy is improved significantly.

The basic methodology of this multisensor analytical localization method is explained as follows: 6 sensors are selected randomly from *m* triggered sensors to combine groups of analytical solutions. Figure 3 shows an example of the combination process, where *m* and *n* are both equal to 7. Then, the logistic probability density function is applied to fit the whole solutions; the coordinates of the AE source are exactly the abscissas corresponding to the maximum value of the logistic probability density function. Figure 4 clarifies the localization process and highlights the proposed 3D analytical localization method.

The basic methodology of this multisensor analytical localization method is explained as follows: 6 sensors are selected randomly from *m* triggered sensors to combine groups of analytical solutions. Then, the logistic probability density function is applied to fit the whole solutions; the coordinates of the AE source are exactly the abscissas corresponding to the maximum value of the logistic probability density function.

##### 3.2. Geiger Method

Current mostly used localization methods are usually based on the classic method proposed by Geiger in 1910 [12, 13], while until 1970s, with the rapid development of the computer technology, the Geiger method was applied to microseismic/AE locating field [51]. The main idea of Geiger method is to linearize the problem, make the residual of arrivals reach the minimum or a certain accuracy (by using iterative computations), and finally obtain the coordinates of the source. The specific steps are listed as follows:

For the arrivals data received by the AE sensors, construct the following function:where is the residual between the observed arrivals and calculated arrivals . is the quadratic sum of the arrivals residual which represents the fitting degree between the hypothetical source and true source, a smaller Φ means that the data we set are more fitted to the actual measured data.

In addition to the Geiger method, many scholars also proposed the localization methods based on the least squares method. Twenty-two significant AE events during rock failure were located using the least squares method based on the time difference between the arrival of S-waves and the six sensors [52]. The location results were used to analyze the propagation of microcracks during rock loading. Fedorow [53] described the least squares method, introduced a consistent estimator with normal distribution and solved it iteratively, and applied this method to AE localization.

##### 3.3. Optimized Algorithms Based on the Geiger Method

The Geiger method of seismic source localization is an iterative method; the iterative process usually converges rapidly unless the data are badly configured or the initial guess is very far from the true solutions. However, it also happens that the solution converges to a local minimum and this would be difficult to detect in the output unless the residuals are very bad. Many scholars have proposed optimized algorithms based on the Geiger method: Buland [54] combined the QR algorithm with the Geiger method to solve the solution. Aki and Lee [55] proposed a series of improved schemes based on the Geiger algorithm, including parameter separation, joint inversion of the 3D velocity structure, and the source and coupling velocity with source resolution to solve the problem. Thurder [56] expanded the residuals formula with the second-order Tailor expansions and solved the coordinates of the epicenter. The second-order Tailor expansions were used to improve the convergence and stability of the algorithm, but the amount of computation also increased.

In 1975, Lee [14] released a series of computer programmers for determining hypocenter, magnitude, and first motion pattern of local earthquakes named HYPO71 and HYPO78-81 written in FORTRAN. Instead of carrying out the traditional procedure, a stepwise multiple regression in used in the HYPO71. A statistical analysis is first performed to see which independent variable should be included in the regression, and the normal equations are then set up for only those significant variables. Therefore, the adjustment vector is obtained by solving a matrix which is never ill-conditioned. Furthermore, convergence to a final hypocenter solution is also more rapid. The HYPO71 program also requires considerable efforts: accurate station coordinates, reasonable crustal structure model, and reliable P and S arrivals. Herrmann et al. [57] relocated the Denver earthquakes of 1967-1968 using HYPO71.

Byerlee [58] applied the Geiger method to the fluid injection test of rocks where six sensors were used to receive the AE signals and establish the 3D location equations with multiparameters.

Klein [15] developed the HYPOINVERSE program to obtain the coordinates and the magnitude of the hypocenter. Lienert et al. [16] combined features of the two well-known methods HYPO71 and HYPOINVERSE, with a new technique—adaptive damping—and proposed a localization method—HYPOCENTER. Each column of the linearized condition matrix *T*, which relates changes in arrival time to changes in the hypocenter position, is centered and scaled to have zero mean and a norm of one, solved iteratively by adding a variable damping factor, , to their diagonal terms before inversion. If the residual sum of squares increases, return to the previous iteration, increase , and then try again. This procedure, which we called adaptive damping, always results in residuals which are less than or equal to the HYPO71 and HYPOINVERSE residuals. HYPOINVERSE fails to converge when the true depth is greater than about 40 km, while HYPO71 fails at depths of greater than 50 km. HYPOCENTER has a better performance on the difference in true versus calculated coordinates, especially depth.

Nelson and Vidale [59] improved the HYPOCENTER method and presented a new method for locating in a region with arbitrarily complex three-dimensional velocity structure called QUAKE3D. The method is a grid search method which searches all possible hypocenters and occur times and finds the global minimum travel-time residual location within the volume. It is obvious that the criterion is much important in the grid search method. QUAKE3D employs the L1 criterion (absolute value) instead of the commonly used L2 criterion (least squares method) based on the fact that L1 criterion is more robust than L2 criterion while the station coverage is sparse [60]. The HYPOINVERSE and QUAKE3D are applied to locate the source of earthquakes occurred in Bear Valley from 1977 to 1986, which suggests that the QUAKE3D has more accuracy in locating the result than HYPOINVERSE.

##### 3.4. Nonlinear Localization Methods

The above localization methods are all based on linear methods. In recent years, the nonlinear method theory has become the frontier in the fields of natural science. Since most of the geophysical problems are nonlinear problems, nonlinear methods tend to be more realistic than the linear methods while solving the geophysical problems. Various nonlinear optimization methods have been developed rapidly, including the methods based on derivatives such as the steepest descent method, Newton method, and conjugate gradient method and the nonlinear methods based on nonderivatives including Monte Carlo method, genetic algorithm (GA), simulated annealing (SA), and random search and simplex search algorithm. Moreover, the waveform-based source localization method based on wave equation is also a nonlinear localization method, which can obtain both the source location and velocity inversion information [61]. Due to the rapid development of computer technology, nonlinear methods have played an important role in geophysical inversion.

The essence of nonlinear localization methods is solving the least squares problem based on equation (1). The principles of some nonlinear localization methods are listed in the following.

###### 3.4.1. Localization Methods Based on Powell’s Method

Powell’s method [62] is a direct method of searching for the minimum value of the objective function. The method does not require partial derivatives or inverse matrices, has low requirements on the initial iteration value, and has good adaptability. The basic principle is to divide the whole calculation process into several stages. Each stage (one iteration) consists of one-dimensional search. At each stage, we search first along the known *n* directions to get the best point and then search along the line connecting the initial point and the best point of this stage to find the best point. After this, the last search direction is used to replace one of the first *n* directions to start the next stage until the calculated residual value is less than the given allowable error or the number of iterations has reached the restraint. Many scholars have applied Powelll’s method for locating the earthquake [63–65].

###### 3.4.2. Localization Methods Based on the Monte Carlo Method

Monte Carlo method was first proposed by Metropolis and Ulam [66]; the method does not search thoroughly in the model space but searches randomly compared to the exhaustion method. Practice shows that if we randomly select the model in the model space and find the global minimum of the objective function, we will save a lot of time compared with planning the space of model space to calculate the global minimum of the model, but the workload is still huge and cannot guarantee the minimum value we found is the global minimum value. Meanwhile, the inherent randomness of the Monte Carlo method can lead to the failure of the calculation. Billings et al. [67] applied the Monte Carlo method to investigate the effect of picking errors.

###### 3.4.3. Localization Methods Based on the Simplex Method

In 1962, Spendley et al. [17] proposed the simplex method which is widely used as a mathematical tool in a variety of engineering fields. As a geometric search method, the simplex method has a property that each iteration is better than the previous one, so the optimal solution can be obtained only by repeated iterations. Meanwhile, if there is no optimal solution, the simplex method can also be used to judge. Nelder and Mead [68] proposed a search iterative method based on the work of Spendley in 1965.

Prugger and Gendzwill [69] introduced the simplex method into the seismic location and got a satisfactory localization result. This method avoids the derivative operation and the matrix transpose operation which greatly reduces the computational complexity and adapts L1 norm to reduce the influence of residual on arrival. Ge [70] evaluated Prugger’s method and refined it. Zhao et al. [71] used the simplex method to locate the earthquakes in Tibet.

###### 3.4.4. Localization Methods Based on the Genetic Algorithm

Genetic algorithm [18, 72] is also a nonlinear global optimization method; its basic idea is imitating the genetic process of the biological world. The first step in using a genetic algorithm is to determine the encoding of the problem parameters, usually encoding of the parameters in binary. For earthquake location, the parameters are (). The upper and lower bounds of the parameters are coded, and a set of individuals are randomly generated, called population. The residual between the calculated time and the actual observed time is used as the fitness function. The smaller the residual is, the higher the survival probability of the individual is and so is the probability of being the parent. Offspring can be created by crossover, and a certain probability of mutation is introduced to enrich the diversity of the population. The above process is repeated for the offspring obtained until the stop rule is satisfied and the individual with the highest fitness function is obtained, which is the optimal hypocenter parameter.

The advantage of genetic algorithm lies in that the process of solving is only related to the object. It only needs to carry out simple operations such as crossover and mutation without complicated mathematical operations such as calculating derivatives and has good global optimization ability. The genetic algorithm has been widely used in the field of source location and geophysical inversion [73–79].

##### 3.5. Microseismic Source Location Method without the Premeasured Wave Velocity (MSLM-WV)

It is discussed in Section 2.2 that the premeasured wave velocity can result in a large error in source locating, which indicates that the localization method without premeasured P-wave velocity can obtain a more accurate result. Dong et al. [80] presented a microseismic/AE (MS/AE) source location method without the need for a premeasured wave velocity. Tests of both the pencil lead break and the thermal fracture in granite were carried out to estimate the applicability and accuracy of the method. Results show that the location accuracy of the proposed method is significantly improved, which is superior to the results of the traditional location method (TM) using premeasured wave velocity. The merit and demerit of the TM and MSLM-WV are listed in Table 1.

##### 3.6. Multistep Localization Method without the Premeasured Velocity (MLM)

Based on the localization function with the model of arrival-time difference, Dong et al. [81] proposed a multistep localization method without premeasured velocity, named MLM. The velocity interval can be narrowed and optimized continuously according to the obtained minimum and maximum velocities, which can be solved in each localization process. The optimal velocity can be determined until the velocity differences are less than the threshold, where the localization results corresponding to this velocity interval have the highest location accuracy. MLM can improve the location accuracy and computation efficiency in the complex environment compared to both the traditional localization method (TLM) and TD method [38]. Figure 5 is the comparison of the MLM, TLM, and TD method.

##### 3.7. Collaborative Localization Method Using Analytical and Iterative Solutions (CLMAI)

The abnormal arrival of a certain sensor is another main factor that affects the localization accuracy. Note that a stable solution with high precision can be obtained using the analytical localization method when the accurate data are provided. Dong et al. [82] proposed a collaborative localization method using analytical and iterative solutions to solve the source coordination, in which the analytical localization method without premeasured wave velocity was applied to remove the abnormal arrivals and the iterative localization method without premeasured wave velocity was used to solve the source coordination with the clear arrivals. It can effectively eliminate not only the error caused by the abnormal arrivals but also the error caused by the difference between the measured wave velocity and the dynamic wave velocity. The proposed CLMAI was verified using the recorded data of blasts and microseismic events in Kaiyang Phosphorous Mine, China, where most of the mining works mainly concentrated in six mining areas [83], and a 32-channel microseismic monitoring system, including 26 single-component sensors and 2 three-component sensors, was established to avoid the destructive disasters [82]. The locating results of the CLMAI are compared with the results of LM-PV and TD. The compared result is shown in Figure 6. It shows that the localization results of CLMAI are obviously better than the results of the LM-PV and TD.

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. Joint Inversion Methods

The localization methods of single event shown in Section 2 mainly concern about the optimization of the method while ignoring the influence of the arrivals and velocity structure. It is noteworthy that the simplification of the wave velocity structure is an important factor affecting the locating accuracy. The effect of arrivals can be corrected by adding a station correction parameter to equation (1). The error caused by the simplified wave velocity structure can be solved by joint inversion of source parameters and wave velocity structure. The joint inversion problem is no longer about the location of a single source, but the location of the earthquake area. Two main methods of joint inversion are described in the following.

##### 4.1. Joint Inversion of Source Location and Station Correction

Cleary and Hales [84] found that equation (1) should consider a term for the correction of the station to describe the difference between the observed travel time and the travel time recorded in travel time table. A joint hypocenter location method–joint epicenter determination (JED)—was first proposed by Douglas in 1967 [19]. The travel time residual can be obtained by applying the station correction to the arrivals:

The coordinates of the events and the correction of the stations can be obtained by applying the standard least squares technique to equation (10).

Dewey [20] proposed a modification of Douglas’ method called joint hypocenter determination (JHD) and applied it to relocate the earthquakes happened in western Venezuela [85]. For a locating system, equations can be obtained where *m* is the number of events and *n* is the number of stations. It should be noticed that the number of equations is always so large that the efficiency of the JHD can be very low. To solve the problem of oversized matrices due to the large number of AE events and stations, the PMLE method using parameter separation based on JHD was proposed by Pavlis and Booker in 1983 [86]. The PMLE requires less computation and still has strong stability. Pujol [87] discussed the relationship between the JHD method and the rough solution of Frohlich [88] to optimize the PMLE.

The JHD method and parameter separation method was applied to the monitoring data of the Kunming seismic network to obtain the corrected P-wave arrivals [89]; the result shows that an accurate locating result can be obtained by using the corrected arrivals. Ratchkovsky et al. [90] applied the JHD to relocate 1604 earthquakes occurred from 1988 to 1993 in the Wadati-Beinioff zone. Guo et al. [91] used the Geiger method and the JHD to locate the earthquakes of Jiashi area in Xinjiang Region, which suggests that the JHD has a much higher locating accuracy than the Geiger method.

##### 4.2. Joint Inversion of Source Location and Velocity Structure

The joint inversion of source location and wave velocity structure was first proposed by Crosson [21]. This method uses the optimized least squares method [92] and inverts the wave velocity structure as an unknown parameter with source parameters simultaneously to reduce the errors caused by the simplified wave velocity structure. The parameters of an earthquake, as well as the information about the velocity structure, can be obtained by this method. Even if there are low-velocity zones in the velocity structure, the inversion results can still be reasonable.

#### 5. Relative Localization Methods

Absolute localization method can get accurate locating results when the structure of wave velocity is known or the velocity is under good condition. However, when the geological structure or the structure of the wave velocity is complex, solving the solution with a uniform wave velocity structure will cause great locating errors. Therefore, under such conditions, it is necessary to eliminate the influence of the complex velocity structure. In seismology, the relative localization method is an accurate method under the condition of complex wave speed.

##### 5.1. Master Event Localization Method

The master event localization method was developed from the JED. The principle of the method is that when the distance between two events is far less than the distance between the event and the station, it can be considered that the time difference between the two events and the stations is determined by the relative distance and the wave speed between them. Therefore, it can eliminate the influence of complicated wave velocity structure between source and station. Spence [93] put forward the master event localization method based on multievent joint location method. In this method, one specially well-located event called master event is selected and the location of a group of earthquakes around it is calculated to determine the source location of these events. The method does not need iterative computation; however, the accuracy of the method depends on the selection of the master event.

Ma [94] used the master event localization method to relocate the focal position of the Huoshan earthquake swarm in 1973 in Anhui province. Poupinet et al. [95] used the method to locate the earthquake of the Calaveras Fault in California. Hu et al. [96] developed a relative localization method for AE in a laboratory. This method has less dependence on velocity structure and is of great significance when applying to the AE locating in rock tests which is nonintegral and anisotropy media.

##### 5.2. Double-Difference Method

If the focal mechanisms of the two earthquakes are similar and closely spaced, the propagation paths and the waveforms recorded on the same station are similar. Through the use of waveform cross correlation, the time difference can be accurate to the millisecond travel time and the relative errors between the two earthquakes can be reduced to tens of meters [97]. Waldhauser and Ellsworth [22] proposed a double-difference method (DD) and developed a localization program—hypoDD [98]—based on this method in 2000. The method was applied to the locating of earthquakes in the Northern Hayward Fault of California.

The double difference between the station *k* and event *i* can be defined as

The following equation can be obtained by applying the first order of Taylor expansion:

Equation (12) can be written as

For all stations and events, the following matrix can be obtained:where *G* is a matrix of size (*p* is the number of double-difference observations and *q* is the number of events) containing the partial derivatives, *d* is the data vector containing the double differences, *m* is a vector of length *4q*, contains the changes in hypocenter parameters we wish to determine, and *W* is a diagonal matrix to weight each equation. The location of the events can be obtained by solving equation (14).

The difference between the double-difference method and other methods is that no station correction term is needed to eliminate the effect of the velocity structure. The outstanding advantage of the double-difference method is that it can use the cross-correlation analysis of waveform to pick the arrival time of the event and greatly improve the accuracy of the data. The double-difference localization method inverts the relative position of each earthquake in a cluster of earthquakes relative to the centroid of the cluster, which is much different from the master event localization method with greatly improved applicability. In addition, the anti-interference and robustness of the double-difference localization method are also strong. The double-difference method has been widely used in the locating of earthquakes [99]. Another multievent localization method is PMEL [100], which is a grid search method. The JHD, HDC, DD, and PMEL were applied to the locating of the earthquakes in Izmit/Duzce, and the results were compared [101].

#### 6. Localization Methods in Spatial Domain

The above methods are all localization methods based on the time domain. There are many deficiencies using the localization method in the time domain. For example, the moment of occurrence is always unknown, the joint significant problems exist between the parameters of the source, and the localization results seriously depend on the velocity structure and the distribution of the network. In order to overcome the above shortcomings, the localization methods based on the spatial domain have been developed. In 1977, Lomnitz [102] proposed a method by replacing the arrivals residual with the distance residual. The equation in this method only involves the location of the epicenter, the depth of the source, and the time of occurrence, which avoids the mutual influence of the parameters. Carza et al. [103] proposed a source localization method without travel time data and used this method to locate the teleseism. The results show that the error of this method is 8–20 km. The method does not need to make a priori assumption of travel time and can be used for the inversion of structures inside the Earth.

In 1957, Romney [104] proposed a method for the locating of the near earthquake and established the distance-residual equation by using the arrival-time difference and surface average apparent velocity of two stations adjacent to each other. The number of needed conditions is small, and it is easy to solve the equation. The locating results have little dependence on the structure. However, the determination of focal depth and time is not well solved. Zhao and Zeng [105] made improvements on this point by using the slope of the arrival curve to determine focal depth and using the intercept of the arrival curve on the time axis to determine the occurrence time.

#### 7. Suggestions and Conclusions

This review of research since the pioneering work of the 1910s is intended to present the state of the localization methods in both microseismic and AE sources. It appears that from the initial single-event, linear algorithm to the later multievent, nonlinear algorithm, the size, complexity, and locating accuracy of hypocenters are constantly improving. This section presents concluding remarks and some prospects requiring investigation.

For single-event localization methods, the velocity is the main effect that influences the locating accuracy. It means that the analytical method shall be applied to the system where the structure of P-wave velocity is isotropy (e.g., metal and intact rock). Meanwhile, the analytical methods can be applied to filter out the abnormal arrivals on account of that the analytical methods are sensitive to the input of arrivals.

The nonlinear localization methods can obtain more accurate results than that based on the linear methods, while the main constraint of its application is the calculation efficiency. Note that the development of computer technology plays a very important role in promoting the evolution of the localization method; more-efficient nonlinear localization methods would be proposed in the future.

For the multievent localization methods (i.e., joint inversion methods and relative localization methods), the coordinates of source swarm can be obtained; however, the defect of the methods is that the methods can only be applied when the events are close to each other so that they can be regarded as a swarm.

As for the localization methods in spatial domain, it provides a new idea for source locating but the accuracy is the main constraint of the methods. More investigation is needed to develop the methods and expand its application.

It should be noted that various localization methods are formed in different engineering backgrounds, which leads to a big difference in the way of the construction, processing, and evaluation of the objective function, so the appropriate localization method should be selected based on the actual background for experiments and engineering rather than the accuracy.

Through the above comprehensive review of microseismic/AE source localization, it can be realized that the location accuracy is significantly affected by some factors, which mainly includes the identification and picking of P/S-wave arrivals, P-wave velocity structure, and P-wave travel path.

As for the P/S-wave arrivals, it is evident that their identifying accuracy needs to be improved, which is the most direct and effective approach to guarantee the location accuracy. Therefore, we can merge some more-advanced algorithms into the current algorithms, to develop several joint identification and picking algorithms, which are beneficial for improving identifying accuracy and computational efficiency. On the other hand, the comprehensive verification of P/S-wave arrivals through multiple identifying methods is also another way to ensure the location accuracy. It is common that an identifying method may be unstable in some situations, which means that it is necessary to apply other identifying methods to verify the previous identifying results. Thus, only the P/S-wave arrivals satisfy all the identifying methods can be inputted for source localization.

The P-wave velocity structure has close relationship with P-wave travel path. In this paper, many localization methods can improve location accuracy, whereas they still assume the P-wave travel path as a straight line. This assumption will only adapt to the monitoring environment containing single propagation media, which is obviously improper for the actual engineering. In the complex rock mass structure, the authentic P-wave travel paths cannot simply be a straight line. For example, the underground goafs generated in the mining process can significantly change the P-wave velocity structure of rock mass and then influence the P-wave travel paths of microseismic sources. Therefore, the location accuracy can be improved if the P-wave travel paths can be indicated by curves, multisegment lines, and the combination of curves and multisegment lines, where the differences between the authentic and calculated P-wave travel paths will be much smaller. For example, we can use the ray-tracing methods to continuously optimize the P-wave travel paths since they can divide the original straight ray into multiple rays according to the actual P-wave velocity structure.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors acknowledge the support of projects (51774327, 51504288) of the National Natural Science Foundation of China, Natural Science Foundation of Hunan Province in China (2018JJ1037), the Young Elite Scientist Sponsorship Program by Hunan Province of China (2018RS3001), and National Basic Research Program of China (2015CB060200).