Stationary phase approximation in the ambient noise method revisited

1.   Introduction

In recent years, it is shown that the inter-receiver Green's function could be extracted from cross correlation of ambient noise recorded at the receivers. This method was first reported in ultrasonic experiments (Weaver and Lobkis, 2001, 2002) and then has been widely applied to crustal seismology to get high resolution surface wave tomography (Shapiro and Campillo, 2004; Shapiro et al., 2005). The scales of studied regions range from regional (Sabra et al., 2005; Kang and Shin, 2006; Yao et al., 2006; Villasenor et al., 2007; Cho et al., 2007) to continental (Yang et al., 2006). However, there are still a lot of controversies about why cross correlation of ambient seismic noise leads to Green's function. Some scientists argued that either equipartition of modes or perfect diffusivity of wave field is necessary (Weaver and Lobkis, 2001, 2002; Shapiro and Campillo, 2004; Shapiro et al., 2005), which is indeed the case in ultrasonic experiments, but might not be realistic for the ambient seismic noise (Snieder, 2004). Instead, a derivation based on time-reversal invariance (Derode et al., 2003a, b) seems to be more physically attractive when the wavefield is not diffuse or the medium is infinite, though attenuation will break this time-reversal invariance (Snieder, 2007). In an effort to overcome the theoretical difficulties mentioned above, Snieder (2004) took a nice mathematical approach involving stationary phase approximation for homogeneous elastic medium. His derivation relaxes the global requirement of equipartitioning of normal modes to the local requirement that the scattered waves propagate isotropically near the receivers. Moreover, Snieder (2007) generalized this idea to inhomogeneous attenuating acoustic waves, which is much closer to the situation of seismic waves in a realistic world. In these two papers, he demonstrated that the emergence of Green's function is a direct consequence of constructive interference only of those scattered waves that propagate along the line connecting receivers. When the scatters are not near the receiver line, their contribution to cross correlation seems to disappear if the stationary phase approximation is assumed.

However, we know that for a homogeneous medium, the sources located on any hyperbola with its foci at the two receivers, not only the receiver line, produce stable phase difference at the two receivers. So we would like to investigate whether the stationary phase approximation has been applied appropriately, in order to understand the physics of emergence of Green's function from ambient noise.

We first examine Snieder's derivation by integrating the contribution from scatters distributed on a line perpendicular to the receiver line. Following that, we compute the cross correlation function in polar coordinates and compare the difference of phase and amplitude between cross correlation and Green's function. Finally, we discuss inaccuracy of ambient noise method due to the deviation of cross correlation from the real Green's function.

The rest of this paper is divided into three parts. In section 2, we will analyze and formulate the problems, with and without the stationary phase approximation. And then in section 3, we will present and discuss the deviation of the approximate solution given by Snieder (2004) from the exact solution from numerical computation. Section 4 is the conclusion.

2.   Numerical integral of NCF

At regional scale, numerous studies have shown that the Green's function extracted from ambient seismic noise is dominated by surface waves (e.g., Shapiro and Campillo, 2004; Shapiro et al., 2005; Sabra et al., 2005; Yao et al., 2006), though there are also some reports of body waves retrieved (e.g., Roux et al., 2005; Draganov et al., 2009; Zhan et al., 2010). In this paper, we only discuss the wave equation in a 2-D attenuating homogeneous medium. For brevity, the simplest case of scalar waves (acoustic waves) is chosen to illustrate our opinions. The geometry for our analysis is similar to that of Snieder (2004) as shown in Figure 1. The two receivers, A and B, are separated by a distance of R and located on the x axis. The origin of our coordinate system is chosen at receiver A. The dashed line L₀ in Figure 1 separates the left and the right half-plane. The scatterers act as secondary sources of singly and multiply scattered waves. Scatterer numbered s emits a signal S_s(t) that is due to all the waves that impinge upon that scatterer. Apart from the scatterers, the medium's phase velocity c and the quality factor Q are assumed to be spatially constant. In this case, the Green's function G(r, r₀) in the frequency domain, that is the displacement response at r to a force δ(r-r₀), is given by

