Method for Detection of First Phase in Seismogram Based on Akaike Information Criterion (AIC) Function

Method for Detection of First Phase in Seismogram Based on Akaike Information Criterion (AIC) Function

Arabind Kumar Yadav

Research Scholar, Department of Geology and Geophysics.

Indian Institute of Technology Kharagpur.

Kharagpur, West Bengal India

Pankaj K. Mishra

Research Scholar, Department of Geology and Geophysics.

Indian Institute of technology Kharagpur.

Kharagpur, West Bengal, India

Abstract A method has been developed to detect and pick-up the first phase arrival in seismogram from the vertical component by using two properties of the Akaike Information Criterion (AIC) function in wavelet domain i. e. 1) the phase arrival point belongs to global minima in the AIC function, and

1. the gradient of incremental after the phase arrival is larger than that before. Present method has been tested on earthquake waveform data with varying signal to noise ratio and different selection window. Consideration of the additional property of the AIC function in present method increase the performance and outperform the presently available methods based on AIC function, that uses property of AIC function having global minimum at the phase arrival in seismogram.

   Keywords- Akaike Information Criterion, particle motion analysis, seismic phase, objective function.

   1. INTRODUCTION

      Accurate detection and pick-up the phase arrival in seismogram is importance for earthquake epicenter location, source characterization, and also essential for other applications such as passive seismic tomography [1]. The identification of seismic phase is traditionally done by experience seismologist. However in the era of digital when volumes of data are available manual phase selection is unreliable. Several methods have been proposed to detect seismic phase using properties of the seismic signal in various domain i. e. frequency, polarization [2-4], amplitude [5-11], and in combination of these [12]. Some of the current method includes artificial neural network [13-14], maximum likelihood method [15-16], fuzzy logic [17], and autoregressive technique [18-21].

      The autoregressive technique also been used to detect the first arrival in microseismic data [22-23] and first arrival time of acoustic emission [24]. Autoregressive based technique use Akaike Information Criterion (AIC) function which has been proposed by two authors: Sleeman and Eck [19], and Maeda [18]. Sleeman and Eck [19] used autoregressive model of a fixed order while Maeda [18] consider seismogram itself in calculate of the AIC function. Sleeman and Eck [19] decided the order for the autoregressive model by trial and error method.

      In AIC based methods the AIC function is computed from part of the seismogram selected by rectangular window. The AIC function has minimum value at the point of phase arrival,

      which is used for phase detection and pick-up. Moreover, in each windowed seismogram there will be a global minimum irrespective to whether true phase arrival is present or not. Moreover seismogram also includes noise with seismic signal which affects the characteristic of AIC function. Zhang et al.,

      [25] calculated AIC function into wavelet domain up-to three levels and use the consistency of the global minimum as the P phase arrival from the vertical component of ground motion. In present method two properties of the AIC function have been used in combination to detect and pick-up phase arrival using AIC function in wavelet domain.

   2. METHODOLOGY

      Various techniques have been proposed to detect and pick-up arrival time of important phases in the seismogram, implemented on single and/or three components recording. The method proposed in present study is an autoregressive type, exploiting the characteristics of AIC function for the phase detection. In AIC function the phase point belong to the global minimum. Zhang et al., (2003) calculated AIC function from windowed seismogram in wavelet domain up to three levels (or scale) using method given by Maeda (1985). There will be a global minim in each time window of the seismogram weather phase arrival is present or not. Therefore, in order to get true pick-up of the P wave/first phase arrival, the selection of data window is important.

      To overcome this problem, in the present method additional property of AIC function has been used with property of global minima at phase arrival. The two properties are: 1) phase arrival belongs to the global minimum in AIC function and, 2) incremental gradient in AIC after the global minimum is larger than that before.

   3. AIC FUNCTION

      In standard Autoregression (AR-AIC) technique, it is assumes that seismogram can be divided into locally stationary segments as an autoregressive process and intervals before and after the onset are two different stationary processes Sleeman and Eck [19]. And order of the AR model or its coefficient or both change at data point when next section of the data is distinct from the following section. AIC is generally used to determine order of the autoregressive model when fitting with time series data [26]. Two methods have been given to

      2,max

      calculate AIC function by Sleeman and Eck [19] and by Maeda [18]. Sleeman and Eck [19] considered fixed order of autoregressive model to and the AIC function given as,

      trial and error, and then AR coefficient can be determined by the Yuke-Walker equations [27].

      1,max

      AIC(k) (k M ) log( 2

      ) (N M k) log( 2

      ) C2

      (1)

      Maeda [18] calculates the AIC function directly from

      where M is order of an Autoregressive (AR) process fitting the

      seismogram, without using the AR coefficients and given as,

      1,max

      data, C2, is constant, and 2

      2

      and

      2,max

      are the variance

      AIC(k) k.log{var( x[1, k])]} (N k 1).log{var( x[k 1, N])}

      (2)

      of the seismogram in the two intervals not explained by AR process. The order of the AR process must be specified by var is the variance. At the point belonging to the phase arrival time variance for each of the two part of the seismogram will be small therefore AIC function will have minimum value i.e. the phase arrival time belong to the point at the global minimum in AIC function. However each section of seismogram will have global minimum therefore, it is necessary to choose the time window for data which includes phase arrival.

   4. PRESENT METHOD

      In present method AIC function has been calculated from windowed seismogram using the method given by Maeda [18] described by (2). Since earthquake is abrupt release of stress energy in form of acoustic wave, therefore on phase arrival, amplitude increases in seismogram rapidly then decays slowly with times. Part of the recorded waveform before phase arrival and afterward are statistically distinct and variance of the two parts is minimum. Therefore value of the AIC function is minimum at the time of the phase arrival and increment after the phase arrival is larger than that before.

      In the present study two properties of the AIC function has been used to detect the phase arrival in the seismogram i. e.1) phase arrival belong to global minimum in AIC has and, 2) increment after global minimum is steeper than that before.

      where, k ranges through the entire seismogram samples, and

   5. MAIN STEPS IN CALCULATION

      Pre-processing

      Preprocessing of waveform data includes two steps: i) the instrument response correction, to remove the effect of the seismograph, and ii) the baseline correction to make mean of the seismogramat zero. Signal to noise ratio (SNR) is energy of the signal with energy of the noise added.

      Step 1: Select the part of seismogram using rectangular window in time domain.

      Step 2: Transform the windowed seismogram into wavelet domain (Daubechies wavelet [28]).

      Step3: Compute the Akaike Information Criterion (AIC) function with the approximation coefficient of wavelet transform at level 1.

      Step 4: Gradient calculation: find the global minimum point in AIC function computed in step 4. Compute the slope of the gradient on both sides of global minimum with the data points using one dimensional linear fit as given in Fig. 1.

      Step 5: Decision making: if incremental gradient (the slope) for the succeeding point is (0.1 in the present study) greater than that for preceding points, point associated in the global minima of the AIC function belong to phase arrival time. Gradient difference (Gd) 0.1 choose by hit and trial method, however is dependent on SNR of the signal. With lower SNR the optimum value decreases, however very small value may not be able to detect the absence of the phase in the windowed seismogram.

      Fig. 1: Calculation of gradient at global minima in AIC function in wavelet domain at level 1. First row represent the approximate coefficient of windowed seismogram at level 1. Second row is AIC function calculated from the approximate coefficient. Third row show the gradient on both sides of the global

      minimum.

   6. RESULTS AND DISCUSSION

      Present method has been applied on vertical component of strong ground motion data with varying signal to noise ratio and different selection window and compared with other AIC based method given by Zhang et al., [25]). I have used rectangular window in time domain for selection of seismogram. Random noise of different magnitude have been added to the seismogram and given in term of signal to noise ratio (SNR) in db. Both the method has been applied of two waveform data from different region. In present method gradient difference () value 0.1 has been used

      Firstly, both the methods have been applied on windowed seismogram including phase arrival in selected part of the seismogram. Fig. [2A 2B] show the performance of the Zhang et al., [25]) and present method respectively, applied on waveform data from Japan earthquake of Magnitude Mw 5.0 in 2004 recorded at the station code FKSH01 at an epicenter distance of 94 km. As depicted in Fig. [2A 2B] both methods correctly pick-up the phase arrival point in the signal having signal to noise ratio (SNR) varying from 1.00 to 0.25 db (SNR

      1.0 represent the observed seismogram before addition of random noise).

      Fig. 2A. Detection of Phase arrival using global minimum in wavelet domain up-to three levels (Zhang et al,. [25]). First row represent the window selected part of vertical component of seismogram. Second to fourth row represent AIC function from the approximate coefficient at level 1 to 3. Green dot is representing correct picking of phase.

      Fig. 2B. Present method. First row is vertical component of seismogram. Second row the approximate coefficient in wavelet domain at level 1 of the windowed seismogram which include phase. Third row is AIC function calculated from approximate coefficient at level 1. Green dot in figure is representing correct picking of phase.

      3, Zhang et al., [25] did not detect phase when it end of the

      To test effect of phase position in windowed seismogram, both the methods have been applied on same data (Japan earthquake) with SNR 0.35 db having phase relatively close to one end of the window as depicted in Fig. 3. As shown in Fig.

      window while the present method pick-up successfully.

      Fig. 3. Effect of selection window on detection of phase arrival. (a) by method given by Zhang et al., [25]) (black vertical line represent the phase arrival), and (b) by Present method (dot represent phase arrival time).

      Fig. [4A 4B] showed that performance of the both algorithm when part of the seismogram selected by the window does not include the phase arrival. As it is depicted in Fig. [4A 4B]

      Zhang et al. [25] method wrongly detect phase in seismogram for signal with SNR 0.35 db, while present method did not detect any phase arrival.

      Fig. 4A. Application of method by Zhang et al,. [25]). First row represent the window selected part of vertical component of seismogram. Second to fourth row represent AIC function from the approximate coefficient at level 1 to 3. Red dot in figure is representing wrong picking of phase.

      Fig. 4B. Present method. First row is vertical component of seismogram. Second row the approximate coefficient in wavelet domain at level 1 of the windowed seismogram which did not include phase. Third row is AIC function calculated from approximate coefficient at level 1.

      Both the methods have been applied on other waveform data of Jammu-Kashmir and Himanchal Pradesh border earthquake of Magnitude Mw 5.8 in 2013 recorded at station code JMU at epicenter distance 80 km. When part of seismogram selected by window includes the phase arrival, both the present method and Zhang et al., [25]) method detect and pick-up the phase

      arrival in signal with SNR 0.35 db while with SNR value 0.25 only present method detect but accuracy decreases as depicted in Fig. [5A 5B].

      Fig. 5A. Detection of Phase arrival using global minimum in wavelet domain up-to three levels (Zhang et al,. [25]). First row represent the window selected part of vertical component of the seismogram. Second to fourth row represent AIC function from approximate coefficient at level 1 to 3. Green dot in figure is representing correct picking of phase.

      Fig. 5B. Present method. First row is vertical component of seismogram. Second row the approximate coefficient in wavelet domain at level 1 of the windowed seismogram which include phase. Third row is AIC function calculated from approximate coefficient at level 1. Green/Yellow dots are representing accurate/inaccurate picking of phase respectively.

      In the case when part of the seismogram selected by window did not include the phase arrival, present method correctly failed to detect any phase for SNR to a minimum 0.35 while

      Zhang et al., [25] method detecting phase arrival in seismogram as depicted in Fig. [6A 6B].

      Fig. 6A. Detection of Phase arrival using global minimum in wavelet domain up-to three levels (Zhang et al,. [25]). First row represent the window selected part of vertical component of seismogram. Second to fourth row represent approximate coefficient at level 1 to 3. Red dot in figure is representing wrong picking of phase.

      Fig. 6B. Present method. First row is vertical component of seismogram. Second row the approximate coefficient in wavelet domain at level 1 of the windowed seismogram which did not include phase. Third row is AIC function calculated from approximate coefficient at level 1.

      variance increment after the phase arrival. In the present

      Comparing performance of both the method on strong ground motion data from different region, it is observed that present method is less sensitive to position of phase arrival in windowed seismogram i. e. selection window, and successfully detect and pick-up the phase arrival in signal having signal to noise ratio greater than or equal to 0.35 db. Moreover present method overcame the limitation of Zhang et al., [25] method of having phase arrival in windowed seismogram, and successfully reject when phase in not present in the windowed seismogram.

      The AIC based methods use global minima property of AIC function at phase arrival. This is equivalent to applying the magnitude of statistical variable (variance) in phase arrival detection i. e. at the point of the phase arrival the two part of the signal are statistical distinct and have minimum variance, therefore has minimum value of AIC function at phase arrival. Slope of increment (gradient) before and after phase arrival is relate to the trend characteristic of statistical variable i. e.

      method combining the two properties of AIC function is equivalent to combining the magnitude and the trend characteristic of the statistical variable; therefore improve the performance of AIC function in detection and pick-up of the phase arrival.

   7. CONCLUSION

