Identifying P phase arrival of weak events: The Akaike Information Criterion picking application based on the Empirical Mode Decomposition

Computers & Geosciences 100 (2017) 57–66

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

Research paper

Identifying P phase arrival of weak events: The Akaike Information Criterion picking application based on the Empirical Mode Decomposition

MARK

⁎

Xibing Lia,b, Xueyi Shanga, , A. Morales-Estebanc, Zewei Wanga a b c

School of Resources and Safety Engineering, Central South University, China Hunan Key Lab of Resources Exploitation and Hazard Control for Deep Metal Mines, China Department of Building Structures and Geotechnical Engineering, University of Seville, Spain

A R T I C L E I N F O

A BS T RAC T

Keywords: P phase arrival picking Weak events and microseisms Empirical Mode Decomposition (EMD) Akaike Information Criterion (AIC) picker Discrete Wavelet Transform (DWT)

Seismic P phase arrival picking of weak events is a difficult problem in seismology. The algorithm proposed in this research is based on Empirical Mode Decomposition (EMD) and on the Akaike Information Criterion (AIC) picker. It has been called the EMD-AIC picker. The EMD is a self-adaptive signal decomposition method that not only improves Signal to Noise Ratio (SNR) but also retains P phase arrival information. Then, P phase arrival picking has been determined by applying the AIC picker to the selected main Intrinsic Mode Functions (IMFs). The performance of the EMD-AIC picker has been evaluated on the basis of 1938 micro-seismic signals from the Yongshaba mine (China). The P phases identified by this algorithm have been compared with manual pickings. The evaluation results confirm that the EMD-AIC pickings are highly accurate for the majority of the micro-seismograms. Moreover, the pickings are independent of the kind of noise. Finally, the results obtained by this algorithm have been compared to the wavelet based Discrete Wavelet Transform (DWT)-AIC pickings. This comparison has demonstrated that the EMD-AIC picking method has a better picking accuracy than the DWT-AIC picking method, thus showing this method's reliability and potential.

1. Introduction Micro-seismic monitoring is an advanced and effective tool for rock burst and other dynamic disaster monitoring. It has been widely used in mining, oil and gas exploitation, slope stability and tunneling. The major technical issues of this tool are monitoring planning, data processing and micro-seismic event location ((Ge, 2005)). The P phase arrival identification is one of the key issues in data processing and the fundamental step to calculate the event's time, magnitude, location, and other seismic parameters ((Gou et al., 2011; Alvarez et al., 2013; Yue et al., 2014)). Although such a task can be done visually by analysts, manual picking is time consuming and its quality depends on the analysts’ experience and subjectivity. Therefore, a large volume of micro-seismic data requires a reliable and automatic identifying method. Nevertheless, background and stationary noises often impact adversely on the micro-seismograms, making the identification tough and unreliable. In this paper, an algorithm based on a brand new time frequency analysis means-Empirical Mode Decomposition (EMD) and the wellknown Akaike Information Criterion (AIC, (Akaike, 1971; Akaike, 1974)) picker- to identify P phase arrival is proposed. It has been

⁎

called the EMD-AIC picker. In contrast to other algorithms which depend on either the algorithm itself, the nature of the noise or the wavelet coefficients, EMD is a novel self-adaptive method that works better. The performance of the EMD-AIC picker has been evaluated on the basis of 1938 micro-seismic signals from the Yongshaba mine (China). The results obtained with this algorithm have been compared with the Discrete Wavelet Transform (DWT, (Mallat, 1989))-AIC pickings. The results show that the EMD-AIC picking method efficiently provides a good estimate of the P phase arrival picking, especially for low Signal to Noise Ratio (SNR) signals. 2. State of the art Many studies have been done to improve the P phase arrival picking performance and various methods have been proposed. The most widely P phase arrival picking method used is the Short Term Average (STA) - Long Term Average (LTA) ratio proposed by (Allen, 1978). In this method, a Characteristic Function (CF) that considers the changes both in amplitude and frequency when P phase comes is introduced. However, this method is not sensitive to a seismic signal that has a low SNR. (Baer and Kradolfer, 1987) modified the envelop function of

Corresponding author. E-mail addresses: [email protected] (X. Li), [email protected] (X. Shang), [email protected] (A. Morales-Esteban), [email protected] (Z. Wang).

http://dx.doi.org/10.1016/j.cageo.2016.12.005 Received 13 May 2016; Received in revised form 9 November 2016; Accepted 5 December 2016 Available online 07 December 2016 0098-3004/ © 2016 Elsevier Ltd. All rights reserved.

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

(a) Amplitude

