On the stability of some imaging functionals

This work is devoted to the stability/resolution analysis of several imaging functionals in complex environments. We consider both linear functionals in the wavefield as well as quadratic functionals based on wavefield correlations. Using simplified measurement settings and reduced functionals that retain the main features of functionals used in practice, we obtain optimal asymptotic estimates of the signal-to-noise ratios depending on the main physical parameters of the problem. We consider random media with possibly long-range dependence and with a correlation length that is less than or equal to the central wavelength of the source we aim to reconstruct. This corresponds to the wave propagation regimes of radiative transfer or homogenization.


Introduction
Imaging in complex media has a long history with many applications such as non-destructive testing, seismology, or underwater acoustics [16,28,3]. Standard methodologies consist in emitting a pulse in the heterogeneous medium and performing measurements of the wavefield or other relevant quantities at an array of detectors. Depending on the experimental setting, this may give access, for instance, to the backscattered wavefield (and its spatio-temporal correlations), or to the wave energy. The imaging procedure then typically amounts to backpropagating the measured data appropriately. When scattering of the wave by the medium is not too strong, measurements are usually migrated in a homogeneous medium neglecting the heterogeneities. This is the principle of the Kirchhoff migration and similar techniques, and is referred to as coherent imaging [11]. When scattering is stronger, the interaction between the wave and the medium cannot be neglected and a different model for the inversion is needed. We will only consider weak heterogeneities in this work, so that a homogeneous model is enough to obtain accurate reconstructions. We refer to [10,8] for a consideration of stronger fluctuations and transport-based imaging.
We are interested here in comparing two classes of imaging methods in terms of stability and resolution: one based on the wavefield measurements (and thus linear in the wavefields), such as Kirchhoff migration; and the other one based on correlations of the wavefield (and therefore quadratic in the wavefields), such as coherent interferometry imaging [14]. Here, "stability" refers to the stability of the reconstructions with respect to changes in the medium or to a measurement noise.
Stability and resolution of coherent imaging functionals in random waveguides were addressed in [13]; for functionals based on topological derivatives, see also e.g. [2]. We refer to [19,20,21,12] and references therein for more details on wave propagation in random media and imaging. Our analysis is done in the framework of three-dimensional acoustic wave propagation in a random medium with correlation length c and fluctuations amplitude of order σ 0 . The quantities c and σ 0 are fixed parameters of the experiments. The random medium is allowed to exhibit either short-range or long-range dependence. Our goal is to image a source centered at the origin from measurements performed over a detector array D. Denoting by λ the central wavelength of the source, we will assume that it is larger than or equal to the correlation length, and that many wavelengths separate the source from the detector, so that we are in a high frequency regime. The case λ = c is referred to as the weak coupling regime in the literature (for an amplitude σ 0 1), while the case λ c is the stochastic homogenization regime. The ratio c /λ is a critical parameter as it controls, along with σ 0 , the strength of the interaction of the sound waves with the heterogeneities. The regime λ c is addressed in [1] and is known as the random geometrical optics regime. The stability/resolution analysis is somewhat direct in that case, since a simple expression for the heterogeneous Green's function is available by means of random travel times. This is not the case in our regimes of interest, where the interaction between the waves and the underlying medium is more difficult to describe mathematically.
We will work with a simplified configuration and will define reduced imaging functionals that retain the main features of the functionals used in practice while offering more tractable computations. In such a setting, we obtain optimal estimates for the stability of the functional in terms of the most relevant physical parameters. We furthermore quantify the signal-to-noise ratio (SNR). The statistical stability is evaluated by computing the variance of the functionals at the source location. For the wave-based functional (WB), this involves the use of stationary phase techniques for oscillatory integrals, while computations related to the correlation-based functional (CB) involve averaging the scintillation function of Wigner transforms [23,25] against singular test functions. It is now well-established that correlation-based functionals enjoy a better stability than wave-based functionals at the price of a lower resolution [14]. This is due to self-averaging effects of Wigner transforms that we want to quantify here.
The paper is organized as follows: in Section 2, we present the setting and our main results. We define the wave-based and correlation-based functionals and describe our models for the measurements in Section 3. Section 4 is devoted to the analysis of the resolution and stability of the correlation-based functional and Section 5 to that of the wave-based functional. The proofs of the main results are presented in Section 6 and a conclusion is offered in Section 7.

