Automatic Arrival Time Detection for Earthquakes

Automatic Arrival Time Detection for Earthquakes Based on Modified Laplacian of Gaussian Filter Omar M. Saada, Ahmed Shalabya,b, Lotfy Samyc, Mohammed S. Sayeda,d a

ECE Department, Egypt-Japan University of Science and Technology (E-JUST), Alexandria, Egypt. b Department of Computer Science, Faculty of Computers and Informatics, Benha University, Egypt. c National Research Institute of Astronomy and Geophysics, Helwan, Egypt. d ECE Department, Zagazig University, Zagazig, Egypt

Abstract— Precise identification of onset time for an earthquake is imperative in the right figuring of earthquake 's location and different parameters that are utilized for building seismic catalogues. P-wave arrival detection of weak events or microearthquakes cannot be precisely determined due to background noise. In this paper, we propose a novel approach based on Modified Laplacian of Gaussian (MLoG) filter to detect the onset time even in the presence of very weak signal-to-noise ratios (SNRs). The proposed algorithm utilizes a denoising-filter algorithm to smooth the background noise. In the proposed algorithm, we employ the MLoG mask to filter the seismic data. Afterward, we apply a Dual-threshold comparator to detect the onset time of the event. The results show that the proposed algorithm can detect the onset time for micro-earthquakes accurately, with SNR of −12 dB. The proposed algorithm achieves an onset time picking accuracy of 93% with a standard deviation error of 0.10 seconds for 407 field seismic waveforms. Also, we compare the results with short and long time average algorithm (STA/LTA) and the Akaike Information Criterion (AIC), and the proposed algorithm outperforms them. Keywords—Arrival time of earthquake (P-wave), Laplacian of Gaussian filter (LoG), Akaike Information Criterion (AIC), Automatic time picks, short and long time average (STA/LTA) algorithm.

I. INTRODUCTION Earthquakes are considered to be the major geohazards around the world. One critical safety issue is to safeguard people from the chemical and the nuclear radiation leaks induced by earthquakes. Therefore, the chemical sources, nuclear sources, and the petroleum transmission lines should be secured during the earthquake intervals. Earthquake Early Warning Systems (EEWSs) are used to alarm people exposed to the seismic threats. Ideally, this alarm should be released shortly after the earthquake occurrence, before the arrival of more destructive shear (S) and surface (Rayleigh and/or Love) waves. In the EEWS, the P-arrival time of the earthquake is necessary to be picked accurately to issue the alert. Therefore; real-time auto-picking algorithm with high accuracy is vital to support the fast decision in EEWS. One of the essential steps in the seismic analysis is detecting the onset time of the earthquake. According to the detection, earthquake location and several parameters are determined. However, the manual picking became more difficult with the increase of seismogram database size. Hence; automatic picking algorithms for the onset time of earthquakes become necessary. Also, one of the important concerns in picking the onset time is the existence of seismic noise. Micro-earthquakes data are often characterized by low signal-tonoise ratios (SNRs). For low SNR data, the manual picking of the onset time becomes more challenging and unreliable. In such environment, several false alarms may be produced or inaccurate onset time picking may occur because the discrimination between background noise and the seismic signal is problematic. Accordingly, many algorithms have been proposed to improve the accuracy of automatic picking. Short term averaging / long time averaging (STA/LTA) and the Akaike Information Criterion (AIC) have been widely used algorithms for event detection and picking algorithms [1]-[2]. STA/LTA algorithm is based on the identification of the abrupt change in characteristic functions. It has drawbacks due to the difficulty of setting up the parameters and the inaccuracy of the onset time detection, especially, with the existence of high background noise. For AIC method, the seismic data are modeled as an Autoregressive (AR) process where the minimum AIC value is set as the onset time of the seismic event. However, in some cases, the global minimum is determined wrongly due to the high noise level. Also, incorrect time window leads the onset time to be missed. Recently, a new method has been proposed to deal with low SNR seismic data based on the Fuzzy C-Means (FCM) clustering [3]. FCM depends on clustering the seismic data to seismic signal and seismic noise groups by obtaining the membership degree matrix. Data points which have high membership degree matrix are assigned to the seismic signal group. The initial time of the seismic signal group is set as the onset time of the seismic event. Moreover, some of the researchers addressed the onset time picking by applying the spectral ratio on time-frequency sub-band [4]. However, this algorithm was replaced by employing the wavelet filter bank on the same dataset and 1