Akaike Information Criterion (AIC) based method has been given for detection and pick-up the first phase arrival in seismogram using two characteristic of the AIC function calculated from seismogram in wavelet domain i. e. phase arrival time/point belong to global minimum in the AIC function and incremental gradient after the minima is greater than that before. Present method has been applied on strong ground motion waveform data from earthquake in different regions with varying signal to noise ratio (SNR) and selection window, and the results are compared with other AIC based method. It has been established that present method is less sensitive with selection window and correctly pick-up the

phase arrival for signal with SNR minimum up to 0.35 db. Also the present method successfully rejects when windowed seismogram did not include phase arrival. In present method combining the additional properties of the AIC function improve the performance to detect and pick-up the phase arrival in the windowed seismogram and outperform the presently available AIC based method that use the global minimum property of the AIC function.

ACKNOWLEDGEMENT

I would like to thank Mr. Hamid Sana to provide me waveform data of Jammu-Kashmir and Himanchal Pradesh border earthquake of Magnitude Mw 5.8 in 2013 recorded at station code JMU and National Research Institute for Earth Science and Disaster Prevention Japan for making the strong motion waveform data of earthquake in Japan available online for research.

REFERENCES

1. Tselentis G-A, N Martakis, P Paraskevopoulos, and Lois 2011b, High Resolution passive seismic tomography (PST) for 3D velocity, poisson ratio v and P- wave quality QP in the Delvina hydrocarbon field, southern Albanai, Geophysics, 76 no. 3, B1-B24.

