| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Short Note |
1 Weston Geophysical Corporation
181
Bedford Street, Suite 1
Lexington, Massachusetts
02420
ileana{at}westongeophysical.com
brittonj{at}westongeophysical.com
 |
Abstract |
|---|
 Top Abstract IntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
|
Introduction |
|---|
|
TopAbstractIntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
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 |
|---|
|
TopAbstractIntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
|
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) |
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) |
Lg-P is the
Lg and P arrival-time difference. The epicentral distance
formula is determined by using a two-layer model: layer thickness
h1 = 20 km, compressional seismic velocity
v1 = 5.7 km/sec, and h2 =
20 km, v2 = 6.6 km/sec. The layer below the Moho has
v3 = 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 fk3C method, adapted from MatSeis (Harris and Young, 1997). Our modifications to the fk3C 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 |
|---|
|
TopAbstractIntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
|
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 |
|---|
|
TopAbstractIntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
|
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 fk3C
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 fk3C 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 |
|---|
|
TopAbstractIntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
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 |
|---|
|
TopAbstractIntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
Manuscript received March 18, 2005
|
References |
|---|
|
TopAbstractIntroductionMethodDatabaseResultsConclusionsAcknowledgmentsReferences |
|---|
Allen, R. (1982). Automatic phase pickers: their present
use and future prospects, Bull. Seism. Soc. Am.72
,S225
–S242.
Baer, M., and U. Kradolfer (1987). An automatic phase
picker for local and teleseismic events, Bull. Seism. Soc.
Am. 77,1437
–1445.
Baker, G. E., and J. L. Stevens (2004). Back azimuth estimation reliability using surface wave polarization, Geophys. Res. Lett. 31,L09611 , doi 10.1029/2004GL019510.[CrossRef]
Bonner, J. L., D. B. Harris, S. Rieven, I. M. Tibuleac, and J. M. Britton (2002). Development of a mining explosion database for Southern Asia, in Proceedings of the 24th Seismic Research Review—Nuclear Explosion Monitoring: Innovation and Integration, LA-UR_02_5048, vol. 1,211 –219.
Britton, J. M., D. Harris, I. M. Tibuleac, and J. L. Bonner (2003). Detection of mining explosions for a southern Asia seismic database, in Proceedings of the 25th Seismic Research Review—Nuclear Explosion Monitoring: Building the Knowledge Base, LA-UR_03_6029, vol. 1,211 –219.
Chael, E. P. (1997). An automated Rayleigh-wave
detection algorithm, Bull. Seism. Soc. Am.87
,157
–163.
Cichowicz, A. (1993). An automatic S-phase picker, Bull. Seism. Soc. Am.83 ,189 .
Daubechies, I. (1992). Ten Lectures on Wavelets, Society for Industrial and Applied Mathematics, Philadelphia, Pennsylvania, 117.
Der, Z. A., and R. H. Shumway (1999). Phase onset time estimation at regional distances using the CUMSUM algorithm, Phys. Earth Planet. Interiors 113,227 –246.[CrossRef]
Engdahl, E. R., R. van der Hilst, and R. Bulland (1998).
Global teleseismic earthquake relocation with improved travel time and procedure
for depth determination, Bull. Seism. Soc. Am.88
,722
–743.
Goforth, T., and E. Herrin (1981). An automatic signal
detection algorithm based on the Walsh transform, Bull. Seism. Soc.
Am. 71,1351
–1360.
Harris, D. B. (1991). A waveform correlation method for
identifying quarry explosions, Bull. Seism. Soc. Am.81
,2395
–2419.
Harris, M., and C. Young (1997). MatSeis: a seismic GUI and toolbox for MATLAB, Seism. Res. Lett.68 ,267 –269.
Kafka, A. (1990). Rg as depth discriminant for
earthquakes and explosions: a case study in New England, Bull. Seism.
Soc. Am. 80,373
–394.
Leonard, M., and B. L. N. Kennett (1999). Multi-component autoregressive techniques for the analysis of seismograms, Phys. Earth Planet. Interiors113 ,247 –263[CrossRef]
Levshin, A. L., and M. H. Ritzwoller (2001). Automated detection, extraction, and measurement of regional surface waves, Pure Appl. Geophys.158 ,1531 –1545.[CrossRef]
Mardia, K. V. (1972). Statistics of Directional Data, Academic Press, New York.
Misiti, M., Y. Misiti, G. Oppenheim, and J. Poggi (1997). Wavelet Toolbox, For Use with Matlab, © MathWorks, Inc., Natick, Massachusetts.
Oonincx, P. J. (1998). Automatic Phase Detection in Seismic Data using the Discrete Wavelet Transform, Report PNA-R9811, ISSN 1386- 3711, Centrum woor Wiskunde en Informatica (CWI), Amsterdam, The Netherlands.
Pinsky, V. (1999). An approach to automatic location of regional events, Phys. Earth Planet. Interiors113 ,293 –302.[CrossRef]
Reiter, D., W. Rodi, and M. Johnson (2005). Development
of a tomographic upper-mantle velocity model beneath Pakistan and northern India
for improved regional seismic-event location, Bull. Seism. Soc.
Am. 95,926
–940.
Ruud, B. O., C. D. Lindholm, and E. S. Husebye (1993).
An exercise in automating seismic record analysis and network bulletin
production, Bull. Seism. Soc. Am.83
,660
–679.
Selby, N. D. (2001). Association of Rayleigh waves using
back azimuth measurements: application to test ban verification,
Bull. Seism. Soc. Am.91
,580
–593.
Storchack, D. A., J. Schweitzer, and P. Bormann (2003). The IASPEI standard seismic phase list, Seism. Res. Lett. 74,761 –772.
Takanami, T., and G. Kitagawa (1993). Multivariate time-series model to estimate the arrival times of S-waves, Comput. Geosci. 19, no. 2, 295–301.
Tibuleac, I. M., and E. T. Herrin (1999). An automatic method for determination of Lg arrival times using wavelet methods, Seism. Res. Lett.70 ,577 –595
Tibuleac, I. M., and E. T. Herrin (2001). Detection and location capability at NVAR for events on the Nevada test site, Seism. Res. Lett.72 ,97 –107.
Tibuleac, I. M., J. Britton, D. B. Harris, T. Hauk, H. Hooper, and J. L. Bonner (2004). Detection methods for mining explosions in Southern Asia, in Proceedings of the 26th Seismic Research Review—Nuclear Explosion Monitoring: Building the Knowledge Base, Orlando, Florida, 21–23 September,367 –376.
Zhang, H., C. Thurber, and C. Rowe (2003). Wave arrival
detection and picking with multiscale wavelet analysis for single-component
recordings, Bull. Seism. Soc. Am.93
,1904
–1912.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
