Chapter 5: Coherence of Seismic Body Waves
1. Introduction
The most common technique used to study details of earthquake sources involves deploying a network of seismic stations around a known region of active seismicity. This is usually done by emplacing portable recorders in the aftershock area of a large earthquake or by setting up a permanent array that encompasses an active fault zone. In both types of experiments the station separations range from one to several tens of kilometers. Even at this close spacing, this corresponds to many spatial wavelengths for those frequencies above a few hertz, making it hard to resolve the desired details.
Several studies have examined the variability of source parameter estimates determined from bodywave spectra calculated at different stations [e.g., Fletcher et al. 1984, Cranswick et al. 1985]. This variability can be caused by the seismic radiation originating preferentially from different parts of the focal sphere, by the propagation of seismic waves through a complicated geological structure, or by localized site and topographic effects. In recent years there has also been considerable interest in the variation of seismic waveforms measured at inter-station spacings of tens to hundreds of meters. The response of shallow sediment-filled valleys was carefully measured by King and Tucker [1984] and Tucker and King [1984], who placed an array with this range of station spacings across the Chusal and Yasman Valleys in the Soviet Union. Tucker et al. [1984] examined the differences between stations located on hard rock which included ridge, slope, and tunnel sit
In this paper, results are presented from a study of the variability of seismic spectra and the spatial coherence of the seismic wavefield over horizontal distances of 50 to 500 meters. These data were recorded on a nine-station seismic array operated on a fairly level granitic terrain for eight months in 1985. The criteria for choosing the array site is discussed in section 2. The experimental design and the events recorded by this array are summarized in sections 3 and 4 respectively. The data analysis procedures, which include multitaper spectral and coherence estimates, are presented in section 5. Examples of analysis results are shown in section 6. Section 7 examines the stability of earthquake source parameters calculated from body wave spectra. Section 8 contains a more detailed discussion, and section 9 lists the conclusions. Appendix A contain
2. Site Choice
Site selection for this seismometer array was based on several criteria. First, the array was to be located in an area where the seismic signals are not expected to be distorted by geological structures near the array. Second, the region of the site must have little topographic variation which can possibly amplify or attenuate seismic signals [Geli et al., 1988]. Finally, the site should be in an active seismic region so enough events can be recorded.
Piñon Flat, which is located in the northernmost part of the Peninsular Ranges batholith in southern California, is a site that satisfies these conditions. The site is situated on the northern part of the Santa Rosa Mountain Pluton (also known as Haystack Pluton) [Dibblee, 1981]. This pluton has a granodiorite composition [Parcel, 1981] and is approximately 30 km long and 15 km wide in a NW-SE orientation.
The array was installed at the Piñon Flat Observatory (PFO) at an elevation of 1280 meters above sea level. Figure 1 [from Wyatt, 1982] shows a map of the area with the array located inside the heavy black ring. Pinyon Flat is bounded on the north by Asbestos Mountain (elev. 1585 m) and on the south by the Santa Rosa Mountains. Two canyons, Palm Canyon and Deep Canyon, bound Pinyon Flat on the west and east respectively. The general array area has little topographic relief within several kilometers of the array. The metasedimentary and metamorphic rocks found in the area of these canyons are roof pendants and screens formed during the emplacement of the plutons.
The near surface geology of the PFO site has been studied in some detail. Wyatt [1982] found the top one meter of ground to be almost totally decomposed granodiorite. This grades into grus, then corestones, becoming jointed granodiorite at a depth of 25 meters or less. In the same study, results of a shallow seismic refraction survey at PFO were presented. This survey showed that the compressional velocities were confined by the extremal bounds from 2.0 km s-1 to 3.3 km s-1 at 18 m depth. P-wave logging of a 300-meter borehole was conducted by Fletcher et al. [1989]. They determined compressional velocities of 4 km s-1 at 65 m, rising to over 5 km s-1 by 80 m depth. These measurements suggest that this site has an extremely thin weathered layer overlying a rigid basement.
The seismicity in the region is quite high since the array is located about 12 km northeast of the San Jacinto fault zone and about 25 km southwest of the San Andreas fault zone [Wyatt, 1982]. The seismicity of magnitude ML ≥ 3 events in the years 1980-1984 is shown in Figure 2, along with the major fault traces.
3. Experimental Design
For this experiment, it is important to ensure that the transfer functions of the sensors and the recording instruments at the sites were very nearly identical so that the differences between the recorded signals are caused solely by the differences in ground motion at each recording site. All individual recording sites were similar in construction and in equipment set up.
The array was made up of nine stations. At each site was a triaxial set of Mark Products L-22 seismometers which have a natural frequency of 2 Hz. The output of the seismometers was recorded by a Sprengnether DR-100 digital cassette recorder which sampled each channel at 200 samples/sec with a 12-bit analog-to-digital converter. The anti-aliasing filters had a 5-pole Butterworth response with a 50 Hz corner frequency.
The array was set up as a set of three nested equilateral triangles (Figure 3) in February 1985. The center point of the array was located on top of the southeast vault of the northwest-southeast laser strainmeter at PFO. Each sensor was precisely located using a Nikon NTD-2 infrared electronic distance measuring device. One problem with using a topographically flat area for an array is finding instrument locations with equipment requiring line-of-sight; local obstructions such as scattered Piñon pines can make this quite difficult. The first station sited was AZG because it was near the laser strainmeter and it was certain there would be no obstruction over a 300 m distance. The other 300 m radius stations (AZH and AZI) were located at clockwise angles of 120° and 240°, respectively, from the reference line of the centerpoint to AZG. The middle triangle of stations AZD, AZE, and AZF were located at a distance of 100 m from the center point and 60&de
Each station consisted of a battery-powered digital cassette recorder and a triaxial seismometer which was shallowly buried to improve coupling to the ground. The digital recorders were modified to run off a single external 1 MHz clock to synchronize the digitizers and eliminate any relative timing drifts between individual recorders. Each digitizer operated with its own independent trigger. The equipment was serviced at 10 to 14 day intervals. After the data were recorded in the field the data cassettes were returned to the USGS labs at Menlo Park and played back onto a computer. Each event from each digitizer was plotted and checked for recording problems or false triggers. Events that were well recorded on five or more stations were chosen for detailed analysis.
4. The Data Set
During the operation of this array from late January to September 1985, there were eight events that triggered five or more recorders simultaneously. These events had local magnitudes ML between 2.6 and 3.9, as estimated from the southern California seismic network array data [Norris, 1986; Given, 1987]. The events, marked with asterisks in Figure 2, have epicentral distances ranging from 14 to 42 km and depths between 4 and 16 km. As none of these events were in the same place, there were seven independent azimuths for this set of earthquakes. Table 1 shows various source parameters including locations, epicentral distance and azimuth from PFO, local magnitude, seismic moment, and stress drop. The moment and stress drop estimates were calculated from data collected on the Anza seismic array using the method described by Fletcher et al. [1987].
An example of a seismic time series from this experiment is shown in Figure 4. The figure displays the transverse components of event 8, located 17 km southwest of the array. Note that the visual correlation between each series is apparently quite low for all parts of the seismogram except for the large low frequency pulse at the beginning of the shear wave arrival. In general, the seismograms for the three components for a given event show some similarity to each other during the P and S phases while the coda for each phase appears to be random in character.
5. Data Analysis Techniques
Since we are dealing with short time series and have all the attendant problems of spectral leakage and variability for the estimates, we decided to use the multiple taper technique devised by Thomson [1982]. The application of this technique for calculation of spectra from seismic data has been described in Chapter IV. These equations in a slightly modified form are used for the calculation of the spectra described here. A deviation of these modified equations from Vernon et al. (manuscript in preparation) is given in the Appendix, with the pertinent ones also reproduced in the text. In practice, spectra calculated by both methods from the same data are similar in appearance.
The multitaper power spectral method attempts to minimize the bias due to spectral leakage outside a chosen central bandwidth. The implicit assumptions for this technique are that the data are a Gaussian distributed, evenly sampled, unaliased, ergodic, stationary time series. These tapers {vt(k)(N,W)}k=1K, are functions of the time series length N and the specified "inner" bandwidth W. The tapers are the set of eigenvectors which are solutions of a matrix eigenvalue problem [Park et al., 1987; Chapter IV] and have associated eigenvalues 1 > λk (N,W) > 0. The data x t are multiplied by each of the K tapers in turn, and the discrete Fourier transform is calculated, as presented in Chapter IV, eqn. (4.8),
y k (f) = sum from t=0 to N-1 ~ v sub t sup (k) ^(N,^W) ^x sub t ^ e sup {-i2 pi ft} ~.
The adaptive spectral estimate S(f) is a weighted average of the K y k (f)2 s. The weights dk (f) depend on frequency, and on the taper order k, since the bias resistance for each taper decreases as k increases.
S(f) = A over K ~ { sum from k=0 to K-1 ~ λ k ^ d k 2 (f) |y k (f) | 2} over {sum from k=0 to K-1 ^ d k 2 (f)}
The coefficient
A = sum from k=0 to K-1 ~ λ k-1
and the weights dk (f) are calculated using the recursion relation in the Appendix.
As pointed out by Thomson [1982], the orthogonal properties of the multitaper spectral analysis can be extended to calculate coherences. If there are two series xti and xtj which have equally sampled data, then (5.2) can be generalized to
S ij (f) = A over K ~ {sum from k=0 to K-1 ~ λ k ^ d k sup i (f) (y k sup i (f)) sup star ^d k j (f) ^ (y k j (f)) } over {left [ ^ sum from k=0 to K-1 ~ d k sup i (f) 2 ^ right ] sup 1/2 ^ left [ ^ sum from k=0 to K-1 ~ d k j (f) 2 ^ right ] sup 1/2 }
If i = j, then (5.3) is the autospectrum. If i ≠ j, then Sij(f) is the cross spectra between xti and xtj, where the weights d kj(f) are the same ones used for the Sjj(f) . The magnitude squared coherence γij(f) is then found directly from
γ ij 2 (f) = {|S ij (f)| 2} over {S ii (f) S jj (f)}
and the phase φ (f) is just the phase of Sij(f) . Multitaper coherence estimation has been used by Lanzerotti et al. [1986] where they choose dkj (f) ~ 1. Here we have a larger dynamic range, so it is necessary to use the dkj (f) explicitly.
One property of the multitaper technique is the orthogonality of each tapered estimate in both the frequency and time domains. The independence of the yk(f) implies that estimates of the mean, bias, and variance for the power spectra and the magnitude squared coherence can be made. Unfortunately the underlying probability distribution is not always known for S(f) and γ2(f). Furthermore, there are a limited number of samples for the estimation of S(f) and γ2(f). Under these constraining conditions, an unknown probability distribution and the small number of samples, it is necessary to use an empirical method of estimating bias and variance. Lanzerotti et al. [1986] first introduced the use of jackknife statistics with multitaper estimates to evaluate the variance of the pseudo-autocorrelations of magnetic field data and cable voltages. Lindberg and Park [1987] also used a jackknifed multitaper estimate for the errors in peak frequency measurements of the earth's normal modes.
The development and background of the jackknife are presented by Miller [1974] and Efron [1982], with an outline of the pertinent equations following in this text. The basic jackknife starts with the assumption that there are n independent and identically distributed random variables Zl (1 ≤ l ≤ n) from an unknown probability distribution. For each random variable Zl there corresponds an observation of that variable, Zl. The quantity of interest is estimated by the statistic &θ; = &θ; (z 1, z2, ..., z n). The jackknife makes use of a set of deleted estimates θ (i) where the i th observation zi is not used in the estimate:
θ hihat sub (i) = θ hihat (z sub 1 ,^ z sub 2 ,^ ...,^ z i-1 ,^ z i+1 ,^ ...,^ z sub n )
The jackknifed estimate of the mean of $θ hihat$ is θ hihat sub {( cdot )} = 1 over n ~ sum from i=1 to n ~ θ hihat sub (i) ~.
Tukey [1958] showed that an estimate of the variance of θ is
var ^ "{" θ hihat "}" = n-1 over n ~ sum from i=1 to n ~ ( θ hihat sub (i) ^-^ θ hihat sub {( cdot )} ) 2 ~.
Miller [1974] points out that it is beneficial, if possible, to transform the statistic θ so that its underlying probability distribution becomes more Gaussian. This can be important when estimating standard errors on quantities such as variance, which has a range of [0,∞ ), or magnitude squared coherence, which has a range of [0,1]. The transformations used to convert to a more normal distribution for power spectra are θ (f) = log (S (f )) [Miller, 1968]; and for magnitude squared coherence θ (f) = arc tanh ( sqrt {γ2(f)}) [Jenkins and Watts, 1968].
To create the k deleted estimates $θi with i = 1,..., K for the power spectra estimation, (5.2) is modified to read
S i (f) = {A-{λ i sup -1}} over {K-1} ~ { sum from {k=0,k ≠ i} to K-1 ~ λ k ^ d k 2 (f) |y k (f) | 2} over {sum from {k=0,k ≠ i} to K-1 ^ d k 2 (f)}
Equations (5.3) and (5.4) are modified in a similar manner to determine the deleted estimate of magnitude squared coherence.
6. Data Analysis
The analysis of the body wave phases for each event starts with the rotation of the time series into radial and tangential components followed by the calculation of the power spectrum of each component. For this data set a time-bandwidth product (NW) of 6.5 with 12 tapers was used in each computation. This gives at most 24 degrees of freedom for S(f) at each f. The time series have been aligned to the P wave arrivals. Each spectrum used two seconds of the time series centered on the particular phase of interest, as shown in Figure 4. Next, the spectra from each station were combined to estimate a average of the log spectrum, SA(f), using the following equation:
S A (f) = exp ^ left [ ^ 1 over N ~ sum from i=1 to N ~ ln ^ (S i (f) ) right ]
where Si (f) is the spectrum at station i and n is the total number of stations. The 95% confidence limits C ± (f) are found using the following expression:
C ± (f) = α sup {+- ^1} (f) S A (f)
where α sup {+-^1} (f) = exp ^ [ 2 ~ sqrt {VAR ( ln (S i (f)))} ^ ]
and VAR (β) is the variance of β.
The gross characteristics of the averaged log spectra are all quite similar (Figure 5). Since the original time series are recordings of ground velocity, the spectra decrease in power as frequency approaches 0 Hz. The velocity power peak for all phases are found to be between 5 and 22 Hz. A remarkable feature of the averaged log spectra is their variance. The 5% to 95% confidence limits α ± (f), determined from the variance of the log spectra between stations, has a mean value of approximately 4, but can vary from 1.4 to over 13.0 as a function of frequency. For reference, the ratio of the 5 to 95% quantities of a χ242 distribution is 2.83. The averaged log spectrum SA(f) from the P wave recorded on the vertical component for event 1 is shown in Figure 5a. On closer examination of the individual spectra, they can be divided into low and high frequency parts by their variance. The individual spectra used to form the averaged log spectra in Figure 5a are shown in Figure 5b. From 0 to 15 Hz all the spectra are similar, differing by a scale factor. The detailed structure of the spectra is not correlated in the band above 15 Hz for these P waves. The S waves show a similar phenomenon except that decorrelation starts at approximately 10 Hz. This feature cannot be explained by the error estimates of the individual spectra. The average double-sided 95% confidence limit for each individual spectrum is shown on the plot, and is much smaller than the differences between spectra. Clearly the details of the seismic body wave spectra above 15 Hz for P waves and above 10 Hz for S waves are controlled by variable localized site effects.
The second part of the analysis was an examination of the coherences between each station pair for every event. The coherences were calculated using data aligned to maximize the cross-correlation between waveforms, and using the same spectral parameters (12 tapers, time-bandwidth product NW = 6.5) as in the power spectra estimation. Figure 6 displays an example of the magnitude squared coherence (MSC) estimates from the transverse component of the shear wave time series in Figure 4.
The MSC estimates shown in Figure 6 gives the general characteristics of all the MSC estimates which were calculated. Each estimate has large random variations superimposed on a general trend of high MSC at low frequency decaying towards zero as the frequency increases. The highest MSC estimates averaged over frequency are found for the P waves measured on the vertical components where MSC values of greater than 0.7 can be found out to over 40 Hz. The detailed structure of the MSC estimate does not correlate from station to station on the same event, nor from event to event for fixed station pairs. The variability of these estimates makes it difficult to characterize expected MSC properties using only single sets of station pairs or data from just one event.
The process of finding a general set of properties using all the MSC estimates for each component of the P and S waves (115 and 140 estimates respectively), started by examining the signal to noise properties of the power spectra. Wyatt and Berger [1980] investigated the coherence between the outputs of two adjacent tiltmeters. They assumed that the output of each instrument xti consisted of a common mode signal ut with an added noise signal vti. If the power signal-to-noise ratio is defined
σ 2 (f) = U(f) / V(f)
where U(f) and V(f) are the power spectra of the common signal and the independent noise, respectively, then the magnitude squared coherence can be estimated by
γ if 2 (f) = 1 over {(1^+^ 1 over {σ i 2 (f)} ) ( 1 ^+^ 1 over {σ j 2 (f) } ) }
where the subscripts i and j represent different sensors. If σ i 2 = σ j 2 = 20$, then γ if 2 = 0.907 by (5.9) In our experiment we approximate
σ i 2 (f) ~ {S i (f)} over {N i (f) }
where Si (f) and Ni (f) are the power spectra of the signal of interest and the pre-event noise, respectively, and exclude any coherence estimates from further processing when σ 2 (f) < 20.
The estimates of γij2 (f0 ) were examined for all station pairs at 2.5 Hz increments for 2.5 Hz ≤ f0 ≤ 50 Hz. For each frequency f0 , all estimates $γ ij 2 (f0 ) for all events of a particular body wave and component were collected, converted to a more normal distribution ($arc tanh ~ sqrt {γ ij 2 (f0 )} ), and sorted by distance between measuring points. All estimates for the same distance were averaged. A weighted cubic spline fit to these average estimates of γ ij 2 (f0 ) was calculated as a function of distance. The results of the spline fits were combined for all f0 . These are presented in contour plots for the three components fo the P wave (Figure 7) and the S wave (Figure 8). These plots give an average value of MSC as a function of frequency and distance. Values for which γ 2 (f, d ) < .3 are not contoured because they are not significantly different from zero at the 95% confidence level. Each of the contour plots shows, as expected, that γ 2 (f, d ) increases as f and d decrease. The graph at the right side of each contour plots gives the signal-to-noise ratio for each event.
The most coherent signals are the vertical components of the compressional waves (which are nearly vertically incident). These MSC estimates have average values of γ 2 (f, d ) ≥ 0.6 when f = 30 Hz and d = 50 meters and γ 2 (f,d ) ≥ 0.6 when f ≤ 10 Hz and d = 500 meters. The least coherent signals are the transverse component of the compressional waves. Since the transverse component of the compressional wave is primarily scattered energy, the average measurements never exceed a value of γ 2 = 0.5 for the whole range of frequency and distance observed in this experiment. The radial component of the compressional wave should measure part of the direct wave as well as some scattered energy. The averaged coherence plot for this component does indeed have similar but lower values than the data from the vertical components.
The shear waves are also nearly vertically incident and should theoretically be separable into SH and SV polarized components. The SH polarization is observed on the transverse component, while the SV polarization is observed on the radial and vertical components. The highest averaged coherences for the S waves are found on the vertical component. The next most coherent component of the S waves is the transverse component, for which γ2 (f,d ) is significantly different from zero for f ≤ 20 Hz at 100 meters. The least coherent component of the shear waves is the radial component. It is important to note that at distances beyond 200 meters and at frequencies greater than 10 Hz, the average coherence is not significantly different from zero coherence for all components of the shear waves.
7. Earthquake Source Parameter Estimation
The primary use of spectra in seismic source studies is to determine various parameters to characterize properties of the seismic source. The far field displacement amplitude spectrum of any reasonable kinematic model of earthquake sources will be bounded by a constant value at low frequencies and is proportional to a negative power at high frequencies. This type of model can be described by three parameters [Aki and Richards, chapter 14, 1980]. The low frequency displacement amplitude level, Ω0, is proportional to the seismic moment of the event. The corner frequency fc is used to estimate a seismic source dimension. The high frequency rolloff rate, given by the exponent N , is dependent on the type of source model. One possible way to model the displacement amplitude spectrum, Ω (f ), of seismic body waves is to use the following equation
Ω (f) = {Ω 0} over {1 + ( f / { f c })N}
which is the same form as a Butterworth filter response [Oppenheim and Schafer, 1975].
One obvious problem that needs to be considered is how does the variability of the spectrum found in Section 6 affects the stability of the spectral parameter estimation. An estimate of seismic moment M can be calculated by M = Ω0F, where F depends on the source-receiver distance, orientation of fault plane, and other physical properties of the medium (see Chapter III). In studies of the aftershocks of the Oroville earthquakes, Boatwright [1984] observed a coefficient of variation Cv = (standard deviation/mean ) of M to be 0.25 to 0.55 for selected events, while Fletcher et al. [1984] found an average Cv of 0.6. In studies of the July 1983 Miramichi earthquakes, Cranswick et al. [1985] found a Cv of about 0.4 for their moment estimates. To test the variability of the spatial parameters from the stations in the small aperture array, all the power spectra and their error estimates were converted into displacement amplitude spectra. The value of Ω0 was estimated from the low frequency part of the displacement amplitude spectrum. Once Ω0 is known then (5.10) can be linearized by taking the log of both sides of the equation, where fc and n are to be determined using a weighted least-squares-fit to the spectrum. The weights used are the jackknifed variance estimates of the spectrum. The results of these fits for the P wave vertical component are presented in Table 2 and for the S wave transverse component in Table 3.
In both tables, the estimates of Ω0 and fc are log averages of the individual estimates for each station, so the errors presented are as coefficients of variation determined by one standard deviation of the log of Ω0 and fc . The errors in the estimates of Ω0 for both the P and S waves are generally 20% to 30% of the mean value, and the errors of fc are between 10% and 20% of the mean. For both phases, the higher signal-to-noise ratio of the spectra gave the more stable estimates of the spectral parameters. The spectral parameters of the P wave of event 5 were not calculated because several recorders only triggered on the S wave and there were too few data.
Recently Snoke [1987] derived an alternative way of determining a corner frequency fs , using the integral of the velocity amplitude spectrum and assuming the Brune source model value N = 2. This method is claimed to be more stable for determining stress drops than the traditional methods of picking Ω0 and fc. The log-averaged estimates of fs have the same stability as the log-averaged estimates of fc but the mean fs value tends to be biased to lower values than the mean value of fc. A misfit of the model, mfit1 in Tables 2 and 3 to the log spectrum was calculated using a Euclidean norm with the parameters Ω0, fc , and N in (5.10). The same model misfit criterion is used for the estimates fs , and Ω0 with N = 2. The model which solves for fc and N gives, on average, a three times smaller misfit to the measured displacement amplitude spectrum than the Snoke method which has the implicit assumption of N = 2. An example of a displacement amplitude spectrum for the transverse shear wave component, from event 8 recorded on station AZE, is shown in Figure 9 along with the ±95% confidence limits.
The moment estimates for the previously mentioned aftershock studies were calculated using handpicked parameters. When the spectral parameters for the current data set were handpicked, using the same procedure as Fletcher, et al. [1987], the Cv for Ω0 ranged from 0.2 to 0.8 and for fc the Cv was from 0.1 to 0.5. It is possible that the instabilities in the moment estimates of previous aftershock studies are due to site variations and parameter estimation errors.
8. Discussion
The spectral and coherence estimates presented are similar to measurements made in other locations. Tucker et al. [1984] compared spectra by calculating displacement amplitude spectral ratios between stations. They recorded earthquakes with different azimuths and angles of incidence in the frequency band 1 to 30 Hz, and found that the spectral ratios vary by no more than a factor of 3 for their two surface hard rock sites. King and Tucker [1984] also calculated spectral ratios for their data with stations across a shallow valley filled with sediments overlying hard rock. In that experiment, amplifications as high as ten were found between hard rock sites and sediment sites, while adjacent sites in the sediments with less than 100 meter spacing showed spectra ratios of up to x5.
To make a direct comparison between the prior results and the current data set, amplitude spectral ratios were calculated for the vertical components of the P waves and the transverse components of the S waves. Each trace in Figures 10a and 11a shows the average spectral ratio as a function of frequency for one event. In Figures 10b and 11b, each trace plots the maximum spectral ratio for each event as a function of frequency. The average spectral ratios for this array can be as high as a factor of two, while the individual peak maximum values reach almost a factor of seven. The station pair with the maximum spectral ratio for an event can be different for different frequencies. This is caused by the intermixing of the spectra which is shown in Figure 5b. Careful examination of all events show that the station pair with the maximum spectral ratio at a specific frequency changes from event to event. The spectral ratio also varies from event to event for each particular station pair. The average spectral ratios are comparable to the hard rock site results of Tucker et al. [1984], but the peak values meet or exceed the peak ratios which King and Tucker [1984] found for adjacent sediment sites. For reference, if this were a true stationary homogeneous process then the expected 1% quantile is approximately 2.7 where the maximum ratio is known as Cochran statistic [Pearson and Hartley, 1969].
McLaughlin et al. [1983] recorded a nuclear explosion on an eight station seismic array and examined the spatial coherence of the seismic waves. Their data show a high degree of variability in the coherences as a function of frequencies. There was significant incoherence energy above 6 Hz for P waves and above 4 Hz for S waves. They found that the correlations between stations is higher in the direction transverse to propagation. McLaughin [1982] installed another array with 12 stations near the San Andreas fault in central California which measured two local earthquakes. In this experiment a coherency maximum was found in the 8 to 13 Hz band for stations spaced 170 to 250 meters apart. Outside this frequency band, the coherency was variable. It is possible that the variability of coherence below 8 Hz is being affected by the response of the 4.5 Hz geophones that were used.
The data presented in this paper have the advantage of a good distribution of source azimuths so that the average properties of the site can be examined. The MSC estimates have a high degree of variability as a function of frequency, superimposed on a general trend of decreasing MSC, as the interstation spacing and the frequency increase. Examination of all P wave MSC estimates on the vertical components shows that station pairs with γ2 (f ) ≥ 0.8 could be found for all f ≤ 45 Hz, and f ≤ 20 Hz for transverse S waves. Unlike the results of McLaughlin et al. [1983], the MSC estimates were not correlated with radial or transverse distance. For each event there might be a small apparent correlation, but when averaged over all events, the correlations vanish. The MSC estimates also showed no correlation with azimuth or station pairs and are independent of the properties of the earthquake source focal mechanism when the data from all events is examined.
The results presented here have strong implications about the ability to resolve features of local seismic sources. This study area was picked because the geology and topography should not introduce large distortions to the seismic waves that propagate through it. Other studies show that this area of the northern Peninsular Ranges has low attenuation [Hough and Anderson, 1988a,b] and that the Piñon site has a high fmax value [Fletcher et al., 1987]. However, in spite of these properties, the observed body waves are incoherent at 500 meters station separation for frequencies greater than 15 Hz for P waves and greater than 10 Hz for S waves. The only coherent (γ 2 > 0.8) frequencies for both body wave phases are well below the observed displacement amplitude corner frequencies of each respective phase. The breakdown in the coherence with increasing frequency is caused by amplitude variations in their power spectra as shown in Figure 5b. The high frequency part of the spectra have similar properties in a statistical sense. The consistency of the source parameter estimates calculated from the body wave spectra show that the spectra have general properties, (fc , fs , Ω0 , N ), which are stable over the spatial extent of the array. The variations of these estimates are smaller than those in experiments where the stations are placed kilometers apart. The differences in source parameter estimates observed in other experiments could be attributed to the sampling of different parts of the source's focal sphere, path or site effects, or the mispicking of Ω0 and fc .
In this study, the earthquakes have source dimensions about 0.1 to 0.3 km, and are small compared to the 15 to 40 km source-receiver distance. If the geology of the region is homogeneous, the seismic body wave will propagate effectively as plane waves across the array. The body waves, however, have measurable variability in the power spectra, and lose coherence when the station spacing is greater than 300 meters for wavelengths shorter than 300 to 400 meters. These observations can be explained by some types of scattering mechanism.
One type of scattering model that is widely used is a random velocity perturbation applied to a homogeneous velocity field. The perturbation field has some correlation length a and might have, for example, Gaussian, exponential or a self-similar distribution [Chernov, 1960; Frankel and Clayton, 1986]. Studies of the velocity structure of the southern California region suggest scale lengths on the order of tens of kilometers are present, with velocity perturbations less than 10% [Hearn and Clayton, 1986a,b]. If the scale length of 10 km is assumed, then this allows certain tests of the appropriate scattering model. McLaughlin and Anderson [1987] showed that a Gaussian medium has the property that for fixed interstation distances rij << a the γij2(f ) will increase with frequency. Chernov [1960] found that the transverse autocorrelation for amplitude and phase approaches unity for station spacings less than the scale length of the medium. Frankel and Clayton [1986] showed that the exponential medium has similar properties to the Gaussian medium for rij << a. These theoretical and numerical results are clearly inconsistent with the data presented here.
Frankel and Clayton [1986] also studied a Von Karman medium, which has the property of being "self similar" in the sense that the standard deviation of the medium remains constant over equal logarithmic intervals of wavenumber. For a two-dimensional model with a = 10 km and a standard deviation of 5%, their average band passed cross-correlation estimates at 200 meter spacing are .7 at 15 Hz and .58 at 30 Hz. The band-passed cross-correlation is equivalent to the coherence estimates. However, if this type of medium represents the earth in the locality of this array, then a larger standard deviation is needed to fit the data.
Another type of scattering model is an undulating interface separating the weathering layer from a homogeneous crystalline rock basement. Theoretical models of low velocity materials overlying curved boundaries have been studied by Bravo et al. [1988] and Sanchez-Sesma et al. [1988] among others. These calculations show similar variations in amplitude, as a function of frequencies, as were observed in this experiment. Bravo et al. [1988] show examples of synthetic seismograms that show characteristics similar to those observed on this array. This type of model is perhaps the preferable interpretation.
We examined the coherence between a set of borehole seismometers located at a site 25 km away. This site, KNW, has very similar geological and topographical characteristics to those found at PFO. The borehole seismometers are located at depths of 300 (KN1) and 150 (KN2) meters. There is also a set of surface instruments (KN3). The data consisted of compressional waves from 13 local earthquakes measured on the vertical components. The values for MSC (f), averaged for all events, are shown in Figure 12. Also plotted is the average MSCij(f) from the PFO array for station spacing Γ ij = 127 meters. At high frequencies, the vertically separated sensors have greater coherence than surface sensors with the same spacing. The coherence between KN1-KN2 is on average greater than that for KN2-KN3 which indicates that most of the signal distortion is occurring in the top 150 meters. Neither the source of the low coherence between KN1 and KN2 at 3 Hz nor the very pronounced scalloping effect in the KN2-KN3 coherence has been resolved yet. These effects are probably caused by the reflections of the compressional waves off the free surface.
9. Conclusions
During the operation of this small aperture seismic array at Piñon Flat, California, eight earthquakes were well recorded. The site was chosen for its homogeneous granitic geology and its planar topography. The sensors and analog filters in the data loggers were all calibrated and had identical response functions. Power spectra were calculated for each body wave arrival. The spectra for each seismic phase and sensor component were log averaged together. A two standard deviation confidence limit was calculated, which is, on average, a multiplicative factor of three times the mean spectral level. Magnitude squared coherences were calculated for all recordings at different sites. The magnitude squared coherence estimates were highest for the P wave arrivals on the vertical component and lowest for the P wave recorded on the transverse component. Coherence for both the P and S waves decreases rapidly as frequency increases and slowly as distance increases. The source parameter estimates Ω0 and fc for these data are more stable that estimates from previous studies with larger station spacings.
The results here suggest that, even for optimal sites, the measured seismic wavefield can be distorted over scale lengths of 500 meters. For sensors located within this distance, the general properties of the spectra are consistent but the details are uncorrelated at higher frequencies. Even though the surface geology and topography are constant, a weathering layer with variable depth could cause this phenomenon. One possible way to avoid this problem at this site is to place the sensors in boreholes. If the sensors are constrained to be located on the surface of the earth, then the study of seismic source properties should be approached from a statistical perspective. High resolution imaging of the seismic source can only be accomplished if the effects of scattering and attenuation of the propagation medium can be removed from the measured data.
Appendix: Multiple Window Spectrum Analysis
(from Vernon et al., manuscript in preparation, 1989)
Prolate Spheroidal Sequences and Wavefunctions
It is convenient to separate the bias of a spectrum into two types: broad-band bias, which occurs in the principal frequency domain (- 1/2 , 1/2 ] outside of a smaller domain (-W,W), and local bias, which is confined to the region (-W,W). W is an arbitrary frequency, although generally W < < 1/2. A good criterion for the choice of a data window is minimization of broad-band bias, or, equivalently, confinement of bias to a limited and specified spectral region. It can be shown that of all possible functions of unit energy which are limited to a finite time, the prolate spheroidal wavefunctions minimize the energy outside of a finite bandwidth, and hence minimize broad-band bias. The properties of these functions were extensively examined in a series of papers in the Bell Systems Technical Journal in the early 1960's [Slepian and Pollak, 1961; Landau and Pollak, 1961,1962; Slepian, 1964]. The first two of these contain useful background on the continuous-time prolate spheroidal wavefunctions and their role in signal theory. A more recent series of papers by Thomson [1977a,1977b] details the properties of the prolate functions as data windows. The superior properties of this family of windows is now well-known, and they should be used routinely when the section-averaging method is applied.
In the sequel, we will use the following notation. The data x(n) are assumed to consist of N contiguous samples of a stationary, real, ergodic, zeromean, Gaussian time series. Some of these conditions can be relaxed (in particular, the requirement of a real and zeromean time series), and non-Gaussianity rarely stops the data analyzer. However, it should be noted that the Gaussian assumption is most critical when inferences are drawn from a spectral estimate (e.g., a confidence limit is placed on a power spectrum). The data sample interval is unity, so that frequency is measured in Nyquists and occupies the principal domain (-1/2,1/2], with f = 1/2 corresponding to the Nyquist frequency. It is assumed that adequate pre-filtering has been performed to eliminate aliasing. The true power spectrum of the sampled time series will be denoted by S(f), while an estimate of it is S(f).
The discrete Fourier transform (DFT) is defined by
y ( f ) = sum from n=0 to {N-1} ~ e sup {-i2 pi f(n-(N-1)/2)} ^ x(n)
where time is written in centered form to keep with the standard notation for the discrete prolate functions [Slepian, 1978]. >From (A1), it is clear that time is a discrete variable due to the sampled nature of the data, but frequency may be a continuous variable on the principal domain (-1/2,1/2]. The usual FFT algorithms yield frequency on an even grid with a spacing of N -1, but this is a computational convenience, not a mathematical requirement.
For a stationary process, the inverse Fourier transform associated with (A1) may be written in terms of the stochastic or Cramer representation
x(n) = int from {- 1/2} to 1/2 ~ e sup {i2 pi v (n-(N-1)/2)} ^ dZ(v)
Equation (A2) is in the form of a Stieltjes integral, and dZ(f) is an orthogonal increment measure, in the sense that dZ(f) and dZ(v) are statistically uncorrelated if f and v are distinct frequencies. However, they may be dependent, since dZ(f) = dZ* (-f) for a real process. Extensive discussions of the statistical basis for the Cramer representation may be found in Doob [1953] and Brillinger [1980]. If the process has zeromean, we have
E [dZ(f)] = 0
E [dZ(f)dZ sup {*} (f)] = S(f)^df
where E denotes the expected value. This means that dZ(f) is a measure of the true spectrum, which can only be known exactly if infinite data are available. Its moments, and especially the second moment or power spectrum, are of immediate interest in spectral analysis.
Combining (A1) with (A2) and interchanging the orders of summation and integration yields
y(f) = int from {- 1/2} to {1/2} ^ sum from {n=0} to {N-1} ^ e sup {i2 pi (v-f)(n-(N-1)/2)} ^ dZ(v)
where
sum from {n=0} to {N-1} ^ e sup {i2 pi (v-f)(n-(N-1)/2)} = {sin N pi (v-f)} over {sin pi (v-f)}
is called the Dirichlet kernel, so that the equation reduces to
y(f) = int from {- 1/2} to {1/2} ^ {sin N pi (f-v)} over {sin pi (f-v)} ^ dZ(v)
This is the fundamental equation of spectral analysis. It describes a projection operation from the infinite stationary process represented by dZ(f) onto the finite sample y(f) using the Dirichlet kernel, which serves as a window through which the true spectrum is viewed. The kernel is very broad for small Nf, so that only a smeared version of the true spectrum is available, but becomes increasingly delta function-like as Nf -> ∞ , and a "good" spectrum is achievable for time series which are long in the usual sense. Note that (A3) does not have a unique inverse; it is like an underdetermined matrix problem. More formally, it is equivalent to a linear integral equation of the first kind, and can be solved approximately in a number of ways. It is the solution of (A3) which leads to multiple window spectral analysis. It should be remembered that it is the statistics of dZ(f) that are of interest, not those of the DFT y(f).
An approach to the spectral analysis problem is suggested by the standard methods used to solve first kind integral equations [e.g., Mikhlin, 1964]. The integral equation
y(x) = int from a to b ^ dx ′ ^ K(x,x ′ ) ^ Z(x ′ )
where the symmetric kernel $K(x, x ′ )$ has a set of eigenvalues $"{" λ sub {k} "}"$ and orthonormal eigenfunctions $"{" psi sub {k} "}"$ given by the solutions of
λ sub {k} psi sub {k} (x) = int from a to b ^ dx ′ ^ K(x,x ′ ) ^ psi sub {k} (x ′ )
yields Z hihat (x) = sum from k to nothing ^ λ sub {k} sup {-1} ^ [ int from a to b dx ′ y(x ′ ) psi sub {k} (x ′ ) ] ^ psi sub {k} (x)
where the term in brackets is the k-th Fourier-Bessel expansion coefficient. However, the eigenvalues are known to decay exponentially for large k, so that higher order terms grow rapidly. It is usually necessary to approximate the solution by truncation and/or by using weights.
In order to solve (A3), the characteristics of the eigenfunctions of the Dirichlet kernel must be known. These were examined thoroughly by Slepian [1978], extending earlier and closely related work on the continuous-time prolate wavefunctions. Slepian [1978] has shown that the Dirichlet kernel eigenfunctions, which are called discrete prolate spheroidal wavefunctions (DPSW), and sometimes Slepian functions, are doubly orthogonal, satisfying the eigenvalue problems
int from {-W} to W ^ dv ^ {sin N pi (f-v)} over {sin pi (f-v)} ^ U sub {k} (N,W;v) mark = λ sub {k} (N,W) ^ U sub {k} (N,W;f)
and
int from {- 1/2} to 1/2 ^ dv ^ {sin N pi (f-v)} over {sin pi (f-v)} ^ U sub {k} (N,W;v) lineup = U sub {k} (N,W;f)
where
int from {-W} to W ^ df ^ U sub {j} (N,W;f) U sub {k} (N,W;f) mark = λ k delta sub {jk}
int from {- 1/2} to 1/2 ^ df ^ U sub {j} (N,W;f) U sub {k} (N,W;f) lineup = delta sub {jk}
The parameter W is a frequency, typically $O(N sup {-1} )$, which defines the local or interior domain; the area outside of (-W,W) is called the exterior domain. Note that W is a free parameter, so that a complete set of eigenvalues and eigenvectors is obtained for each choice of the inner bandwidth 2W, or equivalently, the time-bandwidth product NW. The eigenvalues λ k in (A5a) are distinct, real, and finite in number because the kernel is degenerate, and are ranked so that
1 > λ sub {0} > λ sub {1} > ... > λ sub {N-1} > 0
The first 2NW eigenvalues are very nearly 1, and they decay exponentially to zero at larger order. The eigenvalues give the fractional energy concentration of their corresponding eigenvectors in the interior domain. Of all possible functions which are the Fourier transform of a time-limited sequence, the DPSW of lowest order U0 (N,W;f ) has the largest energy concentration in (-W,W).
The Fourier transforms of the DPSW are called discrete prolate spheroidal sequences (DPSS) or Slepian sequences
v sub {n} sup {(k)} (N,W) = 1 over {epsilon sub {k} λ sub {k} (N,W)} int from {-W} to W ^df ^U sub {k} (N,W;f) ^ e sup {i2 pi f(n-(N-1)/2)}
for k = 0, ..., N -1, where Εk is 1 for k even and i for k odd. An equivalent form exists on the complete interval (- 1/2 , 1/2 ]. The DFT of the DPSS yields the DPSW
U sub {k} (N,W,f) = epsilon sub {k} ^ sum from {n=0} to {N-1} ^ e sup {-i2 pi f(n-(N-1)/2)} ^ v sub {n} sup {(k)} (N,W)
Note that the DPSW and DPSS are both real because of the time-centered form of the Fourier transform in (A1) and (A2).
Computation of the DPSS and DPSW is relatively straightforward. Thomson [1982, appendix] discusses an approximate solution for the DPSS using a relation between the discrete and continuous time prolate functions that works well with double precision arithmetic on a 32 bit computer to time-bandwidth 5. Another approach using a Toeplitz eigenvalue problem satisfied by the DPSS is discussed in Thomson [1982], and yields adequate results for moderate values of N. Much better results have been achieved by treating an analogous tridiagonal eigenvalue problem which yields the DPSS and is derived in Slepian [1978]. For this case, the eigenvalues bear no relation to the ones in (A5)-(A6), so they are obtained by integration of (A7) over the inner and outer bandwidths. In all instances, the DPSW follow from (A7).
Univariate Multiple Window Spectral Analysis
Due to the dual orthogonality properties of the prolate functions, and motivated by the bias concentration properties of the DPSW, we choose to solve the integral equation (A3) on the interior domain (-W,W). This means that the interior bandwidth 2W (or, equivalently, the time-bandwidth product NW) is a free parameter which may be chosen from consideration of the desired resolution and which ultimately determines the variance of the estimate. Since first kind integral equations do not have unique inverses, the existence of free parameters should not cause concern. The solution dZ (f) to (A3) will be expanded in DPSWs over the interior domain (-W,W) using weight functions ak
dZ hihat (f) = sum from k to nothing ^ a sub {k} (fo ) ^ U sub {k} (f-fo )
This solution is valid on fo -W ≤ f ≤ fo +W, and the upper limit to the sum will be specified later. We will define the Fourier-Bessel expansion coefficients of y(f) to be
yk (fo) mark = {λ sub {k}} sup {-1} ^ int from {fo -W} to {fo +W} ^ df ^ y(f) ^ U sub {k} (f-fo)
= {λ sub {k}} sup {-1} ^ int from {-W} to W ^ ds ^ U sub {k} (s) ^ y(s+fo)
The normalization in (A9a) is such that the true spectrum S(f) will be obtained for a white noise process, where bias is not a consideration. To see this, substitute (A3) for y(f) in (A9a) to yield
yk (fo) = int from {- 1/2} to 1/2 ^ U sub {k} (f) ^ dZ(f+fo)
so that
bold E ^ [|yk (f)| sup {2} ] = bold E ^ [|dZ(f)| sup {2}] = S(f) ^ df.
The expansion coefficients ak in (A8) are analogous to the Fourier-Bessel coefficients in (A9a). A solution of (A3), called the high resolution estimate and valid for fo - W ≤ f ≤ fo +W or f - W ≤ fo ≤ f + W, may be written
{dZ hihat} sub {h} (f;fo ) = sum from k to nothing ^bk ^ yk (fo ) ^ U sub {k} (f-fo )
where the bk will be chosen to yield the correct spectrum for a white noise process. For the moment, we are not weighting the terms in a data adaptive manner to minimize bias.
The {yk} are computed numerically using the DFT (A1) after windowing with the DPSS. To verify this, substitute (A1) into (A9a) to get
yk (fo ) mark = {λ sub {k}} sup {-1} ^ int from {fo -W} to {fo +W} ^ df ^ U sub {k} (f-fo ) ^ sum from {n=0} to {N-1} ^ e sup {-i2 pi f(n-(N-1)/2)} ^ x(n)
= {λ sub {k}} sup {-1} ^ sum from {n=0} to {N-1} ^ e sup {-i2 pi fo (n-(N-1)/2)} ^ x(n) ^ int from {-W} to W ^ ds ^ U sub {k} (s) ^ e sup {-i2 pi s(n-(N-1)/2)}
Since the DPSS and DPSW are real when the time-centered form of the Fourier transform is used, the last part of this expression is the complex conjugate of (A5)
int from {-W} to W ^ ds ^ U sub {k} (s) ^ e sup {-i2 pi s(n-(N-1)/2)} = λ sub {k} ^{v sub {n}} sup {(k)} (N,W) ^/^ epsilon sub {k}
so that
yk (fo ) = {epsilon sub {k}} sup {-1} ^ sum from {n=0} to {N-1} ^ e sup {-i2 pi f(n-(N-1)/2)} ^ {v sub {n}} sup {(k)} (N,W) ^ x(n)
The expansion coefficient (A12) will be called the k-th eigentransform, while its square will be called the k-th eigenspectrum and is denoted by Sk (f). It is a direct estimate and is χ 22 distributed.
From the bias properties of the prolate functions, there is little to be gained by using more than the first K = 2NW DPSS as data windows. This determines the upper limit to the sum in (A11). A high resolution spectral estimate follows by squaring (A11) and normalizing, implicitly invoking the ergodic hypothesis
{S hihat} sub {h} (f;fo ) = 1 over N ^ | sum from {k=0} to {K-1} ^ bk ^ yk (fo ) ^ U sub {k} (N,W;f-fo ) | sup {2}
Note that the high resolution estimate has a new free parameter, since it is valid both for fo - W ≤ f ≤ fo + W and for f - W ≤ fo ≤ f + w. This property will be exploited later. It is also a statistically inconsistent estimate, since its variance does not decrease as N increases, and is χ 22 distributed. To obtain statistical consistency, we take its average over the interior domain to get the averaged spectrum S(fo )
S hihat (fo ) mark = 1 over {2W} ^ int from {fo -W} to {fo +W} ^ df ^ {S hihat} sub {h} (f;fo )
= 1 over K ^ sum from {k=0} to {K-1} ^ λ sub {k} ^ {bk} sup {2} ^ |yk (fo ) | sup {2}
where the orthogonality properties of the DPSW are exploited to perform the integration. The coefficients bk follow from the requirement of correct normalization for a white noise process. The eigenspectra are already normalized so that
$bold E ^ [{S hihat} sub {k} (f)] = S(f)$, and $bk = sqrt {{λ sub {k}} sup {-1}}$
to yield that result for the averaged spectrum. This spectral estimate has 2K degrees of freedom, since it is the average of K eigenspectra with 2 degrees of freedom apiece.
Due to the poor bias properties of higher order eigenspectra when the spectrum is (relatively) small, the straight averaging used in (A14b) is not recommended. Instead, we may generalize (A13) to
{S hihat} sub {h} (f;fo ) = A over N {| sum from {k=0} to {K-1} ^ dk (f) yk (fo ) U sub {k} (f-fo )| sup {2}} over {sum from {k=0} to {K-1} ^ |dk (f)| sup {2}}
where the weights dk(f) are computed adaptively from the data in a manner which minimizes the broad-band bias. The scale factor A will then be chosen to yield the correct result for a white noise process in the usual way. The denominator in (A15) is required because the weights are now frequency dependent; we could equivalently insist that the sum of squares of the weights be 1. Considering the definition of the averaged spectrum (A14a) and assuming that the spectrum does not change too rapidly in the interior domain about fo , so that dk(f) ~dk (fo ), we obtain
S hihat (fo ) = A over K {sum from {k=0} to {K-1} ^ λ sub {k} {dk} sup {2} (fo ) |yk ( fo )| sup {2}} over {sum from {k=0} to {K-1} {dk} sup {2} (fo)}
The weights are computed by minimizing the squared difference between the true and estimated Fourier-Bessel coefficients. The latter are given by (A10), while the true Fourier-Bessel coefficients (which cannot be observed) are given by the definition as
dZ sub {k} (fo ) = sqrt {{λ sub {k}} sup {-1}} ^ int from {-W} to W ^ U sub {k} (v) ^ dZ(v+fo )
The normalization of (A17) is again correct for white noise, as can be verified by taking its expected value. The weights will be chosen to minimize
{epsilon sub {k}} sup {2} = |dZ sub {k} (f) ^-^ dk (f) yk (f)| sup {2}
at each frequency f. Using (A10) and (A17), this may be re-expressed as
{epsilon sub {k}} sup {2} = |( sqrt {{λ sub {k}} sup {-1}} ^-^ dk (f)) ^ int from {-W} to W ^ U sub {k} (s) ^ dZ(s+f) ^-^ dk (f) ^ int from {- 1/2} to 1/2 ^ U sub {k} (s) ^ dZ(s+f) | sup {2}
where the cut integral excludes the interior domain (-W,W). This may be expanded, and since the cross terms vanish, the expected values are
E [| int from {-W} to W ^ U sub {k} (s) ^ dZ(s+f)| sup {2}] = λ sub {k} ^ S(f) ^ df
E [| int from {- 1/2} to 1/2 ^ U sub {k} (s) ^ dZ(s+f)| sup {2}] = bk (f) ^ df ~ ≤ ~ σ sup {2} (1- λ sub {k} ) df
where σ2 is the data variance given by
σ2 = 1 over N sum from {n=0} to {N-1} ^ x sup {2} (n)
The second term is the broad-band bias (by definition), and the bound may be obtained with Schwartz's inequality. It is an equality for a white noise process, suggesting the use of prewhitening before data processing. The squared error reduces to
{epsilon sub {k}} sup {2} ~ approx ~ ( sqrt{{λ sub {k}} sup {-1}} ^-^ dk (f)) sup {2} ^ λ sub {k} ^ S(f) ^+^ {dk} sup {2} (f)^ bk (f)
Minimizing this in the usual way yields the approximate weights
dk(f) = { sqrt {λ sub {k}} S(f)} over {λ sub {k} S(f) + bk (f)}
Note that dk(f) = bk in the absence of broad-band bias (Bk(f) =0), recovering the earlier result (A14b). Considering (A16), this suggests taking
A = sum from {k=0} to {K-1} {λ sub {k}} sup {-1}
As the broad-band bias rises, the weights decrease, downweighting portions of the eigenspectra which are dominated by bias. To compute the weights, the true spectra in (A18) are replaced by the estimate (A16), yielding the recursive equation
sum from {k=0} to {K-1} {λ sub {k} ^ ( S hihat (f) - λ sub {k} A over K {S hihat} sub {k} (f))} over {( λ sub {k} S hihat (f) + {B hat} sub {k} (f)) sup {2}} = 0
This may be solved iteratively using standard methods, and the weights are then computed using (A18). An estimate of the equivalent degrees of freedom is given by
v(f) = 2 sum from {k=0} to {K-1} λ sub {k} {dk} sup {2} (f)
Note that this yields the maximum possible value, 2K, only in the absence of bias, and will generally be smaller, especially at low points in the spectrum. A useful indicator of the extent of bias in the estimate can be obtained by averaging (A20) over frequency. If the average is substantially lower than 2K, then either the choice of interior domain frequency W should be increased or the data should be prewhitened more carefully.
The development of the multiple prolate expansion method given here differs slightly from that of Thomson [1982] due to adherence to the correct normalization for white noise processes throughout. In particular, the coefficients bk in (A11)-(A14) are not identical to those in section 3 of Thomson [1982], and the procedure used to get the spectrum recursively in (A19) is also slightly different. In general, these discrepancies are minor, and will not yield different results if proper care is used in data processing.
Multivariate Multiple Window Spectral Analysis
These procedures may easily be extended to the bivariate and multivariate cases. Consider a set of p time series xj(n), n=0, ... N -1, j =1, ... p. We first compute the univariate spectrum and adaptive weights using the methods already described. This includes reshaping of significant spectral lines, as the coherences and transfer functions are ordinarily only defined for the stochastic portion of the spectrum (although new types of coherences and transfer functions could be defined for the line component). The averaged power of the j-th series analogous to (A16) is given by
S hihat sub {jj} ( Ω ) = A over K sum from {i=0} to {K-1} λ i ( e i j ( Ω ) ) 2 | y i j ( Ω ) | 2
where e i j ( Ω ) = d i j ( Ω ) / [ sum from {i=0} to {K-1} ( d i j ( Ω ) ) 2 ] sup {1 over 2}
is a normalized adaptive weight. The subscript on the terms in (A21) refers to the prolate order, while the superscript refers to the series number. An obvious extension of (A21) to the cross-spectrum is
S hihat sub {jm} ( Ω ) = A over K sum from {i=0} to {K-1} λ i e i j ( Ω ) e i sup m ( Ω ) ( y i j ( Ω )) sup * y i sup m ( Ω )
using the adaptive weights from the respective autopowers. The simple coherence for the case p=2
follows immediately as
γ sub {jm} 2 = { | S hihat sub {jm} ( Ω )| 2} over {| S hihat sub {jj} ( Ω ) | | S hihat sub {mm} ( Ω )|}
while the phase is just that of (A22).
References
Aki, K., and P.G. Richards, The seismic source: kinematics, in Quantitative Seismology Theory and Methods, edited by Allan Cox, pp. 932, W.H. Freeman and Company, San Francisco, 1980.
Boatwright, J., Seismic estimates of stress release, J. Geophys. Res., 89, 6961-6968, 1984.
Bravo, M. A., F. J. Sanchez-Sesma, and F. J. Chanez-Garci, Ground motion on stratified alluvial deposits for incident SH waves, Bull. Seism. Soc. Am., 78, 436-450, 1988.
Brillinger, D.R., 1980. Time Series, San Francisco: Holden-Day.
Chernov, L., Wave propagation in a random medium, McGraw-Hill, New York, 1960.
Cranswick, E., R. Wetmiller, and J. Boatwright, High frequency and source parameters of microearthquakes recorded at hard-rock sites, Bull. Seism. Soc. Am., 75, 1535-1569, 1985.
Dibblee, T. W., Geology of the San Jacinto Mountains and adjacent areas, Geology of the San Jacinto Mountains, South Coast Geological Society, 1-47, 1981.
Doob, J.L., Stochastic Processes, New York: Wiley, 1953. .
Efron, B., The jackknife, the bootstrap, and other resampling plans, CBMS-NSF Regional Conference Series in Applied Mathematics #38, Society for Industrial and Applied Mathematics, Philadelphia, 1982.
Fletcher, J., J. Boatwright, L. Haar, T. Hanks, and A. McGarr, Source parameters of aftershocks of the Oroville earthquake, Bull. Seism. Soc. Am., 74, 1101-1124, 1984.
Fletcher, J., L. Haar, T. Hanks, L. Baker, F. L. Vernon III, J. Berger, and J. N. Brune, The digital array at Anza, J. Geophys. Res., 92, 369-382, 1987.
Fletcher, J. B., T. Fumal, H. Liu, and R Porcella, Near-surface velocities and attenuation at two bore-holes near Anza, CA from logging data (judice), 1989.
Frankel, A., and R.W. Clayton, Finite Difference Simulations of Seismic Scattering: Implications for the propagation of short-period seismic waves in the crust and models of crustal heterogeneity, J. Geophy. Res., 91, 64659mi6489, 1986.
Geli, L., P. Bard, and B. Jullien, The effect of topography on earthquake motion, Bull. Seism. Soc. Am., 78, 42-63, 1988.
Hearn, T. M., and R. W. Clayton, Lateral velocity variations in Southern California, I. Results for the upper crust from $P sub g$, Bull. Seism. Soc. Am., 76, 495-510, 1986a.
Hearn, T. M., and R. W. Clayton, Lateral velocity variations in Southern California, II. Results for the lower crust from $P sub n$, Bull. Seism. Soc. Am., 76, 511-520, 1986b.
Hough, S. E., and J. G. Anderson, Attenuation near Anza, California, Bull. Seism. Soc. Am., 78, 672-691, 1988a.
Hough, S. E., and J. G. Anderson, High-frequency spectra observed at Anza, California: Implications for $Q$ structure, Bull. Seism. Soc. Am., 77, 692-707, 1988b.
Jenkins, G. M., and D. G. Watts, Spectral Analysis and Its Applications, Holden-Day, 1968.
King, J. L., and B. E. Tucker, Observed variations of earthquake motion across a sediment-filled valley, Bull. Seism. Soc. Am., 74, 153-165, 1984.
Landau, H.J., and H.O. Pollak, Prolate spheroidal wave functions, Fourier analysis, and uncertainty-II, Bell Syst. Tech. J., 40, 65-84, 1961.
Landau, H.J., and H.O. Pollak, Prolate spheroidal wave functions, Fourier analysis, and uncertainty-III, Bell Syst. Tech. J., 40, 1295-1336, 1962.
Lanzerotti, L. J., D. J. Thomson, A. Meloni, L. V. Medford, and C. G. Maclennan, Electromagnetic study of the Atlantic continental margin using a section of a transatlantic cable, J. Geophys. Res., 91, 7417-7427, 1986.
Lindberg, C., and J. Park, Multiple-taper spectral analysis of terrestrial free oscillations: Part II, Geophys. J. R. Astron. Soc., 91, 795-836, 1987.
McLaughlin, K. L., Spatial Coherency of Seismic Waveforms, unpublished Ph.D. dissertation, 1982.
McLaughlin, K. L., and L. M. Anderson, Stochastic dispersion of short period P waves due to scattering and multipathing, Geophys. J. R. Astr. Soc., 89, 933-963, 1987.
McLaughin, K. L., L. R. Johnson, and T. V. McEvilly, Two-dimensional array measurements of near-source ground accelerations, Bull. Seism. Soc. Am., 73, 349-375, 1983.
Mikhlin, S.G., Integral Equations and their Applications to some Problems of Mechanics, Mathematical Physics, and Engineering, 2nd ed., Oxford: Pergamon, 1964.
Miller, R. G., Jackknifing variances, Ann. Math. Stat., 39, 567-582, 1968.
Miller, R. G., The jackknife -- a review, Biometrika, 61, 1-15, 1974.
Oppenheim, A. V., and R. W. Schafer, Digital Signal Processing, Prentice-Hall, Inc., New Jersey, 1975.
Park, J., C. R. Lindberg, and F. L. Vernon III, Multitaper spectral analysis of high-frequency seismograms, J. Geophys. Res., 92, 12,675-12,684, 1987.
Parcel, R. F., Structure and Petrology of the Santa Rosa Shear Zone in the Piñon Flat Area, Geology of the San Jacinto Mountains, South Coast Geological Society, 139-150, 1981.
Pearson, E. S., and H. O. Hartley, Biometrika Tables for Statisticians, 3rd ed., Cambridge University Press.
Sanchez-Sesma, F. J., F. J. Chavez-Garcia, M. A. Bravo, Seismic response of a class of alluvial valleys for incident SH waves, Bull. Seism. Soc. Am., 78, 83-95, 1988.
Slepian, D., Prolate spheroidal wave functions, Fourier analysis, and uncertainty-IV, Bell Syst. Tech. J., 43, 3009-3057, 1964.
Slepian, D., Prolate spheroidal wave functions, Fourier analysis, and uncertainty-V: the discrete case, , Bell Syst. Tech. J., 57, 1371-1429, 1978.
Slepian, D. and H.O. Pollak, Prolate spheroidal wave functions, Fourier analysis, and uncertainty-I, Bell Syst. Tech. J., 40, 43-64, 1961.
Snoke, J. A., Stable determination of (Brune) stress drops, Bull. Seism. Soc. Am., 77, 530-538, 1987.
Thomson, D. J., Spectrum estimation techniques for characterization and development of WT4 waveguide-I, f2Bell Syst. Tech. J., 56, 1769-1815, 1977a.
Thomson, D. J., Spectrum estimation techniques for characterization and development of WT4 waveguide-II, Bell Syst. Tech. J., 56, 1983-2005, 1977b.
Thomson, D. J., Spectrum estimation and harmonic analysis, IEEE Proc., 70, 1055-1096, 1982.
Tucker, B. E., J. L. King, D. Hatzfeld, and I. L. Nersesov, Observations of hard-rock site effects, Bull. Seism. Soc. Am., 74, 121-136, 1984.
Tucker, B. E., and J. L. King, Dependence of sediment-filled valley response on input amplitude and valley properties, Bull. Seism. Soc. Am., 74, 153-165, 1984.
Tukey, J., Bias and confidence in not quite large samples, abstract, Ann. Math. Stat., 29, 614, 1958.
Wyatt, F., Displacements of surface movements, J. Geophys. Res., 87, 979-989, 1982.
Wyatt, F., and J. Berger, Investigations of tilt measurements using shallow borehole tiltmeters, J. Geophys. Res., 85, 4351-4362, 1980.