Allen´s (1978) method. To do so, an incorporated dynamic threshold that increases the sensitivity to P phase arrival was proposed. (Earle and Shearer, 1994) calculated the STA/LTA time series of a smoothed envelop function generated from a Hilbert transform of the seismogram. They demonstrated that it could identify the primary phases as well as the secondary phases. (Ross and Ben-Zion, 2014) applied filters to remove P-wave energy from seismograms and utilized the STA/LTA and PAI-K detectors in tandem to determine the P phase arrival. Besides the above automatic P phase arrival picking methods, the autoregressive (AR) model is the one most widely used. This method is based on the AIC and assumes that the intervals before and after P phase onset have a worse model fit. (Sleeman and van Eck, 1999) joined the AR modeling of the noise and the fixed seismic signal. Later, they applied the AIC as a P phase arrival onset parameter. The so-called AR-AIC picker was compared with the well-established picker of (Baer and Kradolfer, 1987), obtaining a much better result for P phase onsets based on 1109 seismic pickings. (Zhang et al., 2003) developed a P phase detection and picking algorithm based on the AIC picker and wavelet transform. This takes advantage of those wavelet coefficients which retain P phase arrival well at high resolutions. Then, the AIC picker was applied to different scales. The wavelet-based AIC picker achieves good results even for noisy seismograms. (Sedlak et al., 2009; Sedlak et al., 2013) used a specific characteristic function to improve the sensitivity of the AIC picker to changes in frequency as well as amplitude, obtaining reliable results. Time frequency transformations, such as Short Time Fourier Transform (STFT) and wavelet transform, choose a basis that can represent the seismogram locally in both the time and frequency domains. These have been widely used in seismic P phase arrival identification. (Hafez et al., 2009) proposed an automatic P phase arrival picking algorithm which applies the spectrum-ratio on the timefrequency sub-band, and used the spectrum-ratio as parameter to identify the P phase arrival and to avoid false alarms. (Galiana-Merino et al., 2008) combined higher order statistics and Stationary Wavelet Transform (SWT) to identify P phase arrival. They obtained good results even for a very low SNR signal. (Hafez et al., 2010; Hafez et al., 2013) introduced two wavelet-based algorithms grounded on the Multi Resolution Analysis (MRA) of the Haar wavelet and the MRA of Maximal Overlap Discrete Wavelet Transform (MODWT), respectively. From these wavelet details, the first stage detail was selected to detect P phase arrival as the arrival is clear in that detail. (Karamzadeh et al., 2013) proposed a CF by using the stack of the Continuous Wavelet Transform (CWT) coefficients that cover the frequency content of the desired P phases. Numerous other approaches have been proposed, including higher order statistics ((Saragiotis et al., 2002; Saragiotis et al., 2004; GalianaMerino et al., 2008; Küperkoch et al., 2010; Baillard et al., 2014; Li et al., 2016)), polarization analysis ((Vidale, 1986; Magotra et al., 1987; Kulesh et al., 2007)), the damped predominant period (Tpd) method ((Lockman and Allen, 2005; Hildyard et al., 2008; Hildyard and Rietbrock, 2010)), fuzzy logic theory ((Chu and Mendel, 1994)), maximum likelihood methods ((Christoffersson et al., 1988; Roberts et al., 1989)), cumulative sum-based algorithms ((Der and Shumway, 1999)), damping energy-based method ((Kalkan, 2016)), multiband frequency analysis ((Alvarez et al., 2013; Romero et al., 2016; García et al., 2016)), artificial neural networks-based picking procedure ((Zhao and Takano, 1999; Gentili and Michelini, 2006)), and waveletbased picking methods ((Galiana-Merino et al., 2008; Hafez et al., 2010; Hafez et al., 2013; Karamzadeh et al., 2013)). A review of P phase arrival picking can be found in (Withers et al., 1998), (Tselentis et al., 2012) and (Akram and Eaton, 2016).

x 10−6 2 1317

0 −2 0

500

1000 1500 2000 2500 3000 3500 4000

Samples (b) Amplitude

8

x 10−4 8

6

6

Enlarged

4

4

Power interference of 50 Hz

2

2 0

x 10−4

0 0

0

500

20

1000

40

60

80

1500

100 120 140 160 180 200

2000

2500

3000

Frequency (Hz) Fig. 1. Micro-seismic signal and its amplitude spectrum. (a) Micro-seismic signal. The vertical line corresponds to the manual P phase arrival picking and the number indicated by an arrow corresponds to the P phase arrival point; (b) amplitude spectrum of the micro-seismic signal. The amplitude spectrum of frequency band in [0,200] Hz is enlarged.

picker are explained in detail. 3.1. The EMD The EMD is a fast, effective and self-adaptive time-frequency analysis method for non-linear and non-stationary time series. It has been developed from the simple assumption that any signal can be decomposed into a number of Intrinsic Mode Functions (IMFs). Each of which must satisfy the following two definitions: 1) in the whole data set, the number of extreme and zero-crossing points must be either equal or be different by 1 at the most; 2) at any point, the mean value of the maximum envelope and minimum envelope is zero. The signal x (t ) can be decomposed by EMD as follows ((Huang et al., 1998; Huang et al., 1999)): (1) Identify all the local extrema. Then, connect all the local maxima by a cubic spline-line to form the upper envelope xmax(t ). Repeat this procedure for all the local minima to obtain the lower envelope xmin(t ). (2) The mean of the upper and lower envelopes is marked as m1(t ). Next, separate m1(t ) from x (t ) to obtain the first component h1(t ):

m1(t ) = (xmax(t ) + xmin(t ))/2

(1)

h1(t ) = x (t ) − m1(t )

(2)

(3) Ideally, h1(t ) should be an IMF. Then, h1(t ) is the first IMF component of x (t ). Nonetheless, h1(t ) rarely satisfies the IMF property, which can be measured by the Standard Deviation (SD) between h1(k −1)(t ) and h1k (t ) . More commonly, h1(t ) is treated as the original signal, then Steps (1) and (2) are repeated k times until h1k (t ) satisfies the IMF property:

⎧ h11(t ) = h1(t ) − m11(t ) ⎪ ⎪ h12(t ) = h11(t ) − m12(t ) ⎨ ⎪⋮ ⎪ h (t ) = h 1(k −1)(t ) − m1k (t ) ⎩ 1k

3. Theoretical fundamentals

(3)

where the standard deviation (SD) can be obtained by equation (t ) − h1k (t ) 2 ⎤ T ⎡h SD = ∑t =0 ⎢ 1(k −1)2 ⎥, and a typical value for SD can be set h1(k −1)(t ) ⎣ ⎦

In this section the theoretical fundamentals that support the method proposed in this paper are exposed. The EMD and the AIC 58

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

Fig. 2. Scales of EMD and DWT (blue lines) and their AIC time series (pink lines). The vertical and vertical-dashed lines represent manual pickings and AIC pickings, respectively. The numbers indicated by an arrow correspond to the P phase arrival points. (a) IMFs (c1, c2 ,…, c8 ) and the residual component r8 of EMD and their AIC time series; (b) detail components (D1, D2 ,…, D8 ) and the approximate component A8 of DWT and their AIC time series. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

between 0.2 and 0.3. After, h1k (t ) is designated as the first IMF component of x (t ):

c1(t ) = h1k (t )

