# Two-Dimensional Far Field Source Locating Method with Nonprior Velocity.

1. Introduction

Microseismic monitoring, which was first studied by Obert [1], has been successfully applied in mining [2-5], oil exploitation [6-10], water conservancy construction [11-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:

[mathematical expression not reproducible] (1)

where x, y, and z are coordinates of seismic source and [t.sub.s] is the initial time of source; [x.sub.i], [y.sub.i], and [z.sub.i] are coordinates of the ith sensor; [v.sub.i] is the average wave velocity measured by the ith sensor; [t.sub.i] is arrival time measured by ith sensor. It is obvious that 4 or more sensors at least are needed to solve for (x, y, z, [t.sub.s]) if [v.sub.i] 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-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]. Geigers 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-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-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.

2. 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.

([x.sub.1], [y.sub.1]), ([x.sub.2], [y.sub.2]), ([x.sub.3], [y.sub.3]), and ([x.sub.4], [y.sub.4]) are coordinates of equidistant placed sensors, (x, y) are coordinates of seismic source, v is the average wave velocity, [t.sub.12] is the arrival time difference between ([x.sub.1], [y.sub.1]) and ([x.sub.2], [y.sub.2]), and [t.sub.23] and [t.sub.34] are similar to [t.sub.12]. The unknowns are (x, y) and v. If v is known, (x, y) can be solved by the following equation set:

[mathematical expression not reproducible]. (2)

If v is unknown, the three unknowns can be solved by the following equation set:

[mathematical expression not reproducible]. (3)

Obviously, v in equation set (3) can be eliminated and equation set (3) will change to the following equation set:

[mathematical expression not reproducible]. (4)

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):

[mathematical expression not reproducible], (5)

[mathematical expression not reproducible]. (6)

According to Li's method, the objective function of equation set (3) can be defined as

[([D.sub.12] - [D.sub.23]).sup.2] + [([D.sub.23] - [D.sub.34]).sup.2], (7)

where

[mathematical expression not reproducible]. (8)

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.

Assume the coordinates of sensors (units of coordinates are in meters unless particularly stated in this paper) are (970, 0), (980, 0), (990, 0), and (1000, 0); the average wave velocity is 3000 m/s and the coordinate of seismic source is (500, 500). The graph of objective function (6) when x [member of] (100, 1000) and y [member of] (100, 1000) is shown in Figure 2.

Enlarged view of Figure 2 near (500, 500) is shown in Figure 3.

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 x coordinate is better than y 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.

3. 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 coordinates of seismic source (x, y) satisfy hyperbolic equation (9), shown as follows:

[mathematical expression not reproducible], (9)

where ([x.sub.1], [y.sub.1]) and ([x.sub.2], [y.sub.2]) are coordinates of sensors, v is the average wave velocity, [t.sub.0] is the sampling period of A/D, and N[t.sub.0] is the arrival time difference between ([x.sub.1], [y.sub.1]) and ([x.sub.2], [y.sub.2]); for convenience sake, assume N[t.sub.0] > 0.

The asymptotic equation of (9) is expressed in

[mathematical expression not reproducible]. (10)

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

[mathematical expression not reproducible]. (11)

Define [N.sub.12][t.sub.0] = [absolute value of [t.sub.12]] and [N.sub.23][t.sub.0] = [absolute value of [t.sub.23]]; the solutions of (11) are shown as follows:

[mathematical expression not reproducible], (12)

[mathematical expression not reproducible], (13)

[mathematical expression not reproducible], (14)

[mathematical expression not reproducible], (15)

where

[mathematical expression not reproducible]. (16)

All the y coordinates are equal when sensors are placed in a line as shown in Figure 1. Therefore, [A.sub.12] = [A.sub.23], and there are only two x coordinates in (12)-(15) shown as follows:

x = [[B.sub.23][C.sub.23] - [B.sub.12][C.sub.12]/[B.sub.23] - [B.sub.12]], (17)

x = [[B.sub.23][C.sub.23] + [B.sub.12][C.sub.12]/[B.sub.23] + [B.sub.12]]. (18)

