Bayesian inference in epidemics: linear noise analysis

1.   Introduction

Earthquakes are monitored by using seismic stations that record digitally the motion of the ground and represent it as a function of time in seismograms. However, the seismic signals generated by some anthropogenic activities, including explosions, are frequently misinterpreted as natural earthquakes. Several methods to discriminate between earthquakes and explosions have been proposed. Commonly used methods compare in the time-domain the seismic phases recorded in seismograms. For example, in ^[1] the authors compared the phases $Pg$ and $Lg$ and found that in the tectonically stable eastern USA, the mean $Pg/Lg$ ratios are $0.5$ for earthquakes and $1.25$ for explosions. In southern Russia, the $Pg/Lg$ ratio for earthquakes was found to be about $1.3$, and for explosions $3.2$. However, current discriminating methods are not always successful. A data from the Powder River Basin in the western USA and the Altai-Sayan region in Russia suggest that a simple time-domain $Lg/P$ amplitude ratio analysis does not separate signals originating from earthquakes and explosions consistently enough to provide a reliable discrimination ^[2]. Also, the $S/P$ amplitude ratio does not reliably separate explosions from locally recorded natural seismicity in Southern California ^[3].

In 2014, a $5.2$ magnitude earthquake struck in Arizona USA, triggering a large number of aftershocks. A few years before, a sequence of mining explosions were carried out in the same region, where subsequently the $5.2$ magnitude earthquake and its aftershocks occurred. A nearby seismic station located about $160$ km away, recorded all these events, thus providing an important data set of earthquake and explosion seismograms that makes it possible to apply and investigate different seismic discriminating techniques (see $e.g.$ ^[7]). In this study, we use a data set originating from these explosions and the aftershocks of the $5.2$ magnitude earthquake, in order to investigate the Dynamic Fourier Analysis (DFA) as a discriminating method.

The paper is organized as follows. In section 2, we present the background of the analyzed data. Section 3 is devoted to an introduction of the Dynamic Fourier Analysis and its application for discriminating signals originating from natural earthquakes and mining explosions. Results and a discussion on the suitability of the Dynamic Fourier Analysis technique are presented in section 4. Conclusions are given in section 5.

2.   Data Background

The earthquakes in this study represent a set of aftershocks of the June 26,2014, $5.2$ magnitude intraplate earthquake. That is, this earthquake located between the sates of Arizona and New Mexico in the USA (Figure 1) occurred far from active tectonic boundaries. The explosions analyzed were carried out in a large surface copper mine, as part of quarry blasts activities. We selected earthquakes and explosions with similar local magnitude, in the range $3.0-3.3$, and occurring close to each other, within a $10$ km radius. The magnitude used in this study is measured by the Richter scale (M$_L$) that assigns a number to quantify the size of an earthquake or an explosion. We collected seismograms of earthquakes and explosions recorded at seismic station TUC of the IU network. We also downloaded broadband seismograms from the Incorporated Research Institutions for Seismology Data Management Centers (IRIS DMC). The instrument response was removed and windows of data were selected using the Seismic Analysis Code (SAC) software. We analyzed the seismograms representing the displacement of the ground (in nm) in the vertical direction (Z-component). Information about the seismic station and the selected events is presented in Tables 1 and 2, respectively. Figure 2 shows examples of seismograms recorded at the TUC stations from one earthquake and one explosion.

3.   Dynamic Fourier Analysis (DFA)

A physical process can be described either in the time domain, by the values of some quantity $h$ as a function of time $t$, $e.g.$, $h(t)$, or in the frequency domain, where the process is specified by giving its amplitude $H$ (generally a complex number) as a function of frequency $f$, that is $H(f)$, with $-\infty < f < \infty$. For many purposes it is useful to consider $h(t)$ and $H(f)$ as two different representations of the same function. It is possible to go back and forth between these two representations by means of the Fourier transform formulas:

The Fourier transform $H(f)$ converts waveform data in the time domain into the frequency domain. The Inverse Fourier transform $h(t)$ converts the frequency domain components back into the original time-domain signal. A frequency-domain plot shows how much of the signal lies within each given frequency band over a range of frequencies.

