Typical methods
The sharp onset of high frequency P-waves on a seismogram can be picked by eye. The differential traveltime between two records is simply the difference between the two arrival times. Though simple in concept, this is a tried and true method which c
ommonly resolves traveltime difference of 0.025 seconds (a typical sampling period for teleseismic networks)
Long period S-wave signals (10-20 sec. periods) are more problematic. Their arrival is not as sharp due to their longer period. The standard way of determining differential traveltimes for these signals is by cross correlation. The lag which yields the highest correlation is considered the true time delay.
Recent interest in long period body waves
As broadband three-component instruments have come into wide use, S-waves have attracted increased attention. Unlike compressional waves, shear waves have a polarization. This polarization makes S-waves sensitive to features missed by compressional
waves. A velocity which varies with polarization can be indicative of anisotropy at the crystalline level. Anisotropy can be a manifestation of the stress regime or flow pattern. Anisotropy has proven to be a very powerful tool in determining the struc
ture of both the deep Earth and the crust. This project is concerned with methods used at the crustal scale.
The problem
Simple cross correlation works fine for determining traveltime delays between two distant stations. However, anisotropy studies require greater precision. "Extreme" anisotropy, on the crustal scale, results in delays on the order of one second. Typ
ically, delays are much smaller. Correlating long period waves to greater than 0.1 second precision is dangerous if the waves are not identical. Many factors can alter the shape of a wave during its passage through the crust. The most influential facto
r is the crustal reflection. A signal reaching the surface will have part of it's energy reflected. A portion of this energy can then be reflected back up to the surface. This results in a lower amplitude mirror image being superimposed on the signal at
some later time (there is a phase shift at the Moho reflection.)
figure 1
In a study of S waves, the first step is to isolate the shear wave packet. A time window must be chosen to insure that only the shear wave of interest is considered. Figure 1 shows two teleseismic signals from an event in the Andreanof Isl ands recorded near Lamont. Note how the maximum correlation is good for windows up to 35 seconds long. After this point, the correlation drops off rapidly. At wider windows we are only including more noise. A wider window only includes more noise.
This study makes extensive use of a routine named "correlate.m". It crosscorrelates two signals and determines the number of lags which provides the highest correlation. The number of lags translated into time using the specified sampli
ng rate. It also returns the value of the maximum normalized correlation. It is important to realize that this is not necessarily the zero-lag correlation.
figure 2
The frequency window can be chosen as well. A correctly applied lowpass filter can decrease ambient Earth noise called microseisms. Figure 2 shows the maximum correlation of the same two signals after aplying a variety of lowpass filters.
As frequencies from 0.5 down to 0.3 Hz are filtered out, the correlation r changes little. This implies that it is the lower frequencies which covary. As the filter is made more restrictive and progressively cuts out 0.3 to 0.1 Hz, the correlation incr
eases. This is precisely the microseism range. As they are eliminated, the signals covary more. This can be seen in figure 3. The characteristic signal in the top plot is still largely intact in the second. But as the filter becomes more restrictive
the quality of the signal degrades. Eventually the filter becomes so restrictive that it cuts out the signal itself and the correlation drops off. Beyond 0.05 Hz, the signal is gone. The high correlation of the long period roll is, in fact, a reason to
apply a highpass filter below 0.05 Hz. This is analogous to detrending a dataset.
figure 3
The "prelim.m" procedure selects two time series, crops them to the specified time window and normalizes them by subtracting the mean and divinding by the standard deviation.
The filtering routine "bpf.m" transforms the data to the fourier domain and applies a boxcar filter before returning to the time domain. This approach does not minimize leakage. The time series considered here have the bulk of their var iance in the middle of the signal so this problem is minimized. The advantage of this approach in this study is a clear definition of what frequencies are being passed by the filter.
To gauge the effect of this convolution on a seismogram, one can correlate the pure signal with its crustally altered version. A synthetic signal was used for this test. No single seismogram can adequately represent the variety of styles o
f S wave signals. While a noise-free signal may not be representative of all seismograms, it should do a better job of illustrating the processes which can affect all signals. A synthetic signal is used throughout this study, except in tests where we wi
sh to consider the effects of noise such as in the first section. Figure 4 is perhaps the most insightful diagram of this project. The top plot shows the signal before and after a typical crustal convolution. The middle is a plot of the time offset, ta
u, which results in the highest correlation, seen on the bottom plot. Convolution with the crust drags a signal backwards or forewards depending on the thickness of the crust and the wavelength of the signal. If the reflection arrives perfectly in-phase
or out-of-phase, only the amplitude is affected. At all other times, the wave is shifted as seen in the top plot. This points out tau's dependence on the ratio of the traveltime to the wave period.
figure 4
It is intuitive that the greater the reflection, the greater the time shift tau. This is substantiated in figure 5. A typical filter is shown with tau and r for different reflection strengths. This test was run for a 45 km crust. Figure
4 shows this to be a particularly vulnerable thickness.
figure 5
These tests reveal the sensitivity of tau to the crustal convolution for the test signal. The same test, performed on real seismograms shows similar trends though less variance. This may be due to the their wider frequency range. If a sig nal is spread between more frequencies, there are fewer which can be in-phase at any one time.
The second technique is to devise an inverse crustal filter to counter the effects of the crust. The third technique for comparing a signal altered in two different places is to convolve each signal with the opposite location's forward crus tal filter. The inverse filter required is addressed in this section.
It was expected that the inverse filter to the crustal effects would be highly unstable. Or that it would do a poor job deconvolving the signal. Such is the nature of inverse filters. For the simple crust model however, it worked surprisi ngly well. It was so predictable, based on the forward filter, that regression allowed an analytic expression to be developed.
The method first attempted for finding the inverse filter was in the frequency domain. The script invfilter.m divides the spectrum of the delta function by the spectrum of the filter. The result is, in theory, an inverse filter. Wh
en convolved with the original filter however, it did not produce anything resembling a delta function so it was thrown out. The least squares method was then used. This method in based on inverting an autocovariance matrix of the original filter. It i
s a laborious approach used in invfind.m because the autocovariance matrix for a 60 km crust filter can be 1475x1475 if the sampling rate is 40 Hz. Alternatively, it can be solved via Levinson recursion because of its symmetry. Figure 6 shows how
the convolution of the filter and inverse filter correlate with a delta function. As one might expect, the quality of the inverse filter decreases with larger secondary reflections. But even with a 0.7 reflection, the correlation is over 0.9.
figure 6
The shape of the least squares inverse filter is astonishingly simply though rather intuitive. It is a two spike filter of the same length. A positive spike of roughly unity occurs at the first position and a second spike occurs at the end . The major difference from the forward filter is that the second spike is positive. Conceptually this filter can be thought of as adding back the secondary signal which the crust filter subtracts out.
The magnitude of the two spikes is entirely independent on the crust thickness and depends only on the relative magnitude of the two spikes in the crust filter. The values of these two spikes can be seen in figure 7. Their smooth distribut
ion is easily fit with analytic expressions.
figure 7
A least squares non-weighted regression is performed in coefficients.m to create expressions for the origin spike and the secondary spike. They are functions of the ratio of the two spikes of the crust filter. Typically, when doing a regression one wants to limit the number of terms in the resulting expression. But in this situation we are not trying to decipher the underlying physics, we simply want an expression which will eliminate the need to invert the filter's auotcovariance matrix. Thus it is acceptable to use as many terms as we desire. From a practical standpoint however, we need only use enough terms to capture a significant portion of the variance. Regressions of order zero to four were performed. Each is listed belo w with its corresponding rms.
| order of fit | 0 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|---|
| rms error (origin spike) | 0.0409 | 0.0182 | 0.0032 | 0.0006 | 0.0002 |
| rms error (secondary spike) | 0.1275 | 0.0337 | 0.0038 | 0.0011 | 0.0005 |
The rms decreases by almost an order of magnitude when upgrading from a first to a second order fit. This is no surprise given the gentle curve of the two plots. A third order fit was chosen to be used in this study though it is probably o verkill. The resulting expressions are:
origin spike = 0.999 + 0.015R + 0.006R2 - 0.447R3
secondary spike = -0.002 + 1.058R - 0.325R2 - 0.519R3
Deconvolution
figure 8
Armed with an algorithm for the inverse crustal filter, the deconvolution approach can be tested. The tau-r plots in figure 8 demonstrate how well the deconvolved signal correlates with the unaffected signal for different reflection coefficien
ts. As in figure 6 this method is more effective for smaller amplitude reflections. Figure 9 assumes a constant reflection coefficient of 0.3 and evaluates the procedure for different crustal thicknesses. Recall that for a given reflectio
n coefficient, the amplitude of the inverse filter's spikes are fixed. However, the inverse filter is clearly more effective for thicker crust. In a continental setting, where the crust typically exceeds 30 km, the deconvolution method may be adequate.
figure 9
Cross convolution
The final method for correlating two signals affected by different crusts is to convolve each signal with the opposite location's forward crustal filter. The advantage of this approach is that it eliminates the use of an inverse filter approximation.
The disadvantage is that it further distorts two signals which have already been altered by the crust.
The shortcomings of this method are appearent in the the frequency domain. It will tend to amplify frequencies whose period is twice the filter length. Similarly it will diminish frequencies of the filter's length. Two signals with slight ly different frequency components will be altered in different ways eventhough they will have been subjected to the same filters. In test of seismograms, the cross convolution method still altered time delays by up to 0.2 s.
Comparison of Methods
A sample test of the two methods yielded the following results.
| Correlation of... | tau | r |
|---|---|---|
| raw signals | +0.1250 | 0.5812 |
| crustally altered signals (reflection coefficient = 0.3) | -0.0500 | 0.5421 |
| deconvolved signals | +0.1000 | 0.5844 |
| cross convolved signals | +0.1500 | 0.5568 |
For this example, ironically, the result is a tie. The initial passages through the crust shifts the best correlation by 0.175 seconds. Each method of remediation comes within 0.025 seconds of the original correlation though they err on op posite sides of the real value. There are two severe limitations to this test. First, this test assumes that the forward crustal filter is known exactly. Second, the crustal filter used here happens to have a high quality inverse filter. Recently, a method for the calculation of more detailed crustal filters has been developed. The filters could be substituted into these algorithms. By testing numerous seismic records, the remediation precedures could be tuned and even provide error bars. This wil l ultimately determine which procedure is superior. That is the next step.
Prepared for Quantitative Methods of Data Analysis course, 12/97