High frequency analysis of imaging with noise blending

We consider sensor array imaging for simultaneous noise blended sources. We study a migration imaging functional and we analyze its sensitivity to singular perturbations of the speed of propagation of the medium. We consider two kinds of random sources: randomly delayed pulses and stationary random processes, and three possible kinds of perturbations. Using high frequency analysis we prove the statistical stability (with respect to the realization of the noise blending) of the scheme and obtain quantitative results on the image contrast provided by the imaging functional, which strongly depends on the type of perturbations.


Introduction
An amazing fact in the analysis of imaging functionals, that has been recently pointed out and is currently under investigation, is that specific kinds of noise can improve the image quality or drastically reduce computational costs associated to the evaluation of the imaging functional. A strong motivation is the observation that time reversal refocusing is enhanced when the medium is randomly scattering. A time reversal experiment is based on the use of a special device called a time reversal mirror (TRM), which is an array of transducers that can be used as receivers and transmitters. A typical time reversal experiment consists in two steps. In a first step, a point source emits a short pulse that propagates through the medium and is recorded by the TRM used as an array of receivers. In a second step, the recorded signals are time reversed and reemitted by the TRM used as an array of transmitters. The waves then refocus on the original point source location. The striking observation is that refocusing is enhanced when the medium is randomly heterogeneous and scattering compared to the situation in which the medium is homogeneous. Moreover, the refocused pulse is statistically stable in the sense that it does not depend on the realization of the random medium, but only on the statistical distributions of the fluctuations of the random medium [FGPS].
In the context of sensor array imaging, a similar technique is employed. The typical experiment still consists in two steps. The first step is the experimental data acquisition: a point source emits a wave into the medium, the wave is reflected by the singularities in the propagation speed of the medium and is recorded by an array of receivers. The second step is a numerical processing of the recorded data: the recorded signals are time reversed and resent in a numerical simulation into a model medium to locate the singularities in the propagation speed. However, in contrast with physical time reversal, the fictitious medium employed in the imaging process cannot in general capture complex heterogeneities of the original medium. Therefore, research has focused on other approaches to improve the imaging process, such as the use of random sources [DFGS12], [GP10], [HSH08], [SWH11], [VB11], [WNT12].
The classical approach to the imaging problem consists in performing a large number of experiments sounding each time a different source. For each experiment, the signal recorded at each receiver is stored and time reversed. This produced the data matrix. For each experiment, the time reversed data are reemitted (numerically) into the model medium in a new simulation, and the images obtained by each simulation are stacked. While this method provides a very good "image"of the medium, it involves gathering, storing and processing huge amounts of data [Be09], [MDB11].
In the noise blending approach noisy sources are used and they are all sounded simultaneously in one experiment. In this case, the time reversed recorded signals from the physical experiment are stored in a data vector and resent simultaneously into the model medium in a single simulation. This approach allows for considerable savings in both the data gathering, storing and processing stages. However, special care must be put in the choice of the noisy sources in order to ensure that the cross talk terms are very small and do not compromise image quality. This technique can be successfully applied also to the time reversal approach and its analysis bears similarities with techniques for passive imaging, which exploits ambient noise sources to recover travel times from correlations in between recordings at different stations [GP09] [BSC07], to petroleum prospecting [CGH06] and medicine [FCD00].
As remarked, the crucial step for the noise blending approach lies in the choice of the noisy sources. In [DFGS12] it was suggested the use of stationary random Gaussian sources or of randomly delayed source pulses: for these choices of random sources, fourth moment computations show that the algorithm is statistically stable. In the present work we develop some further investigations on this setting and show how to obtain quantitative results on image quality and statistical stability of this algorithm in the high frequency regime, when the goal is to image singular perturbations in the speed of propagation. The result strongly depends on the type of perturbation. We therefore consider three types of perturbations, supported respectively on small balls, thin tubes and thin discs. With a slight abuse of notation we will refer to them as point, line and surface singularities, as they can be thought of as approximations of singular perturbations of the velocity of propagation supported on subspaces of lower dimension.
For each kind of perturbation we first analyze the average image contrast seen between the center of the perturbation and a point far from it: this will provide an hint on the level of difficulty to correctly image these perturbations. An even more interesting result follows: it concerns the quantitative analysis of the statistical stability of this functional, providing the typical contrast seen for the three perturbations. The question of stability of the imaging functional has already been addressed in [DFGS12]: no quantitative analysis was carried out there, but it was shown that a condition for the statistical stability is that the recording time interval T must be large. With a careful analysis of fluctuations produced by the random sources we show that, in the high frequency regime, the typical contrast is actually much better than just of order √ T ; but again the exact order of amplitude of the contrast depends on the shape of the perturbation. Point singularities are easy to observe, while surface type perturbations are the hardest to locate.
The paper is organized as follows. A short presentation of the model and the imaging functional used will complete this introductory section. In section 2 we analyze the average (with respect to the realization of the random time delays used in the blending process) sensitivity of the method to the three types of singularities. In section 3 we study fluctuations due to the stochastic nature of the result, and from the analysis of the typical behavior we obtain conditions ensuring the possibility to accurately image the perturbations. Finally, all the results obtained are collected and discussed in section 4.
Notation. We use boldfaced characters to denote vectors: for example x ∈ R, but x = (x, y, z) ∈ R 3 . B r will denote the ball of radius r in R d for d equals 2 or 3: B r = {x ∈ R d | |x| ≤ r}. The notation A = O(ε) means that the quantity A is exactly of order ε; in order to say that it is of order ε or smaller we will write A ≤ O(ε).

