# An Undersea Mining Microseism Source Location Algorithm Considering Wave Velocity Probability Distribution

Authors: Furui Du; Nailian Hu; Yuling Xie; Guoqing Li
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The traditional mine microseism locating methods are mainly based on the assumption that the wave velocity is uniform through the space, which leads to some errors for the assumption goes against the laws of nature. In this paper, the wave velocity is regarded as a random variable, and the probability distribution information of the wave velocity is fused into the traditional locating method. This paper puts forwards the microseism source location method for the undersea mining on condition of the probability distribution of the wave velocity and comes up with the solving process of Monte Carlo. In addition, based on the simulated results of the Monte Carlo method, the space is divided into three areas: the most possible area (area I), the possible area (area II), and the small probability area (area III). Attached to corresponding mathematical formulations, spherical models and cylindrical models in different areas are, respectively, built according to whether the source is in the sensor arrays. Both the examples and the actual applications show that (1) the method of microseism source location in this paper can highly improve the accuracy of the microseism monitoring, especially for the source beyond the sensor arrays, and (2) the space-dividing method based on occurrence possibilities of the source can recognize and sweep the hidden dangers for it predicts the probable location range of the source efficiently, while the traditional method cannot.


---

## Body

## 1. Introduction

Jiaodong Peninsula is the most important gold production area in China and it produces about 25% gold in China with 0.2% national territorial area [1]. Xinli minefield of Sanshan Island gold mine is the only domestic seabed mining operation mine that located in Laizhou Bay of Jiaodong Peninsula [2]. Orebody of Xinli minefield and seawater are separated by a few-meter-thick sea mud confining bed. Rock mass stability damage caused by large-scale blasting operation may lead to rock deformation and fracture and further form water passages, which may cause seawater pouring and other serious hazards [3].It will make some vibrations spreading around during the rock fracture and crushing processes. According to this principle, microseism monitoring has become a common and effective location technique. The development of foreign microseism monitoring technology has turned the monitoring of mine microfracture from a wild wish into an integral part of mining safety management. On the other hand, more and more domestic mines have established microseism monitoring system. In 2004, Fankou lead-zinc mine established 16-channel microseism monitoring system to monitor the high ground stress [4]. In 2005, Dongguashan copper mine introduced ISS microseism monitoring system to monitor the rock burst [5]. In 2008, Wutongzhuang coal mine adopted high-precision monitoring system with independent intelligent property right for water-burst monitoring [6].Microseism source location is a core technology for microseism monitoring and it directly affects whether the microseism monitoring is a success or not. With the rapid development of computers, Geiger’s location thought [7] has been increasingly applied and developed and some improvements have been successively proposed at the same time [8–10], especially in the solution stability and the influence of the decreasing arrival time residual. Besides the classical Geiger location thought, multiple-event location, travel time difference localization, and other location methods [11, 12] are also developed. Some nonlinear means are used to solve the difficulties of Geiger linear method [13]. Most of the above-mentioned methods are established on the condition that the wave velocity has been known and kept steady in different propagation paths. However, numerous practices have proved that it is difficult to accurately measure the wave velocity and the complexity of rock composition in different propagation paths may lead to different wave velocities. Some studies have demonstrated that [14] the tiny velocity error may cause more than 100 m location error. Dong et al. regarded wave velocity as an unknown parameter and proposed one location method without measurement of seismic wave propagation velocity [14, 15], which was an efficient method to avoid the wave velocity measurement error. However, this also is based on wave velocity uniformity hypothesis in space. Due to danger and the severity of consequence, undersea mining not only needs accurate location method, but also needs screening of the safety hazard according to the location results. In this paper, the wave velocity in space is viewed as a random variable and the wave velocity probability distributions were introduced into microseism source location to simulate the actual condition of nonuniform wave velocity, hoping to divide the space by aid of the occurrence probability of seismic source location. In this way, the potential safety hazards can be identified and screened.


## 2. Location Method of Uniform Wave Velocity

