The detection of low magnitude seismic events using array-based waveform correlation

Summary

It has long been accepted that occurrences of a known signal are most effectively detected by cross-correlating the incoming data stream with a waveform template. Such matched signal detectors have received very little attention in the field of detection seismology because there are relatively few instances in which the form of an anticipated seismic signal is known a priori. Repeating events in highly confined geographical regions have been observed to produce very similar waveforms and good signals from events at a given site can be exploited to detect subsequent co-located events at lower magnitudes than would be possible using traditional power detectors. Even greater improvement in signal detectability can be achieved using seismic arrays; running correlation coefficients from single sensors can be stacked over an array or network to result in a network correlation coefficient displaying a significant array gain. If two events are co-located, the time separating the corresponding patterns in the wave train as indicated by the cross-correlation function is identical for all seismic stations and this property means that the correlation coefficient traces are coherent even when the waveforms are not. We illustrate the power of array-based waveform correlation using the 1997 August 16 Kara Sea event. The weak event that occurred 4 hr after the main event was barely detected using an STA/LTA detector on the SPITS array but is readily detected by signal matching on a single channel. The main event was also recorded by the far more distant NORSAR array but no conventional detection can be made for the second event. A clear detection is, however, made when the correlation coefficient traces are beamformed over all sensors of the array. We estimate the reduction in detection threshold of a test signal on a regional seismic array using waveform correlation by scaling down a master signal and immersing it into seismic noise. We show that, for this case, waveform correlation using a single channel detects signals of approximately 0.7 orders of magnitude lower than is possible using an STA/LTA detector on the array beam. Waveform matching on the full array provides an additional improvement of approximately 0.4 magnitude units. We describe a case study in which small seismic events at the Barentsburg coal mine on Spitsbergen were detected using the signals from a major rockburst as master waveforms. Many spurious triggers occurred in this study whereby short sections of signal exhibited coincidental similarity with unrelated incoming wave fronts. We demonstrate how such false alarms can almost always be identified and screened out automatically by performing frequency–wavenumber analysis upon the set of individual correlation coefficient traces.

detection threshold, matched signal, seismic array, waveform correlation

1 Introduction

The most effective method of detecting a known signal in a potentially noisy time-series is to cross-correlate a waveform template with successive time segments of incoming data. Any segments of the continuous data-stream which display a high degree of similarity to the template or master waveform will result in a high value of the correlation function. Such a procedure is referred to as a matched filter or a matched signal detector (see e.g. van Trees 1968). The detection by waveform correlation of a synthetic signal buried in noise was demonstrated by Anstey (1966) in a paper that provided a very comprehensive literature review of the early applications of signal correlation in exploration geophysics and other fields.

Such methods have so far received very little attention in the field of detection seismology for the simple reason that they are extremely sensitive to the form of the master waveform applied, and the overwhelming majority of signals detected by the worldwide network of seismometers come from unknown sources and consequently have unknown waveforms. Unknown signals are typically detected by a power detector of the kind proposed by Freiberger (1963) whereby the power over a short time-window (the short-term average, STA) is compared with the power over a long time-window (the long-term average, LTA) for waveforms filtered in a range of different frequency bands. A signal detection occurs when the STA/LTA ratio (the signal-to-noise ratio, SNR) exceeds a predetermined threshold. Sudden increases in the energy over a given frequency range can also be detected by multiscale (wavelet decomposition) methods (see e.g. Oonincx 1999).

The detection of weaker signals often requires array-processing techniques whereby a systematic delay-and-stack (beamforming) of traces from closely spaced instruments increases the SNR through a simultaneous summation of coherent signal and cancellation of incoherent noise (e.g. Schweitzer 2002; Rost & Thomas 2002). Cross-correlation of coherent signals from adjacent instruments allows for the estimation of highly accurate relative delay times from which we can estimate phase propagation parameters (VanDecar & Crosson 1990; Cansi 1995; Almendros 2004) and epicentre locations (see e.g. Almendros 1999).

The semblance of waveforms from a single source recorded at two closely spaced instruments depends upon the distance between the sensors and the homogeneity of the underlying rock; the frequency above which signals become incoherent decreases rapidly as the interstation distance increases. The NORSAR array in southern Norway (Bungum 1971) was designed to detect low-yield underground nuclear explosions at teleseismic distances. The current array consists of 42 sites divided into seven subarrays each consisting of six sites; adjacent instruments within the subarrays are typically separated by approximately 2 km and the full array has a diameter of approximately 80 km. The spacing of the instruments was chosen such that the coherence of the microseismic background noise would be a minimum allowing an optimal SNR-gain for teleseismic signals with frequency of 1–2 Hz for beams with suitable steering parameters. The array provides excellent resolution for parameter estimation of teleseismic phases at the expense of signal coherency for regional signals which, dominated by frequencies above 4 Hz, are not even coherent over the subarrays. The emphasis in nuclear explosion monitoring is currently on the detection, identification and association of regional phases using a denser network of regional arrays. These arrays typically consist of between nine and 25 sites within a diameter of not more than about 2 km. This allows for coherency for signals of far higher frequency, albeit at the expense of resolution capability for teleseismic phases.

Reciprocally, two waveforms recorded at the same instrument are likely to show great similarity if and when the corresponding events occurred within very close proximity of each other and were associated with similar source mechanisms. Geller & Mueller (1980) recognized a set of events that resulted in almost identical waveforms when filtered in a sufficiently low-frequency band, concluding that the source locations of the events could not be separated by more than a quarter of the dominant wavelength. Waveform correlation has been used to measure very small differential traveltimes from event pairs recorded at a network of stations to provide highly improved relative location estimates of nuclear explosions (Shearer & Astiz 1997; Thurber 2001; Phillips 2001; Waldhauser 2004) and earthquakes (e.g. Poupinet 1984; Shearer 1997; Astiz & Shearer 2001; Shearer 2003; Du 2004). The double difference (DD) techniques (e.g. Waldhauser & Ellswoth 2000; Waldhauser & Ellsworth 2002; Schaff 2002; Zhang & Thurber 2002) have vastly improved event location and tomographic modelling through the combined use of absolute traveltime measurements and highly accurate cross-correlation differential traveltimes. Menke (2001) points out that the correlation coefficient itself, and not only correlation-derived delay times, can be useful in constraining earthquake locations.

The degree of waveform similarity between two events has been used repeatedly to classify events according to source location. Israelsson (1990) identified clusters of mining events that indicated high correlation values at frequencies above 15 Hz using the NORES array. The quarter wavelength argument of Geller & Mueller (1980) would constrain these event epicentres to within 100 m, far closer than the accuracy that could be anticipated from standard array processing. Harris (1991), Rivière-Barbier & Grant (1993) and Schulte-Theis & Joswig (1993) all subsequently used cross-correlation techniques to identify events from repeating industrial sources from mines and quarries; the automatic classification of such events with a high degree of confidence is central to the topic of nuclear explosion monitoring in order that analyst time is not wasted upon their identification. However, due to complicated ripple-firing mechanisms (and therefore highly varying source-time functions) and strong heterogeneity within the sites, seismic signals from subsequent mining explosions have often been found to correlate poorly in spite of the close proximity of sources (Bonner 2003; McLaughlin 2004). Poor correlation has motivated the development of signal-subspace detectors for events in a given source-region (Harris 1997), generalizing the principle of rank-1 matched signal detectors (correlators) to multirank matched subspace detectors (Scharf & Friedlander 1994).

