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.

1. Introduction

Microseismic monitoring, which was first studied by Obert [1], has been successfully applied in mining [25], oil exploitation [610], water conservancy construction [1113], and so forth. One of the core technologies of microseismic monitoring is source locating method, which can be reduced to (1) in homogeneous medium:where , , and are coordinates of seismic source and is the initial time of source; , , and are coordinates of the th sensor; is the average wave velocity measured by the th sensor; is arrival time measured by th sensor. It is obvious that 4 or more sensors at least are needed to solve for if 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 [1517], 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 [2628]. 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 [2932], 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.

, , , and are coordinates of equidistant placed sensors, are coordinates of seismic source, is the average wave velocity, is the arrival time difference between and , and and are similar to . The unknowns are and . If is known, can be solved by the following equation set:

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

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

Equation sets (3) and (4) are both nonlinear equations, and they can be solved by Simplex Method, GA, and other global optimization methods.

According to Prugger’s method, the objective function of (3) can be defined as (5) and the objective function of (4) can be defined as (6):

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

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 , , , and ; the average wave velocity is 3000 m/s and the coordinate of seismic source is . The graph of objective function (6) when and is shown in Figure 2.

Enlarged view of Figure 2 near 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 ) is shown in Figure 5.

As shown in Figures 25, two concluding points can be made when the seismic source is far away from sensors as shown in Figure 1; that is, the locating accuracy of coordinate is better than coordinate and the objective function, whether (6) or (7), changes little near true value. As a comparison, the graphs of (6) and (7) when seismic source locates at 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 satisfy hyperbolic equation (9), shown as follows:where and are coordinates of sensors, is the average wave velocity, is the sampling period of A/D, and is the arrival time difference between and ; for convenience sake, assume .

The asymptotic equation of (9) is expressed in

The asymptote will be very close to the hyperbola when the seismic source locates far from sensors as shown in Figure 1. If the error between asymptote and hyperbola is less than predefined locating error, (10) can be used to locate instead of (9) and we will obtain an approximate equation set instead of (2) as shown in

Define and ; the solutions of (11) are shown as follows:where

All the coordinates are equal when sensors are placed in a line as shown in Figure 1. Therefore, , and there are only two coordinates in (12)–(15) shown as follows:

Equation (17) is the coordinate of far field source. Substitute , , , and into (17) and get

Equation (19) shows that if , , and (coordinates of sensors) and and (arrival time differences) are known, ( coordinate of seismic source) can be regarded as a function of (average wave velocity), denoted as . Obviously, there is a one-to-one correspondence between and since in practice; that is to say, exists and is unique when is known.

An approximate graph can be plotted by (19).

Take the derivative of (19) with respect to :where and .

Equation (20) shows that the sign of , namely, the monotonicity of , depends on the sign of 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 ① shows in Figure 8, , ; therefore, , and is a monotonous increasing function of .

As ② shows in Figure 8, , ; therefore, , and is a monotonous decreasing function of .

As ③ shows in Figure 8, , ; therefore, , and is a monotonous increasing function of .

As ④ shows in Figure 8, , ; therefore, , and is a monotonous decreasing function of .

According to (19), when ,

That is, the coordinate of seismic source tends to a constant when , 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 , and therefore or .

Equation (19) also shows that and ; therefore, the supremum of the domain of is finite, denoted as . The value of depends on the smaller value between and ; 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 , the coordinate of seismic source tends to the midpoint of the closer sensors from the seismic source, which is .

So, the graph of can be approximately plotted as shown in Figure 9.

There are two real cases of as shown in Figure 10.

The parameters of (a) and (b) in Figure 10 are shown as follows.

The coordinates of sensors are and the interval of sensors is 1 meter. The average wave velocity is  m/s. And the coordinate of seismic source is . The line nearly parallel to the horizontal axis in (b) is the graph of whose sensor coordinates are and while the other curve is the graph whose sensor coordinates are and .

The parameters of (c) and (d) in Figure 10 are shown as follows.

The coordinates of sensors are and the interval of sensors is 1 meter. The average wave velocity is  m/s. And the coordinate of seismic source is . The line nearly parallel to the horizontal axis in (d) is the graph whose sensor coordinates are and while the other curve is the graph whose sensor coordinates are and .

There is an inference we can make from Figure 10, which is that the intersection point of all the graphs exists and is unique; the reasons are as follows.

Existence. There must be at least one intersection point of the graphs; the horizontal coordinate of the intersection is the real average wave velocity and the vertical coordinate is the real horizontal coordinate of seismic source.