2. Vidale J E, 1986, Complex polarization analysis of particle motion, Bulletin of the Seismological Society of America, 76, 13931405.

3. DAuria L, F. Giudicepietro, M. Martini, M. Orazi, R. Peluso, and G. Scarpato Polarization Analysis in the Discrete Wavelet Domain 2010, An Application to Volcano Seismology Bulletin of the Seismological Society of America April 100 :670-683.

4. Lois A, Sokos E, Martakis N, Paraskevopoulos P, and Tselentis 2013, G-A A new automatic S- onset detection technique: Application in local earthquake data, Geophysics, Vol. 78, No. 1 P. Ks1Ks11,9 Figs.,1Table.10.1190/Geo 201205 01

5. Allen R V, 1978, Automatic earthquake recognition and timing from signal traces: Bulletin of the Seismological Society of America, 68, 152132.

6. Allen R, 1982, Automatic Phase Pickers: Their Present Use and Future Prospects, Bulletin of the Seismological Society of America, Vol. 72, No. 68, pp. S225-S242.

7. Baer M, and Kradolfer U, 1987, An automatic phase picker for local and teleseismic events: Bulletin of the Seismological Society of America, 77, 14371445.

8. Earle P S and Shearer P M ,1994, Characterization of global seismograms using an automatic-picking algorithm, Bulletin of the Seismological Society of America, vol. 84, no. 2, pp. 366376, Apr.

