Global acoustic daylight imaging in a stratified Earth-like model

TL;DR

This paper analyzes acoustic daylight imaging in a stratified Earth-like model, showing how to recover eigenfrequencies and the scattering operator from surface measurements using autocorrelation functions.

Contribution

It establishes a relation between autocorrelation functions of surface measurements and the scattered wave field, enabling eigenfrequency and scattering operator recovery.

Findings

Eigenfrequencies can be recovered from autocorrelation.

Complete scattering operator can be extracted under specific geometric conditions.

Provides a theoretical foundation for acoustic daylight imaging in stratified media.

Abstract

We present an analysis of acoustic daylight imaging in an Earth-like model assuming a random distribution of noise sources spatially supported in an annulus located away from the surface. We assume a situation with scalar wave propagation and that the measurements are of the wave field at the surface. Then, we obtain a relation between the autocorrelation function of the measurements and the trace of the scattered field generated by an impulsive source localized just below the surface. From this relation it is, for example, clear that the eigenfrequencies can be recovered from the autocorrelation. Moreover, the complete scattering operator can be extracted under the additional assumption that the annulus is close to the surface and has a thickness smaller than the typical wavelength.

Equations

$$\frac{\partial }{\partial r}{r}^2\frac{\partial }{\partial r}\hat{p}+\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\sin \theta \frac{\partial }{\partial \theta }\hat{p}+\frac{1}{{\sin }^2\theta }\frac{{\partial }^2}{\partial {\phi }^2}\hat{p}+\frac{{\omega }^2{r}^2}{c^2(r)}\hat{p}=\hat{f}(\omega ,r,\theta ,\phi ),$$

$$\hat{p}(\omega ,r,\theta ,\phi )=$$

$${\hat{p}}_{l,m}(\omega ,r)$$

$${\hat{\mathcal{S}}}_l(\omega ,R_o)={\left|\frac{1+{\hat{\mathcal{R}}}_l(\omega ,R_o)}{1-{\Gamma }_{R_o}{\hat{\mathcal{R}}}_l(\omega ,R_o)}\right|}^2,$$

$${\hat{p}}_{l,m}^{\mathrm{sym}}(\omega ,R_o)={\hat{p}}_{l,m}(\omega ,R_o)-{\overline{\hat{p}}}_{l,m}(\omega ,R_o)={\hat{F}}_{l,m}^{e1}(\omega )(1-{\Gamma }_{R_o}{)}^2{\hat{\mathcal{S}}}_l(\omega ,R_o),$$

$$\displaystyle\mathbb{E}\big{[}{f}(t,r,\theta,\varphi){f}(t^{\prime},r^{\prime},\theta^{\prime},\varphi^{\prime})]=F({t-t^{\prime}})K(r)\delta(r-r^{\prime})\sin(\theta)^{-1}\delta(\theta-\theta^{\prime})\delta(\varphi-\varphi^{\prime}).$$

$$C_{l,m}^T(t,R_o)=\frac{1}{T}{\int }_0^T{p}_{l,m}(t^{\prime },R_o)p_{l,m}(t^{\prime }+t,R_o)d{t}^{\prime },$$

$${\hat{C}}_{l,m}^T(\omega ,R_o)={\hat{F}}_{l,m}^{e2}(\omega )(1+{\Gamma }_{R_o}{)}^2{\hat{\mathcal{S}}}_l(\omega ,R_o),$$

$$\hat{p}(\omega ,r,\theta ,\phi )=\int p(t,r,\theta ,\phi )e^{i\omega t}dt$$

$$f(t,r,\theta ,\phi )=f(t)g(\theta ,\phi )\delta (r-R_s),$$

$$\displaystyle\mathbb{E}\big{[}{f}(t,r,\theta,\varphi)]=0,$$

$$\displaystyle\mathbb{E}\big{[}{f}(t,r,\theta,\varphi){f}(t^{\prime},r^{\prime},\theta^{\prime},\varphi^{\prime})]=F(t-t^{\prime})K(r)\delta(r-r^{\prime})\sin(\theta)^{-1}\delta(\theta-\theta^{\prime})\delta(\varphi-\varphi^{\prime}),$$

$$\displaystyle\big{(}\partial_{r}\hat{p}-r^{-1}\hat{p}\big{)}_{r=R_{o}}=0.$$

$$Y_{l,m}(\theta ,\phi )=\sqrt{\frac{(2l+1)(l-m)!}{4\pi (l+m)!}}{P}_l^m(\cos \theta )e^{im\phi },$$

$${\int }_0^{\pi }{\int }_0^{2\pi }\overline{Y_{l,m}(\theta ,\phi )}{Y}_{l^{\prime },m^{\prime }}(\theta ,\phi )\sin \theta d\phi d\theta ={\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},$$

$$\sum _{l,m}\overline{Y_{l,m}(\theta ,\phi )}{Y}_{l,m}({\theta }^{\prime },{\phi }^{\prime })\sin {\theta }^{\prime }=\delta (\theta -{\theta }^{\prime })\delta (\phi -{\phi }^{\prime }).$$

$$\frac{\partial }{\partial r}{r}^2\frac{\partial }{\partial r}{\hat{p}}_{l,m}-l(l+1){\hat{p}}_{l,m}+\frac{{\omega }^2}{c^2(r)}{r}^2{\hat{p}}_{l,m}={\hat{f}}_{l,m}(\omega ,r),$$

$${\hat{f}}_{l,m}(\omega ,r)$$

$${\partial }_r{\hat{p}}_{l,m}(\omega ,R_o)-R_o^{-1}{\hat{p}}_{l,m}(\omega ,R_o)=0.$$

