Analyzing the coda from correlating scattered surface waves Open Access

Numerical and laboratory experiments have shown how cross correlations of field fluctuations recorded at two points in open and closed systems provide an estimate of the anti-causal and causal Green’s functions between these points.^1–5 Wapenaar and Fokkema⁶ showed how a reciprocity theorem relaxes the necessity for sources of the field fluctuations to be present throughout the volume to sources on a surface surrounding these points. These insights have made a large impact on numerous ocean acoustic ambient noise studies (e.g., Refs. 7–9), but actually date back to Aki.¹⁰ Dispersive ballistic surface waves extracted by cross correlation of ambient noise measurements are now a major tool for imaging the Earth’s crust.^11–13 Moreover, coda waves derived from these correlations have been used for subsurface monitoring.^14,15 We describe here how the coda can be exploited in this correlation technique to locate individual scatterers, and we investigate the accuracy of the recovered coda using an approximate acoustic method rather than the complete elastodynamic Green’s function retrieval method.¹⁶ 

We estimate the ultrasonic Rayleigh-wave Green’s function for surface waves propagating in an aluminum block with 15 scatterers at the surface. In the far field, the acoustic Green’s function G for a medium with constant velocity c and density ρ can be approximated by summing cross correlations of the measured pressure field u [Eq. (32) of Wapenaar and Fokkema⁶]

where G(x_A, x_B, ω) is the vertical component of the Rayleigh-wave Green’s function at x_A from a source at x_B and ℐ denotes the imaginary part. S(ω) is the power spectrum of the sources on the closed contour ∂D and u* denotes the complex conjugate. Strictly speaking, expression (1) is only valid for acoustic media; however, Halliday and Curtis¹⁷ show that G is equivalent to the G_zz component of the elastic surface wave Green’s tensor when c is the Rayleigh-wave speed and u is the vertical component of the displacement. We refer to expression (1) as the acoustic approximation for Green’s function retrieval. Transforming expression (1) to the time domain states that the sum of the cross correlations of the wave fields from sources on the boundary ∂D is proportional to the causal and anti-causal Green’s function that describes wave propagation between x_A and x_B.

One side of an aluminum block (280 mm × 230 mm × 215 mm) contains cylindrical holes 1 mm in diameter and 3 cm deep [Fig. 1(a)]. A high-powered pulsed Nd:YAG (neodymium-doped yttrium aluminum garnet) laser generates ultrasonic waves by briefly (15 ns) heating a 1 mm wide circular spot on the surface. The heating causes thermoelastic expansion and generates broadband ultrasonic waves having a central frequency of 600 kHz. The Rayleigh-wave velocity in the aluminum is c = 2.9 mm/μs, corresponding to a dominant wavelength of ∼5 mm. We measure the vertical component of the displacement as a function of time at a point on the surface of the cube with a laser interferometer, based on a constant-wave 250 mW Nd:YAG laser at 532 nm. The receiver uses two-wave mixing in a photorefractive crystal and is calibrated to measure the absolute out-of-plane displacement field. The measured response is flat between 20 kHz and 20 MHz, with sensitivity on the order of parts of angstroms (see Blum et al.¹⁸ for a complete description).

With the source laser at x_A [star, Fig. 1(a)], we record the wave field at each receiver position (triangles) for 60 μs at a sample rate of 20 MS/s. There are 412 receivers in total, spaced ∼1 mm apart surrounding the source and scatterers. At each receiver position, we repeat the source excitation 128 times and average the recordings to increase the signal-to-noise ratio. The displacement wave fields are presented in Fig. 1(b), after we apply a cosine taper to the first 5 μs to suppress electronic noise from the laser recording system and a band-pass filter between 100 and 2000 kHz. We mute reflections from the edge of the block, when these arrive before t = 60 μs. The wave fields are displayed as a function of receiver azimuth (0 ≤ θ < 360°), defined with respect to the center of the receiver array, with positive θ in the clockwise direction and θ = 0 in the upper left corner of the square. The arrival with the largest amplitude is the direct Rayleigh wave, followed by scattered Rayleigh waves with the opposite polarity caused by scattering at the air-filled holes. The four sharp variations in the coherent arrivals are due to the rectangular nature of the underlying acquisition geometry.

According to expression (1), the Green’s function between two sources can be estimated from receivers on ∂D based on reciprocity of the wave equation.^19,20 Here we estimate the vertical component of the full displacement wave field—both direct and scattered waves—between sources rather than receivers, because with our acquisition setup it is easier to occupy many receiver positions than source positions. We begin by estimating the displacement wave field at a single point using the auto correlation. The auto correlation of the wave field for each receiver as excited by the source at x_A is shown in Fig. 2(a), which henceforth we call a correlation gather. The sum of these auto correlations is shown in Fig. 2(b). As we have not taken the time derivative as defined in expression (1), it represents a phase-shifted causal and anti-causal response for a coinciding source and receiver at x_A. In the auto correlation gather, the “direct” wave is instantaneous (arriving at t = 0), while scattered waves vary smoothly in time as a function of receiver position, θ. The extrema of these curves define the stationary-phase points, where scatterer, source, and receiver are in-line.

The time of a given stationary-phase point in the auto correlation gather (and knowledge of the velocity c) defines the distance r between the source and scatterer. Figure 2(c) illustrates that the stationary receiver at θ = 34° and the arrival time t ≈ −26 μs correspond to a scatterer at a distance r = 37.5 mm from x_A, which agrees with scatterer 1. Because the auto correlation is symmetric in time, there is a similar stationary-phase point at t ≈ +26 μs. Another stationary-phase point exists at θ ≈ 34° with an earlier arrival time of t ≈ ±10 μs, corresponding to scatterer 4. This means the scatterer is in the same direction as scatterer 1, but closer to x_A.