Other recent studies have indicated that the applicability of matched signal or correlation detectors may be greater than previously assumed, both for natural and man-made seismicity. Schaff & Richards (2004a, b) undertook an extensive cross-correlation study of over 14 000 seismic events in China as listed in the Annual Bulletin of Chinese Earthquakes (ABCE). They found that 1301 (or approximately 10 per cent) were repeating events with epicentres not more than 1 km apart; the majority of the repeating events occurred within one month of each other, indicating a causal relationship. Surprisingly, even the highly scattered Lg wave was found to correlate well for many event clusters. The resulting database was used to achieve an impressive lateral location precision of about 100–300 m. Of arguably even greater importance, an automated cross-correlation detector was able to detect all of these 1301 events and possibly many more of lower magnitude than those included in the ABCE. Many of the 1301 events were missed in the both the Reviewed Event Bulletin (REB) of the prototype International Data Center (pIDC) of the CTBTO and the bulletin of the International Seismological Centre (ISC).

NORSAR has recently applied waveform correlation techniques with specific reference to seismic arrays in order to detect small, cavity-decoupled, underground explosions in Central Sweden (Gibbons & Ringdal 2004; Stevens 2004). The largest of these events were recorded with a high SNR on each of the NORES, NORSAR and Hagfors seismic arrays at a distance of approximately 150 km; one such event was selected and the resulting signals were filtered in a frequency band for which the SNR was optimal. A template waveform was extracted for each individual trace and these templates were subsequently correlated with array data over time-periods in which the undetected events were assumed to have taken place. The correlation coefficient channels from the single traces were stacked to give an array correlation trace and, using an ordinary power detector upon this ‘beam’, all eight of the reported explosions were identified. Only one of these events had been large enough to be detected using the standard processing and the smallest detonation (500 Kg TNT in a 1000 cubic meter chamber, a charge density 20 times smaller than that of the master event) was not visible in the waveform data in any available frequency band; the event was only detectable by a summation of multiple cross-correlation traces. Having correlated the master event template with over 2 yr of array data, it could be confirmed that the detector had not triggered once without there having been a confirmed explosion at this site. Given that every explosion at the site that the authors were aware of, with a charge density exceeding 0.5 kg m ⁻³, had been detected using this method, we conclude that this method constitutes a very sensitive, source-specific, detector with an exceedingly low false alarm rate.

The events detected by Gibbons & Ringdal (2004) had very carefully controlled, simple, explosive sources and were known to have taken place within one of three chambers all located within 200 m of each other. In the current paper, we detail and expand upon the method employed by Gibbons & Ringdal (2004) and demonstrate its application to the detection of more general low-magnitude seismic events. In Section 2, we formulate the method for both signal detection and, in the case of a good correlation, the estimation of signal amplitude. We demonstrate the detection using waveform correlation of the 1997 August 16 Kara Sea event aftershock, using the signal from the main event as a template, on stations where traditional array-processing methods either failed to detect or only barely detected the signal. In Section 3, we estimate the gain in detection performance using waveform correlation as opposed to a standard power detector on a seismometer array by scaling down a signal into seismic noise and examining how effectively the different methods manage to detect the submerged signal. In Section 4, we explain a common source of spurious triggers encountered when attempting this procedure on a small-aperture regional array and demonstrate a method by which the majority of such false alarms can be identified and filtered out automatically. In Section 5, we present the findings of a study to detect small-scale seismic events occurring within a very confined region surrounding the source of a rockburst at the Barentsburg coal mine on the Island of Spitsbergen in the summer of 2004. The purpose of this study was to identify as many events as possible which were likely to have originated close to the site of the main rockburst. Finally, Section 6 presents our conclusions from the study and discusses various aspects of applying the array-based correlation technique in practical applications.

2 Waveform Correlation as a Tool for Detection

2.1 Formulation

We use here the notation w_N,Δt(t₀) to denote the vector of N consecutive samples of a non-zero time-series w(t), where t₀ is the time of the first sample and Δt is the spacing between samples:

1

The inner product between v_N,Δt(t_v) and w_N,Δt(t_w) is defined by

2

and the fully normalized cross-correlation coefficient by

3

The coefficient C will always lie in the interval [−1, 1] with the extreme values occurring only when one of the time-series is an exact multiple of the other, i.e.

4

with the sign of C being the same as the sign of α. With real-world data, ± unity will almost never be achieved and the correlation coefficient will fall somewhere between; a high value of C indicates a high degree of waveform similarity and a low value indicates little similarity.

If w(t) denotes the data recorded on a single seismometer channel and t_M the starting time of a data window containing a signal of interest (a template or master waveform) then

5

is a well-defined function of time which measures the similarity between the waveform immediately following time t and that of the template waveform. That C_w(t)_N,Δt can be used by a detector for waveforms that resemble the master signal is clear; what constitutes a detection on C_w(t) is less straightforward. The significance that can be attached to a high correlation coefficient depends upon the nature, primarily the length and complexity, of the template waveform. If w(t) is bandpass filtered in a narrow frequency band and the template signal is short and monochromatic (i.e. has a low time-bandwidth product) then it is likely that C_w(t) will frequently and trivially reach a high value. If the template waveform has a high time-bandwidth product, then it is likely that a high correlation coefficient will be of somewhat greater significance.

The value of the correlation coefficient is not necessarily high when a waveform corresponding to the template is detected; a low correlation coefficient can be significant if it is substantially greater than the values obtained using the same template on surrounding segments of data. To facilitate the automatic detection of significant values of the correlation coefficient, we define a function (analogous to an SNR) which measures the correlation at a given time t relative to the values obtained at times in the vicinity of t. Throughout this paper, we refer to this as a scaled correlation coefficient and apply the definition

6

where C_RMS(t) is given by

7

and the interval I_t is defined as the union of two time windows, one prior to and one following time t

8

with 0 < a < b. The above definition makes C′_w(t)_N,Δt insensitive to the value of the correlation coefficient in the immediate vicinity of t; strong ‘side-lobes’ may be observed in the correlation coefficient since a waveform may resemble itself shifted by a small number of cycles far more than an entirely unrelated segment of data. The values for a and b are rather empirically determined on a case by case basis, in the same way as the window lengths for the STA and LTA are chosen to suit the signal frequency content for a power detector. The value of a must be sufficiently large that the interval I_t for the time of maximum correlation does not include significant peaks in the autocorrelation function. The value of b is set subsequently to provide a data window long enough to provide a relatively stable estimate of the variation anticipated in the correlation function in the absence of close matches with the waveform template; this will vary according to the frequency content of the signal. Typical values of the constants a and b used in this paper are a= 1.0 s and b= 2.5 s. A simple STA/LTA detector could be used upon the correlation coefficient traces, but using a definition like that given in eq. (6) avoids an unnecessary smoothing of the correlation trace.

The single channel cross-correlation coefficient defined in eq. (5) can be generalized to an array or network of M different channels. If the data on channel j is denoted w_j(t) then the fully normalized cross-correlation coefficient for channel j is given by

9

The reference time for the master event, t_R, is fixed for all channels in the array or network, which we will denote by . For the sake of generality, we provide each channel with a separate time offset, τ_j, such that the master waveform for channel j begins at a time (t_R+τ_j). It must be emphasized that the offsets τ_j can be quite arbitrarily set; small changes in τ_j should not make significant changes to the correlation coefficient. (This is in contrast to traditional beamforming where the steering of the beam is controlled entirely by the specified time delays.) We also allow the length of the master time window and the sampling frequency to vary from channel to channel, hence the subscript j on the sample count, N, and the sampling interval, Δt. The formulation specified in eq. (9) means that even when the master event time-windows at the different stations are separated by long time-intervals, all of the correlation functions C_j(t) should attain a maximum at the same time, t, if a signal originating from the same location as the source of the master signal is detected. We can define an array or network correlation coefficient

