## SUMMARY

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.

## INTRODUCTION

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:

*r*is radial distance,

*c*

_{R}is Rayleigh wave velocity,

*c*

_{L}is Love wave velocity,

*R*is the ratio of the horizontal-to-vertical motion of the Rayleigh waves,

*P*

^{R}is the power spectrum of the Rayleigh waves,

*P*

^{L}is the power spectrum of the Love waves and

*J*

_{0},

*J*

_{1}and

*J*

_{2}are Bessel functions of the zeroth, first and second orders, respectively. Note that the correlation matrix for the Love waves lacks non-zero off-diagonal components in this case. Thus, Aki's theory for Love waves is complete for isotropic incidence. Moving beyond an isotropic wavefield, we also derive the Rayleigh and Love wave correlation coefficient matrices for generally anisotropic incidence (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 *u*_{x} and vertical *u*_{z} displacements associated with the Rayleigh wave satisfy

*c*is the wave speed,

*R*is the ratio of the horizontal-to-vertical motion and ρ

_{n}= 2π

*n*/

*L*is the wavenumber associated with the

*n*th lateral mode over a finite interval of length

*L*. Eq. (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

*u*

_{x}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

*u*

_{x}is different from

*u*

_{z}in the following ways:

- (i)
It generally has a different amplitude, represented by the ellipticity

*R*. - (ii)
It is phase-shifted by 90° in time.

- (iii)
The sense of its motion changes whether the Rayleigh wave is propagating to the left or right.

Since *u*_{x} and *u*_{z} are real, the factors *A*_{n} and *B*_{n} 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):

*B*

_{n}. The symbol δ is the Kronecker delta function and the horizontal bar denotes averaging. Since they are uncorrelated random variables, the different lateral modes are taken to be in a state of white noise.

We first consider the spatial correlation of the vertical components, which was previously analysed by Aki (1957):

Based on eq. (5) above, the spatial correlation is given by where we have used the property to cancel the two terms with dependence on time as sin (*c*ρ

_{n}

*t*)cos (

*c*ρ

_{n}

*t*). Note that this is different in form than the analogous property invoked in eq. (6) of Aki (1957), which was simply . In subsequent uses of these quantities, we omit the horizontal bar but the averaging is assumed. Thus, Re (

*A*

_{n}

*B*

_{−n}) = Re (

*A*

_{−n}

*B*

_{n}) = 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 *A*_{n}*B*_{−n} and *A*_{−n}*B*_{n} 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.

As pointed out by Aki (1957), the time dependence in eq. (7) can be removed if we assume that

This is the equipartition condition, with*G*being a deterministic variable. The equipartition condition physically means that the different modes contribute an equal amount of energy (Snieder 2010). In this way, eq. (7) loses its dependence on time and becomes As shown by Aki (1957), we can transition from the discrete wavenumber ρ

_{n}to a continuous variable ρ and change the summation to an integral The crux of Aki's method is that we can map the wavenumber ρ integral to a frequency ω integral by a change of variables ρ =ω/

*c*Aki (1957) further argues that the term |

*G*(ω/

*c*)|

^{2}dρ/dω is equal to the spectral power Φ(ω). As a result For a narrowband signal centred at ω = ±ω

_{0}, Φ (ω) =π

*P*(ω

_{0}) δ (ω −ω

_{0}) +π

*P*(−ω

_{0}) δ (ω +ω

_{0}). If

*P*(ω

_{0}) =

