Research Article | Open Access

Volume 2014 |Article ID 805267 | https://doi.org/10.1155/2014/805267

Furui Du, Nailian Hu, Yuling Xie, Guoqing Li, "An Undersea Mining Microseism Source Location Algorithm Considering Wave Velocity Probability Distribution", Mathematical Problems in Engineering, vol. 2014, Article ID 805267, 7 pages, 2014. https://doi.org/10.1155/2014/805267

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

Academic Editor: Gen Qi Xu
Revised31 Aug 2014
Accepted01 Sep 2014
Published28 Sep 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.

#### 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 [810], 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 to wave’s fastest propagation velocity among the seismic waves and its easily recognizable arrival time, the method of wave location is adopted in most cases. The method of calibration could help predetermine the velocity of wave, namely, , and is generally considered as a constant. Set the position of hypocenter as , origin time as , the coordinates of the stations which have received seismic signals as (, , ) , and the time of receiving seismic signals as . 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: In the above formula, the value of which represents the travel time from the hypocenter to Station could be illustrated in three-dimensional space as follows: where, .

Target function (1) contains four variables, the spatial coordinate parameters of hypocenter and the origin time . 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 The travel time in this equation can be expressed as Compared with (1), (3) contains one more unknown 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 and cumulative distribution function is noted as .

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 , the overall probability density function of wave velocity is , and unknown to-be-estimated parameters are . Therefore, the likelihood function can be noted as The estimated values of enable the likelihood function to reach the maximum value. The estimator of the maximum likelihood will be the solution of the following equation: The distributed fitting quality after estimated parameter needs to be inspected by K-S or , 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  . After various repetitions, sets of the most optimal solutions would be obtained. Using (7), the hypocenter location on average can be obtained: 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.

#### 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 clustered near one point; see Figure 2(a); as for the hypocenter out of the array of sensors, 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.

##### 4.1. Globular Model