10

which should achieve high values when and only when the individual C_j(t) channels interfere constructively. We can similarly form a scaled correlation coefficient, , analoguous to that defined in eq. (6).

Because we have respected the possibility of arbitrarily large time-differences between arrivals on the different channels, an array or network correlation coefficient could in principle be calculated for a highly heterogeneous network. However, in such cases, great care must be taken to ensure that the individual terms in eq. (10) are evaluated at precisely the same time; this is clearly easiest to enforce when the single channel correlation traces, C_j(t) , are evaluated for an array with central timing and uniform sampling rate.

Amplitude measurements of signals with very low SNR are likely to be poor due to the high contribution from the background noise. If such a signal is detected due to a high correlation coefficient with a high-SNR signal template, it may be possible to obtain a far better estimate of the amplitude by finding an optimal scaling factor between the master and detected signals. A more general form of eq. (4) is

11

where x is the data vector comprising the signal templates for all of the channels used, y is the corresponding data vector at the time of maximum cross-correlation and ɛ is the component of y which is not the direct result of a signal of the form x. ɛ vanishes in the limit of perfect correlation. (Note that we assume that our template waveform accurately describes the signal we are attempting to detect and do not make an explicit allowance for a noise component of the vector x.) The crudest estimate of the optimal scaling factor, α, is given by

12

where we have assumed that the noise, ɛ, has a zero correlation with the signal template, x. However, such a formulation does not take into account the quality of each measurement; the lower the amplitude of a sample in the waveform template, the more likely it is that the corresponding sample in the detected waveform is dominated by noise.

If our signal template uses data from M channels, with channel j containing N_j samples, then the length of the vectors x, y and (the unknown quantity) ɛ is . Since we know neither the form of the noise, ɛ, or the scaling factor, α, then eq. (11) is underdetermined. (α= 0 is a perfectly valid solution, in which the vector of observed data, y, is independent of the template, x.) However, assuming that y is deterministically caused by a signal of the specified form and that ɛ is Gaussian noise, eq. (11) is equivalent to an overdetermined system of N linear equations that can be written in the form

13

where we adopt the notation of Gubbins (2004). The vector of model parameters, m, (see Gubbins 2004, for definitions) is our vector of unknowns; it contains only the scalar α and therefore has length 1. The measurement vector, d, is equivalent to our data y, the N× 1 coefficient matrix A is equivalent to our signal template, x, and the unknown vector of errors, e, corresponds to our noise vector, ɛ. The normal equations that minimize the vector e, with equal weighting given to each measurement, lead to the solution

14

which is, on inspection, identical to the expression in eq. (12).

Gubbins (2004) demonstrates how, if the errors in eq. (13) have zero mean and an N×N covariance matrix C_e, the least-squares model parameter vector, m, is given by

15

For our system, this gives α in terms of our master signal and detected waveform using

16

where W is a matrix of weights. A natural assumption is that W_ij= 0 for i≠j, meaning that W is a diagonal matrix with W_ii simply indicating the significance of measurement i in the inversion. W_ii could, for example, be set proportional to the absolute value of x_i; this would ensure that the parts of the master signal with the highest amplitude would be given greatest weight in the inversion. However, we settled on an iterative procedure where an initial estimate α⁽⁰⁾ was obtained using W=I, i.e.

17

and subsequent estimates of α were obtained using a modified weighting matrix:

18

The weighting matrix for iteration k is determined from the residual from the previous evaluation of α with

19

where the residual r^(k−1)_i is given by

20

and ε is a very small number chosen to prevent division by zero or numerical rounding errors. Successive estimates, α^(k), are evaluated using the recipe above until the iteration either converges to a specified precision or a predetermined number of iterations is exceeded without convergence.

This is an example of an iteratively reweighted least squares or IRLS procedure and details and references are provided by Gubbins (2004). The procedure clearly gives the most meaningful estimates of a scaling factor when the residuals are small (i.e. when the correlation coefficient is high) and the inability of the iteration to converge on a value is a sure indication that the correlation is spurious or simply too poor for such an estimation of signal amplitude. All scaling factors quoted in this article are calculated using this procedure.

2.2 Example: The 1997 August 16 event near Novaya Zemlya

In this section, we illustrate the ability of waveform correlation to detect signals from a weak seismic event that has occurred in the near vicinity of a seismic disturbance for which master waveforms exist. We also demonstrate the ability of a multichannel matched filter detector implemented on a large aperture seismic array to detect a signal that would certainly have been missed had only a single channel at the same distance been available.

The example illustrated here has special significance in the field of the seismic monitoring of of the Comprehensive Nuclear Test Ban Treaty (CTBT). Early on 1997 August 16, a small seismic event took place in the Kara Sea, approximately 100 km from the former Soviet nuclear test site on the island of Novaya Zemlya (see Fig. 1). The close proximity of this event to the test site led to initial concerns that the event could have been a small clandestine nuclear explosion in violation of the CTBT that had been signed by Russia 11 months previously. The event has been the subject of many publications (e.g. Richards & Kim 1997; Hartse 1998; Ringdal 1997; Asming 1998; Ringdal 2002; Bowers 2001; Kremenetskaya 2001a; Bowers 2002; Schweitzer & Kennett 2002) and the generally accepted conclusion, based upon location, spectral characteristics and other observations, is that the event was in fact an offshore earthquake.

Figure 1

Map indicating the locations of the Spitsbergen and NORSAR seismometer arrays relative to the location of the 1997 August 16 seismic disturbance in the Kara Sea, approximately 50 km to the east of the island of Novaya Zemlya. The locations of the ARCES array, which was not operational at the time of the event, and the station at Amderma, Russia, are also marked.

One of the most compelling pieces of evidence for the classification of the event as an earthquake was the occurrence of a small event (presumed to be an aftershock) approximately 4 hr following the main event (see, in particular, Richards & Kim 1997). Careful analysis of regional seismograms of both events, especially from the nearby Amderma station (Ringdal & Kremenetskaya 1999; Ringdal 2002), indicate that the events were probably from the same location. Further evidence that the two events were co-located has been presented by Gibbons & Ringdal (2005b) who showed, by expanding upon the waveform correlation procedure in the present paper, that the time delay between P-phase onsets for the two events was the same (to within 0.1 s) for each of the three stations SPITS, NORSAR and Amderma. Bowers (2002) found that the main event produced signals consistent with a double-couple source. A source mechanism for the aftershock has not been determined due to the small number of observations and the low SNR of the signals where they were observed.

The main event, of magnitude m _b= 3.5 (Ringdal 1997), was recorded with a high SNR for both P- and S-phases at the 9-site Spitsbergen array. Only a weak P-phase from the aftershock was detected at SPITS; this was in fact the only detection of this event on any station of the International Monitoring System (IMS) for the CTBTO. Waveforms from the master event were extracted and filtered in the 4.0–8.0 Hz frequency band, resampled to a frequency of 80 Hz, and a-minute-long section of data was cut with a starting time 02.13.44.915. It is reasonable to extract an identical time-window for every channel given that, over the small aperture of the SPITS array, all seismic phases will arrive at all sites within less than a second of each other. The time delays between the phase arrivals at the different sites do not constitute a problem; as we are attempting to detect another seismic event from the same source location, the sought signal will be associated with identical delay times at the same receiver sites. This waveform template was correlated with SPITS data over several hours both prior to and following the main event and a maximum of the array correlation coefficient beam was achieved at a time 1997-228:06.21.55.815, obtained from a spline interpolation of the discrete time-series (see Fig. 2). The correlation results indicate that the origin time of the second event was 14890.9 s following the origin time of the main event.

