Rapid estimation of relative and absolute delay times across a network by adaptive stacking

Summary

1 Introduction

In a recent experiment with a heterogeneous set of short-period recorders we were seeking an effective means of providing a uniform and objective set of picks across an array of 72 stations. After a number of different trials we produced a simple scheme that exploits stacking based on a physical model of the propagation process (e.g. the traveltimes of a reference model), followed by determination of the best fit of each trace to the stack using a direct search. Iterative update then yields a direct estimate of the traveltime residuals from the original model predictions. At this stage we realized that we had developed a reincarnation of the method of adaptive stacking which was applied to array beamforming in the early days of digital seismology.

Jansson & Husebye (1966) and Gangi & Fairborn (1968) were among the first to develop and apply adaptive stacking techniques to data from seismic arrays. The relative delay times of the arriving signal across the stations of the array are given by the time-shifts (or steering delays) that optimally align the traces. In these early works, an initial stacked trace is obtained by roughly aligning all traces using the traveltimes from the Jeffreys–Bullen tables. After calculating the time-shifts that maximize the cross-correlation between each trace and the stacked trace, the traces are realigned and then restacked. This process is applied iteratively until the time-shifts converge, a process that can be described as adaptive stacking or adaptive beamforming. Bungum & Husebye (1971) describe possible sources of error associated with the adaptive stacking technique, and King et al. (1973) successfully apply the method to data from the Warramunga seismic array in Australia. Since this time, other methods which exploit waveform coherence have been developed that have achieved more common usage.

The aim of this paper is to show that adaptive stacking provides a fast and robust means of estimating relative and absolute delay times across a local or regional array, and is therefore well suited to dealing with large data sets of variable quality, such as those encountered in seismic tomography. Our modifications from the earlier work are the use of a more accurate reference model and a direct minimization of an L₃ measure of trace misfit to obtain the relative time-shifts required to achieve good trace alignment. Before describing our adaptive stacking approach in more detail, we briefly examine the characteristics of several other methods for exploiting waveform coherency across an array.

VanDecar & Crosson (1990) describe a method, often called multi-channel cross correlation (MCCC), for determining relative phase arrival times and uncertainty estimates for teleseismic events recorded by a local or regional network. For each pair of traces, a search is made to locate the peak of the cross-correlation function, which gives the relative time-shifts. The best fitting station delay times are then calculated by finding the least-squares solution of a system of over-determined linear equations formulated from the set of relative time-shifts. Timing uncertainties are based on the agreement between the least-squares station delay times and the relative time-shifts computed for each trace pair. Although this approach does not require an initial model, the user must provide a set of preliminary picks of the phase onset, which may be difficult in the presence of significant noise. Cycle skipping is also a potential problem when locating the cross-correlation maximum between trace pairs: a large outlier in the set of relative time-shifts can have a significant effect on the outcome of a least-squares minimization. To counter this problem, VanDecar & Crosson (1990) compute new relative time-shifts for trace pairs that produce very large equation residuals.

Least-squares inversion methods have also been used in beamforming. Mao & Gubbins (1995) develop a scheme that simultaneously inverts for time delays and stack weights, the latter of which reduces the effects of noisy traces on the resulting beam. The quantity that is minimized is the difference between each trace waveform and the corresponding beam formed by stacking the traces from all other stations. Noisy or incoherent data may result in cycle skipping problems and, like other beamforming methods, convergence may fail if the initial alignment contains many traces with large time-shifts.

Global optimization techniques can also be used to form array beams. Chevrot (2002) uses simulated annealing to minimize a cost function that measures the difference between each recorded trace and a reference waveform which can be time delayed. Despite a large number of unknowns equal to the number of samples in the reference trace plus the time delay value for each station, the method is apparently efficient enough to be used in routine analysis of seismic data. Absolute delay times can be found in addition to relative delay times by picking the onset of the waveform from the reference trace.

The underlying assumption in procedures which use cross-correlation between trace pairs to obtain time-shifts is that the signal shape is similar at each station. However, the adaptive stacking approach is more tolerant of significant waveform variability, since the reference trace will tend to average out these differences, and hence can work well with a heterogeneous station network. In the next section we describe our implementation of the adaptive stacking method, and then illustrate the technique by application to teleseismic data recorded at the TIGGER array in Tasmania, southeast Australia. In the presence of noisy data, it can be difficult to ascertain whether the apparent alignment of the traces is reliable. We therefore use events from nearby source regions, but with rather different waveform characteristics; the close correspondence of the residual patterns across the full array verifies the success of the method. Even when large time-shifts are introduced into noisy records, trace alignment is achieved with rapid convergence of the iterative procedure and excellent recovery of the imposed shifts.