Traditional hypocenter location methods are the ones with uniform and known wave velocity, generally based on the microseism location theory proposed by Geiger. Essentially, this traditional method is to evaluate the minimum value of target function. Due toP wave’s fastest propagation velocity among the seismic waves and its easily recognizable arrival time, the method of P wave location is adopted in most cases. The method of calibration could help predetermine the velocity of P wave, namely, V
P, and V
P is generally considered as a constant. Set the position of hypocenter as (
x
,
y
,
z
), origin time as t, the coordinates of the stations which have received seismic signals as (x
i
,
y
i
,
z
i) (
i
=
1,2
,
…
n
), and the time of receiving seismic signals as t
i
(
i
=
1,2
,
…
n
). Since the emphasis of this paper is to introduce the method of location, it is assumed that arrival time has been obtained automatically or manually, and thus the problem can be simplified to evaluate the minimum value of the below target function:
(1)
min
⁡
ϕ
(
t
,
x
,
y
,
z
)
=
∑
i
=
1
n
[
t
i
-
t
-
T
i
(
x
,
y
,
z
)
]
2
.
In the above formula, the value of T
i
(
x
,
y
,
z
) which represents the travel time from the hypocenter to Station i could be illustrated in three-dimensional space as follows:
(2)
T
i
(
x
,
y
,
z
)
=
l
i
V
P
.
where, l
i
=
(
x
-
x
i
)
2
+
(
y
-
y
i
)
2
+
(
z
-
z
i
)
2.Target function (1) contains four variables, the spatial coordinate parameters of hypocenter x
,
y
,
z and the origin time t. Generally speaking, the positioning requires the data from at least 4 stations. Dong et al. put forward a hypocenter location method with no need to measure the wave velocity in advance in [14], which properly solves the problem of inaccurate measurement of wave velocity. Compared with the traditional method, the similarity is to assume the wave is even, but this method is to regard wave velocity as unknown. The target function is
(3)
min
⁡
ϕ
(
t
,
x
,
y
,
z
,
V
P
)
=
∑
i
=
1
n
[
t
i
-
t
-
T
i
(
x
,
y
,
z
,
V
P
)
]
2
.
The travel time in this equation can be expressed as
(4)
T
i
(
x
,
y
,
z
)
=
l
i
V
P
.
Compared with (1), (3) contains one more unknown V
P to be determined, so the location by (3) needs the data from at least 5 stations.


## 3. Location Method Based on Distribution of Wave Velocity Probability

Whether (1) or (3) is adopted, it is assumed that the wave velocity from hypocenter to each station is the same, namely, the hypothesis of uniform velocity. However, due to the complexity of rock constituent and jointed rock mass on wave propagation path, wave velocity in different paths usually differs. This paper considers wave velocities on different paths as a random variable, following the same probability distribution. The probability density function is noted as f
(
T
P
) and cumulative distribution function is noted as F
(
T
P
).The probability distribution of wave velocity can be determined by the parameters from large amounts of the observed samples obtained by the method of calibration. The parameter of probability distribution can be determined by maximum likelihood method. Suppose the observed value of wave velocity sample is(
v
1
,
v
2
,
…
.
,
v
m
), the overall probability density function of wave velocity is f
(
v
i
,
θ
1
,
θ
2
,
…
,
θ
k
), and k unknown to-be-estimated parameters are θ
1
,
θ
2
,
…
,
θ
k. Therefore, the likelihood function can be noted as
(5)
L
(
θ
1
,
θ
2
,
…
,
θ
k
)
=
∏
i
=
1
m
f
(
v
i
,
θ
1
,
θ
2
,
…
,
θ
k
)
.
The estimated values of θ
1
,
θ
2
,
…
,
θ
k enable the likelihood function L
(
θ
1
,
θ
2
,
…
,
θ
k
) to reach the maximum value. The estimator of the maximum likelihood will be the solution of the following equation:
(6)
∂
L
(
θ
1
,
θ
2
,
…
,
θ
k
)
∂
θ
i
=
0
(
i
=
1,2
,
…
,
k
)
.
The distributed fitting quality after estimated parameter needs to be inspected byK-Sor χ
2, or another probability distribution will be selected. In practice, the frequently used model of probability distribution such as normal distribution or log-normal distribution could be adopted. A simple binary equation set could be obtained after the probability distribution model is substituted in (6). Use MATLAB to obtain the solution directly.After the proper probability distribution is determined, the theory of Monte Carlo could be used for location. Specifically speaking, a set of wave velocities can be obtained from the probability distribution determined by (5) and (6). Substitute the wave velocities into (2), and then substitute (2) into (1) to get the minimum value and obtain a set of the most optimal solutions  (
t
′
,
x
′
,
y
′
,
z
′
). After various repetitions, N sets of the most optimal solutions (
t
i
′
,
x
i
′
,
y
i
′
,
z
i
′
) would be obtained. Using (7), the hypocenter location on average can be obtained:
(7)
x
0
=
(
x
1
′
+
x
2
′
+
⋯
+
x
N
′
)
N
,
y
0
=
(
y
1
′
+
y
2
′
+
⋯
+
y
N
′
)
N
,
z
0
=
(
z
1
′
+
z
2
′
+
⋯
+
z
N
′
)
N
.
Since the location method of probability distribution of velocity requires repeated calculations, more calculation time would be needed than traditional method. However, the local optimization algorithm would be adopted in each calculation and the amount of calculation is little, and therefore the calculating time of this method is acceptable. The flow of the location algorithm that considers the probability distribution of wave velocity is depicted in Figure 1.Figure 1
The flowchart of location algorithm presented in this paper.