The wave equation
We consider the solution u of the wave equation in a three-dimensional inhomogeneous medium where c(x) is the velocity of propagation of waves in the medium and x = (x, y, z) ∈ R 3 . We rewrite the velocity in the form , where c 0 (x) is the known smooth background velocity (for simplicity we assume it to be constant) and δc −2 (x) is the velocity perturbation that we want to estimate, whose spatial support is contained in some domain Ω ⊂ R 3 . We will detail these perturbations in the following. To simplify the geometry of the model, we shall take Ω = B R to be a ball with a large radius R. We also assume to have N s point sources located at points (y s ) s=1,...,Ns laying on the surface ∂B R . They can either emit (almost) simultaneously the same short pulse waveform, but randomly delayed in time, or independent stationary random signals. Following [DFGS12], we will refer to the first case as (noise) blended sources and to the second as stationary random sources.
For the first case, the source term n(t, x) is of the form The (short) pulse function (f (t)) t∈R is deterministic. Its carrier frequency is ω 0 and its bandwidth is b. The time delays (τ s ) s=1,...,Ns are zero-mean independent and identically distributed random variables with probability density function p τ (t).
For the second case, the source term is given by where the random functions (n s (t)) t∈R , s = 1, . . . , N s , are independent, zero-mean, stationary Gaussian processes with autocorrelation function The direct and inverse problems can be formulated in terms of the background Green's function, the fundamental solution of the wave equation (1.1). For an homogeneous medium with constant background velocity c 0 , in the Fourier domain the Green's function is given bŷ Here, the Fourier transform of a function f (t) is defined bŷ

Direct and inverse problems
We introduce here the scattering operator, that is, the mapping from velocity perturbations to the data, in the Born approximation [BCS]. We assume that sources (located at points y s , s = 1, ..., N s ) are disposed on the surface of the ball B R containing the perturbations, and are dense enough (ideally, closer than half of the central wavelength) so that a continuum approximation can be used. Signals are observed at the passive sensor array (x r ) r=1,...,Nr for some large time interval [−T /2, T /2]. For noise blended sources, the recording time T should be much larger than the typical travel time, so as to guarantee that the backscattered signals are completely recorded. For stationary random sources it must be taken much larger that the inverse of the bandwidth of the noise sources (i.e. much larger than the decoherence time).
The recorded data consists of the vector d(t) = (d(t, x r )) r=1,...,Nr of the signals recorded by x r for t ∈ [−T /2, T /2]. These data are modeled by the scattering operator F : In the Fourier domainQ This is the formulation of the direct problem: the expression of the data set in terms of the velocity perturbation.
The imaging problem (inverse problem) aims at inverting the map F in order to reconstruct the velocity perturbation δc −2 (x) x∈BR from the data set d(t) t∈[−T /2,T /2] . The usual (least-square) approach would consist in applying the operator (F * F ) −1 F * to the data set d, where the adjoint of the scattering operator F is given by (1.4) However, the full least-square inversion is in practice too complicated and the normal operator F * F is usually dropped in the inversion process. In [DFGS12] it was shown that for T large the normal operator is statistically stable (i.e. its fluctuations are smaller than its expectation) and that its statistical average is close to the identity operator (more precisely, the kernel F * F (x, x ′ ) concentrates near the diagonal x = x ′ ), proving that this procedure can indeed provide a reasonable estimate of the velocity perturbation.
In the following sections we perform a detailed analysis of this imaging functional in the high frequency regime, obtaining quantitative results on the statistical stability of the method for different types of perturbations.
Since the main advantage of our approach lays in the drastic reduction in storage and computational costs, let us stress that the computation of the adjoint operator F * can be done quite easily.
Remark 1. The adjoint operator F * can be computed solving only two wave equations as follows. First, compute the wave u(t, x) emitted by the original source, which is to say solve the wave equation with source n(t, x) and background velocity c 0 : Second, compute the anti causal solution v(t, x) of the wave equation with source term Correlating the two wave solutions at a point x in the search window produces the imaging functional and using the definition (1.2) we obtain

