Two-Dimensional Far Field Source Locating Method with Nonprior Velocity

1 IOT Perception Mine Research Center, China University of Mining and Technology, Xuzhou 221008, China 2School of Mathematics & Physics Science, Xuzhou Institute of Technology, Xuzhou 221008, China 3School of Information and Electrical Engineering, China University of Mining and Technology, Xuzhou 221008, China 4School of Mechanical Engineering, Nanjing University of Science & Technology, Nanjing 221000, China


Introduction
Microseismic monitoring, which was first studied by Obert [1], has been successfully applied in mining [2][3][4][5], oil exploitation [6][7][8][9][10], water conservancy construction [11][12][13], and so forth.One of the core technologies of microseismic monitoring is source locating method, which can be reduced to (1) in homogeneous medium: where , , and  are coordinates of seismic source and   is the initial time of source;   ,   , and   are coordinates of the th sensor; V  is the average wave velocity measured by the th sensor;   is arrival time measured by th sensor.It is obvious that 4 or more sensors at least are needed to solve for (, , ,   ) if V  is known.A set of nonlinear equations as expressed by (1) can be solved either iteratively or noniteratively.The well-known noniterative methods, such as Inglada method [14,15] and USBM (United States Bureau of Mines) method [15][16][17], may not be applicable for current microseismic monitoring system.The Inglada method uses only a minimum number of sensors that are mathematically required for source locating, and because of this requirement, no optimization method can be applied to the algorithm.The USBM method introduces least squares method so that all available arrival time information can be used simultaneously for source locating calculation.However, for its matrix inversion, the method will be unstable while the measured values have gross errors, which is common in a real-world situation.
The main iterative methods for source locating include Geiger's method [18,19], Simplex Algorithm [20,21], and Genetic Algorithm [22].Geiger's method has heavy computing burden and may be divergent because of the foundation of the method-first-order Taylor expansion and least square method.Many scholars studied these problems and have proposed some solutions [23,24].Simplex Algorithm [25], a robust geometry search algorithm, is one of the dominant methods for source locating [26][27][28].Simplex is a geometric figure which contains one more vertex than the solution space's dimension, and the optimization process is simple which is moving the worst vertex till the optimal solution meets the preset conditions.GA (Genetic Algorithm), an optimization method, simulates nature selection in which only the "fittest" solutions survive so that they can create even better answers in the process of reproduction.Although the theory of GA is imperfect, the algorithm was still introduced in earthquake source locations in 1992 [29][30][31][32], and it has shown attractive prospects in parallel computing, such as SPARK [33], because of its intrinsic parallelism.According to the research from Yun and Xi [34] and He et al. [35], EGA (Elitist preserved GA) is global convergent, which provides guidance for GA based algorithms.
All of the algorithms mentioned above have the same hypothesis, which is that the velocity structure is known before location calculation, which is obviously impossible in practice.Therefore, the location accuracy of algorithms mentioned above will be seriously affected by prior velocity.According to Dong et al. 's research [36], the accuracy of methods with prior velocity in homogeneous medium will decrease seriously when the velocity error reaches 1%-5%, while the accuracy of methods with nonprior velocity will not be affected.Consequently, Dong and Li [37] proposed a new location method with nonprior velocity based on arrival time of PS waves; however, the method was not used when P-wave and S-wave could not be distinguished clearly.Li et al. [38] proposed a location method with nonprior velocity based on Simplex Algorithm and the essence of its method is a new objective function without velocity, which could also be used by other optimization methods.
Actually, Li's method is a concrete realization of Prugger's method [21], which will give a correct result when the seismic source is surrounded by sensors and the measured arrival times have no error.However, in practice, the seismic source is not always surrounded by sensors for a variety of reasons, and Li's objective function changes very little near true value which will decrease the accuracy of location by reason of finite word length effect.To solve this problem, this paper proposed a new objective function, which can locate faster and more stable.

