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:  

(1)
formula
 
(2)
formula
where ω is angular frequency, r is radial distance, cR is Rayleigh wave velocity, cL is Love wave velocity, R is the ratio of the horizontal-to-vertical motion of the Rayleigh waves, PR is the power spectrum of the Rayleigh waves, PL is the power spectrum of the Love waves and J0, J1 and J2 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 ux and vertical uz displacements associated with the Rayleigh wave satisfy  

(3)
formula
 
(4)
formula
where 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 nth 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 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:

  • (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 ux and uz are real, the factors An and Bn are the complex conjugates of An and Bn. 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):  

(5)
formula
with an analogous relation for Bn. 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):  

(6)
formula
Based on eq. (5) above, the spatial correlation is given by  
(7)
formula
where we have used the property  
(8)
formula
to cancel the two terms with dependence on time as sin (cρnt)cos (cρnt). Note that this is different in form than the analogous property invoked in eq. (6) of Aki (1957), which was simply graphic. In subsequent uses of these quantities, we omit the horizontal bar but the averaging is assumed. Thus, Re (AnBn) = Re (AnBn) = 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 graphic 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 AnBn and AnBn 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  

(9)
formula
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  
(10)
formula
As shown by Aki (1957), we can transition from the discrete wavenumber ρn to a continuous variable ρ and change the summation to an integral  
(11)
formula
The crux of Aki's method is that we can map the wavenumber ρ integral to a frequency ω integral by a change of variables ρ =ω/c 
(12)
formula
Aki (1957) further argues that the term |G(ω/c)|2 dρ/dω is equal to the spectral power Φ(ω). As a result  
(13)
formula
For a narrowband signal centred at ω = ±ω0, Φ (ω) =πP0) δ (ω −ω0) +πP (−ω0) δ (ω +ω0). If P0) =P (−ω0), the correlation coefficient can then finally be written as  
(14)
formula
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, 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 correlation  

(15)
formula
Similar to eq. (7), we write the expression for the correlation coefficient by exploiting the fact that the different modes are uncorrelated random variables  
(16)
formula
where now the equipartition condition in eq. (9) has removed terms with dependence on time as sin (cρnt) cos (cρnt). We previously discussed that the condition AnBn = −AnBn removed the terms with dependence on time as sin (cρnt) cos (cρnt) from the expression for φZZ. Here, we use the relation AnBn = −AnBn to remove the time-dependence from φZX. We arrive at  
(17)
formula
For the term containing the product of AnBn in the above equation, we substitute the following relationship  
(18)
formula
This relationship agrees with the equipartition condition, eq. (9), in that  
(19)
formula
which is the same as  
(20)
formula
as it should be. Furthermore, the condition in eq. (18) reflects the fact that the terms AnBn and AnBn 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  

(21)
formula
which, as discussed before, can be approximated in the limit as an integral  
(22)
formula
Changing the variable of integration to frequency instead of wavenumber gives  
(23)
formula
Making use of the result by Aki (1957) that the term |G(ω/c)|2 dρ/dω is equal to the spectral power Φ(ω) yields  
(24)
formula
where we set sgn (ω) = sgn (ω/c) since the phase velocity c is always positive. For a narrowband signal centred at ω = ±ω0, Φ(ω) =πP0)δ(ω −ω0) +πP (−ω0)δ(ω +ω0) and, when P0) =P (−ω0), the correlation coefficient becomes  
(25)
formula
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 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:  

(26)
formula
In the next section, we discuss the relationship between this matrix and the Green's function matrix for Rayleigh waves in 1-D.

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  