*P*(−ω

_{0}), the correlation coefficient can then finally be written as This equation shows that in 1-D the correlation coefficient is given by the cosine function (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,

*R*

^{2}.

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 correlation

Similar to eq. (7), we write the expression for the correlation coefficient by exploiting the fact that the different modes are uncorrelated random variables where now the equipartition condition in eq. (9) has removed terms with dependence on time as sin (*c*ρ

_{n}

*t*) cos (

*c*ρ

_{n}

*t*). We previously discussed that the condition

*A*

_{n}

*B*

_{−n}= −

*A*

_{−n}

*B*

_{n}removed the terms with dependence on time as sin (

*c*ρ

_{n}

*t*) cos (

*c*ρ

_{n}

*t*) from the expression for φ

_{ZZ}. Here, we use the relation

*A*

_{n}

*B*

_{−n}= −

*A*

_{−n}

*B*

_{n}to remove the time-dependence from φ

_{ZX}. We arrive at For the term containing the product of

*A*

_{n}

*B*

_{−n}in the above equation, we substitute the following relationship This relationship agrees with the equipartition condition, eq. (9), in that which is the same as as it should be. Furthermore, the condition in eq. (18) reflects the fact that the terms

*A*

_{n}

*B*

_{−n}and

*A*

_{−n}

*B*

_{n}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.

Substituting eq. (18) into the summation in eq. (17) results in a Riemann sum in terms of the discrete wavenumber

which, as discussed before, can be approximated in the limit as an integral Changing the variable of integration to frequency instead of wavenumber gives Making use of the result by Aki (1957) that the term |*G*(ω/

*c*)|

^{2}dρ/dω is equal to the spectral power Φ(ω) yields where we set sgn (ω) = sgn (ω/

*c*) since the phase velocity

*c*is always positive. For a narrowband signal centred at ω = ±ω

_{0}, Φ(ω) =π

*P*(ω

_{0})δ(ω −ω

_{0}) +π

*P*(−ω

_{0})δ(ω +ω

_{0}) and, when

*P*(ω

_{0}) =

*P*(−ω

_{0}), the correlation coefficient becomes This equation shows that in 1-D the correlation coefficient for the ZX correlation is given by the sine function. An important aspect of this derivation is the presence of the absolute value for the frequency |ω

_{0}| in 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

*A*

_{n}and

*B*

_{n}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 *J*_{0}. By analogy between cosine and *J*_{0}, we can expect the ZR correlation to be dependent on the first order Bessel function *J*_{1}, since *J*_{1} 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 is

where*r*

_{1}and

*r*

_{2}are values of the horizontal and vertical components of the Rayleigh wave eigenfunction at the surface,

*U*is the group velocity and

*I*

_{1}is a integral over depth of density weighted by the eigenfunctions at depth. Note the presence of the Hilbert transform operator in the off-diagonal components. For our purposes, we take the normalization as , which is slightly different than the normalization used by Snieder (2002), 1/(8

*cUI*

_{1}) = 1. With this normalization,

*r*

_{1}/

*r*

_{2}is equal to the ratio of the horizontal-to-vertical motion,

*R*. Thus the Green's matrix becomes Comparing eq. (26) to eq. (28), Nakahara (2006) found the relationship between the correlation coefficient and the Green's function in 1-D to be for the unpolarized ZZ component. Note that, for this relation to hold, the normalization

*P*(ω) = 1 is used in 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 that

For negative frequencies ω < 0, we also find such that, in both cases, φ_{ZX}= (2ω/

*c*) × Im (

*G*

_{ZX}) = −

*R*sin (|ω|

*x*/

*c*) and the expression for φ

_{ZX}in eq. (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

*G*

_{ij}can be written as We now turn to expressions for the correlation coefficient matrix in the time domain. As shown by Nakahara (2006), the expression for φ

_{ZZ}in the time domain follows as: where

*F*

^{−1}denotes the inverse Fourier transform. In the context of Rayleigh waves, this result is applicable to a homogeneous half-space since

*c*generally varies as a function of frequency. However, the result is applicable to Rayleigh waves within a narrow frequency band defined by the power spectrum

*P*(ω

_{0}), within which

*c*does not vary greatly. For the φ

_{ZX}component, we need to compute where similarly to φ

_{ZZ}, we have assumed

*R*and

*c*are frequency independent. We proceed by breaking the integral in eq. (34) into two integrals over positive and negative frequencies We then recombine the integrals into two integrals over the entire frequency range From these expressions, it is clear that we are taking the inverse Fourier transform of the Hilbert transform operator −isgn (ω). This transform is known as the Hilbert transform

*H*of the delta function,

*H*[δ(

*t*)] = −1/π

*t*. We refer to this function in the shorthand notation

*H*[δ(

*t*)] =δ

^{H}(

*t*). As a result, the φ

_{ZX}component in the time domain can be written as Whereas 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 be

*P*(ω

_{0}) ≠

*P*(−ω

_{0}). If we set

*P*(ω

_{0}) = 2

*A*and

*P*(−ω

_{0}) = 2

*B*and transform to the time domain, we find, in agreement with Nakahara (2006), that In other words, the delta functions at negative and positive time-lags do not have the same amplitude. The presence of anisotropic incident energy in 1-D does not cause there to be acausal ‘ghost’ waveforms, as exist for anisotropic incident energy in 2-D (Nakahara 2006). By analogy with φ

_{ZZ}, we find that the general form for φ

_{ZX}(

*x*, ω

_{0}) in 1-D is given by Again setting

*P*(ω

_{0}) = 2

*A*and

*P*(−ω

_{0}) = 2

*B*and transforming to the time domain, we find the general form to be Just as seen with φ

_{ZZ}, the right- and left-going waveforms in this case are not equal in amplitude. Before leaving the 1-D case, note that the difference between the expressions in eqs (13) and (24) is the presence of the term

*Ri*sgn(ω) 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 satisfy

where and ψ = tan^{−1}(

*y*/

*x*) are the radial distance and azimuthal angle describing the relative positions of the two seismometers. In contrast to the 1-D case, the displacements are now the result of two summations over the radius ρ

_{n}and polar angle θ

_{m}for the wavenumber. Another difference with the 1-D case is that the factor sgn(ρ

_{n}) does not appear in eqs (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.

We first consider the spatial correlation of the vertical components (ZZ) before moving on to the mixed-component ZR correlation

Similar to the 1-D case, we find that The equipartition condition in 2-D relates the*A*

_{nm}and

*B*

_{nm}coefficients as Inserting this into eq. (46), we obtain the following expression As shown previously in 1-D, we can transition from the discrete wavenumber ρ

_{n}and discrete polar angle θ

_{m}to the continuous variables ρ and θ. This transition changes the summation to an integral From Cox (1973), the exponential in this equation can be expressed as a summation over cylindrical harmonics Inserting this infinite series into eq. (49) gives We further assume that

*G*can be split into separate terms dependent on ρ and θ as |

*G*(ρ, θ)|

^{2}=|

*G*(ρ)|

^{2}

*p*(θ). Inserting this into the above expression results in where the angular variation determines the following factor γ

_{m}As in the 1-D case, we implement the change of variables ρ =ω/

*c*to obtain As shown by Aki (1957), we find that the spectral power is given by and upon substituting this expression into eq. (54), we find that Finally, we take the wavefield to be narrowband by assuming Φ(ω) =π

*P*(ω

_{0}) δ (ω −ω

_{0}), which yields Since

*p*(θ) is real, and the two-sided infinite series can be rewritten as where ε

_{0}= 1 and ε

_{m}= 2 for

*m*≠ 0. For an isotropic wavefield, one in which energy propagates equally in all directions,

*p*(θ) = 1, γ

_{0}= 1 and γ

_{m}= 0 for

*m*≠ 0. In this case, which is the well-known result of Aki (1957) for unpolarized waves.

Using the methodology demonstrated above for the ZZ correlation, we now consider the ZR correlation given by

The derivation in this case follows similar logic as the above derivation for φ_{ZZ}and therefore we only point out the differences for φ

_{ZR}. Due to the term cos (θ −ψ) in eq. (43), the factor γ

_{m}is different from the φ

_{ZZ}case and is given by As pointed out in the 1-D case, φ

_{ZR}is different from φ

_{ZZ}in that the term

*R*isgn(ω) acts as an additional frequency filter for φ

_{ZR}. This results in the following expression for φ

_{ZR}Again, we take Φ(ω) =π

*P*(ω

_{0}) δ (ω −ω

_{0}), resulting in Since

*p*(θ) is real, and the two-sided infinite series can be rewritten as where as before, ε

_{0}= 1 and ε

_{m}= 2 for

*m*≠ 0.

For an isotropic wavefield, *p*(θ) = 1, γ_{1} = 1/2 and γ_{m} = 0 for *m*≠ 1. In this case,

*J*

_{m}(

*x*) is an odd function of

*x*if

*m*is an odd number, with

*J*

_{m}(

*x*) an even function otherwise. This result shows that the correlation function for the ZR component is given by the first-order Bessel function

*J*

_{1}in the case of an isotropic wavefield as opposed to the zeroth-order Bessel function for the ZZ component. Although results for SPAC are often given as normalized correlation functions, with the normalization at

*r*= 0 (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 by

and Note that Aki (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*J*

_{2}(

*x*) = −

*J*

_{0}(

*x*) for large

*x*. Thus, in the far field, φ

_{TT}= 0 and φ

_{RR}is proportional to φ

_{ZZ}by a factor of

*R*

^{2}, that is φ

_{RR}=

*R*

^{2}

*P*(ω)

*J*

_{0}(ω

*r*/

*c*

_{R}). 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.

As shown previously, we can summarize the results for Rayleigh waves in 2-D in the presence of isotropic incidence with the complete correlation matrix

where*P*

^{R}(ω) is the power spectrum of the Rayleigh waves. As shown in 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 is

where*r*

_{1}and

*r*

_{2}are values of the radial and vertical components of the Rayleigh wave eigenfunction at the surface,

*U*is the group velocity and

*I*

_{1}is a integral over depth of density weighted by the eigenfunctions at depth. Note that the Green's matrix in Aki & 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,

*r*

_{1}/

*r*

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

*H*denotes the Hilbert transform operator . As in the 1-D case, the normalization

*P*(ω) = 1 is used for φ

_{ZZ}in this expression. We point out that there is a difference in sign between 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 kind

For large values of ω*r*/

*c*, this can be expressed in the far-field as the sum of exponentials For ω > 0, we use this expression for

*m*= 1 and find that For the Green's function, we take the following expression from the matrix in eq. (70) for ω > 0 to find that 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 holds for all components of the correlation matrix and the Green's function matrix.

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

*T*

_{m}(

*ct*/

*r*) and

*U*

_{m − 1}(

*ct*/

*r*) are the

*m*th order Chebyshev polynomials of the first kind and (

*m*− 1)th order Chebyshev polynomials of the second kind, respectively. Nakahara (2006) showed that for a restricted angular range Δ of incident energy in a cone centred about azimuth φ

_{0}, the coefficient γ

_{m}is given by This result is also given in 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

- (i)
Hilbert transformed compared to the ZZ correlation.

- (ii)
Multiplied by the ratio of the radial-to-vertical motion

*R*. - (iii)
Subject to a different γ

_{m}coefficient given by eq. (61) instead of eq. (53).

Taking these differences into account, we find that the time-domain expressions for the ZR correlation and its Hilbert transform are given for |*t*| < *r*/*c* as

_{m}derived from eq. (B2) in Appendix B and given as In Fig. 2, we evaluate and plot the theoretical expressions for the ZZ and ZR correlations and their Hilbert transforms in eqs (78), (79), (81) and (82) by taking terms in the summations up to

*m*= 10 000. The plots are shown for a noise azimuth of φ

_{0}= 75° and four types of directional noise: isotropic incidence (Δ = 180°), incidence from a half-plane (Δ = 90°), a cone of angles given by Δ = 30° and a narrow cone of Δ = 5°. Comparison of the ZZ and ZR correlations shows that the ZR correlation is less susceptible than ZZ to acausal ghost arrivals that develop for Δ = 30° and Δ = 5°. The resiliency to ghosts for ZR correlations is due to the γ

_{m}coefficient in eq. (61). The polarization of the waves in the ZR case acts as a spatial filter that dampens the out-of-plane ghost arrivals. This makes the ZR correlations more desirable when considering relatively short time windows for the correlations, when the condition of isotropic incidence has a smaller likelihood of being satisfied. In practice, the length of the time window in time-lapse noise correlation studies involves a tradeoff between convergence to the ideal Green's function and the desire to have adequate time resolution. For example, in Haney (2009), correlations of infrasound noise averaged over 10 hr were able to detect day-long variations in sound speed due to changing atmospheric conditions. The robustness demonstrated in Fig. 2 means the ZR correlations are a better candidate for use in detecting time-lapse changes with seismic noise than ZZ correlations. Fig. 3 shows the case when the waves propagate perpendicular to the line of receivers, when φ

_{0}= 90°. The dampening effect of the acausal ghost arrivals is observed to be even more powerful for ZR compared to ZZ. Note that in the presence of both Rayleigh and Love wave noise, the ZZ and ZR correlations are unaffected by the Love waves, even in the near field when the two stations are within a wavelength of each other. Therefore, these results hold for general ambient noise conditions, when the noise is comprised of both Rayleigh and Love waves.

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

*Z*,

*R*,

*T*) coordinate system can be done on-the-fly after computing the correlations. Also note that, had the correlations not been done in such a way as to balance the three-component amplitudes, the rotation to the (

*Z*,

*R*,

*T*) coordinate system would not have been accurate.

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 *k*th zero

*r*= 8.7 km is the interstation distance and

*z*

_{k}is the

*k*th zero of the zeroth or first order Bessel function, depending on whether ZZ or ZR is being analysed. It is clear from Fig. 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

*p*= −

*b*/

*m*.

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 *J*_{0} for ZZ or *J*_{1} 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 *J*_{1} 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.

## CONCLUSION

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.

### ACKNOWLEDGMENTS

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.

## REFERENCES

### Appendices

#### APPENDIX A: GENERAL FORM FOR ANISOTROPIC INCIDENCE

For Rayleigh waves, the expression for the correlation function for generally anisotropic incidence is the following:

where*P*

^{R}(ω

_{0}) is the power spectrum of the Rayleigh waves, the indices

*i*,

*j*take on values of either Z, R or T,

*M*

_{ij}is a multiplier given by and the angular integrals are given by

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 *M*_{RR}, *M*_{RT}, *M*_{TR} and *M*_{TT} since the radial and transverse motions are not elliptically polarized, as is the case for the vertical and either of the horizontal motions.

For Love waves, the vertical displacement is zero and the horizontal displacements are

and Applying the analysis to the RR, RT, TR and TT correlations gives the following general form where*P*

^{L}(ω

_{0}) is the power spectrum of the Love waves, the indices

*i*,

*j*take on values of either R or T, and the angular integrals are given by

#### 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 ψ.

#### APPENDIX C: CONNECTION TO GREEN's FUNCTION FOR ω < 0 IN 2-D

For the Green's function, we take the following expression from the matrix in eq. (70) to find that Eqs (C1) and (C3) can be made equal to each other if eq. (C1) is multiplied by a factor of*i*. This factor is the Hilbert transform operator −isgn (ω) for the case we are currently considering, ω < 0.