9. Montalbetti J F, and Kanasewich E R, 1970, Enhancement of teleseismic body phases with a polarization filter: Geophysical Journal of the Royal Astronomical Society, 21, 119129, doi: 10.1111/j.1365- 246X .1970.tb01771.x.

10. Magotra N, Ahmed N, and Chael E, 1987, Seismic event detection and source location: Bulletin of the Seismological Society of America, 77, 95871.

11. Ruud B O, and Husebye E S, 1992, A new three-component detector and automatic single station bulletin production: Bulletin of the Seismological Society of America, 82, 221237.

12. Bai C Y, and Kennett B L N, 2000, Automatic phase-detection and identification by full use of a single three-component broadband seismogram: Bulletin of the Seismological Society of America, 90, 187198, doi: 10 .1785/0119990070.

13. Zao Y and Takano K, 1999, An articial neural network approach for broadband seismic phase picking, Bulletin of the Seismological Society of America, vol. 89, no. 3,pp. 670680.

14. Wang J and Teng T, 1995, Articial neural network-based seismic detector, Bulletin of the Seismological Society of America, vol. 85, no. 1, pp. 308319.

15. Robert R G, Christoffersson A and Cassidy F 1989, Real-time event detection, phase identication and source location estimation using

    single station three-component seismic data, Geophys. J., vol. 97, pp. 471480.