The cumulative probability function of hypocenter distance in globular model can be defined as The probability function of is The value of represents the proportion of points in Monte Carlo simulation results whose space distances to the point are no more than . According to the law of large numbers, when the simulation times are enough, is close to the probability that the distance of the actual hypocenter to the calculated one is less than . When the value of is 60%, the corresponding means the probability that the distance of the hypocenter to the point is no more than is 60%, noted as . could be defined likewise. According to the distance between space point and , the space is divided into I, II, and III. I represents the area whose probability is less than , and the probability of the existence of hypocenter is the highest in this area; II represents the area whose probability is between and , representing the probability of the existence of hypocenter which is relatively higher in this area; III represents the location whose probability is more than , 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 in space can be written as below: Cumulative probability function of hypocenter distance in cylinder model can be defined as The in the equation is The value of represents the proportion of points in simulation results whose distances to the fitting straight-line (10) are no more than . Similar to globular model, , , 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 (258, 336, 580), (680, 290, 559), (190, 610, 380), and (789, 459, 280), respectively, while two are out of the array of sensors, and their coordinates are (308, 689, 1200) and (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 (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 wave recorded by sensors are shown in Table 1.

 Sensors Trigger times of seismic source/ms 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.961

The 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 , 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 , 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.

 Method Velocity Seismic source Seismic source m×ms−1 /m /m /m /m /m /m /m /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 source Seismic source m×ms−1 /m /m /m /m /m /m /m /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 source Seismic source m×ms−1 /m /m /m /m /m /m /m /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.857

Space partition based on the probability of seismic source is utilized to divide seismic sources (Table 3). It can be seen that all the seismic sources fall in I, except seismic source in II. This indicates that the partition method in this paper is reasonable. It is noteworthy that the location error of seismic source 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.

 Seismic source Partition model Real distance Partition area of seismic source m m m — Globular model 21.036 32.869 20.829 I Globular model 20.085 37.790 32.659 II Globular model 20.027 33.186 17.584 I Globular model 24.115 38.256 11.078 I Cylinder model 45.698 87.269 10.258 I 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 Figure 3.

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

 Number Sensor’s coordinate/m Arrival time of event/ms 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 Nontriggered

This 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 Table 4. 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- (1.545, 0.216) (passing K-S inspection, 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 Figure 4.

2nd event: I, II, and III divided by partition method based on the probability of seismic source in this paper are shown in Figure 5. 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.

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

#### Conflict of Interests

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

#### Acknowledgments

This research was supported by the Fundamental Research Funds for the Central Universities (no. FRF-SD-001A) and the Program for Changjiang Scholars and Innovative Research Team in University (IRT0950).

#### References

1. M. Zhai, J. Yang, and W. Liu, “Large clusters of gold deposits and large-scale metallogenesis in the Jiaodong Peninsula, Eastern China,” Science in China, Series D: Earth Sciences, vol. 44, no. 8, pp. 758–768, 2001. View at: Publisher Site | Google Scholar
2. K. Peng, X.-B. Li, C.-C. Wan, S.-Q. Peng, and G.-Y. Zhao, “Safe mining technology of undersea metal mine,” Transactions of Nonferrous Metals Society of China, vol. 22, no. 3, pp. 740–746, 2012. View at: Publisher Site | Google Scholar
3. Z. Liu, W. Dang, and X. He, “Undersea safety mining of the large gold deposit in Xinli District of Sanshandao Gold Mine,” International Journal of Minerals, Metallurgy and Materials, vol. 19, no. 7, pp. 574–583, 2012. View at: Publisher Site | Google Scholar
4. S.-L. Li, X.-G. Yin, W.-D. Zheng, and C. Trifu, “Research of multi-channel microseismic monitoring system and its application to Fankou lead-zinc mine,” Chinese Journal of Rock Mechanics and Engineering, vol. 24, no. 12, pp. 2048–2053, 2005. View at: Google Scholar
5. L. Tang, C. Yang, and C. Pan, “Optimization of microseismic monitoring network for large-scale deep well mining,” Chinese Journal of Rock Mechanics and Engineering, vol. 25, no. 10, pp. 2036–2042, 2006. View at: Google Scholar
6. F. Jiang, G. Ye, C. Wang, D. Zhang, and Y. Guan, “Application of high-precision microseismic monitoring technique to water inrush monitoring in coal mine,” Chinese Journal of Rock Mechanics and Engineering, vol. 27, no. 9, pp. 1932–1938, 2008. View at: Google Scholar
7. L. Geiger, “Probability method for the determination of earthquake epicenters from the arrival time only (translated from Geiger's 1910 German article ),” Bulletin of St. Louis University, vol. 8, no. 1, pp. 56–71, 1912. View at: Google Scholar
8. G. D. Nelson and J. E. Vidale, “Earthquake locations by 3-D finite-difference travel times,” Bulletin of the Seismological Society of America, vol. 80, no. 2, pp. 395–410, 1990. View at: Google Scholar
9. T. Theunissen, Y. Font, S. Lallemand, and S. Gautier, “Improvements of the maximum intersection method for 3D absolute earthquake locations,” Bulletin of the Seismological Society of America, vol. 102, no. 4, pp. 1764–1785, 2012. View at: Publisher Site | Google Scholar
10. N. Li, E. Wang, M. Ge, and Z. Sun, “A nonlinear microseismic source location method based on Simplex method and its residual analysis,” Arabian Journal of Geosciences, pp. 1–10, 2013. View at: Publisher Site | Google Scholar
11. R. S. Crosson, “Crustal structure modeling of earthquake data: 1. Simultaneous least squares estimation of hypocenter and velocity parameters,” Journal of Geophysical Research, vol. 81, no. 17, pp. 3036–3046, 1976. View at: Publisher Site | Google Scholar
12. T. Garza, C. Lomnitz, and C. R. de Velasco, “An interactive epicenter location procedure for the RESMAC seismic array: I,” Bulletin of the Seismological Society of America, vol. 67, no. 6, pp. 1577–1586, 1977. View at: Google Scholar
13. C. H. Thurber, “Nonlinear earthquake location: theory and examples,” Bulletin of the Seismological Society of America, vol. 75, no. 3, pp. 779–790, 1985. View at: Google Scholar
14. L. Dong, X. Li, L. Tang, and F. Gong, “Mathematical functions and parameters for microseismic source location without pre-measuring speed,” Chinese Journal of Rock Mechanics and Engineering, vol. 30, no. 10, pp. 2058–2067, 2011. View at: Google Scholar
15. X. B. Li and L. J. Dong, “Comparison of two methods in acoustic emission source location using four sensors without measuring sonic speed,” Sensor Letters, vol. 9, no. 5, pp. 2025–2029, 2011. View at: Publisher Site | Google Scholar

Copyright © 2014 Furui Du et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Related articles

No related content is available yet for this article.