Chapter 6: Frequency Dependent Polarization Analysis of High-frequency Seismograms

• Title
• Abstract
• Chapter 1
• Chapter 2
• Chapter 3
• Chapter 4
• Chapter 5
• Chapter 6
• Figures
• Chapter 1
• Chapter 2
• Chapter 3
• Chapter 4
• Chapter 5
• Chapter 6

ABSTRACT

We present a multitaper algorithm to estimate the polarization of particle motion as a function of frequency from three-component seismic data. This algorithm is based on a singular value decomposition of a matrix of eigenspectra at a given frequency. The right complex eigenvector z (hat) corresponding to the largest singular value of the matrix has the same direction as the dominant polarization of seismic motion at that frequency. The elements of the polarization vector z (hat) specify the relative amplitudes and phases of motion measured along the recorded components within a chosen frequency band. The width of this frequency band is determined by the time-bandwidth product of the prolate spheroidal tapers used in the analysis. We manipulate the components of z (hat) to determine the apparent azimuth and angle of incidence of seismic motion as a function of frequency. The orthogonality of the eigentapers allows one to calculate easily uncertainties in the estimated azimuth an

1. Introduction

The polarization of particle motion as measured by a three-component seismometer has been studied by a number of straightforward methods, most simply by tracing the projection of the motion as a function of time onto a chosen plane of reference. Although useful to illustrate the particle motion of simple arrivals, this practice is qualitative and less useful with complicated signals.

The problem of extracting a particular type of wave (e.g., P, SH, Rayleigh) from a noisy background has been studied by correlation techniques and special filters [e.g., Kanasewich, 1981; Archambeau and Flinn, 1965; Vidale, 1986]. Most of these techniques are designed for time domain analysis and implicitly assume that the waveform has essentially the same polarization over all or most frequencies. Samson [1977, 1983a,b,c] describes a method of estimating the polarization as a function of frequency. This is important for the analysis of seismic records. The seismic waveforms of local and regional distance events are often superpositions of direct, refracted, reflected, and scattered waves, with no guarantee that the polarization or phase is constant in frequency. In the presence of strong scattering, one might not expect a respectable "pure state" polarization at any frequency. Alternatively, coherent addition of scattered waves within the crustal waveguide will produce traveling modes whose signature in extended body wave codas may be a well-defined polarization and phase that varies with frequency. The distinct spectral peaks seen in Chapters IV and V in seismic spectra observed on the Anza Network suggest that waveguide modes may be evident in the complex waveforms of events at eπcentral distances of 100-250 km. Inhomogeneities in the crustal waveguide can lead to scattering and coupling of these propagating modes (see, e.g., Kennett [1986] and Odom [1986] for a description of these effects) which will, in general, cause frequency dependent scattering. In such cases, it is more useful to determine the type of seismic motion from its polarization signature, as in the study of Vidale [1986], than to attempt to isolate phases.

In this paper we develop and demonstrate another algorithm for determining the frequency dependence of the polarization of high-frequency seismic records. We have used multitaper spectral analysis [Thomson, 1982] to estimate the spectral density matrix S(f) of Samson [1983a]. This has several advantages. By employing prolate spheroidal wave functions as tapers (instead of cosine or boxcar tapers) to obtain direct spectral estimates, the elements of the estimated spectral density matrix will be less biased [Lindberg, 1986; Park et al., 1987; Chapter IV]. It is also not necessary to apply a moving average to the density matrix estimate to smooth it; smoothing is obtained by summing the eigenspectra of each component of motion (see equation (6.3)). Using multitapers to estimate the spectral density matrix is more suitable for very short records, such as those which include a single seismic phase. This is because data are not discarded by applying a single bell-shaped taper to the record as discussed in Chapter IV. (A similar method has been independently developed and applied to magnetometer data by Lanzerotti et al. [1986].)

We analyze a number of three-component records of seismic codas. In these observations the source pulse has been dispersed and scattered within the crust. In an idealized picture the shape of the source spectrum is retained in the shape of the coda spectrum, but the spectral phase is randomized by scattering effects. Despite this randomized phase, one might expect the particle motion to retain the polarization behavior of the type of wave motion dominant within a selected frequency band. Polarization analysis in the frequency domain offers an opportunity to characterize the signal better. With three-component data we have potentially three independent polarizations. If scattering is not great, a single polarization will predominate. This assumption is often true for the P wave coda. If, for instance, interaction with crustal structure decouples SH and SV motion, there may be two principal polarizations in the S wave coda. The algorithm we describe in this paper offers a quantitative criterion for identifying the single dominant polarization.

In section 2, our multitaper polarization analysis method is described. We apply the algorithm to a synthetic pulse example in section 3. In section 4 we show examples from the P wave codas of data observed on the Anza Seismic Telemetered Network. Section 5 summarizes our findings. Uncertainty estimates for polarization angles and phases are derived in the appendix.