(27)
formula
where r1 and r2 are values of the horizontal and vertical components of the Rayleigh wave eigenfunction at the surface, U is the group velocity and I1 is a integral over depth of density weighted by the eigenfunctions at depth. Note the presence of the Hilbert transform operator graphic in the off-diagonal components. For our purposes, we take the normalization as graphic, which is slightly different than the normalization used by Snieder (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  
(28)
formula
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  
(29)
formula
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  

(30)
formula
For negative frequencies ω < 0, we also find  
(31)
formula
such that, in both cases, φZX = (2ω/c) × Im (GZX) = −Rsin (|ω|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 Gij can be written as  
(32)
formula
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:  
(33)
formula
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 P0), within which c does not vary greatly. For the φZX component, we need to compute  
(34)
formula
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  
(35)
formula
We then recombine the integrals into two integrals over the entire frequency range  
(36)
formula
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  
(37)
formula
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  

(38)
formula
where P0) ≠P (−ω0). If we set P0) = 2A and P (−ω0) = 2B and transform to the time domain, we find, in agreement with Nakahara (2006), that  
(39)
formula
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  
(40)
formula
Again setting P0) = 2A and P (−ω0) = 2B and transforming to the time domain, we find the general form to be  
(41)
formula
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 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 satisfy  

(42)
formula
 
(43)
formula
 
(44)
formula
where graphic 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  

(45)
formula
Similar to the 1-D case, we find that  
(46)
formula
The equipartition condition in 2-D relates the Anm and Bnm coefficients as  
(47)
formula
Inserting this into eq. (46), we obtain the following expression  
(48)
formula
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  
(49)
formula
From Cox (1973), the exponential in this equation can be expressed as a summation over cylindrical harmonics  
(50)
formula
Inserting this infinite series into eq. (49) gives  
(51)
formula
We further assume that G can be split into separate terms dependent on ρ and θ as |G(ρ, θ)|2 =|G(ρ)|2p (θ). Inserting this into the above expression results in  
(52)
formula
where the angular variation determines the following factor γm 
(53)
formula
As in the 1-D case, we implement the change of variables ρ =ω/c to obtain  
(54)
formula
As shown by Aki (1957), we find that the spectral power is given by  
(55)
formula
and upon substituting this expression into eq. (54), we find that  
(56)
formula
Finally, we take the wavefield to be narrowband by assuming Φ(ω) =πP0) δ (ω −ω0), which yields  
(57)
formula
Since p(θ) is real, graphic and the two-sided infinite series can be rewritten as  
(58)
formula
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,  
(59)
formula
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  

(60)
formula
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  
(61)
formula
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 
(62)
formula
Again, we take Φ(ω) =πP0) δ (ω −ω0), resulting in  
(63)
formula
Since p(θ) is real, graphic and the two-sided infinite series can be rewritten as  
(64)
formula
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,  

(65)
formula
where the absolute value follows from the fact that Jm(x) is an odd function of x if m is an odd number, with Jm(x) an even function otherwise. This result shows that the correlation function for the ZR component is given by the first-order Bessel function J1 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 P0).

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  

(66)
formula
and  
(67)
formula
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 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 (ω) J0r/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.

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  

(68)
formula
where PR(ω) 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  

(69)
formula
where r1 and r2 are values of the radial and vertical components of the Rayleigh wave eigenfunction at the surface, U is the group velocity and I1 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 graphic. With this normalization, r1/r2 is equal to the ratio of the horizontal-to-vertical motion, R. Thus the Green's matrix becomes  
(70)
formula
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  
(71)
formula
where H denotes the Hilbert transform operator graphic. 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  

(72)
formula
For large values of ωr/c, this can be expressed in the far-field as the sum of exponentials  
(73)
formula
For ω > 0, we use this expression for m = 1 and find that  
(74)
formula
For the Green's function, we take the following expression from the matrix in eq. (70) for ω > 0  
(75)
formula
to find that  
(76)
formula
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  
(77)
formula
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  

(78)
formula
and  
(79)
formula
where we take the azimuth between the stations to be 0°. In these expressions, Tm(ct/r) and Um − 1(ct/r) are the mth 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  
(80)
formula
This result is also given in eq. (B1) in Appendix B along with other angular integrals for incident energy within a range of angles.

Figure 1

Angles used to demonstrate the effect of anisotropic incidence on correlations.

Figure 1

Angles used to demonstrate the effect of anisotropic incidence on correlations.

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  

(81)
formula
and  
(82)
formula
with the coefficient γm derived from eq. (B2) in Appendix B and given as  
(83)
formula
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.