Figure  1.  Definition of the geometric variables for our analysis. The black triangles denote the receivers A and B. The dashed line L₀ is the perpendicular bisector of the segment AB. We compute integral I(Y) numerically over the vertical dotted line L₁, x=X₀.

DownLoad: Full-Size Img PowerPoint

where k is the wavenumber, R is the distance between r and r₀, and L=Q/k is the attenuation length of the medium (Snieder, 2004). Equation (1) is a far-field approximation of the Hankel function in the 2-D surface wave Green's function. All expressions in this work are formulated in the frequency domain. The wave field at the two receivers A and B can be written as a superposition of the waves radiated by the scatterers:

where ω is the angular frequency and S_s(ω) is the Fourier transform of S_s(t), that is the frequency spectrum of scatterer s. Equation (2) implies that the waves emitted by scatterer s are isotropic, i.e., independent of directions. The summation is over all the scatterers on the plane. As shown by Snieder (2004), the scatterers distributed on the left half-plane (to the left of the dashed line L₀ in Figure 1) contribute to the causal Green's function G(r_B, r_A), while the scatterers on the right half-plane give the acausal Green's function G^*(r_B, r_A). Hereafter, only the scatterers on the left half-plane are considered and the summation is restricted to the left half-plane.

The cross correlation function between displacement records at the two receivers can be written as

Following Snieder (2004), the double sum ∑_{s, s′} could be split into a sum over diagonal terms ∑_s=s′ and a sum ∑_s≠s′ over cross terms. The sum of cross terms could be ignored by averaging over sufficient time/sources (Snieder, 2004). So

As there are many scatterers per wavelength, the resulting sum ∑_s(···) could be replaced by a surface integral

where n is the scatterer density per unit surface area and the ∑ stands for the left half-plane. The geometry of variables r₁ and r₂ is shown in Figure 1. Equation (5) is just the equation (9) in Snieder (2004), but with attenuation taken into account.

2.1   Validity of stationary phase approximation

There are different ways to solve the integral in equation (5). Snieder (2004) evaluated the integral over the transverse coordinate y in the stationary phase approximation, and then evaluated the simpler integral over the coordinate x. These steps lead to a simple but important result:

where is the average power spectrum of scatterers. This equation means that if we know we can extract the Green's function G(r_B, r_A) between A and B from the cross correlation function C(ω).

Obviously, the point y = 0 is the stationary point of the phase term (k(r₂ - r₁)) in equation (5) with the transverse coordinate y as variable for given x = x₀(L₁). This is the only stationary point that Snieder considered in his derivation. However, the phase term approaches zero at y = ±∞ and these stationary points should be taken into account. Though the contribution from ±∞ vanishes for attenuating media, we do not know the exact behavior of the competition between weaker contribution due to attenuation at longer propagation distances and stronger contribution due to constructive interference when y becomes very large. Therefore, we conduct numerical integration along the transverse line L₁ for different quality factor Q and frequency to investigate the behavior.

We calculate the integral (equation (5)) on the line L₁ with X₀ = R = 100 km, with Gaussian quadrature because of its higher efficiency. In Figure 2, we present the results for two Q factors (300 and 3 000) and three frequencies (0.2 Hz, 0.5 Hz, 2 Hz). The case Q = 3 000 describes an unrealistic weak attenuation for the Earth, while Q = 300 is a reasonable value for shield regions (Cong and Mitchell, 1998).

Figure  2.  Integral I(Y) defined by equation (7). In each panel, the solid lines and dashed lines are integrals I(Y) for numerical computation and stationary phase approximation, respectively (red for real and blue for imaginary).

DownLoad: Full-Size Img PowerPoint

In each panel of Figure 2, the solid lines and dashed lines are integrals I(Y) for numerical computation and stationary phase approximation, respectively (red for real and blue for imaginary). For Y close to zero, the two approaches match very well, suggesting that the stationary phase approximation works near the stationary point Y = 0. At high frequency (≥2 Hz), the match is even good for large Y. However, at lower frequencies (0.2 Hz and 0.5 Hz), when Y becomes larger and larger, curves from the two approaches are different from each other, so do their limits.