Wave propagation and measurement setting
The propagation of three-dimensional acoustic waves is described by the scalar wave equation for the pressure p: x ∈ R 3 , t > 0, supplemented with initial conditions p(t = 0, x) and p t (t = 0, x). We suppose for simplicity that the density is constant and equal to one; ρ = 1. We also assume that the compressibility is random and takes the form where σ 0 measures the amplitude of the random fluctuations and c their correlation length. We suppose that V is bounded and σ 0 is small enough so that κ remains positive. The sound speed c = (κρ) − 1 2 thus satisfies c = 1 + O(σ 0 ) 1 and the average sound speed is c 0 = 1. The random field V is a mean-zero stationary process with correlation function where E denotes the ensemble average over realizations of the random medium. We assume R to be isotropic so that R(x) = R(|x|). We will consider two types of correlation functions: (i) integrable R, which models random media with short-range correlations; and (ii) nonintegrable R corresponding to media with long-range correlations. The latter media are of interest for instance when waves propagate through a turbulent atmosphere or the earth upper crust [17,29]. Such properties can be translated into the power spectrumR, the Fourier transform of R, by supposing thatR has the form where S is a smooth function with fast decay at infinity. A simple dimensional analysis shows that R(x) behaves likes |x| δ−3 , δ < 3, as |x| → ∞, which is not integrable for large |x|. The case δ = 0 corresponds to integrable R since the function S is regular. We assume that δ < 2 so that the transport mean free path (see section 3.3) is well-defined. Propagation in media such that 2 < δ < 3 is still possible but requires an elaborate theory of transport equations with highly singular scattering operators that are beyond the scope of this paper. Our setting of measurements is the following: we assume that measurements of the wavefield are performed on a three dimensional detector D and at a fixed time T . We assume that the detector D is a cube centered at x D = (L, 0, 0) of side l < L; see figure 2.1. Our goal is then to image an initial condition localized around x 0 = (0, 0, 0). Practical settings of measurements usually involve recording in time of the wavefield on a surface detector. The two configurations share similar three-dimensional information: spatial 3D for our setting, and spatial 2D + time 1D for the standard configuration. Using wave propagation in a homogeneous medium, it is also relatively straightforward to pass from one measurement setting to the other. Our choice of a 3D detector was made because it offers slightly more tractable computations for the stability of the imaging functionals while qualitatively preserving the same structure of data.
As waves approximately propagate in a homogeneous medium with constant speed c 0 , it takes an average time t D = (L − l/2)/c 0 for the wave to reach the detector. We will therefore suppose that T > t D , and even make the assumption that T = L/c 0 so that the wave packet has reached the center of the detector. Note that in such a measurement configuration, the range of the source is then known since the initial time of emission is given. We therefore focus only on the cross-range resolution.
We consider an isotropic initial condition obtained by Fourier transform of a frequency profile g (a smooth function that decays rapidly at infinity): where |k 0 | is a non-dimensional parameter, λ is the central wavelength, and Bc −1 0 the bandwidth.
Rescaling variables as x → xL, t → tL/c 0 , introducing the parameters η = c /λ, ε = λ/L, and still denoting by p the corresponding rescaled wavefield, the dimensionless wave equation now reads We quantify the bandwidth of the source in terms of the central frequency by setting B = (Lεµ) −1 , where µ is a given non-dimensional parameter such that µ 1. We suppose that µ 1/ √ ε so that the initial condition can be considered as broadband. More details about the latter condition will be given later. The rescaled initial condition then has the form The normalization is chosen so that p ε 0 (0) is of order O(1). Ifĝ denotes the Fourier transform of g, we can write Above, the symbol means equality up to negligible terms and S 2 is the unit sphere of R 3 . The initial condition is therefore essentially a function with support εµ oscillating isotropically at a frequency |k 0 |/ε. As explained in the introduction, we assume that c ≤ λ, which implies η ≤ 1. The case η ∼ 1 leads to the radiative transfer regime in the limit ε → 0 [27]. The case η 1 corresponds to the homogenization regime [5] and a propagation in an effective medium (which is here the homogeneous medium since σ 0 is small). The opposite case η 1 gives rise to the random geometric regime addressed in [1].
For the asymptotic analysis, we recast the scalar wave equation (2.1) as a first-order hyperbolic system on the wavefield u = (v, p), where v is the velocity: augmented with initial conditions p(t = 0, x) = p ε 0 (x) and v(t = 0, x) = 0. The latter system is rewritten as where V = diag(0, V ) ∈ R 4×4 , D j mn = δ m4 δ nj +δ n4 δ mj is a 4×4 symmetric matrix for 1 ≤ j ≤ 3. Above and in the sequel, we use the Einstein convention of summation over repeated indices.

