1 Introduction
The Green's function of the wave equation in an inhomogeneous medium can be estimated by cross correlating signals emitted by ambient noise sources and recorded by a passive sensor array (Bardos et al., 2008; Colin de Verdière, 2009; Shapiro et al., 2005; Wapenaar and Fokkema, 2006; Weaver and Lobkis, 2001). The main result is that the cross correlation of the signals u (t,x1) and u (t,x2), recorded at two sensors x1 and x2 defined by:
(1) |
(2) |
Cross correlations of signals emitted by ambient noise sources and recorded by a passive sensor array can also be used for imaging of reflectors imbedded in the medium (Garnier and Papanicolaou, 2009, 2010; Gouédard et al., 2008). The data that is used for imaging is the cross correlation matrix between pairs of sensors of the passive array (xj)j = 1,...,N. In Section 4 we consider imaging functionals that use the cross correlation matrix to image the reflectors. The forms of these functionals are motivated by the analysis of the cross correlation matrix in the asymptotic regime where the coherence time of the sources is smaller than the typical travel times to be estimated and used. The imaging functionals can then be applied in general, in a smooth or randomly scattering medium.
For travel time estimation and reflector imaging we pay particular attention to the statistical stability of the estimation. There are two types of statistical stability in these problems:
- • First there is the issue of statistical stability with respect to the distribution of the noise sources. This means that the empirical cross correlation CT depends on the realizations of the signals generated by the ambient noise sources. However CT is self-averaging (or statistically stable) in the sense that the time average in (1) over an interval of duration T tends to a statistical average when T is large enough, as shown in Section 2. Therefore, statistical stability relative to the distribution of the noise sources can be controlled through the choice of a long enough recording time-window or by stacking techniques;
- • Second there is the issue of statistical stability with respect to the distribution of the scattering medium. This means that the cross correlation depends on the realization of the randomly scattering medium. It is not in general a self-averaging quantity. Indeed we will see in this note that fluctuations of the cross correlation due to the randomly scattering medium can have a large standard deviation that depends on the spectrum of the noise sources and on statistical properties of the scattering medium. Furthermore, statistical stability with respect to the distribution of the random medium cannot be controlled in general. In Sections 3.2 and 4.3 we analyze the trade-off between resolution enhancement and signal-to-noise ratio (SNR) reduction due to scattering by a random medium. We show in Section 3.3 and 4.4 that the use of iterated cross correlations can improve the SNR. We also carry out some simple numerical simulations to illustrate the results. They confirm that the theoretical predictions obtained by asymptotic analysis can be seen in realistic configurations.
2 Stability with respect to the noise source distribution
We consider the solution u of the scalar wave equation in a three-dimensional smooth medium with speed of propagation c(x):
(3) |
(4) |
The time distribution of the noise sources is characterized by the correlation function F(t2–t1), which is a function of t2–t1 only by time stationarity. The Fourier transform of the time correlation function F(t) is a nonnegative, even, real-valued function proportional to the power spectral density of the sources (Bochner's theorem):
(5) |
The spatial distribution of the noise sources is characterized by the autocovariance function δ(y1 − y2)K(y1). The process n is delta-correlated in space and K characterizes the spatial support of the sources. It is possible to consider a more general form for the spatial autocovariance function (Bardos et al., 2008). The analysis then requires the use of semi-classical methods but the results do not change qualitatively.
We denote by G(t,x,y) the time-dependent outgoing Green's function. It is the fundamental solution of the wave equation
(6) |
The empirical cross correlation of the signals recorded at x1 and x2 for an integration time T is defined by (1). It is a statistically stable quantity, in the sense that for a large integration time T, the empirical cross correlation CT is independent of the realization of the noise sources. We have the following results (the first two items are proved in Garnier and Papanicolaou, 2009):
- • The expectation of the empirical cross correlation CT (with respect to the distribution of the sources) is independent of T:
(7) |
(8) |
- • The empirical cross correlation CT is a self-averaging quantity:
(9) |
- • The covariance of the empirical cross correlation CT is:
(10) |
(11) |
This shows that:
- • all noise sources participate in the fluctuations of the empirical cross correlation (since the volume integral of the source function K appears in (11));
- • the standard deviation of the fluctuations decays as (BT)–1/2 (the square root of the integration time T times the noise bandwidth B) and as L–2 (the square distance from the sources to the sensors).
Therefore errors may occur when the averaging time is not sufficient to ensure that time averages approximate statistical mean values (Gouédard et al., 2008). However, the integration time is usually not a limiting factor, therefore this error can be brought to an arbitariliy small value and we will neglect it in the following.
3 Travel time estimation
If we consider the expression (8) of the cross correlation when the noise sources surround the region of interest that contains the sensors, then the term in the square brackets can be computed by the Helmholtz-Kirchhoff identity: this identity yields that this term is proportional to the imaginary part of the Green's function and we then find that the derivative of the cross correlation is proportional to the symmetrized Green's function convoluted with the time correlation function of the noise sources (Garnier and Papanicolaou, 2009, Section 4.4):
(12) |
Travel time estimation is then straightforward as the singular component of the Green's function G(τ, x1, x2) is concentrated at time lag τ equal to the travel time . When the noise sources are spatially localized the Helmholtz-Kirchhoff identity cannot be used and we address in this section travel time estimation in this situation.
3.1 Analysis in the high-frequency regime
We assume from now on that the decoherence time of the noise sources is much smaller than the typical travel time that we want to estimate. If we denote by ɛ the (small) ratio of these two time scales, then we can write the time correlation function Fɛ of the noise sources in the form
(13) |
- • in applications, such as surface wave tomography, noise records are first bandpass-filtered and then cross correlated (Shapiro et al., 2005). If the central frequency ω0 of the filter is high enough so that the corresponding wavelength λ0 is much smaller than the typical travel distance d, then we have ɛ = λ0/d ≪ 1. As we will see below, the resolution by cross correlation is inversely proportional to the bandwidth, so the hypothesis ɛ ≪ 1 turns out to be natural in order to get some resolution (i.e. an estimate of the travel time with an accuracy smaller than the total travel time through the region);
- •
the Fourier transform of the time correlation function of the sources has the form , so that the cross correlation (8) involves a product of Green's functions evaluated at high frequencies:
(14)
When ɛ ≪ 1 we can use a geometric asymptotics approximation for the Green's function
(15) |
(16) |
(17) |
Note that is not zero only if the ray starting from x2 and going through x1 extends into the source region. In other words, sources located along the ray starting from x1 with the direction of x1–x2 (resp. x2–x1) give rise to a singular component at (resp. ). In Fig. 1a–b one can see a peak at because the ray going through x1 and xj intersects the source region. In Fig. 1c–d one cannot see any peak at because the ray going through x1 and xj does not intersect the source region.
To summarize, when ɛ is small, the cross correlation C(1)(τ, x1, x2) has singular components (at ) if and only if the ray going through x1 and x2 reaches into the source region, that is, into the support of the function K. The results (16), (17) also show that:
- • only the noise sources located in a small tube around the ray joining x1 and x2 contribute to the singular components of C(1)(τ, x1, x2) (this can be seen from the line integral (17));
- • the widths of the peaks are determined by the bandwidth of the noise sources;
- • the heights of the peaks do not depend on the distance from the sources to the sensor array.
This last property follows from the stationary phase analysis and is a consequence of two competing phenomena that cancel each other: on the one hand the geometric decay of the amplitude of the Green's function as a function of the distance from the sources to the sensors, and on the other hand the increase of the diameter of the small tube around the ray that contributes to the singular components.
3.2 Signal-to-noise ratio reduction and enhanced resolution due to scattering
In order to analyze the cross correlation technique in a scattering medium, we first introduce a model for the inhomogeneous medium. We assume that the medium has a homogeneous background speed c0 and small and weak fluctuations responsible for scattering:
Here E stands for the expectation with respect to the distribution of the randomly scattering medium, σs is the standard deviation of the fluctuations of the speed of propagation, ls is the correlation length of the fluctuations, and the function ρ(x) characterizes the spatial support of the scattering region. Note that we have assumed that the correlation length is small enough so that we can consider that the process V is delta-correlated.
To leading order in the scattering strength the average cross correlation is the sum of the unperturbed cross correlation (i.e. the cross correlation (8) in the absence of scattering) and of an additional term of the form:
(18) |
(19) |
Note that the support of the function Kρ is the support of ρ, i.e. the scattering region. Kρ(z) is proportional to the total power reemitted from z by scattering. Eq. (18), which has the same form as (8) but with Kρ instead of K, shows that the random scatterers play the role of secondary sources. The average cross correlation (18) possesses additional peaks at (resp. ) due to the scattered waves provided that there are rays issued from the scattering region that goes through x1 and then x2 (resp. through x2 and then x1). Indeed in the high-frequency regime we find
(20) |
- • all noise sources but only the scatterers along the ray joining x1 and x2 participate in the additional peaks;
- • the heights of the additional peaks decay with the square distance from the sources to the scattering region.
To leading order in the scattering strength and in the high-frequency regime, we find using again stationary phase arguments that the variance of the fluctuations of the cross correlation is:
(21) |
This shows that:
- • all noise sources and all scatterers participate in the fluctuations of the cross correlation due to the randomly scattering medium (and not only the ones along a particular ray);
- • the standard deviation of the fluctuations decay with the square distance from the sources to the scattering region and with the square distance from the scattering region to the sensors;
- • in terms of the noise bandwidth B, the standard deviation scales as B1/2 while the amplitudes of the peaks scale as B, which shows that the relative fluctuations decay with the bandwidth.
The analysis of the mean and variance of C(1) therefore shows that scattering can enhance the directional diversity of the wave fields recorded by the sensors, which can help in travel time estimation, but it also increases the fluctuations of the cross correlation, which may make the additional peaks difficult to detect. The use of fourth-order cross correlations discussed in the next subsection is a way to increase the SNR of travel time estimation in randomly scattering media.
3.3 Use of fourth-order cross correlations
It is possible to estimate the travel time between two sensors x1 and x2 in a scattering medium by looking at the main peaks of a special fourth-order cross correlation (C(3) stands for Correlation of Coda of Correlation) (Stehly et al., 2008). This fourth-order cross correlation uses the data recorded by an array of auxiliary sensors xa,k, k = 1,...,Na, and it is evaluated as follows:
- • calculate the cross correlations between x1 and xa,k and between x2 and xa,k for each auxiliary sensor xa,k as in (1);
- • calculate the coda (i.e. the tails) of these cross correlations:
- • cross correlate the tails of the cross correlations and sum them over all auxiliary sensors to form the coda cross correlation between x1 and x2:
This algorithm is parameterized by three important times:
- • the time T is the integration time and it should be large so as to ensure statistical stability with respect to the distribution of the noise sources;
- • the time Tc1 should be large enough so that the Green's functions and limited to [Tc1,Tc2] do not contain the contributions of the direct waves. This means that Tc1 depends on the index of the auxiliary sensor k and should be a little bit larger than ;
- • the time Tc2 should be large enough so that the Green's functions and limited to [Tc1,Tc2] contain the contributions of the incoherent scattered waves. This means that Tc2 should be of the order of the power delay spread.
The coda cross correlation is a self-averaging quantity and it is equal to the statistical coda cross correlation C(3) as T → ∞:
The coda cross correlation C(3) differs from the cross correlation C(1) in that the contributions of the direct waves are eliminated and only the contributions of the scattered waves are taken into account (note that some of the contributions of scattered waves are also eliminated, but only those which correspond to small additional travel times, which are also those which induce small directional diversity). Since scattered waves have more directional diversity than the direct waves when the noise sources are spatially localized, the coda cross correlation C(3)(τ, x1, x2) usually exhibits a stronger peak at lag time equal to the inter-sensor travel time (Garnier and Papanicolaou, 2009). In particular, in contrast with the cross correlation C(1), the existence of a singular component at lag time equal to the travel time does not require that the ray joining x1 and x2 reaches into the source region, but only into the scattering region.
We illustrate these results in Fig. 2 in which the five sensors are aligned perpendicularly to the energy flux coming from the noise sources and the cross correlation C(1)(τ, x1, xj) does not have a peak at lag time equal to the travel time between the sensors x1 and xj, j = 1,...,5. However the ray going through the sensors x1 and xj intersects the scattering region and the coda cross correlation C(3)(τ, x1, xj) has a peak at lag time equal to the travel time between the sensors. The comparison between Fig. 2c, d also shows that the auxiliary sensors are necessary to get the peak of the coda cross correlation.
4 Reflector imaging
In this section we show that cross correlations of signals emitted by ambient noise sources and recorded by a sensor array can also be used for imaging of reflectors. The data to be used for imaging is the matrix of cross correlations between pairs of sensors of an array . The objective is to image the reflectors buried in the medium. Here we consider the case in which there is only one reflector located at zr. Note that it is often necessary to have available data sets and , with and without the reflector, respectively, so that we can compute the differential cross correlations ΔC = C − C0 and migrate them. In general, but not always, the primary data set C cannot be used directly for imaging because peaks in the cross correlations due to the reflector may be very weak compared both to the peaks of the direct waves, at lag times equal to the inter-sensor travel times, as well as to the non-singular components due to the directionality of the energy flux (Garnier and Papanicolaou, 2009). In many applications where we want to image localized changes in the environment, such as in reservoir or volcano monitoring, both data sets are usually available.
4.1 Analysis in the high-frequency regime
By a stationary phase analysis and a Born approximation for the reflector the singular components of the differential cross correlations C(1)(τ, x1, x2) can be identified (Garnier and Papanicolaou, 2010). This is important because this information will in turn allow us to determine the appropriate imaging functional that should be used to migrate the cross correlations. There are two main types of configurations of sources, sensors, and reflectors:
- • the noise sources are spatially localized and the sensors are between the sources and the reflectors (Fig. 3a). We call this the daylight illumination configuration. In this configuration the singular components of the differential cross correlation are concentrated at lag times equal to (plus or minus) the sum of travel times ;
- • the noise sources are spatially localized and the reflectors are between the sources and the sensors (Fig. 4a). We call this the backlight illumination configuration. In this configuration the singular components of the differential cross correlation are concentrated at lag times equal to the difference of travel times .
In the backlight imaging configuration, when ɛ is small, the differential cross correlation ΔC(1) has a unique singular contribution at lag time equal to the difference of travel times and it has the form:
(22) |
(23) |
Note that is not zero only if the ray going from xj to zr extends into the source region, which is the backlight illumination configuration, while is not zero only if the ray going from zr to xj extends into the source region, which is the daylight illumination configuration. Eq. (23) shows that:
- • only the sources located along the rays joining xj and zr contribute to the singular components;
- • the widths of the peaks are determined by the bandwidth of the noise sources;
- • the heights of the peaks are inversely proportional to the square of the distance from the reflector to the sensor array, and they do not depend on the distance from the sources to the reflector or the sensor array.
4.2 Migration of cross correlations
As shown in the previous subsection, the cross correlation may have additional peaks at lag times that depend on the reflector location. This suggests that the reflector can be imaged by migrating the cross correlation matrix. This form of passive sensor array imaging depends in an essential way on the illumination configuration, that is, the relative positions of the sensors, the noise sources, and the reflector.
To form an image in a daylight illumination configuration (Fig. 3a), we propose to use the daylight imaging functional defined for a search point zS in the image domain by:
(24) |
This functional is built by evaluating each element of the cross correlation matrix at lag time equal to the sum of the travel times and by summing the migrated matrix elements over all pairs of sensors (it is therefore a Kirchhoff-type migration functional for the cross correlation matrix). The resolution analysis of the daylight imaging functional is carried out in (Garnier and Papanicolaou, 2010). The cross range resolution for a linear sensor array with aperture a is given by λ0a/Lr. Here Lr is the distance between the sensor array and the reflector and λ0 is the central wavelength. The range resolution for broadband noise sources is equal to c0/B where B is the bandwidth of the noise sources. The range resolution for narrowband noise sources is .
To form an image in a backlight illumination configuration (Fig. 4a), we propose to use the backlight imaging functional defined for a search point zS in the image domain by:
(25) |
This functional is built by evaluating each element of the cross correlation matrix at lag time equal to the difference of the travel times and by summing the migrated matrix elements over all pairs of sensors. The cross range resolution is λ0a/Lr while the range resolution is whatever the bandwidth, which means that the range resolution is very poor compared to the daylight imaging functional, because it exploits only differences of travel times, which are much less sensitive to the range than the sums of travel times (Garnier and Papanicolaou, 2010).
4.3 Signal-to-noise ratio reduction and enhanced resolution due to scattering
We here discuss the trade-off between resolution enhancement and SNR reduction due to scattering.
When the primary energy flux only gives a backlight illumination of the reflector and when there is no scattering, only the backlight imaging functional (25) can be used, which produces an elongated peak at the reflector location with poor range resolution. When there is scattering, the scatterers can play the role of secondary sources and they can provide a daylight illumination for the reflector. This happens provided there are rays issuing from the scattering region and going to the sensors and then to the reflector. Scattering can therefore lead to the appearance of a singular contribution in the cross correlation at lag time equal to (plus or minus) the sum of travel times and then the daylight imaging functional (24) can have a peak at the reflector location with good range resolution. However, the scattered waves also involve fluctuations in the cross correlations that can be larger than the additional peak exhibited here above. As a consequence, the peak in the daylight imaging functional (24) at the reflector location that we have just mentioned is usually buried in the contributions of the non-singular random components. This can be studied following the lines of Subsection 3.2 and this happens in the setup of Fig. 5, which we discuss in the next subsection.
4.4 Use of fourth-order cross correlations
We have noticed that the peaks produced by the scattered waves in the differential cross correlation ΔC(1)(τ, xj, xl) and that are relevant to imaging of reflectors can be buried in fluctuations. This happens when the SNR of the peaks at lag time equal to (plus or minus) the sum of travel times is low. In this section we propose to use an iterated cross correlation technique that masks the contributions of the direct waves and increases the effective SNR of the peaks produced by the scattered waves. This technique was shown to be efficient for inter-sensor travel time estimation in Section 3.3. In the following we describe this technique for reflector imaging with secondary daylight illumination from scattering.
It is possible to form a special fourth-order cross correlation matrix between sensors from the differential cross correlations ΔCT(τ, xj, xl) obtained from the recorded data. This is done as follows:
- • calculate the coda (i.e. the tails) by truncation of the differential cross correlations:
- • cross correlate the tails of the differential cross correlations and sum them over all complementary sensors in the array to form the coda cross correlation between xj and xl:
(26) |
The roles of the three parameters T, Tc1, and Tc2 are described in Subsection 3.3. The differential coda cross correlation is a self-averaging quantity and it is equal to ΔC(3) as T → ∞:
The time windowing is very important because it selects the contributions that we want to use for imaging the reflector at zr. The asymptotic analysis of the functional ΔC(3) can be carried out in the high-frequency regime using stationary phase arguments. It shows that the differential coda cross correlation ΔC(3) has singular components at lag time equal to (plus or minus) the sum of travel times , and has less additional terms than the usual differential cross correlation studied in Subsection 4.3. As a consequence migration of the differential coda cross correlation using the daylight migration functional can produce an image of the reflector with a higher SNR.
In Fig. 5 we consider a situation in which the primary energy flux is backlight and scattering generates a secondary daylight illumination. The daylight imaging functional applied with ΔC(1) does not exhibit a clear peak at the reflector location (Fig. 5c) because of the strong fluctuations of the cross correlation at lag time equal to the sum of travel times. However, the daylight imaging functional applied with ΔC(3) (Fig. 5d) gives a much better image. The overall result is that the backlight imaging functional IB used with ΔC(1) has good cross-range resolution and SNR while the daylight imaging functional ID used with ΔC(3) has good range resolution and SNR.
5 Further results on cross correlation in scattering media
In this note we have shown that imaging with ambient noise sources in a random medium is statistically stable with respect to the noise source distribution, but it is usually not statistically stable with respect to the distribution of the random medium. It is, however, possible in some situations to exploit the enhanced diversity of the scattered waves in order to improve the resolution in travel time tomography or in reflector imaging without reducing too much the signal-to-noise ratio. This is especially true when using fourth-order cross correlations which tend to increase the signal-to-noise ratio.
These results were obtained here in a weakly scattering medium, in a high-frequency regime, but other propagation regimes can be analyzed. In Garnier and Sølna (2010b) the role of random scattering in the Green's function estimation is studied in the radiative transport regime. This work provides a bridge between the situation with a directional wave field and a fully equipartitioned field. In the radiative transport regime, scattering improves Green's function estimation by enhancing the directional diversity of the waves recorded by the sensors under the following conditions: first the sensors need to be close to each other (that is, closer than the mean free path), while the observation region can be far from the source region. Resolution is enhanced when the propagation distance from the sources to the observation region is larger than the mean free path. In Garnier and Sølna (2010a) and de Hoop and Sølna (2009) other regimes of propagation are considered (randomly layered media and random paraxial regime) which lead to similar conclusions. In Garnier and Sølna (2010a,b), de Hoop and Sølna (2009) as in the weakly scattering high-frequency regime discussed in this note, the important practical issue is the lack of statistical stability of the cross correlations with respect to the distribution of the random medium. In general, averaging over the medium is not possible but it may be possible to smooth or average over mid-points and offsets and to stabilize the data in order to invert for medium parameters in the case of simple, low-dimensional medium models. For instance a robust algorithm to detect an interface in the medium is proposed in Garnier and Sølna (2010a,b), de Hoop and Sølna (2009). The optimal way to combine data in space and frequency to obtain stable and high-resolution estimates remains an open problem.