$$\frac{\partial }{\partial r}{r}^2\frac{\partial }{\partial r}{\hat{p}}_{l,m}-l(l+1){\hat{p}}_{l,m}+\frac{{\omega }^2}{c^2(0)}{r}^2{\hat{p}}_{l,m}=0,$$

$$\displaystyle j_{l}\big{(}\frac{\omega}{c(0)}R_{\delta}\big{)}\frac{\partial}{\partial r}\hat{p}_{l,m}(\omega,R_{\delta})-\frac{\omega}{c(0)}j_{l}^{\prime}\big{(}\frac{\omega}{c(0)}R_{\delta}\big{)}\hat{p}_{l,m}(\omega,R_{\delta})=0.$$

$$\frac{1}{c^2(r)}=\frac{1}{c_o^2(r)}+V(r),$$

$$\frac{\partial }{\partial r}{r}^2\frac{\partial }{\partial r}{A}_l-l(l+1)A_l+\frac{{\omega }^2}{c_o^2(r)}{r}^2{A}_l=0,$$

$$r^2\frac{\partial A_l}{\partial r}{B}_l-r^2{A}_l\frac{\partial B_l}{\partial r}=1,$$

$${\hat{\alpha }}_{l,m}(\omega ,r)$$

$${\hat{\beta }}_{l,m}(\omega ,r)$$

$${\hat{p}}_{l,m}(\omega ,r)=A_l(\omega ,r){\hat{\alpha }}_{l,m}(\omega ,r)+B_l(\omega ,r){\hat{\beta }}_{l,m}(\omega ,r).$$

$${\partial }_r{\hat{\alpha }}_{l,m}$$

$${\partial }_r{\hat{\beta }}_{l,m}$$

$$A_l(\omega ,r){\partial }_r{\hat{\alpha }}_{l,m}(\omega ,r)+B_l(\omega ,r){\partial }_r{\hat{\beta }}_{l,m}(\omega ,r)=0.$$

$$\frac{\partial {\hat{\alpha }}_{l,m}}{\partial r}(\omega ,r)$$

$$\frac{\partial {\hat{\beta }}_{l,m}}{\partial r}(\omega ,r)$$

$$\displaystyle\big{(}\partial_{r}A_{l}(\omega,R_{o})-R_{o}^{-1}A_{l}(\omega,R_{o})\big{)}\hat{\alpha}_{l,m}(\omega,R_{o})$$

Full text

††: Inverse Problems

1 Introduction

The emergence of the Green’s function from cross correlations of noisy signals, in the context of seismic exploration, was first pointed out by Claerbout [5, 6, 17]. He modeled the Earth’s crust as a half-space and considered the situation when waves are generated by unknown sources in the crust and recorded at the surface. He showed that the autocorrelation function of the such signals is the same as the signal reflected by the Earth’s crust when an impulsive source is used at the surface. The latter signal is indeed the one recorded in reflection seismology. He thereby established a formal connection between reflection seismology and wave field correlations. The process of autocorrelating signals on the surface, signals that are generated by noise sources in the interior, has come to be known as the daylight (imaging) configuration.

The physical explanation of why daylight imaging is equivalent to reflection seismology was simple and based on flux conservation. It is possible to give a mathematical proof in the case of a one-dimensional half space with radiation condition [10, 11].

Imaging of the structure of the Earth, from the core to the surface, remains a major challenge. On the one hand, the standard seismic methods using earthquake signals are limited by the uneven spatial distribution of earthquakes, along the major tectonic plate boundaries. On the other hand, seismic interferometry using ambient noise records has been successful. We here refer to seismic interferometry as using cross correlations of signals to gain useful information about the medium. Claerbout in his book [4] asks the reader in an exercise to prove that the temporal autocorrelation of a transmission seismogram for a layered model with a source underground is equivalent to a reflection seismogram, a problem he treated in [4] for a Goupillaud medium. Such a “daylight” imaging approach has been subsequently validated experimentally and mathematically for more general microstructures than layered, see for instance [11]. Seismic interferometry in general has mostly been used to probe the uppermost layers of the Earth through surface waves [2, 19]. However, recently, body-wave path responses from probing the deepest part of the Earth were obtained from noise records [3]. This motivates the question: can the daylight imaging approach as proposed by Claerbout be extended to a global image procedure for the Earth? One caveat is that in global seismology the model domain is a ball, whence the radiation condition that played a crucial role in the analysis of daylight imaging cannot be used anymore. If the noise sources were distributed uniformly, then we would have equipartition of energy amongst the normal modes, and we could invoke this argument to establish the relation between the autocorrelation function of the noise signals and the trace of the scattered field generated by localized and impulsive sources just below the surface [22]. Unfortunately, the noise source distribution is not uniform, and we cannot invoke the diffusion or ergodicity properties of the Earth to claim that the observed field is equipartitioned as in [1, 12] because the scattering is not strong enough. In fact, we will show below that in the propagation regime of interest, the relation between the autocorrelation function of the noise signals and the mentioned scattered field is not as simple as in the case of a one-dimensional half space. In the case of a spherically symmetric Earth model, we show that the point spectrum of the Earth associated with the low angular orders can be extracted from the correlation functions of the signals recorded at the surface and emitted by unknown noise sources localized away from the surface. If the noise sources satisfy some specific conditions in terms of their spatial localization then the correlations of noisy signals can be used to extract the standard scattering operator. These conditions, that the source distribution is spatially localized on a thin annulus with thickness comparable to or smaller than the typical wavelength, could be satisfied in practice. We also include weak angular variations of the wavespeed in the analysis, which confirms that global acoustic daylight imaging is possible for an essentially spherically symmetric Earth.