Results
Our main results concern the signal-to-noise ratio at the point x = 0 defined by where I C stands for the CB functional and I W for the WB functional, which are defined further in section 3. We will also quantify the support of Var{I}(z). We will obtain optimal estimates for the CB and WB functionals in terms of the main parameters that define them as well as ε, η, δ and µ. Our measurements of the pressure are assumed to have the form x) models the statistical instabilities in the Born approximation and σ n n ε p (t, x) is an additive mean zero noise with amplitude σ n . Note that by construction, coherentbased functionals will perform well only in the presence of relatively weak heterogeneities. Hence, modeling the statistical instabilities at first order is both practically relevant and mathematically feasible. Taking into account second order interactions is considerably more difficult mathematically; see [9] for a stability analysis of Wigner transforms (and not of the more complicated imaging functionals) in the paraxial approximation. The term σ n n ε p models an additive noise at the detector and also takes into account in a very crude way the higher order interactions between the wave and the medium that are not included in δp ε . We suppose that n ε p oscillates at a frequency ε −1 (so that a simple frequency analysis cannot separate the real signal from the noise) and that it is independent of the random medium. The variance can then be decomposed as Var{I C } = V C + V C n , where V C is the variance corresponding to the single scattering term and V C n the noise contribution. We also write Var{I W } = V W + V W n . Note that the measurement noise can actually be much larger than the average E{p}. Indeed, a simple analysis of E{p}(t, x) shows that it is of order ε when |x| = 1 (omitting the absorption). It is the refocusing properties of the functionals that lead to a reconstructed source of order one from measurements of order ε.
We formally rescale the single scattering instabilities by e −c 0 Σt/2 , where Σ −1 is the transport mean free path defined in section 3.3, so that the first-order interaction between the average field E{p} (proportional to e −c 0 Σt/2 ) and the medium has a comparable amplitude to E{p}. The fact that the single scattering instabilities are exponentially decreasing as the ballistic part can be proved in simplified regimes of propagation, where a closed-form equation for their variance can be obtained; see e.g. [4,7].
We need to introduce another important parameter, which is the typical length over which correlations are calculated in the CB functional. In dimensionless units, we define it in terms of the wavelength by N C ε. We assume that the detector cannot perform subwavelength measurements, which implies N C ≥ 1. We suppose for simplicity that correlations are calculated isotropically. Accounting for anisotropic correlations would add additional parameters and technicalities, but presents no conceptual difficulties. Since the resolution and the stability of the CB functional are mostly influenced by N C ε and not by the size of the detector as N C ε ≤ l/L by construction (see section 4.1 below), we will systematically suppose that the detector is sufficiently large compared to the wavelength so that its effects on the stability can be neglected in a first approximation. The parameter N C is crucial in that it controls the resolution of the CB functional (shown in section 4.1 to be L/N C ) and its stability. Small values of N C yield a good stability for a poor resolution, while large values lead to less stability with a resolution comparable to that of the WB functional.  Our main results are the following: Theorem 2.1 Denote by V C (z) (resp. V C n ) and V W (z) (resp. V W n ) the variances of the correlation-based and wave-based functionals defined in section 3 for the single scattering contribution (resp. the noise contribution). Let Σ −1 be the transport mean free path defined in (3.9) in section 3.3 below. Then we have: Moreover, V C , V C n and V W are mostly supported on |z| ≤ µε. Above, the notation ∼ means equality up to multiplicative constants independent of the main parameters, and a ∧ b = min(a, b), a ∨ b = max(a, b) . Denote now by SN R C (resp. SN R C n ) and SN R W (resp. SN R W n ) the corresponding signal-to-noise ratios. Then: We below comment on the results of the theorem.
Comparison CB-WB. Let us consider first the case of short-range correlations δ = 0. We neglect absorption at this point (i.e. e −ΣL = 1) and will take it into account when considering the SNR of the external noise. Let us start with the single scattering contributions SN R C and SN R W . Suppose first that N C is small enough and µ large enough so that (L/N C λ) ∧ µ = µ. Then, Hence, the SNR increases when the following quantities increase: λ/ c (the dynamics gets closer to the homogenization regime where the homogenized solution is deterministic), µ (smaller bandwidth, which results in a loss of range resolution) and σ −1 0 (weaker fluctuations). The fact that the SNR increases as µ stems from the self-averaging effects of the quadratic functional. There are no such effects for V W and we find SN R C (0) = SN R W (0)(µλ/ c ) so that the CB functional is more stable than the WB functional in this configuration of small N C , whether in the radiative transfer regime (λ = c ) or in the homogenization regime (λ c ). It is shown in section 4.1 that the resolution of the CB functional is L/N C , so that the best resolution is achieved for the largest possible value N C , where correlations are calculated over the largest possible domain, namely N C = l/λ. Hence, the best resolution is λL/l, which is the celebrated Rayleigh formula and the same cross-range resolution as the WB functional. Now, for this best resolution and when lµ > L, then (L/N C λ) ∧ µ = (L/N C λ) so that SN R C (0) = SN R W (0)(λ/ c ) (equality here means equality up to multiplicative constants independent of the parameters λ, N C , c , σ 0 ). This means that in the absence of external noise, and in a weak disorder regime where multiple scattering can be neglected, the SNR of the CB and WB functionals differ only by the factor (λ/ c ) for a similar resolution. This is a term of order one in the radiative transfer regime, and a large term in the homogenization regime. In the radiative transfer regime, in order to significantly increase the CB SNR compared to WB, one needs to decrease N C and therefore to lower the resolution. This is the classical stability/resolution trade off as was already observed in [1] for the random geometrical regime. Note also that the statistical errors for both functionals are essentially localized on the support (of diameter µλ) of the initial source.
We consider now the noise contributions V C n and V W n . A first important observation is that V C n (x) is mostly localized on the support of the source, while V W n (x) is essentially a constant everywhere. This stems from the fact that in the CB functional, the noise is correlated with the average field so that the functional keeps track of the source direction, while in the WB functional the noise is correlated with itself, see sections 3.3.2 and 5. Moreover,the SNR are As for the single scattering SNR, for identical resolutions (i.e. when N c = l/λ), the SNR are comparable. In order to obtain to better SNR for CB, one needs to lower the resolution (and therefore decrease N C ) so that (L/(N C λ))∧µ = µ and SN R C Minimal wavelength for a given SNR. We want to address here the question of the lowest central wavelength λ m of the source (and therefore the best resolution) that can be achieved for a given SNR. In order to do so, we need to take into account the frequency-dependent absorption factor. We first consider the WB functional. We assume for concreteness that the noise contribution is larger than the single scattering one, but the same analysis holds for the reversed situation. Let us fix the SNR at one. Hence, λ m has to be such that e −ΣL/2 = σ n . Writing σ n = e −τn , this yields ΣL = 2τ n . We will see in section 3.3.1 that Σ admits the following expression .
Hence, as can be expected, λ m decreases as the fluctuations of the random medium and the intensity of the noise decrease (i.e. σ 0 ↓ and τ n ↑.) Remark as well that λ m decreases as c does, and in a faster way compared to σ 0 or τ n . This is also expected since the limit c → 0 corresponds to the homogenization limit in which the wave propagates in a homogeneous medium provided σ 0 is small. The measurements are therefore primarily coherent and both the CB and WB functionals perform well. Let us consider now the CB functional and let δ = 0. Since for (L/(N C λ)) ∧ µ 1, SN R C is greater than SN R W , we may expect in principle to find a lower minimal wavelength. A way to exploit this fact is to consider a central wavelength of the source λ m /α, for the λ m above and some α > 1, and to compute correlations over a sufficiently small domain in order to gain stability, but a domain not too large so that the resolution is still better than λ m . The resolution of the CB functional being L/N C , we choose N C = Lβ/λ m , with β > 1 so that λ m /β is smaller than λ m . The prefactor in the SNR is then the square of (Lα/(N C λ m )) ∧ µ = (α/β) ∧ µ that we suppose is equal to (α/β) 2 . We find and we look for α > 1 and β > 1 such that σ 4α−1 n α/β = 1. Since σ n ≤ 1, this is possible only when α is not too large (so that the central wavelength of the source cannot be too small, otherwise absorption is too strong) and when σ n is not too small either. When σ n is below the threshold, then the minimal wavelength is the same for the CB and WB functionals. Hence, compared to the WB functional, the averaging effects of the CB functional can be exploited in order to improve the optimal resolution only when the noise is significant. Finding the optimal N C is actually a difficult problem, and is addressed numerically in [15].
Effects of long-range correlations. The strongest such effect is seen in Σ since Σ → ∞ as δ → 2. Hence, as the fluctuations get correlated over a larger and larger spatial range, the mean free path decreases and the amplitude of the signal becomes very small. This means that coherent-based functionals are therefore not efficient in random media with long-range correlations, and inversion methodologies based on transport equations [10,8] that rely on uncoherent information should be used preferrably.
There is also a loss of stability in the variance of the single scattering term for the CB functional, which for sufficiently long correlations (i.e. δ > 1) becomes larger than in the short-range case. In such a situation, for the SNR of the CB functional to be larger than that of the WB functional, one needs to decrease the resolution by a factor greater than in the short-range case. Notice moreover that there is an effect of the long-range dependence on the term ( c /λ) 3−δ that measures the distance to the homogenization regime. As the medium becomes correlated over larger distances, the variance increases, which suggests that the homogenization regime becomes less accurate.

