# Bulletin of the Seismological Society of America

- Copyright © 2006 Bulletin of the Seismological Society of America

## Abstract

We describe an automated short-period Rayleigh wave (*Rg*) detector designed to work on local (<2.5° epicentral distance) events recorded at three- component stations. The detector was modeled after an automatic 17- to 22-sec Rayleigh-wave detection method; however, we have modified the algorithms for local distance and short-period applications. We have tested the detector on a well-located cluster of mining events from central India and on a set of ground-truth events in the area. The *Rg* detector was also integrated into a semiautomatic event detection and location algorithm and applied on continuous data. Fourier and wavelet-based methods are evaluated for prefiltering. We observe that sample standard deviations of backazimuth estimates using the *Rg* detector, after wavelet prefiltering, are comparable to *fk*3C *P* backazimuth estimates from event clusters. Our results indicate that using the *Rg*-phase backazimuths for event location is a promising alternative to using small signal-to-noise ratio first-arrival backazimuths. We recommend wavelet prefiltering versus Fourier prefiltering because it is more consistent for the detection of low signal-to-noise ratio events at local distances.

## Introduction

The detection and location of mining explosions, and their discrimination from shallow earthquakes and nuclear explosions in regions of nuclear monitoring concern, are important problems for seismologists. One useful monitoring discriminant is the presence of fundamental-mode short- period (<2.5 sec) crustal Rayleigh waves or *Rg* (Storchack *et al.*, 2003) on waveforms, which indicates shallow (less than 3–4 km) source depth (Kafka, 1990). The repeated detection of *Rg* (≤3 km/sec group velocity) during daytime hours is often an indication of mining operations. Our goal was to develop and test an effective automatic method to detect, identify, and characterize short-period *Rg* phases from small magnitude events observed at local stations.

A typical procedure for automatic shallow-event identification is the use of first-arrival detection algorithms (Goforth and Herrin, 1981; Allen, 1982; Baer and Kradolfer, 1987; Der and Shumway, 1999; Zhang *et al.*, 2003) followed by event location and characterization, with, as an example, waveform correlation methods (Harris, 1991; Bonner *et al.*, 2002). Using this type of method surface-wave arrivals can only be reported after the events have been located and processed by an analyst.

An alternate approach is to continuously scan the data for *Rg* arrivals without reference to known events. A large number of automatic surface-wave detectors have been developed (Cichowicz, 1993; Ruud *et al.*, 1993; Takanami and Kitagawa, 1993; Leonard and Kennett, 1999; Pinsky, 1999). These detectors have been developed at regional-to-teleseismic distances for 8- to 40-sec period Rayleigh waves on the basis of dispersion, polarization (Chael, 1997; Selby, 2001; Baker and Stevens, 2004) and phase-matched filters (Levshin and Ritzwoller, 2001). To our knowledge, our study is the first attempt to develop an automatic short-period (<2.5 sec) *Rg* detection algorithm. The detector presented in this study was inspired by the Chael (1997) automatic teleseismic Rayleigh-wave detection method. The algorithm has the advantage of continuously and independently monitoring all directions for *Rg* arrivals. It uses both amplitude and polarization information for *Rg* detection and determines the *Rg* backazimuth with a precision comparable to single-station *P* backazimuth estimates.

In this study we describe an automatic *Rg* detection algorithm and discuss Fourier and wavelet-based prefiltering methods; present results of detection tests on well-located events recorded at the hyb seismic station (Hyderabad, central India); compare estimates of *Rg* and *P* backazimuth and evaluate the automatic *Rg* detector performance when applied on continuous data.

## Method

### Rg Detection

The *Rg* detector algorithm is a modified version of the Chael (1997) detection method, developed for intermediate- period Rayleigh-wave arrivals (17–22 sec). Modifications are made to account for the shorter periods (<2.5 sec) characteristic of the *Rg* phase from local events. These modifications include (1) prefiltering with a Meyer (Fig. 1) continuous wavelet transform centered on an *Rg*, period characteristic to the study region; (2) calculating the covariance matrix for rotated waveforms; (3) estimating an optimal constant for the sta and lta recursive filters where sta is the short-term average and lta is the long-term average (Allen, 1982); and (4) using new weights of the detection function for planarity verification.

We developed the detection algorithm in two steps.

#### Step 1:

Prefiltering. We have observed that prefiltering is essential for accuracy in backazimuth estimation. We tested Fourier and wavelet filters on analyst-detected *Rg* arrivals from well-located clusters of events. Our goal was to find a simple, fast, reliable, and window-independent prefiltering procedure, suitable for an automatic algorithm at local distances. For this purpose, for central India we chose a Butterworth zero-phase (three poles, two passes) filter, between 0.4 and 1.3 Hz, determined to be optimal for *Rg* arrivals (J. L. Bonner, personal comm., 2002). Alternatively, we used the coefficients of a continuous wavelet transform (cwt) with the Meyer wavelet (Daubechies, 1992; Misiti *et al.*, 1997). The Meyer wavelet scale is chosen to correspond to a center period of 1 sec, as shown in the top and middle plots of Figure 1. The wavelet center period was chosen to be 1 sec, because in central India, Britton *et al.* (2003) observed that, for a distance range between 160 and 200 km, the *Rg* arrivals had periods between 0.9 and 1.7 sec.

The two prefiltering methods do not filter the waveforms in the same frequency range: the wavelet filter is more localized in frequency than the Fourier filter (Fig. 1, middle and lower plots). To make an exact comparison between results using wavelet and Fourier prefiltering, the filtered signal should have the same frequency content. The frequency content should be determined every time the cwt coefficients are calculated, because their values are proportional to the degree of similarity between the signal and the chosen wavelet. Because of its similarity to a sinc function, the frequency domain representation of a Meyer wavelet is a boxcar, starting at 0.55 Hz and ending at 2 Hz, centered on the analyst-chosen period of 1 sec (as shown in Fig. 1, middle plot). A Butterworth zero-phase (three poles, two passes) filter, whether between 0.4 Hz and 1.3 Hz (dotted line in Fig. 1, lower plot) or between 0.55 Hz and 2 Hz (continuous line in Fig. 1, lower plot) is not the best prefiltering choice for an exact comparison because it is less localized in frequency than the Meyer wavelet. For an exact comparison, a better choice of a Fourier-type filter would be, as an example, an autoregressive filter (Der and Shumway, 1999). Such filters could be used to duplicate the wavelet prefiltering results, when appropriate windows of noise and signal are chosen by an analyst or a picking algorithm. However, finding these filters was beyond the scope of our present project, as was an exact evaluation of one prefiltering method against the other. As shown below, both prefiltering methods produce good results for local events, however, they have different levels of consistency for the distance range in our study.

#### Step 2:

*Rg* Detection and Backazimuth Estimation. The *Rg* detection algorithm is applied to the three components of the data in a moving time window three times the predicted *Rg* period (1 sec in this India study), with a half- second step move-up, over a swath of backazimuths ranging between 0° and 360°, using an analyst-chosen step (2° in this study). A detection function is calculated as: 1 where *W _{i}* represent weights calculated from the sta values in each window.

For each backazimuth and timestep, the two covariance elements representing the cross-powers between vertical: radial (*σ*_{zr}) and vertical:transversal (*σ*_{zt}) are used as sta values in the moving time window. lta is calculated by using a similar formula in a window 10 times the predicted *Rg* period. The first weight, W1, is the square of the weight used by Chael (1997) to assess correlation between the vertical and radial motions. The second weight, W2, is applied to eliminate the time windows where the vertical and transverse components are correlated: 2

The third weight (W3) is the square of a planarity function as in Oonincx (1998) used to measure the degree of polarization of the three components in the vertical:radial plane. W3 is also applied to avoid cases when the amplitude of the transversal motion is very large.

For every detected *Rg* phase, values larger than 10% of the detection function's maximum value are plotted as a function of time and backazimuth. The detection surface is scanned for peaks exceeding 90% of the largest detection function value. The backazimuth is calculated by using a procedure described as the direction of the center of mass (Mardia, 1972), with squared values larger than 90% of the detection function as weights. Detections are declared when the maximum value of the detection function is larger than a threshold empirically estimated for each station. We choose the detection threshold at each station, for each prefiltering method, such that it has a lower value than the detection function value of at least 80% of the analyst-confirmed *Rg* arrivals.

### Event Location

The *Rg* detector is incorporated into an ensemble of routines for automatic analysis of continuous data, which also includes local and near-regional event location. Event epicentral distance is estimated from the *P* and *Lg* arrival- time difference using an empirical, station-specific formula: 3 where *d* is distance and Δ_{Lg-P} is the *Lg* and *P* arrival-time difference. The epicentral distance formula is determined by using a two-layer model: layer thickness *h*_{1} = 20 km, compressional seismic velocity *v*_{1} = 5.7 km/sec, and *h*_{2} = 20 km, *v*_{2} = 6.6 km/sec. The layer below the Moho has *v*_{3} = 8.1 km/sec. The *Lg* velocity is considered to be 3.5 km/sec. The velocity values are extracted from the WINPAK3D velocity model (Reiter *et al.*, 2005).