Analysis of the imaging functional
Let us start by the analysis of a kernel which will appear in the imaging functional we have to study. In our model, sources and receivers are located on the surface of the ball B R containing the perturbations and are dense enough so that a continuum approximation can be used. We can therefore rewrite the kernel where ρ(x r ) is the surface density of receivers or sources. We also assume that perturbations are located near the center of the ball; this means that the medium is homogeneous outside of a ball B r of radius r ≪ R. Then, we have the approximate identity [GP09,Proposition 4 which follows from Green's identity and the Sommerfeld radiation condition. This result can be viewed as a version of the Helmholtz-Kirchhoff integral theorem. Using the function sinc(x) = sin(x)/x, the right hand side of the above equation can be rewritten as and assuming that receivers and sources have constant density on the surface ∂B R we obtain We return to the analysis of the imaging functional, which estimates the velocity perturbations, given by (1.4). Using equation (1.3) and the definition of the vector d(t) of recorded signals, we can write the imaging functional as where δc −2 (x) is the estimated velocity perturbation. We will use · to denote the statistical average with respect to the distribution of the random time delays or the stationary random source signals.
The statistical average of the kernel K (which is not the same kernel as the K ω used earlier) is given by for noise blended sources, and by for the stationary random sources, see [DFGS12]. We stress that the focus of this paper is on the high frequency analysis of this functional. We will analyze its performances in localizing singularities in the background velocity of propagation, in the high frequency regime η ≪ 1, where Note that the wavelength is 2πη.

Expected contrast of the estimated perturbation
Even if the estimated perturbation δc −2 (x) provided by the imaging functional does not have the exact shape of the real perturbation, due to the different approximations used, it still shows a peak on the actual location of the velocity perturbation. We study in this section the expected (average in the statistical sense) contrast seen in the estimated velocity perturbation between the location of the real perturbation and points far from it. The expected estimated perturbation for noise blended sources is given by When the bandwidth of the source pulse is smaller than the central frequency, it is enough to study the behavior of the spatial integral I 1 (x) at the central frequency. For simplicity we drop the dependence on ω. It turns out that T · I 1 (x) is the quantity one has to study also in the case of stationary random sources. All the computations presented in the remaining part of this section are performed in the noise blended case; to obtain the different orders of amplitude provided by the imaging functional for stationary random sources it suffices to multiply the results by T .
To simplify the presentation of some computations, we assume that the support of the perturbation is centered on the center of the ball B R on the surface of which are located the sources and receivers, and choose this center as the origin of our system of coordinates.
For the purposes of this section it would not be necessary to distinguish between the scale of the wavelength of the signals (η) and that of the size of the perturbation (ε). However, this distinction turns out to be crucial in the analysis of fluctuations carried out in section 3.

Point singularities
Let us start by considering a point singularity, that is to say, a singularity whose support has a very small diameter. We model it with a perturbation of the velocity of propagation that is supported on a ball of radius ε: (2.1) To simplify computations, we take α to be constant. Then, one observes that it only enters formulas as a multiplicative constant (squared in section 3): since it is of no relevance to our scopes, we set it equal to 1, also for the other types of perturbations.
Changing variables, we have: At the center of the perturbation, x = 0, we have For ε ≃ η we have that I 1 (0) = O(εη 2 ). But if ε ≪ η, the order becomes O(ε 3 ). For |x| = O(1), what gives the order of amplitude of I 1 (x) is the decay of the sinc(x) function, which goes approximately as 1/x. We have that Remark that the bound found is sharp, since there are x for which I 1 is exactly of order ε 3 η 2 . We see that the difference in amplitude between the centre of the perturbation and a point far from it is significant: it is of order η −2 when ε ≪ η, ε −2 if ε ≃ η.

Line singularities
We consider now line-type singularities, which is to say an almost one-dimensional perturbation of the velocity of propagation. We model it by a perturbation supported on a cylinder of radius ε: (2.2) We have At x = 0 this term reduces to which is rapidly oscillating in z ′ . In the integral in z ′ we can therefore approximate sin 2 by its mean: Far from the perturbation the integral is of order ε 2 η 2 . For example, for x 2 + y 2 = C = O(1) we have Therefore, the difference in amplitude seen between the centre of the line perturbation and a point far from it is of order ε −1 .

Plane singularities
Let us consider now singularities that are approximately two-dimensional: we call them surfacetype singularities and model them by a perturbation of the velocity which is supported on a disc of thickness ε: (2.3) With the notation introduced above we have because the sin 2 in the first line is rapidly oscillating in r.

For x = O(1) we have instead that
For this type of singularities, the difference in amplitude seen between the centre of the perturbation and a point far from it is very weak: it is only of order |ln(ε)|. This already hints to the fact that, with the imaging functional we consider, surface-type singularities are harder to locate than the other types.

Fluctuations
We have obtained in the previous section the average contrast of the imaged perturbation. We want now to find confidence intervals for the typical contrast observed during the experiment. To do so, we must analyze the fluctuations in the result. They are given by the standard deviation of the estimated perturbation δc −2 , which at a point x is given by (3.1) We need to compute the standard deviation at the location of the perturbation and at points far from it, and compare them with the expected amplitude of the estimated perturbation. Using (1.5) to write explicitly (3.1) we get In [DFGS12] a formula was obtained for the variance of the kernel K. It is possible to carry out the same computations for the covariance: for noise blended sources one obtains ω, x, y s ) 2Ĝ ω, x ′ , y s )Ĝ ω, x ′ , y s ) .
Here, T τ = p 2 τ (t) dt −1 and p τ (t) is the probability density function of the random time delays.
We will see that this quantity should be large so that the kernel F * F (x, x ′ ) is statistically stable. Using the method of Lagrange multipliers one can show for example that the maximal value of T τ amongst all probability density functions p τ compactly supported in [−τ max , τ max ] is obtained for the uniform density over the interval and gives T τ = 2τ max . Since τ max must be at most of the order of the recording time τ max ≃ T /2, to obtain a large value of T τ one should take T large too (recall that in section 1.2 we had already assumed T to be large).
Observe that in the above equation for the covariance, in each of the two terms on the right hand side we are summing N s terms with a minus sign and N 2 s with a plus sign, and they are all of the same order. The contribution of the terms with a minus sign is therefore small, and we can use the results of section 1.3 to simplify the above equation into Similar computations for the case of stationary random sources, where the terms with a minus sign do not even appear, leads to Let us continue the computations in the case of noise blended sources. Putting everything together, we find that the standard deviation of the imaging functional at a point x is given by Again, we focus on the spatial integral I 2 (and to simplify notations we drop the dependence on ω): the quantity we need to study is I 2 (x)/T τ 1/2 for the noise blended sources and T · I 2 (x) 1/2 for the stationary random sources. The computations presented below are performed in the noise blended sources setting. To obtain the corresponding standard deviation for stationary random sources it suffices to substitute the factor 1/T τ (or 1/ √ T τ ) by T (or √ T ). For comparison with the average amplitude, recall that for stationary random sources the results obtained in the previous section have to be multiplied by a factor T .
We can rewrite the integral I 2 (x) as the sum of the two integrals which is to say the double integral of the two kernels applied to the perturbation δc −2 (x ′ )δc −2 (x ′′ ) , for every one of the three types of perturbation studied above. However, since J 2 (x) = I 2 1 (x) , we will only need to analyze J 1 .
For the following estimates it is important to keep separated the scale of the dimension of the perturbation (ε) from the scale of the wavelength of the sources (η). Their relative amplitude will be specified, but to help the reader to keep track of the different orders, we stress that we will always have ε ≤ η.