Based on Figure 2, we can conclude that it is inappropriate to evaluate the integral over the transverse coordinate y by the stationary phase approximation, either for realistic (Q=300) or low-attenuation media. So in the next section, we compute the integral numerically over the whole plane and discuss the errors incorporated by stationary phase approximation.

2.2   Numerical evaluation in polar coordinate

An alternative way is to compute the integral in equation (5) numerically, without using any approximation. As there is singularity at receiver A because of r₁ at the denominator, we need to transform equation (5) to from which we can see that all quantities in equation (5) that depend on x and y now depend on r₁ and θ, which are shown in Figure 1. The cutoff θ₀ is introduced to ensure that the integral is over the left half-plane.

For brevity, we assume that the power spectrum of ambient noise is spatially constant. In accordance with Snieder (2004), we also assume that the scatterer density n is constant in space. However, it is still difficult to deal with the infinite integral over r₁. Fortunately, the integrand's amplitude is controlled by exp[-(r₁+r₁)/2L], which attenuates exponentially as r₁ increases. So we can replace the upper limit of r₁ with a cutoff value, for example, 20L (L is the attenuation length), for which the amplitude is only 2×10^-9 of its initial amplitude. Numerical tests show that larger cutoff values do not change the results. Then equation (8) could be simplified as

Now the double integral is ready to be computed numerically. We partition the integration intervals of both r₁ and θ into many small enough subintervals, apply the 6-point Gauss-Legendre integration on each subinterval, and sum the results up to get the final C(ω).

3.   Comparison between cross correlation and Green's function

We assume the phase velocity c to be 3.5 km/s, quality factor Q to be 300. Receiver distance R is set to be 50 km, 100 km or 200 km, which are three typical station distances in ambient seismic noise studies. For each R, we evaluate the approximate and numerical cross correlation function C(ω) (denoted by C_ap(ω) and C_nu(ω) in the following, for brevity) by equation (6) and equation (10), respectively, for frequency from 0.01 Hz to 1 Hz, with step length of 0.001 Hz. As C_ap(ω) and C_nu(ω) are complex, we present the phase difference and amplitude difference, respectively.

3.1   Phase deviation of C_ap(ω) from C_nu(ω) and its influence on phase velocity measurement

Figure 3 displays the phase difference δφ between C_ap(ω) and C_nu(ω), defined by

Figure  3.  Phase difference between the approximate and the exact solution of cross correlation function C(ω). The vertical axis is phase difference (in radian) over π. Curves of different color are for different receiver distance R, as shown by the legend. N is the number of wavelengths between the two stations, R/λ. The dots with different color and sizes denote the ends of the curves, respectively.

DownLoad: Full-Size Img PowerPoint

It should be noticed that the x axis of Figure 3 is not frequency or angular frequency, but the number of wavelength between the two receivers, defined by

Curves of different color are for different receiver distance R (black, blue, and red are for R = 50, 100, 200 km, respectively). The dots with different colors and sizes denote the ends of the curves, respectively. The common parts of the three curves are identical, which is expected by the scale analysis. It is very obvious from Figure 3 that there is indeed phase difference between C_ap(ω) and C_nu(ω). The deviation oscillates with smaller amplitude as N increases. When N < 1, the deviation could exceed π/10 and when N is big enough (N>30), the deviation will be below π/50. Figure 4 is the zoomed view of the N < 10 part in Figure 3. In many studies, only the station pairs whose distance is above several wavelengths are used, and then the phase difference is up to π/20 according to Figure 4.

Figure  4.  Zoomed view of the N < 10 part in Figure 3.

DownLoad: Full-Size Img PowerPoint

It is straightforward to translate the phase deviation into errors of phase velocity. From equation (11) we know that the correct formula to measure the phase velocity is

where T is the period, t is the travel time of the surface wave. However, if C_ap(ω) is taken as exact Green's function, the formula of phase velocity measurement will be written as

