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 [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 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.
(a)
(b)
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.
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.
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.
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.
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).