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

Authors: Qing Chen; Xiaowen Liu; Juanjuan Li; Manyi Wang; Peng Liu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Relative position of seismic source and sensors has great influence on locating accuracy, particularly in far field conditions, and the accuracy will decrease seriously due to limited calculation precision and prior velocity error. In order to improve the locating accuracy of far field sources by isometric placed sensors in a straight line, a new locating method with nonprior velocity is proposed. After exhaustive research, this paper states that the hyperbola which is used for locating will be very close to its asymptote when seismic source locates in far field of sensors; therefore, the locating problem with prior velocity is equivalent to solving linear equations and the problem with nonprior velocity is equivalent to a nonlinear optimization problem with respect to the unknown velocity. And then, this paper proposed a new locating method based on a one-variable objective function with respect to the unknown velocity. Numerical experiments show that the proposed method has faster convergence speed, higher accuracy, and better stability.


---

## Body

## 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:(1)
x
-
x
i
2
+
y
-
y
i
2
+
z
-
z
i
2
=
v
i
t
i
-
t
s
,where x, y, and z are coordinates of seismic source and t
s is the initial time of source; x
i, y
i, and z
i are coordinates of the ith sensor; v
i is the average wave velocity measured by the ith sensor; t
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
s
) if v
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]. 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–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 Figure1.Figure 1
Model of two-dimensional far field source locating.The parameters of Figure1 are shown as follows.(
x
1
,
y
1
), (
x
2
,
y
2
), (
x
3
,
y
3
), and (
x
4
,
y
4
) are coordinates of equidistant placed sensors, (
x
,
y
) are coordinates of seismic source, v is the average wave velocity, t
12 is the arrival time difference between (
x
1
,
y
1
) and (
x
2
,
y
2
), and t
23 and t
34 are similar to t
12. The unknowns are (
x
,
y
) and v. If v is known, (
x
,
y
) can be solved by the following equation set:(2)
x
1
-
x
2
+
y
1
-
y
2
v
-
x
2
-
x
2
+
y
2
-
y
2
v
=
t
12
,
x
2
-
x
2
+
y
2
-
y
2
v
-
x
3
-
x
2
+
y
3
-
y
2
v
=
t
23
.Ifv is unknown, the three unknowns can be solved by the following equation set:(3)
x
1
-
x
2
+
y
1
-
y
2
v
-
x
2
-
x
2
+
y
2
-
y
2
v
=
t
12
,
x
2
-
x
2
+
y
2
-
y
2
v
-
x
3
-
x
2
+
y
3
-
y
2
v
=
t
23
,
x
3
-
x
2
+
y
3
-
y
2
v
-
x
4
-
x
2
+
y
4
-
y
2
v
=
t
34
.Obviously,v in equation set (3) can be eliminated and equation set (3) will change to the following equation set:(4)
x
1
-
x
2
+
y
1
-
y
2
-
x
2
-
x
2
+
y
2
-
y
2
x
2
-
x
2
+
y
2
-
y
2
-
x
3
-
x
2
+
y
3
-
y
2
=
t
12
t
23
,
x
1
-
x
2
+
y
1
-
y
2
-
x
2
-
x
2
+
y
2
-
y
2
x
3
-
x
2
+
y
3
-
y
2
-
x
4
-
x
2
+
y
4
-
y
2
=
t
12
t
34
.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):(5)
x
1
-
x
2
+
y
1
-
y
2
-
x
4
-
x
2
+
y
4
-
y
2
v
-
t
12
-
t
23
-
t
34
,
(6)
x
1
-
x
2
+
y
1
-
y
2
-
x
2
-
x
2
+
y
2
-
y
2
x
2
-
x
2
+
y
2
-
y
2
-
x
3
-
x
2
+
y
3
-
y
2
+
x
1
-
x
2
+
y
1
-
y
2
-
x
2
-
x
2
+
y
2
-
y
2
x
3
-
x
2
+
y
3
-
y
2
-
x
4
-
x
2
+
y
4
-
y
2
-
t
12
t
23
-
t
12
t
34
.According to Li’s method, the objective function of equation set (3) can be defined as(7)
D
12
-
D
23
2
+
D
23
-
D
34
2
,where(8)
D
12
=
x
1
-
x
2
+
y
1
-
y
2
-
x
2
-
x
2
+
y
2
-
y
2
t
12
,
D
23
=
x
2
-
x
2
+
y
2
-
y
2
-
x
3
-
x
2
+
y
3
-
y
2
t
23
,
D
34
=
x
3
-
x
2
+
y
3
-
y
2
-
x
4
-
x
2
+
y
4
-
y
2
t
34
.If the seismic source is far away from sensors, as shown in Figure1, 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
∈
(
100,1000
) and y
∈
(
100,1000
) is shown in Figure 2.Figure 2
45-degree view of (6) when seismic source locates at (
500,500
).Enlarged view of Figure2 near (
500,500
) is shown in Figure 3.Figure 3
45-degree/X
Z
/
Y
Z enlarged view of (6).
(a)
45-degree enlarged view of (6) near (
500,500
)
(b)
X
Z enlarged view of (6) near (
500,500
)
(c)
Y
Z enlarged view of (6)The graph of objective function (7) under the same conditions is shown in Figure 4.Figure 4
45-degree view of (7) when seismic source locates at (
500,500
).Enlarged view of Figure4 near (
500,500) is shown in Figure 5.Figure 5
45-degree/X
Z
/
Y
Z enlarged view of (7) near (
500,500
).
(a)
45-degree enlarged view of (7) near (
500,500
)
(b)
X
Z enlarged view of (7) near (
500,500
)
(c)
Y
Z enlarged view of (7)As shown in Figures2–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.Figure 6
45-degree/X
Z
/
Y
Z view of (6) when seismic source locates at (
1000,500
).
(a)
45-degree view of (6) when seismic source locates at (
1000,500
)
(b)
X
Z view of (6) when seismic source locates at (
1000,500
)
(c)
Y
Z view of (6) when seismic source locates at (
1000,500
)Figure 7
45-degree/X
Z
/
Y
Z enlarged view of (7) when seismic source locates at (
1000,500
).
(a)
45-degree view of (7) when seismic source locates at (
1000,500
)
(b)
X
Z view of (7) when seismic source locates at (
1000,500
)
(c)
Y
Z view of (7) when seismic source locates at (
1000,500
)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:(9)
x
-
x
1
+
x
2
/
2
2
N
2
t
0
2
v
2
/
4
-
y
-
y
1
+
y
2
/
2
2
x
2
-
x
1
/
2
2
-
N
2
t
0
2
v
2
/
4
=
1
,where (
x
1
,
y
1
) and (
x
2
,
y
2
) are coordinates of sensors, v is the average wave velocity, t
0 is the sampling period of A/D, and N
t
0 is the arrival time difference between (
x
1
,
y
1
) and (
x
2
,
y
2
); for convenience sake, assume N
t
0
>
0.The asymptotic equation of (9) is expressed in (10)
y
=
y
1
+
y
2
2
±
x
2
-
x
1
2
-
N
2
t
0
2
v
2
N
t
0
v
x
-
x
1
+
x
2
2
.The asymptote will be very close to the hyperbola when the seismic source locates far from sensors as shown in Figure1. 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 (11)
y
=
y
1
+
y
2
2
±
x
2
-
x
1
2
-
N
12
2
t
0
2
v
2
N
12
t
0
v
x
-
x
1
+
x
2
2
,
y
=
y
2
+
y
3
2
±
x
3
-
x
2
2
-
N
23
2
t
0
2
v
2
N
23
t
0
v
x
-
x
2
+
x
3
2
.DefineN
12
t
0
=
t
12 and N
23
t
0
=
t
23; the solutions of (11) are shown as follows:(12)
x
=
A
12
-
A
23
+
B
23
C
23
-
B
12
C
12
B
23
-
B
12
,
y
=
A
12
B
23
-
B
12
A
23
+
B
12
B
23
C
23
-
B
12
B
23
C
12
B
23
-
B
12
,
(13)
x
=
A
23
-
A
12
+
B
23
C
23
+
B
12
C
12
B
23
+
B
12
,
y
=
A
12
B
23
+
B
12
A
23
+
B
12
B
23
C
23
-
B
12
B
23
C
12
B
23
+
B
12
,
(14)
x
=
A
12
-
A
23
+
B
23
C
23
+
B
12
C
12
B
23
+
B
12
,
y
=
A
12
B
23
+
B
12
A
23
-
B
12
B
23
C
23
+
B
12
B
23
C
12
B
23
+
B
12
,
(15)
x
=
A
12
-
A
23
-
B
23
C
23
+
B
12
C
12
B
12
-
B
23
,
y
=
A
12
B
23
-
B
12
A
23
-
B
12
B
23
C
23
+
B
12
B
23
C
12
B
23
-
B
12
,where(16)
A
12
=
y
1
+
y
2
2
,
B
12
=
x
2
-
x
1
2
-
N
12
2
t
0
2
v
2
N
12
t
0
v
,
C
12
=
x
1
+
x
2
2
,
A
23
=
y
2
+
y
3
2
,
B
23
=
x
3
-
x
2
2
-
N
23
2
t
0
2
v
2
N
23
t
0
v
,
C
23
=
x
2
+
x
3
2
.All they coordinates are equal when sensors are placed in a line as shown in Figure 1. Therefore, A
12
=
A
23, and there are only two x coordinates in (12)–(15) shown as follows:(17)
x
=
B
23
C
23
-
B
12
C
12
B
23
-
B
12
,
(18)
x
=
B
23
C
23
+
B
12
C
12
B
23
+
B
12
.Equation (17) is the x coordinate of far field source. Substitute B
12, C
12, B
23, and C
23 into (17) and get (19)
x
=
x
3
-
x
2
/
N
23
t
0
2
-
v
2
x
2
+
x
3
/
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
x
1
+
x
2
/
2
x
3
-
x
2
/
N
23
t
0
2
-
v
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
.Equation (19) shows that if x
1, x
2, and x
3 (coordinates of sensors) and N
12
t
0 and N
23
t
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:(20)
d
x
d
v
=
v
/
2
B
-
A
B
+
A
x
1
+
x
2
-
x
2
-
x
3
A
B
B
-
A
2
,where A
=
(
(
x
2
-
x
1
)
/
N
12
t
0
)
2
-
v
2 and B
=
(
(
x
3
-
x
2
)
/
N
23
t
0
)
2
-
v
2.Equation (20) shows that the sign of d
x
/
d
v, namely, the monotonicity of x
=
f
(
v
), depends on the sign of (
B
-
A
)
(
x
1
-
x
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 Figure8.Figure 8
Relative positions of seismic source and sensors.As ① shows in Figure8, x
1
-
x
3
<
0, (
x
2
-
x
1
)
/
N
12
t
0
2
>
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
>
0, and x is a monotonous increasing function of v.As ② shows in Figure8, x
1
-
x
3
<
0, (
x
2
-
x
1
)
/
N
12
t
0
2
<
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
<
0, and x is a monotonous decreasing function of v.As ③ shows in Figure8, x
1
-
x
3
>
0, (
x
2
-
x
1
)
/
N
12
t
0
2
>
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
>
0, and x is a monotonous increasing function of v.As ④ shows in Figure8, x
1
-
x
3
>
0, (
x
2
-
x
1
)
/
N
12
t
0
2
<
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
<
0, and x is a monotonous decreasing function of v.According to (19), when v
→
0,(21)
lim
v
→
0
⁡
x
3
-
x
2
/
N
23
t
0
2
-
v
2
x
2
+
x
3
/
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
x
1
+
x
2
/
2
x
3
-
x
2
/
N
23
t
0
2
-
v
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
=
x
3
-
x
2
/
N
23
t
0
2
x
2
+
x
3
/
2
-
x
2
-
x
1
/
N
12
t
0
2
x
1
+
x
2
/
2
x
3
-
x
2
/
N
23
t
0
2
-
x
2
-
x
1
/
N
12
t
0
2
.That is, thex 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 (
x
3
-
x
2
)
/
N
23
t
0
≈
(
x
2
-
x
1
)
/
N
12
t
0
≈
v, and therefore x
→
+
∞ or x
→
-
∞.Equation (19) also shows that (
x
3
-
x
2
)
/
N
23
t
0
2
-
v
2
≥
0 and (
x
2
-
x
1
)
/
N
12
t
0
2
-
v
2
≥
0; therefore, the supremum of the domain of x
=
f
(
v
) is finite, denoted as v
m
a
x. The value of v
m
a
x depends on the smaller value between (
x
3
-
x
2
)
/
N
23
t
0 and (
x
2
-
x
1
)
/
N
12
t
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
m
a
x, the x coordinate of seismic source tends to the midpoint of the closer sensors from the seismic source, which is x
=
(
x
m
+
x
m
+
1
)
/
2.So, the graph ofx
=
f
(
v
) can be approximately plotted as shown in Figure 9.Figure 9
Approximate graph ofx
=
f
(
v
).
(a)
(b)There are two real cases ofx
=
f
(
v
) as shown in Figure 10.Figure 10
Real graphs ofx
=
f
(
v
).
(a)
X
=
f
v image of (−10, 0)~(−200, 0)
(b)
X
=
f
v image of (
-
10,0
) and (
-
110,0
)
(c)
X
=
f
v image of (10, 0)~(200, 0)
(d)
X
=
f
v image of (
10,0
) and (
110,0
)The parameters of (a) and (b) in Figure10 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 Figure10 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 Figure10, 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
1
,
y
1
), (
x
2
,
y
2
), and (
x
3
,
y
3
) are coordinates of a set of sensors and t
12 and t
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
2
,
y
2
), (
x
3
,
y
3
), and (
x
4
,
y
4
) are coordinates of a set of sensors and t
23 and t
34 are the corresponding arrival time differences; another equation will be got based on (19), denoted as x
′
=
f
(
v
).Step 3.
Definex
-
x
′ as the objective function; v
m
i
n which minimizes x
-
x
′ can be considered as true value of average wave velocity.Step 4.
Substitutev
m
i
n 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 Section3.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:(22)
The  normal-form  hyperbolic  equation  is
x
2
a
2
-
y
2
b
2
=
1
.
(23)
And  the  normal-form  asymptotic  equation  is
y
=
±
b
a
x
.With (22), the horizontal coordinate of a point on hyperbola is x
=
±
a
y
2
/
b
2
+
1, and the vertical coordinate is y
=
±
b
x
2
/
a
2
-
1. With (23), the horizontal coordinate of a point on asymptote is x
=
±
a
(
y
/
b
), and the vertical coordinate is y
=
±
b
(
x
/
a
):(24)
Define  the  relative  error  of  horizontal  coordinate  as
Δ
x
=
a
y
/
b
-
a
y
2
/
b
2
+
1
a
y
2
/
b
2
+
1
.
Define  the  relative  error  of  vertical  coordinate  as
Δ
y
=
b
x
/
a
-
b
x
2
/
a
2
-
1
b
x
2
/
a
2
-
1
.Δ
x is a decreasing function of y
/
b and lim
y
/
b
→
∞
⁡
Δ
x
=
0; Δ
y is a decreasing function of x
/
a and lim
x
/
a
→
∞
⁡
Δ
y
=
0. Some examples are shown in Table 1.Table 1
Relative error between hyperbola and asymptote.
y
/
b
Δ
x
x
/
a
Δ
y
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.007905261357939As Table1 shows, when y
/
b
>
7, there is Δ
x
<
1
%, and when x
/
a
>
7, there is Δ
y
<
1
%. Assume x
1
=
-
x
0
(
x
0
>
0
), x
2
=
x
0, and y
1
=
y
2
=
0; (9) will change into (25)
x
2
N
2
t
0
2
v
2
/
4
-
y
2
x
0
2
-
N
2
t
0
2
v
2
/
4
=
1
.In comparison with (22) and (25), the inequation y
/
b
>
7 corresponding to (22) will change into (26) corresponding to (25):(26)
y
x
0
2
-
N
2
t
0
2
v
2
/
4
>
7
.And the inequationx
/
a
>
7 will change into (27)
x
N
t
0
v
/
2
>
7
.Sincex
0
2
-
N
2
t
0
2
v
0
2
/
4
<
x
0
2, (26) can be rewritten to y
>
7
x
0
2
>
7
x
0
2
-
N
2
t
0
2
v
2
/
4; therefore, when y
>
3.5
·
2
x
0, there is Δ
x
<
1
%. Since N
t
0
v
<
2
x
0, (27) can be rewritten to x
>
7
x
0
>
7
(
N
t
0
v
/
2
); therefore, when x
>
3.5
·
2
x
0, there is Δ
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 y
>
2.5
·
2
x
0, there is Δ
x
<
2
%, and when x
>
2.5
·
2
x
0, there is Δ
y
<
2
%.According to the analysis above, the distance range of far field source depends on the demanded locating accuracy. For example, assumex
s and y
s are coordinates of seismic source, x
n
1 and y
n
1 are coordinates of the closest sensor, and x
n
2 and y
n
2 are coordinates of the next-closest sensor; if the demanded locating accuracy is less than 1%, then the seismic sources which satisfy x
s
-
(
x
n
1
+
x
n
2
)
/
2
>
3.5
x
n
2
-
x
n
1 and y
s
-
(
y
n
1
+
y
n
2
)
/
2
>
3.5
x
n
2
-
x
n
1 could be considered as far field sources.