Imaging functionals and models
We introduce in sections 3.1 and 3.2 the WB and CB functionals, and in section 3.3 our models for the measurements.

Expression of the WB functional
We give here the expression of the WB functional for generic solutions of the wave equation, which we denote by p and its time derivative ∂ t p. The solution to the homogeneous wave equation with regular initial conditions (p(t = 0) = q 0 , ∂ t p(t = 0) = q 1 ), for (q 0 , q 1 ) given, reads formally in three dimensions where * denotes convolution in the spatial variables and δ 0 the Dirac measure at zero. In the Fourier space, this reduces to where we define the Fourier transform as Fp(k) = R 3 e −ik·x p(x)dx. From the data (p(T, x), ∂ t p(T, x)) at a time t = T for x ∈ D, the natural expression of the WB functional in our setting is obtained by backpropagating the measurements, similarly to the time reversal procedure [18]. This leads to the functional where 1 D denotes the characteristic function of the detector. Using Fourier transforms, this can be recast as This expression is slightly different from the classical Kirchhoff migration functional because of our different measurement setting. It nevertheless performs the same operation of backpropagation. Assuming that all the wavefield is measured, i.e. D = R 3 , and that q 1 = 0 so Fp(t, k) = cos |k|t Fq 0 (k), one recovers the initial condition perfectly, i.e. I W 0 (x) = q 0 (x). In practice, the entire wavefield is generally not available and diffraction effects limit the resolution. In order to simplify the calculations, we assume without loss of generality that only the pressure p(T, x) is used for this functional. This modifies I W 0 as This significantly reduces the technicalities of our derivations while very little affecting the reconstructions. Indeed, for localized initial conditions, the value of the maximal peak of the functional I W 0 is divided by a factor two compared to the full I W 0 : if D = R 3 , and Fp(t, k) = cos |k|t Fq 0 (k), then The first term above yields 1 2 q 0 (x) while the second one is essentially supported on a sphere of radius 2T far away from the source. We will analyse in section 5 the expected value and the variance of the functional I W for random measured wavefields.