In order to analyze a statistical time series, it must be assumed that the structure of the statistical or stochastic process generating the observations is essentially invariant through time. The conventional assumptions are summarized in the condition of stationarity as follows, if a time series $x_t$ is stationary, its second-order behavior remains the same, regardless of the time $t$. Based on the Spectral Representation Theorem ^[8,9], a stationary series is represented through the superposition of sines and cosines that oscillate at various frequencies. Therefore, a stationary time series can be represented with combination of sine and cosine series.

Next, we introduce definitions that will be used when applying the DFA.

3.1.   Tapering

A Discrete Fourier Transform (DFT) is used for discrete time series with a finite number of samples from a process that is continuous in time. When discontinuity occurs, the signal value abruptly jumps, yielding spectral leakage $i.e.$ the input signal does not complete a whole number of cycles within the DFT time window. Consequently, it is required to multiply finite sampled time series by a windowing function, or a taper that decreases and approaches zero. Tapering is a process that reduces spectral leakage by multiplying the original time sequence by a function of time, such that the time series decreases and approaches zero, thus avoiding sudden discontinuities at the first and last point of the time series.

3.2.   Estimation of spectral density with Daniell kernel

A stationary process ${X_t}$ can be defined by using linear combinations of the form:

where $0 \leq \lambda \leq \frac{1}{2}$ is a fixed constant and $A_j$, $B_j$ for $(j=1, m)$ are all uncorrelated random variables with a mean $=0$ and

We assume $\sum_{j=1}^{m} \sigma_j^2=\sigma^2$, so that the variance of the process ${X_t}$ is $\sigma^2$, where the spectral density $f(\lambda)$ satisfies:

Then for large values of $m$ the process given in (3.1) approximates the stationary process with spectral density $f$.

To estimate the spectral density, we define the estimators as weighted averages of periodogram values ($I$) for frequencies in the range $(j-m)/n$ to $(j+m)/n$. In particular:

where the Daniell kernel $W_m(k)$ with parameter m is a centered moving average which creates a smoothed value at time t by averaging all values between times $t-m$ and $t+m$ (inclusive). We define:

The smoothing formula $\{u_t\}$ for a Daniell kernel with $m = 1$ corresponds to the three weights $(\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$ and is given by:

Applying the Daniell kernel again the smoothed values $\{\hat u_t\}$ results in more extensive smoothing, averaging over a wider time interval

Consequently, applying the Daniell kernel transforms the spectral windows into a form of Gaussian probability density function; for more details see ^[10].

3.3.   Detection of Long memory

In this subsection we discuss a unit root methodology that is often used to test long memory behavior, the so-called Augmented Dickey Fuller Test (ADF).

The geophysical time series discussed in this study has a complicated dynamic structure that is not easily captured by traditional methods, such as a simple first-order linear autoregressive (AR(1)) models (see ^[4]). Therefore, there is a need for models that can capture such a unique dynamic structure. It is possible to augment the basic autoregressive unit root test to accommodate general autoregressive moving-average (ARMA($p, q$)) models with unknown orders ^[5] and their test is referred to as the augmented Dickey Fuller (ADF) test (see ^[6]). The ADF is a very powerful test that can handle more complex models. It tests the null hypothesis that a time series $y_t$ is a unit root against the alternative that it is stationary, assuming that the data dynamics has an ARMA structure. The ADF test is based on estimating the test regression:

where $D_t$ is a vector of deterministic term. The $p$ lagged difference terms, $\Delta y_{t-j}$, are used to approximate the ARMA structure of the errors, and the value of $p$ is set so that the error $\epsilon_t$ is serially uncorrelated. An alternative formulation of (3.3) is

where $\pi = \phi-1$. In practice, the test regression (3.4) is used since the ADF t-statistic is the usual t-statistic reported for testing the significance of the coefficient $y_{t-1}$. An important practical issue for the implementation of the ADF test is the specification of the lag length $p$. If $p$ is too small then the remaining serial correlation in the errors will bias the test. If $p$ is too large, the power of the test will suffer. For the Augmented Dickey Fuller $t$-statistic test, small $p$-values suggest the data is stationary. Usually, small $p$-values are those that are smaller than a $5 \%$ level of significance.

3.4.   Application of Dynamic Fourier Analysis

We begin this subsection by using the ADF tests to test for a unit root in the seismic waves generated by the earthquake and explosion time series. To illustrate the ADF test procedure, we consider the seismograms of earthquake and the explosion time series as recorded by the seismic station TUC of the IU network. Figure 2 shows the seismograms recorded for earthquake 1 and explosion 1 from Table 2.

The ADF test was performed for all events (explosions and earthquakes) listed in Table 2. The summary statistics of the ADF results are shown in Table 3.

The ADF tests the following:

$H_0:$ There is a unit root for the time series, or

$H_a:$ There is no unit root for the time series and the series is stationary.

Because the computed $p$-values are lower than the significance level $\alpha = 0.05$, we reject the null hypothesis $H_0$ for all $8$ events, and accept the alternative hypothesis $H_a$. Thus all studied events represent stationary time series.

As observed from the ADF test (Table 3) and the analyzed seismograms, the seismic signals are stationary. Thus the recorded seismograms represent a stationary process and therefore a dynamic stationary Fourier analysis is needed.

The purpose of the dynamic series analysis is to depict the spectral behavior of the seismic signals as they evolve over time. The Fourier Transform is used to extract this information by calculating a base 2 logarithm ($\log_2$); in that case it is required that the length or range of the time series to be evaluated contains a total number of data points equal to $2^m$, where $m$ is a positive integer. If the total number of data points is not a power of two, it is also possible to express the Fourier transform of length $n$ in terms of several shorter Discrete Fourier Transforms (DFT). The Fast Fourier Transforms (FFT), an algorithm for calculation of the DFT, is applied in the present analysis. We obtain the spectral series from a short segment of the data using the DFT. After characterization of the first segment, the spectral analysis is performed on a new segment that is shifted from the first segment. This process is repeated until the end of the data. We assume that $x_t$ represents the series of interest, where $t=1, ......, 6016$. The analyzed data segments are $\{x_{t_k+1}, ......, x_{t_k+256}\}$, for ${t_k = 128k}$, where $k = 0, 1, ......, 45$. Therefore, the first segment analyzed is $\{x_1, ....., x_{256}\}$, the second segment analyzed is $\{x_{129}, ......., x_{384}\}$, and so on. Each segment of 256 observations was tapered using a cosine bell, and spectral estimation is performed using a repeated Daniell kernel with weights $\frac{1}{9}\{1, 2, 3, 2, 1\}$. These coefficients are obtained with the $R$ software by applying the kernel option.

One interpretation of the results of the time-frequency analysis presented above is to consider it based on local transforms of the data $x_t$, given by:

where

In the above formula $h_t$ is a taper and $m$ is a fraction of $n$. In this study, the number of observations is $n = 6016$, window size $m = 256$, overlap $t_k = 128k$ for $k = 0, 1, ...., 45,$ and $h_t$ is a cosine bell taper over 256 points. In equations (1) and (2) above, $j$ indexes the frequency $\omega_j = j/m$, for $j = 1, 2....... [m/2]$ and $k$ indexes the location, or time shift of the transform. We recall that the transforms are based on tapered cosines and sines that are zeroed out over various time intervals.

4.   Results and Discussion

The DFA technique was applied to the recorded seismograms of the events listed in Table 2, and the results are presented in Figures 3 and 4 for the natural earthquakes and the explosions, respectively. The results are also presented numerically as percentages of the total power spectra for $50-$ seconds time windows, and $2-$ Hz frequency windows in Table 4.

The power spectra were evaluated by first estimating the periodograms ^[10], which indicate large and small periodic components of seismic signals in the respective frequencies. This estimation was performed using DFT with a $256$ window size. In order to obtain a stable estimation, we computed the smoothed periodograms by using the Daniell kernel window (see subsection 3.2). Furthermore, we used the tapering process to avoid the spectral leakage of the power spectrum for earthquakes and explosions used in this study ^[11].

There is a considerable difference in the percentage of total power for the time series from earthquakes and explosions; see Figures 3, 4, 5 and Table 4. The results are shown in terms of estimated spectra for frequencies up to $10$ Hz (the folding frequency is $20$ Hz) for starting (time) $t_k = 128k$ with $k = 0, 1........... 45$. The spectral analyses were performed using $R$ and Python.

The power spectrum of the earthquake $S$-waves starts from low frequencies and may remain strong for a long time window. On the other hand, the explosion power spectrum is strong at lower frequencies and gradually tapers off at high frequencies, especially at those greater than about 2 to 4 Hz; see Figures 3, 4, 5 and Table 4.

Specifically, for both earthquakes and explosions, $P$-waves and $S$-waves are identifiable from the power spectra (see Figures 3 and 4) as sequences of picks occurring at around $25$ seconds for the $P$-waves and $50$ seconds for the $S$-waves. For both type of events, the $P$-waves are with very-low amplitudes at all frequencies. There is no evidence that is possible to discriminate between earthquakes and explosions based on $P$-wave spectra. On the other hand, the $S$-waves exhibit a high-amplitude pick between $2$ and $4$ Hz for both type of events. More interestingly, the $S$-waves show characteristic peak frequencies ($\sim$ 6-8 Hz) that is present for earthquakes but not for explosions. In order to further explore this characteristic peak in the $S$-wave frequency content, the values corresponding to the percentages of power in $2-$Hz frequency windows were plotted in Figure 5. It can be noticed that for frequencies between $2$ and $4$ Hz explosions present higher values than earthquakes (Figure 5), and that tendency is reversed for frequencies between $6$ and $8$ Hz, where earthquakes present higher values than explosions.

Identifying the physical mechanisms that cause this notable difference in the frequency content of the $S$-waves of earthquakes and explosions is beyond the scope of this paper. However, we observe that the $6-8$ Hz $S$-wave's peak detected for earthquakes but not for explosions, could be the a direct result from different source mechanisms. In general, it is still not well understood how explosions generate $S$-waves; while natural earthquakes simply shear slip on a fault plane. Explosions are modeled as pressure pulses on a sphere, which in theory should generate only $P$-waves ^[12]. The source depth could also play for this difference; while earthquakes occur deep underground, explosions are generally carried out at much shallower depths, which has an effect on how seismic waves propagate and interfere.

A simple, but very effective discriminant is the power spectrum spike. Figures 3, 4 and 5 illustrate that the typical earthquake series always has more than one spikes in the power spectrum, whereas the explosion series do not have consistently for more than one spike.

The algorithm implemented in this paper is different from "spectogram" and "fft" functions used with the SAC seismological program and MATLAB software. SAC computes a spectrogram using all of the data in memory. The spectra are calculated from a truncated autocorrelation function by using the maximum likelihood, maximum entropy, or Power Density Spectral statistical methodologies. The fft (in MATLAB) is a faster implementation of the DFT which reduces an $n$-point Fourier transform to about $(n/2)\log_{2}n$ complex multiplications from $n^2$ multiplications, providing an increase of speed for large data sets. However, the fft is not suitable for non periodic time series or non commensurate frequency (see subsection 3.1). The Dynamic Fourier methodology is a more reliable approach, since tapering is utilized to fix the spectral leakage by minimizing the effect of the discontinuity between the beginning and the end of the time series.

5.   Conclusion

In this study we investigate the use of the Dynamic Fourier method as a tool to discriminate between seismic signals generated by natural earthquakes and mining explosions. The earthquakes analyzed correspond to a group of small-magnitude aftershocks generated by a $5.2$ magnitude intraplate earthquake in Arizona. The mining explosions were selected such that they occurred very close to the earthquakes and were cataloged with similar magnitude as the earthquakes. We used seimograms from earthquakes and explosions recorded by a seismic station located about $160$ km west from the epicenters of these events. Our results suggest that $S$-waves frequency content can be used as a discriminant tool. In particular, a high-amplitude peak in the $S$-waves spectra for frequencies $\sim$ 6-8 Hz is present for the earthquakes but not for the explosions. Further analysis is suggested in order to investigate how epicentral distance and the magnitude of the events affect these results.

Fourier transform is widely used in digital signal processing for frequency analysis of both stationary and non stationary signals. In the case of transient signals, the characteristics of the signals to be analyzed are strongly time-dependent. This imposes restrictions on the use of the Fourier transform, since the recorded earthquakes and mining explosions do not necessarily fulfill the stationary condition. To overcome this issue we used a short-time Fourier transform where the input signal is a Fourier transform analyzed in a short time window which shifts along the input signal length with some overlap. In addition to dynamic Fourier analysis as a method to overcome the restriction of stationarity, researchers have also sought wavelet methods to discriminate seismic events ^[7].

Acknowledgments

The earthquake data used in this study were provided by the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC).

The authors would like to thank the reviewers for the careful reading of the manuscript and the fruitful suggestions that helped to improve this work.

Conflict of Interest

All authors declare no conflicts of interest in this paper.