16. Christoffersson A E, Husebye S and Ingate S F, 1988, Wavelet decomposition using ML-probabilities in modeling single-site 3- component records, Geophys. J., vol. 93, pp. 197213.

17. Chu C K, and Mendel J, 1994, First break refraction event picking using fuzzy logic systems, IEEE Trans. Fuzzy Syst., vol. 2, no. 4, pp. 255266 November.

18. Maeda N, 1985, A method for reading and checking phase times in autoprocessing system of seismic wave data, Jishin, vol. 38, pp. 365379.

19. Sleeman R and Van Eck T, 1999, Robust automatic P-phase picking: an on-line implementation in the analysis of broad band seismogram-recordings, Phys. Earth Planet. Interiors 113, 265275.

20. Leonard M and Kennett B L N, 1999, Multi-component autoregressive techniques for the analysis of seismograms, Phys. Earth Planet.Interiors 113, 247264.

21. Leonard M, 2000, Comparison of manual and automatic onset time picking, Bulletin of the Seismological Society of America, 90, 13841390.

22. Onge A S, 2011, Akaike Information Criterion Applied to Detecting First Arrival Times on Microseismic

23. Han L,Wong J and Bancroft J C, 2010, Time picking on noisy microseismograms. GeoCanada.

24. Sedlak P, Yuichiro Hirose, Manabu Enoki, And Josef Sikula 2008, Arrival Time Detection In Thin Multilayer Plates, J. Acoustic Emission, 26.

25. Zhang H, Thurber C and Rowe C, 2003, Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single- component recordings, Bulletin of the Seismological Society of America, vol. 93, no. 5, pp. 19041912.

26. Akaike, H, 1973, Information theory and an extension of the maximum likelihood principle. In Proceedings of the Second International Symposium on Information Theory, B. N. Petrov and

    F. Csaki (eds), 267-281. Budapest: Akademiai Kiado.

27. Haykin S, 1996, Adaptive Filter Theory, 3rd ed. Englewood Cliffs, Prentice Hall, NJ.

28. Daubechies I, 1988, "Orthonormal Bases of Compactly Supported Wavelets." Comm. Pure Appl. Math. 41, 909-996.