Figure 2

Detection of an aftershock from the 1997 August 16 Kara Sea event using waveform correlation on the short-period vertical channels of the Spitsbergen array. Each channel was bandpass filtered between 4.0 and 8.0 Hz, and a 60 s long data segment was extracted from the master event signals, labelled ME above, beginning at a time 1997-228:02.13.44.915. Data containing the presumed aftershock was filtered in the same band and a trace of fully normalized correlation coefficients was calculated for each channel. Above the single-channel correlation coefficient traces is the zero-delay array correlation beam and the top trace is the scaled beam, . A clear peak is observed on the beams at a time 1997-228:06.21.55.815. Waveforms from the main event and aftershock are displayed (for three channels only) aligned according to the maximum value of the correlation beam.

The value of the correlation coefficients in Fig. 2 is high despite the low SNR of the aftershock. The correlation coefficient for each individual channel lies between 0.5 and 0.7, with the array correlation beam, , attaining a maximum value of 0.615. The scaled correlation coefficient, , attains a value of 15.1 at this time, compared with typical values in the approximate range ±4 in the absence of a clear correlation. The maximum correlation coefficient is achieved simultaneously for all channels and the peaks in the single channel correlation functions are clear maxima which could be detected by a simple threshold rule or power detector. We conclude that the use of a seismic array was not necessary for the detection of this event at the Spitsbergen site; using waveform correlation for a single sensor would have sufficed. However, the synchronicity of the correlation peaks is a useful observation since it supports the assumption that the slowness vectors for phase arrivals from the two events are identical (Almendros 2004). The master-event signals recorded at the different sites over this small-aperture array show a high degree of semblance. This is the reason that the P-phase from the aftershock could be detected by simple beamforming; the appropriate delay-and-stack of traces provides a significant improvement in the signal-to-noise ratio. Nevertheless, a weak P-phase detection at one regional array is usually not sufficient to provide a confident location estimate. It is also worth noting how much more accurately the arrival time can be estimated as a result of the correlation with the master event waveform.

An even more impressive observation is that the small aftershock could also be detected by waveform correlation at the NORSAR array at a distance of approximately 2300 km. This array had a clear detection of the main event; the signals recorded at three of the array sites are displayed in Fig. 3. The nature of the signal is very different at this much larger distance; due to the attenuation of energy at higher frequencies, the optimal SNR is obtained in a somewhat lower frequency band. Note also, for the frequency band displayed, that the waveforms at the different sites are highly dissimilar. The signal from the aftershock is not visible for any of the single channels of the NORSAR array and no detection was made. Segments of the waveforms from the master event were extracted and cross-correlated with the filtered data from the time of the aftershock. Inspection of the single-channel correlation traces in Fig. 3 indicates no clearly discernible peak values. However, a zero-delay stacking of all of the available correlation channels results in an array correlation beam with a clear maximum at a time 1997-228:06.23.49.999, 14890.9 s following the reference time for the main event. It is noteworthy that this time delay is, to the 0.001 s accuracy aimed for in the interpolation procedure, the same as the time delay between the two events previously estimated by using SPITS array data.

Figure 3

Detection of the Kara Sea aftershock by waveform correlation using the NORSAR array. The frequency band applied in this calculation is 2.5–8.0 Hz. The time windows containing the master event signal are staggered by several seconds to account for the significant time delays across the array; the first master event time window begins at 1997-228:02.15.39.087 for the instrument NC301. Data were present for 36 out of a possible 42 channels and the 36 waveform correlation traces were stacked as shown; the beam displays a clear peak at a time 1997-228:06.23.49.999.

The maximum value of the NORSAR array correlation beam is 0.086, hardly a convincing indication of waveform semblance. However, the scaled coefficient (displayed in the top trace of Fig. 3) attains a value of 11.88, confirming a significantly higher degree of waveform similarity at this particular time. Not only is the signal from the aftershock (estimated magnitude m_b= 2.6) buried deep within the noise at NORSAR, but the single-channel correlation traces do not, in general, register even a local maximum value at this time. The failure of conventional beamforming over arrays to result in an SNR gain is almost invariably the result of incoherence between the signals at the different receiver sites. Given the selection of appropriately defined waveform templates [see eqs (9) and (10)], and the occurrence of signals from two co-located seismic sources, the ‘true cross-correlation peak’ will occur simultaneously for all sites regardless of whether or not the waveforms themselves are coherent at the different receivers. Consequently, beamforming the correlation coefficient traces over the large-aperture NORSAR array results in a spectacular SNR gain; local maxima of the single-channel correlation traces which result from coincidental similarity of unrelated seismic noise cancel out under the stacking operation, leaving only a superposition of the correlation maxima which result from the same deterministic waveform similarity.

The examples presented in this section clearly show the power of the array-based correlation technique. In conventional array processing, a requirement for array gain is incoherent noise and coherent signal. In contrast, array-based waveform correlation requires only incoherent noise in order to be successfully applied. The requirement for waveform coherency across the array is replaced by a requirement for coherency between the master waveform and the target waveform. The ‘signals’ in this case are the cross-correlation traces for each array sensor, and the peaks of these traces occur at exactly the same time when the master event and the target event are co-located. Therefore, an infinite-velocity array beam based on the correlation traces can be calculated with no loss due to missteering or lack of signal coherence.

3 A Comparison Between Detection using Waveform Correlation and Standard Regional Array Processing: A Case Study using Scaled-Down Signals Immersed in Noise

In the previous section it was demonstrated that, given a suitable waveform template, a low SNR signal can be detected easily using waveform correlation on a single channel for an event that was only barely detected by a standard STA/LTA detector on a full array beam (Fig. 2). It was also demonstrated that the same event could be detected on a far more distant array by the beamforming of the individual channel correlation traces when the signal had not been detected by either standard array processing or by single-channel waveform correlation (Fig. 3). In the current section, we aim to quantify the improvement in detection capability for a typical regional seismic array facilitated by the use of full-waveform matching.

We extract a high quality regional signal recorded at a chosen seismic array and simulate signals from hypothetical events of lower magnitude, from the same location, by scaling down the master signal and immersing the resulting waveforms into seismic noise. For a given scaling factor (i.e. magnitude of the hypothetical event) we seek to determine the probability of detecting the event using an STA/LTA detector on an optimally steered beam (using the full array configuration), and also using waveform correlation for many different station configurations including single channel, 3-component (single station), full-array and various logical subsets of sensors.

A linearly scaled waveform only approximates the signal that would be generated by a co-located event of lower magnitude; the power spectra of the generated waveforms do not scale uniformly with event magnitude. However, a linear scaling is likely to provide a reasonable approximation for the lower end of the magnitude range over which an event is detectable at a given station provided that the events are co-located and share a sufficiently similar source mechanism. For a given scaling factor, the probability of detecting a signal by any given method and sensor configuration will depend primarily upon the strength and form of the seismic noise at the time of arrival at the array.

The selected seismic array is ARCES in the far north of Norway and the selected signal is that generated by the tragic explosion onboard the Russian submarine Kursk on 2000 August 12, at a distance of approximately 475 km from ARCES. The ARCES array comprises instruments at 25 seismometer sites within an aperture of approximately 3 km; 3-component stations are located at four of the sites. The reason for selecting ARCES for this study is that its configuration provides excellent array-gain and slowness resolution for regional signals under traditional array processing and the high number of sensors means that there are many different subsets of instruments with which we can test the effectiveness of different array-configurations for detection by waveform correlation. The array design, with sensors arranged in concentric rings, provides natural subsets of sensors with different properties of noise and signal coherence. Although an essentially arbitrary choice of master signal, there are several reasons why the Kursk signal is a sensible choice. The location and nature of the source is well constrained and it provided a strong signal for each of the dominant regional phases; unlike many of the large mining explosions regularly recorded by ARCES, the source-time function was simple and so secondary phases are not obscured by complicated coda waveforms. Most importantly, all the channels of ARCES were operational at the time and so we can test the detection capability over the whole array. The signal is shown together with the location of this event in Fig. 4.

