Automatic detection of wave arrival

Yannick Apfel, summer semester 2019

Introduction

The automatic detection of wave arrival or automatic wave picking is an important issue in many fields of science and involves the algorithmic determination of the onset of a transient signal. As this topic is most commonly addressed in many fields of seismology (Seismic methods), e. g. earthquake detection, the following article is based mainly on findings in that area of science. [1]

Concept

The first step in a seismic analysis is usually the picking of the P-wave arrival. The term picking describes the reliable and accurate estimation of the onset time of a definite seismic arrival [2], whereas detectors are processes for detecting the presence of a specific seismic phase type (e. g. P- or S-wave) [3]. Before the 70s, this task was done completely manually, yielding results based on the pattern recognition capabilities of experienced seismologists. Considering ever-increasing amounts of data in seismology and the tedious nature of the task, many different approaches to automatically identify incoming wave fronts have been done. [4] [1]

Once a wave arrival picking algorithm detects a wave, the detection time is assigned to a specific geophone. Once a certain number of geophones in a seismic monitoring system are triggered, an event is considered as detected and the system would, in the case of an earthquake, proceed to, for instance, compute the location of a hypocenter. [5]

The automatic pickers have drastically improved data acquisition, long-term continuous recording and archiving of seismic signals for seismic stations because they have reduced the data that needs to be stored. Although seismic stations or networks still monitor and process all incoming data in real-time, they only permanently store data acquired during a triggered state, meaning that a trigger algorithm has detected a typical seismic signal against a constant seismic noise backdrop by its first arrival and started the recording process which is then stopped when the algorithm declares the end of the seismic event. In addition to the signal being recorded during ‘trigger active’ times, seismic networks usually add a certain amount of data to the event file before and after the ‘trigger active’ time-window, called the pre- and post-event-times. [6]

The onset time of a signal is generally defined as the point where the signal starts to differ from the background noise of a seismologic recording [1]. A very crucial measure in this regard is the signal-to-noise ratio. It compares the power level of a desired signal to that of undesired background noise. The most common definition, using decibels to express the power, is as follows whereby A signifies the root-mean-square of the signal amplitude. [7] [8]

$\begin{array}{l}\displaystyle SNR=\frac{P_{Signal}}{P_{noise}}=10\log_{10} \left( \frac{A_{S,RMS}}{A_{N,RMS}} \right)\end{array}$

The most prominent difficulties for designing algorithms for automatic wave detection are, in general, the fact that the shapes of arriving waves usually differ strongly depending on the seismic event that caused them and that the reconstruction of their travel paths is difficult. Another problem is the fluctuation of the signal-to-noise ratio due to different seismic noise levels. This strongly depends on the location at which seismic recordings are done. [4]

Figure 1 shows an example of a recording of a seismic event with the P- and S-wave arrival picks.

Figure 1: Example of seismic P- and S- waves and the time picks for their first arrival

Source: extracted from [9] under CC BY-SA 4.0, see DOI: 10.1190/GEO2014-0500.1

Classification

There is a wide range of algorithms that have been proposed and implemented as wave arrival pickers. A common way to classify the methods is to distinguish between:

Single-level and multi-level (=array processing) algorithms and
hybrid and non-hybrid algorithms.

Single-level algorithms are based on single-component or multi-component recordings from an individual receiver level. E. g. the recording of the black line in Figure 1 represents the sum of three single-component recordings of a geophone, thus making it a multi-component recording. As a result, single-level algorithms don’t make use of data from other locations within a receiver array. As a side note, e. g. to localize the hypocentre of an earthquake, it is necessary to receive seismic data from more than one location around the earthquake. Single-level algorithms can be further divided into window-based and non-window-based. While the window-based algorithms derive parameters for arrival-time calculation based on data from a time window of pre-defined starting point and duration within the recordings, non-window-based algorithms don’t require a preselected time window but instead scan the whole signal. [9]

In contrast, multi-level algorithms operate on information on multiple receiver levels within the array. That means that they exploit the similarity of waveforms from nearby events.

Hybrid algorithms or strategies combine two or more algorithms to improve on accuracy and strive to combine their benefits as no single algorithm is optimal in all conditions. [9]

Another way to distinguish between the picking algorithms is whether they work in the time or the frequency domain of a signal.

The following is a list of common single-level algorithms: (apart from the amplitude-threshold picker and the AIC algorithm, all are window-based) [9]