the predetermined value of substantial consequence, or rn(t ) becomes a monotonic function. By summing up Eqs. (5) and (6), the signal x (t ) can be reconstructed by:

(4)

n

(4) Separate c1(t ) from x (t ) to get the rest of the data

x (t ) =

∑ cj(t ) + rn(t ) j =1

r1(t ) = x (t ) − c1(t )

The IMFs c1(t ), c2(t ),…, cn(t ) correspond to different frequency bands ranging from high to low. They change with the variation of signal x (t ), while rn(t ) is the mean trend of signal x (t ).

(5)

Then, r1(t ) is treated as the original signal, and the above steps are repeated to obtain the second IMF component c2(t ) of signal x (t ). The above steps are repeated n times to obtain n-IMFs of signal x (t ):

⎧ r2(t ) = r1(t ) − c2(t ) ⎪ ⎪ r3(t ) = r2(t ) − c3(t ) ⎨ ⎪⋮ ⎪ r (t ) = r (n −1)(t ) − cn(t ) ⎩n

(7)

3.2. The AIC picker The AIC picker supposes that the intervals before and after the P phase arrival are two different stationary time series, and the AIC values of the two-interval model of the micro-seismic signal x (t ) (t=1, 2, …, N) are given as the function of the merging point k ((Sleeman and van Eck, 1999)):

(6)

AIC (k ) = (k − M )⋅ log(σ1,2 max ) + (N − M − k )⋅ log(σ2,2 max ) + C2

This procedure will be accomplished when cn(t ) or rn(t ) is less than 59

(8)

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

(a)

4. The EMD-AIC picker −4

4

x 10

c4

Amplitude

2 0 1 0.5 0

x 10 −5

1

c3

0.5

c2

0 180 500

0 x 10−3

500 1

c5

c6

1000

0 60

3000

2500

3000

c8

150

100

2500

2000

c7

100

50

2000

1500

r8

0.5

0

1500

1000 x 10 −5

The EMD-AIC picking algorithm presented in this paper takes advantage of the IMF of the micro-seismic signal and the AIC picking method. In Fig. 1, an example of micro-seismic signal with a relatively clear P phase arrival at point 1317 and its amplitude spectrum calculated with the Fast Fourier Transform (FFT) are shown. The amplitude spectrum demonstrates that the main information of the P phase is distributed below 80 Hz. However, the P phase is contaminated with the power frequency interference, which makes the automatic P phase onset hard work. The EMD possesses a unique ability for separating high frequency components and residual trend signals. Therefore, the EMD has been applied as an attempt to clear the P phase arrival. In addition, the DWT combined with the AIC picker ((Maeda, 1985)), called DWT-AIC picker, has been selected for comparison and to test the performance of the method proposed. The EMD-AIC and the DWT-AIC picking algorithms operate as follows. Given the micro-seismic signal x (n )(n = 1, 2, ... , N ) with a sampling frequency of 6000 Hz, the EMD and the DWT have been applied to obtain the 8-IMFs (c1, c2 ,…, c8) and the 8 detail components (D1, D2 , …, D8 ). Then, the AIC time series calculated with Eq. (9) has been used to confirm the P phase arrival time, which has been shown in Fig. 2. For a better comparison, their amplitude spectra have been depicted in Fig. 3. In general, the components include different range frequencies from high to low and the scales of the P phase, having a frequency band approximately between 10 and 200 Hz. On the one hand, c1, c2 , D1, D2 and D3 correspond to high frequencies that are very noisy and have some transients which may produce a false P phase onset. On the other hand, c7, c8, r8 , c8, D8 and A8 correspond to low frequencies that have unclear P phase arrivals. Therefore, the components discussed above should not be taken into account for the P phase arrival picking. In order to select components that retain the P phase arrival well, the onset results of the rest of the components based on the EMD-AIC and on the DWT-AIC picking methods have been listed in Table 1. The residuals were calculated by comparing the automatically picked arrivals and the manually picked arrivals, which provides a qualitative evaluation, as the results by analysts are more reliable ((Sleeman and van Eck, 1999)). It is obvious that the c3 and D4 obtained the worst results with picking errors of 158 and 98 samples, while the others have a maximum residue of 50 samples. In addition, the EMD-based method is better than the DWT-based method, having average residues of 10.33 and 44 samples for the other components. Therefore, the c4 , c5, c6 , D5, D6 and D7 have been selected for the AIC picker, thus improving the P phase arrival detection accuracy and reducing the computational time.

c1

200

250

150

200

300

250

300

Frequencey (Hz) (b)

−5

1.5

x 10

−6

6 x 10

Amplitude

1

D4

2 0 300 500

0.5 0 1

0 x 10−3

500

D1

D2 1000

1000

1500

2000

1500

2500

2000

3000

2500

3000

−4

D8

0.5 0

D3

4

A8

1.5 x 10

D7

1 0 50

0

D6

0.5

50

100

100

D5 150

200

150

250

200

300

250

300

Frequencey (Hz) Fig. 3. Amplitude spectra of the EMD and the DWT components, the inner graphs are the amplitude spectra with narrower frequency ranges. (a) Amplitude spectra of the EMD components (c1, c2 ,…, c8 and r8 ); (b) amplitude spectra of the DWT components (D1, D2 ,…,

D8 and A8 ). Table 1 Onset results of the EMD- and the DWT-based AIC picking methods. Method

AIC picking

AIC picking residual

EMD

c3 1475

c4 1331

c5 1328

c6 1323

c3 158

c4 14

c5 11

c6 6

DWT

D4 1415

D5 1350

D6 1268

D7 1267

D4 98

D5 33

D6 −49

D7 −50

where M is the order of an AR to fit the data, σ1,2 max and σ2,2 max indicate the variance of the prediction errors of the time series in intervals [M + 1, k ] and [k + 1, N − M ], and C2 is a constant. In order to obtain the AIC values, the AR coefficients and the order of M must be calculated. These have a high degree of computational complexity. In contrast to the above AIC picker, (Maeda, 1985) calculates the AIC function directly from the seismogram without using the AR coefficients, and the AIC function is defined as:

AIC (k ) = k⋅ log{var(x[1, k ])} + (N − k − 1)⋅ log{var(x[k + 1, N ])}

5. Experimental results and comparisons In this section the dataset used in this study is described. Then, the method used to compare the accuracy is defined. Next, the SNR performance has been tested. Finally, the performance of the method has been evaluated.

(9) 5.1. Dataset

where var(x[1, k ]) is the variance of time series x(1), x(2),…, x (k ), and var(x[k + 1, N ]) is the variance of time series x (k + 1), x (k + 2),…, x (N ). Considering the above two AIC picking methods, the latter is selected in this paper for its efficiency, and the AIC picker chooses the global minimum as the P phase onset. For a seismic signal with a clear P phase, the AIC method is likely to be very accurate, while for a relatively low SNR seismic signal the AIC picker may have large errors. Therefore, there is a great need to apply the AIC method to multiscale analysis and reduce the effect of noises.

The tested micro-seismic signals have been collected from the Institute of Mine Seismology (IMS) system at the Yongshaba mine. This is located in the Guizhou province, China (26°38´N, 106°37´E), and its framework is shown in Fig. 4. The IMS system contains one micro-seismic server, one data exchange center, four master stations and four deputy stations, two triaxial geophones and twenty-six uniaxial geophones. The geophones are distributed across the main exploitation stopes at the 930 m level,

60

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

Fig. 4. Framework of the Institute of Mine Seismology (IMS) system installed at Yongshaba mine. (a) Horizontal view of the micro-seismic monitoring system; (b) the configuration of the micro-seismic network used in the IMS system. The number near the geophone is the geophone number. The uniaxial and triaxial geophones link with their administering station by 4 and 8 core signal cables, respectively. The deputy station, the master station, the data exchange center and the micro-seismic server are connected by a 4 core optical fiber.

61

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

lower hinge

1000

(a)

upper hinge

35

600

lower adjacent value

mean value median value upper adjacent value

29.10% 32.10%

30

400

accuary≥30 25≤accuary<30 20≤accuary<25 15≤accuary<20 10≤accuary<15 ≤5≤accuary<10 accuary<5

25

Overall %

Residuals

800

200 0 -200 -400

20 30.45%

15 37.08%

10 EMD

-600

c4

c5

36.36%

DWT

c6

D5

D6

39.66%

5

D7

AIC picker

0 0

2

4

6

8

10

>12

Fig. 5. Box-plot of the EMD-AIC and the DWT-AIC picking results.

SNR the 1080 m level, and the 1120 m level, with an average network of 250 m×250 m. The uniaxial and triaxial geophones link with their administering station by 4 and 8 core signal cables, respectively. The deputy station, the master station, the data exchange center and the micro-seismic server are connected by a 4 core optical fiber. A geophone records the continuous voltage signal caused by vibrations such as rock fracture, blasting and mine drawing and transports the signal to its administering NetADC (network analogue‐to‐digital converter) through signal cable. Then, the NetADC samples the analog signal with a 6000 Hz frequency. The NetSP (network seismological processor) is an embedded computer that performs tasks such as data capturing from the NetADC, triggering, pre-trigger filtering and data transfering to the data exchange center. Finally, the signal is transported to the micro-seismic server for data storing and calculation. IMS Synapse, Trace, Ticker 3D and Jdi are the main data processing softwares and the IMS Trace uses the STA/LTA ((Allen, 1978)) to pick the P phase arrival. Nonetheless, the micro-seismic signals are often contaminated with various noises, such as truck transportation, mine drawing, blower vibration, blasting vibration and electrical signals, which make it hard to automatically identify the P phase arrival.

(b) 35

19.29% 20.46%

30

accuary≥30 25≤accuary<30 20≤accuary<25 15≤accuary<20 10≤accuary<15 ≤5≤accuary<10 accuary<5

Overall %

25 20 20.83%

15 23.60%

10

19.58% 27.59%

5 0 0

2

4

6

8

10

>12

SNR Fig. 6. Accuracy of the c5- and D6 - based picking methods in different SNR bands. Fill patterns represent the percentage of pickings in specified accuracy bands of the total of events. The numbers above the bar mean the percentages of picking accuracy less than 5 ms in that SNR band. (a) The accuracy of thec5- based picking method; (b) the accuracy

5.2. Accuracy comparison The EMD- and the DWT- based AIC picking methods have been tested and compared with the manual pickings of 1938 micro-seismic signals. The manual pickings are the average results of three experienced analysts’ pickings. The collected micro-seismic signals have different local magnitudes from −3 to 0.5 and occurred at different times from January 1st. 2014 to April 30th. 2014 to ensure that the method proposed sustains different seismic noise conditions. The EMD-AIC and the DWT-AIC picking residuals are shown in the boxplot of Fig. 5.

of the D6 - based picking method.

The c5 - and D6 - based picking methods have the smallest average residual and the c5 - based picking result is better, having a smaller upper hinge residual. The c4 - and D5- based picking methods have a nearly equal median residual. Nevertheless, the c4 - based picking method has a smaller lower hinge residual and a larger upper hinge residual. The c6 - based picking method has a smaller median residual

Table 2 The statistical results of picking residuals and their total costs. Method

Component

Statistical results of different picking residuals

Total cost

0 to 5 ms

5 to 10 ms

10 to 15 ms

15 to 20 ms

20 to 25 ms

25 to 30 ms

> 30 ms

EMD-AIC

c4 c5 c6

681 622 316

320 423 320

163 259 232

98 128 175

76 92 149

57 58 93

543 356 653

1120.3 930.6 1453.5

DWT-AIC

D5 D6 D7

509 403 220

398 359 254

235 308 238

126 255 223

83 165 209

52 80 160

535 368 634

1170.1 1112.0 1558.0

62

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

Table 3 P phase picking performance of 15 micro-seismic signals from the dataset. No.

SNR

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