## 4. Space Partition of Hypocenter Location Based on Probability

To partition space into several areas according to the probability of the existence of hypocenter is helpful to search hypocenter precisely and identify the types of hypocenter. The large amount of simulations manifests; as for the hypocenter in the array of sensors, the results of Monte Carlo simulations(
x
i
′
,
y
i
′
,
z
i
′
) clustered near one point; see Figure 2(a); as for the hypocenter out of the array of sensors, (
x
i
′
,
y
i
′
,
z
i
′
) spreads along a line in band distribution; see Figure 2(b). Therefore, globular model shows the probability distribution of hypocenter in the array of sensors, and cylinder model shows the probability distribution of hypocenter out of the array of sensors. The mathematical versions of these two models are as follows.Space partition based on the probability of seismic focus: (a) globular model and (b) cylinder model.
(a)
(b)


### 4.1. Globular Model

The cumulative probability functionφ
S
(
r
) of hypocenter distance in globular model can be defined as
(8)
φ
S
(
r
)
=
∑
i
=
1
N
H
S
i
(
r
)
N
.
The probability function of H
S
i
(
r
) is
(9)
H
S
i
(
r
)
=
{
1
(
x
i
′
-
x
0
)
2
+
(
y
i
′
-
y
0
)
2
+
(
z
i
′
-
z
0
)
2
≤
r
0
else
.
The value of φ
S
(
r
) represents the proportion of points in Monte Carlo simulation results whose space distances to the point (
x
0
,
y
0
,
z
0
) are no more than r. According to the law of large numbers, when the simulation times are enough, φ
S
(
r
) is close to the probability that the distance of the actual hypocenter to the calculated one is less than r. When the value of φ
S
(
r
) is 60%, the corresponding r means the probability that the distance of the hypocenter to the point (
x
0
,
y
0
,
z
0
) is no more than r is 60%, noted as R
60. R
95 could be defined likewise. According to the distance between space point and (
x
0
,
y
0
,
z
0
), the space is divided into I, II, and III. I represents the area whose probability is less than R
60, and the probability of the existence of hypocenter is the highest in this area; II represents the area whose probability is between R
60 and R
95, representing the probability of the existence of hypocenter which is relatively higher in this area; III represents the location whose probability is more than R
95, representing the probability of the existence of hypocenter which is 5%, namely, the small probability event; see Figure 2(a).


### 4.2. Cylinder Model

Similar to the problem of linear least square fitting on surface, the obtained fitting straight-line(
x
i
′
,
y
i
′
,
z
i
′
) in space can be written as below:
(10)
x
-
x
0
m
=
y
-
y
0
n
=
z
-
z
0
p
.
Cumulative probability function φ
C
(
r
) of hypocenter distance in cylinder model can be defined as
(11)
φ
C
(
r
)
=
∑
i
=
1
N
H
C
i
(
r
)
N
.
The H
C
i
(
r
) in the equation is
(12)
H
C
i
(
r
)
=
{
1
,
|
i
j
k
x
0
-
x
i
′
y
0
-
y
i
z
0
-
z
i
′
m
n
p
|
m
2
+
n
2
+
p
2
≤
r
,
0
,
else
.
The value of φ
c
(
r
) represents the proportion of points in simulation results whose distances to the fitting straight-line (10) are no more than r. Similar to globular model, R
60, R
95, I, II, and III can be defined; see Figure 2(b).


## 4.1. Globular Model

The cumulative probability functionφ
S
(
r
) of hypocenter distance in globular model can be defined as
(8)
φ
S
(
r
)
=
∑
i
=
1
N
H
S
i
(
r
)
N
.
The probability function of H
S
i
(
r
) is
(9)
H
S
i
(
r
)
=
{
1
(
x
i
′
-
x
0
)
2
+
(
y
i
′
-
y
0
)
2
+
(
z
i
′
-
z
0
)
2
≤
r
0
else
.
The value of φ
S
(
r
) represents the proportion of points in Monte Carlo simulation results whose space distances to the point (
x
0
,
y
0
,
z
0
) are no more than r. According to the law of large numbers, when the simulation times are enough, φ
S
(
r
) is close to the probability that the distance of the actual hypocenter to the calculated one is less than r. When the value of φ
S
(
r
) is 60%, the corresponding r means the probability that the distance of the hypocenter to the point (
x
0
,
y
0
,
z
0
) is no more than r is 60%, noted as R
60. R
95 could be defined likewise. According to the distance between space point and (
x
0
,
y
0
,
z
0
), the space is divided into I, II, and III. I represents the area whose probability is less than R
60, and the probability of the existence of hypocenter is the highest in this area; II represents the area whose probability is between R
60 and R
95, representing the probability of the existence of hypocenter which is relatively higher in this area; III represents the location whose probability is more than R
95, representing the probability of the existence of hypocenter which is 5%, namely, the small probability event; see Figure 2(a).


## 4.2. Cylinder Model

Similar to the problem of linear least square fitting on surface, the obtained fitting straight-line(
x
i
′
,
y
i
′
,
z
i
′
) in space can be written as below:
(10)
x
-
x
0
m
=
y
-
y
0
n
=
z
-
z
0
p
.
Cumulative probability function φ
C
(
r
) of hypocenter distance in cylinder model can be defined as
(11)
φ
C
(
r
)
=
∑
i
=
1
N
H
C
i
(
r
)
N
.
The H
C
i
(
r
) in the equation is
(12)
H
C
i
(
r
)
=
{
1
,
|
i
j
k
x
0
-
x
i
′
y
0
-
y
i
z
0
-
z
i
′
m
n
p
|
m
2
+
n
2
+
p
2
≤
r
,
0
,
else
.
The value of φ
c
(
r
) represents the proportion of points in simulation results whose distances to the fitting straight-line (10) are no more than r. Similar to globular model, R
60, R
95, I, II, and III can be defined; see Figure 2(b).


## 5. Example Verification

Reference [14] offers an example of location system, and this paper takes this example to verify the location method considering wave velocity probability distribution. The example is a cubic location system with 8 sensors, sitting at 8 peaks of the cube, respectively, and their code numbers range from A to H. The length of the cube is 800 m. There are 6 events occurring. Four of the microseismic sources are located in the array of sensors, and their coordinates are O(258, 336, 580), P(680, 290, 559), Q(190, 610, 380), and R(789, 459, 280), respectively, while two are out of the array of sensors, and their coordinates are S(308, 689, 1200) and T(860, 910, 1008). In [14], the real wave velocity of the example is 5.6 m/ms, while in fact the wave velocity of each event on different paths varies. Therefore, the example in this paper took the mean velocity from each seismic source to sensor, 5.6 m/ms, and it followed the normal distribution N(5.6,0.1). Suppose the trigger time of source is 10 ms. Each travel time is added with one white noise to simulate the measurement error and interpretation error. The trigger times of P wave recorded by sensors are shown in Table 1.Table 1
Trigger times of seismic source recorded by sensors.
Sensors
Trigger times of seismic source/ms
O
P
Q
R
S
T
A
140.217
176.360
138.273
182.867
269.607
303.767
B
163.213
123.338
177.213
102.220
282.374
268.584
C
172.747
144.828
139.224
86.881
248.827
185.928
D
154.237
196.822
91.167
172.020
238.873
246.680
E
95.277
151.204
146.806
202.678
167.463
238.378
F
128.200
80.245
181.033
127.277
175.813
169.173
G
140.214
120.084
141.824
121.567
124.405
51.801
H
108.838
179.073
97.952
190.028
101.942
166.961The new location method proposed in this paper is the traditional method and the method in [14] which are adopted for location. The location results are shown in Table 2. It can be seen from Table 2 that the location precision of the method proposed in this paper is better than the other two methods, especially the location of the seismic focus out of the array of sensors. In the location of seismic focus S, the error of the traditional method reaches 176 m and the error of the method in [14] is more than 40 m also. However, the error of the method proposed in this paper is merely 12.7 m. In the location of seismic focus T, the errors of the other two methods are over 200 m, while the accuracy of the method in this paper is improved by nearly 100 m, and the stability is enhanced greatly.Table 2
Result comparisons of the new method and traditional methods.
Method
Velocity
Seismic sourceO
Seismic sourceP
m×ms−1
X/m
Y/m
Z/m
D
error/m
X/m
Y/m
Z/m
D
error/m
New method
—
276.070
345.249
584.666
20.829
681.532
274.086
531.022
32.659
Method in Literature [14]
—
258.780
337.458
612.879
32.920
687.257
271.532
532.471
33.129
Traditional methods
5.600
268.941
341.609
597.019
20.995
701.997
265.732
538.570
38.602
Method
Velocity
Seismic sourceQ
Seismic sourceR
m×ms−1
X/m
Y/m
Z/m
D
error/m
X/m
Y/m
Z/m
D
error/m
New method
—
195.787
605.793
375.907
8.242
794.406
448.448
284.529
11.078
Method in Literature [14]
—
178.624
615.321
368.125
17.284
780.302
445.235
268.369
20.010
Traditional methods
5.600
194.751
615.690
374.316
9.341
823.755
452.297
276.689
35.550
Method
Velocity
Seismic sourceS
Seismic sourceT
m×ms−1
X/m
Y/m
Z/m
D
error/m
X/m
Y/m
Z/m
D
error/m
New method
—
319.658
683.683
1195.697
12.742
920.053
995.742
1110.136
140.355
Method in Literature [14]
—
315.320
699.656
1238.86
40.952
955.599
1040.270
1153.380
217.359
Traditional methods
5.600
300.376
749.555
1365.34
176.243
978.467
1063.330
1175.090
255.857Space partition based on the probability of seismic source is utilized to divide seismic sources (Table3). It can be seen that all the seismic sources fall in I, except seismic source P in II. This indicates that the partition method in this paper is reasonable. It is noteworthy that the location error of seismic source T is 140.355 m (Table 2), while the distance from seismic source to the axis of cylinder model is only 31.271 m. This indicates that the zonal distribution of out-of-array seismic source caused by velocity nonuniformity is an important factor of significant location error. This also may be the reason for the failure of traditional methods. Space partition model proposed in this paper can efficiently solve this problem.Table 3
Space partition for the probability of seismic source locations.
Seismic source
Partition model
R
60
R
95
Real distance
Partition area of seismic source
m
m
m
—
O
Globular model
21.036
32.869
20.829
I
P
Globular model
20.085
37.790
32.659
II
Q
Globular model
20.027
33.186
17.584
I
R
Globular model
24.115
38.256
11.078
I
S
Cylinder model
45.698
87.269
10.258
I
T
Cylinder model
38.652
68.691
31.271
I


## 6. Project Application and Analysis

Aiming at the requirements in which Xinli minefield needs high-precision and large-area monitoring, the microseismic monitoring system adopts the systematic framework combining the collective style and distribution style, namely, the arrangement of distribution style for areas and collective style between areas. It includes mainly 4 downhole monitoring substations, 1 downhole monitoring master-station, and ground data processing system. The system architecture is shown in Figure3.Figure 3
Microseism system architecture in Xinli minefield.The types of sensor are 4.5 HZ low-frequency high-sensitivity magnetic-electric sensor and 60 HZ high-frequency high-sensitivity magnetic-electric sensor, which could not only monitor the rock mass stability in small-scale quarry but also monitor the fracture and motion of hanging-wall rock mass caused by mining. Sensors are distributed at −135, −165, −200, and −400 and their coordinates are shown in Table4.Table 4
Coordinates of sensors and triggered times by sensors.
Number
Sensor’s coordinate/m
Arrival time of event/ms
X
Y
Z
1st event
2nd event
M1
94478.3
40514.9
−134
Nontriggered
Nontriggered
M2
95118.6
40784.6
−164
Nontriggered
4076
M3
95001.3
40717.4
−164
Nontriggered
4066
M4
94642.1
40555.1
−164
4081
Nontriggered
M5
94357.2
40427.6
−164
4094
Nontriggered
M6
94262.1
40392.8
−164
4103
Nontriggered
M7
93767.6
40316.0
−164
Nontriggered
Nontriggered
M8
93877.8
40287.5
−164
Nontriggered
Nontriggered
M9
95104.8
40732.8
−199
Nontriggered
4063
M10
95014.8
40680.8
−199
Nontriggered
4056
M11
94968.6
40658.9
−199
Nontriggered
4062
M12
94556.8
40474.8
−199
4062
Nontriggered
M13
94452.4
40433.5
−199
4056
Nontriggered
M14
94319.6
40384.6
−199
Nontriggered
Nontriggered
M15
94420.7
40256.0
−399
Nontriggered
Nontriggered
M16
94420
40262.0
−399
Nontriggered
Nontriggered
M17
94134.9
40328.1
−199
Nontriggered
NontriggeredThis paper offers the locations of two events. The first event is the mine production blasting at 14 : 52 : 59, on January 4, 2013. The known coordinate of source is (94489, 40441, −231). It is used for location accuracy comparisons of the new method and traditional methods. The second event is a mine earthquake measured accidentally. It is used for verifying the application of the partition method for the probability of source locations in screening and identifying potential safety hazard. The arrival times of these two events are listed in Table4. Vibration monitoring data of several production guns during 2013 are used for wave velocity calibration and the mean of one set of wave velocity is 5.1 m/ms. The probability distribution is logarithmic normal distribution log-N (1.545, 0.216) (passingK-S inspection, P value of 0.68). The location results of these two events are illustrated below.1st event: the coordinate obtained from traditional methods is (94534.2, 40432.1, −241.369), with deviation of 47.2 m. The coordinate obtained from the new method is (94505.8, 40448.5, −236.1), with deviation of 19.1 m. The location accuracy of the new method is greatly improved. The waveforms of event 1 are shown in Figure4.Figure 4
Waveforms of the 1st event.2nd event: I, II, and III divided by partition method based on the probability of seismic source in this paper are shown in Figure5. The geophysical prospecting results show that there are four abnormal areas nearby, which may be the area of rack mass fracture, respectively, with coordinates of P1 (95074.5, 40669, −216), P2 (95041, 40642.5, −243.5), P3 (95085.5, 40657, −283.5), and P4 (95057, 40666, −257). Among them, P3 is located in I area, P2 is located in II area, and P1 and P4 are located in III area. Therefore, it can be inferred that this event is most probable to be the rack fracture at P3 and less probable to be the rack fracture at P2. The probability of fracture at P1 and P4 can be screened.Figure 5
Space partition based on the probability of microseism source of the 2nd event.


## 7. Conclusions

A microseism source location algorithm considering the randomness of wave velocity and a space partition based on the probability of source are proposed. The new location method views the propagation wave velocity of seismic wave in different path as a random variable obeying certain probability distribution to further introduce the random probability distribution information of wave velocity into microseism source location algorithm successfully, which reduces the error caused by wave velocity uniformity hypothesis. The space partition makes space distribution statistics on Monte Carlo simulation results to build sphere model (applying to seismic source in sensor array) and cylinder model (applying to seismic source out of sensor array) to describe this distribution. The space is divided into I, II, and III areas, according to possibilities (60%, 95%, and less than 5%) of the location of microseism source.The example of the algorithm and field applied examples in Xinli minefield verified that the location accuracy of the new method in this paper is significantly superior to that of traditional methods, especially for the microseism source out of sensor array. Meanwhile, it also verified that space partition proposed in this paper can accurately predict the probability of microseism source in space. Combined with actual condition of mine, the new method also can identify and screen the potential safety hazards.



---

*Source: 805267-2014-09-28.xml*