Amplitude threshold
STA (=short time average) to LTA (=long time average) ratio
MER (=modified energy ratio)
MCM (=modified Coppens’ method)
AIC (=Akaike information criterion)
PAI (=phase arrival identification) -S (=skewness) and -K (=kurtosis)
SK (=short-term kurtosis) to LK (=long-term kurtosis)

One common multi-level approach are cross-correlation techniques.

Newer or less common approaches include methods based on fractals, neural networks, digital image segmentation and adaptive multiband processing.

Algorithms

Among the myriad of available algorithms, some are looked at more closely in the following section.

Amplitude threshold trigger

The simplest trigger algorithm is the amplitude threshold trigger. It detects any amplitude of a seismic signal exceeding a pre-set threshold amplitude [6]. However, the method proves to be less effective for small amplitude signals and/or signals with high noise level [1]. Thus, this very simple approach is rarely used in weak-motion seismology but is seen as a standard in strong-motion seismic instruments. A common variation of the amplitude threshold trigger (it makes use of the instant signal amplitude) is an RMS-threshold trigger which uses the RMS-values of the amplitude within a short time window, rendering it less sensitive to spike-like seismic noise. [6]

STA/LTA

General description

Despite their age, the STA/LTA pickers remain to be among the most widely used algorithms, e. g. in weak-motion seismology [10]. The following workflow is an example of how this kind of algorithm can be implemented:


Figure 2: STA/LTA ratio algorithm [11]

The algorithm continuously processes filtered seismic signals in two moving time windows. A short one, for the STA (=short time average) and a significantly longer one, for the LTA (=long time average). For both time windows the average absolute amplitude is calculated. Given the different lengths of the time windows this results in the effect that the STA measures the ‘instant’ amplitude of the seismic signal and the LTA measures the current average seismic noise amplitude. Next, the ratio of both parameters is calculated and continuously compared to a STA/LTA trigger threshold level (e. g. 10) pre-defined by the user. If the ratio exceeds the predefined threshold, a channel trigger is declared. Channel de-triggering happens in a similar manner, when the current STA/LTA ratio falls below a predefined de-trigger threshold (e. g. 2). In reality, the STA/LTA triggers are usually slightly more complicated, however, they essentially work like described above [6].

The general idea behind this strategy is that during times of only noise the ratio will usually stay relatively constant. As a signal arrives, e. g. the wave front of a P-wave, the STA captures the change much more quickly than the LTA, resulting in a rising STA/LTA-value. In a way, the algorithm represents a kind of signal-to-noise-ratio. [10]

Generalized expressions for the STA and LTA at the i-th time sample can be formulated as follows:

$\begin{array}{l}STA_{i}=\frac{1}{ns} \sum _{j=i-ns}^{i}CF_{j} \\ LTA_{i}=\frac{1}{nl} \sum _{j=i-nl}^{i}CF_{j}\end{array}$

Where ‘ns’ and ‘nl’ represent the number of samples in the corresponding time windows. ‘CF’, as the characteristic function, can represent any kind of seismic wave information and has traditionally been subjected to experimentation by researchers. E. g. formulations of CF in the past have included the amplitude or wave energy and their derivatives or variance combined in weighted sums. Others have used envelope functions of the signal instead. [9] [10]

Using the algorithm, with amplitudes as the characteristic function, on the signal depicted in Figure 1, one yields the following:

Figure 3: STA/LTA ratio for signal in Figure 1 arrival

Source: extracted from [9] under CC BY-SA 4.0, see DOI: 10.1190/GEO2014-0500.1

Parametrization

The effectiveness of the algorithm strongly depends on the meaningful setting of certain parameters which always represents a tradeoff between sensitivity to seismic events and the minimization of false triggers. Thusly, the goal of searching for optimal parameter settings is to achieve the highest possible seismic station sensitivity for a given type of seismic signal at a tolerable number of false triggers. [6]

The setting of parameters depends strongly on the goal of the application, seismic noise levels present and sensor hardware used [6].

Modern versions of the algorithm continuously adjust the seismic station sensitivity to changes in background noise. [6]

Essentially, the following four parameters must be set: [6]

STA window duration
LTA window duration
STA/LTA trigger threshold level
STA/LTA de-trigger threshold level

As a rule of thumb, the STA window is usually longer than a few periods of the typically expected seismic signal while the LTA should be longer than a few periods of seismic noise fluctuations. As an example, one could set the STA time window to 2-3 times the dominant period of the signal and the LTA time window to 5-10 times the LTA time-window. [5] [9]