Two-Dimensional Source Locating Theory
In this paper, we will discuss only one scenario that frequently happened during mining monitoring as shown in Figure 1.
The parameters of Figure 1 are shown as follows.
( 1 ,  1 ), ( 2 ,  2 ), ( 3 ,  3 ), and ( 4 ,  4 ) are coordinates of equidistant placed sensors, (, ) are coordinates of seismic source, V is the average wave velocity,  12 is the arrival time difference between ( 1 ,  1 ) and ( 2 ,  2 ), and  23 and  34 are similar to  12 .The unknowns are (, ) and V.If V is known, (, ) can be solved by the following equation set: ( If V is unknown, the three unknowns can be solved by the following equation set: Obviously, V in equation set (3) can be eliminated and equation set (3) will change to the following equation set: Equation sets (3) and ( 4) are both nonlinear equations, and they can be solved by Simplex Method, GA, and other global optimization methods.
According to Prugger's method, the objective function of (3) can be defined as (5) and the objective function of (4) can be defined as (6): √ ( According to Li's method, the objective function of equation set (3) can be defined as where If the seismic source is far away from sensors, as shown in Figure 1, neither (6) nor (7) can lead to correct results due to small changes of objective functions near true value.Here is an example.
The graph of objective function (7) under the same conditions is shown in Figure 4.
Enlarged view of Figure 4 near (500, 500) is shown in Figure 5.
As shown in Figures 2-5, two concluding points can be made when the seismic source is far away from sensors as shown in Figure 1; that is, the locating accuracy of  coordinate is better than  coordinate and the objective function, whether (6) or (7), changes little near true value.As a comparison, the graphs of ( 6) and ( 7) when seismic source locates at (1000, 500) are shown in Figures 6-7.
Obviously, the locating accuracy of methods, which take (6) or (7) as objective function, will decline seriously when the seismic source is far away from sensors.This paper presents a new method based on solution space downsizing, and the core of the new method is a new one-variable objective function.Numerical experiments will show that the new method has faster convergence speed, higher accuracy, and better stability.

Two-Dimensional Far Field Source Locating Method Based on Solution Space Downsizing
3.1.Methodology.When locating with nonprior velocity, the locating problem is equivalent to a three-dimensional optimization problem if we choose (5) as objective function, while it is equivalent to a two-dimensional optimization problem if we choose (6) as objective function.Usually, because of dimension reduction of solution space, method taking (6) as objective function will calculate faster than (5) in the same condition.Therefore, if the solution space's dimension of ( 3) or (4) can be reduced to one, the calculation should also be faster.
The asymptotic equation of ( 9) is expressed in The asymptote will be very close to the hyperbola when the seismic source locates far from sensors as shown in Figure 1.If the error between asymptote and hyperbola is less than predefined locating error, (10) can be used to locate instead of ( 9) and we will obtain an approximate equation set instead of (2) as shown in where All the  coordinates are equal when sensors are placed in a line as shown in Figure 1.Therefore,  12 =  23 , and there are only two  coordinates in ( 12)- (15) shown as follows: Mathematical Problems in Engineering 7 Equation (19) shows that if  1 ,  2 , and  3 (coordinates of sensors) and  12  0 and  23  0 (arrival time differences) are known,  ( coordinate of seismic source) can be regarded as a function of V (average wave velocity), denoted as  = (V).Obviously, there is a one-to-one correspondence between  and V since V > 0 in practice; that is to say,  exists and is unique when V is known.
An approximate graph can be plotted by (19).Take the derivative of ( 19) with respect to V: where Equation (20) shows that the sign of /V, namely, the monotonicity of  = (V), depends on the sign of (−)( 1 −  3 ) and then depends on the relative position of seismic source and sensors.There are four possible relative positions of seismic source and sensors as shown in Figure 8.
As B shows in Figure 8, and  is a monotonous decreasing function of V.
As C shows in Figure 8, and  is a monotonous increasing function of V.
According to (19), when V → 0, That is, the  coordinate of seismic source tends to a constant when V → 0, which is determined by coordinates of sensors and arrival time differences.In addition, when the distance between seismic source and sensors tends towards infinity, there is Equation ( 19) also shows that therefore, the supremum of the domain of  = (V) is finite, denoted as V max .The value of V max depends on the smaller value between |( 3 −  2 )/ 23  0 | and |( 2 −  1 )/ 12  0 |; in other words, it depends on sensor interval and arrival time difference of the more distant sensors from seismic source.That is to say, when V → V max , the  coordinate of seismic source tends to the midpoint of the closer sensors from the seismic source, which is  = (  +  +1 )/2.
So, the graph of  = (V) can be approximately plotted as shown in Figure 9.
There are two real cases of  = (V) as shown in Figure 10.
The parameters of (a) and (b) in Figure 10 are shown as follows.
The coordinates of sensors are (−10, 0)-(−200, 0) and the interval of sensors is 1 meter.The average wave velocity is V = 300 m/s.And the coordinate of seismic source is (0, 10).The line nearly parallel to the horizontal axis in (b) is the graph of  = (V) whose sensor coordinates are (−10, 0) and (−11, 0) while the other curve is the graph whose sensor coordinates are (−110, 0) and (−120, 0).
The parameters of (c) and (d) in Figure 10 are shown as follows.
The coordinates of sensors are (10, 0)-(200, 0) and the interval of sensors is 1 meter.The average wave velocity is V = 300 m/s.And the coordinate of seismic source is (0, 10).The line nearly parallel to the horizontal axis in (d) is the  = (V) graph whose sensor coordinates are (10, 0) and (11, 0) while the other curve is the graph whose sensor coordinates are (110, 0) and (120, 0).
There is an inference we can make from Figure 10, which is that the intersection point of all the  = (V) graphs exists and is unique; the reasons are as follows.
Existence.There must be at least one intersection point of the  = (V) graphs; the horizontal coordinate of the intersection is the real average wave velocity V and the vertical coordinate is the real horizontal coordinate of seismic source.Uniqueness.Suppose there is another intersection, named (  , V  ), of two  = (V) graphs.As mentioned above, when V is known,  (the horizontal coordinate of seismic source) exists and is unique; that is to say, when the average wave velocity is V  , all of the  = (V) graphs will uniquely intersect at one point (  , V  ).However, as presented in (b) and (d) in Figure 10, the sensors away from seismic source and the sensors close to seismic source will intersect at only one point, the real horizontal coordinate of seismic source and the real average wave velocity; there is not another intersection named (  , V  ).Therefore, there is only one intersection of  = (V) graphs.
According to the analyses above, the procedure of the new two-dimensional far field source locating method is as follows.
Step 1. Assume ( 1 ,  1 ), ( 2 ,  2 ), and ( 3 ,  3 ) are coordinates of a set of sensors and  12 and  23 are the corresponding arrival time differences; get an equation relating the horizontal coordinate of seismic source and the average wave velocity based on (19), denoted as  = (V).
Step 3. Define | −   | as the objective function; V min which minimizes | −   | can be considered as true value of average wave velocity.
Step 4. Substitute V min and other known quantities into (12) and (15) to get the coordinate of seismic source.

Theoretical Error.
The reason why the new locating method proposed in Section 3.1 can reduce a twodimensional solution space to a one-dimensional solution space is that the hyperbola is replaced by asymptote when the seismic source is far away from sensors, but this replacing will unquestionably cause theoretical error while locating: The normal-form hyperbolic equation is And the normal-form asymptotic equation is  = ±   . ( Δ is a decreasing function of |/| and lim |/|→∞ Δ = 0; Δ is a decreasing function of |/| and lim |/|→∞ Δ = 0. Some examples are shown in Table 1. As Table 1 shows, when |/| > 7, there is |Δ| < 1%, and when |/| > 7, there is |Δ| < 1%.Assume  1 = − 0 ( 0 > 0),  2 =  0 , and  1 =  2 = 0; (9) will change into In comparison with ( 22) and ( 25), the inequation |/| > 7 corresponding to (22) will change into (26) corresponding to (25): And the inequation |/| > 7 will change into Since 27) can be rewritten to || > 7 0 > 7( 0 V/2); therefore, when || > 3.5 ⋅ 2 0 , there is |Δ| < 1%.That is to say, when the distance between the vertical coordinate of seismic source and the vertical coordinate of midpoint of the two sensors reaches 3.5 times larger than sensor interval, the relative error of horizontal coordinate of locating result is less than 1%, and when the distance between the horizontal coordinate of seismic source and the horizontal coordinate of midpoint of the two sensors reaches 3.5 times larger than sensor interval, the relative error of vertical coordinate of locating result is less than 1%.Similarly, when || > 2.5 ⋅ 2 0 , there is |Δ| < 2%, and when || > 2.5 ⋅ 2 0 , there is |Δ| < 2%.
According to the analysis above, the distance range of far field source depends on the demanded locating accuracy.For example, assume   and   are coordinates of seismic source,  1 and  1 are coordinates of the closest sensor, and  2 and  2 are coordinates of the next-closest sensor; if the demanded locating accuracy is less than 1%, then the seismic sources which satisfy |  − ( 1 +  2 )/2| > 3.5| 2 −  1 | and |  − ( 1 +  2 )/2| > 3.5| 2 −  1 | could be considered as far field sources.

Numerical Experiments and Discussion
4.1.Optimization Method.GA (Genetic Algorithm) is a naturally parallel method of optimization, which can be conveniently migrated to parallel environments to improve computing speed.In this paper, we use SGA (Simple Genetic Algorithm) to optimize the objective function, proposed in Section 3.1.The program is developed by a toolbox called GATBX (published by UK's University of Sheffield).The parameters of SGA used in this paper are shown in Table 2.

Validity Experiments.
Assume the coordinates of sensors are (970, 0), (980, 0), (990, 0), and (1000, 0), and the average wave velocity is 3000 m/s.To test the validity of the new method proposed in this paper, we move the seismic source where RE denotes the relative error, IR denotes the inversion result, and TV denotes the true value.
The LEFT (means the  coordinate of seismic source satisfies  ∈ (60, 960)) relative errors of V, , and  are shown in Figure 11, and the RIGHT (means the  coordinate of seismic source satisfies  ∈ (1010, 2010)) errors are shown in Figure 12.
All the comparatively larger relative errors highlighted in Figures 11 and 12 are shown in Table 3.
As shown in Table 3, the comparatively larger relative errors will arise when seismic source comes close to sensors which is consistent with the analysis in Section 3.2.For example, when the largest relative error of LEFT  is −0.1014% and the corresponding  coordinate of seismic source is 10, the distance between the  coordinate of seismic source and the  coordinate of sensors is equal to sensor interval.One more example is that when the largest relative error of LEFT  is −1.72%, the corresponding  coordinate of seismic source is 960, and the  coordinate is 10, it is obvious that the distance between the  coordinate of seismic source and the  coordinate of sensors is equal to sensor interval and the minimum distance between the  coordinate of seismic source and the  coordinate of the midpoint of sensors is slightly greater than sensor interval.According to the numerical experimental results and analyses above, the locating result can be considered unauthentic when the minimum distance between the located result and sensors is less than 2.5 times the sensor interval.

Contrast Experiments.
Assume the coordinates of sensors are (970, 0), (980, 0), (990, 0), and (1000, 0), the coordinates of seismic source are (600, 1000), and the average wave velocity is 3000 m/s; the corresponding theoretical arrival time differences are  21 ≈ 1.1704 ms,  32 ≈ 1.1976 ms, and  43 ≈ 1.2246 ms after numerical simulation.Table 4 shows 20 contrast experiment results based on the coordinates of sensors and the theoretical arrival time differences calculated above, where NEW denotes the method taking | −   | as objective function and OLD denotes the method taking (6) as objective function.The optimization methods of NEW and OLD are both SGA and the parameters of SGA are shown in Table 2.  Comparison chart of convergence speed is shown in Figure 13.
Obviously, Figure 13 shows that the NEW method converges faster than the OLD method and also the NEW method's fluctuation of iterative times is smaller than the OLD method.
The stability of optimal solution is shown in Figure 14, and clearly the optimal solution of NEW method is more stable than OLD method.
According to Table 4, neither optimal solution of the two methods converges to true value, which is caused by the poor timing accuracy.If we substitute arrival time differences with higher accuracy into the NEW method, such as  21 = 1.170400 ms,  32 = 1.197628 ms, and  43 = 1.224583 ms, then the optimal solution of V will be 3002.4. Figure 15 shows errors between calculated results and true values where the calculated results were obtained from NEW and OLD method by 1 dB-200 dB white Gauss noise added arrival time differences (arrival time differences are  21 = 1.170400227006 × 10 −3 s,  32 = 1.197627771071 × 10 −3 s, and  43 = 1.224582829460 × 10 −3 s). Figure 16 shows enlarged view of Figure 15.As shown in Figures 15  and 16, the calculation accuracy of OLD method is better than NEW method when SNR (Signal-to-Noise Ratio) is lower than about 80 dB, and the calculation accuracy of NEW method will be higher than OLD method's when SNR becomes higher than about 85 dB.Obviously, OLD method has better antinoise ability than NEW method, but NEW method will give more accurate results with higher timing precision.

Conclusions and Discussions
This paper proposes a fast and stable method for twodimensional far filed source locating with nonprior velocity and the new method has been proved valid by numerical experimentations.The new method uses asymptote instead of hyperbola to locate seismic source which reduces the dimension of solution space from two to one.Numerical experiments show that the new method works faster and more stable than usual locating method at expense of acceptable locating accuracy.However, the method proposed in this paper is still imperfect.
(1) There is no strict theoretical proof for the uniqueness of graphs' intersection point of  = (V) in Section 3.1.
(2) The relative errors shown in Figures 11 and 12 and Table 3 are smaller than the theoretic errors got from Section 3.2.The maximum relative error of  coordinate is −0.1014% (LEFT) and −0.08024% (RIGHT), respectively; they are both smaller than the theoretic error (1%) deduced in Section 3.2.
(3) NEW method has weaker antinoise ability than OLD method.The three questions mentioned above will need to be solved in future research.

Figure 1 :
Figure 1: Model of two-dimensional far field source locating.

Figure 8 :
Figure 8: Relative positions of seismic source and sensors.

Figure 13 :
Figure 13: Comparison chart of convergence speed.

Table 1 :
Relative error between hyperbola and asymptote.

Table 3 :
Values and coordinates of comparatively larger relative errors.

Table 4 :
Results of contrast experiments.