Figure 2

Comparison of ghosts in ZZ and ZR correlations for φ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°. The correlations φZZ and φZR are given as dashed lines, with their Hilbert transforms shown as solid lines. Time has been normalized by distance and phase velocity as in Nakahara (2006). A normalized time of unity is therefore equal to the arrival time of the direct wave between the stations.

Figure 2

Comparison of ghosts in ZZ and ZR correlations for φ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°. The correlations φZZ and φZR are given as dashed lines, with their Hilbert transforms shown as solid lines. Time has been normalized by distance and phase velocity as in Nakahara (2006). A normalized time of unity is therefore equal to the arrival time of the direct wave between the stations.

Figure 3

Comparison of ghosts in ZZ and ZR correlations for φ0 = 90° and 4 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°. The correlations φZZ and φZR are given as dashed lines, with their Hilbert transforms shown as solid lines. Time has been normalized by distance and phase velocity as in Nakahara (2006). A normalized time of unity is therefore equal to the arrival time of the direct wave between the stations.

Figure 3

Comparison of ghosts in ZZ and ZR correlations for φ0 = 90° and 4 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°. The correlations φZZ and φZR are given as dashed lines, with their Hilbert transforms shown as solid lines. Time has been normalized by distance and phase velocity as in Nakahara (2006). A normalized time of unity is therefore equal to the arrival time of the direct wave between the stations.

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  