2 Residual Estimation From Adaptive Stacking

We start by selecting a segment of each trace around the phase of interest and then achieve approximate alignment of the suite of traces by applying time-shifts, relative to a reference point, derived from a specific propagation model. For a local array it may be appropriate to adopt a simple plane wave front model, but for a network spread out over an extended region specific allowance needs to be made for the dependence on epicentral distance. We have found that the phase times for the ak135 model (Kennett et al. 1995) work well at far-regional to teleseismic distances, both for the selection of a suitable interval from each trace, typically 20–40 s long, and for approximate trace alignment. Any suitable phase can be selected, e.g. P, PP, pP, S, ScP, provided there is evidence of signal.

After the initial alignment using the set of time-shifts (moveout corrections) {t^c_i} derived from the specific propagation model, the linear and quadratic stacks of the records are calculated for a segment of the window that includes the desired portion of the traces associated with the selected phase. For N stations with records {u_i(t)}, the linear stack V_l(t) across the full suite of stations is defined by:

(1)

and the quadratic stack V_q(t) by:

(2)

The linear stack V_l(t) represents an estimate of the typical waveform across the array and the quadratic stack V_q(t) is an indicator of the spread in alignment between stations.

For each moveout-corrected trace u_i(t), the optimum match with the stacked trace V_l(t) is obtained by using a direct search over time-shift τ to minimize an L_p measure of misfit, defined by:

(3)

where M is the number of samples in the trace segment. τ is restricted to a specified interval (e.g. −1 < τ < 1 s) so that the computational effort of evaluating the misfit for each discrete time-shift is not large. The full search means that the minimum of the misfit measure for the ith trace in this interval is definitely located, with associated time-shift τ_i.

Once the time-shifts {τ_i} have been estimated for all the traces, the composite time corrections {t^c_i+τ_i} are applied to each trace to improve alignment. The linear and quadratic stacks are then recalculated with the revised time adjustments. The new stacked trace V′_l(t) represents an improvement on the initial stack. The alignment procedure is therefore repeated for each station trace, using the new stack, to produce an improved estimate of the residuals {τ_i} from values predicted by the propagation model {t^c_i}. The process is then iterated until accurate and stable trace alignment is achieved.

We have found that the choice of p= 3 in eq. (3) is particularly effective in achieving trace alignment for teleseismic data. The use of an L₃ measure may appear to be an unusual choice for measuring trace misfit with the stacked trace, due to its intolerance to outliers. Experiments with different choices for p suggest that the benefits of larger p come from the strong penalties imposed on even slight time offsets from the minimum of . This leads to more rapid convergence of the iterative process with very stable results achieved in two to three iterations.

The result of the adaptive procedure is to produce a set of residuals {τ_i} for each of the stations, relative to the predictions of the propagation model. When, for example, we use the ak135 predictions for teleseismic phases, we recover the pattern of residuals from ak135 across the array of stations from the shifts needed to achieve alignment of all the signals.

The high signal-to-noise for the stacked trace V_l(t), aids the picking of the onset of motion corresponding to the reference point and, hence, the absolute time for the particular phase at that point. The absolute times for the full set of traces can then be recovered by applying the composite time corrections {t^c_i+τ_i} for each of the stations.

When the waveform for a trace and the stacked signal are very similar (e.g. high signal-to-noise ratio) then the minimum at the optimum time-shift is pronounced; when they are dissimilar (e.g. low signal to noise ratio), the minimum is broader and more shallow. This behaviour suggests how one might measure the accuracy of traveltime residuals. Let represent the minimum misfit for a particular trace n. The uncertainty in the estimated residual is then defined as the smallest difference |T_n−τ_n| that results in is a calibration factor chosen a priori, and principally controls the absolute error level. For example, ε= 1.25 will result in an error estimate that equals the time-shift required to increase by 25 per cent from its minimum value. How one might choose an appropriate value for ε is examined more closely in the next section, but in general we have found that ε= 1.25 gives good results with the TIGGER data.