Point singularities
We return to the case of point perturbations introduced in the previous section and modeled by (2.1). In this case, even a very rude estimation is sufficient to obtain a bound which guarantees that this perturbation can be imaged. Since |sinc| ≤ 1, changing variables we get where the constant K = 4/3π comes from Cauchy-Schwarz inequality. Since the sinc function is bounded, we have obtained a bound of order ε 6 . This is only a rough upper bound, but it is not necessary to look for an improvement since it is already of the same order of the integral of the second kernel, for which we have Using the decay of the sinc function, we can find near the perturbation (|x| ≃ ε) a bound for I 2 of the same order. Oscillations are therefore of order ε 3 / √ T τ . Recall that the average value observed on the peak is of order ε 3 for ε ≪ η, so that the typical value observed remains of the same order due to the large factor T τ . The same results holds true also for ε ∼ η, namely the typical value observed is of the same order of the average value.
Far from the perturbation the integrals of the two kernels decrease. Using the bound (3.2), the integral of the first kernel can be bounded for |x| = O(1) by With some more work, one can show that this bound is sharp, at least for ε ≪ η. This can be done using the Fourier representation of the sinc function written in spherical coordinates u = u(r, θ, φ) ∈ R 3 : where S 2 = {x ∈ R 3 : |x| = 1} is the unitary sphere in R 3 . We can write where we have used the assumption ε ≪ η. Simplifying this equation and using again (3.3), we get for |x| = O(1) As for the integral of the second kernel, the bound we have is of order ε 6 η 4 . For ε ≪ η we have therefore In the general case ε ≤ η, the above equation becomes an upper bound.
Recall that the statistical average of the imaging operator far from the perturbation is of order ε 3 η 2 . Assuming that T is large, but still 1 ≪ T τ ≤ 1/η 2 , the above result implies that the typical value observed is at most of order ε 3 η/ √ T τ (it is exactly of this order for ε ≪ η). Therefore, the typical contrast is still at least of order √ T τ /η, allowing for a precise location of the perturbation (both when ε ≪ η and ε ∼ η).