which is without the phase deviation δφ. Then we define the relative error as

This means that the relative error introduced by the phase deviation decreases very fast as N increases. The result of δc/c′ is shown in Figure 5, also with the x axis of N < 10. When N is small (N < 2), the relative error is substantial (>3%). However, for the receiver distances used in most of seismological studies (N>5), the relative error is less than 1%.

Figure  5.  Relative error of phase velocity introduced by the approximation. N is the number of wavelength between receivers, and the vertical axis is the ratio between phase velocity error and the true phase velocity.

DownLoad: Full-Size Img PowerPoint

3.2   Amplitude deviation of C_ap(ω) from C_nu(ω)

The amplitude of Green's function (equation (1)) contains the information of geometrical spreading and attenuation of the medium. The amplitudes of C_ap(ω) and C_nu(ω) (denoted as A_ap and A_nu respectively in the following) for different values of R are shown in Figure 6. The red curves stand for approximate solutions and the blue ones for the numerical solutions. Although they agree with each other fairly well, the numerical amplitude oscillates around the approximate amplitude. To quantitatively present the amplitude differences between C_ap(ω) and C_nu(ω), we define the relative error of amplitude, just as what we do with the phase velocity.

Figure  6.  Amplitude of approximate and exact solution of cross correlation function C(ω).

DownLoad: Full-Size Img PowerPoint

The result is shown in Figure 7. The x axis is the number of wavelength between the two receivers, A and B. The dots with different colors and sizes denote the ends of the curves, respectively. The relative error of amplitude has similar pattern with the phase difference shown in Figure 3, which is oscillating around zero and with decreasing magnitude. When N < 3, the mismatch is more than 10%. But for R used in most studies, the amplitude difference is less than 10%.

Figure  7.  Relative difference of amplitude between the approximate and the exact solution of cross correlation function C(ω). The vertical axis is the ratio between amplitude error and the true amplitude. The dots with different colors and sizes denote the ends of the curves, respectively.

DownLoad: Full-Size Img PowerPoint

4.   Discussion and conclusions

The theoretical framework (Snieder, 2004, 2007) is very helpful for understanding emergence of Green's function from cross correlation of ambient noise. His approach is more appropriate for application of ambient noise to seismology, as it neither require equipartition of energy, nor time-reversal invariance of wave equation in elastic medium, which seems to be absent for our Earth. However, a stationary phase approximation was introduced to deal with the integral equation (5) over the whole plane in the derivation, which seems to work well in high frequency band. To study the error resulted from the approximation, we numerically evaluated the integral equation (5) to get a more accurate result and compare it with the real Green's function. We find that the approximation is effective as the exact and approximate solution agree with each other generally. However, it indeed introduces some errors to equation (1). There is an oscillating deviation of both amplitude and phase in the approximate solution from the exact solution, which is around zero and attenuating when the frequency increases. This means that the stationary phase approximation is not appropriate when the frequency is relatively low (that is, the number of wavelength between receivers are very small), but suitable for the high frequency condition.

Although theoretically the relationship shown in equation (1) is not accurate, it does not affect most applications of ambient seismic noise method, as in almost all the researches only the station pairs whose distances are above several wavelengths are used to compute the cross correlation. Our computation implies that if the receiver distances are above several wavelengths, the error introduced by the approximation would be negligible. We recommend that only the station pairs whose distances are above 5 wavelengths can be used.

In this research, it is still assumed that the scattered waves propagate isotropically near the receivers. As we treat the scatterers as isotropic, mathematically this assumption is expressed by the condition that the scatterer density n is constant in space. However, many practical researches found that this might not be the actual case. Sabra et al. (2005) found that the SNR of cross correlation function is higher for station-pairs oriented perpendicular to the coast line. Pederson and Kruger (2007) found that in Baltic shield the ambient seismic noise is very close to a plane wave in some frequency band. Some researchers even use this property to study the characteristic of ambient seismic noise (Stehly et al., 2006). All of these imply that it still need more theoretical and practical researches on the relationship between cross correlation function and the Green's function.