# A Multistep Method for Automatic Determination and Optimization of Microseismic P-Phase Arrival Times in a Coal Mine

Authors: Quanjie Zhu; Xiaohui Liu; Xiaoyun Liu; Haitao Chai; Jianjun Shi
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

We present an algorithm for P-phase arrival time automatic determination and optimization for microseismic events of a coal mine in steps. For P-phase detection, we apply four modules: single denoising (band-filter), fast searching (STA/LTA and peak value), accurately picking (AIC), and optimization selection (time difference of arrival time and envelope circle judgment). This method is applied for a precise determination of P-phase arrival time; in addition, an automatic quality estimate is obtained to optimize and select the effective stations and used for source location. Through the analysis and validation of seismic data in a mine located in Hebei province, the result shows that this method can pick up the first arrival time of MS signal fast and accurately and automatically judge and remove the error station. Comparing automatic picks with the manual method in terms of speed and accuracy, computation cost (0.828 s) is far less than manual method (15 s per set), precision of manual and automatic represented by picking error is up to 71.7% (between −5 ms and 5 ms, 190 sets), 83.4% (between −10 ms and 10 ms, 221 sets), and 90.57% (between −20 ms and 20 ms, 240 sets). This study shows that the method presented in this paper is feasible to precise automatic P-phase determination and is a basis on automatic and precision location.


---

## Body

## 1. Introduction

For a long time in China, deep marching into the earth is a strategic and technological problem that has been solved. The forecast of rockburst disaster derived from deep mining and geotechnical engineering is one of the urgent technological challenges that need to be realized urgently. Such disasters are characterized with strong suddenness, short process, and strong destructive force. It will lead to a significant increase in occurrence frequency and intensity of deep disasters, which seriously threaten the production safety of project sites. Microseismic (MS) is an effective technique for rockburst disasters forecast. How to improve the rapid and accurate analysis abilities (such as automatic positioning, unattended operations) of microseismic technology and the accurate pickups at first arrival (PK-FA) signal is a very important precondition. Consequently, development of an algorithm for PK-FA is the key for microseismic data rapid processing and advantageous to the disaster emergency treatment.Microseismic monitoring in mining engineering is different to the earthquake, and it is difficult to distinguish P and S waves using a single component sensor. Also, the pickup of characteristic P wave is vital to monitor mine/rock microseisms. However, attributes to effects from operating conditions and environment, the microseismic signal component is quite complex and its frequency even tends to high level. Different from the simple mechanical vibrations, the actual microseismic signal is not in single and periodic forms but accompanied by external interferences. The interferences mainly originate from artificial operations, transmission of microseisms system, heterogeneity of transmission medium, etc. Therefore, to meet the requirements of timeliness and position accuracy, the vast amounts and various kinds of triggered data from engineering site will bring great difficulties on field processing and analysis. As reported in the literature, the theoretical study method (TSM) on ideal or high signal-to-noise signals have shown satisfactory results. However, there are certain limitations for TSM whose identification precision is too poor to process and analyse the low SNR signals being in mine/rock microseisms.In present, there are mainly four methods for the vibrational signal PK-FA: time domain method (TDM), frequency domain method (FDM), time-frequency method (TFM), and comprehensive analysis method (CAM). Among these methods, TDM could be roughly categorized in Short- and Long-Time Average (STA/LTA) ratios, digital image (DI) [1], cross-correlation [2], fractal theory [3, 4], and autoregression model. The STA/LTA algorithm [5] proposed by Allen is a classical TDM that takes advantages of simple in model and high running speed. However, the accuracy of STA/LTA algorithm was limited by the size of time-window. By modifying for proposing new feature function of the STA/LTA algorithm, many methods for vibrational signal pickups were developed [6], such as Modified Energy Ratio (MER) [7], Modified Coppens’ Method (MCM) [8], etc. These modified methods have their advantages on pickups of characteristic function’s composition, commonly using PAI-S/K and AIC algorithms. In the PAI-S/K algorithm [9, 10], the PK-FA was based on statistical analysis of waveform data to obtain the skew and kurtosis of signal. But the PAI-S/K algorithm was limited by the signal length during the meshing of geometrical configuration. At present, Akaike Information Criterion (AIC) [11] is a widely used statistical criterion for the pickups, being in features on high speed and precision.The algorithms based on the frequency domain (FD) and time domain (TD) have been proposed [12] for the PK-FA of seismic waves by means of techniques of time-frequency domain energy ratio [13], time-frequency analysis (TFA) [14], neural network [15], etc. The TFA realized the PK-FA by means of Fourier transform, Wavelet transform, etc. to map the one-dimensional seismic signal to two-dimensional time-frequency plane. The time and frequency characteristics of signal is highlighted during mapping. In addition to the using of single algorithm, the PK-FA can be realized via comprehensive analysis method (CAM) which combined the above methods [16, 17]. As reported in literatures, the CAM comprehensively took signal frequency, energy, and vibration track into account, showing a well results on PK-FA in engineering site test [18, 19].Attributes to the heterogeneity in transmission medium, the positioning accuracy is often affected by the station distribution. Since the characteristics of surrounding rocks are quite different from the platform, the real time is different from that of other stations at initial time (nonlinear). If the station is selected to participate in positioning calculation, a large error will be caused. Our previous study presented the concepts of waveform correlation and time-difference method, which showed good results but with low efficiency on positioning calculations of channels screen via artificial method involved [20].Based on the aforementioned background, a model of multistage detection/optimization for PK-FA of P wave is developed, aiming to process the microseismic data in automation and high precision. The modules of waveform denoising, on-time fast searching, accurate picking, and efficient channel judgments are proposed. Our model is expected to lay a foundation for the subsequent high-precision positioning calculation of PK-FA.


## 2. Common Methods and Their Principles

### 2.1. STA/LTA Method

It could be concluded that the more complex the algorithm, the larger computations and longer time. In present, PK-FA of microseismic waveform on earthquakes, petroleum, rocks/mines mainly by means of STA/LTA, MER, MCM, PAI-K, S/L-Kurt, AIC and JER, etc. Akram statistically studied the pickups of field data (1000 groups) by using these methods [21]. From the detailed results of comparisons included time-consuming ratio and error size of the algorithms, it indicated that AIC method is in high picking precision and stable calculation precision on PT-FA. And the STA/LTA [22] as well as its extended algorithms has been widely utilized for the PT-FA of seismic wave. The algorithm of STA/LTA is expressed as(1)Ri=STAiLTAi≥λ.The STA and LTA is, respectively, calculated as(2)STAi=∑j=i−ns+1nsCFjns,(3)LTAi=∑j=i−nl+1nlCFjnl,where λ is the trigger threshold of feature function; nl and ns are the size of long time-window (LTW) and short time-window (STW), respectively; j refers to jth moment sampling point in [i − nl + 1, i] range; and CF(j) denotes the eigenvalue (such as the amplitude, amplitude squared, mean variance, etc.) of microseismic signal at jth moment.In the STA/LTA method, the STW was considered to be including LTW. The eigenvalue of LTA and STA is then calculated by the sliding established LTW (equation (2)) and STW (equation (3)), which represents the mean value of background noise and effective microseismic signal, respectively. Once the time-window develops from the beginning to arrival time, STA/LTA ratio value R will manifest a significant increase. The field engineering application revealed that the STA/LTA algorithm was characterized as simple in model and fast in calculation speed. It could meet the demands in field monitoring system for rapid detection of microseismic events. However, its precision on PT-FA is limited by the over-dependent time-window size and threshold settings. Therefore, the fast and efficient features of STA/LTA can be used to quickly locate the general position of PT-FA.