3 Examples Of Residual Estimation By Adaptive Stacking

We illustrate the application of our modified adaptive stacking procedure, using teleseismic data from the TIGGER experiment in southeastern Australia. TIGGER is a multifaceted seismic study of the crust and upper mantle beneath Tasmania and southern Victoria, and involved the deployment of 64 short-period and 17 broad-band recorders in 2001 and 2002. All short-period recorders and eight of the broad-band recorders formed a 72 station network with a nominal 15 km spacing in northern Tasmania (see Fig. 1). A feature of this set of stations is that it included a variety of instrumentation. There were L4A, L4C, Wilmore IIIA and IIIC vertical component short-period seismometers, Guralp 40T and 3ESP three-component broad-band seismometers, and three different recording systems: ANU short-period, Orion broad-band and Reftek broad-band. Although we attempt to unify the response function of the seismometers in pre-processing, some relative distortions of the waveform are unavoidable.

Figure 1

Location of the 72 stations that form part of the TIGGER array in northern Tasmania. Black dots denote short-period recorders and grey dots denote broad-band recorders. The surface expression of several geological features are also shown: Arthur lineament (AL), Mt. Read volcanics (MRV), Tamar fracture system (TFS).

Open in new tabDownload slide

Two events from nearby source regions (δΔ≈ 6°), but with different signal-to-noise ratios and contrasting P character, are chosen to test the method. The first example is an event in the Fiji Islands which occurred on 2002 April 20 (mb 6.0, depth 33 km). The P arrival comprises a relatively short sharp pulse with good signal-to-noise ratio on most (but not all) stations. The traces aligned using the ak135 arrival times for the Preliminary Determination of Epicentres (PDE) location are shown in Fig. 2(a): the imperfect alignment reflects the presence of lateral structure beneath the array. Note that a low pass filter has been applied to remove frequencies above 5.0 Hz.

Figure 2

Adaptive stacking example using an event from the Fiji Islands (located at 173.26°E, 16.38°S and 33 km depth): (a) initial trace alignment given by ak135; (b) final trace alignment achieved after five iterations; (c) pattern of residuals (in seconds) obtained from the difference between the initial and final stack. In both (a) and (b), the top two traces labelled zssl and zscp represent the linear and quadratic stack, respectively. Flat traces denote stations which did not record data [shown as crosses in (c)].

Open in new tabDownload slide

After five iterations of the adaptive stacking procedure, the traces are now very well aligned (Fig. 2b), and the time-shifts that needed to be applied at the individual stations form a coherent pattern of anomalies across the array (Fig. 2c). Although the adaptive scheme only requires two iterations to converge, the rapid computation time of about 0.25 s per iteration across all the 72 stations (on a SunBlade 150) means that a few extra iterations can be routinely applied to ensure the stability of the solution. It is interesting to note that even the very noisy traces (e.g. stations tsn1 and tsj4) appear to be correctly aligned. The onset of the stacked trace at the final iteration (Fig. 2b) is easy to pick due to the enhanced signal-to-noise ratio. A comparison of the quadratic stack at the initial (Fig. 2a) and final (Fig. 2b) stages, clearly indicates the improved alignment.

The second example consists of data from an event in the Vanuatu Islands on 2002 April 11 (mb 5.9, depth 10 km). Compared to the Fiji event shown in Fig. 2, the waveform for P is much longer and less distinct for this shallower event because of interference with the depth phases (pP, sP), and the signal-to-noise ratio is not as high. Fig. 3(a) shows the traces across the TIGGER network approximately aligned using the ak135 arrival times for the PDE location. As in the Fiji example, five iterations of the adaptive stacking procedure result in a set of exceptionally well-aligned traces (Fig. 3b), but convergence is achieved after only two iterations. The first motion of the optimally aligned stacked trace can be easily picked, so accurate relative and absolute delay times can be obtained.

Figure 3

Adaptive stacking example using an event from the Vanuatu Islands (located at 167.69°E, 14.39°S and 10 km depth). Refer to the Fig. 2 caption for a description of each plot.

Open in new tabDownload slide