The paper is organized as follows. In Section 2 we give a summary of the main results of the paper. In Section 3 we introduce the wavefield decomposition and coupling in a radial Earth and the relevant asymptotic analysis. The scaling in our analysis corresponds to modes that essentially propagate in the radial direction. In Section 4 we describe the scattering operator in a stratified spherically symmetric Earth. In Section 5, we present global acoustic daylight imaging and we establish the relationship between the scattering operator and the autocorrelation functions of the ambient noise signals. We also consider the special case where the spatial support of the noise sources is a thin annulus located below the surface with thickness smaller than the typical wavelength, providing a simplified result. In Section 6, we show that our results are robust with respect to the properties of source and receiver distributions, with respect to small angular undulations in the Earth’s parameters (rather than a purely radial Earth), and with respect to measurement noise. In particular, it follows that even with small angular variations the observed point spectrum of the Earth can be reasonably accurately explained by a radial model.

2 Summary of Modeling and Main Results

In this section we summarize the main modeling assumptions, quantities of interest, and main result. We compare the symmetrized field (5) measured at the Earth’s surface and transmitted by an impulsive point source just below the surface with the empirical autocorrelation function (7) of the field measured at the Earth’s surface and generated by a distributed noise source distribution. We establish a correspondence between them and we briefly discuss aspects of robustness of the results to the modeling assumptions.

2.1 Wave Decomposition in the Radial Earth and Main Result

We consider here scalar wave propagation in the spherical Earth so that the harmonic wave field $\hat{p}=\hat{p}(\omega,r,\theta,\varphi)$ satisfies

[TABLE]

for $(r,\theta,\varphi)\in(0,R_{o})\times(0,\pi)\times(0,2\pi)$, where $c(r)$ is the radial velocity model and $\hat{f}(\omega,r,\theta,\varphi)$ is the source term in the frequency domain. We assume here that $\hat{f}$ is supported away from the origin in $\omega$ and a boundary condition at the surface $r=R_{o}$ that takes into account boundary dissipation, see Eq. (38). We moreover assume a high-frequency situation, see Section 3.3 for a discussion of the WKB context that this entails and which is central to our analysis. The wave field is then decomposed as

[TABLE]

with the spherical harmonics $Y_{l,m}$ defined in Eq. (12) and

[TABLE]

We consider the two following source scenarios.

2.1.1 Point Source.

This corresponds to a “classic” seismology configuration with a point source just below the surface. Let $\hat{\cal R}_{l}(\omega,R_{o})$ be the fundamental scattering function associated with Earth with no surface reflection, note that this is independent of $m$ as the velocity model is radial (see Proposition 1). Then the symmetrized surface Earth response function is (see Proposition 2)

[TABLE]

and this is the function that encapsulates information about the Earth’s interior. The 1 in the numerator corresponds to a direct transmission from the source and the $\hat{\cal R}_{l}$ term is the reflection from the Earth’s interior while the denominator produces multiples from reflections at the Earth’s surface with a reflectivity of $\Gamma_{R_{o}}$.

The symmetrized field measured at the Earth’s surface is then of the form

[TABLE]

with $\hat{F}^{\rm e1}_{l,m}$ the mode-dependent effective source trace and the factor $(1-\Gamma_{R_{o}})^{2}$ is an effective transmittivity factor through the Earth’s surface (see Proposition 2).

2.1.2 Distributed Random Source Field.

This is our main configuration and the main result is that under “ideal” circumstances we can recreate the surface Earth response function in Eq. (4) via forming correlations. Specifically, we model the source field as a random field delta-correlated in space located below the surface:

[TABLE]

We then have for the empirical correlation

[TABLE]

that for $T$ large:

[TABLE]

with again $\hat{F}^{\rm e2}_{l,m}$ an effective source trace (see Proposition 3).

This generalizes the classic daylight imaging result to the context of the spherical Earth. By looking at (7) we can immediately deduce that the eigenfrequencies -for which $\hat{\cal S}_{l}(\omega,R_{o})$ essentially blows up- can be extracted from the empirical correlation $\hat{C}^{T}_{l,m}(\omega,R_{o})$. Moreover, if one wants to reconstruct the Earth scattering operator $\hat{\cal S}_{l}(\omega,R_{o})$ quantitatively, then some additional assumptions have to be made on the noise source distribution so that one can extract $\hat{\cal S}_{l}(\omega,R_{o})$ from the product $F^{\rm e2}_{l,m}(\omega)\hat{\cal S}_{l}(\omega,R_{o})$. Such an additional assumption is proposed and discussed at the end of Section 5: If the spatial support of the noise source is localized in a small annulus below the surface, then the autocorrelation is directly related to the scattering operator via a classical seismic interferometry formula and one can extract $\hat{\cal S}_{l}(\omega,R_{o})$ from $\hat{C}^{T}_{l,m}(\omega,R_{o})$ up to a multiplicative frequency-dependent function (see Proposition 4).

2.2 Robustness of Main Result with Random Source Field

The assumptions that the Earth is radial, that the noise source field is statistically homogeneous in angles (Eq. (6)), and that the wave field is recorded over the whole surface may seem quite restrictive. However, as we discuss in detail in Section 6 the analysis indicates that our results are surprisingly robust with respect to these assumptions. These observations are consistent with the success of helioseismology where the “source” waves are generated by the turbulence in the convection zone immediately beneath the Sun’s surface [13] and the eigenfrequencies can be relatively robustly observed and giving information about the radial variation of the Sun’s parameters [15, 8].

We briefly comment on the robustness and refer to Section 6 for a detailed discussion.

• •

Full surface measurement aperture: In fact, exploiting spherical symmetry our result is essentially unchanged with very limited measurements, see Section 6.1.

• •

Full surface source aperture: The assumption of noise sources whose angular distribution is uniform can be slightly relaxed, see Section 6.2.

• •

A radial Earth: We can allow for small undulations in the velocity model beyond the radial case, see Section 6.3.

• •