2. Polarization Analysis With The Multitaper Algorithm

Polarization analysis involves determining the eigenstructure of the spectral density matrix S(f). Suppose one has three-component data recorded in the time domain of the form

x (t) = (x ¹ (t), x ² (t), x ³ (t))

t = n τ ; n = 0, 1, ..., N -1

where τ; is the sampling interval, N τ is the length of the time series, the coordinate system is right-handed, and x¹(t) is the vertical component. If the j th record x^j(t) has the frequency domain representation z^j(f), the spectral density matrix S(f) has components

S_jk(f) = E "{" (z^j(f) ^* z sup k (f) "}"

where E denotes the expectation operator. Samson [1983a] forms an estimate of the spectral density matrix, S(hihat) (f), with components

S hihat_jk(f) = (y^j(f)) ^* y sup k (f)

j,k = 1, 2, 3

where y^j(f) = 1 over {N τ} ~ sum from n=0 to N-1 ~ w_nx^j(n τ ) e sup {-i2 π fn τ}

is a discrete Fourier transform of the j th component of x(t) and {w_n}_{n = 0}^{N - 1} is a chosen data taper. The matrix S(hihat) (f) is then smoothed in the frequency domain by applying a moving average, and the eigenvectors and eigenvalues of the smoothed matrix are found.

To apply the multitaper algorithm to the estimation of S(f), one employs a set of K prolate spheroidal wave function "eigentapers" v_n^(k) (N,W ); k = 0, 1, ..., K -1, which are optimally resistant to spectral leakage from outside a chosen frequency band of width 2W [Thomson, 1982; Lindberg, 1986; Park et al., 1987]. For k = 0, 1, ..., K -1 the spectral estimates (from (4.8))

y _k ^(j) (f) = 1 over {N τ} ~sum from n=0 to N-1 ~ v_n^(k) (N,W)x^j(n τ ) e sup {-i2 π fn τ}

of each component of x(t) can be made. Then a multitaper estimate of the spectral density matrix is

1 over K ~bold M ^H (f) ~cdot~ bold M (f) where superscript H denotes conjugate transpose and

M (f) = left [ ~ matrix { ccol { y ₀ ⁽¹⁾ (f) above y ₁ ⁽¹⁾ (f) above 3dot above y sub K-1 ⁽¹⁾ (f) } ccol { y ₀ ⁽²⁾ (f) above y ₁ ⁽²⁾ (f) above 3dot above y sub K-1 ⁽²⁾ (f) } ccol { y ₀ ⁽³⁾ (f) above y ₁ ⁽³⁾ (f) above 3dot above y sub K-1 ⁽³⁾ (f) } } ~ right ]

The value of K, the number of eigenspectra used, depends on 2W, the width of the frequency band in which the spectral energy at frequency f is concentrated. The K = 2NW τ - 1 lowest order eigentapers possess sufficient spectral leakage resistance to be useful [Slepian, 1983]. To investigate the eigenstructure of S(f), we perform a singular value decomposition

M (f) = U cdot D cdot V ^H,

where U is a K xK unitary matrix of left eigenvectors of M, V is a 3 x 3 unitary matrix of right eigenvectors v_j of M, and D is a K x 3 matrix with D_jj= d_j, j = 1, 2, 3, the singular values of M, and D _ij = 0 for i ≠ j [Golub and Van Loan, 1983]. The polarization vector z (hat) is the right eigenvector corresponding to the largest singular value of the matrix M. It specifies the direction of particle motion at frequency f which contains the largest fraction of seismic energy [Samson, 1983b]. The components of z(hat) can be complex, allowing for phase lags between components. Phase lags between components represent elliptical particle motion. Our ability to identify z(hat) with the principal polarization of motion at f can be qualitatively assessed by comparing the singular values d₁ ≥ d₂≥ d₃. If d₁ >> d₂, d₃, the polarization z(hat) = v(hat)₁ is well determined. We can use the ratio of the singular values to estimate the uncertainty in z(hat) and any quantities we calculate from it. The estimation of the polarization uncertainty follows the derivation of Park and Chave [1984] and is outlined in the appendix. If d₁ ~ d₂>> d₃, there is a strong possibility that coherent seismic motion at f₀ exists at two separate polarizations. The dot product v(hat)₁* cdot v(hat)₂= 0 by virtue of the singular value decomposition, but this orthogonal relationship need not carry over into the seismic polarizations. In an S wave arrival, one expects SV and SH motion to be orthogonal to first order in most situations, but the superposition of other signals (e.g., reflected P arrivals) need not have orthogonal polarizations.

If d₁ >> d₂, d₃, the three-component particle motion x (t) in the neighborhood of frequency f can be represented by the real part of

R z(hat) e sup {+i2 π ft},

where R is the amplitude of motion. We can adjust the phase of z(hat) so that R is real. If there exists a phase &φ; such that z(hat) e ^{i &φ;} is purely real, then all motion described by z(hat) lies along a single line in three space. More generally, particle motion will follow an ellipse confined to the plane spanned by the two real vectors Re ( z(hat) ) and Im ( z(hat) ). If this ellipse is strongly elongated along its major axis, reasonable horizontal and vertical azimuths can be found. If the wave type is known, such as a P wave, then the propagation direction can be determined. Strongly elliptical polarization suggests modelike particle motion (for example, a Rayleigh wave) with a poorly defined angle of incidence.

We can project the particle motion described by the complex unit vector z(hat) onto an ellipse in the horizontal plane which is defined by z_H = z(hat) - ( e (hat)₁ cdot z(hat) ) e(hat)₁, where e(hat) ₁ = (1,0,0). The major axis of this horizontal ellipse is taken to be the principal direction of horizontally polarized motion. To find the azimuth of the major axis, we determine the point of greatest displacement for the projection z_H in the horizontal plane by finding the maximum value of

| Re ( bold z _H e sup {i2 π f t} ) | ²

If the components (z₁ , z₂, z₃ ) of z(hat) are expressed in the form z_j= |z_j| e ^i&φ;j, this is equivalent to finding the maxima of

|z₂| ² cos ² (2 π f t + φ₂) + |z ₃ | ² cos ² (2 π f t + φ ₃ )

The extremes of this expression, remembering

| z₂| ² sin 2 φ₂ + | z ₃ | ² sin 2 φ ₃ = Im [ z₂² + z ₃ ² ]

are found when the phase angle θ defined as θ = 2 π ft takes the values

θ _l = - half roman arg [ z₂² + z ₃ ² ] + { l π } over 2

where l is an integer. Let l be the integer closest to zero which minimizes (6.5), the horizontal displacement, and for which Re (z₁ ) < 0. Define the phase angle θ_H to be the value of θ _l for this l . Once θ _H has been determined, the horizontal azimuth of the major axis Θ_H measured counterclockwise from e(hat)₂= (0,1,0) can be defined as

Θ _H = tan ^-1 left [ {Re (z ₃ e sup {-i θ _H} )} over {Re (z₂e sup {-i θ _H} )} right ] = Re ( tan ^-1 (z ₃ / z₂) )

The range of the arctangent function is 0 < Θ _H ≤ 180 if Re (z ₁ z ₃^* ) < 0 and -180 < Θ _H ≤ 0 if Re (z ₁ z ₃^* ) ≥ 0. If the particle motion is P like, Θ _H can be interpreted as pointing in the direction of the wave source. A representation of an elliptical motion for which Θ _H < 0 is shown in Figure 1.

Another useful quantity is φ ₃ - φ₂= φ _HH, the phase difference between the horizontal components of particle motion. If φ₂- φ₃ ~ 0 or 180, the particle motion is predominantly linear. The value φ₂- φ₃ ~ 90 represents elliptical motion with the major and minor axes oriented along the axes of the instruments. If z₂= ± iz₃, the particle motion is circular, with no definable azimuth. In this case, the uncertainty in Θ _H, given in the appendix, goes to infinity as it is proportional to |z₂² + z₃² | ^-1.

The expressions relating horizontal to vertical motion are similar. We want to find the angle Θ_V made with the vertical by the major axis of the ellipse defined by Re ( z(hat) e^{+2 π f t}). Define the phase angles

θ_m= 2 π f t = - 1/2 arg [z₁² + z_H² ] + { m π } over 2

where m is an integer and z_H² = z₂² + z₃². The phase angle θ _V is the value of θ_m at an m for which the particle motion displacement is maximized. The angle of incidence is

Θ _V = tan ^-1 left [ left | {Re [z ₁ e ^-iθv ]} over {Re [ z _H e ^-iθv]} right | right ]

where Im z_H ≥ 0. The absolute value is taken to restrict Θ_V to lie between 0° and 90°, the usual convention for the angle of incidence (Figure 1). The phase lag between vertical and horizontal motion can also be defined. Define φ_VH = θ_H - φ₁. Since the end points of the major axis of the horizontal motion ellipse correspond to θ_H and θ_H ± π, we can restrict the range of φ_VH to (-90°, 90°).

3. A Synthetic Example

We first illustrate the definitions of Θ_H, Θ_V, φ _HH, and φ _VH in a synthetic example. We constructed a three-component record (Figure 2) from a sum of cosinusoids:

matrix { rcol { x ¹ (n τ ) above x ² (n τ ) above x ³ (n τ ) } ccol { = above = above = } lcol { sum from f=1 to 72 ~ cos ~ left ( π over 80 ~ f over 2 right ) ~sin~ left ( {2 π fn τ} over 2 ~-~ π over 50 ~ f over 2 right ) above sum from f=1 to 72 ~cos~ left ( π over 20 ~ f over 2 right ) ~sin~ left ( π over 80 ~ f over 2 right ) ~sin~ (2 π n τ ~f over 2 ) above sum from f=1 to 72 ~sin~ left ( π over 20 ~ f over 2 right ) ~sin~ left ( π over 80 ~ f over 2 right ) ~sin~ (2 π n τ ~ f over 2 ) } }

where n = 0, 1, ..., N -1 and the sampling interval is τ = 0.004 s. The polarization vector of this signal can be written immediately as

z(hat) (f) = left ( cos left ( {π f} over 80 right ) e sup {i φ _VH} , cos left ( {π f} over 20 right ) sin left ( {π f} over 80 right ) , sin left ( {π f} over 20 right ) ~cdot~ sin left ( {π f} over 80 right ) right )

where φ _HH = 0 and φ _VH = (- π f)/50. Figure 3 shows the results of a multitaper polarization analysis for frequencies 0 ≤ f ≤ 30 Hz. The uncertainties are plotted as one standard deviation error bars in this and succeeding figures. Figure 3b shows the three scaled singular values as a function of frequency. The principal polarization appears well determined. The amplitude spectra for the three components are plotted in Figure 3a. The angles Θ_H and φ_HH are plotted in Figures 3c and 3d. The angle φ_HH is not well determined near zero frequency, as the horizontal signal amplitude is dwarfed by vertical component energy. The apparent horizontal azimuth Θ_H"wraps around" from 180° to -180° at 20 Hz and jumps 180° at 25 Hz. The former jump is obvious; the latter is an artifact of φ_VH passing through 90°. The phase angle φ_HH, estimated from the synthetic record, has a value of 0° or ±180°, to observational accuracies. These values correspond to rectilinear motion and are dependent on the quadrant where the horizontal azimuth is directed. The phase lag φ_VH between vertical and horizontal components is well determined everywhere except very near zero frequency where the horizontal component amplitude vanishes. The ellipticity of particle motion disrupts the linear trend in Θ_V, as shown in Figure 3e. At 25 Hz, φ _VH = 90° and the particle motion is an ellipse with major and minor axes oriented horizontally and vertically, respectively. Therefore Θ _V = 90° at 25 Hz. At higher frequencies, φ _VH > 90, the relative sign of vertical and horizontal motion reverses, and the particle motion ellipse "tips" in an opposite manner relative to its orientation for φ _VH < 90°. This causes the observed 180° jump in apparent horizontal azimuth Θ_H. This example suggests that one should use caution in interpreting the angles Θ_H and Θ_V wherever the particle motion is nearly fully elliptical, i.e., when φ_HH or φ_VH is within 20° of ± 90°.

4. Data Examples

We illustrate this method of determining the polarization as a function of frequency with several examples. We analyzed several waveforms which were recorded on the Anza array after an earthquake that occurred at 0521:39.5 UT, September 9, 1982, with hypocenter positioned at 32.93°N, 115.85°W, and depth 4.2 km. The magnitude M_l was determined to be 4.4. The event was located near Superstition Mountain, California, on the western edge of the Imperial Valley. The earthquake was recorded on only four stations in the array (PFO, KNW, FRD, and CRY; see Berger et al. [1984] or Chapter II for the definitions of these three-letter acronyms) as the event occurred prior to the completion of the array. The hypocenter was roughly 100 km southeast of the array. The m = 1 component is the vertical seismometer output with positive motion defined as up. We choose the m = 2 component so that positive motion points 45° east of north. Positive motion along the m = 3 axis is directed 45° west of north, forming a right-handed coordinate system. Let the angle Θ_H be measured counterclockwise from the primary horizontal axis (m = 2). If the wave propagation is along a straight line connecting the source and receiver, -85° > Θ _H > -98° for the four stations. The first 30 s of recorded motion for this event are shown in Figure 4, along with range and azimuth information (azimuth is measured counterclockwise from N45E). Both S and P arrivals are extended wave trains, although the S energy is more concentrated in time. An interesting feature of this event is the small precursor to the main P arrival, shown in the enlarged detail for stations FRD and CRY in Figures 5a and 5b. This waveform corresponds to a lower crustal phase.

Polarization analysis reveals that the first arrivals have complicated polarization signatures. The time window taken is short (1.6 s), corresponding to a Rayleigh resolution 1/(N τ ) of 0.625 Hz. Analysis using seven 4 π prolate tapers averages energy over a band of width 8/(N τ ), so that all of the estimates shown represent an average over a 5-Hz frequency band. If the true polarization varied significantly over this bandwidth, one would expect Θ_V, Θ_H, φ_VH, and φ_HH to be relatively poorly determined. The results for FRD are shown in Figure 6. The singular values d₂ and d₃ displayed in Figure 6b show local maxima at several places in the spectrum from 0 to 30 Hz. Maxima at 2.5, 7.5, and 14 Hz correspond to boundaries between distinct spectral features (Figure 6a). All the maxima below 25 Hz correspond to frequencies at which one or more of the polarization angles change rapidly. Horizontal motion is roughly rectilinear below 13 Hz, but its azimuth is variable and significantly different from the nominal azimuth of -87°. In fact, the largest amplitude signal, from 8 to 13 Hz, is oriented clockwise 125° from the primary component, a deflection of nearly 40° from the nominal P wave arrival azimuth. The phase lag between horizontal and vertical motion is alternately positive and negative in adjacent frequency bands but is never more than partially elliptical. The angle between vertical and horizontal motion, which can be interpreted in this case as the nominal angle of incidence, varies smoothly with frequency in Figure 6e, with Θ _V ~ 25° - 30° for f < 10 Hz, and Θ _V ~ 15° above 13 Hz.

Figure 7 shows an analysis of the small amplitude P precursor observed at station CRY. The variation of the largest singular value d₁ with frequency shows four frequencies (2.5, 7, 12, and 16 Hz) at which the principal polarization vector is poorly determined, and there is a peak in d₂. Each of these peaks in d₂ occurs where there is an abrupt change in the three-component spectra and in one or more of the polarization angles. Although the estimated uncertainties are larger than those in the last example, the variability among frequency bands is clearly visible in Figures 7c-f. Motion in the horizontal plane is dominantly elliptical below 14 Hz, but particle rotation proceeds in opposite senses in the two frequency bands 2.5 Hz ≤ f ≤ 7 Hz and 7 Hz ≤ f ≤ 14 Hz. The azimuthal angle Θ_H hovers near the value expected for the epicenter (-85°), but our synthetic example in Figure 3 suggests that this may be due to the ~ 90° phase lag between component motions. At higher frequencies, including the substantial spectral peak at 18-20 Hz, the observed horizontal azimuth of particle motion is roughly transverse to the arrival azimuth, as though the energy at these frequencies were SH in character. A better interpretation is in terms of side-scattered P energy, as the vertical azimuth of particle motion Θ_V remains in the 20°-40° range across all frequencies in Figure 7e.

Similar behavior is observed on stations PFO and KNW. The nature of this polarization behavior is quite puzzling. It is unlikely that instrument calibrations are at fault. A timing error among components would result in a linear drift in the relative phase angles, similar to that shown in Figure 3f. There are no poles or zeroes in the instrument response over the frequency region shown. A perturbation in the response filter characteristics would have difficulty mimicking the apparent boundaries between spectral processes. Moreover, we show below that the relative polarization shift from frequency band to frequency band varies greatly within the P coda. This argues for a signal-generated effect rather than an instrument effect. This behavior may reflect the modal structure of an intercrustal head wave in a stratified crust. Another interpretation is in terms of resonant vibrational modes in the earth structure near the receiver. Structure of scale lengths 100-200 m could account for the higher-frequency resonances observed in Figures 6 and 7.

We performed experiments to see if such resonant behavior could be found in the P codas for this event. When the entire coda was used for polarization analysis, the results were poor. The three-component seismogram recorded at station KNW is shown in Figure 8. Figure 9 presents polarization data from the 14-s P wave coda. There appear to be competing signals at nearby frequencies, creating either rapid variations in the polarization, which are difficult to interpret, or else large uncertainties in the polarization. Likewise, the presence of both SV- and SH-polarized energy in the S arrivals made the identification of a "principal" polarization uncertain.

We chose, therefore, to analyze the P codas of these records in successive 2-s (500 sample) segments. We observed what appear to be resonances over 4-6 Hz frequency bands and variations in polarization over time that suggest the arrival of P energy which has been scattered within the crust. The results of a polarization analysis of the first, fourth, and sixth 2-s time segments of the P wave coda recorded at station KNW are shown in Figures 10-11-12. The growth of the "noise" singular values d₂ and d₃ as the time window moves through the coda suggests an increase in scattered energy. The most prominent features in the spectra of the principal polarization components are the spectral peaks near 5 and 14 Hz. Comparison of the values of Θ_H in the time windows indicates that there is a boundary between two distinct spectral processes at 7-7.5 Hz. The 7-14 Hz process is characterized by dominantly rectilinear horizontal motion and steeply vertical particle motion. The relative phase angles φ_VH and φ_HH for the lower-frequency process exhibit more variability. Within a 2-s time window the horizontal azimuth varies only slightly within the 0-7 Hz frequency band, with more shallow vertical angles. Figure 12e shows that Θ_V > 60° in this frequency band, which may indicate SV-converted motion. Particle motion at frequencies greater than 15 Hz bears little relation to the higher-amplitude low-frequency signal and often cannot be interpreted in terms of P-, SV- or SH-polarized motion traveling directly from source to receiver.

The similar frequency dependences of Θ_H and Θ_V in these 2-s time windows contrast with the absence of a clear pattern in the larger time analysis shown in Figure 9. Similar effects are found when records from the other three stations for this event are analyzed. This is not surprising when one notes the large variation of polarization among the three time windows shown in Figures 10-11-12. The azimuth of the epicenter has Θ_H = -92° (i.e., clockwise) from the second component. The horizontal azimuth Θ_H of particle motion is, for 7.5 Hz < f <14 Hz, always oriented more to the south, with values that vary among time windows by 40 or more. At f ≤ 7 Hz, several of the time windows tested were consistent with -92° relative azimuth, but the fourth and sixth segments, shown in Figures 11 and 12, show particle motion whose horizontal orientation is nearly pure east-west. We take this variation as evidence for the arrival of scattered off-azimuth P energy.

A detailed interpretation of these results is beyond the scope of this paper, but we can draw parallels with recent studies of high frequency seismic spectra. Sereno and Orcutt [1985] have shown that the extended P_n wave train observed in ocean bottom seismic data can be modeled by reverberations in the oceanic sediment layer and overlying water column, buttressing their comparison by demonstrating a simple pattern of spectral peaks corresponding to leaky vibrational modes. Bard and Bouchon [1985] have shown spectra from seismic events for which the retrieval of simple source parameters like corner frequency and high-frequency roll-off is contaminated by a high-frequency resonance which they model as a reverberation in the low-velocity surface layer. The apparent polarization resonances observed in the P wave codas of the September 1982 Superstition Mountain event probably argue for an even more complex structure than was postulated in these studies.

The interpretation of the coda using resonance models may offer a more direct method for characterizing near-receiver structure than time-domain models of scattered waves [e.g., Sato, 1984]. If the resonances of the structure beneath one's receivers are known, we can hope to determine better the spectral shape of the original seismic source. If we model the response R (f , z(hat) ) of the crustal structure local to a receiver to waves traveling in the lithospheric wave guide with frequency f and polarization z(hat), we expect observed three-component amplitude spectra U (f ) to be found by integrating

U (f ) = int from Ω ~ R (f, z(hat) ) s(f, z(hat) ) d Ω

where s (f , z(hat) ) is the amplitude of the impinging signal. We integrate z(hat) over the lower half of the unit sphere in order to account for energy arriving from all vertical azimuths and out of plane scattering. In the example of Sereno and Orcutt [1985], R (f , z(hat) ) was calculated for a simple layered model. For arrays (such as Anza) positioned atop a heterogeneous medium, constraints on R (f , z(hat) ) can be found empirically using a number of events at different azimuths. Determination of R (f , z(hat) ) may be helpful in evaluating the earthquake hazards of a potential building site, especially as polarization analysis specifies both seismic amplitude and particle motion at the recording site. More research is necessary to determine if such a project is feasible. The above examples suggest that s (f , z(hat) ) varies significantly within the coda, complicating the determination of the near-receiver resonant structure.

5. Conclusions

We have devised a multitaper algorithm to determine the polarization of particle motion as a function of frequency and applied it to data recorded on the Anza Seismic Telemetered Array. We form a matrix of eigenspectra of three-component records and perform a singular value decomposition to estimate the complex-valued unit vector z(hat) whose components specify the sense of particle motion in the plane defined by the two real vectors Re ( z(hat) ) and Im ( z(hat) ). We manipulate the components of z(hat) in order to specify four angles. The angle φ_HH represents the relative phase between the components of horizontal motion. The angle φ_HH = 0° or ± 180° if the particle motion is rectilinear in the horizontal plane, and φ_HH ~ ± 90° if the motion is elliptical and oriented along the component axes. The phase angle φ_VH is the relative phase between horizontal and vertical motion. The apparent azimuth Θ_H is defined by the maximum displacement of the horizontal projection of the particle motion ellipse. It is measured in the counterclockwise direction from the first horizontal component. Finally, an angle of incidence Θ_V of the particle motion is estimated. The uncertainties in these polarization angles can be estimated from the singular value decomposition used to obtain z(hat) (appendix).

The variability of the spectra and polarization over 0 ≤ f ≤ 30 Hz suggest that the P coda observations can be separated into several distinct varieties of seismic motions, each occupying a separate frequency band. This behavior suggests that in some cases it may be more appropriate to model the P wave coda as a set of resonant modes caused by near-receiver structure rather than a number of randomly scattered compressional pulses. Evidence for scattered energy is not lacking, however, as the principal polarization accounts for a smaller proportion of the total seismic energy late in the P coda, accounting for only 60-65% in some frequency bands. We also observe that the apparent P wave arrival azimuth can vary by up to 50°, both between adjacent frequency bands and in adjacent time windows. Both rectilinear and elliptically polarized signals are found, often coexisting in the same time window in adjacent frequency bands. We find that the apparent modal structure of the signal polarization breaks down if the length of the time window is much greater than 2 s, suggesting incoherent excitation by direct and scattered seismic waves.

We are currently investigating the polarization behavior of the data recorded at each site in the Anza array. We want to use the polarization information to obtain better estimates of the seismic source spectrum. Such an endeavor requires that one be able to identify the factors causing the apparent jumps in polarization, both as a function of frequency and time.

Acknowledgments. This chapter has appeared in the Journal of Geophysical Research, vol. 92, pp. 12,664-12,674, 1987 under the title "Frequency dependent polarization analysis of high-frequency seismograms," by J. Park, F. L. Vernon III, and C. R. Lindberg. Frank Vernon participated in all phases of the development of this research except for the programming of the computer algorithms.

Appendix: Formal Uncertainty of Polarization Estimates

We estimate the uncertainties in the angles Θ_V, Θ_H, φ_VH, φ_HH from uncertainties δ z(hat) (f ) in the unit eigenvector z(hat) (f ), which represents the principal polarization of particle motion at frequency f. The derivation of the rms expectation of δ z(hat) can be found in the work by Park and Chave [1984]. We only define the problem and state the results here. The vector z(hat) = v(hat)₁, the right eigenvector of M (defined in (6.3)) associated with largest singular value d₁. The uncertainty σ is estimated from the two smaller singular values

σ ² = K over K-1 ~ (d₂² + d ₃ ² ) / 2

where K is the number of eigenspectra used in forming A. The covariance matrix for the first-order uncertainty δ z(hat) has expectation value

langle δ z(hat) citimes ( δ z(hat) ) ^* rangle = {sigma ²} over {Kd ₁} ~ sum from j=2 to 3 ~ bold v hat_jcitimes ( bold v hat_j) ^*

It is also true that l angle δ z(hat) citimes δ z(hat) rangle = 0. The $citimes$ symbol denotes the tensor (outer) product of two vectors. Since

langle δ z(hat) citimes ( δ z(hat) ) ^* rangle _ij = langle δ z sub i δ z_j^* rangle

we have complete information on the formal uncertainties of the components of the principal polarization. Note that since ( δ z(hat) ) ^* cdot z(hat) = 0 as z(hat) is a vector of unit length, δ z(hat) is composed of v(hat)₂ and v(hat)₃, the right eigenvectors associated with the "noise" singular values d₂, d₃.

Given (A2), we can determine the formal first-order uncertainty of any well-behaved function of z(hat) = (z₁, z₂, z₃ ); φ_HH = φ₂- φ₃ = arg (z₂) - arg (z₃) to within an additive constant, Since

δ φ_j= {i( δ z_j^* z_j- δ z_jz_j^* )} over {2|z_j| ²} = Im left ( {δ z} over z right )

langle | δ φ _HH | ² rangle = half left [ {langle | δ z₂| ² rangle} over {|z₂| ²} ~-~ 2 Re left ( {langle δ z₂δ z ₃ ^* rangle } over {z₂z ₃ ^*} right ) + {langle | δ z ₃ | ² rangle} over {|z ₃ | ²} right ]

Since Θ_H = Re ( tan^-1 (z₃/z₂) ),

δ Θ _H = Re left [ {z₂²} over {z₂² + z ₃ ²} ~cdot~ δ left ( {z ₃} over {z₂} right ) right ] = Re left [ {z₂δ z ₃ - z ₃ δ z₂} over {z₂² + z ₃ ²} right ]

so that

langle | δ Θ _H | ² rangle = half | z₂² + z ₃ ² | sup -2 (|z ₃ | ² langle | δ z₂| ² rangle

- 2 Re [ z₂^* z ₃ langle δ z₂δ z ₃ ^* rangle ] + |z₂| ² langle | δ z ₃ | ² rangle )

Note that l angle | δ Θ _H | ² rangle -> ∞ as z₂-> iz₃, i.e., circular polarization.

The uncertainties of the vertical polarization angles Θ_V and θ_VH are similar. With Θ_V given by (6.8), where z^~2 = z₂² + z₃², we use the relation δ z ~ = z ~ ^-1 (z₂δ z₂+ z ₃ δ z ₃ ) to find

δ Θ _V = Re left [ z ~ ^-1 left ( {z ~ ² δ z ₁ - z ₁ z₂δ z₂- z ₁ z ₃ δ z ₃} over {z ₁ ² + z₂² + z ₃ ²} right ) right ]

and

langle | δ Θ _V | ² rangle = half left [ {X ( z(hat) , δ z(hat) )} over {|z₂² + z ₃ ² | | z₂² + z ₃ ² + z ₁ ² | ²} right ]

where

X( z(hat) , δ z(hat) ) ~ mark =~ | z ~ | sup 4 langle | δ z ₁ | ² rangle + | z ₁ z₂| ² langle | δ z₂| ² rangle + | z ₁ z ₃ | ² langle | δ z ₃ | ² rangle

- 2 Re ( z ~ ² z ₁ ^* z₂^* langle δ z ₁ δ z₂^* rangle ) - 2 Re ( z ~ ² z ₁ ^* z ₃ ^* langle δ z ₁ δ z ₃ ^* rangle )

+~ 2 Re (| z ₁ | ² z₂z ₃ ^* langle δ z₂δ z ₃ ^* rangle )

The restriction of the argument of the arctangent to be positive in the definition of Θ_V does not alter its uncertainty. Following (A3) and (A4), the uncertainty of φ_VH is

δ φ_VH = Im left ( {δ z ~} over {z ~} ~-~ {δ z ₁} over {z₁} right )

and

langle | δ φ _VH | ² rangle = half left [ {X( z(hat) , δ z(hat) )} over {|z ₁ | ² | z₂² + z ₃ ² | ²} right ]

where X is given in (A9).

References

Archambeau, C. B., and E. A. Flinn, Automated analysis of seismic radiations for source characteristics, Proc. IEEE, 53, 1876-1884, 1965.

Bard, P., and M. Bouchon, The two-dimensional resonance of sediment-filled valleys, Bull. Seismol. Soc. Am., 75, 519-542, 1985.

Berger, J., L. M. Baker, J. N. Brune, J. B. Fletcher, T. C. Hanks, and F. L. Vernon III, The Anza array: A high-dynamic range, broadband, digitally radio telemetered seismic array, Bull. Seismol. Soc. Am., 74, 1469-1481, 1984.

Golub, G. H., and C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, Md., 1983.

Kanasewich, E. R., Time Sequence Analysis in Geophysics, 3rd ed., University of Alberta Press, Edmonton, Canada, 1981.

Kennett, B. L. N., Wavenumber and wavetype coupling in laterally heterogeneous media, Geophys. J. R. Astron. Soc., 87, 313-332, 1986.

Lanzerotti, L. J., D. J. Thomson, C. G. Maclennan, and L. V. Medford, Intermediate bandwidth polarization analysis using multiple-window methods, Eos Trans. AGU, 67 (44), 873, 1986.

Lindberg, C. R., Multiple taper spectral analysis of terrestrial free oscillations, Ph.D. thesis, 182 pp., Univ. of Calif., San Diego, La Jolla, 1986.

Odom, R. I., A coupled-mode examination of irregular waveguides including the continuum spectrum, Geophys. J. R. Astron. Soc., 86, 425-454, 1986.

Park, J., and A. D. Chave, On the estimation of magnetotelluric response functions using the singular value decomposition, Geophys. J. R. Astron. Soc., 77, 683-709, 1984.

Park, J. J., C. R. Lindberg, and D. J. Thomson, Multiple taper spectral analysis of terrestrial free oscillations: Part I, Geophys. J. R. Astron. Soc., 91, 755-794, 1987.

Park, J., F. L. Vernon III, and C. R. Lindberg, Multitaper spectral analysis of high-frequency seismograms, J. Geophys. Res., 92, 12,675-12,684, 1987.

Samson, J. C., Matrix and Stokes vector representations of detectors for polarized waveforms: Theory, with some applications to teleseismic waves, Geophys. J. R. Astron. Soc., 51, 583-603, 1977.

Samson, J. C., Pure states, polarized waves, and principal components in the spectra of multiple, geophysical time-series, Geophys. J. R. Astron. Soc., 72, 647-664, 1983a.

Samson, J. C., The spectral matrix, eigenvalues and principal components in the analysis of multichannel geophysical data, Ann. Geophysicae, 1 (2), 115-119, 1983b.

Samson, J. C., The reduction of sample-bias in polarization estimators for multichannel geophysical data with anisotroπc noise, Geophys. J. R. Astron. Soc., 75, 289-308, 1983c.

Sato, H., Attenuation and envelope formulation of three-component seismograms of small local earthquakes in randomly inhomogeneous lithosphere, J. Geophys. Res., 89, 1221-1241, 1984.

Sereno, T. J., and J. A. Orcutt, Synthetic seismogram modelling of the oceanic P_n phase, Nature, 316, 246-249, 1985.

Sleπan, D., Some comments on Fourier analysis, uncertainty and modelling, SIAM Rev., 25, 379-393, 1983.

Thomson, D. J., Spectrum estimation and harmonic analysis, IEEE Proc., 70, 1055-1096, 1982.

Vidale, J. E., Complex polarization analysis of particle motion, Bull. Seismol. Soc. Am., 76, 1393-1405, 1986.