Figure 4

The explosion on the Russian submarine Kursk, 2000 August 12, as recorded by the ARCES array. The waveforms are filtered between 2.0 and 8.0 Hz and the Signal to Noise Ratio trace is calculated from the Pn-steered beam from all of the vertical component instruments. The f–k plots to the right indicate the slowness measured in windows of length 3.0 s starting at the stated times.

Instead of attempting to generate synthetic forms for noise, we superimpose the scaled-down signal upon many segments of actual ARCES data; the raw data were bandpass-filtered in the same frequency band chosen for the master signal prior to the combination of signal and noise. 2376 segments of data were selected by an arbitrarily defined formula that specified the starting point of each segment at regularly spaced intervals throughout 33 selected days. The only condition applied to this selection was that every channel had to contain real data; days on which one or more channel was subject to data-outages or systematic faults were not used. The presence of ‘signals’ in the chosen data segments was not checked for; as far as the detection of a given signal is concerned, an unrelated signal is merely an unfortunate occurrence of correlated noise.

All waveforms used in the study were filtered between 2.0 and 8.0 Hz; the SNR for the master signal is high for this frequency range and the signal is reasonably coherent over the full ARCES array. The first 60.0 s following the initial Pn-arrival (data segment starting 2000-225:07.31.46.25) were used as the master signal with filtering of the data performed prior to the cutting of the segment. An identical time window was cut for each channel, hence not taking account of the small time differences separating the Pn-arrivals at the different sites of ARCES. The interval chosen contains the same length of data that was used for the Kara Sea event example and, for the Kursk signal at ARCES, this includes the Pn, Pg and Sn phases but not the Lg phase and subsequent coda. The part of the wave train not included is that which would probably be susceptible to the greatest change given a co-located event of lower magnitude.

For the detection by the STA/LTA detector, an SNR trace was calculated for a beam steered for an optimal Pn-signal (Azimuth = 86.17°, Apparent Velocity = 8.29 km s ⁻¹: see Fig. 4). A detection was deemed to be made if an SNR of 3.2 was attained in the 8.0 s long interval [t_E− 3.0 : t_E+ 5.0], where t_E is the anticipated arrival time of the Pn phase on the beam. This excludes instances where the scaled signal was detectable but only for the later phase arrivals, but would include any instances where an unrelated signal in the data coincidently arrives at approximately the same time as our constructed signal. Such instances could, in principle, be identified by a measurement of slowness immediately following the detection but the evaluation of such measurements would require the definition of further sets of rules of acceptance or rejection of slowness measurements (typically made by the frequency–wavenumber or f–k analysis of the kind proposed by Capon 1969). This was deemed to be an unnecessary additional complication.

For the detection by waveform correlation, the single channel correlation traces were calculated for each sensor in the selected configuration. These channels were stacked as described in the preceding section and a detection defined by a scaled beam correlation coefficient greater than or equal to 6.0.

We first examine how the ability to detect the immersed signal changes with the imposed scaling factor for a single segment of noise. The top row of panels in Fig. 5 shows the Kursk signal superimposed without downscaling onto the noise at ARCES recorded at 03.15.00 UTC on 2005 May 2. No other signal was detected in the time window shown and the level of noise is comparable to that observed immediately prior to the arrival of the Kursk signal in 2000 August. As a result, the arrival on the Pn-beam is clearly visible and is detected with an SNR of 132. The panels to the right show that the signal is also detected clearly by waveform correlation on both a single channel and the full array; the scaled correlation coefficient exceeds 10.0 in both cases. In the second row of Fig. 5, the scale of the imposed signal is reduced by a factor of 20 and the SNR on the Pn beam is now only 6.7, a reduction by a factor of approximately 20. By contrast, the scaled correlation coefficients for both single channel and array are changed little. In the third row, with an imposed scaling of 0.025, the STA/LTA detector is essentially at the detection threshold in the frequency band selected with an SNR of 3.48 being recorded at the time of the Pn-arrival on the beam; the correlators still show clear detections for both single channel and full array.

Figure 5

ARCES signal from the Kursk event scaled down into an arbitrarily selected segment of array noise. For each row, the specified multiple of the master signal (60.0 s starting at time 2000-225:07.31.46.25) is immersed into ARCES data starting at the time 2005-122:03.15.00; this time corresponds to 0.0 on the time scales and the Pn-arrival at the central array element is at 1.0 s. The waveform shown is the Pn-steered beam and the panels to the right show the scaled correlation coefficients for the ARA0_sz instrument and full array (10.0 s displayed prior to and after the time of anticipated maximum correlation). Pn-beam waveforms are displayed on different axes in each of the panels whereas the same axes are used for all the scaled correlation coefficient traces.

With the Kursk signal reduced by a factor of 200.0 (Fig. 5, fifth row) the signal is now buried in the noise; although a clear detection on the full array correlation beam, the single channel correlation coefficient is considerably reduced, although still in excess of the threshold of 6.0. With an imposed scaling factor of 0.0025, the single channel correlation coefficient displays a local maximum at the time anticipated for the submerged signal, but, with a scaled coefficient below 5.0, it is such that it could probably not be detected by an algorithm with an acceptable false alarm rate. The detection on the array correlation beam is still above the set threshold. With the signal reduced by a factor of 1000.0, it is not detectable by any of the methods examined. This single example suggests that the array-based waveform correlation can detect a signal an order of magnitude weaker than the weakest detectable by a standard power detector.

The 2376 segments of noise chosen for the current study display different levels of background noise (including the possibility of high-SNR signals) and it is therefore by no means guaranteed that the situation indicated in Fig. 5 will apply more generally. Fig. 6 indicates the percentage of these submerged signals which were detected according to the criteria described above. The STA/LTA detector on the Pn-beam detects essentially all occurrences of the signal down to a scaling factor of approximately 0.07. With a scaling factor of 0.037, slightly over 90 per cent of the immersed signals are detected and, with a scaling factor of 0.016, fewer than 10 per cent of the signals are detected with the required signal-to-noise ratio. At the lower scaling factor of 0.010, essentially all of the signals are still detected using the correlators. The single-channel correlator detects a similar proportion of the signals as the STA/LTA detector at approximately 0.7 mag units lower; the full array correlator (with 25 channels) detects a similar proportion of events at approximately 0.4 mag units lower than the single-channel correlator.

Figure 6

Percentage of scaled signals submerged into ARCES noise detected by the methods indicated. The Kursk signal shown in Fig. 4 was superimposed with the scalings shown onto 2376 arbitrarily selected segments of data. ‘C ring’ refers to the seven short-period vertical channels at sites ARC1 up to ARC7, each at approximately 900 m distance from the central element. ‘sz’ refers to the short-period vertical channels of the array.

Fig. 6 shows the percentage of signals detected by waveform correlation for four different subsets of the ARCES sensors. Predictably, the best performance was achieved using the full array and the poorest correlator was the single-channel. Some single sensor correlators performed better than others, moving the detection curves slightly to the left when the rate of detection was improved; this is probably a function of the signal-to-noise ratio at each receiver site. The three-channel correlators provided by each of the four three-component sites at ARCES all performed slightly better in this experiment than the three-channel correlator provided by the short-period vertical channels at the ARA1, ARA2 and ARA3 sites. It may, however, be the case that three channels from distinct sites with a greater separation provide a more effective correlator than three channels from the same site due to the decorrelation of noise. The relative performance of three-component and vertical-only correlators is also likely to be a function of the partitioning of energy in the different directions. Further case studies will be required to test the performance of different configurations for different detection scenarios.