6.20 1.22 3.57 4.00 2.78 3.64 4.92 4.63 11.30 6.49 9.79 2.69 1.67 8.01 7.93

Manual picking

1330 1221 1462 1363 1310 1247 1218 1379 1116 1177 1212 1031 1403 1470 1086

AIC picking

Picking residual

c5

D6

c5

D6

1349 1224 1547 1400 1334 1252 1453 1397 1158 1625 1207 1024 1409 1473 1113

1303 1346 1434 1344 2320 2313 1162 1268 1150 1608 1216 1098 1370 1406 1340

19 3 85 37 24 5 235 18 42 448 −5 −7 6 3 27

−27 125 −28 −19 1010 1066 −56 −111 34 431 4 67 −33 −64 254

∑ CFi i

⎧ 0.0, ⎪ ⎪ 0.2, ⎪ 0.4, ⎪ CFi = ⎨ 0.6, ⎪ 0.8, ⎪ ⎪1.0, ⎪1.5, ⎩

accuracy < 5ms 5ms ≤ accuracy < 10ms 10ms ≤ accuracy < 15ms 15ms ≤ accuracy < 20ms 20ms ≤ accuracy < 25ms 25ms ≤ accuracy < 30ms accuracy ≥ 30ms

Clear P arrival with little blower vibration and high frequency noises Not clear P arrival with mine drawing and high frequency noises Clear P arrival Not clear P arrival with shock and power frequency noises Not clear P arrival with power and high frequency noises Not clear P arrival with power frequency noise Not clear P arrival with high frequency noise Clear P arrival Very noisy environment, poor P arrival Poor clear P arrival with mine drawing, power and high frequency noises Clear P arrival Poor P arrival, noisy environment with blower vibration noise Clear P arrival with little high frequency noise Clear P arrival with blasting vibration noise Noisy environment, not clear P arrival

shown in Fig. 6. The picking performance is measured by the accuracy percentage of the total of events in each SNR band. The majority of SNRs are in the bands of [0,2] and (2, 4], with 32.09% and 29.26% of the total of events, respectively. For both methods, the higher the SNR, the more accurate is the onset, generally speaking. For example, the accuracy less than 5 ms of the c5- based picking method has made up to 29.10%, 32.10%, 30.45%, 37.08%, 39.66% and 36.36% at each SNR band. The c5- and D6 - based picking methods obtain good accuracies even when the SNRs are low and have a similar percentage of picking accuracy of less than 30 ms. Nonetheless, the c5- based picking method gets much better results in accuracy of less than 5 ms. In summary, the c5- based picking method adapts to micro-seismic signals of different SNRs very well.

than the D7- based picking method, while the D7- based picking method has a smaller upper hinge residual. Therefore, it is difficult to compare c4 - and D5-, and -c6 and D7- based picking methods. For a better comparison of the P phase onset quality based on different components, a cost minimization approach has been proposed. This method takes advantage of the cost of each picking residual and the minimum combination of cost functions corresponds to the best picking method. If a P phase accuracy of less than 30 ms is picked, then this will be assigned as an acceptable cost. Otherwise, it will be penalized by a high cost. The Total Cost-Function (TCF) and the CostFunction (CFi ) are defined as:

TCF =

Analyst comments

(10)

5.4. Performance of the method

The cost minimization approach has been applied to the 1938 picking residuals and the statistical results of certain picking residuals and their total costs have been listed in Table 2. It is obvious that there is a large number of picking residuals below 5 ms for the c4 - and c5based picking methods. This indicates that these two components retain P phase arrival information well. Moreover, the c5- based picking method has the smallest number statics of picking residuals of more than 30 ms. The c5- based picking method has the smallest total cost (930.6). The c4 -, D5- and D6 - based picking method have similar total costs (1120.3, 1170.1 and 1112.0, respectively). The c6 - based picking method has the worst picking result. Yet it is a little better than the D7based picking residuals. In conclusion, the EMD-based picking method is better than the DWT-based picking method, and the c5-AIC picker obtains the best picking result.

To determine the performance of the method proposed, some picking results have been listed in Table 3. For a better understanding of the waveform adaptation, six typical examples have been selected from Table 3. Their waveforms (Top), c5-AIC and D6 -AIC time series (Middle), and amplitude spectra of waveforms (blue lines), c5 (red lines) and D6 (green lines) (Bottom) have been depicted in Fig. 7. It can be noticed that both the c5- and D6 - based picking methods obtain good results in most cases. For example, a good picking accuracy is shown for clear P arrival signals (e.g., cases 1, 3, 8, 11 and 14). This is due to them being able to retain a good approximate P phase arrival for clear P arrival signals, as can be observed in Fig. 7(e). However, there can be large errors for unclear P phases arrival and poor onset phases (e.g., cases 7 and 10 for c5- based picking method, and cases 5, 6, 10 and 15 for D6 - based picking methods), as can be noticed in Fig. 7(d). What is more, the c5- based picking method is better than the D6 - based picking method given that it has a small average picking residual and a better adaption to different waveforms. This may be due to the fact that the c5 usually obtains a frequency which is more likely to the main frequency of original signals (e.g., Fig. 7(a), (b) and (f)).

5.3. The SNR effect

6. Discussion

SNR applicability is necessary to test the performance of automatic picking results. Yet the SNR definition varies among studies. In this paper, the SNR definition by (Gentili and Michelini, 2006) has been selected. The SNR is defined as the ratio of the mean absolute amplitude after and before the P phase onset. The mean value has been calculated in two 100 ms length windows in this paper. The accuracy of c5- and D6 - based picking methods in different SNR bands is

The above analysis has shown that the c5 of the EMD-based picking method is better than the D6 of the DWT- based picking method. It has a greater accuracy and adaptability to different waveforms and SNRs. It is known that the essence of DWT is the Fourier transform of an adjustable window. Once the wavelet basis function is selected, the main frequency band is only related to the sampling frequency of the signal. Furthermore, it has a low frequency resolution for high

(11)

63

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