For reasons of causality the STA time window should always lead the LTA time window. To avoid statistical independence the time windows shouldn’t overlap. [9]

Other parameters include: trigger filters, PEM (=pre-event time), PET (=post-event time), trigger voting scheme, etc. [6]

Fields of application

The STA/LTA trigger is most beneficial at seismically quiet sites (e. g. natural seismic noise like marine noise) or sites with continuous but changing man-made seismic backdrops. The algorithm has proven less effective in the presence of irregular, high amplitude man-made seismic noise. [6]

Similar algorithms

The MER and MCM algorithms are two examples of algorithms similar to the STA/LTA ratio. The MER, for instance, is an extension of the STA/LTA with equally sized pre- and post-sample windows. With ‘w’ as the window length and ‘x_i’ as the input time-series, the MER is defined as follows:

$\begin{array}{l}\displaystyle MER_{i}= \left( \frac{ \sum _{j=i}^{i+w}x_{j}^{2}}{ \sum _{j=i-w}^{i}x_{j}^{2}}\ast \vert x_{i} \vert \right) ^{3}\end{array}$

Because the ratio of energies is computed using post- and pre-sample windows, the time index with the maximum MER represents the arrival time pick.

The MCM algorithm is very similar to the MER. The latter has a slightly different scheme for the window size selection. Additionally, to make the algorithm more robust, a stabilization constant ‘β’ is introduced which reduces rapid fluctuations of the MCM-curve, thus reducing the number of false picks. [9]

PAI-K

The PAI-K and PAI-S as well as SK/LK represent algorithms using higher order statistics on sliding time windows. The PAI-K algorithm uses kurtosis values to form the characteristic function. The kurtosis can be defined as follows:

$\begin{array}{l}\displaystyle K=\frac{E \left( \left( X-E \left( X \right) \right) ^{4} \right) }{E \left( \left( X-E \left( X \right) \right) ^{2} \right) ^{2}}=\frac{m_{4}}{m_{2}^{2}}\end{array}$

‘E(X)’ denotes the expected value of ‘X’ (the expectation of a continuous distribution) while ‘m_k’ represents the central statistical moment of order ‘k’. k takes different values depending on the type of statistical distribution. [9]

The PAI-S, which takes advantage of the skewness of a distribution, can be calculated in a similar manner: [12]

$\begin{array}{l}\displaystyle S=\frac{E \left( \left( X-E \left( X \right) \right) ^{3} \right) }{E \left( \left( X-E \left( X \right) \right) ^{2} \right) ^{\frac{3}{2}}}=\frac{m_{3}}{m_{2}^{\frac{3}{2}}}\end{array}$

For a symmetrical distribution the value of ‘S’ becomes 0. [12]

AIC

The AIC algorithm is based on the concept that microseismic signals are non-stationary and can be approximated by dividing a signal into locally stationary segments where each segment is treated as an autoregressive process (noise and signal). The method assumes that variance changes occur only at the arrival time (between the two stationary windows): [9] [13]

$\begin{array}{l}\displaystyle AIC \left( k \right) = \left( k-M \right) log \left( \sigma _{1,max}^{2} \right) + \left( N-M-k \right) log \left( \sigma _{2,max}^{2} \right) +C\end{array}$

‘AIC(k)’ is the AIC at the k-th data sample of a microseismic waveform of length ‘N’. ‘M’ denotes the order of the autoregressive model, while the ‘σ’ coefficients are the variances in the two windows not covered by the autoregressive process. The order of the autoregressive model is estimated using the window containing the noise. The AIC function computed provides a measure of the model fit. Optimal separation of the two stationary time series is indicated by the time index associated with the minimum value of AIC (AIC defines the onset point as a global minimum) [9]. In essence, the point where the AIC function is minimized is regarded as the wave arrival time. There exists a variation of the algorithm where the AIC is directly obtained from the time series without using the autoregressive model coefficients. In that case the AIC is calculated as follows:

$\begin{array}{l}\displaystyle AIC \left( k \right) =k\log \left( var \left( x \left( 1,k \right) \right) \right) + \left( N-k-1 \right) \log \left( var \left( x \left( k+1,N \right) \right) \right)\end{array}$

Analogous to the first expression, ‘k’ is the sliding sample that divides the two different variances of the time series.

Waveform cross-correlation

Waveform cross-correlation is the dominant example for multi-level array-based algorithms. These take advantage of similar characteristics of waveforms across a receiver-array or multiple events on a single receiver. Other examples of these include image-processing techniques, global-optimization-based techniques and beamforming (delay and stack) of waveforms. [9]