achieved better accuracy in [5]. Using the same dataset, the accuracy of the onset time picking hit the peak using maximal overlap discrete wavelet transform [6]. Also, in [7]-[10], the authors proposed an automatic phase picker based on wavelet transformation. Furthermore, unsupervised machine learning technique was proposed to detect the arrival time in a low SNR environment [11]. The authors proposed to use Fuzzy clustering algorithm as an intelligent arrival picking algorithm [11]. Meanwhile, in [12], the authors proposed an algorithm to obtain the kurtosis value of small time windows, and if this value exceeded a pre-defined threshold, it would consider as the P-arrival time. In [13], the authors proposed automatic picking algorithm to pick the P-phase arrival based on the energy-frequency by obtaining the local maxima (LM) in earthquake seismograms. In [14], an artificial neural network was proposed to detect the P- and S-phases. Meanwhile, the presence of the arrival time was obtained by lifting the seismogram into a highdimensional space, and graph-based methods were used [15]. On the other hand, there are several seismic denoising algorithms were proposed to deal with seismic data. Chen et al. proposed a novel approach for denoising based on f-x Empirical-Mode Decomposition (EMD) predictive filtering [16]. Meanwhile, the conventional denoising operator with applying a weighting operator was proposed to predict the signal-leakage energy and retrieve it from the noise section [17]. A novel algorithm based on modified Multichannel Singular Spectrum Analysis (MSSA) was proposed to attenuate random noise with a new formula of low-rank reduction, which was named as the damped MSSA algorithm [18]. Also, an Ensemble Empirical Mode Decomposition (EEMD) combined with adaptive thresholding was proposed for microseismic denoising [19]. In [20], a multi-scale morphological method was proposed to detect the arrival time for weak micro-seismic signals. The authors proposed to decompose the data into multi-scale components based on the mathematical morphology theory [20]. Moreover, Mathematical Morphological Filtering (MMF) method was proposed to suppress low-frequency noise and preserve the low-frequency components of the signal at the same time [21]. The authors extracted empirical relation between morphological scale and frequency band so that the mathematical morphology-based approach can be conveniently used in low-frequency noise attenuation [21]. Huijian Li et al. proposed a mathematical morphology theory to suppress low-frequency noise in microseismic by discrimination between signals and background noise [22]. Meanwhile, Short Time Fourier Transform (STFT) with neighboring thresholding was proposed to reduce the background noise in microseismic data [23]. Other algorithms for denoising were proposed based on Time-Frequency domain with thresholding techniques to detect the micro-seismic event [24]-[25]. The Gaussian filter has been widely used in numerous image processing applications like the edge detection field [26]. The Laplacian of Gaussian (LoG) has gotten much attention since it was proposed by Marr [27] as a physiological model of the early human visual framework. In LoG technique, the image is filtered using the second derivative of the Gaussian mask to find zero-crossings in an image which represent the edges. The Gaussian mask is used to smooth the input signal, while the second derivative is obtained to enhance the image [28]. This algorithm is sensitive to sudden changes in the image hence; it can detect edges easily. In this paper, we propose a novel algorithm for picking the onset time based on filtering the seismic data using Modified Laplacian of Gaussian (MLoG) mask; then the filtered signal passes through a Dual-threshold system to detect the onset time. In the proposed algorithm, the background noise is well smoothed (denoised) such that the transition (P-arrival time) between the background noise and seismic signal can be detected using the Dual-threshold comparator. Two types of datasets are used to check the robustness of the proposed algorithm. First, synthetic signals with a different signal to noise ratios [3]. Second, field seismic data recorded by Egyptian National Seismic Network (ENSN) stations [4]-[6]. The given Receiver Operating Characteristic (ROC) curves and experimental results indicate that the proposed algorithm achieves higher accuracy compared to the common algorithms in the literature [1]-[6]. This paper is organized as follows: In Section II, the theory and the proposed algorithm for onset time detection are presented. Section III shows the experimental results. Section IV presents the implementation of the proposed algorithm on Field Programmable Gate Array (FPGA). Section V concludes this paper.

II. THE PROPOSED ALGORITHM In Image processing, the Laplacian of Gaussian (LoG) filter is used to detect the edges in the image. The edge can be defined as a sudden change in the image, and similarly, the P-arrival time in seismic data can be considered as a sudden change from background noise to the seismic signal. Therefore; we propose the Laplacian of Gaussian filter to detect the sudden change in seismic data by smoothing the background noise and make the P-arrival time more clear to be picked. The coefficients of LoG filter can be obtained by determining the second derivative of the Gaussian filter. First, Equation (1) shows the 1D Gaussian filter. 𝐺(𝑛) =

1 √2Π𝜎

.𝑒

−𝑛2 ) 2𝜎2

(

(1),

2

where 𝜎 represents the standard deviation of the Gaussian filter and n is the Gaussian index. Second, the first derivative of the Gaussian filter is shown in Equation (2).

𝐺

′ (𝑛)

=

−1

𝑛

.

√2Π 𝜎3

.𝑒

−𝑛2 ) 2𝜎2

(

(2)

Third, the LoG filter is the second derivative of the Gaussian filter, so that, the derivative of Equation (2) is obtained to reach the formula of LoG filter as shown in Equation (3). 𝐺 ′′ (𝑛) =

1 √2Π

.

1 𝜎2

.(

𝑛2 𝜎

2 − 1) . 𝑒

−𝑛2 ) 2𝜎2

