key: cord-0058496-ogi7izyk
authors: D’Angelo, Nicoletta; Adelfio, Giada; D’Alessandro, Antonino; Chiodi, Marcello
title: A Fast and Efficient Picking Algorithm for Earthquake Early Warning Application Based on the Variance Piecewise Constant Models
date: 2020-08-26
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58820-5_65
sha: 800a7ee593f937b370fc456f5420268445dc7a87
doc_id: 58496
cord_uid: ogi7izyk

An earthquake warning system, or earthquake early warning system, is a system of accelerometers, seismometers, communication, computers, and alarms that is devised for notifying adjoining regions of a substantial earthquake while it is in progress. This is not the same as earthquake prediction, which is currently incapable of producing decisive event warnings. The implementation of efficient and computationally simple picking algorithm is necessary for this purpose, as well as automatic picking of seismic phases for seismic surveillance and routine earthquake location for fast hypocenter determination. In this paper, a picking method, based on the detection of signals changes in variance, is proposed, taking advantage of a generalized linear model formulation of the investigated problem. An application to simulated data is provided.

The term "Earthquake Early Warning" (EEW) is used to describe real-time earthquake information systems that have the potential to provide a warning before significant ground shaking. Warning times range from a few seconds to a little more than a minute and are primarily a function of the distance of the user from the earthquake epicentre. In EEW system, generally, two approaches could be followed, called an onsite and regional warning. The principle of onsite or single-station warning [7] is to detect seismic energy at a single location and provide warning of coming ground shaking at the same location, i.e., detect the P-wave and predict the peak shaking. The regional warning makes use of a seismic network and typically combines information derived by several stations. Generally speaking, the EEW systems first detects earthquakes, in particular, Pwave first arrival and then transmits a useful warning. Given that the strongest ground shaking usually arrives at the time of, or after, the S-wave arrival, using the P-wave to provide warning has the potential to increase the warning time everywhere, to reduce the radius of the blind zone, and to potentially provide a warning at the epicentre. It is therefore important for a robust EEW system, the implementation of efficient and computationally simple picking algorithm. Automatic picking of seismic phases is also important in seismic surveillance and routine earthquake location for fast hypocenter determination. To be suitable for both application, a false alarm must be avoided and time picking must be as accurate as possible.

In this work, we propose a new automatic picking algorithm suitable for the implementation of EEW and in seismic surveillance. The algorithm, based on changes in variance, is tested on synthetic seismograms. This paper is structured as follows: first, in Sect. 2, an overview about the most widespread automatic picking algorithms is reported; secondly, the description of the proposed methodology is reported in Sect. 3; finally, in Sect. 4, the variance piecewise constant models are applied to seismic waveforms; the last section is devoted to conclusions and final remarks.

For correct early warning, it is fundamental to recognize early the beginning of a seismogram employing fast and robust algorithms. This is important because the algorithm must run in real-time and false alarm or lost event should be avoided. First arrival times on seismograms coincides with the arrival of the first P-wave. The time of the phase-detectionT i at a station i interpreted as the first P-phase arrival time, which is, of course, afflicted with an error i .T i may be written as:

where T 0 is the source time and t i is the travel time of a P-wave to station i. The coincidence trigger detects an event, if for any combination of a minimum number of stations (typically three or four) the condition

is met. is the maximum allowed difference between trigger times at neighbouring stations. This coincidence trigger works satisfying for local networks, where the number of stations and the aperture of the network is not large. For regional and global networks this simple event detection algorithm has to be modified. [12, Chapter 16] review the most widespread automatic picking algorithms and analyze their properties. Here we report a brief overview. It is also worth to notice that comparative works among different pickers have been carried out in literature [4, 11, 19] . [5, 6] introduce the concept of characteristic function (CF), where the 'character' of the seismic trace is specified. The CF is obtained by one or several non-linear transformations of the seismogram and should increase abruptly at the arrival time of a seismic wave. In addition to the calculation of the CF the next steps of a picking algorithm are the estimation of the arrival time from the CF and the quality estimation.