Line singularities
Consider the case of line singularities, modeled by (2.2). Using rude estimations similar to the ones presented above, we could only bound the integral of the (absolute value of the) first kernel near the origin with something of order ε 4 η 2 ln 2 (ε). This means that, in order to be sure to able to image the perturbation, we would need to have ε| ln(ε)| ≪ η. But we can do better.
Assuming simply ε ≪ η, we can approximate It is then possible to use the Fourier representation of the sinc function to obtain the amplitude of oscillations. At x = 0 we can rewrite the integral of the first kernel as an integral over C 1 = B 1 × [−1, 1] ⊂ R 2 × R, use (3.4) and integrate in x ′ , y ′ , x ′′ , y ′′ : and using the Fourier representation introduced above we get Since the function s → s 0 sinc(u) du is uniformly bounded in s, we get that For the second kernel, we have J 2 (0) = I 2 1 (0) = O(ε 2 η 4 ) .
Remark that for |x| ≃ ε we can still bound fluctuations in the same way, because (3.4) still holds, and the integral We see that fluctuations near x = 0 are of order εη 2 / √ T τ . Since the average value observed at x = 0 for the imaging functional is of order εη 2 , it is thanks to the large factor T τ ≫ 1 that we get the statistical stability of the operator. This means that the typical value observed on the perturbation remains of order εη 2 . When x is far from the perturbation, oscillations are even smaller. Indeed, we can proceed as in (3.2) to find a bound for the integral of the absolute value of the first kernel. We get Since the integrand is bounded, the bound we get is of order ε 4 η 2 . The second kernel is of a higher order, J 2 (x) O(ε 4 η 4 ). This bound implies that the typical value observed far from the perturbation is of order at most ε 2 η/ √ T τ , so that the contrast is at least of order η √ T τ /ε.