Expression of the CB functional
We define here the CB functional for a measured wavefield u. This requires us to introduce the Wigner transform of u [25,23] and to decompose it into propagating and vortical modes. Since the initial velocity is identically zero, the amplitude of the latter modes remains zero at all times, see [27]. The projection onto the propagating mode is done with the eigenvectors of the dispersion matrix L(k) = k i D i where D i was defined in (2.3). These vectors are given by b ± (k) = (k, ±1)/ √ 2, withk = k/|k| and we define in addition the matrices B ± = b ± ⊗ b ± , where ⊗ denotes tensor product of vectors. The full matrix-valued Wigner transform of u is defined by The quantity W ε is a real-valued matrix. In our setting, the field u is measured at a time t = T = 1. We compute the correlations over a domain C ⊂ D, so that we do not form the full Wigner transform but only a smoothed version of it given by We only consider the propagating mode associated with the positive speed of propagation c 0 , the other one can be recovered by symmetry. This mode corresponds to the vector b + and the associated amplitude is given by where (W ε D ) T denotes the matrix transpose of W ε S and Tr the matrix trace. The CB functional is then defined in our setting by Above, we suppose implicitely that if (x + c 0 Tk) ∈ D, then C is such that (x + c 0 Tk) ± εy 2 ∈ D. The functional is slightly different from the classical coherent interferometry functional of [14] because of our measurement setting. However, the two functionals qualitatively perform the same operation, that is the backpropagation of field-correlations along rays. The difference resides in how these correlations are calculated: in our configuration, we have access to 3D volumic measurements at a fixed time so that the 3D spatial Wigner transform is available and is the retropropagated quantity; in [14], measurements are performed on a 2D surface and recorded in time so that the 3D spatial Wigner transform is replaced by a 2D spatio-1D temporal Wigner transform.
For the stability and resolution analysis, it is convenient to recast I C in terms of the amplitude a associated with the full Wigner transform (i.e. computed for C = R 3 ). We then obtain the expression, using the rescaled variables, c 0 = T = 1: Note that I C is real-valued.

Models for the measurements
We introduce in this section our different models for the measurements. We will define a model for the mean of the measurements, and a model for the statistical instabilities. The latter is obtained by using the single scattering approximation (Born approximation). We start with the mean.