The Allen picker is a fast and robust algorithm, which also accounts for automatic quality assessment. However, as this algorithm is amplitude based only, it might miss emergent P-onsets. A comparative study by [11] shows that this algorithm tends to pick somewhat early compared to what an analyst would pick.

Another widely used picking algorithm is the one proposed by [8] . This algorithm is frequently applied, e.g. by 'Programmable Interactive Toolbox for Seismological Analysis' (PITSA, [17] ) and the picking system MannekenPix [4] .

In contrast to Allen's squared envelope function, this CF is sensitive to changes in amplitude, frequency as well as in phase.

The Baer and Kradolfer picker is also very fast and robust and quite userfriendly, as this algorithm only needs 4 input parameters. A shortcoming of this algorithm is the missing automated quality assessment. Several comparative studies [4, 11, 19] show how this picking algorithm tends to be somewhat late compared to manual P-picks.

When an earthquake signal occurs, the statistical properties of a seismogram change abruptly. Therefore, the measurement of statistical properties in a moving window are suitable for the determination of a CF and subsequent estimation of arrival times. The statistical properties of the seismogram might be characterized by its distribution density function and by parameters like variance, skewness and kurtosis. The latter two are parameters of higher order statistics (HOS) and are defined by [10] .

Though just amplitude-based, higher order statistics are quite sensitive even to emergent P-onsets. In combination with a sophisticated picking algorithm (e.g. [11] ), which exploits the entire information provided by the determined CF, it yields excellent results. If precisely tuned, the automated quality assessment proposed by [11] gives similar weights as the analysts. However, choosing the parameters for this sophisticated algorithm is quite difficult and needs a great experience.

The so called autoregressive-Akaike-Information-Criterion-piker (AR-AIC) proposed by [19] is based on the work by [2, 3, 13] and [21] .

A longer time series is divided into two locally stationary segments each modelled by an autoregressive (AR) process. The first segment represents noise, the second segment contains the signal. After estimating the two sets of AR parameters, two prediction errors are computed and then the minimum of the two-model Akaike-Information-Criterion (AIC) indicates the arrival time.

The AR-AIC picker is a highly more sophisticated algorithm based on information theory. The algorithm is computationally quite expensive and hence much slower than the other reviewed pickers. [1] considers the case of changepoint detection procedure for changes in variation, assuming that the variance function can be described by a piecewise constant function with segments delimited by unknown changepoints. Let y i be the outcome and x i be the observed sample, for i = 1, 2, . . . , n occasions. Let us assume that y i = μ i + i , where μ i is for instance a sinusoidal function representing the observed signal and i ∼ N (0, σ 2 i ) is an error temr. In this context, σ 2 i is a variance function approximated by a piecewise constant regression function with K 0 + 1 segments. An example is shown in Fig. 1 . For simplicity, the model for changes in variance after the k * th observation is

with λ,λ, and k * unknown and

Taking advantage of a generalized linear model formulation of the investigated problem, the test for stepwise changes in variance of a sequence of Gaussian random variables may be transformed equivalently to the case of testing for changes in mean of the squared residuals from an estimated linear model that accounts for the mean behaviour of the observed signal. The estimation of the mean signalμ can be carried out by using a common smoothing procedure, e.g., fitting a cubic smoothing spline to the data. Following a suggestion in [20] , a gamma generalized linear model (GLM) is fitted with a log-link function, with response given by the squared studentized residuals s i = (y i −ŷ i ) 2 /w i , witĥ y =μ and weights w i = 1 − h i , where h i is the ith diagonal element of the hat matrix H. According to this approach, testing H 0 against H 1 means that we are looking for a change in the mean of the residuals from a fitted linear model. The proposed approach can be considered as a wider version of the cumSeg models proposed in [16] for independent normally distributed observations with constant variance and piecewise constant means to detect multiple changepoints in the mean of the gene expression levels in genomic sequences by the least squares approach. The authors assume that the datum y i , ∀i is defined as the sum of the signal μ i and noise i ∼ N (0, σ 2 i ) and that μ i is approximated by a piecewise constant regression function with K 0 + 1 segments, that is:

Here, I(·) is the indicator function, such that I(x) = 1 is x is true, ψ represents the K 0 locations of the changes on the observed phenomenon, β 1 is the mean level for x i < ψ 1 , and δ is the vector of the differences in the mean levels at the change points. The authors proceed to take the cumulative sums of the jump-points model to get a convenient modelling expression that faces the discontinuities at the changepoints ψ k assuming a piecewise linear or segmented relationship. Therefore, looking for changes in variance, the model is specified as

where

and it has the advantage of an efficient estimating approach via the algorithm discussed in [14, 15] , fitting iteratively the generalized linear model:

. The parameters β 1 and δ are the same of Eq. (1), while the γ are the working coefficients useful for the estimation procedure [15] . At each step the working model in Eq. (2) is fitted and new estimates of the changepoints are obtained viâ ψ k =ψ k +γ k δ k iterating the process up to convergence. K * (< K) values are returned, producing the fitted model g(θ * i ) =β 1 +δV i1 + . . . +δ K * V iK * , where V ik = I(x i >ψ k ) for k = 1, 2, . . . , K * and the squared residuals are modelled as the response of a gamma GLM with logarithmic link function. Selecting the number of significant changepoints means selecting the significant variables among V 1 , . . . , V k , where K * is the number of estimated changepoints from model (1) . The author solves the model selection problem by using the lars algorithm by [9] . Thus, the fitted optimal model withK ≤ K * changepoints is selected by the generalized Bayesian Information Criterion (BIC Cn ), defined by:

where L is the likelihood function, edf is the actual model dimension quantified by the number of estimated parameters (including the intercept, the δ and ψ vectors), and C n is a known constant.

The choice of C n . The first issue concerns the value of C n to used in the BIC Cn criterion to select the changepoints. Therefore, by simulation, we assess the performance of different specifications of C n . Assuming for simplicity x i = i, i = 1, . . . , 700, the true variance used for the simulations is σ 2 i = 0.5+10I(i > 175), while the true mean is not specified, as the mean behaviour of the signal is assumed to be unknown. We consider BIC Cn with different values of C n and we report the results in Table 1 . Among the different examined specifications of C n , simulations reveal that C n = log log n has the best performance. Thus, we use this value for the provided analysis.

To test the performance of the algorithm, a dataset consisting of 100 waveforms over 60 s is simulated. The changepoints are three, equal to 10, 12, and 22 s. The first changepoint represents the arrival of the first P-wave, the second one represents the arrival time of the fist S-wave and the last one corresponds to the end of the seismic event [18] .

The simulated dataset is shown in Fig. 2 . Each signal has the same changepoints but different variances corresponding to each phase.

First, the algorithm is applied to the entire dataset. For each waveform, we obtain a different set of changepoints, estimated along with the signal. We report the output of the procedure applied to the first waveform of the dataset in Fig. 3 .

The dashed red lines are placed in correspondence of the changepoints estimated through the algorithm. As we may notice, the proposed algorithm succeeds in the identification of the arrival of both the P-wave and S-wave. The straight lines are three changepoints closer to the true values (10, 12 and 22 s), and therefore they represent the changepoints, among all the estimated ones, that a further algorithm should pick as the three relevant changepoints.

In Table 2 the summary statistics of the three relevant changepoints estimated along the whole dataset are reported. We may see that the mean values are close to the true values and that the smallest mean squared error is the one of the first changepoint, that is, the arrival of the first P-wave. Indeed, the estimates of the other two changepoints are less precise. This effect can be related to the occurrence of the first changepoint in correspondence of the most abrupt change in variation of the signal, making it easier to detect. Simulated Trend. Then, to test the performance of the algorithm in different settings, further analysis is proposed.