(84)
formula
where θ is the azimuth (in degrees counter-clockwise from east) between the stations. The importance of this type of processing is that the angular rotation can be performed after the cross-correlation step. This means that for an array of stations, the rotations into the (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.

Figure 4

A regional map showing the location of Akutan Volcano (upper panel) and a local map of Akutan Island (lower panel). The locations of the two stations used to compute cross correlations are given on the local map. Bold contours are plotted at 0 and 850 m elevation; light contours are shown between 150 and 675 m elevation in intervals of 175 m. The highest contour line, at 850 m, clearly marks the location of the summit caldera at Akutan Volcano.

Figure 4

A regional map showing the location of Akutan Volcano (upper panel) and a local map of Akutan Island (lower panel). The locations of the two stations used to compute cross correlations are given on the local map. Bold contours are plotted at 0 and 850 m elevation; light contours are shown between 150 and 675 m elevation in intervals of 175 m. The highest contour line, at 850 m, clearly marks the location of the summit caldera at Akutan Volcano.

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 zero  

(85)
formula
where r= 8.7 km is the interstation distance and zk is the kth 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.

Figure 5

The real part of the frequency spectrum for the average time-domain correlation between stations AKMO and AKLV at Akutan Volcano for ZZ (top panel) and ZR (bottom panel). The spectrum is tapered to zero below 0.1 Hz. The zeros of the real part of the frequency spectrum are shown for both ZZ and ZR by vertical lines. Rayleigh wave phase velocities, measured at the locations of the zero crossings of ZZ and ZR, are shown in Fig. 6.

Figure 5

The real part of the frequency spectrum for the average time-domain correlation between stations AKMO and AKLV at Akutan Volcano for ZZ (top panel) and ZR (bottom panel). The spectrum is tapered to zero below 0.1 Hz. The zeros of the real part of the frequency spectrum are shown for both ZZ and ZR by vertical lines. Rayleigh wave phase velocities, measured at the locations of the zero crossings of ZZ and ZR, are shown in Fig. 6.

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  

(86)
formula
where the derivatives are computed according to the relation p = −b/m.

Figure 6

Rayleigh wave phase velocities at Akutan Volcano measured from the first four ZZ (red squares) and ZR (red circles) zero-crossings. The phase velocities are obtained by associating the ZZ zero-crossings with J0 and the ZR zero-crossings with J1. Uncertainties in the phase velocities are also plotted as described in the text. The blue curve is the Rayleigh wave phase velocity dispersion curve for a 1-D model at nearby Okmok Volcano (Masterlark 2010). The Okmok dispersion curve shows that the phase velocities obtained at Akutan volcano are consistent with structure inferred at a nearby volcano in the Aleutian Arc.

Figure 6

Rayleigh wave phase velocities at Akutan Volcano measured from the first four ZZ (red squares) and ZR (red circles) zero-crossings. The phase velocities are obtained by associating the ZZ zero-crossings with J0 and the ZR zero-crossings with J1. Uncertainties in the phase velocities are also plotted as described in the text. The blue curve is the Rayleigh wave phase velocity dispersion curve for a 1-D model at nearby Okmok Volcano (Masterlark 2010). The Okmok dispersion curve shows that the phase velocities obtained at Akutan volcano are consistent with structure inferred at a nearby volcano in the Aleutian Arc.

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.

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

1
Aki
K.
,
1957
.
Space and time spectra of stationary stochastic waves, with special reference to microtremors
,
Bull. Earthq. Res. Inst. Tokyo Univ.
 ,
25
,
415
457
.
2
Aki
K.
,
1965
.
A note on the use of microseisms in determining the shallow structures of the Earth's crust
,
Geophysics
 ,
30
,
665
666
.
3
Aki
K.
Richards
P.
,
1980
. ,
W. H. Freeman and Company
,
San Francisco, CA
.
4
Aki
K.
Richards
P.
,
2002
. ,
W. H. Freeman and Company
,
San Francisco
.
5
Campillo
M.
Paul
A.
,
2003
.
Long-range correlations in the diffuse seismic coda
,
Science
 ,
299
,
547
549
.
6
Chávez-García
F.J.
Luzón
F.
,
2005
.
On the correlation of seismic microtremors
,
J. geophys. Res.
 ,
110
,
B11313
, doi:10.1029/2005JB003671.
7
Chávez-García
F.J.
Rodríguez
M.
Stephenson
W.R.
,
2005
.
An alternative approach to the SPAC analysis of microtremors: exploiting stationarity of noise
,
Bull. seism. Soc. Am.
 ,
95
,
277
293
.
8
Cho
I.
Tada
T.
,
2006
.
A generic formulation for microtremor exploration methods using three-component records from a circular array
,
Geophys. J. Int.
 ,
165
,
236
258
.
9
Chouet
B.
Tilling
R.I.
,
1996
.
New methods and future trends in seismological volcano monitoring
, in , pp.
Springer
,
New York, NY
.
23
97
.
10
Chouet
B.
De Luca
G.
Milana
G.
Dawson
P.
Martini
M.
Scarpa
R.
,
1998
.
Shallow velocity structure of Stromboli Volcano, Italy, derived from small-aperture array measurements of strombolian tremor
,
Bull. seism. Soc. Am.
 ,
88
,
653
666
.
11
Chouet
B.
Dawson
P.
Arciniega-Ceballos
A.
,
2005
.
Source mechanism of Vulcanian degassing at Popocatépetl Volcano, Mexico, determined from waveform inversion of very long period signals
,
J. geophys. Res.
 ,
110
,
B07301
, doi:10.1029/2004JB003524.
12
Cox
H.
,
1973
.
Spatial correlation in arbitrary noise fields with application to ambient sea noise
,
J. acoust. Soc. Am.
 ,
54
,
1289
1301
.
13
Dixon
J.P.
Stihler
S.D.
Power
J.A.
,
2008
. ,
U.S. Geol. Surv. Data Ser., No. 367, U.S. Geological Survey
,
Reston, VA
.
14
Dorman
J.
Ewing
M.
,
1962
.
Numerical inversion of seismic surface wave dispersion data and crust-mantle structure in the New York-Pennsylvannia area
,
J. geophys. Res.
 ,
67
,
5227
5241
.
15
Durand
S.
Montagner
J.P.
Roux
P.
Brenguier
F.
Nadeau
R.M.
Ricard
Y.
,
2011
.
Passive monitoring of anisotropy change associated with the Parkfield 2004 earthquake
,
Geophys. Res. Lett.
 ,
38
,
L13303
, doi:10.1029/2011GL047875.
16
Ekström
G.
Nettles
M.
Abers
G.A.
,
2003
.
Glacial earthquakes
,
Science
 ,
302
(
5645
),
622
624
.
17
Ekström
G.
Abers
G.A.
Webb
S.C.
,
2009
.
Determination of surface-wave phase velocities across US Array from noise and Aki's spectral formulation
,
Geophys. Res. Lett.
 ,
36
, doi:10.1029/2009GL039131.
18
Ferrazzini
V.
Aki
K.
Chouet
B.
,
1991
.
Characteristics of seismic waves composing Hawaiian volcanic tremor and gas-piston events observed by a near-source array
,
J. geophys. Res.
 ,
96
,
6199
6209
.
19
Haney
M.M.
,
2009
.
Infrasonic ambient noise interferometry from correlations of microbaroms
,
Geophys. Res. Lett.
 ,
36
,
L19808
, doi:10.1029/2009GL040179.
20
Haney
M.M.
,
2010
.
Location and mechanism of very long period tremor during the 2008 eruption of Okmok Volcano from interstation arrival times
,
J. geophys. Res.
 ,
115
,
B00B05
, doi:10.1029/2010JB007440.
21
Haney
M.M.
Nies
A.
Masterlark
T.
Needy
S.
Pedersen
R.
,
2011
.
Interpretation of Rayleigh-wave ellipticity observed with multicomponent passive seismic interferometry at Hekla Volcano, Iceland
,
Leading Edge
 ,
30
,
526
531
.
22
Harmon
N.
Rychert
C.
Gerstoft
P.
,
2010
.
Distribution of noise sources for seismic interferometry
,
Geophys. J. Int.
 ,
183
,
1470
1484
.
23
Lin
F.-C.
Moschetti
M.P.
Ritzwoller
M.H.
,
2008
.
Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps
,
Geophys. J. Int.
 ,
173
,
281
298
.
24
Lobkis
O.I.
Weaver
R.L.
,
2001
.
On the emergence of the Greens function in the correlations of a diffuse field
,
J. acoust. Soc. Am.
 ,
110
,
3011
3017
.
25
Lu
Z.
Wicks
C.
Jr
Power
J.A.
Dzurisin
D.
,
2000
.
Ground deformation associated with the March 1996 earthquake swarm at Akutan volcano, Alaska, revealed by satellite radar interferometry
,
J. geophys. Res.
 ,
105
(
B9
),
21 483
21 495
.
26
Masterlark
T.
Haney
M.
Dickinson
H.
Fournier
T.
Searcy
C.
,
2010
.
Rheologic and structural controls on the deformation of Okmok Volcano, Alaska: FEMs, InSAR, and ambient noise tomography
,
J. geophys. Res.
 ,
115
,
B02409
, doi:10.1029/2009JB006324.
27
Muyzert
E.
,
2007
a
Seabed property estimation from ambient-noise recordings. Part 1: compliance and Scholte wave phase-velocity measurements
,
Geophysics
 ,
72
,
U21
U26
.
28
Muyzert
E.
,
2007
b
Seabed property estimation from ambient-noise recordings. Part 2: Scholte-wave spectral-ratio inversion
,
Geophysics
 ,
72
,
U47
U53
.
29
Nakahara
H.
,
2006
.
A systematic study of theoretical relations between spatial correlation and Green's function in one-, two-, and three-dimensional random scalar wavefields
,
Geophys. J. Int.
 ,
167
,
1097
1105
.
30
Okada
H.
,
2003
.
The Microtremor Survey Method
 ,
Geophys. Monogr. No. 12, Society of Exploration Geophysicists
,
Tulsa , OK
.
31
Prieto
G.A.
Lawrence
J.F.
Beroza
G.C.
,
2009
.
Anelastic Earth structure from the coherency of the ambient seismic field
,
J. geophys. Res.
 ,
114
,
B07303
, doi:10.1029/2008JB006067.
32
Rhie
J.
Romanowicz
B.
,
2004
.
Excitation of Earth's continuous free oscillations by atmosphere-ocean-seafloor coupling
,
Nature
 ,
431
(
7008
),
552
556
.
33
Roux
P.
,
2009
.
Passive seismic imaging with directive ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA
,
Geophys. J. Int.
 ,
179
,
367
373
.
34
Roux
P.
Wathelet
M.
Roueff
A.
,
2011
.
The San Andreas Fault revisited through seismic-noise and surface-wave tomography
,
Geophys. Res. Lett.
 ,
38
,
L13319
, doi:10.1029/2011GL047811.
35
Saccorotti
G.
Chouet
B.
Dawson
P.
,
2003
.
Shallow-velocity models at the Kilauea Volcano, Hawaii, determined from array analyses of tremor wavefields
,
Geophys. J. Int.
 ,
152
,
633
648
.
36
Snieder
R.
,
2001
. ,
Cambridge University Press
,
Cambridge
.
37
Snieder
R.
,
2002
.
Scattering of surface waves
, in , pp.
Academic Press
,
San Diego, CA
.
562
577
.
38
Snieder
R.
Fan
Y.
Slob
E.
Wapenaar
K.
,
2010
.
Equipartitioning is not sufficient for Green's function extraction
,
Earthq. Sci.
 ,
23
,
403
415
.
39
Stehly
L.
Campillo
M.
Shapiro
N.
,
2007
.
Travel time measurements from noise correlation: stability and detection of instrumental errors
,
Geophys. J. Int.
 ,
171
,
223
230
.
40
Tada
T.
Cho
I.
,
2007
.
Beyond the SPAC method: Exploiting the wealth of circular-array methods for microtremor exploration
,
Bull. seism. Soc. Am.
 ,
97
,
2080
2095
.
41
Tanimoto
T.
Alvizuri
C.
,
2006
.
Inversion of the HZ ratio of microseisms for S-wave velocity
,
Geophys. J. Int.
 ,
165
,
323
335
.
42
Teal
P.D.
Abhayapala
T.D.
Kennedy
R.A.
,
2002
.
Spatial correlation for general distributions of scatterers
,
IEEE Signal Process. Lett.
 ,
9
,
305
308
.
43
Tsai
V.C.
,
2009
.
On establishing the accuracy of noise tomography travel-time measurements in a realistic medium
,
Geophys. J. Int.
 ,
178
,
1555
1564
.
44
Tsai
V.C.
Moschetti
M.P.
,
2010
.
An explicit relationship between time-domain noise correlation and spatial autocorrelation (SPAC) results
,
Geophys. J. Int.
 ,
182
,
454
460
.
45
Van Wijk
K.
Mikesell
T.D.
Schulte-Pelkum
V.
Stachnik
J.
,
2011
.
Estimating the Rayleigh-wave impulse response between seismic stations with the cross terms of the Green tensor
,
Geophys. Res. Lett.
 ,
38
,
L16301
, doi:10.1029/2011GL047442.
46
Zobin
V.M.
Plascencia
I.
Reyes
G.
Navarro
C.
,
2009
.
The characteristics of seismic signals produced by lahars and pyroclastic flows: Volcán de Colima, México
,
J. Volc. Geotherm. Res.
 ,
179
(
1-2
),
157
167
.

Appendices

APPENDIX A: GENERAL FORM FOR ANISOTROPIC INCIDENCE

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

(A1)
formula
where PR0) is the power spectrum of the Rayleigh waves, the indices i, j take on values of either Z, R or T, Mij is a multiplier given by  
(A2)
formula
and the angular integrals are given by  
(A3)
formula

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.

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

(A4)
formula
and  
(A5)
formula
Applying the analysis to the RR, RT, TR and TT correlations gives the following general form  
(A6)
formula
where PL0) 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  
(A7)
formula

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

(B1)
formula
 
(B2)
formula
 
(B3)
formula
 
(B4)
formula
 
(B5)
formula
 
(B6)
formula

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

For ω < 0, we find that  

(C1)
formula
For the Green's function, we take the following expression from the matrix in eq. (70) 
(C2)
formula
to find that  
(C3)
formula
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.