Equation (17) is the x coordinate of far field source. Substitute [B.sub.12], [C.sub.12], [B.sub.23], and [C.sub.23] into (17) and get

[mathematical expression not reproducible]. (19)

Equation (19) shows that if [x.sub.1], [x.sub.2], and [x.sub.3] (coordinates of sensors) and [N.sub.12][t.sub.0] and [N.sub.23][t.sub.0] (arrival time differences) are known, x (x coordinate of seismic source) can be regarded as a function of v (average wave velocity), denoted as x = f(v). Obviously, there is a one-to-one correspondence between x and v since v > 0 in practice; that is to say, x 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:

[dx/dv] = [(v/2) (B - A)(B + A) ([x.sub.1] + [x.sub.2] - [x.sub.2] - [x.sub.3])/ AB [(B - A).sup.2], (20)

where A = [square root of [(([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.0]).sup.2] - [v.sub.2]] and B = [square root of [(([x.sub.3] - [x.sub.2])/[N.sub.23][t.sub.0]).sup.2] - [v.sub.2].

Equation (20) shows that the sign of dx/dv, namely, the monotonicity of x = f(v), depends on the sign of (B - A)([x.sub.1] -[x.sub.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 (1) shows in Figure 8, [x.sub.1] - [x.sub.3] < 0, [(([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.0]).sup.2] [(([x.sub.3] - [x.sub.2])/[N.sub.23][t.sub.0]).sup.2]; therefore, dx/dv = (B - A)([x.sub.1] - [x.sub.3]) > 0, and x is a monotonous increasing function of v.

As (2) shows in Figure 8, [x.sub.1] - [x.sub.3] < 0, [(([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.0]).sup.2] < [(([x.sub.3] - [x.sub.2])/[N.sub.23][t.sub.0]).sup.2]; therefore, dx/dv = (B - A)([x.sub.1] - [x.sub.3]) < 0, and x is a monotonous decreasing function of v.

As (3) shows in Figure 8, [x.sub.1] - [x.sub.3] > 0, [(([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.0]).sup.2] > [(([x.sub.3] - [x.sub.2])/[N.sub.23][t.sub.0]).sup.2]; therefore, dx/dv = (B - A)([x.sub.1] - [x.sub.3]) > 0, and x is a monotonous increasing function of v.

As (4) shows in Figure 8, [x.sub.1] - [x.sub.3] > 0, [(([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.0]).sup.2] < [(([x.sub.3] - [x.sub.2])/[N.sub.23][t.sub.0]).sup.2]; therefore, dx/dv = (B - A)([x.sub.1] - [x.sub.3]) < 0, and x is a monotonous decreasing function of v.

According to (19), when v [right arrow] 0,

[mathematical expression not reproducible]. (21)

That is, the x coordinate of seismic source tends to a constant when v [right arrow] 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 [absolute value of ([x.sub.3] - [x.sub.2])/[N.sub.23][t.sub.0]] [approximately equal to] [absolute value of ([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.0]] [approximately equal to] v, and therefore x [right arrow] +[infinity] or x [right arrow] -[infinity].

Equation (19) also shows that [(([x.sub.3] - [x.sub.2])/[N.sub.23][t.sub.0]).sup.2] - [v.sup.2] [greater than or equal to] 0 and [(([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.0]).sup.2] - [v.sup.2] [greater than or equal to] 0; therefore, the supremum of the domain of x = f(v) is finite, denoted as [v.sub.max]. The value of [v.sub.max] depends on the smaller value between [absolute value of ([x.sub.3] -[x.sub.2])/[N.sub.23][t.sub.0]] and [absolute value of ([x.sub.2] - [x.sub.1])/[N.sub.12][t.sub.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 [right arrow] [v.sub.max], the x coordinate of seismic source tends to the midpoint ofthe closer sensors from the seismic source, which is x = ([x.sub.m] + [x.sub.m+1])/2.

So, the graph of x = f(v) can be approximately plotted as shown in Figure 9.

There are two real cases of x = f(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 x = f(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 x = f(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 x = f(v) graphs exists and is unique; the reasons are as follows.

Existence. There must be at least one intersection point of the x = f(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 (x', V'), of two x = f(v) graphs. As mentioned above, when V is known, x (the horizontal coordinate of seismic source) exists and is unique; that is to say, when the average wave velocity is v', all of the x = f(v) graphs will uniquely intersect at one point (x', 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 (x', v'). Therefore, there is only one intersection of x = f(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 ([x.sub.1], [y.sub.1]), ([x.sub.2], [y.sub.2]), and ([x.sub.3], [y.sub.3]) are coordinates of a set of sensors and [t.sub.12] and [t.sub.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 x = f(v).

Step 2. Assume ([x.sub.2], [y.sub.2]), ([x.sub.3], [y.sub.3]), and ([x.sub.4], [y.sub.4]) are coordinates of a set of sensors and [t.sub.23] and [t.sub.34] are the corresponding arrival time differences; another equation will be got based on (19), denoted as x' = f(v).

Step 3. Define [absolute value of x - x] as the objective function; [v.sub.min] which minimizes [absolute value of x - x] can be considered as true value of average wave velocity.

Step 4. Substitute [v.sub.min] and other known quantities into (12) and (15) to get the coordinate of seismic source.

3.2. Theoretical Error. The reason why the new locating method proposed in Section 3.1 can reduce a two-dimensional 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 [[[x.sup.2]/[a.sup.2]] - [[y.sup.2]/[b.sup.2]] = 1. (22)

And the normal-form asymptotic equation is y

= [+ or -][b/a] x. (23)

With (22), the horizontal coordinate of a point on hyperbola is x = [+ or -]a[square root of [y.sup.2]/[b.sup.2] + 1], and the vertical coordinate is y = [+ or -]b [square root of [x.sup.2]/[a.sup.2] - 1]. With (23), the horizontal coordinate of a point on asymptote is x = [+ or -]a(y/b), and the vertical coordinate is y = [+ or -]b(x/a):

Define the relative error of horizontal coordinate as [DELTA]x

[mathematical expression not reproducible]. (24)

Define the relative error of vertical coordinate as [DELTA]Ay

[mathematical expression not reproducible].

[DELTA]x is a decreasing function of [absolute value of y/b] and [lim.sub.[absolute value of y/b][right arrow][infinity]] [DELTA]x = 0; [DELTA]y is a decreasing function of [absolute value of x/a] and [lim.sub.[absolute value of x/a][right arrow][infinity]] [DELTA]y = 0. Some examples are shown in Table 1.

As Table 1 shows, when [absolute value of y/b] > 7, there is [absolute value of [DELTA]x] < 1%, and when [absolute value of x/a] > 7, there is [absolute value of [DELTA]y] < 1%. Assume [x.sub.1] = -[x.sub.0] ([x.sub.0] > 0), [x.sub.2] = [x.sub.0], and [y.sub.1] = [y.sub.2] = 0; (9) will change into

[mathematical expression not reproducible]. (25)

In comparison with (22) and (25), the inequation [absolute value of y/b] > 7 corresponding to (22) will change into (26) corresponding to (25):

[mathematical expression not reproducible]. (26)

And the inequation [absolute value of x/a] > 7 will change into

[[absolute value of x]/N[t.sub.0]v/2] > 7. (27)

Since [square root of [x.sup.2.sub.0] - [N.sup.2][t.sup.2.sub.0][v.sup.2.sub.0]/4] < [square root of [x.sup.2.sub.0]], (26) can be rewritten to [absolute value of y] > 7[square root of [x.sup.2.sub.0]] > 7[square root of [x.sup.2.sub.0] - [N.sup.2][t.sup.2.sub.0][v.sup.2]/4]; therefore, when [absolute value of y] > 3.5 x 2[x.sub.0], there is [absolute value of [DELTA]x] < 1%. Since N[t.sub.0]v < 2[x.sub.0], (27) can be rewritten to [absolute value of x] > 7[x.sub.0] > 7(N[t.sub.0]v/2); therefore, when [absolute value of x] > 3.5 x 2[x.sub.0], there is [absolute value of [DELTA]y] < 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 [absolute value of y] > 2.5 x 2[x.sub.0], there is [absolute value of [DELTA]x] < 2%, and when [absolute value of x] > 2.5 x 2[x.sub.0], there is [absolute value of [DELTA]y] < 2%.

According to the analysis above, the distance range of far field source depends on the demanded locating accuracy. For example, assume [x.sub.s] and [y.sub.s] are coordinates of seismic source, [x.sub.n1] and [y.sub.n1] are coordinates of the closest sensor, and [x.sub.n2] and [y.sub.n2] are coordinates of the next-closest sensor; if the demanded locating accuracy is less than 1%, then the seismic sources which satisfy [absolute value of [x.sub.s] - ([x.sub.n1] + [x.sub.n2])/2] > 3.5[absolute value of [x.sub.n2] - [x.sub.n1]] and [absolute value of [y.sub.s] - ([y.sub.n1] + [y.sub.n2])/2] > 3.5[absolute value of [x.sub.n2] - [x.sub.n1]] could be considered as far field sources.

4. 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.

4.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 in x, y plane at interval of 100 meters and calculate the arrival time (in seconds with 16 significant digits) of each point, where x [member of] (60, 960) or x [member of] (1010, 2010) and y [member of] (10, 1010). Define relative error formula of v, x, and y as

RE = [IR - TV/TV] x 100%, (28)

where RE denotes the relative error, IR denotes the inversion result, and TV denotes the true value.

The LEFT (means the x coordinate of seismic source satisfies x [member of] (60, 960)) relative errors of v, x, and y are shown in Figure 11, and the RIGHT (means the x coordinate of seismic source satisfies x [member of] (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 X is -0.1014% and the corresponding Y coordinate of seismic source is 10, the distance between the Y coordinate of seismic source and the Y coordinate of sensors is equal to sensor interval. One more example is that when the largest relative error of LEFT Y is -1.72%, the corresponding X coordinate of seismic source is 960, and the Y coordinate is 10, it is obvious that the distance between the Y coordinate of seismic source and the Y coordinate of sensors is equal to sensor interval and the minimum distance between the X coordinate of seismic source and the X 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.

4.3. 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 [t.sub.21] [approximately equal to] 1.1704 ms, [t.sub.32] [approximately equal to] 1-1976 ms, and [t.sub.43] [approximately equal to] 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 [absolute value of x - x'] 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 [t.sub.21] = 1.170400 ms, [t.sub.32] = 1.197628 ms, and [t.sub.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 [t.sub.21] = 1.170400227006 x [10.sup.-3] s, [t.sub.32] = 1.197627771071 x [10.sup.-3] s, and [t.sub.43] = 1.224582829460 x [10.sup.-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.

5. Conclusions and Discussions

This paper proposes a fast and stable method for two-dimensional 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 x = f(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 x 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.

http://dx.doi.org/10.1155/2016/8540303

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors thank the financial support of the National Key Basic Research Program of China, 973 Program, Grant no. 2014CB046300.

References

[1] L. Obert, "The microseismic method: discovery and early history," in Proceedings of the 1st Conference on Acoustic Emission/ Microseismic Activity in Geologic Structures and Materials, pp. 11-12, June 1977.

[2] H. Wang and M. Ge, "Acoustic emission/microseismic source location analysis for a limestone mine exhibiting high horizontal stresses," International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 5, pp. 720-728, 2008.

[3] M. Sikora, "Induction and pruning of classification rules for prediction of microseismic hazards in coal mines," Expert Systems with Applications, vol. 38, no. 6, pp. 6748-6758, 2011.

[4] M. Cai, P. K. Kaiser, and C. D. Martin, "Quantification of rock mass damage in underground excavations from microseismic event monitoring," International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 8, pp. 1135-1145, 2001.

[5] M. Ge, "Efficient mine microseismic monitoring," International Journal of Coal Geology, vol. 64, no. 1-2, pp. 44-56, 2005.

[6] S. Sasaki, "Characteristics of microseismic events induced during hydraulic fracturing experiments at the Hijiori hot dry rock geothermal energy site, Yamagata, Japan," Tectonophysics, vol. 289, no. 1-3, pp. 171-188, 1998.

[7] H. Zhang, S. Sarkar, M. N. Toksoz, H. S. Kuleli, and F. Al-Kindy, "Passive seismic tomography using induced seismicity at a petroleum field in Oman," Geophysics, vol. 74, no. 6, pp. WCB57-WCB69, 2009.

[8] O. H. Al-Harrasi and A. Al-Anboori, "Seismic anisotropy in a hydrocarbon field estimated from microseismic data," Geophysical Prospecting, vol. 59, no. 2, pp. 227-243, 2011.

[9] G. A. Jones, D. Raymer, K. Chambers, and J.-M. Kendall, "Improved microseismic event location by inclusion of a priori dip particle motion: a case study from Ekofisk," Geophysical Prospecting, vol. 58, no. 5, pp. 727-737, 2010.

[10] S. I. Kaka, "Passive microseismic experiments at King Fahd University of Petroleum and Minerals in Saudi Arabia," Seismological Research Letters, vol. 83, no. 4, pp. 680-685, 2012.

[11] N. W. Xu, C. A. Tang, L. C. Li et al., "Microseismic monitoring and stability analysis of the left bank slope in Jinping first stage hydropower station in southwestern China," International Journal of Rock Mechanics and Mining Sciences, vol. 48, no. 6, pp. 950-963, 2011.

[12] N. W. Xu, F. Dai, Z. Z. Liang, Z. Zhou, C. Sha, and C. A. Tang, "The dynamic evaluation of rock slope stability considering the effects of microseismic damage," Rock Mechanics and Rock Engineering, vol. 47, no. 2, pp. 621-642, 2014.

[13] N.-W. Xu, C.-A. Tang, H. Li et al., "Excavation-induced microseismicity: microseismic monitoring and numerical simulation," Journal of Zhejiang University SCIENCE A, vol. 13, no. 6, pp. 445-460, 2012.

[14] V. Inglada, "Die berechnung der herdkoordinated eines nahbebens aus den dintrittszeiten der in einingen benachbarten stationen aufgezeichneten P-oder P-wellen," Gerlands Beitrage zur Geophysik, vol. 19, pp. 73-98, 1928.

[15] M. Ge, "Analysis of source location algorithms, part I: overview and non-iterative methods," Journal of Acoustic Emission, vol. 21, pp. 14-28, 2003.

[16] F. Leighton and W. Blake, "Rock noise source location techniques," Tech. Rep. RI 7432, US Bureau of Mines, Washington, DC, USA, 1970.

[17] F. Leighton and W. I. Duvall, A Least Squares Method for Improving Rock Noise Source Location Techniques, CDC, 1972.

[18] L. Geiger, "Herdbestimmung bei Erdbeben aus den Ankunftszeiten," Konigliche Gesellschaft der Wissenschaften zu Gottingen, vol. 4, pp. 331-349, 1910.

[19] L. Geiger, "Probability method for the determination of earthquake epicenters from the arrival time only (translated from Geiger's 1910 German article)," Bulletin of St. Louis University, vol. 8, pp. 56-71, 1912.

[20] J. A. Nelder and R. Mead, "A simplex method for function minimization," The Computer Journal, vol. 7, no. 4, pp. 308-313, 1965.

[21] A. F. Prugger and D. J. Gendzwill, "Microearthquake location: a nonlinear approach that makes use of a simplex stepping procedure," Bulletin of the Seismological Society of America, vol. 78, no. 2, pp. 799-815, 1988.

[22] J. H. Holland, Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence, MIT Press, Boston, Mass, USA, 1992.

[23] K. R. Anderson, "Robust earthquake location using Mestimates," Physics of the Earth and Planetary Interiors, vol. 30, no. 2-3, pp. 119-130, 1982.

[24] W. H. K. Lee and S. W. Stewart, Principles and Applications of Microearthquake Networks, Academic Press, 1981.

[25] M. Ge, "Analysis of source location algorithms part II: iterative methods," Journal of Acoustic Emission, vol. 21, no. 1, pp. 29-51, 2003.

[26] J. Li, S.-C. Wu, Y.-T. Gao, L.-J. Li, and Y. Zhou, "An improved multidirectional velocity model for micro-seismic monitoring in rock engineering," Journal of Central South University, vol. 22, no. 6, pp. 2348-2358, 2015.

[27] N. Li, E. Wang, M. Ge, and Z. Sun, "A nonlinear microseismic source location method based on Simplex method and its residual analysis," Arabian Journal of Geosciences, vol. 7, no. 11, pp. 4477-4486, 2013.

[28] J.-P. Liu, Y.-H. Li, and S.-D. Xu, "Estimation of cracking and damage mechanisms of rock specimens with precut holes by moment tensor analysis of acoustic emission," International Journal of Fracture, vol. 188, no. 1, pp. 1-8, 2014.

[29] B. L. N. Kennett and M. S. Sambridge, "Earthquake location-genetic algorithms for teleseisms," Physics of the Earth and Planetary Interiors, vol. 75, no. 1-3, pp. 103-110, 1992.

[30] M. Sambridge and K. Gallagher, "Earthquake hypocenter location using genetic algorithms," Bulletin of the Seismological Society of America, vol. 83, no. 5, pp. 1467-1491, 1993.

[31] J. Sileny, "Earthquake source parameters and their confidence regions by a genetic algorithm with a 'memory'," Geophysical Journal International, vol. 134, no. 1, pp. 228-242, 1998.

[32] W. Kim, I.-K. Hahm, S. J. Ahn, and D. H. Lim, "Determining hypocentral parameters for local earthquakes in 1-D using a genetic algorithm," Geophysical Journal International, vol. 166, no. 2, pp. 590-600, 2006.

[33] M. Zaharia, M. Chowdhury, T. Das et al., "Resilient distributed datasets: a fault-tolerant abstraction for in-memory cluster computing," in Proceedings of the 9th USENIX Conference on Networked Systems Design and Implementation, p. 2, USENIX Association, San Jose, Calif, USA, April 2012.

[34] W. Yun and Y. Xi, "The analysis of global convergence and computational efficiency for genetic algorithm," Control Theory and Applications, vol. 13, no. 4, pp. 455-460, 1996 (Chinese).

[35] L. He, K. Wang, G. Li et al., "The discussion about the paper 'the analysis of global convergence and computational efficiency for genetic algorithm'," Control Theory & Applications, vol. 18, no. 1, pp. 142-145, 2001 (Chinese).

[36] L. Dong, X. Li, L. Tang, and F. Gong, "Mathematical functions and parameters for microseismic source location without pre-measuring speed," Chinese Journal of Rock Mechanics and Engineering, vol. 30, no. 10, pp. 2057-2067, 2011 (Chinese).

[37] L. Dong and X. Li, "A microseismic/acoustic emission source location method using arrival times of PS waves for unknown velocity system," International Journal of Distributed Sensor Networks, vol. 2013, Article ID 307489, 8 pages, 2013.

[38] J. Li, Y. T. Gao, Y. Xie, Y. Zhou, and K. Yang, "Improvement of microseism locating based on simplex method without velocity measuring," Chinese Journal of Rock Mechanics and Engineering, vol. 33, no. 7, pp. 1336-1346, 2015 (Chinese).

Qing Chen, (1,2,3) Xiaowen Liu, (1,3) Juanjuan Li, (1,3) Manyi Wang, (4) and Peng Liu (1,3)

(1) IOT Perception Mine Research Center, China University of Mining and Technology, Xuzhou 221008, China

(2) School of Mathematics & Physics Science, Xuzhou Institute of Technology, Xuzhou 221008, China

(3) School of Information and Electrical Engineering, China University of Mining and Technology, Xuzhou 221008, China

(4) School of Mechanical Engineering, Nanjing University of Science & Technology, Nanjing 221000, China

Correspondence should be addressed to Xiaowen Liu; chenqxzit@gmail.com

Received 22 December 2015; Accepted 30 March 2016

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

Caption: Figure 2: 45-degree view of (6) when seismic source locates at (500, 500).

Caption: Figure 3: 45-degree/XZ/7Z enlarged view of (6).

Caption: Figure 4: 45-degree view of (7) when seismic source locates at (500, 500).

Caption: Figure 5: 45-degree/XZ/YZ enlarged view of (7) near (500, 500).

Caption: Figure 6: 45-degree/'XZ/YZ view of (6) when seismic source locates at (1000, 500).

Caption: Figure 7: 45-degree/'XZ/YZ enlarged view of (7) when seismic source locates at (1000, 500).

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

Caption: Figure 9: Approximate graph of x = f(v).

Caption: Figure 10: Real graphs of x = f(v).

Caption: Figure 11: Relative error of LEFT V/X/Y.

Caption: Figure 12: Relative Error of RIGHT V/X/Y.

Caption: Figure 13: Comparison chart of convergence speed.

Caption: Figure 14: Stability of optimal solution.

Caption: Figure 15: Calculated results by 1-200 dB white Gauss noise added arrival time differences.

Caption: Figure 16: Calculated results by 80-200 dB white Gauss noise added arrival time differences.
```Table 1: Relative error between hyperbola and asymptote.

[absolute     [DELTA]x               [absolute       [DELTA]y
value of                             value of x/a]
y/b]

2              -0.105572809000084    2               0.154700538379252
3              -0.051316701949486    3               0.060660171779821
4              -0.029857499854668    4               0.032795558988644
5              -0.019419324309080    5               0.020620726159658
6              -0.013606076167856    6               0.014185105674220
7              -0.010050506338834    7               0.010362971081845
8              -0.007722123286332    8               0.007905261357939

Table 2: Parameters of SGA.

Number of     Maximum       Variable
individuals   number of     number of
generations   binary
digits

100           100           25

Gap of        Crossover     Mutation
generations   probability   probability

0.95          0.7           0.1

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

Error     Relative     True value of      True value of
name      error (%)   X coordinate (m)   Y coordinate (m)

LEFT V     0.02626          960                 10
LEFT V    -0.08622          960                110
LEFT X     0.08442          960                 10
LEFT X     0.05909           60                110
LEFT X     -0.1014           60                 10
LEFT Y     0.1885           960                110
LEFT Y      -1.72           960                 10
RIGHT V    0.02626          1010                10
RIGHT V   -0.08622          1010               110
RIGHT X   -0.08024          1010                10
RIGHT Y    0.1885           1010               110
RIGHT Y     -1.72           1010                10

Table 4: Results of contrast experiments.

Experiment    Iteration      Optimum v       Error of
number        number of        of OLD       OLD method
OLD method       method       between V
and true
value

1                100         2569.7434      -430.2566
2                100         2665.3755      -334.6245
3                 90         2306.9505      -693.0495
4                100         2237.0807      -762.9193
5                 71         2205.2818      -794.7182
6                 99         2660.2377      -339.7623
7                100          2476.904       -523.096
8                100         2204.2403      -795.7597
9                 67         2474.3871      -525.6129
10                35         2165.6936      -834.3064
11               100         2573.8501      -426.1499
12                75         2475.6581      -524.3419
13                80         1756.5821      -1243.4179
14                96         1759.8416      -1240.1584
15                30          2276.137       -723.863
16                87         1846.0683      -1153.9317
17                41         2781.9136      -218.0864
18                60         1756.4908      -1243.5092
19                64         3280.5975       280.5975
20                99         2343.6613      -656.3387

Experiment    Iteration      Optimum v       Error of
number        number of        of NEW       NEW method
NEW method       method       between v
and true
value

1                 17          2614.70        -385.30
2                 13          2614.70        -385.30
3                 18          2614.70        -385.30
4                 16          2614.70        -385.30
5                 13          2614.70        -385.30
6                 6           2614.70        -385.30
7                 15          2614.70        -385.30
8                 14          2614.70        -385.30
9                 17          2614.70        -385.30
10                16          2614.70        -385.30
11                25          2614.70        -385.30
12                19          2614.70        -385.30
13                31          2614.70        -385.30
14                14          2614.70        -385.30
15                19          2614.70        -385.30
16                14          2614.70        -385.30
17                13          2614.70        -385.30
18                8           2614.70        -385.30
19                13          2614.70        -385.30
20                25          2614.70        -385.30
```