Uniqueness. Suppose there is another intersection, named , of two graphs. As mentioned above, when is known, (the horizontal coordinate of seismic source) exists and is unique; that is to say, when the average wave velocity is , all of the graphs will uniquely intersect at one point . 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 . Therefore, there is only one intersection of graphs.

According to the analyses above, the procedure of the new two-dimensional far field source locating method is as follows.

Step 1. Assume , , and are coordinates of a set of sensors and and 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 .

Step 2. Assume , , and are coordinates of a set of sensors and and are the corresponding arrival time differences; another equation will be got based on (19), denoted as .

Step 3. Define as the objective function; which minimizes can be considered as true value of average wave velocity.

Step 4. Substitute 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:

With (22), the horizontal coordinate of a point on hyperbola is , and the vertical coordinate is . With (23), the horizontal coordinate of a point on asymptote is , and the vertical coordinate is :

is a decreasing function of and ; is a decreasing function of and . Some examples are shown in Table 1.

As Table 1 shows, when , there is , and when , there is . Assume , , and ; (9) will change into

In comparison with (22) and (25), the inequation corresponding to (22) will change into (26) corresponding to (25):

And the inequation will change into

Since , (26) can be rewritten to ; therefore, when , there is . Since , (27) can be rewritten to ; therefore, when , there is . 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 , there is , and when , there is .

According to the analysis above, the distance range of far field source depends on the demanded locating accuracy. For example, assume and are coordinates of seismic source, and are coordinates of the closest sensor, and and are coordinates of the next-closest sensor; if the demanded locating accuracy is less than 1%, then the seismic sources which satisfy and 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 , , , and , 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 , plane at interval of 100 meters and calculate the arrival time (in seconds with 16 significant digits) of each point, where or and . Define relative error formula of , , and as where RE denotes the relative error, IR denotes the inversion result, and TV denotes the true value.

The LEFT (means the coordinate of seismic source satisfies ) relative errors of , , and are shown in Figure 11, and the RIGHT (means the coordinate of seismic source satisfies ) errors are shown in Figure 12.

All the comparatively larger relative errors highlighted in Figures 11 and 12 are shown in Table 3.

As shown in Table 3, the comparatively larger relative errors will arise when seismic source comes close to sensors which is consistent with the analysis in Section 3.2. For example, when the largest relative error of LEFT is −0.1014% and the corresponding coordinate of seismic source is 10, the distance between the coordinate of seismic source and the coordinate of sensors is equal to sensor interval. One more example is that when the largest relative error of LEFT is %, the corresponding coordinate of seismic source is 960, and the coordinate is 10, it is obvious that the distance between the coordinate of seismic source and the coordinate of sensors is equal to sensor interval and the minimum distance between the coordinate of seismic source and the coordinate of the midpoint of sensors is slightly greater than sensor interval. According to the numerical experimental results and analyses above, the locating result can be considered unauthentic when the minimum distance between the located result and sensors is less than 2.5 times the sensor interval.

4.3. Contrast Experiments

Assume the coordinates of sensors are , , , and , the coordinates of seismic source are , and the average wave velocity is 3000 m/s; the corresponding theoretical arrival time differences are  ms,  ms, and  ms after numerical simulation. Table 4 shows 20 contrast experiment results based on the coordinates of sensors and the theoretical arrival time differences calculated above, where NEW denotes the method taking as objective function and OLD denotes the method taking (6) as objective function. The optimization methods of NEW and OLD are both SGA and the parameters of SGA are shown in Table 2.

Comparison chart of convergence speed is shown in Figure 13.

Obviously, Figure 13 shows that the NEW method converges faster than the OLD method and also the NEW method’s fluctuation of iterative times is smaller than the OLD method.

The stability of optimal solution is shown in Figure 14, and clearly the optimal solution of NEW method is more stable than OLD method.

According to Table 4, neither optimal solution of the two methods converges to true value, which is caused by the poor timing accuracy. If we substitute arrival time differences with higher accuracy into the NEW method, such as  ms,  ms, and  ms, then the optimal solution of 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  s,  s, and  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 in Section 3.1.(2)The relative errors shown in Figures 11 and 12 and Table 3 are smaller than the theoretic errors got from Section 3.2. The maximum relative error of coordinate is −0.1014% (LEFT) and −0.08024% (RIGHT), respectively; they are both smaller than the theoretic error (1%) deduced in Section 3.2.(3)NEW method has weaker antinoise ability than OLD method.

The three questions mentioned above will need to be solved in future research.

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.