The normalized cross-correlation of two digital waveforms ‘x_i’ and ‘y_i’ can be expressed as follows:

$\begin{array}{l}\displaystyle \varphi _{xy} \left( \tau \right) _{norm}=\frac{ \varphi _{xy} \left( \tau \right) }{\sqrt[]{ \varphi _{xx} \left( 0 \right) \varphi _{yy} \left( 0 \right) }}\end{array}$

The expression represents the cross-correlation of two digital waveforms at a time lag of ‘τ’. The denominator contains their zero-lag autocorrelation values. A correlation value of 1 indicates a perfect match whereas -1 indicates that the waveforms have opposite polarity. [9]

The following form can be assumed for microseismic data recordings on two different receivers:

$\begin{array}{l}x_{1} \left( t \right) =s \left( t \right) +n_{1}\left( t \right) \\ x_{2} \left( t \right) =a \left( t- \tau \right) +n_{2} \left( t \right)\end{array}$

The two signals have a delay of ‘τ’, with two different levels of background noise ‘n₁(t)’ and ‘n₂(t)’. ‘s(t)’ denotes the signal and ‘a’ denotes the amplitude ratio between them. [9]

In schemes like this a reference (pilot) waveform is usually chosen from one of the receivers to be cross-correlated with the waveforms of all the other receiver levels. This yields the time delay. The workflow can look like this: [9]

Choose pilot waveform based on high S/N-ratio
Guess an initial arrival time (e. g. manually) and stack waveforms from all receiver levels accordingly to yield another (the new) pilot waveform
Iteratively update pilot waveform obtained through stacking after time delay estimation and lag adjustments until time delay converges to value below pre-defined threshold value

An important parameter to be set in this algorithm is the size of the correlation window. [9]

Literature

J. H. Kurz, C. U. Grosse and H.-W. Reinhardt, "Strategies for reliable automatic onset time picking of acoustic emissions and of ultrasound signals in concrete," Ultrasonics 43, pp. 538-546, 2005.
H. Dai and C. MacBeth, "Automatic picking of seismic arrivals in local earthquake data using an artificial neural network," Geophysical Journal International 120, pp. 758-774, 1995.
S. E. J. Nippress, A. Rietbrock and A. E. Heath, "Optimized automatic pickers: application to the ANCORP data set," Geophysical Journal International 181, pp. 911-925, 2010.
A. G. Hafez, M. Rabie and T. Kohda, "Seismic noise study for accurate P-wave arrival detection via MODWT," Computers & Geosciences 54, pp. 148-159, 2013.
K. Munro, "Automatic event detection and picking of P-wave arrivals," CREWES Research Report 16, pp. 1-10, 2004.
A. Trnkoczy. [Online]. Available: http://gfzpublic.gfz-potsdam.de/pubman/item/escidoc:4097:3/component/escidoc:4098/IS_8.1_rev1.pdf. [Accessed 17 July 2019].
N. Karamzadeh, G. J. Doloei and A. M. Reza, "Automatic Earthquake Signal Onset Picking Based on the Continuous Wavelet Transform," IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING 51, No. 5, pp. 2666-2674, 2013.
"Signal-to-noise_ratio," [Online]. Available: https://en.wikipedia.org/wiki/Signal-to-noise_ratio. [Accessed 15 July 2019].
J. Akram and D. W. Eaton, "A review and appraisal of arrival-time picking methods for downhole microseismic data," GEOPHYSICS 81, No.2, pp. KS 71 - KS 91, 2016.
V. M. Pikoulis and E. Z. Psarakis, "A ratio test for the accurate automatic P-wave onset detection," ICSP2010 Proceedings, pp. 2621-2624, 2010.
A. E. Ruano, G. Madureira, O. Barros, H. R. Khosravani, M. G. Ruano and P. M. Ferreira, "Seismic detection using support vector machines," Neurocomputing 135, pp. 273-283, 2014.
L. Küperkoch, T. Meier, J. Lee and W. Friederich, "Automated determination of P-phase arrival times at regional and local distances using higher order statistics," Geophysical Journal International 181, pp. 1159-1170, 2010.
Z. Zhou, R. Cheng, Y. Rui, J. Zhou and H. Wang, "An improved automatic picking method for arrival time of acoustic emission signals," IEEE Access, pp. 75568-75576, 7 June 2019.

Seitenhierarchie