(b)

2000

2500

3000

3500

0

−5.8

1224

1000

1500

2000

2500

3000

3500

0

−5.8 1346

−5 0

Amplitude

4000 x 104 −5.6

500

1000

1500

#2

0 0

2500

3000

3500

20

c5

D6

40

60

80

100 120 140 160 180 200

1500

4000 4 x 10 −5.5 −5.6

1000

1500

2000

2500

3000

3500

0

−5.7 4000 4 x 10 −5.6 −5.7

2320

500

1000

1500

2 1 0 0

2000

2500

3000

3500

−5.8 4000

Samples

x 10−4

c5

#5

20

D6

40

60

80

100 120 140 160 180 200

4000 4 x 10

1453

1000

1500

2000

2500

3000

3500

0

−5.8 4000 x 104 −5.8 −5.9

1162

500

1000

1500

2500

3000

3500

20

40

60

80

100 120 140 160 180 200

1500

1000

1500

2500

3000

3500

4000 4 x 10 −5.6

2500

3000

3500

2500

3000

3500

4000 4 x 10 −5.6 −5.7 −5.8 −5.9 4000

−5.8

1625

2000

1608

−1 0

4000

D6

1000

2000

0

Amplitude

c5

#7

1500

0 −5 0 −6 500 x 10 1

Samples

x 10−4 2 1 0 0

2000

1000

2

500

Samples

x 10−4

c5

#10

1 0 0

Frequency (Hz)

2000

AICc5

3500

c5

3000

D6

2500

AICc5

2000

AICD6

1500

−5.6

0

D6 Amplitude

1000

1177

AICD6

1218

−5 0 −7 500 x 10 5

x 10−7 10 #10 5 0 −5 0 −7 500 x 10 5

Amplitude

x 10 #7

0

c5

20

40

D6

60

80

100 120 140 160 180 200

Frequency (Hz)

(e)

3000

3500

−5.0 1207

1000

1500

2000

2500

3000

3500

0

−5.2 4000 4 x 10 −5.0 −5.2

1216

500

1000

1500

x 10

1 0 0

2500

3000

3500

c5

D6

1000

1500

2000

2500

3000

3500

0

4000 4 x 10 −5.2 −5.4

1113

1000

1500

2000

2500

3000

3500

0

−5.6 4000 4 x 10 −5.6 −5.7

1340

−5 0

Samples

−3 #11

2000

−5.4 4000

−2 0 −6 500 x 10 1

−1 0 −7 500 x 10 5

Amplitude

−2 0 2

4000 4 x 10 −4.8

5

500

1000

2000

2500

3000

3500

−5.8 4000

Samples

x 10−4 #15

1500

AICc5

2500

1086

0

c5

2000

x 10−6 #15

2

D6

1500

AICD6

1000

AICc5

1212

0 −5 0 −6 500 x 10 2

Amplitude

(f) x 10−6 5 #11 0 −5 0 −6 500 x 10 5

AICD6

Amplitude

5

−2 0

Amplitude

3500

(d) −7

−5 0 500 x 10−7 2

c5

3000

Frequency (Hz)

(c)

D6

2500

1334

Frequency (Hz)

Amplitude

2000

0

−5 0

Samples

x 10−4 1

2000

−6.0 4000

1000

−1 0 −7 500 x 10 5

Amplitude

D6

−5 0 500 x 10−7 5

4000 4 x 10 −5.6

1310

AICc5

1500

c5

1000

D6

−5 0 −7 500 x 10 5

x 10−6 1 #5 0 −1 0 −6 500 x 10 1

AICD6

1221

AICc5

0

x 10 #2

AICD6

5

Amplitude

−7

c5

Amplitude

(a)

c5

D6

2.5

20

40

60

80

100 120 140 160 180 200

Frequency (Hz)

0

0

20

40

60

80

100 120 140 160 180 200

Frequency (Hz)

Fig. 7. Six typical examples selected from Table 3 and their waveforms (Top), c5-AIC and D6 -AIC time series (Middle), and amplitude spectra of waveforms: c5 (red lines) and D6 (green lines) (Bottom). The vertical and vertical-dashed lines correspond to manual pickings and AIC pickings, respectively, and the numbers indicated by an arrow correspond to the P phase arrival points. The numbers in the left corner of the waveforms represent the case numbers indicated in Table 3. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

64

Computers & Geosciences 100 (2017) 57–66

X. Li et al.