To determine the accuracy of the coda retrieval, and because it is not trivial to measure the impulse response directly for a coinciding source and receiver position, we repeat the same acquisition described above for a source at x_B [Fig. 1(a)]. In this case, between θ = 60° and 180° the side reflection from the block edge interferes with waves scattered from some of the cylindrical holes. Muting the side reflection also removes some of the scattering coda. Next, we cross correlate the processed wave fields recorded at each receiver from the two sources [Fig. 3(a)] and sum [Fig. 3(b)]. Following expression (1), each trace in the sum is weighted by the difference in angle between the receiver and proceeding receiver.

Toward estimating the wave field that propagates between x_A and x_B, we time-reverse the anti-causal weighted sum and stack with the causal weighted sum to improve the signal-to-noise ratio. As prescribed by expression (1), we then take a time derivative. Because the correlation contains a contribution S from the source wavelet, we deconvolve a tapered version of the direct Rayleigh-wave arrival from our result using a water-level deconvolution.²¹ The tapered wavelet used in the deconvolution is shown in the inset of Fig. 3(c).

We compare our estimate of the wave field propagating between a source x_A and a virtual receiver x_B to a direct measurement for a source at x_A and a real receiver at x_B [Fig. 3(c)]. The amplitudes of the two wave fields are normalized by the maximum of the direct Rayleigh wave. The direct surface wave is extracted with great accuracy. The amplitude of the coda—relative to the direct wave—is accurately estimated, and some arrivals in the coda are well reproduced, while others are not. Differences between the direct estimate of the impulse response and the result based on expression (1) are caused by violations of the assumptions that went into the approximate expression (1). An exact formulation for the extraction of the full elastodynamic Green’s function has been developed.¹⁶ This formulation requires, however, multicomponent point forces and double couples that excite the field fluctuations and is difficult to implement in laboratory and field experiments. Specifically, in this experiment we (1) suppress some scattered waves when muting the side reflection, (2) violate the far-field assumption, (3) apply an approximate deconvolution, and (4) measure the wave field at a finite number of locations on ∂D rather than continuously in space. Halliday and Curtis¹⁷ justified the use of this acoustic formulation based on its accuracy for extracting the direct surface wave, not necessarily on scattered arrivals. One aspect of recovery of the coda in this technique is that we are able to identify which scatterer is responsible for a certain time-section of the coda.

In theory—and previously only shown numerically—the accurate extraction of scattered waves from the correlation of field fluctuations relies on the cancellation of correlated energy with a constant arrival time. These waves come from the correlation of waves traveling from a source to a particular scatterer and onto both receivers; the arrival time of this phase is independent of source location.²² Its cancellation is due to the generalized optical theorem that constrains the scattering coefficient. Figure 3(a) does not provide any evidence of correlations with a constant arrival time because of the dominating amplitude of the direct wave. In Fig. 4(a), the direct wave at each receiver is muted before cross correlation. In other words, it is the correlation gather for scattered waves only. Several events with a constant arrival time are now visible, some of these are highlighted with boxes around t ∼−9 μs and ∼ +5 μs. The scattering coefficient of surface waves depends, in general, strongly on the scattering angle,²³ and the events with a constant arrival time in Fig. 4(a) are strongest when the source, receiver, and scatterer are aligned. Because there are 15 scatterers in the model, there should be 15 waves with a constant arrival time. It is difficult to verify this because the correlation gather is complicated by cross terms of waves radiated by different scatterers (whose arrival time depends on the receiver angle θ).

The wave field from the summed correlation gather of only the scattered waves is shown in Fig. 4(b). Note that in contrast from the wave field extracted from the cross correlation of the full wave field [Fig. 3(b)], Fig. 4(b) does not show the direct wave at all. In fact, the summed wave field is dominated by arrivals t ∼ −9 μs and t ∼ +5 μs that correspond to the sum of the waves in the boxes in the correlation gather with a constant arrival time. These spurious arrivals should not contribute to the estimated impulse response, and theory and numerical examples show that these waves cancel by stationary phase contributions from a cross term of direct and scattered waves.²² Since we mute the direct wave, such cancellation does not take place, with the result that these spurious arrivals dominate in Fig. 4(b). It has been shown theoretically that in order to retrieve the scattered waves, one needs the cross correlation of the direct wave with scattered waves.²⁴ Since the direct waves are muted in the correlation gather of Fig. 4(a), one cannot expect that the waves in Fig. 4(b) accurately represent the scattered waves.

We present laboratory measurements of scattered ultrasonic surface waves that propagate in an aluminum block with cylindrical holes that act as scatterers. An estimate of the wave field between two source locations can be obtained by summing the cross correlated wave fields measured at receiver locations enclosing the sources and scatterers. The correlation gather computed from the recorded wave fields shows distinct arrivals related to individual scatterers, and we show that the stationary-phase points of these arrivals can be used to constrain the location of scatterers.

The accuracy of our Green’s function estimate can be compared by directly measuring the response at the location of one of the two sources. The direct wave is reconstructed quite accurately, but the reconstruction of the coda is less accurate. The acoustic formulation for Green’s function retrieval is justified on arguments that apply to the direct wave only, and we attribute the reduced ability to extract the coda to the inaccuracy of using the acoustic formulation for scattered elastic surface waves.

Finally, theory predicts that for accurate recovery of the wave field, the correlation gather should contain energy with a constant arrival time from the correlation of scattered arrivals from the same scatterer at each receiver. We show the presence of such arrivals in laboratory data, after muting the direct wave before cross correlation of the wave fields. Summing this cross correlation gather does not lead to a reconstruction of either the direct or scattered waves. This agrees with theory, which predicts that the cross correlation of direct and scattered waves are essential in Green’s function extraction by cross correlation.