It is necessary to emphasize some limitations of the experiment described here. First, this represents the ideal case for the correlation detector; the signal detected is exactly the same shape as the template waveform. Scaling a seismic source by three orders of magnitude will not generally lead to negligible changes in the waveforms generated. Also, variations in the source-time function and exact source location will change the resulting waveforms to degrees that will vary according to the location and source type. Secondly, the study only used a single-template waveform with a single length and filtered in a single-frequency band. The length of usable signal will vary greatly from case to case as will the frequency content of the signal and the frequency content of the background noise. The performance of the STA/LTA detector will depend upon the latter two of these factors; the performance of the correlation detector will depend upon all three. Thirdly, it is not possible to conduct a systematic comparison of ‘false alarm rates’ for these different approaches since different definitions of false alarms apply. The STA/LTA detector is a general signal detector which will trigger whenever a sufficiently energetic data segment is encountered, whichever direction it comes from. Such a detection can only be considered to be a false alarm if it is subsequently identified incorrectly. What constitutes a false alarm for a correlation detector is dealt with in Section 4; the likelihood of a given template being triggered by an unrelated source will vary according to the complexity of the signal available, the monitoring network available and to what degree other signals may ‘resemble’ the template waveform.

We showed in Section 2 a method of inverting for the multiple of the waveform template which gives the best fit to the data (in this case comprising noise and an imposed multiple of the master signal). If we denote the imposed scaling factor by ε and the best-fit scaling factor by α then the ratio α : ε measures the quality of our amplitude estimate. The calculations from which Fig. 6 was generated provide a large data set of measurements of α and (the array correlation coefficient) which allow us to examine how the ratio α : ε varies as a function of . A ratio α:ε which differs greatly from unity indicates that our amplitude estimate is misleading and Fig. 7 indicates over which range of , for various array configurations, the amplitude estimate is meaningful.

Figure 7

On each occasion on which a multiple ε of the Kursk signal, immersed into background noise from ARCES, was detected using the master waveform as a template, the array correlation coefficient and the best-fit scaling factor α were calculated. The panels indicate the ratio α:ε as a function of for the ARA1_sz channel only (upper left), the ARC2 3-component station (upper right), the seven C-ring sz channels (lower left) and the full ARCES array (lower right). The plots were compiled by specifying a value for , calculating which 200 observations produced a measurement of closest to this value and then reordering the corresponding 200 values of α:ε. The central dashed line in each panel shows the mean, the solid lines indicate the 5 and 95 percentiles and the outermost dashed lines show the minima and maxima of the ratio α:ε. The large variability of the limits at low reflects the fact that fewer detections are made with low correlation coefficients and the data bins cover a larger range of .

For the single-channel case (ARA1_sz), the amplitude ratio between the data and master signal could be expected to be within 10 per cent of the actual value imposed for 90 per cent of calculations for which the array correlation coefficient was approximately 0.5. However, a correlation coefficient of nearly 0.8 was required before every scaling measurement was within 10 per cent of the imposed value. Adding additional channels reduces the bias in the amplitude estimate caused by noise; when the ARC2 three-component station is used, the ratio α:ε was within 10 per cent for all observations with as low as 0.53. In the 25-channel case (Fig. 7, lower right panel), 90 per cent of calculations returned inverted scaling factors within 10 per cent of the imposed scaling factor with a correlation coefficient as low as 0.2. A of over 0.3 was required before all the observations were within 10 per cent of the correct value.

4 The Automatic Screening of False Alarms

A ‘false alarm’ in the context of a site-specific correlation detector would correspond to the detection of a signal which correlates sufficiently well with the template waveform to exceed a predetermined threshold but which comes from an entirely unrelated source. A full-waveform cross-correlation detector is a highly sensitive form of detector which should, ideally, only result in a detection for a waveform which matches the template waveform extremely closely. However, identical matches do not occur in reality due to the presence of noise and it is exactly instances of low signal-to-noise ratio that we are attempting to detect; as was demonstrated in Fig. 3, many significant detections correspond to quite marginal values of the correlation coefficient. To prevent the non-detection of such cases, we need to set detection thresholds accordingly and we must therefore accept a certain risk of false alarms. In an on-line detector scenario, the number of false alarms should be reduced to an absolute minimum and as many of these as possible should be filtered out automatically.

Section 5 describes a case study in which we attempted to detect low-magnitude seismic events in the immediate vicinity of the site of a major rockburst using the waveforms from the main event as signal templates. A vast number of detections were made and visual inspection of waveform data at these times revealed that many of the more marginal correlations were entirely spurious and the result of completely unrelated signals. Such an example is shown in Fig. 8.

Figure 8

A typical correlation detector ‘false alarm’. A marginal correlation over a short time-window is sufficient to make the scaled array correlation beam exceed the detection threshold. The master signal resulted from a rockburst at the Barentsburg coal mine on Svalbard, approximately 50 km to the South West of the SPITS array. A phase arrival from a similar direction causes an increase in the correlation coefficients for each channel; however, unlike the example in Fig. 2, the maxima are not aligned and an f–k plot indicates a clear peak for a non-zero horizontal slowness. This is to say that, if these correlation traces corresponded to waveform data, this would indicate a plane-wave propagating across the array. All waveforms are filtered between 3.0 and 6.0 Hz.

The most obvious difference between this detection and those resulting from the Kara Sea event aftershock (Figs 2 and 3) is that, although the stacking of the correlation coefficient traces leads to a significant increase in the array correlation beam, the single-channel correlation traces are not aligned at the time of the maximum array coefficient. Performing broadband f–k analysis on the correlation coefficient channels at this point reveals that the delays associated with this misalignment are, in fact, indicative of a plane-wave propagation across the array with a non-zero horizontal slowness. Measuring the plane-wave propagation parameters for both master and detected waveforms, in a time-window for which the correlation is greatest, reveals that the slowness indicated by the correlation coefficient channels is approximately the difference between the slowness vectors for the two correlating incoming wave fronts.

That this should be the case is demonstrated in Fig. 9 where we define our master signal as a hypothetical plane wave front of synthetic wavelets and correlate against ‘data’ consisting of identical wavelets but delayed to simulate an incoming plane wave front from a different direction and with a different apparent velocity. An identical wavelet, w(t), is used for the arrival on each channel and, therefore, each of the cross-correlation traces achieves a value of 1.0. However, due to the disparity in interstation delay-times, this maximum is not reached simultaneously on all of the channels as would occur given an incoming wave front propagating with the same slowness (cf.Fig. 2). Nevertheless, a zero time-delay stacking of the traces does not lead to a simple cancellation of the signals and a detection is made.

Figure 9

A cross-correlation detection over a small aperture array can, in principle, occur when an incoming plane wave is cross-correlated with a plane wave with a different backazimuth (BAZ) and apparent velocity (v_app). The waveforms in the left and centre panels show synthetic signals simulating arrivals at the ARCES array with v_app= 4.5 km s ⁻¹, BAZ = 45° (left) and v_app= 7.2 km s⁻¹, BAZ = 135° (centre). Only three channels out of 25 used are shown. The slowness measurements from broadband f–k analysis are shown below: s^M (left) and s^D (centre). The correlation traces are shown in the upper-right panel. Broadband f–k analysis is performed on these traces in a short time-window centred around the maximum on the correlation sum; the slowness measurement, s^C, is indicated in the lower-right panel. s^C is demonstrably equal to s^D−s^M.