First-arrival backazimuths are estimated using the *fk*3C method, adapted from MatSeis (Harris and Young, 1997). Our modifications to the *fk*3C method for continuous data analysis include *f*-*k* analysis in a sliding 1-sec window with a step of a half-second, performed within 5 sec of the analyst-picked *P* arrival to obtain the most stable backazimuth estimate (Tibuleac *et al.*, 2004).

## Database

To examine the effects of prefiltering, we tested the *Rg* detector on 23 events in central India (Fig. 2) recorded by station Hyderabad, India (hyb; 17.4169° N, 78.5531° E) of the GEOSCOPE network between 21 April 1999 and 20 May 1999. These events were detected and assigned to four clusters using a waveform correlation method (Harris, 1991, Britton *et al.*, 2003).

We also tested the *Rg* detector on a group of four hyb ground-truth events (Table 1) from the EHB bulletin (Engdahl *et al.*, 1998), characterized by a ground-truth gt15 classification (15-km location accuracy), as defined by the International Association of Seismology and Physics of the Earth Interior (IASPEI) Working Group of Reference Events (http://lemond.colorado.edu/∼copgte). The detector was applied on 13 days of continuous hyb data available from GEOSCOPE between 24 April 1999 and 14 May 1999 and locations of 114 larger signal-to-noise ratio (snr) events with automatically detected *Rg* phases are represented as black dots in Figure 2.

## Results

### Prefiltering Test on Well-Located Events

We applied the automatic *Rg* detector after prefiltering the data from four clusters (23 events) recorded at station hyb in central India (Fig. 2), and we estimated the *Rg* backazimuth (Fig. 3) for Fourier prefiltering (*bazF*) and wavelet prefiltering (*bazW*).

Whereas *Rg* is detected as the most probable arrival for all 23 events using wavelet prefiltering, it is detected as the most probable arrival for only 21 of 23 (87%) events using Fourier prefiltering. Next we examine the low *snr* events that were not correctly detected with Fourier methods.

For event 21 (Fig. 3) from cluster 4, the *Rg* phase arriving at about 115 sec was not automatically detected after Fourier prefiltering (Fig. 4, bottom plot). The arrival that corresponds to the largest value of the detection function after Fourier prefiltering was another Rayleigh-type arrival at ∼440 sec, with a period of ∼2.5 sec and a backazimuth of 162° (Fig. 4, bottom plots). In the upper plots in Figure 4 we show that Meyer wavelet prefiltering (scale 14, corresponding to a 1-sec period) increased the *Rg* snr and resulted in the detection of the *Rg* phase.

The other event not detected by Fourier prefiltering was event 9 of cluster 3, for which the maximum of the detection function was below the empirical threshold value of 670 at hyb. Additionally, after Fourier prefiltering of event 16 of cluster 4, an *Rg* arrival was detected at the right time for event 16, but from a direction 40° off the true backazimuth. Wavelet prefiltering improved the detection results for these two events as well (see Fig. 3). However, for 20 events with larger snr, wavelet and Fourier prefiltering had similar results.

The total estimated *Rg* backazimuth sample standard deviation for the 23 events near hyb, using wavelet prefiltering, was 5.2°, equivalent to a maximum location error between 16 km (clusters 1, 3, 4) and 17 km (cluster 2). This standard deviation was better than the estimates using Fourier prefiltering (see very large Fourier prefiltering residuals in Fig. 3) and was close to the standard deviation of the first- arrival backazimuth residuals (3.1°) calculated by using data from a single station and wavelet prefiltering. Given the estimated location error, gt20, for these 24 events, we consider these single-station location results encouraging.

### Detector Performance on Ground-Truth Events

We tested the *Rg* detector on a set of four local gt15 ground-truth events (Table 1). Of these, we assumed that events with similar location and source mechanism would have similar low-frequency components of the vertical waveforms and be less affected by small-scale inhomogeneities in the upper crust. Results of waveform cross-correlation, after prefiltering from 0.01 Hz to 0.2 Hz (Tibuleac and Herrin, 2001), show that events 2, 3, and 4 are colocated (maximum of the normalized cross-correlation of 0.81), whereas event 1 has a different location (Table 1). We observe the most consistent estimates of the backazimuth (lowest sample standard deviation in Table 1) with the *fk*3C method (*P* backazimuth) followed by wavelet prefiltering (*bazW*) and finally, by Fourier prefiltering (*bazF*).

These results demonstrate that for local events, the performance of the *Rg* detector is more consistent when using wavelet decomposition as a prefiltering method and produces backazimuth estimates of a precision comparable to single-station *P* backazimuth estimates. A possible explanation is that the transient character and narrow frequency band of the *Rg* phases is best described by wavelet analysis (Tibuleac and Herrin, 1999). Based on these results, we used wavelet prefiltering as the primary prefiltering method for the *Rg* automatic detection algorithm at local distances in central India.

### Detector Performance on Continuous Data

The *Rg* detector was applied on 10-min waveform segments, with 1 min of overlap. The resulting detections were confirmed by an analyst, because no seismic bulletin was available.

#### Detection Results

The empirical detection threshold was chosen to be 2500 when using wavelet prefiltering. Detection results are presented in Table 2.

All the detected *Rg* arrivals were clearly visible. Nine percent of the automatic detections were false alarms. Most of the false alarms were spikes in the data that trigger the *Rg* detector. The *Rg*-type arrivals in the coda of large events automatically detected have not yet been added to the “false alarm” category, thereby raising to a total of 21% of the detections to be false alarms for the 2500 threshold, because the possibility that some of these phases are real *Rg* arrivals from local shallow events is still under investigation.

For most of the events with analyst-confirmed *Rg* arrivals not automatically detected, the *Rg* frequency was higher than 2 Hz, above the frequency range of our investigation. Also, the algorithm does not pick arrivals within 10 min of a larger, previously detected *Rg* phase, even if the detection function value is larger than the threshold.

Finally, 114 events with larger *snr*, all with automatically detected *Rg* arrivals, were chosen for a more detailed study (Fig. 2). All events were recorded between 1 a.m. and 2 p.m. UTC time, that is, during daylight hours (the local time is UTC plus 5.5 hr). The consistency of the recorded *Rg* arrival times supports the interpretation of most event sources as mine blasts. However, considering natural seismicity patterns in the area, some of the events are probably shallow earthquakes.

A plot of the estimated *fk*3C *P* versus *Rg* backazimuths for 114 events is presented in Figure 5. The mean residual backazimuth (*P* backazimuth − *Rg* backazimuth) is −1.6° with a standard deviation of 21.5°. The events with large differences between *P* and *Rg* backazimuth have low snr first arrivals. The right plot of Figure 5 represents the *Rg* and *P* arrival-time difference as a function of distance. Two regression lines were fit to the data, for both *Pg* and *Pn* first arrivals, and an *Rg* group velocity of 2.9 km/sec was determined to have the best L1 norm fit.

Detections with *Rg* detector function values lower than the threshold were considered “*Rg* noise.” Although the detection function values are approximately normally distributed in backazimuth, we have observed that the backazimuth distribution is not uniform (consistent with observations of Selby [2001] and Baker and Stevens [2004]). In Figure 6 we show that several peaks are observed: a less pronounced peak at 45° backazimuth, to the northeast; a peak to the southeast, with a maximum at 135° backazimuth; and a peak to the south-southwest, with a maximum at 200° backazimuth. We interpret these peaks to be caused by mining activity in the area.

## Conclusions

We have developed and tested an automatic *Rg*- detection method, with promising results, on a database of events with *Rg* arrivals and on four ground-truth (gt) events, all recorded at hyb, in central India. For 23 well-located events clustered near hyb, *Rg* is always detected after wavelet prefiltering. By using Fourier prefiltering between 0.4 and 1.3 Hz, *Rg* is detected as the most probable arrival for the 21 largest snr events of the 23 (87%). All gt event *Rg* arrivals are detected, with best backazimuth estimates obtained after wavelet prefiltering. Wavelet-based prefiltering of the local-event waveforms results in more consistent short-period surface-wave backazimuth estimates when used as an alternative to a Fourier filter.

The estimated *Rg* backazimuth sample standard deviation is 5.2°, comparable with the standard deviation of the first-arrival backazimuth residuals (3.1°). Typically, single- station locations use only the *P* backazimuth estimates; however, this study demonstrates that estimates of *Rg* backazimuth could also be used for location of events with low snr first arrivals and comparatively larger *Rg* amplitudes.

The algorithm is also tested on continuous data at station hyb (Hyderabad, India). Of the 207 analyst-identified *Rg* arrivals, 161 were automatically detected (78%) by using a 2500 detection threshold value. The undetected events were local *Rg* arrivals out of the chosen frequency range or events within 10 min of a detected *Rg* arrival. The false alarm rate was 9%.

Further improvements of the *Rg* detector will be focused on reducing the false alarm rate by decreasing the algorithm's sensitivity to spikes and high-amplitude teleseismic arrivals. Other improvements will include recording of all detections within the search window versus only the detection with the largest detection function value.

## Acknowledgments

We express our gratitude to Jessie Bonner, who had the idea to adapt the Chael (1997) method for *Rg* detection. We thank Delaine Reiter, Heather Hooper, Mark Leidig, and Anastasia Stroujkova for in-depth review and fruitful discussions. We thank James Lewkowicz for continuous support during this project. We express our gratitude to an anonymous reviewer for useful comments and suggestions.

- Manuscript received 18 March 2005.