The analyses are carried out assuming underlying trends equal to μ 1 = sin πx and μ 2 = 0.4 sin 100πx, where x in equally spaced in the range [1 : 6000]. These are superimposed and added to the simulated signals. In Fig. 4 (left panel) the first waveform of the simulated dataset is shown, together with the simulated trends, plotted in red. Therefore, the algorithm can account for a mean behaviour of the observed signal, without any loss of precision. In particular, two different models are considered to take into account the mean behaviour of the signal. The first is linear and the second is obtained fitting a smooth spline with 100 knots. These are represented in Fig. 4 (right panel) in red and green, respectively. In Table 3 the results of the application along with the whole dataset, to which the mean behaviour is added, are shown.

We may notice that the estimates of the relevant changepoints are more precise when taking into account the mean behaviour of the signal through a smooth spline, rather than ignoring it. Indeed, as we may notice from Table 3 , the mean values are closer to the true changepoints and the Mean Square Error Values are smaller. Post Selection. Finally, we propose a further algorithm to detect, among the estimated changepoints, the three ones corresponding to the arrival of the first P-wave, the arrival of the first S-wave and the end of the seismic event. In particular, we compare the ratio between the variances of subsequent phases of the signal and select only the biggest three. As far as the first algorithm, we fit a smoothing spline to account for the mean behaviour of the signal (see Fig. 5 and Table 4 ). Of course, the post-selection adds uncertainty to the estimates obtained through the application of the algorithm, so the Mean Squared Error values obtained are the sum of the Mean Squared Errors from the true values plus the uncertainty due to the post-selection. Therefore, as the Mean Squared Error value referred to the first changepoint is equal to the one of the proposed procedure, this means that the first changepoint is always correctly selected by the post-selection algorithm. Overall, we may conclude that the algorithm succeeds most of the times in picking the first relevant changepoints, that is the one corresponding to the arrival of the first P-wave. Concerning the second changepoint, that is the one corresponding to the arrival of the first S-wave, this is more easily detected if the mean behaviour of the underlying trend is taken into account through a smoothing spline with the right choice of the number of knots. Indeed, some preliminary analyses, not reported here for brevity, show that a smoothing spline with too many knots tends to overfit the signal, hiding that variation of the signal that could make clearer the detection of the changepoints. Finally, we see that in all the considered setting, the last relevant changepoint, that corresponds to the end of the seismic event, is the most difficult to detect. This may related to the less abrupt change in variation at the end of the signal and, therefore, both the proposed algorithms and the post-selection one can not be accurate in this case.

The proposed approach is a highly sophisticated algorithm based on changes in variation.

All the analysis carried out in the present work account for simulated data, showing that the proposed algorithm accounts for unknown trends of the signals. Future work will deal with real data, depending on the different characteristics of the analysed region that may affect the observed waveforms. Moreover, taking into account the trend of the signal, by fitting splines, is important to remove all those changes in the waveform that are not attributable to the seismic event.

Moreover, in future work, we will perform comparison between the proposed method and the other already existing algorithms, as the Allen picker, Baer and Kradolfer picker, higher order statistics and AR-AIC picker, to asses their general performance, also in terms of computational efficiency and robustness.

Change-point detection for variance piecewise constant models. Commun

Markovian representation of stochastic processes by canonical variables

Autoregressive model fitting for control

Toward Three-dimensional Crustal Structure of the Dead Sea Region From Local Earthquake Tomography

Automatic phase pickers: their present use and future prospects

Automatic earthquake recognition and timing from single traces

The status of earthquake early warning around the world: an introductory overview

An automatic phase picker for local and teleseismic events

Least angle regression

Statistik: Lehr-und Handbuch der angewandten Statistik

Automated determination of p-phase arrival times at regional and local distances using higher order statistics

Automated event and phase identification

Automatic detection of onset time of seismic waves and its confidence interval using the autoregressive model fitting

Segmented: an R package to fit regression models with broken-line relationships

Estimating regression models with unknown break-points

Efficient change point detection for genomic sequences of continuous measurements

Programmable Interactive Toolbox for Seismological Analysis

Efficient global matrix approach to the computation of synthetic seismograms

Robust automatic p-phase picking: an on-line implementation in the analysis of broadband seismogram recordings

Exact and approximate REML for heteroscedastic regression

A new efficient procedure for the estimation of onset times of seismic waves