Plane singularities
We turn now to analyze fluctuations in the case of surface-type perturbations, modeled by (2.3). The difficult part is again to obtain good estimates on the integral of the first kernel, for which we use the Fourier representation of the sinc function obtained in (3.3). At x = 0 we have where the approximate equality holds for ε ≪ η and ⊥ denotes the projection on the last two coordinates: We can rewrite also the integrals on S 2 as (twice the) integrals on the projection B 1 ∈ R 2 , compute the integral in dx ′ ⊥ and change variables: Here, J 1 is the Bessel function of the first kind. Let us focus on the integral inside the square brackets. Observe that the origin of our system of coordinates is always inside B 1/η (w ⊥ ). Changing to polar coordinates we have where we denote by φ w the square root term (written in polar coordinates), and the function ρ w (θ) takes values in [1/η − |w ⊥ |, 1/η + |w ⊥ |] ⊂ [0, 2/η]. We claim that the integral term Y is bounded. This can be proved integrating by parts in r. Remark that the square root term is concave (as a function of u ⊥ ), take its maximum over the domain of integration B 1/η (w ⊥ ) at the center of the ball and is zero at the boundary. Therefore, for every fixed (θ, w), the function φ w (θ, ·) is still concave and bounded by 1. Then, one easily obtains that the integral of the absolute value of its derivative in r is bounded by 2. Also, the antiderivative of J 1 (r) is the Bessel function of order zero −J 0 (r), which is bounded (the maximum of its absolute value is taken at r = 0, and J 0 (0) = 1). Putting everything together, we get This proves the claim.
Using again the boundedness of the square root, we can bound the integral in dw ⊥ by π/η 2 . We have therefore obtained a bound for J 1 (0) of order O(ε 2 η 4 ). This is only an upper bound, but there is no need to look for an improvement, since it is already of a smaller order than the integral of the second kernel, for which we have J 2 (0) = I 2 1 (0) = O(ε 2 η 4 ln 2 (ε)) .
Far from the perturbation, oscillations are even smaller. Denote x = (x, x ⊥ ) ∈ R × R 2 and u = (u 1 , u ⊥ ) ∈ R × R 2 ; the same notation is used for v and w. Let us look at the integral of the first kernel; following the computations presented above we have For |x| ≫ 1, the last two integrals above are now much smaller than the corresponding ones for x = 0. This is due to the fact that for |x| large, at least one of the two exponential terms, which have mean zero, is rapidly oscillating with respect to J 1 . We therefore have that J 1 (x) is at most of order ε 2 η 4 . Far from the perturbation, the (sharp) bound we have on the integral of the second kernel is of the same order: J 2 (x) O(ε 2 η 4 ). We have obtained that Thanks to the large factor T τ , fluctuations are therefore smaller than the average value given by the imaging functional, both on the perturbation and far from it. The typical contrast remains therefore of the same order as the average contrast, namely of order | ln(ε)|.

Conclusions and comments
We have analyzed the imaging functional given by (1.4)in the high frequency regime (η ≪ 1) with respect to small perturbations (ε ≪ 1) of the background velocity of propagation. Using a suitable disposition of the sources and receivers, we have been able to obtain quantitative estimates on the (average, with respect to the realization of the random time delays or the stationary random source signals) sensitivity of the imaging functional. The image presents a peak on the location of the perturbation, and the contrast is of order η −2 for point perturbations, of order ε −1 for line perturbations, and only of order | ln(ε)| for surface perturbations. The most interesting result obtained in this paper concerns the quantitative analysis of the statistical stability of this functional, providing the typical contrast seen for the three perturbations considered. The question of stability of the imaging functional has been addressed in [DFGS12]: no quantitative analysis was carried out there, but it was shown that a condition for the statistical stability is that the quantities T (for stationary random sources) and T τ (for noise blended sources) must be large. For random time delays uniformly distributed on the interval [−τ max , τ max ], T τ large means that τ max must be large, which in turn implies that the recording time T ≃ τ max must be large.
An important fact is that the typical contrast found only depends on the type of perturbation one is trying to image, and not on the method used. All results are described below for noise blended sources, but the corresponding contrast for stationary random sources are obtained simply substituting T τ with T .
We have shown that for point perturbations, both when ε ≪ η and ε ∼ η, fluctuations due to the stochastic nature of the method are small, and the typical contrast is at least of order √ T τ /η ( √ T /η for stationary random sources): point perturbations are easy to find. For line perturbation the situation is different. We can image with a satisfactory accuracy and reasonable recording time T ≫ 1 only very thin line perturbations, ε ≪ η. The typical contrast in this case is at least of order √ T τ η/ε. For plane perturbations the average contrast is quite poor, only of order | ln(ε)|. However, for very thin perturbations , ε ≪ η, also the typical contrast is of the same order.
These results are summarized in the following tables, where we list for the three type of perturbations considered the order of the average value given by the imaging functional and of the standard deviation at the center of the perturbation and far from it. Table 1: Noise blended sources: mean δc −2 and standard deviation (S) of the estimated velocity perturbation at the center of the perturbation (x = 0) and far from it (|x| ≫ 1), in the regime ε ≪ η ≪ 1. The cases of point, line and disc singularities are displayed.