For the idealized situation shown in Fig. 9, we write v_i(t) for continuous data recorded at site i of an array and v (t) for the vector of array data, where [v (t)]_i=v_i(t). In particular, v^M (t) denotes the array data of the master signal and v^D (t) the data with which the correlation detection occurred. Given that the waveform recorded at each site is the same wavelet w(t) with an appropriate time shift, τ_i, then we can write

21

and

22

for the master and detected events, respectively. Both v^M(t) and v^D(t) are here defined to be perfect plane wave fronts associated with slowness vectors s^M and s^D respectively, such that, if the coordinates for site i are given by Δ_i= (x_i, y_i), then the site delay times are given by

23

and

24

where τ^M₀ and τ^D₀ are time shifts common to all sites within the array. The un-normalized correlation trace for channel i is given by

25

where T is the length of the extracted template and t_M is the starting time for the master signal [see eq. (5)]. From the definitions of our master and detected wave fronts [eqs (21) and (22)] we can write

26

which, due to the identical wavelets, will achieve an optimal correlation when

27

If we use t_i to denote the time t for which eq. (27) holds, and express the delay times for site i in terms of the array geometry [eqs (23) and (24)], it follows that the peak correlation occurs for time

28

where the slowness s^C and time shift τ_C are common to all sites of the array. Eq. (28) demonstrates that the times of maximum cross-correlation are associated with a slowness vector given by

29

as indicated by the idealized example in Fig. 9 and approximated by the real-data example in Fig. 8.

This result offers a very simple automatic method of eliminating obvious false alarms. Once the detection threshold on the correlation beam is exceeded, broadband f–k analysis is performed upon the single-trace correlation coefficient channels at this time. If the indicated slowness vector is very clearly not zero then it is almost guaranteed that the detected waveform came from a different source to the master waveform. Section 5 provides more details of how such a test was applied during the Barentsburg study.

This form of false alarm is only likely to occur when using small-aperture regional seismic arrays. Only in such circumstances is it likely that incoming wave fronts from unrelated sources are likely to correlate with a phase within the master signal at sufficiently many sites within a sufficiently short time-window to result in a detection on the array correlation beam. On the NORSAR array, where the interstation distances are too large for the effective processing of regional seismic phases, the time delays between any coincidental correlations on individual channels are typically so large that the beamforming of correlation channels effectively eliminates the chances of exceeding the detection threshold. Paradoxically then, while the trend in seismic monitoring using conventional array processing has been towards small-aperture regional arrays, a large-aperture array or sparse network of seismometers is likely to accommodate a far more robust correlation detector. If the only available observations of an event are on a small-aperture regional array, the danger of false alarms can also be mitigated by the extraction of longer time-windows and increasing the band width with which the waveforms are filtered.

5 Application: Monitoring of Seismicity at the Barentsburg Coal Mine on Spitsbergen

On 2004 July 26, a rockburst in the Barentsburg coal mine on the island of Spitsbergen caused the death of mineworker; a magnitude estimate of 1.76 was obtained for the event together with origin time 2004-208:06.42.19.44 and coordinates 77.9375 °N, 14.0703 °E. The event was well recorded by the SPITS array at a distance of approximately 50 km. Many roof-collapses have taken place at the mine (see Kremenetskaya 2001b) and Fig. 10 displays SPITS signals from the July 26 event together with the corresponding signals from three previous events. The waveforms for the four events filtered in the lower frequency band show great similarity and an event location based upon traditional array processing methods would be essentially unable to separate these events, even given the close proximity of the mine and array. The waveforms in the higher frequency band show far greater differences. By calculating a continuous array correlation coefficient, using the July 26 event signal as a template waveform, we aim to detect as many seismic disturbances from the immediate vicinity of the site of this event as possible.

Figure 10

Waveforms on the SPITS central seismometer (SPA0_sz) of four events at the Barentsburg coal mine; traces in the upper-left panel are filtered between 1.2 and 2.5 Hz and traces in the lower-left panel are filtered between 2.0 and 4.0 Hz. In each of the panels, the top trace is from the 2004 July 26 and the lower three traces are aligned according to the best correlation with a 5.0 s long time-window beginning at a time 2004-208:06.42.30. The map shows the location of the mine in relation to the SPITS array and the 3-component station at Kings Bay.

The application of the waveform correlation procedure serves the dual purposes of detecting events of a lower magnitude than would be otherwise possible and excluding events within a greater radius that a traditional array-processing algorithm (e.g. Gibbons 2005) could not resolve. The master waveform for this study consisted of 20.0 s of data beginning at 2004-208:06.42.29.695 for all the short-period vertical channels of SPITS, filtered between 3.0 and 6.0 Hz. Under the quarter-wavelength argument of Geller & Mueller (1980), the use of this frequency band should constrain two well-correlated events to be within approximately 150 m of each other (based upon an Lg velocity of 3.5 km s ⁻¹ and a frequency of 6 Hz). However, the ‘correlation distance’ of two different events is likely to vary greatly with local structure and heterogeneity (see e.g. Nakahara 2004) and measurements at local distances would be required to verify precise location estimates.

The correlation procedure detailed in Section 2 was run on all SPITS data from 2004 January 1 until 2004 August 10, the day on which SPITS was temporarily taken out of service for the installation of new instrumentation. During this period, a scaled array correlation coefficient, , exceeded a threshold of 5.5 on a total of 7292 occasions. However, many of these detections demonstrably did not correspond to an event from Barentsburg; an example of such an instance is shown in Fig. 8. An automatic screening algorithm was implemented by performing f–k analysis on the correlation coefficient channels and excluding any detections for which the slowness exceeded 0.04 s km⁻¹. Such detections are almost certain to correspond to coincidental correlations with wave fronts from different directions as demonstrated in the previous section. However, this condition alone is not sufficient to exclude all of the signals that were on inspection clearly not from Barentsburg.

Two additional conditions were empirically determined and applied as a requirement for a correlation detection to be considered a likely candidate for a Barentsburg event:

• (i)

  The relative power in the correlation beam had to exceed 0.39.

• (ii)
  The ratio had to exceed a value of 0.58. Here, t_M is the time of the maximum array correlation coefficient [defined in eq. (10)] and is defined by

  

  30

  where C^L_j(t_M) is the value of the closest local maximum to the time t_M on the correlation coefficient trace for channel j. This quantity measures the amplitude of the zero-delay correlation beam relative to the amplitude of the individual correlation coefficient traces at the time of detection; it is quite similar to the f–k power but is not constrained to the set of time delays permitted by plane wave front models. A low value of indicates that much energy is lost in forming the correlation beam.

These empirically determined conditions were devised simply from the manual observations of several hundred examples and examination of how these quantities varied when we had a clear false alarm or a very likely Barentsburg detection. A total of 1578 detections passed these tests. The Kola Regional Seismological Center (KRSC) in Apatity, Russia, operates two three-component stations very close to the mine. These stations provide data at local distances for some, although unfortunately not all, of the days in our 223-day period of investigation. The best of these stations, BRBB, was operational at the times of 50 of the detections and colleagues at KRSC have confirmed that, at each of the times provided, a signal is observed that is consistent with an event at the mine. A second station, BRBA, was operational at the times of 168 of the detections although only 78 of these detections could be confidently associated with a signal from a nearby source. It must be emphasized that the noise levels at the BRBA station are far worse than at the BRBB site and signals which were clearly visible at BRBB corresponded to events for which BRBA does not make a detection. It is admittedly unfortunate that so much of the on-site data are not available, but, for the days on which the local data exist, the indications are that these criteria were indeed able to remove spurious detections (Gibbons & Ringdal 2005a).