### 2.2. AIC Method

Based on the Akaike Information Criterion (AIC), AR-AIC algorithm divides the nonstationary signal into several fixed length waveform segments. The divided segments are then imported into a regression (AR) model whose order and coefficient are further solved. Thus, the calculation process of AR-AIC aims to find the local minimum value of AIC function in microseismic signal, expressed as(4)AICi=i−klogσ1,max2+N−i−klogσ2,max2+C,where C is a constant; N the size of microseismic data; k the order of AR process; i the boundary points between two local statistical periods; σ2,max2 and σ1,max2 the errors of two local statistical periods.Maeda [23] developed a modified AR-AIC model (equation (5) by eliminating the autoregressive coefficient of AR model.(5)AICi=ilogvarx1,i+N−i−klogvarxi+1,N+C,where x(i) is the amplitude of a set of waveform data. Figure 1 shows the microseismic wave pickups of PT-FA by the modified AR-AIC model. It could be suggested that the advantages of the modified AR-AIC model lie in fast calculation speed and accurate pickups. But it has the defect of excessive dependence of pickups accuracy on a window selection.Figure 1
Result of the first arrival picking using the AIC method. (a) A typical MS signal; (b) scaling figure from red part of (a) and two vertical lines show the first arrival time (onset time) and terminal time (offset time); (c) local eigenvalue of AIC method; and (d) global eigenvalue in the whole period.


## 2.1. STA/LTA Method

It could be concluded that the more complex the algorithm, the larger computations and longer time. In present, PK-FA of microseismic waveform on earthquakes, petroleum, rocks/mines mainly by means of STA/LTA, MER, MCM, PAI-K, S/L-Kurt, AIC and JER, etc. Akram statistically studied the pickups of field data (1000 groups) by using these methods [21]. From the detailed results of comparisons included time-consuming ratio and error size of the algorithms, it indicated that AIC method is in high picking precision and stable calculation precision on PT-FA. And the STA/LTA [22] as well as its extended algorithms has been widely utilized for the PT-FA of seismic wave. The algorithm of STA/LTA is expressed as(1)Ri=STAiLTAi≥λ.The STA and LTA is, respectively, calculated as(2)STAi=∑j=i−ns+1nsCFjns,(3)LTAi=∑j=i−nl+1nlCFjnl,where λ is the trigger threshold of feature function; nl and ns are the size of long time-window (LTW) and short time-window (STW), respectively; j refers to jth moment sampling point in [i − nl + 1, i] range; and CF(j) denotes the eigenvalue (such as the amplitude, amplitude squared, mean variance, etc.) of microseismic signal at jth moment.In the STA/LTA method, the STW was considered to be including LTW. The eigenvalue of LTA and STA is then calculated by the sliding established LTW (equation (2)) and STW (equation (3)), which represents the mean value of background noise and effective microseismic signal, respectively. Once the time-window develops from the beginning to arrival time, STA/LTA ratio value R will manifest a significant increase. The field engineering application revealed that the STA/LTA algorithm was characterized as simple in model and fast in calculation speed. It could meet the demands in field monitoring system for rapid detection of microseismic events. However, its precision on PT-FA is limited by the over-dependent time-window size and threshold settings. Therefore, the fast and efficient features of STA/LTA can be used to quickly locate the general position of PT-FA.


## 2.2. AIC Method

Based on the Akaike Information Criterion (AIC), AR-AIC algorithm divides the nonstationary signal into several fixed length waveform segments. The divided segments are then imported into a regression (AR) model whose order and coefficient are further solved. Thus, the calculation process of AR-AIC aims to find the local minimum value of AIC function in microseismic signal, expressed as(4)AICi=i−klogσ1,max2+N−i−klogσ2,max2+C,where C is a constant; N the size of microseismic data; k the order of AR process; i the boundary points between two local statistical periods; σ2,max2 and σ1,max2 the errors of two local statistical periods.Maeda [23] developed a modified AR-AIC model (equation (5) by eliminating the autoregressive coefficient of AR model.(5)AICi=ilogvarx1,i+N−i−klogvarxi+1,N+C,where x(i) is the amplitude of a set of waveform data. Figure 1 shows the microseismic wave pickups of PT-FA by the modified AR-AIC model. It could be suggested that the advantages of the modified AR-AIC model lie in fast calculation speed and accurate pickups. But it has the defect of excessive dependence of pickups accuracy on a window selection.Figure 1
Result of the first arrival picking using the AIC method. (a) A typical MS signal; (b) scaling figure from red part of (a) and two vertical lines show the first arrival time (onset time) and terminal time (offset time); (c) local eigenvalue of AIC method; and (d) global eigenvalue in the whole period.


## 3. Principles and Process of Step Picking Method

The drawbacks of common methods on PK-FA lie in time consuming and inaccuracy. It is unable to conduct the analysis and treatment timely, which inevitably affects the emergency response of outburst on the on-site microseismic data in massive scale. Furthermore, common methods can only extract the initial value of the waveform without judgments on the accuracies and errors, which affect the locating calculation being in the time-interval method.


### 3.1. Principles and Process

To realize the fast and accurate pickups, an operation mechanism for quick judgment of arrival time is developed in Figure2. In brief, such mechanism aims to realize the arrival time picking, includes two major sections of accurate picking (AP) and effective judgment/optimization (EJ&O). The whole picking process is carried out step-by-step in four levels: denoising (DNS), fast searching (FSC), AP, and EJ&O, which is, respectively, illustrated as follows:(1)
DNS step: Considering the relatively concentrated feature of frequency range of mine microseismic signal, the original signal is denoised by using Butterworth band-pass filter (10∼200 Hz) [24] which is coded in MATLAB®.(2)
FSC step: The peak point in signals is searched fast via the peak judgment method. The position of initial arrival is then determined by the STA/LTA method by moving 500 sampling points around the peak point.(3)
AP step: The picking of initial arrival is based on the AIC algorithm, which is a conventional method with rapid identification, short time consumption, high precision, and wide application advantages.(4)
EJ&O step: The EJ&O channels can be accomplished through three steps: cluster analysis for reference channel, time-difference judgment eliminates invalid channels, and envelope curve judgment on error channel, respectively.Figure 2
Flow chart of arrival time picking and effective judgment.


### 3.2. Accurately Picking Up

The accurate pickups of initial arrival are the hierarchical testing and optimization core of the whole process which consist two steps of FSC and AP, as shown in Figure2. In the FSC step, the nonoverlapping time-window of 500 ms is firstly used to search within the entire signal (5000 ms length). The last sampling point of the first time-window is treated as the starting point of the second searching. The sum and mean of absolute amplitude value in time-window is calculated to obtain a series of eigenvalues during the whole searching process. Once the signal has been processed, the maximum position of eigenvalue is calculated, which is considered as the center and extends to both sides of time-window, respectively. Then, a larger extensional time-window with 1500 ms length is obtained. As illustrated in Figures 1(a) and 1(b), the red background ranges refer to the time-window position from initial to final arrivals, in which the dash line being the curve of searched eigenvalues. After the FSC, the eigenvalues are used to determine the initial arrival position of microseismic waveform via AIC and MER methods in the AP step. Results of the AP step are shown as the red position of Figure 1(c) that refers to picking up between first arrival and final arrival.


### 3.3. Validity Judgment and Optimization

The EJ&O step of channels includes following three parts:(1)
Cluster analysis. In this part, the difference of initial to final arrivals (DIF) is analysis via spatial cluster (k = 2). The cluster with more individuals is treated as the main cluster, in which point A with the smallest distance from the center of concentration is established as optimal reference point.(2)
DIF judgment method. The reference channel A is selected as a reference channel to calculate the spatial distance and actual time difference between channels, respectively. And the theoretical maximum DIFs of any channel to channel A can be calculated by spatial distance and measured wave velocity. By comparing the theoretical and actual time differences, the judgment threshold is selected and the channel with large time error can be eliminated. The judgment denotes the preliminary screening result.(3)
Envelope curve judgment on error channel. Based on the judgment result of previous step, the first trigger channelC0 (with smallest arrival) is selected as reference channel, which is a standard point to calculate the actual DIF to Ci channels (i is the number of channels other than 0). The radius value Ri is then calculated from the wave velocity and drawn a circle, i.e., the spatial position of microseismic events undoubtedly locates in the circle, while the events outside the circle can be preliminarily located to obtain the distance Ui between the microseismic source and stations. When the sum (Ui) is at its lowest, it is the optimal location, and then reverse-pushRi (the subsequent research in this paper will be introduced). The optimum location being in minimum sum (Ui) is used to reversely derive Ri.The as-location of microseismic source is used to screen the channels that involved in the calculations. By using the simplex method, the located calculations of the model is express asO’(xf, yf, zf, tf). And propagation time between microseismic source and stations can be expressed as(6)ΔTi=xi−xf2+yi−yf2+zi−zf2V=ti−tf.Calculation of the optimal microseismic source actually is to minimize the objective function of(7)λ=∑i=1Nti−tf=∑i=1NRiV,where N is the stations number and Ri distance (radius) between microseismic source and stations. According to this judgment method, the smaller error of time between microseismic source and station indicates the higher reliability of positioning calculation. On the contrary, the larger error means lower reliability.Theoretically, under the condition of homogeneous media, center of the platform is suggested as the circle center. The distance between platform stations to microseismic source is drawn as the radius, whose circular curves (stations) are converged to microseismic source, as illustrated in Figure3 (blue curve denotes the monitoring station, and red is microseismic source).Figure 3
Sketch map of sensors (circles) and source location.However, the heterogeneity of propagation medium and variation of wave velocity indicate the assumptions exist only in theory. Even so, the validity of PT-FA judgment can be determined by considering the above factors, giving a certain surplus ΔR amount and a fictional Δ radius R + R.


## 3.1. Principles and Process

To realize the fast and accurate pickups, an operation mechanism for quick judgment of arrival time is developed in Figure2. In brief, such mechanism aims to realize the arrival time picking, includes two major sections of accurate picking (AP) and effective judgment/optimization (EJ&O). The whole picking process is carried out step-by-step in four levels: denoising (DNS), fast searching (FSC), AP, and EJ&O, which is, respectively, illustrated as follows:(1)
DNS step: Considering the relatively concentrated feature of frequency range of mine microseismic signal, the original signal is denoised by using Butterworth band-pass filter (10∼200 Hz) [24] which is coded in MATLAB®.(2)
FSC step: The peak point in signals is searched fast via the peak judgment method. The position of initial arrival is then determined by the STA/LTA method by moving 500 sampling points around the peak point.(3)
AP step: The picking of initial arrival is based on the AIC algorithm, which is a conventional method with rapid identification, short time consumption, high precision, and wide application advantages.(4)
EJ&O step: The EJ&O channels can be accomplished through three steps: cluster analysis for reference channel, time-difference judgment eliminates invalid channels, and envelope curve judgment on error channel, respectively.Figure 2
Flow chart of arrival time picking and effective judgment.


## 3.2. Accurately Picking Up

The accurate pickups of initial arrival are the hierarchical testing and optimization core of the whole process which consist two steps of FSC and AP, as shown in Figure2. In the FSC step, the nonoverlapping time-window of 500 ms is firstly used to search within the entire signal (5000 ms length). The last sampling point of the first time-window is treated as the starting point of the second searching. The sum and mean of absolute amplitude value in time-window is calculated to obtain a series of eigenvalues during the whole searching process. Once the signal has been processed, the maximum position of eigenvalue is calculated, which is considered as the center and extends to both sides of time-window, respectively. Then, a larger extensional time-window with 1500 ms length is obtained. As illustrated in Figures 1(a) and 1(b), the red background ranges refer to the time-window position from initial to final arrivals, in which the dash line being the curve of searched eigenvalues. After the FSC, the eigenvalues are used to determine the initial arrival position of microseismic waveform via AIC and MER methods in the AP step. Results of the AP step are shown as the red position of Figure 1(c) that refers to picking up between first arrival and final arrival.


## 3.3. Validity Judgment and Optimization

The EJ&O step of channels includes following three parts:(1)
Cluster analysis. In this part, the difference of initial to final arrivals (DIF) is analysis via spatial cluster (k = 2). The cluster with more individuals is treated as the main cluster, in which point A with the smallest distance from the center of concentration is established as optimal reference point.(2)
DIF judgment method. The reference channel A is selected as a reference channel to calculate the spatial distance and actual time difference between channels, respectively. And the theoretical maximum DIFs of any channel to channel A can be calculated by spatial distance and measured wave velocity. By comparing the theoretical and actual time differences, the judgment threshold is selected and the channel with large time error can be eliminated. The judgment denotes the preliminary screening result.(3)
Envelope curve judgment on error channel. Based on the judgment result of previous step, the first trigger channelC0 (with smallest arrival) is selected as reference channel, which is a standard point to calculate the actual DIF to Ci channels (i is the number of channels other than 0). The radius value Ri is then calculated from the wave velocity and drawn a circle, i.e., the spatial position of microseismic events undoubtedly locates in the circle, while the events outside the circle can be preliminarily located to obtain the distance Ui between the microseismic source and stations. When the sum (Ui) is at its lowest, it is the optimal location, and then reverse-pushRi (the subsequent research in this paper will be introduced). The optimum location being in minimum sum (Ui) is used to reversely derive Ri.The as-location of microseismic source is used to screen the channels that involved in the calculations. By using the simplex method, the located calculations of the model is express asO’(xf, yf, zf, tf). And propagation time between microseismic source and stations can be expressed as(6)ΔTi=xi−xf2+yi−yf2+zi−zf2V=ti−tf.Calculation of the optimal microseismic source actually is to minimize the objective function of(7)λ=∑i=1Nti−tf=∑i=1NRiV,where N is the stations number and Ri distance (radius) between microseismic source and stations. According to this judgment method, the smaller error of time between microseismic source and station indicates the higher reliability of positioning calculation. On the contrary, the larger error means lower reliability.Theoretically, under the condition of homogeneous media, center of the platform is suggested as the circle center. The distance between platform stations to microseismic source is drawn as the radius, whose circular curves (stations) are converged to microseismic source, as illustrated in Figure3 (blue curve denotes the monitoring station, and red is microseismic source).Figure 3
Sketch map of sensors (circles) and source location.However, the heterogeneity of propagation medium and variation of wave velocity indicate the assumptions exist only in theory. Even so, the validity of PT-FA judgment can be determined by considering the above factors, giving a certain surplus ΔR amount and a fictional Δ radius R + R.


## 4. Evaluation of Experimental Data

In order to verify the rationality and validity of this method, we used the theoretical method and experimental data to verify above method. The process of this verification includes the synthetic signals of microseismic generation, first arrival time picking, and result evaluation.


### 4.1. Generation of Synthetic Signals

The theory method of synthetic signal generation presented in Miao’s article [25] was applied in this paper. On his view, different phases of seismic waves are generated by an internal point source and propagated in the homogeneous rock mass; considering the unevenness of the propagation medium, the influence of noise is added. Based on the wave mechanics and other relevant theories, the wave equation of mining microseismic can be approximately considered as(8)∂2rfr,t∂r2−1v2∂2rfr,t∂t2=0,fr,t=1rf10,t−rv,where r is the propagation distance of seismic; t propagation time of seismic wave in rock mass; and f(r, t) seismic wave which propagated t time and r distance of the source.Based on the real engineering of mining, the model above is transformed to geodetic coordinates. Assume that the location of seismic source isR0(x0, y0, z0), the seismic wave Ri(xi, yi, zi) in position Ri(xi, yi, zi) which is far from source time ts and distance v∗Δt can be considered as(9)ti=Δt+t,ri=R−R0=xi−x02+yi−y02+zi−z02,fri,ti=gxi,yi,zi,ti.By means of the model above and the microseismic monitoring network of a coal mine in Hebei province, a set of seismic data is generated in this paper (as shown in Figure4). Figure 4(a) shows the seismic model and survey diagram (12 receivers corresponds to 12-channel waveforms). The red dot is the location of real seismic source, and the blue circle shows the monitoring receiver. Figure 4(b) shows the 12-channel received records of the monitoring network.Figure 4
Model and survey diagram and receiving records of monitoring network. (a) Model and survey diagram. (b) Receiving records of monitoring network.
(a)
(b)From this, it can be seen that difference in distance or time corresponds different record (waveform). Also, it has a different first arrival time. Based on the above theoretical model and experimental data, the arrival time pickers of all 12 receivers has been obtain using the model provided in Section3.


### 4.2. Data Preprocess

The first arrival time pickers of receiving records obtained by the present method are shown in Figure5. Blue curves in this figure show the eigenvalues of arrival time calculated and red lines show the onset time and the end time. We can see that the eigenvalues curve is so smooth so that we can pick the key point clearly and accurately, although noise has been added to the theoretical data.Figure 5
Eigenvalues (onset time and end time) of receiving record.


### 4.3. Evaluation Result

Theoretically, the ratio of distance from source and travel time of seismic in isotropic media is a fixed value, that is that propagation velocity of seismic is a constant. The velocity of seismic is 3000 m/s in our theoretical model. In order to illustrate the effectiveness of the proposed method, picking results of all receivers are shown in Figure6. In this figure, blue curve is the waveform of seismic, red circle is the first arrival time of receiving record, and green curve is the line segments of all first arrival times. It can be seen form the figure that green curve is approximately a straight line; this indicates that the slope and the velocity of seismic is consistent. Accordingly, theoretical analysis agrees well with experiment results.Figure 6
Picking results of all receivers showed in the time-distance figure.In order to describe the similarity of theoretical value and calculated value more intuitively, a contrastive analysis on true value and calculated value of seismic first arrival time is studied and shown in Figure7. From the image, we can see that red circle represents the onset time (calculated value) of seismic waveform and red curve represents the fitted curve of calculated value. Blue curve represents the line segments of all first arrival times of theoretical values. Obviously, the blue curve and red curve coincide with each other approximately. This indicates that the method presented in Section 3 is effectively applied to first time picking of seismic waveforms and has obvious characteristics of rationality and validity.Figure 7
Comparison between the calculated value and true (theoretical) value.Theoretical data cannot depict the characteristics exactly and wholly; therefore, we will continue to study and verify the method provided in this paper using field monitoring data.


## 4.1. Generation of Synthetic Signals

The theory method of synthetic signal generation presented in Miao’s article [25] was applied in this paper. On his view, different phases of seismic waves are generated by an internal point source and propagated in the homogeneous rock mass; considering the unevenness of the propagation medium, the influence of noise is added. Based on the wave mechanics and other relevant theories, the wave equation of mining microseismic can be approximately considered as(8)∂2rfr,t∂r2−1v2∂2rfr,t∂t2=0,fr,t=1rf10,t−rv,where r is the propagation distance of seismic; t propagation time of seismic wave in rock mass; and f(r, t) seismic wave which propagated t time and r distance of the source.Based on the real engineering of mining, the model above is transformed to geodetic coordinates. Assume that the location of seismic source isR0(x0, y0, z0), the seismic wave Ri(xi, yi, zi) in position Ri(xi, yi, zi) which is far from source time ts and distance v∗Δt can be considered as(9)ti=Δt+t,ri=R−R0=xi−x02+yi−y02+zi−z02,fri,ti=gxi,yi,zi,ti.By means of the model above and the microseismic monitoring network of a coal mine in Hebei province, a set of seismic data is generated in this paper (as shown in Figure4). Figure 4(a) shows the seismic model and survey diagram (12 receivers corresponds to 12-channel waveforms). The red dot is the location of real seismic source, and the blue circle shows the monitoring receiver. Figure 4(b) shows the 12-channel received records of the monitoring network.Figure 4
Model and survey diagram and receiving records of monitoring network. (a) Model and survey diagram. (b) Receiving records of monitoring network.
(a)
(b)From this, it can be seen that difference in distance or time corresponds different record (waveform). Also, it has a different first arrival time. Based on the above theoretical model and experimental data, the arrival time pickers of all 12 receivers has been obtain using the model provided in Section3.


## 4.2. Data Preprocess

The first arrival time pickers of receiving records obtained by the present method are shown in Figure5. Blue curves in this figure show the eigenvalues of arrival time calculated and red lines show the onset time and the end time. We can see that the eigenvalues curve is so smooth so that we can pick the key point clearly and accurately, although noise has been added to the theoretical data.Figure 5
Eigenvalues (onset time and end time) of receiving record.


## 4.3. Evaluation Result

Theoretically, the ratio of distance from source and travel time of seismic in isotropic media is a fixed value, that is that propagation velocity of seismic is a constant. The velocity of seismic is 3000 m/s in our theoretical model. In order to illustrate the effectiveness of the proposed method, picking results of all receivers are shown in Figure6. In this figure, blue curve is the waveform of seismic, red circle is the first arrival time of receiving record, and green curve is the line segments of all first arrival times. It can be seen form the figure that green curve is approximately a straight line; this indicates that the slope and the velocity of seismic is consistent. Accordingly, theoretical analysis agrees well with experiment results.Figure 6
Picking results of all receivers showed in the time-distance figure.In order to describe the similarity of theoretical value and calculated value more intuitively, a contrastive analysis on true value and calculated value of seismic first arrival time is studied and shown in Figure7. From the image, we can see that red circle represents the onset time (calculated value) of seismic waveform and red curve represents the fitted curve of calculated value. Blue curve represents the line segments of all first arrival times of theoretical values. Obviously, the blue curve and red curve coincide with each other approximately. This indicates that the method presented in Section 3 is effectively applied to first time picking of seismic waveforms and has obvious characteristics of rationality and validity.Figure 7
Comparison between the calculated value and true (theoretical) value.Theoretical data cannot depict the characteristics exactly and wholly; therefore, we will continue to study and verify the method provided in this paper using field monitoring data.


## 5. Example Verification

In order to verify the aforementioned method, we take a previous detected microseismic event (2011-10-07T13:42:17) as an example. The experimental site (parameters listed in Table1) is located in a workface of coal mine in Hebei province. The aim of this algorithm is to monitor the microseismic event in real-time.Table 1
Station layout parameters.
Sensor coordinates (m)
Station number (#)
X axial
Y axial
Z axial
Station location
1
438.82
389.3
−191.49
Haulage roadway
2
446.95
360.04
−184.306
Haulage roadway
3
454.68
333.56
−178.669
Haulage roadway
4
463.56
303.22
−174.108
Haulage roadway
5
471.99
275.47
−166.204
Haulage roadway
6
481.08
247.3
−162.075
Haulage roadway
7
572.29
387.38
−186.378
Ventilation roadway
8
598.85
301.87
−172.128
Ventilation roadway
9
590.08
329.94
−175.486
Ventilation roadway
10
598.85
301.87
−172.128
Ventilation roadway
11
607.09
272.79
−166.429
Ventilation roadway
12
616.86
244.66
−160.446
Ventilation roadway


### 5.1. Bandfilter: Denoising Module

The statistical analysis of field data shows that the frequency range of microseismic signal is affected by many factors, such as transmission distance, surrounding rock lithology, and microseismic source itself. These factors not only constrain the PT-FA, but also affect the accuracy to a great extent. Therefore, it is very necessary to denoise the microseismic data before the PT-FA beginning. It is famous that Butterworth band-pass filter is widely utilized in denoise of seismic signal. Taking the above factors in to consideration, a Bandfilter (code in MATLAB®) is used to process the microseismic signal by setting the effective frequency distribution range of 10∼200 Hz. Figure8 shows the band-pass denoise results of typical microseismic signal. It could be found that the main components and energy of signal distributes in range of 10∼200 Hz. The original signal shows a sparse distribution in the high-frequency region. However, but the main components in the high-frequency region is eliminated. According to the method described in the literature [16], the calculated SNR = 35.6126 db is higher and indicates that the denoising signal SNR has been greatly improved.Figure 8
Waveform and spectrum figures: signal and its frequency spectrum before (a) and after (b) denoising.
(a)
(b)Two typical microseismic signals were processed by the Bandfilter module, and the results of those can be seen in Figure9. As can be seen from the figures, this method can effectively remove the interference components and improve the signal-to-noise ratio (SNR). This is the basis of first arrival time pick automatically.Figure 9
Validation of typical signals: (a) denoising before (left) and after (right) signal 3#; (b) denoising before (left) and after (right) signal 7#.
(a)
(b)


### 5.2. First Arrival Time Picking Module

According to the mentioned calculation process in Section2, the PT-FA of studied microseismic event was carried out and is shown in Figure 10, which shows the “time-amplitude” curve within 12 channels. The red line and Stations 1∼12 denote the PT-FA results and 12 sensors, respectively. It can be obviously found that our algorithm can pick up the signal of Stations 1, 2, 3, 5, 6, 7, 8, 9, and 11. There are no effective waveform fluctuations found in channels 10 and 12 whose calculated initial arrival time is 1 s. Regarding channel 4 with lower SNR value, its signal noise is relatively low with a pickup time of 0.003 s. These results indicate that our algorithm is effective to realize the PT-FA of microseismic event under certain SNR.Figure 10
Waveforms and picked first arrival times acquired by the MS system with 12 sensors (vertical lines show the first arrival times).From Figure10 and Table 2, it could be concluded that the higher the SNR, the more accurate the picking time. On the contrary, it is difficult to accurately pick up or judge the SNR or abnormal signal of initial arrival. In Figure 10, channels of 4, 6, 10, and 12 can be directly observed by the naked eye. But if the algorithm returns values without judgment and optimization, it will seriously affect the accuracy of subsequent positioning. After the pickup of initial arrival, all arrival times are substituted in the positioning model for seismic source calculations. It could be found in Figure 10 that the judgments of channel 4, 10, and 12 are invalid according to the algorithm. The judgment process can be automatically realized through EJ&O step in Section 2.Table 2
First arrival time of every station (sensor).
Channel number
Arrival time (s)
1
0.507
2
0.505
3
0.519
4
0.003
5
0.530
6
0.490
7
0.542
8
0.471
9
0.527
10
1
11
0.529
12
1


### 5.3. Efficient Channel Optimization

The EJ&O of channels is to ensure the accuracy of positioning calculation. The initial arrival values of each time period are screened out, and invalid or error values are eliminated.


#### 5.3.1. Choosing Reference Object

Since the sensors are located in a same working face and arranged at 30–50 m intervals, the time interval between triggering sensors is not large sensor as for the microseismic events happen in network layout. The channels with larger errors can be eliminated by using clustering calculations (k-means clustering algorithm,k = 2) on these trigger moments. The selection of reference channels is realized by means of clustering according to Table 3 and Figure 11, in which 12 initial arrivals are clustered into two types and “∗” signs denotes clustering center (Figure 8).Table 3
Result of first arrival times.
Channel number
Distance of station and cluster A
Cluster number
1
0.0020
1
2
0.0018
1
3
0.0032
1
4
0.2110
1
5
0.0064
1
6
7.67 × 10–4
1
7
0.0064
1
8
7.57 × 10–5
1
9
0.0042
1
10
0.2891
2
11
0.0044
1
12
0.2891
2Figure 11
Clustering of arrival times and choice of refer station. (a) Clustering results using k-means; (b) distance (normalized value) of current station and cluster A.
(a)
(b)It can be observed from Figure11 that all the channels are clustered into one type except Channels 10 and 12. As the elements in cluster A are much more than those in cluster B, class A with larger element number are selected to continue the study. Then, the closest point (station) to the cluster center (mean value at that time in the cluster) in cluster is taken as the reference, and class B is excluded. That is to say, once the number of elements in two clusters equals each other, the centers (two clusters) can be selected as the cluster center in order to find the distance (Ri). Figure 11(b) gives the distance between each element to its cluster center (cluster A), where 8# station having a minimum distance of 7.57E-4 can be used as reference station.


#### 5.3.2. Culling of Dead Channel

After the treatment of first step, it is necessary to solve the problem of incorrect initial arrival time value or large error. Figure12 compares the actual arrival time difference, ΔTA (actual value) between 8# and other stations. The theoretical ΔTT of channels 4, 10, and 12 are lower than the actual ΔTT. By assuming λ = ΔTA/ΔTT, the channels with larger errors can be judged as λ > 10.Figure 12
Comparison histogram of theoretical maximum value and actual value.From the calculated results of relevant parameters of arrival time difference in Table4, the λ values of channels 4, 10, and 12 are not in accordance with actual with far more larger than actual values. Thus, the channels 4, 10, and 12 ought to be eliminated because they are too invalid to be used in subsequent positioning calculations. In addition to removing the abnormal channel, the channels with larger error should also be taken in to account.Table 4
Time difference between station 8 and other stations.
Channel number
ΔTA
ΔTT
λ
1
0.036
0.0393
0.9169
2
0.034
0.0361
0.9414
3
0.048
0.0347
1.3828
4
0.468
0.0351
13.343
5
0.059
0.0372
1.5841
6
0.019
0.0408
0.4662
7
0.071
0.0081
8.7215
9
0.056
0.0084
6.6419
10
1.471
0.0164
89.7603
11
0.058
0.0247
2.3484
12
1.471
0.0329
44.7076


#### 5.3.3. Culling of Big Error Channel

By eliminating the channels 4, 10, and 12, the rest channels are optimized by using envelope curve, whose radius ought to be calculated firstly. RadiusRi is equal to ΔTT multiplied by the propagation wave velocity V. Table 5 and Figure 10 give the calculated results of the envelope curve. It can be obviously found that each envelope curve is closely distributed at the microseismic source (with red “∗” mark). Also, most of the envelope curves are in two intersecting forms.Table 5
Maximum value of circle radius.
Channel number
R
1
133.2
2
125.8
3
177.6
5
218.3
6
70.3
7
262.7
9
207.2
11
214.6By comparing the situation Figure3, and combining with theoretical analysis, it can be concluded that the common area of envelope intersects is the microseismic source distributed area, which can be further verified from Figure 13, i.e., all envelope curves are covered within the source’s area except channel 6. The channel 6 cannot be used in locating calculations, since its waveform is suspected to be an electrical pulse waveform.Figure 13
Envelope curve of MS event potential distribution area (calculated based on-time difference and wave velocity).The above analysis results show that channels 4, 10, and 12 cannot be accurately detected at the initial arrival, and error of channel 6 is even larger. These channels are not suitable for positioning calculations and treated as invalid channels. The channels 1∼3, 5, 7∼9, and 11 are effective channels. Therefore, the above channels are finally selected for positioning calculations of this microseismic event.


### 5.4. Contrastive Analysis

Manual picking method is used for the initial arrival pickup of the aforementioned events. The comparisons between pickup accuracy and time-consumption of manual and automatic methods are summarized in Table6. For actual measurement result, it takes 3 min to pick the above 12-channel waveform at initial arrival, which is very larger than automatic methods of only 0.039 s.Table 6
Comparison of picking methods in terms of speed and accuracy.
Number
Accuracy
Calculating speed
Remarks
Manual pick
6.6 ms (average error)
3 min
This contrast is limited to data of 12 channels
Automatic pick
0.039 s


## 5.1. Bandfilter: Denoising Module

The statistical analysis of field data shows that the frequency range of microseismic signal is affected by many factors, such as transmission distance, surrounding rock lithology, and microseismic source itself. These factors not only constrain the PT-FA, but also affect the accuracy to a great extent. Therefore, it is very necessary to denoise the microseismic data before the PT-FA beginning. It is famous that Butterworth band-pass filter is widely utilized in denoise of seismic signal. Taking the above factors in to consideration, a Bandfilter (code in MATLAB®) is used to process the microseismic signal by setting the effective frequency distribution range of 10∼200 Hz. Figure8 shows the band-pass denoise results of typical microseismic signal. It could be found that the main components and energy of signal distributes in range of 10∼200 Hz. The original signal shows a sparse distribution in the high-frequency region. However, but the main components in the high-frequency region is eliminated. According to the method described in the literature [16], the calculated SNR = 35.6126 db is higher and indicates that the denoising signal SNR has been greatly improved.Figure 8
Waveform and spectrum figures: signal and its frequency spectrum before (a) and after (b) denoising.
(a)
(b)Two typical microseismic signals were processed by the Bandfilter module, and the results of those can be seen in Figure9. As can be seen from the figures, this method can effectively remove the interference components and improve the signal-to-noise ratio (SNR). This is the basis of first arrival time pick automatically.Figure 9
Validation of typical signals: (a) denoising before (left) and after (right) signal 3#; (b) denoising before (left) and after (right) signal 7#.
(a)
(b)


## 5.2. First Arrival Time Picking Module

According to the mentioned calculation process in Section2, the PT-FA of studied microseismic event was carried out and is shown in Figure 10, which shows the “time-amplitude” curve within 12 channels. The red line and Stations 1∼12 denote the PT-FA results and 12 sensors, respectively. It can be obviously found that our algorithm can pick up the signal of Stations 1, 2, 3, 5, 6, 7, 8, 9, and 11. There are no effective waveform fluctuations found in channels 10 and 12 whose calculated initial arrival time is 1 s. Regarding channel 4 with lower SNR value, its signal noise is relatively low with a pickup time of 0.003 s. These results indicate that our algorithm is effective to realize the PT-FA of microseismic event under certain SNR.Figure 10
Waveforms and picked first arrival times acquired by the MS system with 12 sensors (vertical lines show the first arrival times).From Figure10 and Table 2, it could be concluded that the higher the SNR, the more accurate the picking time. On the contrary, it is difficult to accurately pick up or judge the SNR or abnormal signal of initial arrival. In Figure 10, channels of 4, 6, 10, and 12 can be directly observed by the naked eye. But if the algorithm returns values without judgment and optimization, it will seriously affect the accuracy of subsequent positioning. After the pickup of initial arrival, all arrival times are substituted in the positioning model for seismic source calculations. It could be found in Figure 10 that the judgments of channel 4, 10, and 12 are invalid according to the algorithm. The judgment process can be automatically realized through EJ&O step in Section 2.Table 2
First arrival time of every station (sensor).
Channel number
Arrival time (s)
1
0.507
2
0.505
3
0.519
4
0.003
5
0.530
6
0.490
7
0.542
8
0.471
9
0.527
10
1
11
0.529
12
1


## 5.3. Efficient Channel Optimization

The EJ&O of channels is to ensure the accuracy of positioning calculation. The initial arrival values of each time period are screened out, and invalid or error values are eliminated.


### 5.3.1. Choosing Reference Object

Since the sensors are located in a same working face and arranged at 30–50 m intervals, the time interval between triggering sensors is not large sensor as for the microseismic events happen in network layout. The channels with larger errors can be eliminated by using clustering calculations (k-means clustering algorithm,k = 2) on these trigger moments. The selection of reference channels is realized by means of clustering according to Table 3 and Figure 11, in which 12 initial arrivals are clustered into two types and “∗” signs denotes clustering center (Figure 8).Table 3
Result of first arrival times.
Channel number
Distance of station and cluster A
Cluster number
1
0.0020
1
2
0.0018
1
3
0.0032
1
4
0.2110
1
5
0.0064
1
6
7.67 × 10–4
1
7
0.0064
1
8
7.57 × 10–5
1
9
0.0042
1
10
0.2891
2
11
0.0044
1
12
0.2891
2Figure 11
Clustering of arrival times and choice of refer station. (a) Clustering results using k-means; (b) distance (normalized value) of current station and cluster A.
(a)
(b)It can be observed from Figure11 that all the channels are clustered into one type except Channels 10 and 12. As the elements in cluster A are much more than those in cluster B, class A with larger element number are selected to continue the study. Then, the closest point (station) to the cluster center (mean value at that time in the cluster) in cluster is taken as the reference, and class B is excluded. That is to say, once the number of elements in two clusters equals each other, the centers (two clusters) can be selected as the cluster center in order to find the distance (Ri). Figure 11(b) gives the distance between each element to its cluster center (cluster A), where 8# station having a minimum distance of 7.57E-4 can be used as reference station.


### 5.3.2. Culling of Dead Channel

After the treatment of first step, it is necessary to solve the problem of incorrect initial arrival time value or large error. Figure12 compares the actual arrival time difference, ΔTA (actual value) between 8# and other stations. The theoretical ΔTT of channels 4, 10, and 12 are lower than the actual ΔTT. By assuming λ = ΔTA/ΔTT, the channels with larger errors can be judged as λ > 10.Figure 12
Comparison histogram of theoretical maximum value and actual value.From the calculated results of relevant parameters of arrival time difference in Table4, the λ values of channels 4, 10, and 12 are not in accordance with actual with far more larger than actual values. Thus, the channels 4, 10, and 12 ought to be eliminated because they are too invalid to be used in subsequent positioning calculations. In addition to removing the abnormal channel, the channels with larger error should also be taken in to account.Table 4
Time difference between station 8 and other stations.
Channel number
ΔTA
ΔTT
λ
1
0.036
0.0393
0.9169
2
0.034
0.0361
0.9414
3
0.048
0.0347
1.3828
4
0.468
0.0351
13.343
5
0.059
0.0372
1.5841
6
0.019
0.0408
0.4662
7
0.071
0.0081
8.7215
9
0.056
0.0084
6.6419
10
1.471
0.0164
89.7603
11
0.058
0.0247
2.3484
12
1.471
0.0329
44.7076


### 5.3.3. Culling of Big Error Channel

By eliminating the channels 4, 10, and 12, the rest channels are optimized by using envelope curve, whose radius ought to be calculated firstly. RadiusRi is equal to ΔTT multiplied by the propagation wave velocity V. Table 5 and Figure 10 give the calculated results of the envelope curve. It can be obviously found that each envelope curve is closely distributed at the microseismic source (with red “∗” mark). Also, most of the envelope curves are in two intersecting forms.Table 5
Maximum value of circle radius.
Channel number
R
1
133.2
2
125.8
3
177.6
5
218.3
6
70.3
7
262.7
9
207.2
11
214.6By comparing the situation Figure3, and combining with theoretical analysis, it can be concluded that the common area of envelope intersects is the microseismic source distributed area, which can be further verified from Figure 13, i.e., all envelope curves are covered within the source’s area except channel 6. The channel 6 cannot be used in locating calculations, since its waveform is suspected to be an electrical pulse waveform.Figure 13
Envelope curve of MS event potential distribution area (calculated based on-time difference and wave velocity).The above analysis results show that channels 4, 10, and 12 cannot be accurately detected at the initial arrival, and error of channel 6 is even larger. These channels are not suitable for positioning calculations and treated as invalid channels. The channels 1∼3, 5, 7∼9, and 11 are effective channels. Therefore, the above channels are finally selected for positioning calculations of this microseismic event.


## 5.3.1. Choosing Reference Object

Since the sensors are located in a same working face and arranged at 30–50 m intervals, the time interval between triggering sensors is not large sensor as for the microseismic events happen in network layout. The channels with larger errors can be eliminated by using clustering calculations (k-means clustering algorithm,k = 2) on these trigger moments. The selection of reference channels is realized by means of clustering according to Table 3 and Figure 11, in which 12 initial arrivals are clustered into two types and “∗” signs denotes clustering center (Figure 8).Table 3
Result of first arrival times.
Channel number
Distance of station and cluster A
Cluster number
1
0.0020
1
2
0.0018
1
3
0.0032
1
4
0.2110
1
5
0.0064
1
6
7.67 × 10–4
1
7
0.0064
1
8
7.57 × 10–5
1
9
0.0042
1
10
0.2891
2
11
0.0044
1
12
0.2891
2Figure 11
Clustering of arrival times and choice of refer station. (a) Clustering results using k-means; (b) distance (normalized value) of current station and cluster A.
(a)
(b)It can be observed from Figure11 that all the channels are clustered into one type except Channels 10 and 12. As the elements in cluster A are much more than those in cluster B, class A with larger element number are selected to continue the study. Then, the closest point (station) to the cluster center (mean value at that time in the cluster) in cluster is taken as the reference, and class B is excluded. That is to say, once the number of elements in two clusters equals each other, the centers (two clusters) can be selected as the cluster center in order to find the distance (Ri). Figure 11(b) gives the distance between each element to its cluster center (cluster A), where 8# station having a minimum distance of 7.57E-4 can be used as reference station.


## 5.3.2. Culling of Dead Channel

After the treatment of first step, it is necessary to solve the problem of incorrect initial arrival time value or large error. Figure12 compares the actual arrival time difference, ΔTA (actual value) between 8# and other stations. The theoretical ΔTT of channels 4, 10, and 12 are lower than the actual ΔTT. By assuming λ = ΔTA/ΔTT, the channels with larger errors can be judged as λ > 10.Figure 12
Comparison histogram of theoretical maximum value and actual value.From the calculated results of relevant parameters of arrival time difference in Table4, the λ values of channels 4, 10, and 12 are not in accordance with actual with far more larger than actual values. Thus, the channels 4, 10, and 12 ought to be eliminated because they are too invalid to be used in subsequent positioning calculations. In addition to removing the abnormal channel, the channels with larger error should also be taken in to account.Table 4
Time difference between station 8 and other stations.
Channel number
ΔTA
ΔTT
λ
1
0.036
0.0393
0.9169
2
0.034
0.0361
0.9414
3
0.048
0.0347
1.3828
4
0.468
0.0351
13.343
5
0.059
0.0372
1.5841
6
0.019
0.0408
0.4662
7
0.071
0.0081
8.7215
9
0.056
0.0084
6.6419
10
1.471
0.0164
89.7603
11
0.058
0.0247
2.3484
12
1.471
0.0329
44.7076


## 5.3.3. Culling of Big Error Channel

By eliminating the channels 4, 10, and 12, the rest channels are optimized by using envelope curve, whose radius ought to be calculated firstly. RadiusRi is equal to ΔTT multiplied by the propagation wave velocity V. Table 5 and Figure 10 give the calculated results of the envelope curve. It can be obviously found that each envelope curve is closely distributed at the microseismic source (with red “∗” mark). Also, most of the envelope curves are in two intersecting forms.Table 5
Maximum value of circle radius.
Channel number
R
1
133.2
2
125.8
3
177.6
5
218.3
6
70.3
7
262.7
9
207.2
11
214.6By comparing the situation Figure3, and combining with theoretical analysis, it can be concluded that the common area of envelope intersects is the microseismic source distributed area, which can be further verified from Figure 13, i.e., all envelope curves are covered within the source’s area except channel 6. The channel 6 cannot be used in locating calculations, since its waveform is suspected to be an electrical pulse waveform.Figure 13
Envelope curve of MS event potential distribution area (calculated based on-time difference and wave velocity).The above analysis results show that channels 4, 10, and 12 cannot be accurately detected at the initial arrival, and error of channel 6 is even larger. These channels are not suitable for positioning calculations and treated as invalid channels. The channels 1∼3, 5, 7∼9, and 11 are effective channels. Therefore, the above channels are finally selected for positioning calculations of this microseismic event.


## 5.4. Contrastive Analysis

Manual picking method is used for the initial arrival pickup of the aforementioned events. The comparisons between pickup accuracy and time-consumption of manual and automatic methods are summarized in Table6. For actual measurement result, it takes 3 min to pick the above 12-channel waveform at initial arrival, which is very larger than automatic methods of only 0.039 s.Table 6
Comparison of picking methods in terms of speed and accuracy.
Number
Accuracy
Calculating speed
Remarks
Manual pick
6.6 ms (average error)
3 min
This contrast is limited to data of 12 channels
Automatic pick
0.039 s


## 6. Field Data Analysis and Application

In order to evaluate the practical application effect of hierarchical detection and pickup algorithm, microseismic events in mining site (265 sets of waveform data) were selected for calculation on arrivals. The computing environment was as follows: Windows 10 professional edition of the operating system, Intel(R) Core(TM) i7-7700k, CPU 4.20 GHz, and memory 16 GB.


### 6.1. Algorithm Testing

To quantitatively compare the accuracy of several methods, 265 groups of data (time-consuming: 0.828 s) were picked up manually and compared with the proposed method in Figure14, where Figure 14(a) shows the distributed relationship between manual and model calculated results and Figure 14(b) shows arrival errors distribution (within the range of 50 ms errors).Figure 14
P arrival residual plot for the algorithm. (a) Comparison of automatic pick and manual pick methods. (b) Error statistical analysis of automatic pick.
(a)
(b)It is obvious that a large number of automatic picking results (red∗) fall into the blue circle of manual pickups, indicating that the new proposed method is similar to manual pickups. Furthermore, the final error analysis shows the distribution efficiency of errors in −5∼5 ms is 71.7% (190 groups) −10∼10 ms is 83.4% (221 groups), and −20∼20 ms is 90.57% (240 groups). The high efficiencies illustrate that the proposed method in this study has higher precision and can replace the manual method.


### 6.2. Contrastive Analysis of Location Results

In order to verify the effects of channel optimization on final positioning calculations, the former calibration gun is used as a reference to compare the difference between the real value and calculated value before and after optimizations. The corresponding results can be considered as an evaluation criterion to judge the advantages and disadvantages of the proposed hierarchical detection and picking method. The proposed positioning method is a four-four combination positioning method: before the optimization, all channels are selected for calculation; after optimization, optimized channels are selected for positioning calculations. The final calculation results are shown in Table7.Table 7
Comparison of results before and after optimization.
Item
Source location (m)
Remarks
X axial
Y axial
Z axial
Real value
483.62
395.87
−189.14
—
Before
Calculated value
487.76
382.59
−173.25
Manual picking channels: 1, 2, 3, 5, 6, 7, 8, 9, 11
Error
−4.14
13.28
−15.89
After
Calculated value
490.29
389.43
−187.56
Automatic picking channels: 1, 2, 3, 5, 7, 8, 9, 11
Error
6.67
−6.44
1.58It can be seen that the selection of different channels has a great effect on final calculation results. Error between the original and real microseismic sources before optimization is about 21.12 m, which is larger than the result 9.41 m after optimization. There are two major defects in manual pickup arrival and effective channel. One lies in the fact that the manual method takes a longer time. The other one is that it is unable to screen channels with large difference in time value, which will seriously affect the processing rate of microseismic data.In summary, the method proposed in this study can solve the problem of PT-FA in mine microseismic wave to some extent. Our new method can be suggested as a precondition for microseismic positioning calculations. However, some improvements ought to be further developed in future work, such as effective waveforms identification and selection and judgment of external field microseismic events.


## 6.1. Algorithm Testing

To quantitatively compare the accuracy of several methods, 265 groups of data (time-consuming: 0.828 s) were picked up manually and compared with the proposed method in Figure14, where Figure 14(a) shows the distributed relationship between manual and model calculated results and Figure 14(b) shows arrival errors distribution (within the range of 50 ms errors).Figure 14
P arrival residual plot for the algorithm. (a) Comparison of automatic pick and manual pick methods. (b) Error statistical analysis of automatic pick.
(a)
(b)It is obvious that a large number of automatic picking results (red∗) fall into the blue circle of manual pickups, indicating that the new proposed method is similar to manual pickups. Furthermore, the final error analysis shows the distribution efficiency of errors in −5∼5 ms is 71.7% (190 groups) −10∼10 ms is 83.4% (221 groups), and −20∼20 ms is 90.57% (240 groups). The high efficiencies illustrate that the proposed method in this study has higher precision and can replace the manual method.


## 6.2. Contrastive Analysis of Location Results

In order to verify the effects of channel optimization on final positioning calculations, the former calibration gun is used as a reference to compare the difference between the real value and calculated value before and after optimizations. The corresponding results can be considered as an evaluation criterion to judge the advantages and disadvantages of the proposed hierarchical detection and picking method. The proposed positioning method is a four-four combination positioning method: before the optimization, all channels are selected for calculation; after optimization, optimized channels are selected for positioning calculations. The final calculation results are shown in Table7.Table 7
Comparison of results before and after optimization.
Item
Source location (m)
Remarks
X axial
Y axial
Z axial
Real value
483.62
395.87
−189.14
—
Before
Calculated value
487.76
382.59
−173.25
Manual picking channels: 1, 2, 3, 5, 6, 7, 8, 9, 11
Error
−4.14
13.28
−15.89
After
Calculated value
490.29
389.43
−187.56
Automatic picking channels: 1, 2, 3, 5, 7, 8, 9, 11
Error
6.67
−6.44
1.58It can be seen that the selection of different channels has a great effect on final calculation results. Error between the original and real microseismic sources before optimization is about 21.12 m, which is larger than the result 9.41 m after optimization. There are two major defects in manual pickup arrival and effective channel. One lies in the fact that the manual method takes a longer time. The other one is that it is unable to screen channels with large difference in time value, which will seriously affect the processing rate of microseismic data.In summary, the method proposed in this study can solve the problem of PT-FA in mine microseismic wave to some extent. Our new method can be suggested as a precondition for microseismic positioning calculations. However, some improvements ought to be further developed in future work, such as effective waveforms identification and selection and judgment of external field microseismic events.


## 7. Conclusions

Aiming at the problem of large amount and low precision on waveforms PT-FA pickup, a new model is proposed in this paper. Through theoretical analysis and engineering example verification, following conclusions were drawn:(1)
This paper provides an improved method of PK-FA which has four optimization modules: DNS, FSC, AP, and EJ&O module. STA/LTA (FSC module) was used to quickly search for the initial moment of first arrival time; AIC (AP module) provides an accurate pick method and solve the problem about multiple peak values of AIC eigenvalue. This provides some helps for subsequent accurate pickup and greatly reduces the calculated amount and time consuming of data process.(2)
Bandfilter module can be used for mining microseismic signals which have the limit frequency distribution. This method can effectively remove the interference components and improve the signal-to-noise ratio (SNR). This is the basis of first arrival time pick automatically.(3)
EJ&O of channels is to ensure the accuracy of positioning calculation. This module has the judgment criteria: cluster analysis, DIF judgment, and envelop curve. Based on choice of reference object, culling of dead channels, and big error channels, the left channels are finally selected for positioning calculations of the microseismic event(4)
Comparing manual and automatic picks in terms of speed and accuracy using field 265 sets of microseismic data, computation cost (0.828 s) is far less than manual method (15 s per set), precision of manual and automatic represented by picking error is up to 71.7% (between −5 ms and 5 ms, 190 sets), 83.4% (between −10 ms and 10 ms, 221 sets), and 90.57% (between −20 ms and 20 ms, 240 sets). This study shows that the method presented in this paper is feasible to precise automatic P-phase determination and is a basis on automatic and precision location.



---

*Source: 3780875-2019-06-19.xml*