Model for the mean
The wavefield. Denoting by = ∂ 2 t 2 − ∆ the d'Alembert operator and by the ballistic part associated to an initial condition q 0 (with vanishing initial ∂ t p), the solution to (2.1) reads: Obtaining an expression for the expectation of p requires the analysis of the term EV ((·/(ηε))∂ 2 t 2 p which is not straighforward and requires a diagrammatic expansion. Rather, we follow the simpler, heuristic, method of [24], which amounts to iterating the above inversion procedure one more time and getting and to replace p in the last term by p B . Since E{p} p B at first order in σ 0 , p B may be replaced by E{p} to obtain a homogenized equation for E{p}. The result is [24] that a harmonic wavep(ω, x) (Fourier transform in time of p) is damped exponentially, i.e. for |x| = 0, The absorption coefficient γ is obtained by adapting the setting of [24] to ours. Since our initial condition (2.2) is mostly localized around wavenumbers with frequency |k 0 |/ε, we find that the waves are absorbed by a factor A Taylor expansion then leads to the classical |k| 4 dependency of the absorption associated to the Rayleigh scattering: We implicitely assumed above thatR(0) was defined, which holds in random media with shortrange correlations but not in media with long-range correlations. The latter case is addressed below. We can relate the latter expression for γ ε to the mean free path Σ −1 := Σ −1 (|k 0 |) defined in the next paragraph in (3.9) by Hence Σ 2γ(|k 0 |/ε).
Obtaining such a relation in the case of long-range correlations is slightly more involved. To do so, we first notice that where the constant c δ is given by [22], Γ being the Gamma function. The expression of γ ε is then Recall that δ < 2 so that the above integral is well-defined. The inverse of the mean free path is now in the long-range case .
We used again the fact that δ < 2 to make sense of the integral. Since the following relation holds , we recover that Σ 2γ(|k 0 |/ε) in the long-range case. We will therefore systematically replace γ(|k 0 |/ε) by Σ/2 in the sequel. The expression of E{p}(t, x) is obtained by Fourier-transforming back E{p}(ω, x) and using the fact that |x| = c 0 t since the ballistic part is supported on the sphere of radius |x| = c 0 t. This yields This is our model for the average pressure.
The Wigner transform. It is shown in [27] that the Wigner transform W ε of a wavefield u satisfies the following system with Noticing that we can recast the r.h.s of (3.7) as It is then well established [27] that a good approximation of E{W ε } is where the matrices B ± were introduced in section 3.2 and the amplitudesā ε ± (t, x, k) satisfy the following radiative transfer equation Above, δ 0 is the Dirac measure at zero. The equation (3.8) is supplemented with the initial conditionā ± (t = 0) = a ε 0,± = Tr(W ε,0 ) T B ± , where W ε,0 is the Wigner transform of the initial wavefield u(t = 0). Sinceā ε − (k) =ā ε + (−k), we focus only on the modeā ε + and drop both the + lower script and the ε upper script for notational simplicity. We now need to compute the Wigner transform of the initial condition.
We consider an approximate expression that simplifies the analysis of the imaging functionals. The scalar Wigner transform of p ε 0 , denoted by w ε 0 , reads: Since the bandwidth parameter µ is such that µε ε, we can separate the scales of theĝ terms and the oscillating exponentials. We can then state that, as is classical for Wigner transforms, different wavevectors lead in a weak sense to negligible contributions because of the highly oscillating term e i|k 0 |x·(k 1 −k 2 )/ε . The leading term in the initial condition is therefore The effect of the function w ε 0 in the phase space is essentially to localize the variable x around εµ and the variable |k| on a shell of radius |k 0 | with width µ −1 . This is what is expected since the initial condition for the wave equation is isotropic, oscillates at a frequency ε −1 |k 0 | and has a bandwidth of order (εµ) −1 . The fact that µ 1/ √ ε implies that the localization is stronger in the spatial variables than in the momentum variables, which can be seen as a broadband property. For the sake of simplicity, we replace the initial condition by a function that shares the same properties but has an easier expression to handle. We write: where the function w 0 is smooth. With c 0 = 1, we then recast (3.8) as ∂ā ∂t +k · ∇ xā + Σ(k)ā = Q + (ā),ā ± (t = 0, x, k) = µw 0 x εµ , µ(|k| − |k 0 |) (3.9) (η(k − |k|p))dp.
When σ 2 0 = ε, η = 1, we recover the usual equation for the weak coupling regime. Note that since R(x) = R(|x|), we haveR(k) =R(|k|) and Σ(k) = Σ(|k|). Going back to the measured amplitude a D defined in (3.2), an accurate approximation of its mean is therefore where the filter F x ε was defined in (3.3). The integral solution to (3.9) reads, denoting by T tā (x, k) :=ā(x − c 0 tk, k) the free transport semigroup: The averageā(t, x, k) is the sum of a ballistic term and the multiple scattering contribution.

Model for the random fluctuations in the measurements
The wavefield. We follow the Born approximation, which consists in retaining in (3.5) the terms at most linear in V . In order to take into account the absorption as explained in section 2.2, we replace the ballistic term in (3.5) by E{p}. In doing so, both the average E{p} and the fluctuations are exponentially decreasing. The random instabilities on the pressure are then generated by the term We suppose that the additive noise on the wavefield u has the form σ n n ε (x) = σ n (n ε v , n ε p )(x) = σ n (n v ( x ε ), n p ( x ε )) and is independent of the random medium. For simplicity, we assume that the noise entries are real and have the same correlation structure, E{n ε i (x)n ε j (y)} = Φ( x−y ε ), i, j = 1, · · · , 4, where Φ(x) ≡ Φ(|x|) and is smooth. Using (3.6) and (3.12), our model for the measurements is therefore The Wigner transform. Let a ε be the projection of the full Wigner transform W ε onto the + mode, i.e. a ε = Tr(W ε ) T B + , so that, for the a S defined in (3.2), we have a S = F x ε * k a ε . We already know from (3.10) that E{a S } F x ε * kā , withā the solution to the radiative transfer equation (3.9). We subsequently write a ε =ā + δa ε , where δa ε accounts for the random fluctuations. The simplest model for δa ε is obtained for the single scattering approximation which consists in retaining in W ε only terms at most linear in V , so that the related variance will be at most linear inR. This leads to defining the random term W ε with vanishing initial conditions and where As was the case for the pressure, such a definition of W ε 1 does not take into account the absorption factor e −c 0 Σt . We thus formally correct this and set The fluctuations of the amplitude δa ε and the wavefield can be related by writing the field u as E{u} + δu, where δu is obtained in the single scattering approximation. The perturbation δa ε is then the projection on the + mode of the sum of Wigner transforms W [E{u}, δu] + W [δu, E{u}]. Taking into account the external noise σ n n ε (x) and denoting by a ε n the projection of W [E{u}, n ε ] + W [n ε , E{u}] on the + mode as in (3.15), our complete model for the measurements in the single scattering approximation is therefore a ε (t, x, k) =ā(t, x, k) + δa ε (t, x, k) + σ n a ε n (t, x, k).