Der, Z.A., Shumway, R.H., 1999. Phase onset time estimation at regional distances using the CUSUM algorithm. Phys. Earth Planet. Inter. 113, 227–246. http://dx.doi.org/ 10.1016/S0031-9201(99)00053-9. Earle, P.S., Shearer, P.M., 1994. Characterization of global seismograms using an automatic-picking algorithm. Bull. Seismol. Soc. Am. 84, 366–376. Galiana-Merino, J.J., Rosa-Herranz, J.L., Parolai, S., 2008. Seismic P phase picking using a Kurtosis-based criterion in the stationary wavelet domain. Geosci. Remote Sens. IEEE Trans. on 46, 3815–3826. http://dx.doi.org/10.1109/ TGRS.2008.2002647. García, L., Álvarez, I., Benítez, C., Titos, M., Bueno, Á., Mota, S., de la Torre, Á., Segura, J.C., Alguacil, G., Díaz-Moreno, A., Prudencio, J., García-Yeguas, A., Ibáñez, J.M., Zuccarello, L., Cocina, O., Patanè, D., 2016. Advances on the automatic estimation of the P-wave onset time. Ann. Geophys. 59, S0434. http://dx.doi.org/10.4401/ag7087. Ge, M.C., 2005. Efficient mine microseismic monitoring. Int. J. Coal Geol. 64, 44–56. http://dx.doi.org/10.1016/j.coal.2005.03.004. Gentili, S., Michelini, A., 2006. Automatic picking of P and S phases using a neural tree. J. Seismol. 10, 39–63. http://dx.doi.org/10.1016/j.coal.2005.03.004. Gou, X.T., Li, Z.M., Qin, N., Jin, W.D., 2011. Adaptive picking of microseismic event arrival using a power spectrum envelope. Comput. Geosci. 37, 158–164. http:// dx.doi.org/10.1016/j.cageo.2010.05.022. Hafez, A.G., Khan, M.T.A., Kohda, T., 2010. Clear P-wave arrival of weak events and automatic onset determination using wavelet filter banks. Digit. Signal Process. 20, 715–723. http://dx.doi.org/10.1016/j.dsp.2009.10.002. Hafez, A.G., Khan, T.A., Kohda, T., 2009. Earthquake onset detection using spectro-ratio on multi-threshold time-frequency sub-band. Digit. Signal Process. 19, 118–126. http://dx.doi.org/10.1016/j.dsp.2008.08.003. Hafez, A.G., Rabie, M., Kohda, T., 2013. Seismic noise study for accurate P-wave arrival detection via MODWT. Comput. Geosci. 54, 148–159. http://dx.doi.org/10.1016/ j.cageo.2012.12.002. Hildyard, M.W., Nippress, S.E., Rietbrock, A., 2008. Event detection and phase picking using a time-domain estimate of predominate period Tpd. Bull. Seismol. Soc. Am. 98, 3025–3032. http://dx.doi.org/10.1785/0120070272. Hildyard, M.W., Rietbrock, A., 2010. Tpd, a damped predominant period function with improvements for magnitude estimation. Bull. Seismol. Soc. Am. 100, 684–698. http://dx.doi.org/10.1785/0120080368. Huang, N.E., Shen, Z., Long, S.R., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. R. Soc. 454, 903–995. http://dx.doi.org/10.1098/rspa.1998.0193. Huang, N.E., Shen, Z., Long, S.R., 1999. A new view of nonlinear water waves: the Hilbert spectrum. Annu. Rev. Fluid Mech. 31, 417–457. http://dx.doi.org/10.1146/ annurev.fluid.31.1.417. Kalkan, E., 2016. An automatic P-phase arrival-time picker. Bull. Seismol. Soc. Am. 106, 971–986. http://dx.doi.org/10.1785/0120150111. Karamzadeh, N., Doloei, G.J., Reza, A.M., 2013. Automatic earthquake signal onset picking based on the continuous wavelet transform. IEEE Trans. Geosci. Remote Sens. 51, 2666–2674. http://dx.doi.org/10.1109/TGRS.2012.2213824. Kulesh, M., Diallo, M.S., Holschneider, M., Kurennaya, K., Kruger, F., Ohrnberger, M., Scherbaum, E., 2007. Polarization analysis in the wavelet domain based on the adaptive covariance method. Geophys. J. Int. 170, 667–678. http://dx.doi.org/ 10.1111/j.1365-246X.2007.03417.x. Küperkoch, L., Meier, T., Lee, J., Friederich, W., 2010. Automated determination of Pphase arrival times at regional and local distances using higher order statistics. Geophys. J. Int. 181, 1159–1170. http://dx.doi.org/10.1111/j.1365246X.2010.04570.x. Li, X.B., Shang, X.Y., Wang, Z.W., Dong, L.J., Weng, L., 2016. Identifying P-phase arrivals with noise: an improved Kurtosis method based on DWT and STA/LTA. J. Appl. Geophys. 133, 50–61. http://dx.doi.org/10.1016/j.jappgeo.2016.07.022. Lockman, A.B., Allen, R.M., 2005. Single-station earthquake characterization for early warning. Bull. Seismol. Soc. Am. 95, 2029–2039. http://dx.doi.org/10.1785/ 0120040241. Maeda, N., 1985. A method for reading and checking phase times in auto-processing system of seismic wave data. Zisin= Jishin 38, 365–379. http://dx.doi.org/10.1785/ 0120040241. Magotra, N., Ahmed, N., Chael, E., 1987. Seismic event detection and source location using single-station (three-component) data. Bull. Seismol. Soc. Am. 77, 958–971. Mallat, S.G., 1989. A theory for multiresolution signal decomposition: the wavelet representation. Pattern Anal. Mach. Intell. IEEE Trans. on 11, 674–693. Roberts, R.G., Christoffersson, A., Cassidy, F., 1989. Real-time event detection, phase identification and source location estimation using single station three-component seismic data. Geophys. J. Int. 97, 471–480. http://dx.doi.org/10.1111/j.1365246X.1989.tb00517.x. Romero, J.E., Titos, M., Bueno, A., Alvarez, I., Garcia, L., de la Torre, A., Benitez, M.C., 2016. APASVO: a free software tool for automatic P-phase picking and event detection in seismic traces. Comput. Geosci. 90, 213–220. http://dx.doi.org/ 10.1016/j.cageo.2016.02.004. Ross, Z.E., Ben-Zion, Y., 2014. Automatic picking of direct P, S seismic phases and fault zone head waves. Geophys. J. Int. 199, 368–381. http://dx.doi.org/10.1093/gji/ ggu267. Saragiotis, C.D., Hadjileontiadis, L.J., Panas, S.M., 2002. PAI-S/K: a robust automatic seismic P phase arrival identification scheme. Geosci. Remote Sens. IEEE Trans. on 40, 1395–1404. http://dx.doi.org/10.1109/TGRS.2002.800438. Saragiotis, C.D., Hadjileontiadis, L.J., Rekanos, I.T., Panas, S.M., 2004. Automatic P phase picking using maximum kurtosis and k-statistics criteria. Geosci. Remote Sens. Lett. IEEE 1, 147–151. http://dx.doi.org/10.1109/LGRS.2004.828915. Sedlak, P., Hirose, Y., Enoki, M., 2013. Acoustic emission localization in thin multi-layer