(

(3)

The scaling factor in front of the exponential is primarily concerned with ensuring that the area underneath the Gaussian is one. Therefore, the normalized LoG filter can be reached by dividing Equation (3) by the exponential sum as shown in Equation (4).

𝐿𝑜𝐺 𝐹𝑖𝑙𝑡𝑒𝑟𝑁𝑜𝑟𝑚𝑎𝑙𝑖𝑧𝑒𝑑 =

1 1 . √2Π 𝜎2

−𝑛2 ( )

2

𝑛 .( 2 −1) . 𝑒 2𝜎2 𝜎

(4)

−𝑛2 ( ) ∑𝑛 𝑒 2𝜎2

The LoG filter can be considered as a denoising/FIR filtering technique for seismic data. For FIR filter, when the sum of the coefficients of the filter tends to zero, it acts as a high-pass filter. One of the characteristics of high-pass FIR filter that it has no response at zero frequency. Note that, a perfect edge detection filter which looks for sudden change must be a high-pass FIR filter which taps are sum to zero. Therefore; we modify the formula of LoG filter such that the sum of taps is zero as shown in Equation (5). We name the final formula for LoG filter as Modified Laplacian of Gaussian (MLoG). The MLoG filter smooths the background noise such that the transition between the background noise and the signal becomes more clear and sharp to be picked.

𝑀𝑜𝑑𝑖𝑓𝑖𝑒𝑑 𝐿𝑜𝐺 𝐹𝑖𝑙𝑡𝑒𝑟 =

1 1 . √2Π 𝜎2

−𝑛2 ( )

2

𝑛 .( 2 −1) . 𝑒 2𝜎2 𝜎 −𝑛2 ( ) ∑𝑛 𝑒 2𝜎2

1

- 𝐶 ∑𝑛

1 1 . √2Π 𝜎2

2

−𝑛2 ( )

𝑛 .( 2 −1) . 𝑒 2𝜎2 𝜎 −𝑛2 ( ) ∑𝑛 𝑒 2𝜎2

(5)

In Equation (5), C represents the number of taps (filter order or the number of filter coefficients) of the MLoG filter. By applying such a modification, the SNR of the output signal is enhanced by 4 dB. The order of MLoG filter and 𝜎 are set to be 10 and 2.5, respectively. The selection of 𝜎 and the order of MLoG filter is illustrated in section III. The MLoG filter is used to filter the seismic data. Then, the filtered signal is squared to maximize the difference between the seismic signal and the background noise. Afterward, the average over a sliding window with the length of m samples is calculated for smoothing as shown in Equation (6).

𝑊(𝑖) =

1 𝑚

∑𝑖+𝑚−1 𝑌1 (𝑖) 𝑖

(6)

In Equation (6), 𝑌1 represents the squared filtered signal using MLoG filter while m represents the length of the sliding window. After analysis, eight samples per window gives the best onset time accuracy result. Thus we set the value of m to eight. The selection of window length (m) is illustrated in section III. Finally, we set two thresholds to detect the onset time. The first one differentiates between seismic signal and background noise, while the second one is to guarantee that the sample is the onset time. In details, the average of a moving window of eight samples is calculated. Then, comparators are applied. Every two adjacent windows (W(i) and W(i-1)) must exceed the first threshold. The second window (W(i)) must be greater than the first window (W(i-1)) by a factor of the second threshold. If these two conditions met, the time of the first sample in the second window would be the onset time of the event. The selection of the two thresholds is illustrated in section III. Table I represents the pseudocode for the proposed algorithm, and a flowchart for the proposed algorithm is shown in Fig. 1. According to Fig.1 and Table I, the input seismic signal is stored. Afterward, the coefficients of the MLoG filter construct the MLoG mask according to a certain value of σ and filter order which tends to a discrete filter. Then, the stored seismic signal passes through the MLoG mask for denoising process. The output signal from MLoG filter is squared, and finally, a Dual-threshold comparator is applied to detect the arrival time of the event.

3

Reading Seismic Samples

Filtering Using MLoG MASK

Creating MLoG MASK

Squaring

Sliding windows (8samples per each) Obtaining the average for each window

Delay

W(i)

Decision is ZERO

NO

IF W(i) >Threshold 1 And IF W(i-1) >Threshold 1

W(i-1)

YES

NO

IF W(i) >W(i-1)*Threshold 2

YES

Fig.1. Flowchart for the Arrival time detection proposed algorithm.

4

Decision is ONE

Table I. PSEUDO CODE FOR THE PROPOSED ALGORITHM. Set Threshold1 and Threshold2; //( Symmetric FIR Filter 10 Taps (index number (n) ) ) Set n to [-5 -4 -3 -2 -1 1 2 3 4 5]; // (Initialization for FIR Registers) Set R[1 to 10] to zero; // (C[1 to 10] Set the Coefficients of MLoG Filter; // Reading Seismic event (x) Read (x); // (MLoG Filter) For II=1 to length(x) // Shift Register (R) For index=10 to 1 IF index =1 Then R[index] ← x(II); ELSE R[index] ← R[index-1]; End IF End FOR // Computing the output from MLoG filter Y(II) ← C[1]*R[1]+C[2]*R[2]+C[3]*R[3]+C[4]*R[4]+C[5]*R[5]+ C[6]*R[6]+C[7]*R[7]+C[8]*R[8]+C[9]*R[9]+C[10]*R[10]; End FOR // (Squaring the output from MLoG Filter) Y1 ← Y^2; // (Smoothing by taking the mean for each 8 samples) // Initialization (W = Mean (Y1 [1 to 8]) ) W(1) ← (Y1[1] + Y1[2] + Y1[3] + Y1[4] + Y1[5] + Y1[6] + Y1[7] + Y1[8])/8 ; // (Threshold Step) For I=2 to length(Y1) // (Smoothing the squared filtered signal by taking the mean for each 8 samples) ( Mean(Y1[I to I+7])) W(I) ← (Y1[I] + Y1[I+1] + Y1[I+2] + Y1[I+3] + Y1[I+4] + Y1[I+5] + Y1[I+6] + Y1[I+7])/8; IF( W(I) > Threshold1 & W(I-1)> Threshold1 & W(I) > W(I-1) * Threshold2) Then Decision ← 1; Break; End IF I ← I+1; End FOR

III.

EXPERIMENTAL RESULTS AND ANALYSIS

A. Analysis and Test of Proposed Algorithm The proposed algorithm is evaluated using synthetic datasets with different SNRs [3]. In addition, the proposed algorithm is applied on two real datasets; each one has local earthquakes with different SNRs and different magnitude. One of these datasets were used in [4]-[6]. For a fair comparison, the same environment of [3] is established but with our parameters. We want to guarantee that the

5

proposed algorithm is non-sensitive to the noise especially if the event is a micro-earthquake with a small magnitude which can be affected by the noise. Eleven synthetic models are constructed to obtain the best setting parameters of the proposed algorithm; each model has 1000 Ricker wavelets with variously dominated frequencies of 1, 5, 10, 15, 20 and 25 Hz and various amplitudes. Each model has different SNR, most of the possible cases are tested, and the SNR of the 11 models varies from -10 dB to 0 dB. Hence; we test 11000 synthetic signals (11 models * 1000 Ricker wavelets). Various sigma values are attempted starting from 1 to 10 with a step size of 0.05, and several filter orders are tested starting from 1 to 20 with a step size of 1. According to the analysis, the MLoG filter is more sensitive to the seismic signal when the filter order is greater than or equal 10 with 𝜎 of 2.5. Therefore; we set the MLoG filter order and 𝜎 to be 10 and 2.5, respectively. For example, Fig.2(a) shows three Ricker wavelet with dominated frequencies of 1, 5 and 10 Hz and the sampling frequency is 100 Hz. A white Gaussian noise (WGN) is added to the synthetic signal with SNR of -8 dB as shown in Fig. 2(b). In this paper, we define the SNR as shown in Equation (7) [3].

𝑆𝑁𝑅 = 10 𝑙𝑜𝑔10 ∑

∑𝑡 |𝑥(𝑡)|2

(7),

2 𝑡 |𝑠(𝑡)−𝑥(𝑡)|

where, x(t) is the seismic signal and, s(t) is the noise signal (background noise). Noise signal

Ricker Wavelet

sigma = 2.500000 , order = 5

1

1

1

(a)

(b) 0.8

0

0.5

Amplitude

Amplitude

Amplitude

0.5

0

0.6 0.4

-0.5

0.2

2

4

-0.5 0

6

Time (s)

sigma = 10.000000 , order = 5

Amplitude

(e)

4 Time (s) sigma = 2.500000 , order = 15

6

1

0.4

2


0.2

6

(h)

0.6 0.4 0.2

2

4 Time (s)

6

0.4

0 0 3

0.8 Amplitude

0.4

0.6

x 10

2

4

6

Time (s)

-3

2

(g)

0.6

(f)

0.2

1

0.8

6

0.8

0.2

2

4 Time (s)

0.6

0 0

2

sigma = 10.000000 , order = 10

0 0

MloG Mask Values

Amplitude

0.4 0.2

Amplitude

0 0

6

1

0.8

(d)

0.6

0 0


1

1 0.8

2

Amplitude

-1 0

0 0

(c)

1 0 -1

(I)

-2 -3 -4 1

2

4

6

2

3

4

5 6 Index (n)

7

8

Time (s)

Fig.2. (a) Ideal signal. (b) Noisy signal (SNR=-8 dB). (c) – (h) output signal from our proposed algorithm with different choices of 𝜎 (2.5 and 10) and different choices of filter order (5, 10, and 15),and (I) the MLoG Mask with filter order and 𝜎 of 10 and 2.5, respectively.

6

9

10

Fig.2 (c)-(h) show the output signal of the proposed algorithm using various 𝜎 and different MLoG filter orders. Fig. 2 (c) and (d) show the effect of using 5th order filter with different 𝜎 of 2.5 and 10, respectively. Meanwhile, Fig.2 (e) and (f) show the influence of using 10th order filter with different 𝜎 of 2.5 and 10, respectively. Finally, Fig.2 (g) and (h) show the effect of using 15th order filter with various 𝜎 of 2.5 and 10, respectively. Fig.2 (I) shows the proposed MLoG mask with filter order and 𝜎 of 10 and 2.5, respectively. We plot ROC curves for the signals with various SNRs. ROC curve plots the false positive rate versus the true positive rate. The two thresholds in the proposed algorithm are tuned to reach the maximum Area Under the ROC Curves (AUC) value which indicates the accuracy of the proposed algorithm. Higher AUC value indicates a higher accuracy performance of the proposed algorithm. We set two classes for ROC curves, class one if the algorithm detects a seismic event and class zero if no seismic signal found (noise). Two datasets are created to evaluate the proposed algorithm. First, 1000 Ricker wavelets with variously dominated frequencies of 1, 5, 10, 15, 20 and 25 Hz are created, and a WGN is added with different SNRs. The second dataset contains 1000 Ricker wavelets with the same dominated frequencies of the first dataset but, real noise is added with various SNRs. In the proposed algorithm, the two thresholds which maximize the AUC values for the two datasets are 0.2 and 2, respectively. Table II shows the accurate picking, inaccurate picking, AUC values and right seismic signals picked within 0.02 seconds for the previous two datasets. While, Table III shows the comparison between the proposed algorithm, STA/LTA, AIC, and FCM algorithms when the SNR is -5 dB. Fig.3 (a) and (b) show the ROC curves with different SNRs reached by the WGN and real noise for the proposed algorithm, respectively. While Fig.3 (c) and (d) show a comparison of the ROC curves between the proposed algorithm versus STA/LTA, AIC and FCM algorithms when the SNR is -5 dB. According to Table II, Table III, and Fig.3, we can conclude that the proposed algorithm outperforms the other three algorithms. The STA/LTA and AIC cannot pick accurately the onset time when the SNR is -5 dB. At SNR of -5 dB, the proposed algorithm achieves correct onset time accuracy of 99% and an AUC value of 0.9868 while, the STA/LTA, AIC and FCM algorithms produce picking accuracy of 79.1%, 75.5%, and 94.5%, respectively. Meanwhile, the STA/LTA, AIC, and FCM algorithms reach AUC values of 0.781, 0.728 and 0.9791, respectively. Furthermore, 99% of the correctly picked seismic signals which are detected using the proposed algorithm falls within the range of 0.02 seconds at SNR of -5 dB. Also, the proposed algorithm reaches an acceptable accurate onset time when the SNR is -12 dB while, FCM algorithm loses its ability to pick the onset time accurately when the SNR goes below -8 dB. The proposed algorithm reaches correct time picking accuracy of 97%, 84% and 81% for SNR of -8, -10 and -12 dB, respectively. Hence, the proposed algorithm is suitable when the SNR is higher than -12 dB. 1

0.8 0.6 0.4 0.2 0 0

0.2

0.4 0.6 False Positive Rate

(a)

0.8

0.8

0.4 0.2

1

1

0.8

0.8

MLOG FCM STA/LTA AIC

0.6 0.4 0.2 0 0

(c) 0.2


0.8

0.2


0.4

0.8

1

(d)

0.2 0.2


Fig.3. ROC curves for signals with (a) WGN and (b) Real noise for difeerent SNRs. (c) and (d) ROC for the four methods with WGN and Real noise, respectievly.

7

(b)

MLOG FCM STA/LTA AIC

0.6

0 0

1

snr=5dB snr=0dB snr=-5dB snr=-7dB snr=-8dB snr=-10dB snr=-12dB snr-14dB

0.6

0 0

1

True Positive Rate

True Positive Rate

snr=5dB snr=0dB snr=-5dB snr=-7dB snr=-8dB snr=-10dB snr=-12dB snr-14dB

True Positive Rate

True Positive Rate

1

0.8

1

Table II. AUC VALUES FOR THE TWO DATASETS WITH DIFFERENT SNRS USING PROPOSED ALGORITHM The WGN test

The real noise test

SNR (dB)

Number of Signals

Correct picking

Incorrect picking

Picking within 0.02 s

AUC values

Number of Signals

Correct picking

Incorrect picking


AUC values

5 0 -5 -7 -8 -10 -12 -14

1000 1000 1000 1000 1000 1000 1000 1000

1000 1000 999 996 983 851 820 753

0 0 1 4 17 149 180 247

1000 999 993 941 905 680 600 334

0.9999 0.9980 0.9871 0.9743 0.9458 0.8890 0.8221 0.7800

1000 1000 1000 1000 1000 1000 1000 1000

1000 999 995 990 978 840 810 748

0 1 5 10 22 160 190 252

999 995 991 938 900 672 580 320

0.9998 0.9978 0.9868 0.9738 0.9430 0.8871 0.8180 0.7770

Table III. AUC VALUES FOR THE TWO DATASETS FOR DIFFERENT ALGORITHMS. SNR = -5 dB Method

STA/LTA AIC FCM Proposed

The WGN test

The real noise Test

Number of Signals

Correct picking

Incorrect picking


AUC values

Number of Signals

Correct picking

Incorrect picking


AUC values

1000 1000 1000 1000

791 753 948 999

209 247 52 1

121 100 859 993

0.7950 0.7510 0.9822 0.9871

1000 1000 1000 1000

774 726 945 995

226 274 55 5

99 78 848 991

0.7810 0.7280 0.9791 0.9868

Output signal SNR = -10

4

6 Tim e (s) Ricker Wavelet SNR = -8

8

10

Amplitude

0.5 0 -0.5

0

2

4


8

10

Amplitude

0.5 0 -0.5

2

4


8

10

0 -1

0

2

4


8

0 -1

-0.5

0

2

4 6 Tim e (s) (a)

8

10

Amplitude

4

6 Tim e (s) Output signal SNR = -8

8

10

0

2

4


8

10

0.5 0 -0.5

10

1

2

0

Amplitude

1

0

0.5

Amplitude

Amplitude

0

-0.5

Amplitude

2

Amplitude

0

0

0

2

4


8

10

1 0 -1

0

2

4


8

0 -0.5

0

0

2

4 6 Tim e (s) (b)

8

10

2

2.2

2.4

2.6 2.8 3 Tim e (s) Output signal SNR = -8

3.2

3.4

3.6

3.8

4.6

4.8

5

5.2

5.4

6.2

6.4

-1 6.4

6.6

6.8

7.6

7.8

8.6

8.8

9.6

9.8

0.5 0 -0.5 3.4

4 4.2 4.4 Tim e (s) Output signal SNR = -6

0.5 0 -0.5

10

1

-1

0.5

Amplitude

-0.5

Output signal SNR = -10

0.5

Amplitude

Amplitude

0

Amplitude

Amplitude

Amplitude

Amplitude

Ricker Wavelet SNR = -10 0.5

5.6 5.8 6 Tim e (s) Output signal SNR = -4

1 0

7 7.2 7.4 Tim e (s) Output signal SNR = -2

1 0 -1 8.4

9

9.2

9.4

Time (s) (c)

Fig.4. (a) Syntheticseismic signal with various SNR, (b) output signal from the proposed algorithm for synthetic seismic signal, and (c) 1-second zoom for Parival time for synthetic seismic signal

8

Fig. 4(a) shows an example of different synthetic events (5 events) with SNR range of -9 dB to -1 dB. The proposed algorithm success to detect the onset time accurately as shown in Fig.4(b), despite the existence of high background noise. Fig.4(c) shows a 1-second zoom for the P-arrival time. The black stars in Fig.4 represent the manual picking while the red circles represent the onset time according to the proposed algorithm. Moreover, in order to find the optimal parameters for the proposed algorithm, we optimize the filter order, 𝜎, first threshold, second threshold and, the number of samples per window (m) using Particle Swarm Optimization (PSO) algorithm [29]. Particle Swarm Optimization (PSO) has been widely used in many applications due to its outstanding advantages including its simplicity and easiness of implementation [30]-[32]. The basic PSO algorithm consists of three steps generating positions of particles and velocities, velocity update, and position update [33]. Each particle represents a possible solution to the problem that changes its position from one iteration to another based on velocity updates. First, the positions, xi, and velocities, vi, of the initial swarm of particles are randomly generated. The PSO consists of many particles which form a swarm, design space. At each step, each particle updates its velocity and distance according to Equations (8) and (9), respectively [30]. 𝑣𝑖= 𝑣𝑖 + 𝑐1 ∗ 𝑟𝑎𝑛𝑑 ∗ (𝑃𝑖 − 𝑥𝑖 ) + 𝑐2 ∗ 𝑟𝑎𝑛𝑑 ∗ (𝑃𝑔𝑙𝑜𝑏𝑎𝑙 − 𝑥𝑖 ) xi= xi + vi

(8) (9)

In Equation (8), the Pi represents the best previous position and the global best position stored in Pglobal. The velocity update formula includes some random parameters (rand) to ensure excellent coverage of the design space. These random parameters are uniformly distributed. According to [29], the original PSO algorithm used the value of two for both constants C1 and C2. For PSO, several models are tested to reach the optimum parameters for the proposed algorithm. Eleven different models are tested; each model has 1000 Ricker wavelets with variously dominated frequencies of 1, 5, 10, 15, 20 and 25 Hz and various amplitudes. The 11 models have different SNRs between ranges of -10 dB to 0 dB. Hence; we test 11000 synthetic signals (11 models * 1000 Ricker wavelets). The cost function is minimizing the time error between the automatic picking according to the proposed algorithm and the manual picking. In the first iteration, random initialization takes place for the five parameters, and then the cost function is obtained. Afterward, for each iteration, the five parameters are updated according to PSO scheme, then the newly updated parameters are used to obtain the cost function. If the cost function is less than the global minimum, the current parameters will be considered as the optimum parameters. Meanwhile, the cost function will reach zero if the difference between manual picking and automatic picking according to the proposed algorithm is zero (Error = zero). If this condition happens, the PSO will be stopped, but if this condition did not meet, the process would continue till the number of iterations reaches 2000, and the program will be halted. Fig.5 shows the flowchart of PSO process. Finally, the optimum value of the five parameters obtained by PSO matched the value of the parameters obtained in section II and section III which cope with our analysis. Start

Set the PSO parameters Sigma, Filter order, threshold1, threshold2, and m Passing the 11 models (11000 Ricker wavelets) via the proposed algorithm using the previous parameters.

Update the new velocity and distance for each particle in the swarm

Determining the Cost Function =error Between manual Picking and automatic Picking for the 11000 Ricker Wavelets.

No Update the optimum parameters and the global minimum

If the cost function = 0 or the number of iterations =2000

No

If the Cost Function less than the global minimum

yes end

Fig.5 Flowchart for PSO test.

9

yes

B. Field Seismic Data The proposed algorithm is tested using the dataset in [4]–[6]. This dataset consists of 146 local earthquakes from Cairo region recorded by the Egyptian national seismic network (ENSN) and detected by the vertical component (Z-Component) of KOT, HAG, and FYM stations. For this dataset; earthquakes have different SNR to check the robustness of the proposed algorithm. Furthermore, to validate the reliability of the proposed algorithm, 261 waveforms are collected for local earthquakes recorded by KOT and HAG stations between 2014 and 2016. Both datasets include micro/local earthquakes with small/large magnitudes (ML) and low/high SNRs. For both datasets, the earthquake's magnitude varies from 0.5 to 5 ML, While, the range of SNR is -8 dB to 15 dB. Fig. 6 shows the location of these local earthquakes. The total waveforms of the two datasets are 407. For both datasets, the proposed algorithm is applied to every seismic waveform to detect the onset time then we compare the result with the manual picking. Table IV shows an accuracy comparison between the proposed algorithm and the common algorithms in the literature [1]-[6]. The proposed algorithm achieves correct picking time accuracy of 93% with standard deviation error of 0.10 seconds while, STA/LTA, AIC, and FCM algorithms reach picking accuracy of 76%, 72%, and 88%, respectively with standard deviation error of 0.11, 0.12 and 0.11 seconds. Also, the proposed algorithm outperforms the algorithms in [4]-[6] as shown in Table IV.

FYM

▲

Fig.6. Location of local earthquakes detected by KOT and HAG stations during the interval of 2014 and 2016

Table IV. ONSET TIME PICKING ACCURACY COMPARISON. Method [4] [5] [6] STA/LTA [1] AIC [2] FCM [3] Proposed

Number of waveforms

Correct waveforms Pick

Accuracy

within 0.05s

Standard Deviation(s)

407 407 407 407 407 407 407

321 335 357 310 294 359 379

79% 82% 87.5% 76% 72% 88% 93%

228 244 246 178 127 318 346

0.133 0.13 0.13 0.11 0.12 0.11 0.10

10

Fig. 7(a) shows an example of different micro-earthquake events (9 events) with SNR range of -8 dB to -4 dB, and a magnitude range of 0.5ML to 2ML. The proposed algorithm success to detect the onset time accurately despite, the existence of high background noise. Also, in Fig. 7(a), the proposed algorithm can detect multiple onset time in the same waveforms as seen in waveform number 3 and 4. Moreover, Fig. 7(d) shows a high SNR (>2 dB) local earthquake events (9 events) with a magnitude range of 2ML to 4ML for which the proposed algorithm detects the onset time easily. Meanwhile, Fig.7. (b) and (e) show the output signal after passing through the proposed algorithm. Fig.7 (c) and (f) show a 1-second zoom around the P-arrival time. According to Fig.7, the background noise is smoothed, while the transition between the background noise and the seismic signal becomes clear and sharp to be picked. The black stars in Fig.7 represent the manual picking while the red circles represent the onset time according to the proposed algorithm. According to these results, the proposed algorithm can detect the onset time with low and high SNR also; the proposed algorithm can deal with the micro-earthquakes under low SNR conditions. We obtain the Probability Distribution Function (PDF) of the error for the real seismic datasets for different algorithms as shown in Fig. 8. The PDF for the proposed algorithm, FCM, STA/LTA and AIC algorithms are shown in Fig. 8(a), Fig. 8(b), Fig. 8(c) and Fig. 8(d), respectively. According to Fig.8 and Table IV, more than 90% of the correct picking for the proposed algorithm falls within 0.05 seconds from the manual pick. The STA/LTA algorithm picks the P- wave time prematurely while the AIC algorithm picks the onset time belatedly. The correct picking for FCM, AIC, and STA/LTA algorithms are 88.5%, 43% and 57% within 0.05 seconds from the manual pick, respectively. Thus; the proposed algorithm outperforms the other three algorithms. The setting parameters for the proposed algorithm (filter order, σ, the number of samples per window (m), and the second threshold in the proposed algorithm) are the same as mentioned in section II. Notice that, the first threshold is adaptive to the input seismic data since the input field seismic data are recorded using a Digitizer with a sensitivity of 24 bits. Therefore; we obtain the first threshold by determining the average seismic noise for the first dataset [4]-[6] and Particle Swarm Optimization (PSO) is utilized to confirm this value.

I. IMPLEMENTATION OF THE PROPOSED ALGORITHM After the analysis, it is clear that the proposed algorithm can detect the onset time accurately with various SNRs levels and different magnitude while it reduces the false alarm. Also, the proposed algorithm is simple and can be easily implemented as a hardware module in an Early Earthquakes Warning System (EEWS). In this section, we report the computational complexity of the proposed algorithm on real hardware resources. The proposed algorithm is prototyped on Field Programmable Gate Array (FPGA). In [34], onset time picking module was previously implemented on FPGA to support the decision on EEWS. There are two basic techniques for EEWS, the single station (on-site) approach, and the network approach. For on-site warning system, the data are real-time processed at the station [35]. However, using the single station for onset time detection can produce several false alarms. On the other hand, the network approach uses several seismic stations to release the alarm [36]-[37]. This approach depends on sending the seismic data to the main site where the decision is taken centrally. In this approach, the false alarm is decreased while the system takes more time for releasing the alarm because it requires multiple stations to detect the earthquake. Therefore, our approach is to implement the algorithm on FPGA to reduce the processing time in both on-site system and the network approach. When the picking algorithm is implemented in each station, the decision for earthquake detection will be performed in the station, and then this information will be sent to the main site. Also, our proposed algorithm reduces the false alarm in the on-site warning system approach because it has high accuracy for picking the onset time. To implement the algorithm on FPGA, we apply the digital systems design flow: floating-point MATLAB, fixed-point MATLAB, Register Transfer Level (RTL) implementation using Verilog and synthesis over FPGA. We implement a simple FIR filter [38] to represent the MLoG mask which consists of ten coefficients. Furthermore, all the multipliers in the proposed algorithm are implemented using shift and add algorithm [39]. Finally, the proposed design is prototyped on a cheap FPGA kit that is called CMOD A7 [04] to be compatible with EEWS. Table V shows that the proposed design has fewer hardware resources than the common algorithms in the literature [1]–[6]. The STA/LTA algorithm has fewer resources rather than the proposed algorithm, but it fails to detect the onset time under certain conditions.

11

Waveforms Number Waveforms Number

Fig.7. (a) local micro earthquakes with low SNR, (b) output signal from the proposed algorithm for micro earthquakes, (c) 1-second zoom for P-arival time for microearthqaukes, (d) local earthquakes with high SNR (e) output signal from the proposed algorithm for local earthquakes, and (f) 1-second zoom for P-arival time for local earthquakes. Manual pick and picks accoridng to the proposed algorithm are black stars and red circles, respectievly.

12

The Probability Distribution Function

100

0 -0.1

-0.05 0 0.05 Discrepancy

0.1

(b)

Percentage (%)

50

Percentage (%)

Percentage (%)

(a)

50

0 -0.1



0.1

50

40

(c)

30 20 10 0 -0.1


0.1

(d)

40 30 20 10 0 -0.1

(s)

(s)

(s)


50 Percentage (%)


100


(s)

Fig. 8. Probability Distribution Function of error for (a) Proposed Algorithm, (b) FCM, (c) STA/LTA and (d) AIC

Table V. HARDWARE UTILIZED RESOURCES FOR DIFFERENT ALGORITHMS Hardware Resources

Slice LUTs

Slice Registers

CMOD available Resources [4] [5] [6]

20800 21500(> 100%) 14560(70%) 14560(70%)

41600 10100(24.3%) 2100(5%) 2100(5%)

FCM [3] AIC [2] STA/LTA [1]

4161(20%) 11856(57%) 2496(12%) 3016(14.5%)

1051(2.4%) 1372(3.3%) 500(1.2 %) 749(1.84%)

Proposed

II. CONCLUSION An automatic onset time picking algorithm is introduced based on Modified Laplacian of Gaussian (MLoG) filter. For image processing field, Laplacian of Gaussian (LoG) filter is widely used for edge detection. The Modified Laplacian of Gaussian filter (MLoG) is introduced in this paper to automatically detect the onset time of any seismic event despite the existence of high background noise. Firstly, the MLoG mask is utilized to filter the seismic data. Then, the output signal from the MLoG filter is squared and passes through a Dual-threshold module to ensure the onset time is picked accurately. The strength of the proposed algorithm lies behind its ability for smoothing the background noise using MLoG filter. ROC curves show that the proposed algorithm works perfectly with SNR higher than -12 dB. Furthermore, the proposed algorithm is evaluated using local earthquakes and reaches high picking time accuracy when compared to the common algorithms in the literature [1]-[6]. Many micro-earthquakes with low SNR are tested using the proposed algorithm and accurate onset time is obtained successfully. The computational complexity of the proposed algorithm is less than the common algorithms in the literature [1]-[6] except for STA/LTA [1]. A comparison of the computational complexity is demonstrated by prototyping the proposed algorithm on FPGA and estimating the required resources for implementing it. ACKNOWLEDGMENT We would like to thank Egypt-Japan University of Science and Technology (E-JUST) for the continuous support and the Egyptian Ministry of Higher Education for funding this work. Also, we would like to thank Dr. Ali Gamal (National Research Institute of Astronomy and Geophysics - NRIAG) for the suggestion of searching in auto-picking, earthquake early warning system, and for his support and advice.

13

0.1

REFERENCES [1]

R. V. Allen, “Automatic earthquake recognition and timing from single traces,” Bull. Seismol. Soc. Amer., vol. 68, no. 5, pp. 1521–1532, Oct. 1978.

[2]

M. Leonard and B. L. N. Kennett, “Dual-component autoregressive techniques for the analysis of seismograms,” Phys. Earth Planetary Interiors, vol. 113, nos. 1–4, pp. 247–263, Jun. 1999.

[3]

Zhu, Dan, Yue Li, and Chao Zhang. "Automatic Time Picking for Microseismic Data Based on a Fuzzy C-Means Clustering Algorithm." IEEE Geoscience and Remote Sensing Papers 13, no. 12 (2016): 1900-1904.

[4]

Hafez, Ali G., Tahir A. Khan, and Tohru Kohda. "Earthquake onset detection using spectro-ratio on dual-threshold time–frequency sub-band." Digital Signal Processing 19, no. 1 (2009): 118-126.

[5]

Hafez, Ali Gamal, Muhammad Tahir Abbas Khan, and Tohru Kohda. "Clear P-wave arrival of weak events and automatic onset determination using wavelet filter banks." Digital Signal Processing 20, no. 3 (2010): 715-723.

[6]

Hafez, Ali G., Mostafa Rabie, and Tohru Kohda. "Seismic noise study for accurate P-wave arrival detection via MODWT." Computers & Geosciences54 (2013): 148-159.

[7]

Bogiatzis, P., and M. Ishii, 2015, Continuous wavelet decomposition algorithms for automatic detection of compressional- and shear-wave arrival times Bull. Seismol. Soc. Am., 105, pp. 1628–1641

[8]

Mousavi, S. M., C. A. Langston, and S. P. Horton, 2016, Automatic microseismic denoising and onset detection using the synchrosqueezedcontinuous wavelet transform: Geophysics, 81, no. 4, V341 –V355.

[9]

Karamzadeh, N., G. J. Doloei, and A. Reza, 2013, Automatic earthquake signal onset picking based on the continuous wavelet transform: IEEE Transaction on Geoscience and Remote Sensing, 51, 2666–2674.

[10] Zhang, H., C. Thurber, and C. Rowe, 2003, Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings: Bulletin of Seismological Society of America, 93, 1904– 1912. [11] Chen, Yangkang. "Automatic microseismic event picking via unsupervised machine learning." Geophysical Journal International 212, no. 1 (2017): 88-102. [12] Baillard, C., W. C. Crawford, V. Ballu, C. Hibert, and A. Mangeney, 2014, An automatic kurtosis-based P- and S-phase picker designed for local seismic networks: Bulletin of Seismological Society of America, 104, 394–409. [13] Panagiotakis, C., E. Kokinou, and F. Vallianatos, 2008, Automatic P-phase picking based on local-maxima distribution: IEEE Transaction on Geoscience and Remote Sensing, 46, 2280–2287. [14] Gentili, S., and A. Michelini, 2006, Automatic picking of P- and S-phases using a neural tree: Journal of Seismology, 10, 39–63. [15] Taylor, K. M., M. J. Procopio, C. J. Young, and F. G. Meyer, 2011, Estimation of arrival times from seismic waves: A manifold-based approach: Geophysical Journal International, 185, 435–452. [16] Chen, Y., J. Ma, (2014), Random noise attenuation by f-x empirical mode decomposition predictive filtering, Geophysics, 79, pp. V81–V91 [17] Chen, Y., S., Fomel, (2016), Random noise attenuation using local signal-and-noise orthogonalization, Geophysics, 80 (2015), pp. WD1–WD9. [18] Huang,W., R. Wang, Y. Chen, H. Li, and S. Gan, 2016, Damped multichannel singular spectrum analysis for 3D random noise attenuation, Geophysics, vol. 81, no. 4, pp. V261–V270, 2016. [19] Han, J., and M. van der Baan, (2015). Microseismic and seismic denoising via ensemble empirical mode decomposition and adaptive thresholding, Geophysics, 80(6), KS69–KS80. [20] Li, Huijian, Runqiu Wang, Siyuan Cao, Yangkang Chen, Nan Tian, and Xiaoqing Chen. "Weak signal detection using multiscale morphology in microseismic monitoring." Journal of Applied Geophysics 133 (2016): 39-49. [21] Huang, Weilin, Runqiu Wang, Shaohuan Zu, and Yangkang Chen. "Low-frequency noise attenuation in seismic and microseismic data using mathematical morphological filtering." Geophysical Journal International 211, no. 3 (2017): 1318-1340. [22] Li, H., R. Wang, S. Cao, Y. Chen, W. Huang, 2016, A method for low-frequency noise suppression based on mathematical morphology in microseismic monitoring, Geophysics, 81 (3) (2016), pp. V159–V167. [23] Mousavi, S. M., and C. A. Langston (2016). Adaptive noise estimation and suppression for improving microseismic event detection, Journal of Applied Geophysics, 132, 116-124. [24] Mousavi, S. M., and C. A. Langston, (2016). Fast and novel microseismic detection using time-frequency analysis. SEG Technical Program Expanded Abstracts 2016: pp. 2632-2636. [25] Mousavi, S. M., and C. A. Langston (2017). Automatic Noise-Removal/Signal-Removal Based on the General-Cross-Validation Thresholding in Synchrosqueezed domains, and its application on earthquake data, Geophysics. 82(4), V211-V227. [26] Basu, Mitra. "Gaussian-based edge-detection methods-a survey." IEEE Transactions on Systems, Man, and Cybernetics, Part C 32, no. 3 (2002): 252-260. [27] Marr, David, and Ellen Hildreth. "Theory of edge detection." Proceedings of the Royal Society of London B: Biological Sciences 207, no. 1167 (1980): 187217. [28] Singh, Meghna, Mrinal K. Mandal, and Anup Basu. "Gaussian and Laplacian of Gaussian weighting functions for robust feature based tracking." Pattern Recognition Papers 26, no. 13 (2005): 1995-2005. [29] Shi, Yuhui, and Russell C. Eberhart. "Empirical study of particle swarm optimization." In Evolutionary Computation, 1999. CEC 99. Proceedings of the 1999 Congress on, vol. 3. IEEE, 1999.

14

[30] Kennedy, James. "Particle swarm optimization." In Encyclopedia of machine learning, pp. 760-766. Springer US, 2011. [31] Hassan, Rania, Babak Cohanim, Olivier De Weck, and Gerhard Venter. "A comparison of particle swarm optimization and the genetic algorithm." InProceedings of the 1st AIAA multidisciplinary design optimization specialist conference, pp. 18-21. 2005. [32] Selvi, V., and Dr R. Umarani. "Comparative analysis of ant colony and particle swarm optimization techniques." International Journal of Computer Applications (0975–8887) 5, no. 4 (2010). [33] Bai, Qinghai. "Analysis of particle swarm optimization algorithm." Computer and information science 3, no. 1 (2010): 180. [34] Peng, Chaoyong, Xiaoyi Zhu, Jiansi Yang, Bing Xue, and Yang Chen. "Development of an integrated onsite earthquake early warning system and test deployment in Zhaotong, China." Computers & Geosciences 56 (2013): 170-177. [35] Colombelli, S., A. Caruso, A. Zollo, G. Festa, and H. Kanamori. "AP wave‐based, on‐site method for earthquake early warning." Geophysical Research Letters 42, no. 5 (2015): 1390-1398. [36] Böse, M., R. Allen, H. Brown, G. Gua, M. Fischer, E. Hauksson, T. Heaten et al. "CISN ShakeAlert: An earthquake early warning demonstration system for California." In Early Warning for Geological Disasters, pp. 49-69. Springer Berlin Heidelberg, 2014. [37] Picozzi, Matteo, Aldo Zollo, Piero Brondi, Simona Colombelli, Luca Elia, and Claudio Martino. "Exploring the feasibility of a nationwide earthquake early warning system in Italy." Journal of Geophysical Research: Solid Earth 120, no. 4 (2015): 2446-2465. [38] Saramaki, Tapio. "Finite impulse response filter design." Handbook for digital signal processing (1993): 155-277. [39] Mirzaei, Shahnam, Anup Hosangadi, and Ryan Kastner. "FPGA implementation of high speed FIR filters using add and shift method." In2006 International Conference on Computer Design, pp. 308-313. IEEE, 2007. [40] http://www.xilinx.com/products/boards-and-kits/1-f3zdsm.html (Last Access 26-12-2017).

15