Analysis of the CB functional
We compute in this section the mean and the variance of the CB functional.

Average of the functional
Assume for the moment that our measurements have the form, for g a smooth function: Above, F x ε (k) is the filter defined in (3.3) and T t is the free transport semigroup introduced in (3.11). Plugging this expression into the CB functional yields: Let us verify first that when D = R 3 , we recover that I C [g](x) = (2π) 3 R 3 dqg 0 (x, q), so that if g 0 is the scalar Wigner transform of a function ψ, we find I C [g](x) = (2π) 3 |ψ(x)| 2 and the reconstruction is then perfect. This is immediate since from the definition of F x ε (k) given in (3.3) we conclude that F x ε (k) = (2π) 3 δ 0 (k) when C = R 3 . When C is finite as is the case in practical situations, one does not recover R 3 dqg(x, q) but an approximate version of it limited by the resolution of the functional. This point is addressed below.
Amplitude and resolution. For simplicity, we suppose that the domain C is chosen such that F x ε defined in (3.3) is independent of x. We suppose in addition that C is a ball centered at zero of a certain radius. We parametrize C by γ ∈ [0, 1] such that if ± εy 2 ∈ C, then we have |εy| ≤ r 0 ε 1−γ for some r 0 > 0. This means that the diameter of the ball is equal to where the parameter N C was introduced in section 2.2. Since the (rescaled) detector has a side l L < 1, we have necessarily by construction that r 0 < l L . We could easily accommodate for anisotropic domains C with additional technicalities. In order to deal with a regular filter, we smooth out the characteristic function of the unit disk and replace it by some approximate function χ(x) ≡ χ(|x|). Rescaling y as y → yr 0 ε −γ , it comes for the filter F ε : Recall that our measurements read a S = F x ε * a ε , where a ε is given by (3.16). The total functional is denoted by I C [a ε ](x) and its average is given by I C [ā](x). The mean of the measurementsā is the sum of a ballistic term a B (t, x, k) = T t a ε 0 (x, k)e −Σ(|k|)t and a multiple scattering term defined in (3.11). Measurements in a homogeneous medium have the form T t a ε 0 (x, k). The CB functional is tailored to backpropagate the data along the characteristics of the transport equation, so as to undo the effects of free transport. Since the multiple scattering term is smoother than the ballistic term, one can expect the inversion operation to produce a signal with a lower amplitude for the multiple scattering part than for the ballistic part. We therefore neglect the multiple scattering in the computation of the average functional. Moreover, the resolution of I C is mostly limited by the filter F and not by the size of the detector since C is included in D. As a consequence, the limitations due to D on the resolution are negligible and we replace D x by R 3 in the definition of I C . We then find, so that the reconstructed image is essentially obtained by a convolution-type relation in the k variable between a localizing kernel and the source. This is clear when ε γ εµ where We can then define the resolution of the functional to be the scale of the localization, which is ε γ /r 0 = N −1 C , or L/N C in dimension variables. Hence, the larger the domain on which the correlations are calculated, the better is the resolution. When the domain is large, more phase information about the wavefield is retrieved, and this is the most useful information to achieve a good resolution. On the contrary, when γ = 0, the phase is lost and imaging is performed mostly using the singularity of the Green's function, which lowers the resolution.
Regarding the amplitude of the signal, we set x = 0 and compute the value of the peak. If the bandwidth parameter µ = ε −α with α > 1 − γ, we find In this case, the convolution is done at a smaller scale than the spatial support of w 0 , so that there is no geometrical loss of amplitude with respect to w 0 (x = 0). When α < 1 − γ, convolution is done at a larger scale and we have Above, θ 1 and θ 2 are such thatk q + εµ(e 1 θ 1 + e 2 θ 2 ) with e 1 ·q = e 2 ·q = e 1 · e 2 = 0 and |e 1 | = |e 2 | = 1. The case α = 1 − γ follows similarly. As a conclusion of this section, we therefore have

Variance of the functional
We compute in this section the variance of the CB functional with our measurements (3.16). For simplicity, we replace c 0 and T by their actual values, c 0 = T = 1. The total variance is defined by For some function ϕ to be defined later on, let us introduce the following quantity The function J is the single scattering approximation of the scintillation function of the Wigner transform and can be seen as a measure of the statistical instabilities. See [4,6,7] for an extensive use of scintillation functions in the analysis of the statistical stability of wave propagation in random media. We have the following technical lemma, proved in section 6.1, that will be used to obtain an explicit expression for the variance V C : Lemma 4.1 w ε is given by where the function H ε reads in the definition of w ε , then The analysis of the variance is somewhat technical in that the scintillation is integrated against singular test functions (due to the application of the CB functional) while a classical stability analysis amounts to integrating against regular test functions. As before, we assume that the detector D is large enough so that its effects on the variance can be neglected and we replace it by R 3 as a first approximation. We have the following two propositions proved in sections 6.2 and 6.3: The variance of the CB functional for the single scattering contribution satisfies and V C is mostly supported on |z| ≤ µε.
The above result gives the optimal dependency on the parameters η, ε, δ, γ and µ. Regarding the contribution of the noise, we have: The variance of the CB functional for the noise contribution satisfies and V C n is mostly supported on |z| ≤ µε.
As before, the result is optimal and the variance decreases as γ goes to zero and resolution is lost. Notice that both variances are essentially supported on the support of the initial source. We now turn to the WB functional.