The pattern of traveltime residuals across the array for the Vanuatu Islands event (Fig. 3c) is not only coherent, but compares remarkably well with the residual pattern for the Fiji event in Fig. 2(c). Given that the two events are from a similar location, but have different waveforms and signal-to-noise ratios, this result is a strong endorsement of the success of our adaptive stacking procedure in extracting the signal of lateral heterogeneity. Another desirable feature of the method is that it is relatively insensitive to high frequency noise. In both Figs 2 and 3, all traces are low pass filtered at 5 Hz: using smaller values (e.g. 2 Hz) removes more of the high frequency noise, but we found that the residual patterns were virtually unchanged.

Using the approach described in Section 2, the uncertainty associated with each traveltime residual estimate is calculated for the two events (see Fig. 4) with ε= 1.25. A minimum value of 37.5 ms (equal to 75 per cent of the sample interval) is also imposed on the uncertainty estimates: it is unrealistic to expect greater accuracy given the presence of noise and imperfect coherence of the waveform across the array. As noted in the previous section, relative errors are not strongly affected by variations in ε, which is principally a function of absolute error. However, a knowledge of the absolute error is important in many applications, such as seismic tomography, where one needs to know how well a given model satisfies the data. This means that the value given to ε should not be chosen arbitrarily.

Figure 4

Uncertainty estimates for the traveltime residuals obtained by adaptive stacking: (a) Fiji Islands event (see Fig. 2); (b) Vanuatu Islands event (see Fig. 3). For both plots, stations that did not record data have error estimates set to zero.

Open in new tabDownload slide

The approach we use for calibrating the uncertainty estimates shown in Fig. 4 is based on the ability of adaptive stacking to recover synthetic time-shifts applied to all traces. For both events, random time-shifts with a Gaussian distribution and standard deviation of σ= 0.75 s are applied to the aligned data (as shown in Figs 2b and 3b). Adaptive stacking is then used to try and recover the imposed time-shifts: the RMS difference (δ) between the applied and recovered time-shifts is used as an indicator of the absolute error. For the Fiji Islands event, δ= 60 ms; for the Vanuatu Islands event, δ= 73 ms. This difference in δ is consistent with the lower signal-to-noise ratio of data from the Vanuatu Islands event. The appropriate choice for ε is the one that produces error estimates with an RMS value equal to δ. We found that ε= 1.25 produces error estimates for the Fiji Islands and Vanuatu Islands events with RMS values of 61 ms and 73 ms, respectively, which are almost identical to the δ values obtained using the imposed time-shift technique. It should be noted, however, that the value of δ has some dependence on the size and distribution of the imposed time-shifts, so although the calibration scheme appears to work well, it is not completely free from ad hoc assumptions.

A comparison of the error patterns displayed in Fig. 4, with the seismic traces in Figs 2 and 3 shows a satisfyingly strong anticorrelation between the strength of the signal and the size of the uncertainty. For example, stations tsq4 in Fig. 2 and tsp2 in Fig. 3 have little evidence of P and large error estimates result. As previously noted, high frequency noise may give a visual mask to the true signal, but not significantly influence residual estimation. For instance, the error associated with station tso2 is larger than that of tsc4 in Fig. 4(a), even though the latter clearly has a much smaller signal-to-noise ratio (Fig. 2). The reason for this is that the waveform recorded by tso2 has a reduced coherence, due to the presence of noise within the frequency band of the signal, structural effects or a change in instrument response. On the other hand, the signal recorded by station tsc4 is more coherent with the reference stack trace when separated from the high frequency noise that contaminates it.

The magnitude of the ak135 residuals in Figs 2(c) and 3(c) do not exceed 0.5 s, so the initial stacked traces are quite similar to the final stacked traces. In areas of greater upper mantle heterogeneity, the magnitude of residuals will be larger; for example, Seber (1996) observe residuals as large as 1.5 s due to lateral structure beneath the Atlas and Rif mountains in Morocco. To examine the behaviour of adaptive stacking in such circumstances, we generate random time-shifts with a Gaussian distribution and standard deviation of σ= 0.75 s and apply them to very noisy, unaligned data from an event that occurred in the south Vanuatu Islands on 2002 April 10 (mb 5.4, depth 33 km). Inspection of Fig. 5(a) shows that the signal-to-noise ratio of these data are much poorer than those used in Figs 2 and 3, and that time-shifts well in excess of 1 s are required to align many of the traces (where signal is evident). This example clearly represents a significant challenge to any method that relies on waveform coherence to obtain relative arrival times.