Noisy measurements: The fact that we average the observations over different locations makes the scheme robust with respect to additive measurement noise, see Section 6.4.

3 Wave Field Asymptotics in a Stratified Ball

3.1 Helmholtz Equation in the Ball

In this subsection we formulate the Helmholtz-type equation that the mode amplitudes satisfy with the appropriate boundary conditions. In spherical coordinates $(r,\theta,\varphi)$, the time-harmonic wave field

[TABLE]

satisfies Eq. (1) for $(r,\theta,\varphi)\in(0,R_{o})\times(0,\pi)\times(0,2\pi)$, where $c(r)$ is the velocity model and $\hat{f}(\omega,r,\theta,\varphi)$ is the source term in the frequency domain. If we consider the situation when this model comes from the acoustic wave equation, then this means that the density is constant and only the bulk modulus is heterogeneous so that we can view $c(r)$ as the independent variable as we do here. We consider two types of sources in this paper:

• the source corresponding to seismology is point-like and located just below the surface:

[TABLE]

where ${f}(t)$ is a short pulse. We will discuss this case in detail in Section 4.

• the source corresponding to daylight imaging is localized in an annulus and emits random noise that is uncorrelated in space and stationary in time:

[TABLE]

where $F(t-t^{\prime})$ is the time correlation function of the noise sources (its Fourier transform is the power spectral density) and $K(r)$ is a smooth function compactly supported in an annulus $r\in[R_{s},R_{o}]$ below the surface. We will discuss this case in detail in Section 5.

We have Neumann (traction-free) boundary condition at the surface $r=R_{o}$, which corresponds to a Robin-type boundary condition:

[TABLE]

In fact, we will consider in the following a modified version (see Eq. (38)) that takes into account boundary dissipation [16, 20] so that when an interior wave hits the surface not all energy is reflected back into the interior.

At the center $r=0$ the condition is that the field should not have any singularity.

We expand the wave field in spherical harmonics as in Eq. (2). The spherical harmonics

[TABLE]

with $P^{m}_{l}$ being the associated Legendre polynomials, form a complete orthonormal set in the space $L^{2}(\sin\theta d\theta d\varphi)$ and they satisfy:

[TABLE]

The mode amplitudes $\hat{p}_{l,m}(\omega,r)$ defined by (3) satisfy the system of second-order ODEs:

[TABLE]

for $r\in(0,R_{o})$, where

[TABLE]

At the surface $R_{o}$ the mode amplitudes satisfy the Robin condition:

[TABLE]

At the center $r=0$ the condition is that the mode amplitude $\hat{p}_{l,m}$ should not have any singularity. This condition can be made more explicit if we assume that the velocity is homogeneous in a small ball with radius $R_{\delta}$ at the center. Then the mode amplitude $\hat{p}_{l,m}$ satisfies

[TABLE]

for $r\in(0,R_{\delta})$. Since the only regular solution of this second-order ODE is the spherical Bessel function, we get that $\hat{p}_{l,m}$ must be equal to j_{l}\big{(}\omega r/c(0)\big{)}, for $r\in(0,R_{\delta})$, up to multiplicative constant. We consider here $\omega>0$. Therefore an explicit boundary condition is that the mode amplitude should satisfy the Robin condition at $r=R_{\delta}$:

[TABLE]

3.2 Radial Wavefield Decomposition

In this subsection we introduce a mode decomposition that turns out to be useful in the forthcoming high-frequency analysis. Let $c_{o}(r)$ be the smooth component of the speed of propagation:

[TABLE]

where $V(r)$ contains the rapidly varying component responsible for scattering. We introduce two linearly independent solutions $A_{l}(\omega,r)$ and $B_{l}(\omega,r)$ of the second-order ODE:

[TABLE]

for $r\in(R_{\delta},R_{o})$. We do not require the solutions to satisfy any boundary condition, but we require these independent solutions to have a Wronskian equal to one:

[TABLE]

for $r\in(R_{\delta},R_{o})$, which is possible by (17) and indeed ensures linear independence. We will determine the choice for solutions $A_{l}(\omega,r)$ and $B_{l}(\omega,r)$ later on.

We define:

[TABLE]

The mode amplitude $\hat{p}_{l,m}$ can then be written in the form

[TABLE]

Differentiating (19) and (20) we get using (17):

[TABLE]

and thus the amplitudes $\hat{\alpha}_{l,m}$ and $\hat{\beta}_{l,m}$ satisfy

[TABLE]

Therefore, using (13) the mode amplitudes $\hat{\alpha}_{l,m}$ and $\hat{\beta}_{l,m}$ satisfy the system of first-order ODEs:

[TABLE]

for $r\in(R_{\delta},R_{o})$.

At the surface $R_{o}$ we find from Eq. (15) that the mode amplitudes satisfy the linear relation:

[TABLE]

(note that this relation is not trivial as $\partial_{r}A_{l}-R_{o}^{-1}A_{l}$ and $\partial_{r}B_{l}-R_{o}^{-1}B_{l}$ cannot both be zero since the Wronskian of $A_{l}$ and $B_{l}$ is one).

At the surface $r=R_{\delta}$ we find from Eq. (16) that the mode amplitudes satisfy the linear relation:

[TABLE]

(note that this relation is again not trivial since the two coefficients within the big parentheses cannot both be zero since the Wronskian of $A_{l}$ and $B_{l}$ is one and since moreover $j_{l}$ and $j_{l}^{\prime}$ cannot both be zero).

3.3 High-frequency Asymptotics

In this subsection we consider a frequency such that the radius of the Earth is much larger than the corresponding wavelength and we carry out high-frequency asymptotic expansions. More precisely, (in the assumed non-dimensionalized coordinates) we assume that the surface radius $R_{o}$ is an order one quantity and that the source wavelength is small. We therefore here consider a high frequency of the form

