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

## 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, 1983*a,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* [1983*a*]. 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 ^{1} (t), x ^{2} (t), x ^{3} (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 ^{1}(t)* is the vertical component. If the

*j*th record

*x*has the frequency domain representation

^{j}(t)*z*, the spectral density matrix

^{j}(f)**S**

*(f)*has components

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

where *E* denotes the expectation operator. *Samson* [1983*a*] 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_{n}x^{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 2*W* [*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 _{0} ^{(1)} (f) above y _{1} ^{(1)} (f) above 3dot above y sub K-1 ^{(1)} (f) } ccol { y _{0} ^{(2)}
(f) above y _{1} ^{(2)} (f) above 3dot above y sub K-1 ^{(2)} (f) } ccol { y _{0} ^{(3)} (f) above y _{1} ^{(3)} (f) above 3dot above y sub K-1 ^{(3)}
(f) } } ~ right ]

The value of *K*, the number of eigenspectra used, depends on 2*W*, the width of the frequency band in which the spectral energy at frequency *f* is concentrated. The *K*
= 2*NW* τ - 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* x*K* 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*, 1983*b*]. 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_{1} ≥ d_{2}≥ d_{3}. If d_{1} >> d_{2}, d_{3}, the polarization **z**(hat)
= **v**(hat)_{1} 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_{1} ~ d_{2}>> d_{3},
there is a strong possibility that coherent seismic motion at *f*_{0} exists at two separate polarizations. The dot product **v**(hat)_{1}* cdot **v**(hat)_{2}=
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_{1} >> d_{2}, d_{3}, 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)_{1} cdot **z**(hat) ) **e**(hat)_{1}, where **e**(hat) _{1} = (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} ) | ^{2}

If the components (z_{1} , z_{2}, z_{3} ) of **z**(hat) are expressed in the form *z _{j}*= |

*z*| e

_{j}^{i&φ;j}, this is equivalent to finding the maxima of

|z_{2}| ^{2} cos ^{2} (2 π f t + φ_{2}) + |z _{3} | ^{2} cos ^{2} (2 π f t + φ _{3} )

The extremes of this expression, remembering

| z_{2}| ^{2} sin 2 φ_{2} + | z _{3} | ^{2} sin 2 φ _{3} = Im [ z_{2}^{2} + z _{3} ^{2} ]

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

θ _{l} = - half roman arg [ z_{2}^{2} + z _{3} ^{2} ] + { 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_{1} ) < 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)_{2}=
(0,1,0) can be defined as

Θ _{H} = tan ^{-1} left [ {Re (z _{3} e sup {-i θ _{H}} )} over {Re (z_{2}e sup {-i θ _{H}} )} right ] = Re ( tan ^{-1} (z _{3}
/ z_{2}) )

The range of the arctangent function is 0 < Θ _{H} ≤ 180 if Re (z _{1} z _{3}^{*} ) < 0 and -180 < Θ _{H} ≤ 0 if Re (z _{1}
z _{3}^{*} ) ≥ 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 φ _{3} - φ_{2}= φ _{HH}, the phase difference between the horizontal components of particle motion. If φ_{2}- φ_{3}
~ 0 or 180, the particle motion is predominantly linear. The value φ_{2}- φ_{3} ~ 90 represents elliptical motion with the major and minor axes oriented along the axes of the instruments.
If z_{2}= ± iz_{3}, 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_{2}^{2} + z_{3}^{2} | ^{-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_{1}^{2} + z_{H}^{2} ] + { m π } over 2

where *m* is an integer and z_{H}^{2} = z_{2}^{2} + z_{3}^{2}. 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 _{1} 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} - φ_{1}. 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 ^{1} (n τ ) above x ^{2} (n τ ) above x ^{3} (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 φ _{V}H} , 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 3*b* 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 3*a*. The angles Θ_{H}
and φ_{HH} are plotted in Figures 3*c* and 3*d*. 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 3*e*. 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 5

*a*and 5

*b*. 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 _{2}* and

*d*displayed in Figure 6

_{3}*b*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 6

*a*). 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 6

*e*, 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 _{1}* 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

_{2}*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 7

_{2}*c*-

*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 7

*e*.

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 3*f*. 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 _{2}* 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 Θ