## 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:(9)
x
-
x
1
+
x
2
/
2
2
N
2
t
0
2
v
2
/
4
-
y
-
y
1
+
y
2
/
2
2
x
2
-
x
1
/
2
2
-
N
2
t
0
2
v
2
/
4
=
1
,where (
x
1
,
y
1
) and (
x
2
,
y
2
) are coordinates of sensors, v is the average wave velocity, t
0 is the sampling period of A/D, and N
t
0 is the arrival time difference between (
x
1
,
y
1
) and (
x
2
,
y
2
); for convenience sake, assume N
t
0
>
0.The asymptotic equation of (9) is expressed in (10)
y
=
y
1
+
y
2
2
±
x
2
-
x
1
2
-
N
2
t
0
2
v
2
N
t
0
v
x
-
x
1
+
x
2
2
.The asymptote will be very close to the hyperbola when the seismic source locates far from sensors as shown in Figure1. 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 (11)
y
=
y
1
+
y
2
2
±
x
2
-
x
1
2
-
N
12
2
t
0
2
v
2
N
12
t
0
v
x
-
x
1
+
x
2
2
,
y
=
y
2
+
y
3
2
±
x
3
-
x
2
2
-
N
23
2
t
0
2
v
2
N
23
t
0
v
x
-
x
2
+
x
3
2
.DefineN
12
t
0
=
t
12 and N
23
t
0
=
t
23; the solutions of (11) are shown as follows:(12)
x
=
A
12
-
A
23
+
B
23
C
23
-
B
12
C
12
B
23
-
B
12
,
y
=
A
12
B
23
-
B
12
A
23
+
B
12
B
23
C
23
-
B
12
B
23
C
12
B
23
-
B
12
,
(13)
x
=
A
23
-
A
12
+
B
23
C
23
+
B
12
C
12
B
23
+
B
12
,
y
=
A
12
B
23
+
B
12
A
23
+
B
12
B
23
C
23
-
B
12
B
23
C
12
B
23
+
B
12
,
(14)
x
=
A
12
-
A
23
+
B
23
C
23
+
B
12
C
12
B
23
+
B
12
,
y
=
A
12
B
23
+
B
12
A
23
-
B
12
B
23
C
23
+
B
12
B
23
C
12
B
23
+
B
12
,
(15)
x
=
A
12
-
A
23
-
B
23
C
23
+
B
12
C
12
B
12
-
B
23
,
y
=
A
12
B
23
-
B
12
A
23
-
B
12
B
23
C
23
+
B
12
B
23
C
12
B
23
-
B
12
,where(16)
A
12
=
y
1
+
y
2
2
,
B
12
=
x
2
-
x
1
2
-
N
12
2
t
0
2
v
2
N
12
t
0
v
,
C
12
=
x
1
+
x
2
2
,
A
23
=
y
2
+
y
3
2
,
B
23
=
x
3
-
x
2
2
-
N
23
2
t
0
2
v
2
N
23
t
0
v
,
C
23
=
x
2
+
x
3
2
.All they coordinates are equal when sensors are placed in a line as shown in Figure 1. Therefore, A
12
=
A
23, and there are only two x coordinates in (12)–(15) shown as follows:(17)
x
=
B
23
C
23
-
B
12
C
12
B
23
-
B
12
,
(18)
x
=
B
23
C
23
+
B
12
C
12
B
23
+
B
12
.Equation (17) is the x coordinate of far field source. Substitute B
12, C
12, B
23, and C
23 into (17) and get (19)
x
=
x
3
-
x
2
/
N
23
t
0
2
-
v
2
x
2
+
x
3
/
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
x
1
+
x
2
/
2
x
3
-
x
2
/
N
23
t
0
2
-
v
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
.Equation (19) shows that if x
1, x
2, and x
3 (coordinates of sensors) and N
12
t
0 and N
23
t
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:(20)
d
x
d
v
=
v
/
2
B
-
A
B
+
A
x
1
+
x
2
-
x
2
-
x
3
A
B
B
-
A
2
,where A
=
(
(
x
2
-
x
1
)
/
N
12
t
0
)
2
-
v
2 and B
=
(
(
x
3
-
x
2
)
/
N
23
t
0
)
2
-
v
2.Equation (20) shows that the sign of d
x
/
d
v, namely, the monotonicity of x
=
f
(
v
), depends on the sign of (
B
-
A
)
(
x
1
-
x
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 Figure8.Figure 8
Relative positions of seismic source and sensors.As ① shows in Figure8, x
1
-
x
3
<
0, (
x
2
-
x
1
)
/
N
12
t
0
2
>
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
>
0, and x is a monotonous increasing function of v.As ② shows in Figure8, x
1
-
x
3
<
0, (
x
2
-
x
1
)
/
N
12
t
0
2
<
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
<
0, and x is a monotonous decreasing function of v.As ③ shows in Figure8, x
1
-
x
3
>
0, (
x
2
-
x
1
)
/
N
12
t
0
2
>
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
>
0, and x is a monotonous increasing function of v.As ④ shows in Figure8, x
1
-
x
3
>
0, (
x
2
-
x
1
)
/
N
12
t
0
2
<
(
x
3
-
x
2
)
/
N
23
t
0
2; therefore, d
x
/
d
v
=
(
B
-
A
)
(
x
1
-
x
3
)
<
0, and x is a monotonous decreasing function of v.According to (19), when v
→
0,(21)
lim
v
→
0
⁡
x
3
-
x
2
/
N
23
t
0
2
-
v
2
x
2
+
x
3
/
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
x
1
+
x
2
/
2
x
3
-
x
2
/
N
23
t
0
2
-
v
2
-
x
2
-
x
1
/
N
12
t
0
2
-
v
2
=
x
3
-
x
2
/
N
23
t
0
2
x
2
+
x
3
/
2
-
x
2
-
x
1
/
N
12
t
0
2
x
1
+
x
2
/
2
x
3
-
x
2
/
N
23
t
0
2
-
x
2
-
x
1
/
N
12
t
0
2
.That is, thex 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 (
x
3
-
x
2
)
/
N
23
t
0
≈
(
x
2
-
x
1
)
/
N
12
t
0
≈
v, and therefore x
→
+
∞ or x
→
-
∞.Equation (19) also shows that (
x
3
-
x
2
)
/
N
23
t
0
2
-
v
2
≥
0 and (
x
2
-
x
1
)
/
N
12
t
0
2
-
v
2
≥
0; therefore, the supremum of the domain of x
=
f
(
v
) is finite, denoted as v
m
a
x. The value of v
m
a
x depends on the smaller value between (
x
3
-
x
2
)
/
N
23
t
0 and (
x
2
-
x
1
)
/
N
12
t
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
m
a
x, the x coordinate of seismic source tends to the midpoint of the closer sensors from the seismic source, which is x
=
(
x
m
+
x
m
+
1
)
/
2.So, the graph ofx
=
f
(
v
) can be approximately plotted as shown in Figure 9.Figure 9
Approximate graph ofx
=
f
(
v
).
(a)
(b)There are two real cases ofx
=
f
(
v
) as shown in Figure 10.Figure 10
Real graphs ofx
=
f
(
v
).
(a)
X
=
f
v image of (−10, 0)~(−200, 0)
(b)
X
=
f
v image of (
-
10,0
) and (
-
110,0
)
(c)
X
=
f
v image of (10, 0)~(200, 0)
(d)
X
=
f
v image of (
10,0
) and (
110,0
)The parameters of (a) and (b) in Figure10 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 Figure10 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 Figure10, 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
1
,
y
1
), (
x
2
,
y
2
), and (
x
3
,
y
3
) are coordinates of a set of sensors and t
12 and t
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
2
,
y
2
), (
x
3
,
y
3
), and (
x
4
,
y
4
) are coordinates of a set of sensors and t
23 and t
34 are the corresponding arrival time differences; another equation will be got based on (19), denoted as x
′
=
f
(
v
).Step 3.
Definex
-
x
′ as the objective function; v
m
i
n which minimizes x
-
x
′ can be considered as true value of average wave velocity.Step 4.
Substitutev
m
i
n 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 Section3.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:(22)
The  normal-form  hyperbolic  equation  is
x
2
a
2
-
y
2
b
2
=
1
.
(23)
And  the  normal-form  asymptotic  equation  is
y
=
±
b
a
x
.With (22), the horizontal coordinate of a point on hyperbola is x
=
±
a
y
2
/
b
2
+
1, and the vertical coordinate is y
=
±
b
x
2
/
a
2
-
1. With (23), the horizontal coordinate of a point on asymptote is x
=
±
a
(
y
/
b
), and the vertical coordinate is y
=
±
b
(
x
/
a
):(24)
Define  the  relative  error  of  horizontal  coordinate  as
Δ
x
=
a
y
/
b
-
a
y
2
/
b
2
+
1
a
y
2
/
b
2
+
1
.
Define  the  relative  error  of  vertical  coordinate  as
Δ
y
=
b
x
/
a
-
b
x
2
/
a
2
-
1
b
x
2
/
a
2
-
1
.Δ
x is a decreasing function of y
/
b and lim
y
/
b
→
∞
⁡
Δ
x
=
0; Δ
y is a decreasing function of x
/
a and lim
x
/
a
→
∞
⁡
Δ
y
=
0. Some examples are shown in Table 1.Table 1
Relative error between hyperbola and asymptote.
y
/
b
Δ
x
x
/
a
Δ
y
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.007905261357939As Table1 shows, when y
/
b
>
7, there is Δ
x
<
1
%, and when x
/
a
>
7, there is Δ
y
<
1
%. Assume x
1
=
-
x
0
(
x
0
>
0
), x
2
=
x
0, and y
1
=
y
2
=
0; (9) will change into (25)
x
2
N
2
t
0
2
v
2
/
4
-
y
2
x
0
2
-
N
2
t
0
2
v
2
/
4
=
1
.In comparison with (22) and (25), the inequation y
/
b
>
7 corresponding to (22) will change into (26) corresponding to (25):(26)
y
x
0
2
-
N
2
t
0
2
v
2
/
4
>
7
.And the inequationx
/
a
>
7 will change into (27)
x
N
t
0
v
/
2
>
7
.Sincex
0
2
-
N
2
t
0
2
v
0
2
/
4
<
x
0
2, (26) can be rewritten to y
>
7
x
0
2
>
7
x
0
2
-
N
2
t
0
2
v
2
/
4; therefore, when y
>
3.5
·
2
x
0, there is Δ
x
<
1
%. Since N
t
0
v
<
2
x
0, (27) can be rewritten to x
>
7
x
0
>
7
(
N
t
0
v
/
2
); therefore, when x
>
3.5
·
2
x
0, there is Δ
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 y
>
2.5
·
2
x
0, there is Δ
x
<
2
%, and when x
>
2.5
·
2
x
0, there is Δ
y
<
2
%.According to the analysis above, the distance range of far field source depends on the demanded locating accuracy. For example, assumex
s and y
s are coordinates of seismic source, x
n
1 and y
n
1 are coordinates of the closest sensor, and x
n
2 and y
n
2 are coordinates of the next-closest sensor; if the demanded locating accuracy is less than 1%, then the seismic sources which satisfy x
s
-
(
x
n
1
+
x
n
2
)
/
2
>
3.5
x
n
2
-
x
n
1 and y
s
-
(
y
n
1
+
y
n
2
)
/
2
>
3.5
x
n
2
-
x
n
1 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 Section3.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.Table 2
Parameters of SGA.
Number of individuals
Maximum number of generations
Variable number of binary digits
Gap of generations
Crossover probability
Mutation probability
100
100
25
0.95
0.7
0.1