frequency components and a low temporal resolution for low frequency components ((Daubechies, 1992)). Nevertheless, the amplitude spectrum of micro-seism is mainly concentrated in low frequency bands. Besides, the main frequency is not always found in certain frequency bands due to the signal itself, the transmission media and the transmission distance (as can be observed in Figs. 3 and 7). Therefore, the D6 of the DWT-based picking method retains approximate P phase arrival information rather than an accurate P phase arrival. While EMD is based on the local features to decompose signals into several IMFs and the IMFs depend on the information of the signal itself, its decomposition is self-adaptive and adapts to non-linearity and non-stability signal analysis well. The random noise and other incoherent features will disappear in low IMFs and large IMFs, while the main IMFs retain P phase arrival well. Therefore, the c5 of EMD has better P phase arrival information than the D6 of DWT. From the above results and discussion, it can be concluded that the c5 of the EMDbased picking method provides effective means to identify P phase arrival. 7. Results In this paper, a new method for determining the P phase arrival of weak events has been presented. It is based on the EMD and the AIC picker. The performance of this algorithm has been tested by its application to 1938 micro-seismic signals from the Yongshaba mine (China). The DWT-based AIC picker has been selected for comparison. The P phase onsets have been compared to the onsets picked manually from the database. The results show that the EMD-AIC picking method works efficiently for the majority of P phase arrival identifications (even for low SNR micro-seismic signals) and is independent of the kind of noise. In addition, the EMD-AIC picking method has a better picking accuracy than the DWT-AIC picking method. In conclusion, the picker proposed is a useful P phase arrival picking method. Acknowledgements The authors gratefully acknowledge the financial supports of the National Natural Science Foundation of China (41630642 and 11472311) and National Key Research and Development Program of China (2016YFC0600706). The authors would also like to thank the micro-seismic project team members Chuxuan Zhang, Dong Liu and Yongyong Zhou for providing the analysis data picked manually. References Akaike, H., 1971. Autoregressive model fitting for control. Ann. Inst. Stat. Math. 23, 163–180. http://dx.doi.org/10.1007/978-1-4612-1694-0_12. Akaike, H., 1974. Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes. Ann. Inst. Stat. Math. 26, 363–387. http://dx.doi.org/10.1007/978-1-4612-1694-0_17. Akram, J., Eaton, D.W., 2016. A review and appraisal of arrival-time picking methods for downhole microseismic data. Geophysics 81, Ks71–Ks91. http://dx.doi.org/ 10.1190/Geo2014-0500.1. Allen, R.V., 1978. Automatic earthquake recognition and timing from single traces. Bull. Seismol. Soc. Am. 68, 1521–1532. Alvarez, I., Garcia, L., Mota, S., Cortes, G., Benitez, C., de la Torre, A., 2013. An automatic P-phase picking algorithm based on adaptive multiband processing. Geosci. Remote Sens. Lett., IEEE 10, 1488–1492. http://dx.doi.org/10.1109/LGRS.2013.2260720. Baer, M., Kradolfer, U., 1987. An automatic phase picker for local and teleseismic events. Bull. Seismol. Soc. Am. 77, 1437–1445. Baillard, C., Crawford, W.C., Ballu, V., Hibert, C., Mangeney, A., 2014. An automatic kurtosis‐based P‐and S‐phase picker designed for local seismic networks. Bull. Seismol. Soc. Am. 104, 394–409, (DOI: 0.1785/0120120347). Christoffersson, A., Husebye, E.S., Ingate, S.F., 1988. Wavefield decomposition using ML-probabilities in modelling single-site 3-component records. Geophys. J. Int. 93, 197–213. http://dx.doi.org/10.1111/j.1365-246X.1988.tb01996.x. Chu, C.-K.P., Mendel, J.M., 1994. First break refraction event picking using fuzzy logic systems. Fuzzy Syst., IEEE Trans. on 2, 255–266. http://dx.doi.org/10.1109/ 91.324805. Daubechies, I., 1992. Ten lectures on wavelets, in CBMS-NSF Regional Conference Series in Applied Mathematics. SIAM, Philadelphia, PA, USA, p. 357. DOI: http://dx.doi. org/10.1137/1.9781611970104.

65

Computers & Geosciences 100 (2017) 57–66

X. Li et al. plates using first-arrival determination. Mech. Syst. Signal Process. 36, 0888–3270. http://dx.doi.org/10.1016/j.ymssp.2012.11.008. Sedlak, P., Hirose, Y., Khan, S.A., Enoki, M., Sikula, J., 2009. New automatic localization technique of acoustic emission signals in thin metal plates. Ultrasonics 49, 254–262. http://dx.doi.org/10.1016/j.ultras.2008.09.005. Sleeman, R., van Eck, T., 1999. Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings. Phys. Earth Planet. Inter. 113, 265–275. http://dx.doi.org/10.1016/S0031-9201(99)00007-2. Tselentis, G.A., Martakis, N., Paraskevopoulos, P., Lois, A., Sokos, E., 2012. Strategy for automated analysis of passive microseismic data based on S-transform, Otsu's thresholding, and higher order statistics. Geophysics 77, Ks43–Ks54. http:// dx.doi.org/10.1190/geo2011-0301.1. Vidale, J.E., 1986. Complex polarization analysis of particle motion. Bull. Seismol. Soc.

Am. 76, 1393–1405, (DOI: 10.1.1.577.1425). Withers, M., Aster, R., Young, C., Beiriger, J., Harris, M., Moore, S., Trujillo, J., 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection. Bull. Seismol. Soc. Am. 88, 95–106. Yue, B.B., Peng, Z.M., Zhang, Q.H., 2014. Seismic Wavelet Estimation Using Covariation Approach. IEEE Trans. Geosci. Remote Sens. 52, 7495–7503. http://dx.doi.org/ 10.1109/TGRS.2014.2313116. Zhang, H.J., Thurber, C., Rowe, C., 2003. Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings. Bull. Seismol. Soc. Am. 93, 1904–1912. http://dx.doi.org/10.1785/0120020241. Zhao, Y., Takano, K., 1999. An artificial neural network approach for broadband seismic phase picking. Bull. Seismol. Soc. Am. 89, 670–680.

66