As in our scaled signal experiment, we calculated a best-fit scaling factor for each of the correlation detections. Assuming a magnitude m_b= 1.75 for the main event, we can estimate a magnitude for each of the events from this scaling factor. In order to assess the accuracy of these magnitude estimates, we made independent magnitude estimates by using the on-site stations installed by KRSC. It is assumed that the amplitude at these on-site instruments is dominated by the event in the mine and that the event magnitude is proportional to log (STA). Fig. 11 shows the magnitude estimates from the SPITS correlations plotted against the corresponding estimates from the BRBA and BRBB on-site instruments. The pattern appears largely consistent. The panels of Fig. 7 indicate the variability in the scaling factor which can be expected for different correlation coefficients even in the perfect situation in which detected and template waveforms are identical aside from added noise. In this case, there is no guarantee that the events come from exactly the same location and so fundamental differences in the waveforms may be expected which may lead to additional bias in the correlation magnitude estimates. The magnitude estimates from the on-site stations are subject to considerable environmental noise.

Figure 11

Magnitude estimates of events at the Barentsburg mine from the scaling factor inverted for from SPITS correlation detections and from short-term averages obtained from on-site stations BRBA and BRBB. The master event (not recorded by either BRBA or BRBB) was fixed to m_b= 1.75 for the purpose of scaling.

Finally, we need to compare the detection by waveform correlation with that achieved with conventional array processing. Of the 1578 presumed Barentsburg events detected by the correlation method, a total of 304 were detected using conventional processing on the SPITS array. Out of these, only 148 had both P-phase and S-phase detected with the appropriate time difference and consistent slowness estimates; this is the minimum requirement that a location estimate can be obtained. Of the remaining conventional detections, 63 only had a consistent P-phase detection and 93 only had a consistent S-phase detection. Many additional events detected by conventional processing were allocated automatic locations in the vicinity of the Barentsburg mine which did not correlate sufficiently well with the 2004 July 26 event to trigger a detection. It is conceivable that the seismicity in the entire region could be monitored using waveform correlation against a set of master event signals which was dynamically modified according to predetermined algorithms. However, the reliability of such a process would require verification using a far more sensitive on-site monitoring system than is currently available.

6 Conclusions

The principal results from this study are summarized as follows:

1. We have developed a waveform correlation detection algorithm that extends the traditional single-channel matched filter detector to a seismic array or network. This is done by applying array processing techniques to the individual correlation traces. We demonstrate that a simple stacking of the correlation traces can provide a significant beamforming gain relative to single-channel correlation detectors. The actual gain depends upon array geometry, the number of sensors, and the selection of time window and filter frequency band.

2. The basis for this performance is the fact that the correlation traces are coherent across a seismic array even when the individual signal traces are not. The reason for this is that if two events are co-located then the time separating the corresponding patterns in the resulting wave trains is identical for all receiver sites. We show as an illustrative example how the array-based correlation processing can be used to detect a small (m_b= 2.6) aftershock of the well known Kara Sea event of 1997 August 16 at the large aperture NORSAR array at a distance as large as 2300 km.

3. Frequency—wavenumber analysis of the correlation traces on a small aperture array provides an effective method of screening out a certain category of false alarms and can therefore be used to improve detector sensitivity by lowering the threshold for automatic array detection. We have demonstrated mathematically (under some simplifying assumptions) that the slowness indicated by the correlation traces is in fact equal to the difference between the slowness of the incoming signal and that of the master signal.

4. The array-based approach is equally applicable to local or regional seismic networks. Since signal coherency is not required, the method can be applied even at high frequencies across a sparse network.

5. An interesting aspect is the potential for control and verification of timing accuracy and data quality of individual stations in a seismic network. Given two or more co-located events, the cross-correlation traces between the events are expected to line up very accurately if the timing is consistent from one event to the other. Thus, data problems such as timing inconsistencies can be readily detected in such cases (see also e.g. Koch & Stammler 2003).

7 Discussion

In practical monitoring situations the use of the waveform correlation approach has some important limitations. These limitations apply both to the single-channel correlation and array-based correlation, and are directly related to the requirement for very waveform-specific master event templates. First, the master event and the target event must be very nearly co-located (a typical separation should be less than 1–2 wavelengths) in order for the matched filter detector to be effective. Secondly, the source-time functions of the master event and the target event must be very similar; the latter requirement means that large earthquakes should not be used as master events if the purpose is to detect small earthquakes. A low-SNR signal may provide a more suitable waveform template than a high-SNR signal if the high-SNR signal results from an event which is too large to have a comparable source mechanism. On the other hand, if the SNR is too low, the correlation coefficient may be dominated by the noise and a detection may not be made; also, magnitude estimates may subsequently be erroneous. The inability to obtain a suitable waveform template provides the greatest single hindrance to the practical implementation of these methods, although we believe that the examples in this paper demonstrate their potential for greatly improved effectiveness of automatic processing at least in situations such as earthquake aftershocks and rockburst sequences.

While the idealized study presented in Section 3 gives a good indication of the maximum improvement in detection capability for matched signal detectors for various configurations of an array over traditional array-processing, it does not address many fundamental questions. It does not address the variation of detectability with the length of the signal template or the frequency content of the master signal and the background noise; all these factors will vary from case to case in real-world monitoring scenarios. It also fails to address the level of deterministic variation permissible within a signal before it ceases to be detected by the master event waveform; such variation will be caused by differences in the source type and source location. For the full potential of the array-based correlation detector to be exploited in practical applications, there is a need for a station-site specific calibration encompassing the selection of the time window, threshold setting and filter frequency band.

The frequency–wavenumber analysis method described for identifying false-alarms when performing multichannel correlation detection on small aperture arrays becomes less effective when the direction of propagation for the incoming and master signal wave fronts becomes more similar. Further work will be required to establish criteria for identifying spurious correlation detections for low SNR signals.

The extent to which correlation detectors can be applied to the routine processing of data remains to be seen. Withers (1999) have already described an algorithm for signal detection and event location using waveform correlation. However, this system is aimed at providing a rapid and robust generalized detector and is based upon the waveform envelope; such a scheme would not be applicable to the detection of weak signals buried in noise as demonstrated in Fig. 3. On the other hand, the small geographical footprint of the waveform-specific detectors presented here restrict their use to far more specific seismic sources than envelope methods.

Much research is currently underway to determine the proportion of seismic events which are repeating sources (Schaff & Richards 2004a, b) and this may determine how widely correlators are employed for general detection purposes. The increase in computing power and storage media available to seismological institutes is rapidly making real-time cross-correlation detection highly feasible. However, what is certain is that correlation detectors and related methods (e.g. signal subspace detectors) will constitute the most sensitive detectors available for the monitoring of specific sites for which good calibration data are available.

Acknowledgments

This research was sponsored by the United States Army Space and Missile Defense Command and was monitored by AFTAC, Patrick Air Force Base, FL 32925, under contract number F08650-01-C-0055.

The authors are grateful to Vladimir Asming and colleagues at the Kola Regional Seismological Center in Apatity, Russia, for analysing the seismograms from the instruments installed at the Barentsburg mine. We also thank David B. Harris for useful discussions related to signal subspace detectors.

Maps were created using GMT software (Wessel & Smith 1995).

Copies of NORSAR technical reports can be obtained by sending an e-mail to or by writing to NORSAR, PO Box 53, 2027 Kjeller, Norway.