[TABLE]

and accordingly introduce a WKB type parameterization with respect to spherically propagating waves. For any $l$ we parametrize the solutions $A_{l}^{\varepsilon}$ and $B_{l}^{\varepsilon}$ as

[TABLE]

where $\tau(r,R_{o})$ is the travel time (obtained from the eikonal equation):

[TABLE]

and the amplitudes $\tilde{A}_{l}(r)$ and $\tilde{B}_{l}(r)$ satisfy the frequency-independent transport equations:

[TABLE]

which can be integrated as

[TABLE]

The mode $A_{l}^{\varepsilon}$ is a down-going mode, while $B_{l}^{\varepsilon}$ is an up-going mode. In order to satisfy the condition (18) that the Wronskian is one, we can take:

[TABLE]

Remark 1

The mode profiles are independent of $l$ to leading order in the regime $\varepsilon\to 0$ as long as $l$ is of order one. There are, however, corrective terms in the WKB expansion that depend on $l$. Indeed, by computing higher order terms in the WKB expansion, it is easy to check that the modes can be expanded as

[TABLE]

or equivalently

[TABLE]

with

[TABLE]

It is indeed possible to carry out the forthcoming analysis with these refined expressions, but below we shall continue with the leading order phase $\tau(r,R_{o})$ as this is sufficient to characterize the leading behavior for our quantities of interest, but note here that the leading-order terms may take a different form when $l={\mathcal{O}}(\varepsilon^{-1/2})$.

The rapidly varying component of the speed of propagation may vary at the scale $\varepsilon$, so we denote it by $V^{\varepsilon}(r)$ so as to remember it may depend on $\varepsilon$. Then, substituting the mode expressions (29-30) with (31-32-33) into the system (25-26), the system of first-order ODEs for the mode amplitudes reads:

[TABLE]

for $r\in(R_{\delta},R_{s})$.

The mode amplitudes also satisfy jump and boundary conditions:

At the surface $R_{o}$ the mode amplitudes should satisfy the linear relation:

[TABLE]

which corresponds to a perfect traction-free boundary condition and a perfect reflection condition $\hat{\beta}_{l,m}^{\varepsilon}(\omega,R_{o})/\hat{\alpha}_{l,m}^{\varepsilon}(\omega,R_{o})=1$ at the surface. If we assume that there is some attenuation or dissipation at the surface, then the reflection condition is not equal to one exactly, but to some number $\Gamma_{R_{o}}\in(-1,1)$, and the mode amplitudes satisfy the linear relation:

[TABLE]

More exactly, if we consider the dissipation boundary condition [16, 20]:

[TABLE]

instead of the traction free boundary condition (11), then we get the boundary condition (37) with

[TABLE]

We will look at $\kappa=0$ or $\Gamma_{R_{o}}=1$ as a limiting case in the following.

At the surface $r=R_{\delta}$ the mode amplitudes satisfy the linear relation:

[TABLE]

where we have used the asymptotic formulas for the spherical Bessel functions:

[TABLE]

Note that this condition is independent on the value of $R_{\delta}$ (which can be expected as the mode amplitudes $\hat{\alpha}_{l,m}^{\varepsilon}$ and $\hat{\beta}_{l,m}^{\varepsilon}$ are constant within the homogeneous central small ball) and that we can take the limit $R_{\delta}\to 0$:

[TABLE]

Remark 2

When there is no attenuation, i.e. when $\kappa=0$ and $\Gamma_{R_{o}}=1$, then there is existence and uniqueness of the solution provided $\omega$ is not an eigenfrequency. A frequency $\omega$ is an eigenfrequency if the linear system (34-35) for $r\in(0,R_{o})$ admits a non-zero solution that satisfies the two boundary conditions (37-40) when there is no source. When there is attenuation then there is existence and uniqueness of the solution for any $\omega$, this follows from the form of the propagator associated with (34-35), see Chap. 7 in [9].

Across a source at depth $R_{s}$ whose Fourier component at frequency $\omega/\varepsilon$ has the form

[TABLE]

the mode amplitudes satisfy the jump source conditions:

[TABLE]

It is convenient to introduce the propagator of the system (34-35), that is, the $2\times 2$ matrix ${\bf P}^{\varepsilon}_{l}(\omega,r^{\prime},r)$ solution of

[TABLE]

with ${\bf P}^{\varepsilon}_{l}(\omega,r^{\prime},r=r^{\prime})={\bf I}$. The matrix ${\bf P}^{\varepsilon}_{l}(\omega,r^{\prime},r)$ has the symplectic form

[TABLE]

where $(\hat{a}^{\varepsilon}_{l}(\omega,r^{\prime},r),\hat{b}^{\varepsilon}_{l}(\omega,r^{\prime},r))^{T}$ satisfies

[TABLE]

starting from $\hat{a}^{\varepsilon}_{l}(\omega,r^{\prime},r=r^{\prime})=1$, $\hat{b}^{\varepsilon}_{l}(\omega,r^{\prime},r=r^{\prime})=0$. The solution satisfies the energy conservation relation

[TABLE]

3.4 Wave Decomposition in the Radial Earth

In this subsection we summarize the basic high-frequency wave field decomposition in the stratified sphere. We define the scaled Fourier transform $\hat{f}^{\varepsilon}(\omega)$ of a function $f(t)$ as

[TABLE]

In spherical coordinates $(r,\theta,\varphi)$, the time-harmonic wave field $\hat{p}^{\varepsilon}(\omega,r,\theta,\varphi)$ satisfies

[TABLE]

for $(r,\theta,\varphi)\in(0,R_{o})\times(0,\pi)\times(0,2\pi)$, where $c(r)$ is the velocity model and $\hat{f}^{\varepsilon}(\omega,r,\theta,\varphi)$ is the source term in the frequency domain. We assume here that $\hat{f}^{\varepsilon}$ is supported away from the origin in $\omega$. The radially dependent speed parameter $c(r)$ is modeled by