### 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
∈
(
60,960
) or x
∈
(
1010,2010
) and y
∈
(
10,1010
). Define relative error formula of v, x, and y as (28)
R
E
=
I
R
-
T
V
T
V
×
100
%
,where RE denotes the relative error, IR denotes the inversion result, and TV denotes the true value.The LEFT (means thex coordinate of seismic source satisfies x
∈
(
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
∈
(
1010,2010
)) errors are shown in Figure 12.Figure 11
Relative error of LEFTV
/
X
/
Y.
(a)
Relative error of LEFTV
(b)
Relative error of LEFTX
(c)
Relative error of LEFTYFigure 12
Relative Error of RIGHTV
/
X
/
Y.
(a)
Relative error of RIGHTV
(b)
Relative error of RIGHTX
(c)
Relative error of RIGHTYAll the comparatively larger relative errors highlighted in Figures11 and 12 are shown in Table 3.Table 3
Values and coordinates of comparatively larger relative errors.
Error name
Relative error (%)
True value ofX coordinate (m)
True value ofY coordinate (m)
LEFTV
0.02626
960
10
LEFTV
−0.08622
960
110
LEFTX
0.08442
960
10
LEFTX
0.05909
60
110
LEFTX
−0.1014
60
10
LEFTY
0.1885
960
110
LEFTY
−1.72
960
10
RIGHTV
0.02626
1010
10
RIGHTV
−0.08622
1010
110
RIGHTX
−0.08024
1010
10
RIGHTY
0.1885
1010
110
RIGHTY
−1.72
1010
10As shown in Table3, 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
21
≈
1.1704 ms, t
32
≈
1.1976 ms, and t
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 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.Table 4
Results of contrast experiments.
Experiment number
Iteration number of OLD method
Optimumv of OLD method
Error of OLD method betweenv and true value
Iteration number of NEW method
Optimumv of NEW method
Error of NEW method betweenv and true value
1
100
2569.7434
−430.2566
17
2614.70
−385.30
2
100
2665.3755
−334.6245
13
2614.70
−385.30
3
90
2306.9505
−693.0495
18
2614.70
−385.30
4
100
2237.0807
−762.9193
16
2614.70
−385.30
5
71
2205.2818
−794.7182
13
2614.70
−385.30
6
99
2660.2377
−339.7623
6
2614.70
−385.30
7
100
2476.904
−523.096
15
2614.70
−385.30
8
100
2204.2403
−795.7597
14
2614.70
−385.30
9
67
2474.3871
−525.6129
17
2614.70
−385.30
10
35
2165.6936
−834.3064
16
2614.70
−385.30
11
100
2573.8501
−426.1499
25
2614.70
−385.30
12
75
2475.6581
−524.3419
19
2614.70
−385.30
13
80
1756.5821
−1243.4179
31
2614.70
−385.30
14
96
1759.8416
−1240.1584
14
2614.70
−385.30
15
30
2276.137
−723.863
19
2614.70
−385.30
16
87
1846.0683
−1153.9317
14
2614.70
−385.30
17
41
2781.9136
−218.0864
13
2614.70
−385.30
18
60
1756.4908
−1243.5092
8
2614.70
−385.30
19
64
3280.5975
280.5975
13
2614.70
−385.30
20
99
2343.6613
−656.3387
25
2614.70
−385.30Comparison chart of convergence speed is shown in Figure13.Figure 13
Comparison chart of convergence speed.Obviously, Figure13 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 Figure14, and clearly the optimal solution of NEW method is more stable than OLD method.Figure 14
Stability of optimal solution.According to Table4, 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
21
=
1.170400 ms, t
32
=
1.197628 ms, and t
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
21
=
1.170400227006
×
10
-
3 s, t
32
=
1.197627771071
×
10
-
3 s, and t
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.Figure 15
Calculated results by 1–200 dB white Gauss noise added arrival time differences.Figure 16
Calculated results by 80–200 dB white Gauss noise added arrival time differences.


## 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 Section3.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.Table 2
Parameters of SGA.
Number of individuals
Maximum number of generations
Variable number of binary digits
Gap of generations
Crossover probability
Mutation probability
100
100
25
0.95
0.7
0.1


## 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
∈
(
60,960
) or x
∈
(
1010,2010
) and y
∈
(
10,1010
). Define relative error formula of v, x, and y as (28)
R
E
=
I
R
-
T
V
T
V
×
100
%
,where RE denotes the relative error, IR denotes the inversion result, and TV denotes the true value.The LEFT (means thex coordinate of seismic source satisfies x
∈
(
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
∈
(
1010,2010
)) errors are shown in Figure 12.Figure 11
Relative error of LEFTV
/
X
/
Y.
(a)
Relative error of LEFTV
(b)
Relative error of LEFTX
(c)
Relative error of LEFTYFigure 12
Relative Error of RIGHTV
/
X
/
Y.
(a)
Relative error of RIGHTV
(b)
Relative error of RIGHTX
(c)
Relative error of RIGHTYAll the comparatively larger relative errors highlighted in Figures11 and 12 are shown in Table 3.Table 3
Values and coordinates of comparatively larger relative errors.
Error name
Relative error (%)
True value ofX coordinate (m)
True value ofY coordinate (m)
LEFTV
0.02626
960
10
LEFTV
−0.08622
960
110
LEFTX
0.08442
960
10
LEFTX
0.05909
60
110
LEFTX
−0.1014
60
10
LEFTY
0.1885
960
110
LEFTY
−1.72
960
10
RIGHTV
0.02626
1010
10
RIGHTV
−0.08622
1010
110
RIGHTX
−0.08024
1010
10
RIGHTY
0.1885
1010
110
RIGHTY
−1.72
1010
10As shown in Table3, 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
21
≈
1.1704 ms, t
32
≈
1.1976 ms, and t
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 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.Table 4
Results of contrast experiments.
Experiment number
Iteration number of OLD method
Optimumv of OLD method
Error of OLD method betweenv and true value
Iteration number of NEW method
Optimumv of NEW method
Error of NEW method betweenv and true value
1
100
2569.7434
−430.2566
17
2614.70
−385.30
2
100
2665.3755
−334.6245
13
2614.70
−385.30
3
90
2306.9505
−693.0495
18
2614.70
−385.30
4
100
2237.0807
−762.9193
16
2614.70
−385.30
5
71
2205.2818
−794.7182
13
2614.70
−385.30
6
99
2660.2377
−339.7623
6
2614.70
−385.30
7
100
2476.904
−523.096
15
2614.70
−385.30
8
100
2204.2403
−795.7597
14
2614.70
−385.30
9
67
2474.3871
−525.6129
17
2614.70
−385.30
10
35
2165.6936
−834.3064
16
2614.70
−385.30
11
100
2573.8501
−426.1499
25
2614.70
−385.30
12
75
2475.6581
−524.3419
19
2614.70
−385.30
13
80
1756.5821
−1243.4179
31
2614.70
−385.30
14
96
1759.8416
−1240.1584
14
2614.70
−385.30
15
30
2276.137
−723.863
19
2614.70
−385.30
16
87
1846.0683
−1153.9317
14
2614.70
−385.30
17
41
2781.9136
−218.0864
13
2614.70
−385.30
18
60
1756.4908
−1243.5092
8
2614.70
−385.30
19
64
3280.5975
280.5975
13
2614.70
−385.30
20
99
2343.6613
−656.3387
25
2614.70
−385.30Comparison chart of convergence speed is shown in Figure13.Figure 13
Comparison chart of convergence speed.Obviously, Figure13 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 Figure14, and clearly the optimal solution of NEW method is more stable than OLD method.Figure 14
Stability of optimal solution.According to Table4, 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
21
=
1.170400 ms, t
32
=
1.197628 ms, and t
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
21
=
1.170400227006
×
10
-
3 s, t
32
=
1.197627771071
×
10
-
3 s, and t
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.Figure 15
Calculated results by 1–200 dB white Gauss noise added arrival time differences.Figure 16
Calculated results by 80–200 dB white Gauss noise added arrival time differences.


## 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 ofx
=
f
(
v
) in Section 3.1.
(2)
The relative errors shown in Figures11 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.



---

*Source: 8540303-2016-04-26.xml*