Analysis of the WB functional
In this section we compute the mean and the variance of the WB functional.

Average of the functional
Using our model (3.13) and the definition of the functional (3.1), we have The cross-range resolution of I W is the same as the classical Kirchhoff functional and given by the Rayleigh formula λL/l.

Variance of the functional
The variance is defined by As before, we distinguish the contribution of the single scattering, denoted by V W , from the one of the noise, denoted by V W n . We then have the following propositions, proved in sections 6.4 and 6.5: The variance of the WB functional for the single scattering contribution satisfies: and V W is mostly supported on |z| ≤ µε.
As for the CB functional, the results are optimal. Regarding the contribution of the noise, we have: The variance of the WB functional for the noise contribution satisfies for all z: In the latter case, the variance is uniform in z and not just supported on the support of the source. This is due to the fact that, contrary to the CB functional which considers the correlation of the noise with the coherent wavefield (i.e. V C n involves (p B ) 2 (σ n n ε p ) 2 ), the variance of the noise for the WB functional does not involve the average wavefield and the information about the source is lost (i.e. V W n involves only (σ n n ε p ) 2 ). Remark also that V W n and V C n are comparable when γ = 1, namely when both functionals have the same resolution. V C n decreases when the resolution worsens. The proof of theorem 2.1 is then a straightforward consequence of (4.1), propositions 4.2, 4.3, 5.1, and 5.2 after recasting ε, η and r 0 ε −γ in terms of λ, L, c and N C .

Proof of lemma 4.1
We start from (3.14) using the notations of section 3.3. We project the initial condition W ε 0 onto the propagating modes and introduce the related amplitudes a ± (recall there are no vortical modes since the initial velocity is zero): where w µ 0 (x, |k| − |k 0 |) := µw 0 (x/µ, µ(|k| − |k 0 |)). Split then S ε in (3.7) into S ε 1 + S ε 2 with obvious notation. We project (3.14) on the + mode and need to compute Tr(S ε 1 W ε 0 ) T (k)B + (k). Direct calculations yield In the same way, we find and finally, recasting δa ε by a 1 for simplicity of notation: Above, stands for real part. The integral equation for a 1 reads: Let us introduce the product function We then have Using the fact that E{V (ξ)V (ν)} = (2π) 3R (ξ)δ 0 (ξ + ν), we find Therefore Accounting finally for the absorption e −Σ , this allows to obtain the following expression for w ε : We have finally which replaced in H ε yields the expression given in the lemma. This is the final result and ends the proof.

Proof of proposition 5.1
We use here the notation of section 3. With δp ε the random fluctuations defined in (3.12), the variance of the WB functional admits the expression Similarly to the CB functional, we neglect the finite size of the detector as a first approximation, so that V W reads dkdk e i(k−k )·z cos |k| cos |k |E{(Fδp ε )(t = 1, k)(Fδp ε )(t = 1, k )}.
Terms in the expression above involving an oscillating exponential are negligible, which shows that V 1 (0) ∼ σ 2 0 e −Σ , and for the same reason as V 2 , V 1 is mostly supported on |z| ≤ εµ. The leading term is therefore V 1 , which ends the proof.

Proof of proposition 5.2
The proof is straightforward, the variance of the WB functional for the noise contribution admits the expression dkdk e i(k−k )·z cos |k| cos |k |E{(F1 D n ε p )(t = 1, k)(F1 D n ε p )(t = 1, k )}.

Conclusion
This work is concerned with the comparison in terms of resolution and stability of prototype wave-based and correlation-based imaging functionals. In the framework of 3D acoustic waves propagating in a random medium with possibly long-range correlations, we obtained optimal estimates of the variance and the SNR in terms of the main physical parameters of the problem. In the radiative transfer regime, we showed that for an identical cross-range resolution, the CB and WB functionals have a comparable SNR. The CB functional is shown to offer a better SNR provided the resolution is lowered, which is achieved by calculating correlations over a small domain compared to the detector. This is the classical stability/resolution trade-off. We obtained morever that the minimal central wavelength λ m that the functionals could accurately reconstruct were identical in the regime of weak fluctuations in the random medium, and that in the case of larger fluctuations, the CB functional offered a better (smaller) λ m (resolution) than the WB functional.
We also investigated the effects of long-range correlations in the complex medium. We showed that coherent imaging became difficult to implement because the mean free path was very small and therefore the measured signal too weak compared to additive, external, noise. In such a case, transport-based imaging with lower resolution is a good alternative to wave-based (coherent) imaging. This will be investigated in more detail in future works.