Using ambient seismic noise for imaging subsurface structure dates back to the development of the spatial autocorrelation (SPAC) method in the 1950s. We present a theoretical analysis of the SPAC method for multicomponent recordings of surface waves to determine the complete 3 × 3 matrix of correlations between all pairs of three-component motions, called the correlation matrix. In the case of isotropic incidence, when either Rayleigh or Love waves arrive from all directions with equal power, the only non-zero off-diagonal terms in the matrix are the vertical–radial (ZR) and radial–vertical (RZ) correlations in the presence of Rayleigh waves. Such combinations were not considered in the development of the SPAC method. The method originally addressed the vertical–vertical (ZZ), RR and TT correlations, hence the name spatial autocorrelation. The theoretical expressions we derive for the ZR and RZ correlations offer additional ways to measure Rayleigh wave dispersion within the SPAC framework.
Expanding on the results for isotropic incidence, we derive the complete correlation matrix in the case of generally anisotropic incidence. We show that the ZR and RZ correlations have advantageous properties in the presence of an out-of-plane directional wavefield compared to ZZ and RR correlations. We apply the results for mixed-component correlations to a data set from Akutan Volcano, Alaska and find consistent estimates of Rayleigh wave phase velocity from ZR compared to ZZ correlations. This work together with the recently discovered connections between the SPAC method and time-domain correlations of ambient noise provide further insights into the retrieval of surface wave Green's functions from seismic noise.
Seismologists often encounter signals that do not exhibit distinct P- or S-wave phase arrivals. Examples of these signals include oceanic microseisms (Tanimoto & Alvizuri 2006; Muyzert 2007a,b), the Earth's hum (Rhie & Romanowicz 2004), glacial earthquakes (Ekström 2003), volcanic tremor (Ferrazzini 1991; Haney 2010), volcanic explosions (Chouet 2005) and debris flows such as lahars (Zobin 2009). For these cases, the seismologist's standard toolbox, based on picking the arrival times and amplitudes of body-wave phase arrivals, is not applicable. At best, for short duration signals such as regional earthquakes and teleseisms, passing surface waves can be identified and used to infer structure of the crust and upper mantle (Dorman & Ewing 1962). Recent advances in seismology, including ambient noise tomography (ANT), have pioneered the use of extended signals lacking body waves for the purpose of imaging Earth structure. ANT achieves this goal by correlating long-time recordings of ambient seismic noise within the frequency band dominated by the oceanic microseism (<1 Hz). Although the proliferation of ambient noise techniques has mostly occurred during the past decade, the idea of correlating long duration seismic signals from spatially extended sources goes back much further. In fact, the basic approach to interstation correlation originates over 50 yr ago, in seminal papers by Aki (1957, 1965). Aki's spatial autocorrelation (SPAC) method has subsequently been used to study shallow subsurface structure using wavefields composed of microseisms (Okada 2003) and volcanic tremor (Ferrazzini 1991; Chouet 1996; Chouet 1998; Saccorotti 2003). Okada (2003) provided a comprehensive review of the SPAC method, also known as the microtremor method.
Recently, Aki's work has been revisited in light of advances made in the use of ambient noise following the groundbreaking work of Campillo & Paul (2003). Connections between time-domain ambient noise correlations and the SPAC method have been explored by several authors (Chávez-García & Luzón 2005; Chávez-García 2005; Nakahara 2006; Ekström 2009; Prieto 2009; Tsai 2009; Harmon 2010; Tsai & Moschetti 2010). In the presence of an isotropic noise field, Nakahara (2006) connected the Fourier transform of the time-domain ambient noise correlation and the spectral correlation coefficient from the SPAC method. Building on this connection, Ekström (2009) showed that the zeros in the real part of the spectrum of averaged time-domain correlations could be associated with the spectral zeros predicted by Aki's SPAC theory. In the case of an anisotropic noise field, Nakahara (2006) derived expressions for time-domain ambient noise correlations drawing on results derived in the field of acoustics (Cox 1973; Teal 2002). Through these studies, the SPAC method and time-domain correlations have been shown to be different but largely equivalent ways of using ambient noise to gain information on the structure of the subsurface.
Despite the connections forged between the two methods, comparisons between SPAC and time-domain correlations have thus far been limited to discussions of correlations between vertical–vertical (ZZ) components, with some discussion of radial–radial (RR) and transverse–transverse (TT) components (Chávez-García & Luzón 2005). Cho & Tada (2006) and Tada & Cho (2007) have extended the SPAC method under general conditions for three-component recordings. However, mixed components of polarized waves, such as ZR and RZ, have not been considered before. This is due in part to the fact that Aki (1957) only discussed ZZ, RR and TT correlations in the original paper on the SPAC method. In fact, the SPAC method has its name in part because it considers correlations between identical components, and therefore, a comparison between the complete Green's tensor of time-domain correlations and the SPAC method cannot be made. Even so, several authors have presented the full matrix of possible time-domain ambient noise correlations from three-component instruments (Campillo & Paul 2003; Stehly 2007; Roux 2009; Durand 2011; Roux 2011). Recent work by Van Wijk (2011) on the mixed components in the time domain suggests these independent results provide additional information that may be more robust than the diagonal terms. In this work, we extend the SPAC method to the complete matrix of correlation coefficients and compare the result to time-domain correlations between three-component instruments. Furthermore, we quantify the sensitivity of the different components of the correlation matrix to out-of-plane, directional ambient noise. To do so, we proceed in the way shown by Nakahara (2006) for unpolarized waves recorded on a single vertical component.
For isotropic incidence, our results for the off-diagonal components of the Rayleigh φR and Love φL wave correlation coefficient matrices, combined with Aki's results for the main diagonal, can be summarized as follows:Appendix A). For anisotropic incidence, all off-diagonal components are non-zero in general for both Rayleigh and Love waves.
In the following, we present the theory for Rayleigh wave propagation in 1-D and 2-D and make connections with the time-domain Green's functions. We show the robustness of ZR correlations in the presence of out-of-plane ambient noise. With data from two broadband seismometers at Akutan Volcano, Alaska, we use ZR correlations to obtain Rayleigh wave dispersion curve measurements within the SPAC methodology that are complementary to measurements from ZZ correlations.
THE 1-D CASE FOR RAYLEIGH WAVES
We begin by analysing Rayleigh waves in a 1-D setting. Although the results we derive are not the ultimate findings of our study, the 1-D case has many conceptual similarities to the 2-D case. In this way, we can demonstrate the extensions we propose to the original SPAC method of Aki (1957) without some of the complexities of the 2-D case. Our derivation follows closely the development of SPAC in Aki (1957) and we adopt much of the notation used in that paper. We refer to the original manuscript for detailed discussions on the assumptions inherent to the method.
The 1-D case generally consists of Rayleigh waves propagating to the left or right within a single vertical plane. The horizontal ux and vertical uz displacements associated with the Rayleigh wave satisfyEq. (3) has the same form as the solution to the 1-D wave equation presented by Aki (1957). To analyse Rayleigh waves, we have added the coupled equation for ux in eq. (4). The Rayleigh waves considered here are assumed to consist of a single vertical mode, for example, the fundamental mode as in Aki (1957). Note that ux is different from uz in the following ways:
It generally has a different amplitude, represented by the ellipticity R.
It is phase-shifted by 90° in time.
The sense of its motion changes whether the Rayleigh wave is propagating to the left or right.
Since ux and uz are real, the factors An and Bn are the complex conjugates of A−n and B−n. Furthermore, as shown in eq. (4) of Aki (1957), the amplitudes of the different modes are assumed to be uncorrelated random variables (Lobkis & Weaver 2001):
We first consider the spatial correlation of the vertical components, which was previously analysed by Aki (1957):eq. (5) above, the spatial correlation is given by Aki (1957), which was simply . In subsequent uses of these quantities, we omit the horizontal bar but the averaging is assumed. Thus, Re (AnB−n) = Re (A−nBn) = 0 as in Aki (1957) but the imaginary parts of these terms are non-zero and opposite in sign.
Aki (1957) explains that the condition physically means the initial distributions of displacement and particle velocity are independent. It is worth noting that the alternative condition we adopt in this paper, eq. (8), also satisfies this physical requirement. The condition in eq. (8) is less restrictive than the condition from Aki (1957): it implies that the sum of the terms AnB−n and A−nBn is zero instead of each term being equal to zero individually. This distinction was not important for Aki (1957) since mixed-components were not analyzed, but will play a major role here.Snieder 2010). In this way, eq. (7) loses its dependence on time and becomes Aki (1957), we can transition from the discrete wavenumber ρn to a continuous variable ρ and change the summation to an integral Aki (1957) further argues that the term |G(ω/c)|2 dρ/dω is equal to the spectral power Φ(ω). As a result Aki 1957). This result also holds for the correlation coefficient of the horizontal components φXX, except that φXX is multiplied by a factor equal to the squared ellipticity, R2.
We now analyse the correlation coefficient that results from mixing vertical and horizontal components. These were not considered by Aki (1957) in his original paper. First we look at the ZX correlationeq. (7), we write the expression for the correlation coefficient by exploiting the fact that the different modes are uncorrelated random variables eq. (9) has removed terms with dependence on time as sin (cρnt) cos (cρnt). We previously discussed that the condition AnB−n = −A−nBn removed the terms with dependence on time as sin (cρnt) cos (cρnt) from the expression for φZZ. Here, we use the relation AnB−n = −A−nBn to remove the time-dependence from φZX. We arrive at eq. (9), in that eq. (18) reflects the fact that the terms AnB−n and A−nBn are purely imaginary. Finally, the sign convention we use in eq. (18) is consistent with the convention used to define Rayleigh wave ellipticity in Aki & Richards (1980, 2002), in which retrograde motion corresponds to a negative value for R. We further demonstrate this consistency in the following section. Aki (1957) that the term |G(ω/c)|2 dρ/dω is equal to the spectral power Φ(ω) yields eq. (25). As we will see, the absolute value ensures that the correlation function φZX is closely related to the ZX component of the Rayleigh wave Green's function in 1-D. Since the imaginary parts of the two terms with An and Bn in eq. (18) are opposite in sign, the XZ component of the correlation coefficient matrix φXZ is equal to the negative of φZX.
Eq. (25) foreshadows the 2-D case since we know from Aki (1957) that the ZZ correlation is the zero order Bessel function J0. By analogy between cosine and J0, we can expect the ZR correlation to be dependent on the first order Bessel function J1, since J1 is analogous with sine. In a later section, we will rigorously show that this is indeed the case. We can summarize our results for 1-D Rayleigh waves as follows:
RELATION TO THE GREEN's FUNCTION IN 1-D
As in Nakahara (2006), we now show that there is a relationship between the correlation matrix and the Green's function matrix. The Green's matrix for Rayleigh waves is given by eqs (7.108)–(7.110) in Aki & Richards (1980) or eqs (7.109)–(7.111) in Aki & Richards (2002). The matrix for the fundamental mode Rayleigh wave isSnieder (2002), 1/(8cUI1) = 1. With this normalization, r1/r2 is equal to the ratio of the horizontal-to-vertical motion, R. Thus the Green's matrix becomes eq. (26) to eq. (28), Nakahara (2006) found the relationship between the correlation coefficient and the Green's function in 1-D to be eq. (26). This relationship shows that the correlation coefficient is proportional to the time-derivative of the difference between the Green's function and the time-reversed Green's function. Since φZZ and φXX are proportional to each other, this relation holds for the XX component.
We now show that the relationship holds for the off-diagonal components as well. First, consider the ZX component for positive frequencies ω > 0. From eq. (28), we find thateq. (26) is verified. The expression for φXZ can similarly be verified and the connection between any element of the matrix of correlation coefficients φij and the Green's function matrix Gij can be written as Nakahara (2006), the expression for φZZ in the time domain follows as: eq. (34) into two integrals over positive and negative frequencies eq. (33) consisted of two propagating delta functions, eq. (37) indicates that there are two propagating waveforms which are the Hilbert transform of the delta function. Although the Hilbert transform of the delta function is an asymmetric wavelet, the two waveforms are subtracted and so the ZX correlation is symmetric with respect to t = 0. Thus, it is possible to average the positive and negative lags of the correlation in the time domain to obtain the so-called symmetric component (Lin 2008). The Hilbert transform of the delta function is the waveform since, relative to φZZ, the ZX component is phase shifted by 90°. This is a direct result of the elliptical polarization of the Rayleigh wave, as exploited in Van Wijk (2011).
The final issue we discuss for the 1-D case concerns the correlation coefficients when the incident energy is anisotropic. Although we did not point it out at the time, our assumption in the discussion after eq. (13), that the power spectrum was an even function of frequency (or, by ρ =ω/c, an even function of wavenumber), is equivalent to assuming that the incident energy is equal from both directions. When this is not the case, we find the general form for φZZ(x, ω0) to beNakahara (2006), that Nakahara 2006). By analogy with φZZ, we find that the general form for φZX(x, ω0) in 1-D is given by eqs (13) and (24) is the presence of the term Risgn(ω) within the integrand of eq. (24). This means that φZX is related to φZZ by multiplication by R and a Hilbert transform. We utilize this property in the 2-D case to simplify the derivation for the vertical–radial correlation φZR.
THE 2-D CASE FOR RAYLEIGH WAVES
We turn our attention to Rayleigh waves along a 2-D surface. The previous development of the 1-D case serves as a useful guide for the 2-D case. This was also the approach taken by Aki (1957) in the original derivation of SPAC. By analogy with eqs (3) and (4), the three-component displacements associated with Rayleigh waves in 2-D satisfyeqs (43) and (44) as it did in eq. (4). This is because the sign change due to different propagation directions is implicitly handled by the factors cos (θm −ψ) and sin (θm −ψ). Finally, note that the Rayleigh waves in eqs (42)–(44) do not necessarily propagate along the direct path between the receivers and, as a result, the horizontal particle motion is in general not strictly radial. eq. (46), we obtain the following expression Cox (1973), the exponential in this equation can be expressed as a summation over cylindrical harmonics eq. (49) gives Aki (1957), we find that the spectral power is given by eq. (54), we find that Aki (1957) for unpolarized waves. eq. (43), the factor γm is different from the φZZ case and is given by Aki 1957; Chouet 1996), here we can see that such a normalization is not possible for the ZR correlation since φZR is equal to 0 at r = 0. As a result, we have chosen to present correlation functions that are not normalized throughout this paper. The lack of normalization means that all of our expressions contain the spectral power at the reference frequency P (ω0).
As we saw previously for the 1-D case, it can be shown that the RZ correlation is the negative of the ZR correlation. Although the two correlations are equal to each other in magnitude, they do constitute independent measurable quantities since they utilize different components of the wavefield. Thus it is correct to consider both the ZR and RZ correlations as independent quantities. The RR and TT correlations were considered previously by Aki (1957) for an isotropic Rayleigh wavefield and shown to be given byAki (1957) used normalized correlation functions whereas these are not normalized. A surprising result from Aki (1957) is that the TT correlation in the presence of an isotropic Rayleigh wavefield is not zero. It is important to remember that asymptotically J2(x) = −J0(x) for large x. Thus, in the far field, φTT = 0 and φRR is proportional to φZZ by a factor of R2, that is φRR =R2P (ω) J0(ωr/cR). However, in the near field, φTT is non-zero. Due to the connections with the time-domain Green's function, as shown in the next section, this means that in the near-field of a transverse force source, Rayleigh waves are measured on a transverse receiver.
The remaining correlations involve combinations of Z or R components with the transverse component: ZT, TZ, RT and TR. For an isotropic wavefield, all of these terms are equal to zero. In Appendix A, we show the γm for these components in the case of generally anisotropic incidence. For isotropic incidence, the γm for these components are given in Appendix B by setting the angular range variable Δ = 180°. Following this procedure for the ZT, TZ, RT and TR components shows that they are indeed equal to zero for isotropic incidence.Appendix A, anisotropic incidence generally causes the off diagonal terms that are zero for isotropic incidence to become non-zero. Anisotropic incidence also causes the correlation matrix to become complex-valued. For Love waves, all off-diagonal terms are zero in the case of isotropic incidence, meaning that only RR and TT components exist since ZZ is equal to zero. Thus, the results of Aki (1957) for the main diagonal are complete for Love waves. Note that, in the presence of both Rayleigh and Love wave energy, the ZZ, ZR and RZ components are not affected by the Love waves. This is in contrast to the RR and TT coefficients, which are sensitive to both Rayleigh and Love waves, at least within the near-field as described above. In Appendix A, we give expressions for both Rayleigh and Love wave correlation matrices in the case of generally anisotropic incidence.
RELATION TO THE GREEN's FUNCTION IN 2-D
As in the 1-D case, we now show that there is a relationship between the 2-D correlation matrix and the Green's function matrix. The Green's matrix for Rayleigh waves is given in the far-field by eq. (7.146) in Aki & Richards (1980) or eq. (7.147) in Aki & Richards (2002). The matrix for the fundamental mode Rayleigh wave isAki & Richards (1980, 2002) is a 3 × 3 in terms of the three Cartesian coordinates x, y and z. Here, we use a 2 × 2 matrix in terms of the radial and vertical components, r and z, since the Rayleigh waves only exist on those components in the far-field. We again take the normalization as . With this normalization, r1/r2 is equal to the ratio of the horizontal-to-vertical motion, R. Thus the Green's matrix becomes Nakahara (2006) found a relationship between the correlation coefficient and the Green's function in 2-D. For the unpolarized ZZ correlation, the relation between the two is eq. (71) and the result of Nakahara (2006). This stems from the fact that there is a sign difference in the definition of the 2-D Green's function in eq. (18.23) of Snieder (2001), on which the result of Nakahara (2006) was based, and the Rayleigh wave Green's function in eq. (7.100) of Aki & Richards (1980). The sign difference arises because the forcing terms shown in these two references are on opposite sides of the 2-D wave equation. Note that there was no difference in sign for the previously considered 1-D case, eq. (29), and the corresponding result of Nakahara (2006). This is because, in addition to the sign difference, a sign error exists in the expression for the 1-D Green's function in eqs (18.55) and (18.56) of Snieder (2001).
In the following, we show that the relation in eq. (71) holds for the ZR component in 2-D. In contrast to the 1-D case, we rely on asymptotic approximations of Bessel functions to show this and therefore the demonstration applies to the far-field. This is necessary since the Green's matrix eq. (69) is derived in the far-field. From Snieder (2001), we know that the Bessel function of any order can be written as the sum of Hankel functions of the same order of both the first and second kindeq. (70) for ω > 0 Eqs (74) and (76) can be made equal to each other if eq. (74) is multiplied by a factor of −i. This factor is the Hilbert transform operator −isgn (ω) for the case we are currently considering, ω > 0. The equivalence we have just shown for ω > 0 is demonstrated for ω < 0 in Appendix C. Therefore, we have shown that
RESILIENCY TO GHOSTS
The ZR and RZ correlations have advantageous properties compared to the ZZ correlation in the presence of anisotropic ambient noise. We demonstrate these properties theoretically here, but this has also been shown observationally by Van Wijk (2011). To investigate the robustness of the ZR correlation, we proceed in the manner of Nakahara (2006) and study the theoretical time-domain correlation functions within a cone of times prior to the arrival of the direct wave. Since a relationship exists between correlation functions and the Green's function, there should ideally be no high-amplitude arrivals prior to the direct wave. As we show here, such arrivals—called ‘ghost’ arrivals—do occur when the noise wavefield is highly anisotropic in the out-of-plane direction between two seismometers. The ghost arrivals are more pronounced in ZZ correlations compared to ZR correlations and as a result, the ZR correlations can be considered to be more robust.
The situation we consider is shown in Fig. 1, where two stations are subjected to highly directional seismic noise centred about the azimuth φ0 and spanning an angular range of 2Δ. From Nakahara (2006), the time-domain expressions for the ZZ correlation and its Hilbert transform are given for |t| < r/c asNakahara (2006) showed that for a restricted angular range Δ of incident energy in a cone centred about azimuth φ0, the coefficient γm is given by eq. (B1) in Appendix B along with other angular integrals for incident energy within a range of angles.
For the expression for the ZR correlation in the time domain, we take advantage of the fact that the main differences with the ZZ correlation are that the ZR correlation is
Hilbert transformed compared to the ZZ correlation.
Multiplied by the ratio of the radial-to-vertical motion R.
SPAC WITH ZR CORRELATIONS AT AKUTAN VOLCANO
Recently, Ekström (2009) demonstrated the connection between the SPAC method and time-domain correlations of ambient noise between vertical components (ZZ). Namely, the zeros of the real part of the spectrum can be associated with the zeros of the zeroth-order Bessel function. Here, we show that our results allow the same methodology to be applied to ZR correlations as well, by associating the zeros of the real part of the spectrum with the zeros of the first-order Bessel function.
We use data from two broadband, three-component seismometers located at Akutan Volcano, historically one of the most active volcanoes in the Aleutian Island chain (Lu 2000). The two sensors, AKMO and AKLV, are both Guralp CMG-6TD instruments (Dixon 2008) that exist within the network operated by the Alaska Volcano Observatory. The location of Akutan Island and the seismic stations are given in Fig. 4. We compute time-domain correlations between all components (Z, N and E) over 40 d beginning on 2007 November 1. The correlations are done in a manner that maintains the relative amplitudes between the three-components associated with a single station (Haney 2011). This is accomplished by applying identical time-normalization and frequency-normalization filters to each of the three components. For the time-normalization, a moving average window of rms amplitude is computed at each time sample for each of the three components associated with a single station. We use a time-window of 10 s for the moving averages. A single time-normalization filter is formed from the maximum of the three moving-averages at each time sample. The three-component seismograms are then divided by the same value of the time-normalization filter at each time sample. Note that this balancing of amplitude in time normalization is impossible if the sign-bit type of normalization is applied to each of the three components individually. Hour-long time-normalized seismograms are then subjected to frequency normalization by forming smoothed versions of the amplitude spectra for each of the three components. The smoothing is done over a frequency range of 0.025 Hz. A frequency-normalization filter is formed from the maximum of the three smoothed amplitude spectra at each discrete frequency between 0.1 and 0.7 Hz. The spectra for each of the three components are then divided by the same frequency-normalization filter at each discrete frequency. Outside of this frequency band, the frequency spectrum of the signals are tapered to zero. Once the signals have been time and frequency normalized, all combinations of hour-long, three-component seismograms are correlated between the stations. Once we compute the correlation matrix between the stations in the (Z, N, E) coordinate system, we transform to (Z, R, T) by applying a tensor rotation
In Fig. 5, we plot the real part of the frequency spectrum of the time-domain correlations for both ZZ and ZR. It is clear that the zeros of ZZ and ZR do not coincide. This is expected since the zeros of the ZZ spectrum are associated with the zeros of the zeroth-order Bessel function, with the zeros of the ZR spectrum given by the first-order Bessel. Over the frequency band from 0.1 to 0.7, we pick four zeros in both the ZZ and ZR spectra, as shown in Fig. 5. We compute the phase velocities c at these zeros by setting the argument of the Bessel functions equal to their kth zeroFig. 5 that the lowermost frequency, 0.1 Hz, is low enough to capture the first zero crossing on the vertical component. For greater interstation distances, the lowermost frequency must be decreased in order to detect the first zero crossing. If seismic noise does not exist for lower frequencies, there can be some ambiguity in the association of zero-crossings with a particular zero.
We plot the computed phase velocities at Akutan in Fig. 6 along with the Rayleigh wave dispersion curve for a 1-D model of Okmok Volcano developed by Masterlark (2010). The phase velocities at Akutan are in good agreement with the Okmok model, showing that the values are generally consistent with those obtained at a nearby volcano in the Aleutian Arc. Uncertainty estimates are also plotted for the phase velocities in Fig. 6. The errors are obtained by fitting lines to the ZZ and ZR spectra in the immediate vicinity of the zero crossings. Close to the zero-crossings, the spectra should be almost linear, just as a sine function is linear near its zero-crossings. From the linear fit, we obtain error estimates (variances) for the slope and the y-intercept directly as well as the covariance between the slope and y-intercept estimates. Finally, the errors for the x-intercept (i.e. zero-crossing), indicated by the symbol p, are estimated from the variances and covariance for the slope m and y-intercept b by simple propagation of errors
In Fig. 6, the phase velocity estimations from both ZZ and ZR correlations are overall consistent with each other and show a decreasing phase velocity with increasing frequency. Note how the values of phase velocity from the ZR correlations naturally interpolate the values from the ZZ correlation. The agreement between the ZZ and ZR measurements below 0.4 Hz supports the theoretical results derived in this paper. Above 0.4 Hz, there is some discrepancy between the phase velocities estimated from the ZZ and ZR correlations as the estimates from the ZR correlations appear to be slightly larger. This is likely related to the wavelength of the Rayleigh waves at 0.4 Hz approaching the interstation distance. For a nominal phase velocity of 2.5 km s−1, the wavelength at 0.4 Hz is 6.25 km which is comparable to the interstation distance of 8.7 km. The assumption that the subsurface is an effective homogeneous medium between the stations becomes progressively less applicable as the frequency increases and the wavelength becomes smaller than the interstation distance. At higher frequencies, ray-bending and scattering may also become more important. Ray-bending, in particular, poses a problem since the ZR correlation assumes the azimuth and backazimuth of Rayleigh wave propagation between the stations are directed along the interstation azimuth. Another possible cause for the discrepancy at frequencies above 0.4 Hz is if the source distribution is not close to being isotropic and the zeros of the ZZ and ZR correlations are not determined simply from the first order term, either J0 for ZZ or J1 for ZR. These issues deserve further attention in future work.
Additional challenges for future work would involve fitting the observed spectra over the entire frequency band for a source power distribution and a generally anisotropic incident wavefield instead of only using the zeros of the spectra. Alternatively, the J1 dependence would be evident for the ZR correlation within a small-aperture array of several seismometers in which the spatial sampling was dense enough to observe the correlation coefficient as a function of interstation distance.
We have extended the SPAC method to the full matrix of correlation coefficients that exists between three-component seismometers. Applications of the theory include the use of ZR and RZ correlations within the SPAC methodology and the demonstration of increased robustness for ZR and RZ correlations in the presence of anisotropic noise distributions. A clear relationship exists between the time-domain correlations of ambient noise and the results of the SPAC method for the entire matrix of three-component correlation coefficients. Future topics to be explored include the extension of the 1-D and 2-D theories presented here to the 3-D case with P and S waves, inversion of the correlation spectra for the noise distribution, and the use of ZR and RZ correlations for imaging and monitoring the structure of the subsurface.
MH thanks Professor Paul Michaels of Boise State University for initially giving him a copy of the 1957 paper by Aki. This manuscript has benefitted from comments by Bernard Chouet, Lapo Boschi, and an anonymous reviewer. This work was partially funded under NSF award EAR-1142154 and USGS/ARRA award G10AC00016.
APPENDIX A: GENERAL FORM FOR ANISOTROPIC INCIDENCE
In the above expressions, R is the ratio of the horizontal-to-vertical motion. Also, note that the factor isgn (ω0) does not appear for MRR, MRT, MTR and MTT since the radial and transverse motions are not elliptically polarized, as is the case for the vertical and either of the horizontal motions.
APPENDIX B: ANGLE INTEGRALS OVER AN INTERVAL
For incident energy that is constant (p(θ) = 1) over an angular range of 2Δ centred on the angle φ0, the angular integrals are given by the following expressions. The azimuth between the seismometers is denoted by ψ.