[TABLE]

The wave field is then decomposed as

[TABLE]

with $Y_{l,m}$ defined in Eq. (12).

For $0<r^{\prime}<r<R_{o}$ so that there is no source in the interval $(r^{\prime},r)$, the mode amplitudes are related by the propagator as

[TABLE]

with the propagator ${\bf P}^{\varepsilon}_{l}$ defined by (43).

At the surface $r=R_{o}$ the mode amplitudes satisfy the boundary condition (37).

At the center $r=0$ the mode amplitudes satisfy the boundary condition (40).

Assume a source at depth $r=R_{s}$ as

[TABLE]

then the mode amplitudes satisfy the jump conditions in Eq. (42).

4 Scattering in Seismology

We consider a source of the form (8) just below the surface $R_{s}=R_{o}^{-}$ whose emission has a typical frequency of order $\varepsilon^{-1}$:

[TABLE]

where $f$ is a function whose Fourier transform is supported away from the origin. Therefore the source term defined as in (47)-(53) reads

[TABLE]

with

[TABLE]

The mode amplitudes satisfy the jump source conditions:

[TABLE]

We sum Eq. (55) multiplied by $-R_{o}^{2}\partial_{r}B_{l}(\omega,R_{o})$ and Eq. (56) multiplied by $R_{o}^{2}B_{l}(\omega,R_{o})$, and we make use of the Wronskian relation (18) to get

[TABLE]

Similarly, we sum Eq. (55) multiplied by $R_{o}^{2}\partial_{r}A_{l}(\omega,R_{o})$ and Eq. (56) multiplied by $-R_{o}^{2}A_{l}(\omega,R_{o})$ and we obtain

[TABLE]

At the surface we “measure” the spherical harmonics of the field:

[TABLE]

Proposition 1

The measured field is of the form

[TABLE]

where we have defined the modulus one function

[TABLE]

Proof. We have from Eq. (51):

[TABLE]

The mode amplitudes $\hat{\alpha}_{l,m}^{\varepsilon}$ and $\hat{\beta}_{l,m}^{\varepsilon}$ satisfy the boundary condition (37) at $r=R_{o}$, the source jump condition (42) at $r=R_{o}^{-}$, the propagation equation (52) from $r=0$ to $r=R_{o}^{-}$ (using the form (44) of the propagator matrix), and the boundary condition (37) at $r=0$:

[TABLE]

By solving this $6\times 6$ linear system for

[TABLE]

we obtain the expression of $(\hat{\alpha}_{l,m}^{\varepsilon}(\omega,R_{o}),\hat{\beta}_{l,m}^{\varepsilon}(\omega,R_{o}))$ that we substitute into (61), which gives the desired result (59). $\Box$

If we assume that the pulse is an even function and we consider the symmetrized response function

[TABLE]

then its Fourier transform (47) is of the form

[TABLE]

which is equal to

[TABLE]

Using |1+\hat{\cal R}_{l}^{\varepsilon}(\omega,R_{o})|^{2}=2{\rm Re}\big{\{}1+\hat{\cal R}_{l}^{\varepsilon}(\omega,R_{o})\big{\}}, we finally establish the main result of this section.

Proposition 2

The symmetrized response function is of the form:

[TABLE]

where $\hat{\cal R}_{l}^{\varepsilon}(\omega,R_{o})$ is defined by (60).

The scattering operator

[TABLE]

is the quantity of interest that characterizes the medium and that can be extracted from standard seismology. Recall that $\hat{\cal R}_{l}^{\varepsilon}(\omega,R_{o})$ has modulus one. Frequencies such that $\hat{\cal R}_{l}^{\varepsilon}(\omega,R_{o})=1$ correspond to the eigenmodes, and we remark that then

[TABLE]

as $\kappa\to 0$ in view of the relation (39).

5 Daylight Imaging

We consider a noise source as in (9) and (10) with a correlation function of the form

[TABLE]

which means that the power spectrum of the noise source contains frequencies of the order of $\varepsilon^{-1}$. Under these conditions the source term $\hat{f}_{l,m}^{\varepsilon}(\omega,r)$ defined by (47) and (53) is a random process with mean zero and covariance function

[TABLE]

The field at the surface is $p(t,R_{o},\theta,\varphi)$ and we compute its empirical autocorrelation components associated with the spherical harmonics:

[TABLE]

for

[TABLE]

In practice, the signals $p_{l,m}(t,R_{o})$ are computed from the signals $(p(t,R_{o},\theta_{j},\varphi_{j}))_{j=1,\ldots,N}$, recorded by a collection of $N$ receivers located at $((R_{o},\theta_{j},\varphi_{j}))_{j=1,\ldots,N}$ at the surface of the Earth via a quadrature formula. We comment in more detail on robustness with respect to sampling in Section 6.1.

As $T\to\infty$ the correlations converges in probability towards the statistical autocorrelation defined by

[TABLE]

In the following proposition we study its scaled Fourier transform:

[TABLE]

Proposition 3

The statistical autocorrelation is of the form

[TABLE]

with

[TABLE]

Proof. We have

[TABLE]

where $\hat{\alpha}_{l,m}^{\varepsilon}(\omega,r;R_{s})$ and $\hat{\beta}_{l,m}^{\varepsilon}(\omega,r;R_{s})$ are the mode amplitudes corresponding to a unit-amplitude spherical harmonic point source at $r=R_{s}$. We get then in view of (66) and (69):

[TABLE]

The mode amplitudes $\hat{\alpha}_{l,m}^{\varepsilon}(\omega,r;R_{s})$ and $\hat{\beta}_{l,m}^{\varepsilon}(\omega,r;R_{s})$ satisfy the boundary, propagation, and source jump conditions:

[TABLE]

By solving this linear system we get the desired result. To get the desired expression (70), we use the energy conservation relation (46):

[TABLE]

and we also make use of the propagation relation ${\bf P}^{\varepsilon}_{l}(\omega,0,R_{o})={\bf P}^{\varepsilon}_{l}(\omega,R_{s},R_{o}){\bf P}^{\varepsilon}_{l}(\omega,0,R_{s})$ which gives

[TABLE]

$\Box$

Although the autocorrelation function (70) is related to the scattering operator of interest (64), it also depends on the source distribution in a non-trivial way, which may complicate the extraction of the scattering operator from the autocorrelation of the measured data, compared to the standard daylight imaging configuration addressed in [5, 6, 11, 17], in which the bottom condition is a radiating condition. However, with arbitrary noise source distributions, the autocorrelation function (70) and the symmetrized response function (63) have the same denominator, so that the identification of the eigenfrequencies can be carried out with both sets of data with the same accuracy. Moreover, in a realistic source configuration the relationship between the autocorrelation function (70) and the scattering operator (64) becomes much simpler. In fact, as we show in Eq. (73), the autocorrelation may produce the scattering operator up to a frequency-dependent modulation function that depends on the temporal spectrum of the noise source trace. This is the case if we assume that the support of the noise sources (that is to say, the support of the function $K$) is localized below the surface, and its thickness is smaller than the typical wavelength. Then for any $R_{s}$ in the support of $K$, we have $\hat{\cal S}^{\varepsilon}_{l}(\omega,R_{s})\simeq 1$ by (71), $\hat{\cal T}^{\varepsilon}_{l}(\omega,R_{s})\simeq\hat{\cal R}^{\varepsilon}_{l}(\omega,R_{o})$ by (60) and (72), and $\tau(R_{s},R_{o})\simeq 0$, so the square brackets in (70) can be simplified as

[TABLE]

As a consequence we get the following result.

Proposition 4

If the spatial support of the noise source is localized in a small annulus below the surface, then the statistical autocorrelation is of the form

[TABLE]

This result shows that the autocorrelation of the noise signals (70) is directly related to the scattering operator (64). More precisely, the time derivative of the autocorrelation function is proportional to the symmetrized response function (62), provided $\hat{F}(\omega)=\hat{f}(\omega)$ (which imposes in particular that the pulse profile $f$ should be an even function):

[TABLE]

We recover the classical seismic interferometry formula that has been established in many other configurations [11, 18, 21].

6 On Robustness

6.1 Robustness with Respect to Sampling Configuration

In the previous Section 5 we studied the autocorrelation function $C_{l,m}(t,R_{o})$. Note that this is the autocorrelation function of the $(l,m)$-th component associated with the spherical harmonics, so that it is a linear transform of the cross correlation of the field observed at the surface:

[TABLE]

From this it may seem that we need to observe the field everywhere at the surface to estimate this quantity (by further substituting a time average for the expectation). However, we do not need such extensive data and we can relax the hypothesis that we observe the field everywhere at the surface of the Earth. Indeed we have shown that $C_{l,m}(t,R_{o})$ does not depend on $m$. This is due to the spherical symmetry of the Earth model and the noise source distribution. Therefore we can consider

[TABLE]

where actually all terms are equal and given by the expression in Proposition 3. Using the addition theorem of spherical harmonics, we have

[TABLE]

where $P_{l}$ is the Legendre polynomial and $\cos\Omega$ is the angle between two unit vectors oriented at the polar coordinates $(\theta,\varphi)$ and $(\theta^{\prime},\varphi^{\prime})$. The expectation

[TABLE]

depends only on $\cos\Omega$ and it can be expanded as

[TABLE]

The quantity of interest $C_{l}(t,R_{o})$ is given by

[TABLE]

This shows that observation points at the surface of the Earth that are such that the angles $\Omega$ between pairs of points cover the interval $(0,\pi)$ are sufficient to estimate the Earth spectrum (and the scattering operator under the additional hypothesis in Proposition 4): this data set gives estimates of ${\cal C}(t,\Omega)$ for a sufficiently dense grid of $\Omega$; then the quantity of interest $C_{l}(t,R_{o})$ can be obtained by a numerical evaluation (a quadrature formula) of the integral (75); the Earth point spectrum can then be obtained by inspection of the Fourier transform $\hat{C}_{l}(\omega,R_{o})$ whose peaks correspond to eigenfrequencies.

6.2 Robustness with Respect to Source Configuration

Consider the case when the assumption in (6) is generalized as

[TABLE]

thus allowing for general lateral source distribution with density function $G(\theta,\varphi)$. Then

[TABLE]

with

[TABLE]

In this general setting we find that Propositions 3 and 4 still hold true but $\hat{C}_{l,m}^{\varepsilon}(\omega,R_{o})$ is now multiplied by $G_{l,m}$ compared to Eqs. (70) and (73). Since these new expressions depend on $m$ only though $G_{l,m}$, we can sum over $m=-l,\ldots,l$ to get rid of the dependence with respect to $m$ of the source distribution, and there remains only a dependence with respect to $l$ via a positive multiplicative factor $\sum_{m=-l}^{l}G_{l,m}$. This shows that, as long as the main objective is to extract the eigenfrequencies, we have robustness with respect to the assumption of angular homogeneity in the random source distribution. However, the quantitative estimation of the Earth scattering operator is sensitive to the angular source distribution through the multiplicative factor $\sum_{m=-l}^{l}G_{l,m}$ that affects the amplitude of the estimation.

6.3 Robutsness with Respect to Medium Noise or Model Error