Figure 5

Adaptive stacking example using an event from the south Vanuatu Islands (located at 169.29°E, 20.74°S and 33 km depth) with Gaussian noise (σ= 0.75 s) used to offset the initial trace alignment: (a) initial trace alignment; (b) final trace alignment achieved after five iterations; (c) pattern of residuals (in seconds) obtained from the difference between the initial and final stack with imposed time-shifts removed.

Open in new tabDownload slide

After five iterations of the adaptive stacking scheme, most traces appear to be well aligned (Fig. 5b): the change in spread of the quadratic stack between the initial and final stacks supports this claim. However, many of the traces are very noisy so it is difficult to ascertain whether they are all correctly aligned or not. Fortunately, the event we are using is within 7° of the two events we have already examined, so we can compare the ak135 residuals obtained via the adaptive stacking process. For the south Vanuatu Islands event, we remove the imposed time-shifts from the recovered residuals in order to make a comparison. Fig. 5(c) shows the ak135 residuals for the south Vanuatu Islands event. The residual pattern closely resembles those of Figs 2(c) and 3(c), with only a few exceptions. For example, the residuals at stations tsp5 and tso3 differ greatly from those obtained for the Fiji Islands and Vanuatu Islands event. Inspection of these two traces (Fig. 5a) reveals little evidence of signal, and error estimates obtained with ε= 1.25 are of the order of 250 ms. Therefore, if estimated uncertainties are accounted for, the residual pattern is consistent with those obtained for the Fiji Islands and Vanuatu Islands events. As expected, the RMS value of the error estimates (83 ms) is larger than those determined for the events shown in Fig. 2 (61 ms) and Fig. 3 (73 ms).

4 Discussion and Conclusions

The strength of adaptive stacking is that it is simple, fast and works well in the presence of noisy data and some waveform distortion, as illustrated by the teleseismic examples from the TIGGER experiment. The method depends on having an a priori propagation model, but, as we have seen, in teleseismic applications the use of a 1-D Earth model such as ak135 is sufficient for the initial trace alignment. The time-shifts determined from the adaptive stacking procedure then provide direct estimates of the residuals from the ak135 model. The test case displayed in Fig. 5 demonstrates that we can expect adaptive stacking to work with noisy data and in areas of significant heterogeneity that result in less accurate initial trace alignments than were obtained for the TIGGER data. The computational cost of adaptive stacking scales linearly with the number of traces, which means that it can be efficiently applied to data from large seismic arrays.

Compared to the MCCC technique of VanDecar & Crosson (1990), adaptive stacking is expected to be faster, since there is no need to solve a system of linear equations or cross-correlate all trace pairs, both of which scale approximately quadratically with the number of traces. MCCC also requires an initial pick of the arriving waveform, which is not required with adaptive stacking. We would expect that adaptive stacking is more robust than MCCC, and we have demonstrated that our scheme can deal with very noisy traces and with a variety of seismometer responses. The method of Chevrot (2002), based on a simulated annealing inversion, may be more robust as no a priori model or preliminary picks are required and a fully non-linear waveform fitting approach is used to determine the delay time residuals. However, it is not clear how it would perform in the presence of high frequency noise and varying signal character. Adaptive stacking is much simpler to implement than the Chevrot (2002) procedure, and again uses non-linear waveform matching.

Our examples of adaptive stacking have used teleseismic waveforms, but the technique is not limited to this class of data alone. The philosophy of adaptive stacking is to exploit any available models of Earth structure in order to detect unexplained traveltime variations in observed data. This can be done at any scale provided the waveform coherency assumption is valid and reasonable a priori models are available. For example, across a medium aperture array a simple plane wave front model is normally adequate.

Adaptive stacking, first developed several decades ago by Gangi & Fairborn (1968) and others, is a valuable alternative to other methods such as MCCC that are prevalent today (e.g. VanDecar 1995; Frederiksen 1998; Ritsema 1998; James 2001; Tilmann 2001). As we have shown, the introduction of ak135 as a traveltime predictor, and the use of an L₃ norm to measure trace misfit, rather than cross-correlation, results in a robust technique for teleseismic applications.

Our adaptive stacking software is freely available from: http://rses.anu.edu.au/seismology/astack as FORTRAN code, along with instructions on its use.

References