_{3}_{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 12

*e*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)_{1}, the right eigenvector
of **M** (defined in (6.3)) associated with largest singular value *d _{1}*. The uncertainty σ is estimated from the two smaller singular values

σ ^{2} = *K* over *K*-1 ~ (d_{2}^{2} + d _{3} ^{2} ) / 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 ^{2}} over {Kd _{1}} ~ sum from j=2 to 3 ~ bold v hat_{j}citimes (
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)_{2} and **v**(hat)_{3}, the right eigenvectors
associated with the "noise" singular values *d _{2}, d_{3}*.

Given (A2), we can determine the formal first-order uncertainty of any well-behaved function of **z**(hat) = (z_{1}, z_{2}, z_{3} ); φ_{HH} = φ_{2}-
φ_{3} = arg (z_{2}) - arg (z_{3}) to within an additive constant, Since

δ φ_{j}= {i( δ z_{j}^{*} z_{j}- δ z_{j}z_{j}^{*} )} over {2|z_{j}| ^{2}} = Im left ( {δ z} over z
right )

langle | δ φ _{H}H | ^{2} rangle = half left [ {langle | δ z_{2}| ^{2} rangle} over {|z_{2}| ^{2}} ~-~ 2 Re left ( {langle δ z_{2}δ
z _{3} ^{*} rangle } over {z_{2}z _{3} ^{*}} right ) + {langle | δ z _{3} | ^{2} rangle} over {|z _{3} | ^{2}} right ]

Since Θ_{H} = Re ( tan^{-1} (z_{3}/z_{2}) ),

δ Θ _{H} = Re left [ {z_{2}^{2}} over {z_{2}^{2} + z _{3} ^{2}} ~cdot~ δ left ( {z _{3}} over {z_{2}} right )
right ] = Re left [ {z_{2}δ z _{3} - z _{3} δ z_{2}} over {z_{2}^{2} + z _{3} ^{2}} right ]

so that

langle | δ Θ _{H} | ^{2} rangle = half | z_{2}^{2} + z _{3} ^{2} | sup -2 (|z _{3} | ^{2} langle | δ z_{2}| ^{2}
rangle

- 2 Re [ z_{2}^{*} z _{3} langle δ z_{2}δ z _{3} ^{*} rangle ] + |z_{2}| ^{2} langle | δ z _{3} | ^{2}
rangle )

Note that l angle | δ Θ _{H} | ^{2} rangle -> ∞ as z_{2}-> iz_{3}, 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_{2}^{2}
+ z_{3}^{2}, we use the relation δ z ~ = z ~ ^{-1} (z_{2}δ z_{2}+ z _{3} δ z _{3} ) to find

δ Θ _{V} = Re left [ z ~ ^{-1} left ( {z ~ ^{2} δ z _{1} - z _{1} z_{2}δ z_{2}- z _{1} z _{3} δ z
_{3}} over {z _{1} ^{2} + z_{2}^{2} + z _{3} ^{2}} right ) right ]

and

langle | δ Θ _{V} | ^{2} rangle = half left [ {X ( **z**(hat) , δ **z**(hat) )} over {|z_{2}^{2} + z _{3} ^{2}
| | z_{2}^{2} + z _{3} ^{2} + z _{1} ^{2} | ^{2}} right ]

where

X( **z**(hat) , δ **z**(hat) ) ~ mark =~ | z ~ | sup 4 langle | δ z _{1} | ^{2} rangle + | z _{1} z_{2}| ^{2} langle | δ
z_{2}| ^{2} rangle + | z _{1} z _{3} | ^{2} langle | δ z _{3} | ^{2} rangle

- 2 Re ( z ~ ^{2} z _{1} ^{*} z_{2}^{*} langle δ z _{1} δ z_{2}^{*} rangle ) - 2 Re ( z ~ ^{2} z _{1} ^{*}
z _{3} ^{*} langle δ z _{1} δ z _{3} ^{*} rangle )

+~ 2 Re (| z _{1} | ^{2} z_{2}z _{3} ^{*} langle δ z_{2}δ z _{3} ^{*} 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 _{1}} over {z_{1}} right )

and

langle | δ φ _{V}H | ^{2} rangle = half left [ {X( **z**(hat) , δ **z**(hat) )} over {|z _{1} | ^{2} | z_{2}^{2}
+ z _{3} ^{2} | ^{2}} 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, 1983*a*.

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

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, 1983*c*.

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.