We have assumed in the previous sections that the speed of propagation has only radial variations. Here we want to show that the result is robust with respect to certain angular variations of the velocity model. More exactly, we consider an Earth-model with slow angular variations, whose velocity model has the form:

[TABLE]

in which there are small and slow angular velocity fluctuations $V_{2}^{\varepsilon}$. We call them slow because they vary at the scale one with respect to the angles $\theta$ and $\varphi$, while they may vary slowly or rapidly in the radial coordinate $r$. We want to clarify under which circumstances the perturbation $V_{2}^{\varepsilon}$ can be neglected.

The field has the form (50-51) and the coupled system of first-order ODEs for the mode amplitudes now reads:

[TABLE]

for $r\in(0,R_{o})$, instead of (34-35), with

[TABLE]

We consider the case in which the fluctuation term $V_{2}^{\varepsilon}$ may have slow and fast components:

[TABLE]

where $V_{21}$ and $V_{22}$ can be arbitrary functions and $V_{23}$ is a zero-mean random process, that is stationary and mixing in $r$.

We have the following results:

• Provided $a>1$, the term $V_{21}$ gives rise to terms in (77-78) that vanish in the limit $\varepsilon\to 0$.

Provided $b>1$, the term $V_{22}$ gives rise to terms in (77-78) that vanish in the limit $\varepsilon\to 0$.

Provided $c>0$, the term $V_{23}$ gives rise to terms in (77-78) that vanish in the limit $\varepsilon\to 0$.

The first two assertions are trivial by direct inspection of the amplitude of the coupling terms $\varepsilon^{a-1}q(r)$ and $\varepsilon^{b-1}q(r/\varepsilon)$, respectively, but the third assertion is not so trivial as it involves coupling terms of the form $\varepsilon^{c-1}q(r/\varepsilon^{2})$ for a mixing and zero-mean process $q$. However, diffusion approximation theory reveals that this coupling becomes effective only when $c=0$ [9]. Therefore all crossed terms in the equations (77-78) cancel and we get a system similar to (34-35). The energy conservation equation (46), which is the key property, holds true and the conclusion of the previous section holds including robustness with respect to partial source and measurement aperture.

6.4 Robustness with Respect to Measurement Noise

The results presented in the previous sections are based on the behavior of the statistical autocorrelation function $C_{l,m}(t,R_{o})$. In this section we want to explain that this function can be estimated by the empirical autocorrelation of the recorded signals and that it is quite robust with respect to additive measurement noise.

The statistical autocorrelation function $C_{l,m}(t,R_{o})$ defined by (69) is the one that is related to the scattering operator. The empirical autocorrelation function $C_{l,m}^{T}(t,R_{o})$ defined by (67) is the one that is computed from the recorded data. It is easy to check that the expectation of $C_{l,m}^{T}(t,R_{o})$ is exactly $C_{l,m}(t,R_{o})$ for any $T$. By a detailed fluctuation analysis it is also possible to show that the variance of $C_{l,m}^{T}(t,R_{o})$ is proportional to $1/T$ [11, Chapter 2]. This ensures that the empirical autocorrelation function is a good estimate of the statistical autocorrelation function provided the recording time window $T$ is large enough. As noticed in [11, Chapter 2], the required recording time is all the longer as lower frequency components are being investigated.

When the recorded signals are polluted by additive measurement noise, which is independent from the noise sources, the empirical autocorrelation function in Eq. (67), computed via a discretization of Eq. (68), is the sum of the statistical autocorrelation function $C_{l,m}(t,R_{o})$ and of the autocorrelation function $C_{l,m}^{\rm n}(t)$ that comes from the additive noise at the $N$ receivers. Measurement noise is independent from one receiver to the other one, so that, if the $N$ receivers are distributed uniformly at the surface of the Earth and the measurement noise statistics is the same at each receiver we have:

[TABLE]

where $F^{\rm n}(t)$ is the autocorrelation function of the measurement noise ($\hat{F}^{\rm n}(\omega)$ is the power spectral density of the measurement noise). This shows that the impact of the measurement noise decays with the number of receivers. Moreover, measurement noise usually has higher frequencies than the spectral band that is of interest for global seismology. However, if the spectrum of the additive noise intersects the spectral band over which we look for eigenfrequencies and/or the scattering operator, then measurement noise may corrupt the estimation. The estimation of the eigenfrequencies should be quite robust as measurement noise is very unlikely to have spectral peaks similar to the ones that are investigated. The estimation of the scattering operator (which is more sensitive to the estimated amplitude) may be affected.

7 Concluding Remarks

We have shown that global acoustic daylight imaging is possible. The point spectrum of the Earth can be extracted from the correlation functions of the signals recorded at the surface and emitted by unknown noise sources localized away from the surface. Under an additional realistic assumption of the spatial support of the noise sources, the complete scattering operator can be extracted from the correlation functions. The first result (on the extraction of the point spectrum) is robust to receiver distribution, source distribution, and medium and measurement noise, but one should be careful to interpret the second result when the noise sources are not evenly distributed as this may perturb the amplitudes of the recovered signals.

In this paper we have only considered low angular modes (with Legendre number $l$ of the order of one). Indeed we do not need to analyze the high angular modes because there is no coupling between low and high angular modes in our setting. We could consider high angular modes, but then a non-trapping condition would be necessary, such as the Herglotz [14] and Wiechert and Zoeppritz [23] conditions.

We have here, in Eq. (10), considered the situation that the noise sources are delta-correlated in space. This assumption simplifies the analysis but can be relaxed provided the correlation radius is small by using stationary phase methods [11] or semi-classical analysis [7].

8 Acknowledgement

Maarten V. de Hoop was supported in part by the Simons Foundation under the MATH + X program. Knut Sølna was supported in part by AFOSR under grant FA9550-14-1-0197 and NSF under grant 1616954.