J-D Benamou
email: jean-david.benamou@inria.fr

F Collino

S Marmorat

Numerical MicroLocal Analysis of 2-D Noisy Harmonic Plane Wave and Source Point Wavefields

Keywords: plane waves, point source, inverse scattering AMS classification scheme numbers: 78A05, 78A46, 78M35 P p=1

We present a mathematical and numerical analysis of the stability and accuracy of the NMLA (Numerical MicroLocal Analysis) method [3] and its discretization. We restrict to homogeneous space and focus on the two simplest field data cases : 1) Noisy plane wave packets, 2) Noisy point sources solution. A new stability result is obtained through the introduction of an impedant observable. The point source analysis leads to a modified second order (curvature Dependant) correction of the algorithm. As NMLA is local, it can be applied the second order improved version can be applied to general data (heterogeneous media). See [2] for a pedestrian description of the algorithm and an application to source discovery.

Introduction

The use of a formal high frequency asymptotic solution ‡ of the Helmholtz equation can be traced back to a 1911 paper of Debye on "the article of Sommerfeld", see [1]. It marked the beginning of the development of modern Geometric Optics numerical models and methods for the simulation of high frequency wave propagation which have been used in the industry for decades [5] [8]. In [3], we considered the inverse problem of numerically and locally extract the Geometric Optics components from harmonic wavefields. We also introduced a new approach named NMLA for Numerical MicroLocal Analysis to solve it. It is based a family of linear operators (or filters) B k , indexed by the wavenumber k (see (4) for the new version of the filter). The filters are applied to observable harmonic wave data collected on a circle of fixed radius r 0 around an observation point x 0 , in [3] we used the Dirichlet data : U k (ŝ) = u k (x), x = x 0 + r 0 ŝ, ŝ = (cos θ, sin θ).

(

) 1
When B k is applied to an observable constructed from a finite number of plane waves :

u k (x) =
the output is an angular function presenting P picks of amplitude b p at points localized approximatively at the ŝp 's. Both the localization of the picks and the quality of the evaluation of the amplitudes b p increases with the frequency as we move in the high Frequency asymptotic regime, thus allowing in principle to estimate the number of high frequency wavefronts passing through x 0 , their directions and amplitudes. The stability of the original NMLA method was not clear and a L 2 Tychonov regularization was proposed in [3], see also [12] for a L 1 regularization. This paper introduces a new "impedant" observable and associated filters for which rigorous mathematical stability and localization error estimates are obtained in the case of noisy plane wave data. This is the topic of §2. In particular, Estimate (7) shows that the new "impedant filter" is very robust to noise.

The method however still relies on the plane wave assumption. General wave data will therefore carry linearization errors growing with the frequency which cannot be avoided. It can be rigorously analyzed on point source solution ( §3.1). Ideally, second order terms should be identified and would allow to get a good approximation of the high frequency geometric optics amplitudes. It seems to be a quite difficult open problem. The closest related work that we know of, is the generalized Radon transform [4] which can be used to identify non linear patterns in images. This is not the philosophy of NMLA which uses the asymptotic geometric optics model (actually only the Eikonal equation so far) to reduce the complexity of the problem. As of today, we are still unable to incorporate the full second order term of the linearization in the identification process and use the transport equation for the geometric amplitudes.

So, we turned to the simplest solution of the Helmholtz equation after the plane wave, the above mentioned point source fundamental solution (in 2-D the Hankel function). Using classic FMM type asymptotic techniques [10], a second order correction can be analytically calculated (see §3) for the NMLA filters applied to a point source solution in homogeneous space. This is the simplest "constant curvature" non linear wave data. As we use this correction locally, it can be applied to general data, just assuming that instead of the the plane wave linearization we approximate the underlying high frequency solution by a constant curvature wavefront with an unknown virtual source. This correction to the classic NMLA filter is second order accurate, see [2] for a pedestrian presentation and numerical results.

§4 discusses the choice for the various discretization parameters of the NMLA filters.

§5 introduces an additional Gaussian filtering procedure which aims at removing interferences between plane waves when more that one. This is important in view of angular decomposition of the observable as the "constant curvature" correction of §3 only works when the wave data radiates from a single source.

A numerical illustration of these results is presented in §6.

Stable extraction of angles from noisy plane wave packets

A new impedant observable for NMLA

As discussed in the introduction, our new NMLA observable uses the impedant quantity on the same circle of radius r 0 centered at x 0 :

U k (ŝ) = 1 ik ∂ r 0 u k (x) + u k (x) , x = x 0 + r 0 ŝ, ŝ = (cos θ, sin θ). (3)
It is known that the knowledge of U k (ŝ) = U k (θ) and that u k satisfies the Helmholtz equation inside the circle, is enough to reconstruct the whole solution since there is no resonance for an impedant Helmholtz problem : our observable contains hence all the information. As mentioned in the introduction, the observable in [3] was the restriction of u k on the circle (compare with (1)). It does not allow to reconstruct the solution inside the circle for some exceptional values of k (Dirichlet resonance) and this deficiency underlies technical and practical important difficulties which are removed by (3).

Let us give the expression of the NMLA filters B k for the new observable :

       B k U(θ) = 1 2L k + 1 L k ℓ=-L k (F U) ℓ e iℓθ (-i) ℓ (J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 )) with L k = max(1, [kr 0 ], kr 0 + (kr 0 ) 1 3 -2.5 ), (4)
where J ℓ (x) is the Bessel functions of order ℓ and argument x, J ′ ℓ (x) their derivative and (F U) ℓ is the ℓ-th Fourier coefficient of U (F U) ℓ = 1 2π 2π 0 U(θ)e -iℓθ dθ.

Note that the denominator in (4) never vanishes since J ℓ (x) and J ′ ℓ (x) have no common zeros (cf. Wronksian formula for Bessel Functions (A.2)).

Preliminary results for perturbed plane waves packets

We now assume noise perturbation of the plane wave data (2) :

U k (ŝ) = U plane k (ŝ) + δU(ŝ). (5)
The first term corresponds to the plane wave packets, the second one to the perturbation δu k . The filters being linear we get

B k U k (θ) = B k U plane k (θ) + B k δU(θ),
and as already noticed in [3] the part of U k corresponding to the plane waves can be computed analytically; we have

B k U plane k (θ) = P p=1 b p ψ L k (θ -θ p ),
where

ψ L (ϕ) = sin((L + 1 2 )ϕ) (2L + 1) sin ϕ 2 .
The proof relies on the Jacobi Anger Expansion, [11, page 66]; we have

1 ik ∂ r 0 + 1 (e -ikr 0 cos(θ-θp) ) = ∞ ℓ=-∞ 1 ik ∂ r 0 + 1 (-i) ℓ J ℓ (kr 0 )e iℓ(θ-θp) , or 1 ik ∂ r 0 + 1 (e -ikr 0 cos(θ-θp) ) = ∞ -∞ (-i) ℓ (-iJ ′ ℓ (kr 0 ) + J ℓ (kr 0 )) e iℓ(θ-θp) ,
and, therefore,

B k 1 ik ∂ r 0 + 1 (e -ikr 0 cos(θ-θp) ) = 1 2L k + 1 L k ℓ=-L k e iℓ(θ-θp) = ψ L k (θ -θ p ).
To understand how applying B k to the observable provides a filtered output with picks localized at the angles of incidence, we use in the next section the properties of the function ψ L (ϕ) and a rather technical result that shows that one controls the norm of the operator B k independently of the wavenumber k.

We start with function Ψ L 's properties :

(i) Ψ L (ϕ) is even; its maximum is reached at ϕ = 0 sup Ψ L (ϕ) = 1 = Ψ L (0).
(ii) Ψ L (ϕ) take values in the range of its maximum only in the vicinity of 0; we have

|Ψ L (ϕ)| > 0.45 ⇒ |ϕ| < 2 L(L + 1
) .

(iii) Ψ L (ϕ) is "small" far away from 0; we have

Ψ L (ϕ) ≤ π 2L + 1 1 |ϕ| .
Concerning the bound on the norm of the B k operators, we have the following result :

There exists a pure constant B ⋆ (i.e. independent of k) such that

sup θ |B k U(θ)| ≤ B ⋆ √ 2π U L 2 (0,2π) ≤ B ⋆ U L ∞ (0,2π)
This result strongly depends on the choice of the Fourier modes truncation level L k in (4). Numerical experiments also shows that we have

B ⋆ ≤ 1,
and, even sharper B ⋆ ≤ 0.89, if we deal with high frequencies such that kr 0 > 3. These bounds on B ⋆ mean that the filter operator does not amplify too much the noise.

The proof starts with Cauchy-Swartz

|B k U(θ)| ≤   1 2L k + 1 L k ℓ=-L k 1 J 2 ℓ (kr 0 ) + J ′ 2 ℓ (kr 0 ) 1 2   ∞ ℓ=-∞ |(F U) ℓ | 2 1 2 , then Parseval |B k U(θ)| ≤   1 2L k + 1 L k ℓ=-L k 1 J 2 ℓ (kr 0 ) + J ′ 2 ℓ (kr 0 ) 1 2   1 2π 2π 0 |U(ϕ)| 2 dϕ 1 2 , (6)
and finally on a technical results (cf. appendix A) that shows that the term between hooks is bounded independently of k.

Confidence interval for angle extraction

Assuming θ ⋆ denotes the angle for which θ → |B k U k (θ)| is maximum, we expect this will be the ray direction of the most energetic plane wave component. This section gives an error estimate on θ ⋆ . We begin with a single perturbed plane wave data :

B k U k (θ) = b 1 ψ L k (θ -θ 1 ) + B k δU. Assume |θ -θ 1 | > 2π 2L k +1
, we use the third property of function Ψ L k and obtain

|B k U k (θ)| ≤ |b 1 | π 2L k + 1 2L k + 1 2π + B ⋆ δU L ∞ ≤ 1 2 |b 1 | + B ⋆ δU L ∞ , whereas, at θ = θ 1 , |B k U k (θ 1 )| ≥ |b 1 | -B ⋆ δU L ∞ .
Thus, if perturbation due to the noise satisfies

2B ⋆ δU L ∞ < 1 2 |b 1 |, (7)
the maximum of |B k U plane (θ)| will be reached in a confidence interval centered around θ 1 with size π

L k + 1 2
. When the frequency increases, since L k ∼ kr 0 , the width of this interval will diminish and the angle of incidence is recovered with more precision. Remark that Inequality (7) is not very restrictive since a numerical estimate gives 1 4B⋆ ≃ 0.28: if the noise level does not overpass 28 % the angle will be detected within an interval of size ∼ 2π k .

The general case of P plane waves is more technical, for simplicity we only detail the case P = 2. We assume

B k U plane (θ) = b 1 Ψ L k (θ -θ 1 ) + b 2 Ψ L k (θ -θ 2 )δU, with |b 2 | < |b 1 |.
We anticipate the maximum to be reached for θ ∼ θ 1 . At θ 1 , we have

|B k U plane (θ 1 )| > |b 1 | - π 2L k + 1 1 |θ 2 -θ 1 | -B ⋆ δU L 2 , while for θ ∈ I = θ, max p=1,2 |θ -θ p | > 4π 2L k + 1 , |B k U plane (θ ) | < π 2L k + 1 2L k + 1 4π (|b 1 | + |b 2 |) + B ⋆ δU L ∞ < 1 2 |b 1 | + B ⋆ δU L ∞ .
The maximum is therefore be outside

I if 2B ⋆ δU L ∞ + π 2L k + 1 |b 2 | |θ 2 -θ 1 | < 1 2 |b 1 |.
A sufficient condition to ensure that the maximum is reached near θ 1 , the angle of maximal amplitude, is

|b 2 | + |b 1 | π 2L k + 1 1 |θ 1 -θ 2 | + 2B ⋆ δU L ∞ ≤ |b 1 |.
This inequality is asymptotically satisfied when the amplitude difference is not larger than the noise level. When |b 1 | is too close to |b 2 |, the maximum may be erroneously located nearer to θ 2 .

Let us point out that the width of the confidence interval I is twice larger than in the one plane wave case. This is a direct consequence of the possible interferences in the image by the NMLA filter of the two plane wave component if the angles are too close. Adding a Gaussian filter as proposed in §5 allows to consider treat well separated waves.

Source points identification using NMLA

Points sources viewed as local plane waves

After noisy plane wave packets (2), the simplest Helmholtz data are point sources solutions in homogeneous media. It also is the only other case we can handle analytically. So, we assume the data is generated from a P collection of source points x p plus the now usual noise δu k :

u k (x) = P p=1 a p i 4 H (1) 0 (k|x -x p |) + δu k (x). (8)
We have now 3 × P unknowns : the a p 's, θ p 's and R p 's, where

x p = x 0 + R p ŝp . = x 0 + R p (cos θ p , sin θ p ).
When the frequency goes to infinity, we can substitute to the Hankel function its asymptotic development, [11, page 65]

u k (x) = P p=1 a p e ik|xp-x| -ik 8π|x p -x| + δu k (x).
If we want to use NMLA, a natural idea is to linearize the data on the observation circle x 0 + r 0 ŝ around the observation point x 0 to find our usual plane wave packet. If r 0 is "small enough" with respect to all the R p 's, we will get have

u k (x 0 + r 0 ŝ) = P p=1 b p e -ikr 0 cos(θ-θp) + δu k (x), with b p = a p e ikRp -ik 8πR p . (9)
This means that, on the disk of center x 0 and radius r 0 , the high frequency field is the sum of P plane waves plus some noise perturbation. The next step is to apply the previous filter (i.e. our operator B k ) to the impedant quantity (3). However, we need to be careful about the errors accumulated over the successive approximations (see the number of tildes rising over the perturbation term !).

If we just focus on the linearization error, we see that we used one term for the modulus and two terms for the phase in the expansion

k|x p -x| = k|x p -x 0 | -k(x -x 0 ) • ŝp + k 2|x p -x 0 | |x -x 0 | 2 -(x -x 0 • ŝp ) 2 + . . . , or, for x = x 0 + r 0 ŝ k|x p -(x 0 + r 0 ŝ)| = kR p -kr 0 cos(θ -θ p ) + kr 2 0 2R p sin 2 (θ -θ p ) + . . . (10)
For a fixed radius r 0 , since the second order term goes to infinity with k, the NMLA filter is likely to fail. One would like to scale down the radius r 0 with the frequency to control it. However, we saw in §2 that stability and accuracy of the angle extraction depend on the truncation level L k which is homogeneous to r 0 . Somehow a finer analysis is needed.

Analysis of a NMLA filtered source point

We first assume that the data field corresponds to a unique source point

u k (x) = a 1 i 4 H (1) 0 (k|x -x 1 |) + δu k (x).
Then, we use the Graf addition Theorem, [11, page 66]: when r 0 is smaller than

R 1 = |x 1 -x 0 |, we have u k (x 0 + r 0 ŝ) = a 1 i 4 +∞ ℓ=-∞ H (1)
ℓ (kR 1 )J ℓ (kr 0 ) e iℓ(θ-θ 1 ) + δu k (x 0 + r 0 ŝ), where θ (resp. θ 1 ) is the angle associated to ŝ (resp. to

x 1 -x 0 |x 1 -x 0 |
). The impedant quantity on the circle with the usual noise perturbation is

U k (θ) = U pt k (θ) + δU k (θ)
, where

U pt k (θ) = ( 1 ik ∂ r 0 + 1)u k (x 0 + r 0 ŝ) = a 1 i 4 +∞ ℓ=-∞ H (1)
ℓ (kR 1 )(J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 )) e iℓ(θ-θ 1 ) , or, after a change of variable

U pt k (θ) = a 1 e ikR 1 √ -ik 8πR 1 +∞ ℓ=-∞ C ℓ (kR 1 ) (-i) ℓ (J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 )) e iℓ(θ-θ 1 ) , with C ℓ (kR 1 ) = i ℓ ikπR 1 2 e -ikR 1 H (1) ℓ (kR 1 ). ( 11
)
It is important to notice that when kR 1 is large enough C ℓ (kR 1 ) is a perturbation of 1 due to the asymptotic, [13, page 198 ],

C ℓ (kR 1 ) ∼ 1 + ∞ m=1 (ℓ, m) (-2ikR 1 ) m , (12) with
(ℓ, m) = (4ℓ 2 -1 2 )(4ℓ 2 -3 2 ) . . . (4ℓ 2 -(2m -1) 2 ) 2 2m m! .
The zeroth order approximation C ℓ (kR 1 ) ∼ 1 corresponds to a plane wave linearization of the data. As in §2.2, the NMLA filter will map it to a the shape function Ψ L k (θθ 1 ).

Consequently, we split the data between linear and remaining quantities and obtain (cf. definition of b 1 in (9) and the expression of the filter B k in (4))

B k U pt k (θ) = b 1 1 2L k + 1 L k -L k C ℓ (kR 1 ) e iℓ(θ-θ 1 ) , (13)
or

B k U pt k (θ) = b 1 Ψ L k (θ -θ 1 ) + δβ(θ),
where

δβ(θ) = b 1 1 2L k + 1 L k -L k (C ℓ (kR 1 ) -1) e iℓ(θ-θ 1 ) . ( 14
)
δβ is precisely the image of the linearization error by the filter. It can be bound as follows

|δβ(θ)| ≤ |b 1 | 1 2L k + 1 L k ℓ=-L k |C ℓ (kR 1 ) -1| = |b 1 | E(L k , kR 1 ).
We can then use the results of §2.3 , applying the filter to the impedant data linked to a unique noisy source point. The only change is the additional linearization noise :

B k U k (θ) = b 1 Ψ L k (θ -θ 1 ) + δβ(θ) + B k δU k (θ).
The supremum of this is function will correspond to an angle located at a distance less than π

L k + 1 2 of θ 1 , as soon as 2B ⋆ |b 1 | δU L ∞ + E(L k , kR 1 ) ≤ 1 2 .
To go further, we need an estimate of the perturbation; it can be done numerically; for E(L k , kR 1 ) bounded by, let say ǫ = 10 -1 , 10 -2 or 10 -3 , numerical computations show that the L k parameter must approximately satisfy the respective inequalities (with some approximations)

L k ≤ C ǫ kR 1 with C 0.1 = 0.775, C 0.1 = 0.25 C 0.001 = 0.0775. (15)
The √ kR 1 dependency is not surprising when we look at the second order term neglected in the Taylor development of the phase (10).

Note that Inequality [START_REF] Cassereau | Time reversal techniques in ultrasonic nondestructive testing of scattering media[END_REF] is not unconditionally satisfied when L k (= [kr 0 + (kr 0 ) 1 3 -2.5]) grows and k is too large (as discussed in §3.1). To fix the problem one may try to chose the radius of the circle r 0 such that

kr 0 + (kr 0 ) 1 3 -1.5 ≤ C ǫ (kR 1 ) 1 2 .
Asymptotically r 0 must go to 0 like

r 0 ∼ 1 4
R 1 k (we chose ǫ = 10 -2 to simplify). As R 1 is unknown, this formula is not practical, but we can replace R 1 by any estimate from below of this distance. More important, picking a radius which decreases like the square root of the wavenumber will induce a decrease of L k and, consequently, of the localization of the pick, i.e. of the angle of incidence.

In §3.3, we use more terms in the asymptotic series (12) to correct the NMLA filter up to second order.

Second order correction of the filter for a source point

We saw in the previous sections that the good angle localization property obtained for noisy plane wave data does not carry to point source data. We also know that the linearization is responsible for this loss of accuracy. Roughly speaking, the localization is in C k for a plane wave and C ′ √ k for the source point. This section shows that it is possible to identify the curvature from the filtered image of the observable and then use it to improve the angle extraction.

We forget about the noise for now and apply NMLA to the impedant data associated to a source point. We then get (using the maximum test) a first estimate of the true angle

θ 1 = x 1 -x 0 |x 1 -x 0 |
; this first estimation is denoted θ est . The Fourier coefficients βℓ of B k U pt k (θ) are analytically given by (cf. Formula ( 13))

βℓ = 1 2L k + 1 C ℓ (kR 1 ) e -iℓθ 1 .
We shift the phase for convenience :

βℓ e iℓθest = 1 2L k + 1 C ℓ (kR 1 ) e iℓ(θest-θ 1 ) = 1 2L k + 1 C ℓ (kR 1 ) e iℓδθ ,
where δθ = θ estθ 1 is the error we make on the angle. We now use a second order asymptotic from (12),

C ℓ (kR 1 ) ≃ e i ℓ 2 -1 4 2kR 1 1 - ℓ 2 -1 4 (kR 1 ) 2 1 4 , (16)
and obtain

βℓ e iℓθest = 1 2L k + 1 e i " ℓδθ+ ℓ 2 -1 4 2kR 1 « 1 1 - ℓ 2 -1 4 (kR 1 ) 2 1 4 + O( 1 (kR 1 ) 3 ).
If we consider(ℑ stands for imaginary part)

ϕ ℓ = ℑ log( βℓ e iℓθest ) -log( β0 ) ≃ ℓδθ + ℓ 2 2kR 1 , (17)
we see that we get a parabolic drift of the phase. The correction process consists in computing the closest parabola (for instance in the least square sense) of ℓ → ϕ ℓ for |ℓ| ≤ L k . A weighted optimization can be used to take into account that the asymptotic is more accurate for small ℓs. The estimated parabola coefficients δθ est and R est provide a correction on the angle and an estimation of the distance to the source or equivalently of the front curvature. These, in turn can be used to construct a corrected filter :

B cor k U(θ) = 1 2L k + 1 L k ℓ=-L k (F U) ℓ e iℓθ Ĉcor ℓ (-i) ℓ (J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 )) (18)
with Ĉcor

ℓ = e -iℓ " δθest+ ℓ 2 -1 4 2kR est « 1 - ℓ 2 -1 4 (kR est ) 2 1 4
, which gives better localization as the parabolic drift of the phase of βℓ has been removed. Remark that Ĉcor ℓ ≤ 1, and so the uniform stability estimate in k of the norm of B k (7) also holds for B cor k : B cor k ≤ B ⋆ independent of k. This modified corrected filter will yield a better result only if the error induced by approximation ( 16) is small enough to not perturb too much the value of the maximum. Using the approach described in the previous section, we have to make sure that

E ′ (kR 1 , L k ) =    1 2L k + 1 L k ℓ=-L k C ℓ (kR 1 ) - e i ℓ 2 -1 4 2kR 1 1 - ℓ 2 -1 4 (kR 1 ) 2 1 4    ≤ ǫ, (19)
where we used the truncation, now up to second order, of the asymptotic series (12). Again, we performed numerical computations which show that when kR 1 > 4.25 then ( 19) is satisfied when

L k + 1 ≤ C ǫ (kR 1 ) 3 4
, with C ǫ = 1.75, 1.05, 0.59, for ǫ equals to 10 -1 , 10 -2 , 10 -3 . Corrected NMLA will only succeed if

kr 0 + (kr 0 ) 1 3 ≤ C ǫ (kR 1 ) 3 4 .
The radius r 0 must asymptotically behave like

r 0 ∼ C ǫ R 3 4 1 k 1 4
.

The confidence interval

π L k + 1 2
for angle localization (see §3. .

Discretization parameters

We discuss here the discretization of the process. In [3], we proposed to discretize the observable as well as its filtered image by B k using the same 2L k + 1 number of equidistributed samples. In other words, the number of Fourier coefficients of the filter was linked to the discretization. In practice though, we are free to use different numbers. We therefore relax this constraint and discuss here the choice of these parameters. The observable will be sampled with P > 2L k + 1 points while its filtered image will be sampled in Q > 2L k + 1 points.

We show that number P must be chosen large enough for the L k Fourier coefficients to be correctly evaluated, while the number Q must be large enough to capture with sufficient precision the angle corresponding to the maximum of the modulus of the filtered image.

Discretization of the observable

Let P be an integer larger than 2L k + 1, we define

       B P k U(θ) = 1 2L k + 1 L k ℓ=-L k (F P U) ℓ e iℓθ (-i) ℓ (J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 )) où L k = max([kr 0 ] , kr 0 + (kr 0 ) 1 3 -2.5 ), (20)
with

(F P U) ℓ = 1 P P p=1 U(θ p )e -iℓθp , with θ p = 2π(p -1) P .
Operator B P k is a discrete approximation of the filter B k . A straightforward calculation shows that the estimate of the norm of B k remains true for B P k ; indeed, we have by following the inequality of the continuous case

B P k U(θ) ≤   1 2L k + 1 L k ℓ=-L k 1 J 2 ℓ (kr 0 ) + J ′ 2 ℓ (kr 0 ) 1 2   L k ℓ=-L k F P U ℓ 2 1 2
, then, since 2L k + 1 is assumed less than P , the discrete Parseval equality provides

B P k U(θ) ≤   1 2L k + 1 L k ℓ=-L k 1 J 2 ℓ (kr 0 ) + J ′ 2 ℓ (kr 0 ) 1 2   1 P P p=1 |U(θ p )| 2 1 2
, or, bounding the term between hooks by B ⋆ (cf. appendix) and each sample |U(θ p )| by the maximum of |U|

B P k U(θ) ≤ B ⋆ U L ∞ (0,2π) .
To evaluate the number of samples necessary to a good approximation, we split the observable as follows

U(θ) = U P (θ) + U ⊥ P (θ), U P (θ) = M ℓ=-M (F U) ℓ e iℓθ , with M = P -1 2 .
Remark that (as M > L k ) the Fourier coefficients of U P are calculated exactly by our equidistributed quadrature rule: we have (F P U P ) ℓ = (F U P ) ℓ , which implies that B P k U P coincides with B k U P . Thus,

B P k U(θ) = B P k U P (θ) + B P k U ⊥ P (θ) = B k U P (θ) + B P k U ⊥ P (θ) = B k U(θ) + (B P k -B k )U ⊥ P (θ); but since B k U ⊥ P vanishes (as M > L k again), it remains B P k U(θ) = B k U P (θ) + B P k U ⊥ P (θ), whence |B P k U(θ) -B k U(θ)| ≤ B ⋆ sup θ |U ⊥ P (θ)|.
Assuming that P is large enough for the sup of the remaining part to be small, we have that the function

B P k U(θ) is, when U(θ)
is the trace of some perturbed plane waves, the weighted sum of a linear combination of sinc like functions ψ L k (θ-θ p ) plus a perturbation (sum of the noise raised with the error induced by the discretization). It is then possible to proceed again as in the §2.3 and shows that the supremum is reached at a point localized in a confidence interval of size π L k + 1 2 around the the angle of incidence. Given ǫ small, it is therefore sufficient to choose P such that

sup θ |U ⊥ P (θ)| ≤ ǫ|b 1 |.
If U is the observable associated to the plane wave b 1 e -ikr cos(θ-θ 1 ) , this is equivalent to sup θ ℓ>M

(J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 ))i ℓ cos ℓθ ≤ ǫ, M = P -1 2 .
This problem has been the subject of many studies in the framework of the "Fast Multipole Method". For instance, using the techniques developed in [7], it can be shown that when kr 0 is larger than 10 then, taking P > 2M ǫ (kr 0 ) + 1 points with

M ǫ (kr 0 ) = kr 0 + 1 2 3 2 W 2 3πε 2 2 3 (kr 0 ) 1 3 + 2, (21)
The function W (u) is the Lambert function defined as the unique solution of

W (u)e W (u) = u (22)
(this function has a sub-logarithmic behavior when u tends toward infinity). We get

sup ϕ ℓ>M ǫ (kr 0 ) (J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 ))(-i) ℓ e iℓϕ < 2ǫ.
A Numerical study shows that the constant in front of (kr 0 ) 21) is, respectively, 3.24, 49.37 and 156.01 for ǫ = 10 -2 , ǫ = 10 -6 and ǫ = 10 -10 . If we compare with L k in filter [START_REF] De La Bourdonnaye | Numerical simulation of scattering problems with Fourier-integral operators[END_REF], we see that it is always necessary choose M ǫ (kr 0 ) > L k , that is more samples than Fourier coefficients.

1 3 in (
Remark: The same approach can be followed to study the discretization in the point source case. Ones uses in this case truncature error estimates in the Graf'series that have been also studied in the Fast Multipole Method. The estimates for the M ǫ (kr 0 ) are of the same kind but with constants in front of the term in (kr 0 ) 1 3 which depends now on r 0 R , see [6].

Taking the maximum over a finite set of samples

The natural and practical post processing of the NMLA output consist in selecting the supremum of B k U k over a finite number of discretization angles. We are interested here in the evaluation of the number of samples Q needed to capture it accurately. First, we know that before discretization, the size of the confidence interval linked to the truncation to 2

L k + 1 is π L k + 1 2
. So a uniform discretization of size Q should be chosen such that there is least one sample in every interval of this length. We need

Q ≥ 2L k + 1.
Let us proceed for a single plane wave; the discretization of the image of the NMLA operator is

B k U k (θ Q q ) = b 1 Ψ L k (θ Q q -θ 1 ) + δU(θ Q q ), θ Q q = 2π(q -1) Q .
We now discuss and use the function

Ψ L k as in §2.3. For θ Q q outside the interval |θ Q q -θ 1 | ≤ π L k + 1 2
, we have

|B k U k (θ Q q )| ≤ 1 2 |b 1 | + B ⋆ δU L ∞ .
For θ Q q inside this interval, the supremum of our discretization is bounded below by the value taken at the sampled angle the closest to θ 1 , denoted θ Q q 0 . We have

|B k U k (θ Q q 0 )| ≥ |b 1 |Ψ L (δθ) -B ⋆ δU L ∞ , where δθ = θ Q q 0 -θ 1 . We necessarily have |δθ| ≤ π Q .
and the inequality is an equality when the angle θ 1 is exactly at the middle of two samples. When Q is large enough, δθ is close to 0 and so Ψ L (δθ) is close to 1. Indeed, we have

Ψ L (δθ) ≥ sin[(L k + 1 2 )δθ] (L k + 1 2 )δθ ≥ 1 - 1 6 (L k + 1 2 )δθ 2 , which in terms of Q gives Ψ L (δθ) ≥ 1 - 1 6 π(L k + 1 2 ) Q 2 .
Therefore we can choose Q large enough such that , for instance,

Ψ L (δθ) ≥ 3 4 , then, |B k U k (θ Q q 0 )| ≥ 3 4 |b 1 | -B ⋆ δU L ∞ > 1 2 |b 1 | + B ⋆ δU L ∞
and the conclusion follows as §2.3, when the noise level is small enough, see (7).

An extension: the Gaussian weights filter

In the case of multiple plane waves or point sources arriving at the observation with angles θ n , the NMLA output is a sum of functions Ψ L (θθ n ). These functions decrease like C |θ-θn| , which is not fast enough to avoid pollution of the (local) maxima identification for each branch of the data when angles are close.

A possible way to improve this situation is to convolute the result with a smooth window, we used a Gaussian function. It smooths the Ψ L functions, make them decrease faster and preserves the maxima location. The modification is easily incorporated into the filter by a multiplication in the Fourier space

       B g k U(θ) = 1 2L k + 1 L k ℓ=-L k (F U) ℓ ̟ k ℓ e iℓθ (-i) ℓ (J ℓ (kr 0 ) -iJ ′ ℓ (kr 0 )) with L k = max(1, [kr 0 ], kr 0 + (kr 0 ) 1 3 -2.5 ), (23)
with

̟ k ℓ = Ĝk ℓ ν k , ν k = 1 2L k + 1 L ℓ=-L Ĝk ℓ , (24)
where Ĝk ℓ are real positive and even coefficients. We decided on using the Fourier coefficients of a Gaussian function

G k (θ) (τ is a free parameter) Ĝk ℓ = e -1 2 " ℓ κL k " 2 . ( 25
) As Ĝk ℓ < 1, lim L k →∞ ν k = τ 1 τ 0 e -1 2 ξ 2 dξ,
we infer that there exists some number C τ , independent of k such that

sup k max ℓ, |ℓ|≤L k ̟ k ℓ ≤ C τ ,
and get the bound

max θ |B g k U(θ)| < B ⋆ C τ max θ |U(θ)| .
The new Gaussian weighted filter ( 23)-( 24)-( 25) enjoys the same stability properties as the original NMLA filter. We have good control over noise.

If we apply it to a plane wave packet we obtain

B g k U plane k (θ) = P p=1 b p Ψ τ L k (θ -θ p ),
where the new shape function Ψ L k is the convolution of the original one with the Gaussian

Ψ L k ⋆ G k (θ) = Ψ τ L k (θ) =   L k ℓ=-L k Ĝk ℓ   -1 L k ℓ=-L k Ĝk ℓ e iℓ(θ-θp) .
When τ is not too small, each function associated to a particular angle decrease much faster. It allows to separate angle contributions when they are not too close. See the numerical study section or the examples in [2]. A detailed study of function ψ τ L k , as well as the possibility to identify analytically or numerically through non linear optimization several maxima in interfering sums of such functions is on-going work.

Numerical study

Noisy synthetic point source data

We use the exact solution in homogeneous space, the Hankel function. The observation point x 0 is at a distance R = 1 meter from the source point in the direction (cos θ i , sin θ i ). We arbitrarily chose θ i = 157.77 degrees.

In order to add correlated noise, we proceed as follows : We introduce a regular Cartesian grid on top of the observation circle. The grid step is

h = 2π k 1 n g λ
where k is the wave number. We use n g λ = 5 and n g λ = 10 for strong and weak noise.The data values are not taken exactly at the points uniformly sampled on the observation circle but at their closest grid points. We do not interpolate.

The test algorithm for the curvature correction

We assume 50% error on the available D estimate. So we take R := R -= 0.5 meter. We first numerically estimate L 0 as a function of k such that [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF] satisfies

E ′ (kR -, L) ≃ 0.01.
On the wavelength interval we consider we find

L real ≃ 0.93(kR -) 0.763 .
Then we set the radius of the observation circle r 0 to be the unique solution of

L real = kr 0 + (kr 0 ) 1 3 -2.5.
The filter [START_REF] De La Bourdonnaye | Numerical simulation of scattering problems with Fourier-integral operators[END_REF] integer truncation level L k is close to L real . The observation points are uniformly sampled on the circle and then projected on the "noise" grid. We take P = 2M + 2 points where M is given by ( 21) where ǫ = 10 -13 . Then we apply the NMLA (4) discrete filter on this data and select the angle corresponding to the maximum modulus

β p = B P k U k 2π(p -1) P for p = 1, . . . , P
(we take Q = P to simplify).

We then perform the curvature analysis ( §3.3). We compute the phase ϕ ℓ defined by ( 17) being careful to unwrap the possible 2π rotations (we assume the phase cannot jump more than π 2 between two successive ℓ). There are different possible ways to perform the parabolic fit of ℓ → ϕ ℓ . We do as follows : for two values of L e , L e = L k 2 and L e = L k 2 + 1, we find the best least square parabolic fit of ϕ ℓ between -L e and L e . It determines three polynomial coefficients : the zeroth order term is a constant phase ϕ m (L e ), the first order term is the angle correction δθ est (L e ) and the third order term can be used to compute the distance to the source R est (L e ) (see §3.3). We then take the mean value of the two L e , denoted δθ est and R est and use them to correct the filter as in [START_REF] Collino | Analyse microlocale dans le domaine temporel[END_REF]. We thus obtain new corrected value sequence β p . Assuming the maximum in modulus is reached for p = p⋆, the angle θ est = 2π(p ⋆ -1) P + δθ est is the new approximation of the incidence angle while β p ⋆ approximated the asymptotic amplitude of the source point data

β p ⋆ = B est (k) ≃ B exact = e ikD √ - 8iπkD
We present below numerical results comparing this algorithm using NMLA and its Gaussian weight version ( §5) where we use τ = 2 7 .

Numerical results

For fixed k

We fix the wavelength to 2.5 centimeters or kR = 40π ≃ 251.32.. Using the law found in §6. 2 The radius of the observation circle is

r 0 ≃ 14.46 cm.
The grid noise in relative L 2 norm with h = λ 10 is 22%. We apply B P k and its Gaussian weight version (L k = 37) to this data. On Figure 2 we plot β p the output of the filter versus the angle in degrees. The left column is the first NMLA pass, the right is the curvature correction (second NMLA) pass. The top row corresponds to the regular NMLA filter, the bottom row the Gaussian weight version. We notice that all filters peak near the correct angle. Amplitude and angles are not very accurate for the first pass (especially for the non Gaussian weight version) because The plane wave approximation error is too large for this r 0 (81% in relative L 2 norm). The maximum is reached at θ i + 1.05 degree.

The parabolic fit of ( 17) is presented in figure 1: we indeed observe a cloud of points around a parabola. The estimated curvature is R est = 1.1075 meter (10% error). The curvature correction (second pass) on the right is sharper. We get

θ est -θ i = -0.037degree, |B est | -|B exact | |B exact | ≃ -2.2%
while for the Gaussian filter,

θ est -θ i = -0.037degree, |B est | -|B exact | |B exact | ≃ -0.5%,
We remark that the Gaussian weight version (bottom) does a better job on the modulus of the amplitude.

One can improve this result by running the algorithm §6.2 again but with the better estimate we just found for R : R -= 1.10075 meter. The result is shown in figure3. The radius r 0 is then bigger (24.7cm) and the new estimated curvature R est is found to be 1.01 meter (1% error).

Varying k

We study numerically the result of algorithm §6.2 when k takes 300 values such that the source-observation distance R = 1 represents between 10 and 400 wavelengths.

The "noise" grid step is one tenth of the wavelength. The relative noise level on the data varies with k between 8 and 14% (See figure 4,top ).

The curvature estimate R est (k) varies around the exact value R = 1 with an error which decrease with k down to 2% from a maximum at 10% (cf. figure 4, middle). . Middle :

k 2π → |Rest(k)|-R| R . Bottom : k 2π → θ i -θ est (k).
The error on the angle (figure 4 bottom) shows a maximum of 0.2 degree but then decrease to 0.01 degree for high frequencies.

There are no noticeable differences when we use the Gaussian weight filter except on the modulus of the amplitude. Figure 5 show the relative error for this quantity which is underestimated by 0.8% for the standard NMLA (top) and 0.15% for the Gaussian weight filter (bottom). N(L, kr 0 ) defined as

N 2 (L, x) = 1 (2L + 1) 2 L n=-L 1 J 2 n (x) + J ′ 2 n (x)
.

A straightforward calculation gives N 2 (1, 0) = 1 and, due to the Hansen Equality and its variant

∞ n=-∞ J 2 n (x) = 1, ∞ n=-∞ J ′ 2 n (x) = 1 2 ,
Cauchy Swartz provides

N 2 (L, x) > 2 3 .
Now, looking at the curves L → N(L, x) for a different positive values of x, the numerical computations depicted in figure A1 show that this sequence is at first decreasing on some interval [1, L(x)] then increasing for L in [L(x), +∞[, the increasing rate being more than exponential when L overpasses a value slightly larger than x. In Figure A2, we show the numerical behavior of functions x → N(L, x) for different values of L. It appears that the minimum over all L of N(L, x) remains between 0.8 and 1. It is possible to approximate numerically an optimal law L opt (x) for the L which minimizes N L (x). For x large enough (> 5) and less than 8000 we find

L opt (x) ≃ x + x 1 3 -2.5.
So the optimal integer L for a given value of x will be Remark Long technical arguments we skipped here, give evidence that a better approximation of L opt (x) is

L ⋆ (x) = max 1, [x], [x + x 1 3 -2.5] . 0.
L opt (x) ≃ max(1, x + 1 32 W 2 x 4 8 3π 3 1 3 x 1 3 ),
where u → W (u) is the Lambert function, see [START_REF] Hugonnet | High resolution radon: a review[END_REF].

Appendix A.2. The stability result

We prove the Theorem Appendix A.1 Let L be a positive integer and x a positive number; let

Ñ2 (L, x) = 1 (L + 1 2 ) 2 L n=-L 1 J 2 n (x) + J ′ 2 n (x)
, then there exists a pure constant C such that

sup n (x,L)/0≤x≤L≤x+x 1 3 o Ñ2 (L, x) ≤ C.

Proof :

As a preamble, we recall that since J -ℓ (x) = (-1) ℓ J ℓ (x), see [13, (2) page 15], we have as well

Ñ2 (L, x) = 1 (L + 1 2 ) 2 L n=0 ε ℓ 1 J 2 n (x) + J ′ 2 n (x)
, where ε 0 = 1, and ε ℓ = 2 when ℓ = 0. Thus, we can work with Bessel functions of non negative order.

For technical reasons, we consider separately the cases x "small" and x "large". More precisely, let x 0 be large enough to satisfy

x + x 1 3 -x 2 3 + 1 ≤ √ x 2 -1, for all x > x 0 . (A.1)
Numerically, we find x 0 = 4.6237... First, we consider the case x < x 0 ; the function

x → (x + 1 2 )(J 2 n (x) + J ′ 2 n (x)
) is non negative and does not vanish for x > 0. Moreover, it tends to 2 π when x goes to infinity since, cf. [13, page 199]

J n (x) = 2 πx cos x - π 2 (n + 1 2 ) 1 + O(x -1 ) , J ′ n (x) = 1 2 (J n-1 (x) -J n+1 (x)) = 2 πx sin x - π 2 (n + 1 2 ) 1 + O(x -1
) .

We deduce that there are some constants

β n such that 1 (x + 1 2 )(J 2 n (x) + J ′ 2 n (x)) ≤ β n ,
and, consequently

Ñ2 (L, x) ≤ x + 1 2 (L + 1 2 ) 2 L n=0 ε n β n . Now, as x < x 0 , we have if L ≤ [x + x 1 3 ] ≤ [x 0 + x 1 3 0 ], Ñ2 (L, x) ≤ x 0 + 1 2 [ 1 2 ] 2 [x 0 +x 1 3 0 ]+1 n=0 ε n β n = C ste .
If the norm Ñ2 (L, x) diverged over x ≤ L ≤ x + x 1 3 , it can only be for x large. Consider now x large, i.e. x > x 0 . Bounding above being easier than bounding below, we use the identity (Y n (x) is the Neumann Function, see [13])

1 J 2 n (x) + J ′ 2 n (x) = Y 2 n (x) + Y ′ 2 n (x) (J n (x)Y n (x) + J ′ n (x)Y ′ n (x)) 2 + (J n (x)Y ′ n (x) -J ′ n (x)Y n (x)) 2 .
We recognize in the second term inside the denominator the square of the Wronksian

J n (x)Y ′ n (x) -J ′ n (x)Y n (x) = 2 πx , (A.2)
where Y n (x) is the Neumann Function, [13, pp 64]; we then have

1 J 2 n (x) + J ′ 2 n (x) ≤ π 2 x 2 4 Y 2 n (x) + Y ′ 2 n (x) .
To get rid of the derivative, we use

Y ′ n (x) = 1 2 (Y n-1 (x) -Y n+1 (x)) , n = 0, Y ′ 0 (x) = -Y 1 (x) and therefore 1 J 2 n (x) + J ′ 2 n (x) ≤ π 2 x 2 4 Y 2 n (x) + 1 2 Y 2 n-1 (x) + 1 2 Y 2 n+1 (x) ,
from which, we deduce (the constant π 2 is not optimal)

Ñ2 (L, x) ≤ π 2 x 2 (L + 1 2 ) 2 L+1 n=0 Y 2 n (x) (A.3)
We are going bound above the right hand side of (A.3). We split the summation over n into four parts: n = 0, then n = 1, . . . Q, then n = Q + 1, . . . L -1 and at finally n = L, L + 1 where Q is chosen such that

x + x 1 3 -x 2 3 ≤ Q < √ x 2 -1.
Such an integer does exist because the difference between the upper and lower bounds is always larger than 1 when x > x 0 (x 0 was precisely chosen for that). In the following, we will often use the square of the modulus of the Hankel function J n (x) + iY n (x) since one controls rather well its behavior, thanks to the Nicholson formulae, [13, page 444]

J 2 n (x) + Y 2 n (x) = 8 π ∞ 0 K 0 (2x sinh t) cosh(2nt) dt. (i) For n = 0 : as x → J 2 n (x) + Y 2 n (x) is a decreasing function, [13, page 446]. So, for x ∈ [x 0 , ; ∞] Y 2 0 (x) ≤ J 2 0 (x 0 ) + Y 2 0 (x 0 ) = C 0 . (A.4) (ii) For n = 1, . . . Q: we use [13, page 447], when x ≥ n ≥ 1 2 J 2 n (x) + Y 2 n (x) ≤ 2 π 1 √ x 2 -n 2 ,
and therefore

Q n=1 Y 2 n (x) ≤ 2 π Q n=1 1 √ x 2 -n 2 , but, since the sequence n → 1 √ x 2 -n 2 is increasing, we have Q n=1 J 2 n (x) + Y 2 n (x) ≤ 2 π 1 x 2 -Q 2 + Q 0 dν √ x 2 -ν 2 , or Q n=1 J 2 n (x) + Y 2 n (x) ≤ 2 π 1 x 2 -Q 2 + arcsin ν x Q ν=0 ≤ 2 π 1 x 2 -Q 2 + π .
Now, one's recall that the Q was chosen such that Q ≤ √ x 2 -1 and therefore The behavior of function ψ(x) for x large can be deduced from the approximations of Bessel's Functions which argument is closed to its order. From the estimates given in [9], we deduce (painfully) that when κ is given and x is large,

1 √ x 2 -Q 2 ≤ 1; we deduce Q n=1 J 2 n (x) + Y 2 n (x) ≤ 2 π (1 + π) = C 1 . (A.5) (iii) For n = Q + 1, . . . L -1: function ν → J 2 ν (x) + Y 2 ν (x) is increasing (cf. Nicholson formulae), we thus have L-1 n=Q+1 Y 2 n (x) ≤ L-1 n=Q+1 J 2 n (x) + Y 2 n (x) ≤ (L -Q) J 2 L-1 (x) + Y 2 L-1 (x) . But if L ≤ x + x 1 3 + 1, we have L-1 n=Q+1 J 2 n (x) + Y 2 n (x) ≤ (L -Q) J 2 x+x 1 3 (x) + Y 2 x+x 1 3 (x) = (L -Q)ψ(x).
J x+κx 1 3 (x) ≃ 2 x 1 3
Ai(2

1 3 κ), Y x+κx 1 3 (x) ≃ 2 x 1 3
Bi(2

1 3 κ),
Where Functions Ai and Bi are the Airy Functions . Thus (cf. Fig A3 ), there exists a constant C 2 such that sup

x>x 0 x 2 3 J 2 x+x 1 3 (x) + Y 2 x+x 1 3 (x) ≤ C 2 ,
and therefore

L-1 n=Q+1 J 2 n (x) + Y 2 n (x) ≤ (L -Q) x 2 3 C 2 ≤ C 2 , (A.6)
The last equality coming from the bound

L -Q ≤ x + x 1 3 -(x + x 1 3 -x 2
3 ) = x 2 3 . (iv) For the two last terms of the series L, L + 1 : in order to use the previous estimate, we move back twice on the order of the Bessel functions . Since

Y n+1 (x) = -Y n-1 (x) + 2n x Y n (x), we have 1 2 Y 2 n+1 (x) ≤ Y 2 n-1 (x) + 4n 2
x 2 Y 2 n (x), and, by successive bounding, we get

Y 2 n (x) + Y 2 n+1 (x) ≤ 2 1 + 8n 2 x 2 1 + 4(n -1) 2 x 2 max Y 2 n-2 (x), Y 2 n-1 (x) .
and therefore

Y 2 n (x) + Y 2 n+1 (x) ≤ 2 1 + 8n 2 x 2 1 + 4n 2 x 2 Y 2 n-1 (x) + J 2 n-1 (x)
and that, always thanks to the increasing property of the sequence n → J 2 n (x)+Y 2 n (x). We are brought to the previous bound. Thus, we have, when L < x + x

1 3 L+1 n=L Y 2 n (x) ≤ 2 1 + 8L 2 x 2 1 + 4L 2 x 2 (J 2 x+x 1 3 (x) + Y 2 x+x 1 3 (x)),
or also, when

x 0 ≤ x ≤ L ≤ x + x 1 3 , L+1 n=L Y 2 n (x) ≤ 2 1 + 8(x + x 1 3 ) 2 x 2 1 + 4(x + x 1 3 ) 2 x 2 C 2 x 2 3 ≤ C 3 . (A.7)
We gather all our results, when x + x 1 3 > L > x > x 0 , there exists four constants such that

Ñ2 (L, x) ≤ π 2 x 2 (L + 1 2 ) 2 (C 0 + C 1 + C 2 + C 3 ) ≤ π 2 (C 0 + C 1 + C 2 + C 3 )
, and since Ñ2 (L, x) is uniformly bounded over x ≤ x 0 and L ≤ [x 0 + x

1 3 0 ], Q.E.D.
(ii) time reverse this signal and re-emit it in the medium.

Then, physics observes and maths proves that the time reversed signal will exhibit refocusing properties towards the scatterers or the source points depending on the set up of the experiment.

We propose a new interpretation of time reversal using the High Frequency (HF) asymptotic model of wave propagation, also known as Geometric Optics (GO). The HF model only makes sense with smooth deterministic media when the characteristic wavelength of the signal is smaller than the variation scale of the acoustic speed. We will therefore NOT be concerned with the widely studied "super-resolution" effect of the focusing observed in random media [9].

The goal of the paper is twofold :

• Provide a non technical description of the Numerical MicroLocal Analysis algorithm (NMLA) [7] and its recent improved second order version [6]. We apply it to analyze the underlying GO ray directions of the signal recorded in step (i). The NMLA algorithm works in the frequency domain. We therefore use classic pre and post processing techniques to extend it to time domain and also to treat a signal recorded on linear arrays of transducers. We also show that the use of time domain data simplify multiple rays discrimination and that frequency redundant information provides a traveltime estimate from the source.

• Demonstrate numerically that the output ray direction and traveltime are accurate enough, at least with synthetic noisy data, to shoot the ray back to its unknown source location. The notion of "shooting back" implies that the traveltime, the actual propagation time, is reversed in the same way as in step (ii) above.

In order to make the paper self contained, it is organized as follows :

Section 2 recalls the HF model derivation and formulates the HF asymptotic time reversal source discovery inverse problem.

Section 3 reviews the refined NMLA algorithm [6] and its extension to time domain data recorded on a linear array of "transducers".

Section 4 presents numerical source discovery from synthetic data.

2

2 Forward and inverse Problem

High frequency acoustic wave equation

We start from the classical acoustic wave equation in the time domain. The wave field u(x, t) is generated by a finite time wavelet emitting from a source point :

(P )    1 c 2 (x) ∂ 2 tt u(x, t) -∆u(x, t) = δ(x -x s ) W ωs (t), ∀(x, t) ∈ R 2 × [0, T ) u(x, t = 0) = ∂ t u(x, t = 0) = 0, ∀x ∈ R 2 , (1) where
• [0, T ) is the simulation time (one may consider T = +∞)

• x → c(x) is the acoustic speed law function on R 2 . As we deal with HF GO models, we will consider smooth speeds.

• x s is the source point location

• t → W ωs (t) is the source time wavelet.
In all the sequel the radian frequency ω will be the dual variable of time in the Fourier transform denoted .

û(x; ω) = ∞ 0 u(x, t) e iωt dt. (2)
The frequency content of the solution depends on the source wavelet's. HF wave fields have strong HF content. It can be tuned using the "characteristic pulsation (or frequency)" ω s which parametrizes W ωs (t). More precisely, we set

W ωs (t) = W (ω s t),
where W (t) is a shape function with support in [0, 1]. Assuming further that W is absolutely integrable and has zero mean, the following properties can be checked :

Support (W ωs ) = [0, 1 ω s ]. ( 3
)
Ŵωs (ω) = 1 ω s Ŵ ω ω s . ( 4
)
For ǫ > 0, there exists two positive constants γ + ǫ and γ - ǫ such that

ω / ∈I | Ŵωs (ω)| 2 dω ≤ ǫ | Ŵωs (ω)| 2 dω, I = [γ - ǫ ω s , γ + ǫ ω s ] (5)
So when ω s tends to infinity, (3) shows that the source impulse lapse shrinks and (4-5) that the frequency content becomes large and is shifted to infinity.

Frequency domain and ansatz

The HF model uses a formal asymptotic representation of the solution, where the unknowns are independent of the radian "frequency" ω. It is easier to derive and interpret it in the frequency domain :

Assume T = +∞ and Fourier transform (see ( 2)) equation ( 1) :

- ω 2 c 2 (x) û(x; ω) -∆û(x; ω) = δ(x -x s ) Ŵωs (ω). (6)
We consider that c(x) = v is constant outside a compact in Ω, and add the Sommerfeld radiation boundary condition at infinity

lim |x|→∞ ∂ |x| û(x; ω) - iω v û(x; ω) = 0, (7)
to ensure uniqueness [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF].

In the homogeneous case (c(x) = v constant everywhere), the analytic solution of (6-7) is given by û

(x; ω) = i 4 H (1) 0 ω |x -x s | v Ŵωs (ω), (8)
where H

(1) 0 (u) is the Hankel function of first type and order 0. Using the following asymptotic estimate (cf. [START_REF] Watson | avenue de l'Europe -38334 Montbonnot Saint-Ismier Centre de recherche INRIA Lille -Nord Europe : Parc Scientifique de la Haute Borne -40, avenue Halley -59650 Villeneuve d'Ascq Centre de recherche INRIA Nancy -Grand Est : LORIA, Technopôle de Nancy-Brabois -Campus scientifique 615, rue du Jardin Botanique -BP 101 -54602 Villers-lès[END_REF]) 8), we obtain the asymptotic approximation of the solution û(x; ω) ≃ A(x; x s )e iωφ(x,xs) Ŵωs (ω)

H (1) 0 (x) ∼ x→∞ 2 πx e i(x-π 4 ) , in (
√ -iω , (9)
where the large parameter is the number of wavelengths between the source and the observation point (the product ω |x-xs| v

) and A and φ are the GO solution (see next section) :

A(x, y) = v 8π|x -y| , φ(x, y) = |x -y| v . ( 10
)
The "ansatz" idea, apparently going back to a 1911 paper of Debye on "the article of Sommerfeld" (see [1]), consists in generalizing the a priori form (9) to the general heterogeneous case and search for a formal asymptotic solution of the Helmholtz equation in the form ûray (x; ω) ≃ e iωϕ(x)

∞ j=0 A j (x) (i ω) j , (11)
(we assume the source wavelet and the half primitive (symbol √ -iω) have been removed, see section 3.1).

In general, only the first A 0 term of the amplitude series is used. Oscillations are formally enforced in the phase and travel at traveltime ϕ. The radian frequency ω is the large parameter.

Geometric optics -Ray tracing

Assuming a solution taking the ansatz form (11) locally away from the source. One can plug it into ∆û(x; ω) + ω 2 c 2 (x) û(x; ω) = 0, (12) and expand in inverse power of ω. Canceling the coefficients of the first and second order term of the series gives respectively the Eikonal equation and amplitude equations

   |∇ϕ(x)| = 1 c(x) 2∇ϕ(x) • ∇A 0 (x) + A 0 (x)∆ϕ(x) = 0. ( 13
)
These equations are local but global solutions are constructed using a Lagrangian method called Ray Tracing.

The traveltime ϕ is computed along the integral curves of the gradient called rays and denoted s → y(s, y 0 ). Each ray may be parameterized by its initial position y 0 and assuming the traveltime and its gradient are known at this location one solves the following system of ODEs :

        
where • y(s, y 0 ) is the ray shot from y 0 and s is a parameterization of the curve.

• p = ∇φ(y) is called the "slowness" vector and gives the ray direction.

• n = 1 c is called the slowness of the medium. When the velocity is homogeneous (c(x) = v) around a source point, the analytic asymptotic form (10) may be used to define the initial conditions : For τ small enough the level curve W C = {y 0 , s.t. φ(y 0 ) = |y 0 -xs| v = τ } is a circle around the source point and p(0, y 0 ) = ∇ϕ(y 0 ) = (y 0 -xs) v |y 0 -xs| . Letting τ go to 0 we consider that the HF asymptotic form of a source point corresponds to rays leaving the source point x s isotropically with the same amplitude. We then switch the parameterization of the rays to the initial take off angle y 0 → θ 0 . The global HF solution is formed of all rays; this is the solution, for all θ 0 ∈ [0, 2π[, of

           d ds y(s, θ 0 ) = ∇ϕ(y(s, θ 0 )) := p(s, θ 0 ), y(0, θ 0 ) = x s d ds p(s, θ 0 ) = 1 2 ∇n 2 (y(s, θ 0 )), p(0, θ 0 ) = s0 v d ds ϕ(y(s, θ 0 )) = n 2 (y(s, θ 0 )), ϕ(x s ) = 0, (15)
where s0 = (cos θ 0 , sin θ 0 ).

Ray Tracing (14-15) can be extended to compute solutions of the amplitude equation (see [5] for instance).

In smooth heterogeneous media, the Lagrangian traveltime function remains smooth but it is quite common for rays to cross. Then, solutions to the Eulerian equation ( 13) become multivalued. In the simplest case, the HF ansatz remains a superposition/sum of elementary (11) ansatz. See [START_REF] Engquist | Computational high frequency wave propagation[END_REF] [5] for more details and references on multivalued solutions. In this paper and also for the NMLA algorithm we will always assume that the Helmholtz solution has a general HF representation of the form

û(x; ω) ≃ N n=1 A 0,n (x)e iωϕn(x) , ( 16
)
where N is the number of rays going through x and n the subscript labeling the different branches of the solution. We kept only the zeroth order term of the amplitude series of (11).

Source point discovery with HF time reversal

We now propose an HF local interpretation of the time reversal algorithm (i) -(ii) :

(i) HF Record an acoustic signal around some observation point x 0 and for "some" time; this will be our data. Then, identify at least one significant ray going through x 0 , its traveltime and slowness vector. Let us denote them ϕ data (x 0 ) and p data 0 = ∇ϕ data (x 0 ).

(ii) HF Time reverse this signal and re-emit it in the medium using the HF model. This is accomplished by changing ds in -ds [START_REF] Cassereau | Time reversal techniques in ultrasonic nondestructive testing of scattering media[END_REF] and use the initial conditions obtained in (i) HF for this particular ray :

         d -ds y(s) = ∇ϕ(y(s)) := p(s), y(0) = x 0 d -ds p(s) = 1 2 ∇n 2 (y(s)), p(0) = p data 0 d -ds ϕ(y(s)) = n 2 (y(s, θ 0 )), ϕ(x s ) = ϕ data (x 0 ). (17)
The traveltime ϕ(y(s)) will then decrease and we stop integrating [START_REF] Collino | Peut-on obtenir des amplitudes correctes avec les quations paraxiales ?[END_REF] when it gets to 0. If step (i) HF yields a reliable HF component estimate, we have time reverted the ray and y(s) has reached the source x 0 .

There are no theoretical or numerical difficulties with step (ii) HF . The next section presents the NMLA algorithm as a way to invert the ansatz ( 16) in order to perform step (i) HF

Numerical MicroLocal Analysis [7] [6]

We assume frequency domain data û(x, ω) is available locally on a circle around an observation point x 0 . Appendix A explains in details how this format can be generated from time domain data recorded on a linear array1 .

NMLA does not directly invert ( 16) but a linear equation obtained through the successive steps described in the following subsections.

Source deconvolution

When the wavelet W ωs is known, one can perform the following source deconvolution on the data

û(x; ω) ← √ -iω Ŵωs (ω) û(x; ω). (18)
We will keep the same notation for the data. See appendix B on getting rid of this step when the source wavelet is unknown.

Plane wave approximation

We linearize the traveltime around x 0 . In the GO regime, we expect the phase and amplitude to be smooth real with variations on a scale larger than ω. Using a Taylor expansion on the traveltime and within a few wavelength of x 0 we have

ϕ(x) ≃ ϕ(x 0 ) + (x -x 0 ) • ∇ϕ(x 0 ) + 1 2 (x -x 0 ) T Hϕ(x 0 ) (x -x 0 ) + ... (19)
The usual first order "plane wave" approximation gives, for N = 1 in (16),

û(x; ω) ≃ B(x 0 )e iω(x-x 0 )•∇ϕ(x 0 ) , (20)
where

B(x 0 ) = A 0 (x 0 )e iωϕ(x 0 ) .
In the more general case where N rays cross at x 0 we assume û can be identified as the N -rays ansatz

û(x, ω) ≃ N n=1 B n (x 0 )e iω(x-x 0 )•∇ϕn(x 0 ) . (21)

Relaxation

Equation ( 21) is non-linear, remember the unknowns are the (B n , ∇ϕ n ) plus N the number of rays. As the norm of the gradient of the traveltimes is given by the Eikonal equation,

∇ϕ n (x 0 ) = 1 c(x 0 )
for all n.

We instead consider the linear problem of finding the amplitude distribution

{B m , m = 1..M } over a discretization of the set of 2-D unit vectors, { dm = m = 1..M } with dm = (cos θ m , sin θ m ) θ m = m 2 π M , m = 1..M (22)
The . notation will always refer to unit vectors. The relaxed version of ( 21) is

û(x, ω) ≃ M m=1 B m (x 0 )e i ω c(x 0 ) (x-x 0 )• dm . (23)
The size of the discretization M is discussed later.

The observable and the Linear system

In our first paper [7] the data û was sampled on a scaled circle around x 0 using the same discretization [START_REF] Hugonnet | High resolution radon: a review[END_REF].

U α,m = û(x 0 + r dm ; ω) r = αc(x 0 ) ω |s| = 1.. (24)
This formulation had stability issues and a L 2 regularization was used in [7]. A L 1 regularization is proposed in [START_REF] Landa | Discovery if point sources in the hemlhotz equation posed in unknown domains with obstacle[END_REF]. The new version proposed in Figure 1: Geometry of the observable support, a sample point in the direction s [6] still uses the circle geometry but the data is sampled from the impedant quantity c(x 0 ) iω ∂ û ∂r + û. The resulting NMLA filter is stable independently of its parameters. The new observable is

U m = c(x 0 ) iω ∂ û ∂r (x 0 + r dm ; ω) + û(x 0 + r dm ; ω), r = αc(x 0 ) ω . (25)
Modifying [START_REF] Fomel | Applications of plane-wave destructor filters[END_REF] according to this new observable and setting the vectors

U α = {U m } B α = {B m } m = 1, . . . , M.
and the matrix

K α = {( s • dm + 1) e iα d m• dl } l, m = 1, . . . , M. (26)
The linear system to be inverted is

U α = K α B α . (27)

K α pseudo-inverse and stability result

The spectrum of the linear operator K α and its continuous version have been carefully studied in [START_REF] Collino | Analyse microlocale dans le domaine temporel[END_REF]. In particular it can be diagonalized on the Fourier basis {e l (s) = 1 √ 2π e ilθ s } l∈Z and its eigenvalues depend on Bessel functions

K α e l = D l (α) e l , D ℓ (α) = 2πi ℓ (J ℓ (α) -i J ′ ℓ (α)); J ℓ (α)
is the Bessel function of order ℓ and argument α, it decreases more than exponentially for ℓ above a threshold which depends on α see figure 2. K α is not invertible in any Sobolev norm but a stable and normalized pseudo-inverse can be computed as follows

{B m } := 1 2 L(α) + 1 F -1 ({ βℓ }), βℓ = H ℓ F({U m }) ℓ (28)
where :

• F is the M discrete Fourier transform on our discrete set of directions [START_REF] Hugonnet | High resolution radon: a review[END_REF] .

• H ℓ = D -1 ℓ
if |ℓ| < L(α) and 0 else.

• L(α) corresponds to a truncation threshold. It is shown in [6] that if

L(α) = min{α, α + α 1/3 -2.5} (29)
then the norm of the normalized inverse operator defined by ( 28) is bounded by 1, independently of α:

sup m |B m | ≤ sup m |U m |
• The sampling M must be larger than L(α).

• The new impedant observable [START_REF] Landa | Discovery if point sources in the hemlhotz equation posed in unknown domains with obstacle[END_REF] prevents the eigenvalues to accidentally vanish for ℓ < L(α) (note that J ℓ (x) and J ′ ℓ (x) have no common zeros, [START_REF] Watson | avenue de l'Europe -38334 Montbonnot Saint-Ismier Centre de recherche INRIA Lille -Nord Europe : Parc Scientifique de la Haute Borne -40, avenue Halley -59650 Villeneuve d'Ascq Centre de recherche INRIA Nancy -Grand Est : LORIA, Technopôle de Nancy-Brabois -Campus scientifique 615, rue du Jardin Botanique -BP 101 -54602 Villers-lès[END_REF]).

Exact solutions for plane wave (linear) data

When the data is a pure plane wave, i.e. linearization is exact, an analytical form of the output of NMLA is available. We denote d the direction of the plane wave (θ d the corresponding angle) and A the amplitude. The asymptotic HF content consists in one ray going through the observation point x 0 with the same direction. The observable data [START_REF] Giladi | A hybrid numerical asymptotic method for scattering problems[END_REF] is given by

• We see that the inverse ( 28) is normalized such that the B m of maximum modulus matches the amplitude of the plane wave.

• The larger the L(α) the more peaked the Sinc like function S α will be.

• We note that for perfect plane wave data, a perfect ray direction approximation can be achieved by increasing M the number of sampled directions.

• As we will see in the next section this is not the case when the asymptotic wavefront has curvature or just adding noise the maximum of the peaked function is moved and a confidence interval can be estimated thanks to the stability result [6].

Second order curvature correction

As can be seen in Figure 8, L(α) and thus α (see [START_REF] Nicolay | Gaussian beam decomposition of high frequency wave fields[END_REF]) bounds the number of Fourier modes and limits the accuracy as we hope to reconstruct Dirac masses in the correct ray direction. The last plot also shows that we are limited in increasing α because of both the plane wave approximation and the error linked to the second order terms in the Taylor expansion [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. We did not manage to include this term in the NMLA filter nor to use the transport equation (13). That seems to remain an open problem.

Instead, we first assume there is only one ray (N = 1). This is not a very restrictive assumption as both time and direction windows can be used to select a particular ray, see Appendix C. We then make the simplest non linear approximation on the HF asymptotic solution : we look for locally constant curvature as if the wavefield solution were the fundamental solution G(x, x s ) = i 4 H

(1) Then, instead of the plane wave approximation, the NMLA data is first approximated locally by a constant curvature wavefield G(x, x s ). Its argument xx s is decomposed into the sum of the main propagation direction (xx 0 = d d) plus x 0x s a vector turning around x 0 : x 0x s = rs . Figure 4 present this decomposition.

0 (ω |x-xs| v ) (check ( 8
U α (s) ≃ 8π d c(x 0 ) A 0 (x 0 ) e iω(ϕ(x 0 )-d) i 4 H (1) 0 ω c(x 0 ) |d d + rs| . ( 31
)
The high frequency amplitude and phase above have been corrected as the global wavefield is not necessarily globally a point source solution in an homogeneous medium. We use the known HF asymptotics of the Hankel function, see (10). Inspired by asymptotic expansions from the Fast Multipole Method applied to electromagnetics [12] [13], one can then consider r as a perturbation in the argument of the Hankel function and expand the argument. Using technical asymptotic expansions, we show in [6] that a second order accurate approximation of the β Fourier modes of the NMLA filter ( 28) is βℓ ≃ Ae iωϕ(x 0 ) e i(ℓθ d+ (ℓ 2 -1

) 2γ

) ,

where γ = ω d c(x 0 ) is a new large parameter and θ d is the angle associated to d. This is to be compared to the exact solution for one plane wave linear data [START_REF] Hodgkiss | Phase conjugation in the ocean: Experimental demonstration of an acoustic time-reversal mirror[END_REF]. We see that we get a γ (curvature dependent) second order correction. From there we simply remark that

1 i log( βℓ β0 ) ≃ ℓθ d + (ℓ 2 -1 4 ) 2γ ) ( 33
)
is a parabolic function of the mode number ℓ and the coefficients θ d and 1 2γ can be fitted by a least square approximation, see figure 5. We summarize the 2nd order correction of the NMLA algorithm (see figure 9 for a numerical illustration) :

1. Get an estimate θ m ⋆ of θ d over the set of all discretized directions [START_REF] Hugonnet | High resolution radon: a review[END_REF] using [START_REF] Porat | Analysis of the asymptotic relative efficiency of music algorithm[END_REF] output :

m ⋆ such that |B m ⋆ | = max m |B m |. (34)
2. Using (33), we perform a parabolic fitting of

1 i log( βℓ β0 e -iℓθ m ⋆ ) ≃ ℓδθ + (ℓ 2 -1 4 ) 2γ . (35)
The parabola coefficients are δθ = θ dθ m ⋆ , the angle correction and the ℓ quadratic term coefficient which depends on the curvature of the front through γ.

Correct NMLA amplitudes

{B ′ m } := F -1 ( βℓ e -i(ℓδθ+ (ℓ 2 -1 4 ) 2γ ) }), (36)
and get the angle corresponding to the maximum of the modulus of the B ′ m 's. Figure 9 illustrates the improved angle correction.

Traveltime estimation

We first notice that manual or automatic (most energetic for instance) picking of an event on a seismogram gives an estimate of the traveltime, assuming of course that the source W ωs (t) is known. The wavelet is however not a perfect wavefront and may spread over a given period of time (around .1s in our numerical illustration figure 17. NMLA combined with the use of time data can provide an accurate traveltime estimation. The idea is to use the implicit redundancy of the HF traveltime over the frequency range. More precisely we have, again assuming perfect asymptotic one ray ansatz data, and performing the curvature correction of section 3.9 to improve the accuracy :

B m ⋆ (ω) ≃ A 0 exp i ω φ(x 0 ) . (37)
The ω dependence is explicit in the notation since we are going to use the variation with respect to frequency :

φ(x 0 ) ≃ ∂ ω B m ⋆ i B m ⋆ . (38)
The ω derivative is obtained by applying the same NMLA procedure and post processing to the ω derivative of the data. The derivative data can be computed using finite difference or as we derive the frequency domain from the time domain by computing the Fourier transform of t u(t, x 1 , 0).

In figure 6, using the case study of section 4.4, we show that this method is more accurate than picking by sliding the center t 0 of a time window around the correct traveltime and plotting the estimated traveltime τ = φ(x 0 ). As long as the time window contains a significant signal the estimated traveltime is constant and correct and therefore does not depend on the picking time as expected. We generate synthetic data using the following procedure :

• A standard FDTD scheme (second order in time and space) is used to simulate the time domain propagation from one or more point sources.

• A "Ricker" (Gaussian derivative) is the source wavelet.

• We use a perfect reflecting boundary condition ∂ ∂ν u = 0 at the bottom of the domain and Higher Order Absorbing Conditions, [START_REF] Collino | Conditions aux limites absorbantes d'ordre élevé pour les modèles de propagation d'onde. problème des domaines rectangulaires[END_REF], elsewhere.

• We record over time a seismogram (just u) at all discretization points of the domain.

• Once the observation time t 0 and surface position x 0 are picked (see figure 17), a smooth window in space and time can be used to single out the event.

• We Fourier transform all the time lines, then choose a frequency ω.

• The NMLA frequency data is then computed using the one way method described in appendix A.

Homogeneous medium 4.3 one source point

Figure 7 depicts the seismogram recorder over time (vertical axis). Each vertical line is the signal t → u(t, x) recorded at a grid point x on the surface (horizontal axis). The line horizontal deviation is proportional to u thus giving an idea of the succession of waves (events) reaching the surface. In this case there is only one source points but two events due to the reflecting boundary condition at the bottom. Our interface allows to select a time t 0 and a position x 0 (red lines) where we want to perform the NMLA. A window in time can be used to select only one event corresponding to one wavefront and one ray. Figure 8 shows the impact of α value on the basic NMLA filter [START_REF] Porat | Analysis of the asymptotic relative efficiency of music algorithm[END_REF]. If α is too small, we do not have enough Fourier modes to focus, note the angle is still correctly estimated for one ray, however for more than one ray (as in the next section) α = 10 would not allow to observe separate maxima when the ray directions are close. If α is too large, then the curvature error induced by the plane wave approximation pollutes the result. Figure 9 shows the output of the NMLA procedure [START_REF] Porat | Analysis of the asymptotic relative efficiency of music algorithm[END_REF] and also the second order correction of section 3.7. After correction the B ′ amplitudes are better focused around the correct ray angle and the corrected angle is very accurate. Table 1 quantifies the error on the angle and the source point discovery. Local Radon transform (LSS) completely fails while NMLA does a better job than plane wave destructor (PWD) (references on these two methods in section 4.6).

Figure 10 shows the error on angle and traveltime when applying NMLA second order to several position on the surface y = 1.

Method

Two source points

Basic NMLA is able to recover well separated multiple rays. We treat here the case of two but it works equally well on more see the original paper [7].

Figure 11 shows the recorded seismogram for two source points. The picked observation point is at equal traveltime from the source points. Thus, the time window will contain two wavefronts and two rays.

Figure 12 shows the NMLA basic filter result on the left and the application of the Gaussian window convolution (45) (see appendix C) which enables to separate the rays if they are not too close and get the β modes associated to only one ray. Then, the curvature correction and traveltime estimation (section 3.7 and 3.8) can be used.

Homogeneous medium Noisy data -Stability

These results illustrate the stability claim on the algorithm (section 3.5) for the two source point case. The NMLA filter (28) easily filters out white noise (13) as it is precisely built to single out coherent signals. In order to test correlated noise, we added to the NMLA input data also collected on circle around the observation point but with a different radius, see figure 14.

Heterogeneous medium one source point

We now consider the smooth speed velocity distribution c of figure 15. The + sign is the source location and our observation point is the cross x at the surface. Table 2 displays the ray direction found with several methods. NMLA basic and the second order correction. The Plane Wave Destructor Method PWD [START_REF] Fomel | Applications of plane-wave destructor filters[END_REF] and the Local Slant Stack LSS (Radon transform) [10] [22] are popular methods in Geophysics and radar processing but are very sensitive to the necessary windowing. Note that PWD is first order by construction; the Radon transform can be generalized to identify parabolas [8] but we did not tested this against NMLA second order. Both PWD and LSS can be considered as image processing methods and do not use at all the wave models.

Figure 18 illustrates the source discovery through the proposed HF TR method (section 2.4) In order to test NMLA performances. NMLA second order yields a ϕ(x 0 ) = 0.4602s. traveltime and we use this together with the different ray directions of table 2 to shoot the rays backwards. Red is NMLA second order, green LSS and yellow PWD. A red circle centered around the actual source position highlight the better performance of NMLA.

Conclusion

The new second order version of NMLA method is a robust stable local HF wave analysis tool. Unlike the Plane Wave Destructor, the Local Slant Stack methods (Radon) [10] or other signal procession methods like MUSIC [START_REF] Porat | Analysis of the asymptotic relative efficiency of music algorithm[END_REF], it is based on the true wave HF models and robust even to correlated noise. The companion paper [6] also tries to mathematically establish how close in the near field this method is going to work.

It is completely automatic and has no other tuning parameters than α. A mathematical study of the error depending on α and the discretization parameter M is done in [6]. It also proposes an optimal law for α depending on the front cuvature. The source deconvolution is only needed for the traveltime estimation.

We are confident that NMLA can be extended to 3D and other linear models including elastic waves. • GO is simple and cheap when the underlying wave field is not too complicated (no diffraction, reflections, caustics ...). In these more complicated regime, GO has to manage several complex interacting ansatz. As simple propagating waves can be treated locally by GO, coupling GO and full wave fields across physical interfaces or in a multiscale framework make sense [START_REF] Medgyesi-Mitschang | Hybrid methods for analysis of complex scatterers[END_REF] [4] A reliable process to invert the GO ansatz is needed. [START_REF] Nicolay | Gaussian beam decomposition of high frequency wave fields[END_REF], Curvelets [3], diplets, beamlets ... )

• The Direction of Arrival problem (DOA) is a classic signal processing problem [START_REF] Porat | Analysis of the asymptotic relative efficiency of music algorithm[END_REF] where one tries to find the direction of the source of a signal impinging on a, generally linear, array of antennae

• Source point reduced models for far field representations [11] It must be stressed that, currently, none of these methods use traveltime nor the amplitude quantities but only the gradient of the traveltime. This seriously hinders their mathematical and numerical development. data on a circular array around an "observation" point. Regarding real data, circular arrays do not seem to be a common acquisition device either. Interpolation is a solution but it needs data on a volume covering the acquisition circle which typically is several wavelengths large.

Considering the most common acquisition format is a linear array, like for instance seismograms recorded on an ideal flat surface in geophysics, we propose here to use one way wave extrapolation to reconstruct the circular array from a seismogram from this type of data .

We assume for simplicity that :

• We are in 2D and x = (x 1 , x 2 ).

• The data û0 (x 1 ; ω) is given on the surface x 2 = 1.

• Our observation point is on the data surface x 0 = (x o 1 , 1).

• The medium speed is constant c(x) = v in a layer containing the observation circle of radius r (see [START_REF] Giladi | A hybrid numerical asymptotic method for scattering problems[END_REF]).

We recall here the one way extrapolation procedure. More details, references and also the more complicated case of one way extrapolation in heterogeneous medium, can be found in [START_REF] Collino | Peut-on obtenir des amplitudes correctes avec les quations paraxiales ?[END_REF].

First, we Fourier transform the data and Equation ( 12) with respect to (x 1 , x 2 ) and denote the Fourier space variable (k 1 , k 2 ). Fourier transform rules that the partial derivatives ∂ x j corresponds in Fourier space to a multiplication by, respectively, -i k j , that is

(k 2 1 + k 2 2 - ω 2 v 2 ) û(k 1 , k 2 ; ω) = 0 û0 (k 1 ; ω) given on x 2 = 1. ( 39
)
We keep the û notation and drop the ω dependence for simplicity. If we restrict ourselves to real "up-going" k 2 > 0 Fourier mode, this is the one way assumption, and after factorization, equation (39) can be simplified into

(k 2 - ω 2 c 2 -k 2 1 ) û(k 1 , k 2 ) = 0. (40)
The significant real Fourier modes of the solutions have to satisfy the (so called) "dispersion" relation characterizing the paraxial operator :

k 2 = ω 2 c 2 -k 2 1 . (41)
The solution can be computed as a function of (k 1 , x 2 ). Indeed after an inverse Fourier transform in x 2 , we get the ODE (we abusively keep the û notation again)

-i∂ x 2 û(k 1 , x 2 ) = k 2 û(k 1 , x 2 ) û(k 1 , 1) = û0 (k 1 ), ( 42
)
where k 2 is given by (41), with trivial solution The exact solution of (42) at depth x 2 is û(k 1 , x 2 ) = û0 (ω, k 1 ) e i k 2 x 2 .

The solution at any point (x 1 , x 2 ) is recovered by inverse Fourier transform from k 1 to x 1 . we are able to select for instance only the first arrival. In this case we only have one ray and can perform the curvature correction on section 4. The oscillations of S α produced by the truncation of the harmonic series may be a problem to identify the maxima of the B's. To avoid these small spurious oscillations, or to avoid them being mixed with noise, it is possible to convolute B with a smooth Gaussian weight function for instance. This is simply done using a additional filter in Fourier space before the last Fourier transform in [START_REF] Porat | Analysis of the asymptotic relative efficiency of music algorithm[END_REF] :

βℓ ← βℓ ŵℓ { ŵℓ } = F({exp a 2 θ 2 m }) (45)
where a can be used to tune the Gaussian filter. Figure [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF] shows comparison for one source point in homogeneous medium. See also figure 12.

2 )

 2

Figure 1 .

 1 Figure 1. ℓ → ϕ ℓ (red crosses) and parabolic fitting (blue curve)

Figure 2 .

 2 Figure 2. θ (in degree) → |β(θ)| |Bexact| . Left column first NMLA pass, right column second NMLA pass (curvature correction). Top row: standard filter; Bottom row : Gaussian weight filter.

Figure 3 .

 3 Figure 3. Same as figure 2 but with R -= 1.1R exact .

Figure 4 .

 4 Figure 4. Top :k 2π → |Unoisy(k)| L 2 (D) -|Uexact(k)| L 2 (D) |Uexact(k)| L 2 (D)

Figure 5

 5 Figure 5. k 2π → |Best(k)|-|Bexact| |Bexact| . Top : standard NMLA. Bottom : Gaussian weight version.

Figure A2 .

 A2 Figure A2. Functions x → N (L, x) for L = 1, . . . 10. For each x the lowest N L (x, L) can be read as the lowest envelop of all the curves.

Figure A3 .

 A3 Figure A3. Graph of Function x 2 3 ψ(x) over x ∈[4, 3000]

d

  ds y(s, y 0 ) = ∇ϕ(y(s, y 0 )) := p(s, y 0 ), y(0, y 0 ) = y 0 , d ds p(s, y 0 ) =1 2 ∇n 2 (y(s, y 0 )), p(0, y 0 ) = ∇ϕ(y 0 ) given, d ds ϕ(y(s, y 0 )) = n 2 (y(s, y 0 )), ϕ(y 0 ) given,

=Figure 2 :

 2 Figure 2: The stars show |J ℓ (α)i J ′ ℓ (α)| for integer values of ℓ and α = 30.The outputs of the NMLA procedure[START_REF] Porat | Analysis of the asymptotic relative efficiency of music algorithm[END_REF] is  

  )). The virtual source point x s is such that x 0x s = d d where1 d and d are the sought for curvature of the front and direction of the ray, see figure(3).

Figure 3 :

 3 Figure 3: A constant 1 d curvature wavefront going through the observable geometry with ray direction d at the observation point.

Figure 4 :

 4 Figure 4: Observable geometry data in the direction s on the circle is written as d d + rs

Figure 5 :

 5 Figure 5: Blues stars are the (33) values and the red curve the fitted parabola.

Figure 6 :

 6 Figure 6: traveltime for different windows

Figure 7 :

 7 Figure 7: Seismogram and time position picking

Figure 8 :

 8 Figure 8: One source point : |B m | versus θ m . varying α : 10, 20, 40

Figure 9 :

 9 Figure 9: One source point : |B m | versus θ m . Left : Black : NMLA -Blue : NMLA + curvature correction -Red vertical line : exact ray direction. Right : Zoom of the max zone.

Figure 10 :

 10 Figure 10: One source point (0.1, 0.2) : Angle error (solid line) and Traveltime error (dashed line) versus the position of the observation point x 0 on the "surface" y = 1.

Figure 11 :Figure 12 :

 1112 Figure 11: seismogram and time/position picking

Figure 13 :Figure 14 :

 1314 Figure 13: |B m | versus θ m . White noise (20%-40%)

Figure 15 :

 15 Figure 15: Heterogeneous medium : velocity distribution

Figure 16 :

 16 Figure 16: Synthetic data numerical simulation

Figure 18 :

 18 Figure 18: Ray backward propagation with initial conditions specified in table 2.

Figure 19 :

 19 Figure 19: Solid line incorporates Gaussian window. One source point.

Table 1 :

 1 HF ray direction analysis and source point discovery

		Ray angle error Error (distance to source point)
	NMLA basic	0.5988 o
	NMLA 2 nd order	0.0184 o	0.0011
	LSS	failed
	PWD	0.1917 o	0.0071

•

  In Geophysics "Cleaning" data by filtering out low frequencies either to use faster HF models of simply because only high frequencies are really Figure 17: Synthetic seismogram + point/time picking

	Method	Ray angle Error (distance to source point)
	NMLA basic	90.2970 o
	NMLA 2 nd order	90.5503 o	0.0069
	LSS	91.002 o	0.0878
	PWD	90.7418 o	0.0201

Table 2 :

 2 HF ray direction analysis and source point discovery relevant (migration) or nicer for inversion. (use of Gaussian Beams

A common acquisition geometry, seismogram in geophysics for instance.

We established, see (6), that the norm of the filter is bounded by the quantity

Source point discovery through High Frequency Asymptotic Time Reversal

Jean-David BENAMOU -Francis COLLINO -Simon MARMORAT January 18, 2011

Abstract

Given local frequency domain wave data, the Numerical MicroLocal Analysis (NMLA) method [7] in its recent improved version [6] gives a pointwise estimate of the number of rays, their slowness vectors and corresponding wavefront curvatures. With time domain wave data and assuming the source wavelet is given, the method also estimates the traveltime. The paper provides a pedestrian presentation of the NMLA algorithm and a numerical application which can be interpreted as a High Frequency asymptotic version of the classical Time Reversal method [9].

keyword : Linear Waves, Geometrical Optics, Time Reversal

Introduction

Source point or scatterer location identification is a classical inverse problem in computational wave propagation. It is of great practical interest in, for example, non destructive testing of alloys defect, kidney stones location discovery [START_REF] Cassereau | Time reversal techniques in ultrasonic nondestructive testing of scattering media[END_REF], underwater communications [START_REF] Hodgkiss | Phase conjugation in the ocean: Experimental demonstration of an acoustic time-reversal mirror[END_REF]. It is a basic step in inverse synthetic aperture radar [10]. See also [START_REF] Landa | Discovery if point sources in the hemlhotz equation posed in unknown domains with obstacle[END_REF] for the same problem of source discovery.

The popular time reversal (TR) method (see [9] for a mathematical study and references) solves this inverse problem. It is rapidly summarized in the following steps :

(i) record a wave signal on some array of "transducers" for "some" time Appendix A : Seimogram data/ FT / One-way extrapolation of linear array data

The circular geometry of the observable's support plays an important part in the properties of the K α operator [START_REF] Langdon | High frequency scattering by convex curvilinear polygons[END_REF]. This is however a constraint for synthetic data as few numerical discretizations of wave fields directly provide This procedure generates the data on a circular array for up-going rays in constant medium at the price of two 1D Fourier transforms. It could also be used to filter out down-going rays. As we restrict to observation point near the x 2 = 1 depth of the linear array, it is possible to use a x 1 windows which size depends on the finite speed v and the resulting dependence cone.

Appendix B : When the wavelet W ω s is unknown

In that case, the source deconvolution ( 18) is not possible so we just remove the half primitive factor

This implies that the NMLA algorithm expected output is multiplied by the frequency domain wavelet

The curvature correction method of section 3.8 is unchanged since the amplitudes are removed in (35).

The traveltime estimation method (38) does not work since the amplitude B m ⋆ contains now the wavelet with an unknown ω dependant phase.

Appendix C :Time and direction windows

As shown in Figure 17 we assume that, after picking the observation point x 0 = (x 0 1 , 1) on the surface data (the vertical red line), we select an observation time in the seismogram (the horizontal red line). This can be done by the user or an ad hoc automatic signal processing procedure.

The advantage of time picking, is that rays crossing at the observation point with different traveltime may be separated using a time window. In Figure 17 we see that a second later wavefront passes the observation point. If we just Fourier transform in time and then perform NMLA, we expect β α to correspond a two ray ansatz. Conversely, using a smooth time window,


A landscape

come   L'archive ouverte pluridisciplinaire

A landscape approach to sustainable forest management and habitat quality modelling

Luque, S. (1) (1) CEMAGREF, Institute for Agricultural and Environmental Engineering Research, Grenoble, France. sandra.luque@cemagref.fr ABSTRACT This work provides views and examples on how the holistic approach of landscape ecology can be considered towards the interplay between biodiversity value and the needs of forestry activities. Focusing on biodiversity conservation as a proxy for ecological dimensions of sustainable forest management the work proposes a series of requirements for the conservation of habitat networks and ecological processes to be met by forest landscape managers. This is particular challenging nowadays because of the increasing pressures to intensify wood production and timber exploitation. Nevertheless, there are also demands for improving the actions in favour of safeguarding biodiversity, and in a more general way improving forest ecosystems functioning. This twofold requirement of an optimized forest production and an environmental quality improvement represents a true challenge for the years to come. A fast reorganization of the system is needed in order to find the right balance between management within a forest ecosystem services approach at different scales. In particular, we need to consider a valorisation of wood resources and production, thorough knowledge of their vulnerability within an intensification management scenario. Vis-a-vis of such stakes, which concepts, methodologies and tools can be validated on strong scientific grounds that can be proposed to the actors charged to implement policies and actions on the ground? Landscape ecology, proposes a comprehensive and integrative approach from the plot level to the landscape level, considering adaptive management and an analysis of ecological thresholds. Spatial conservation prioritization methods are applied to a national-scale conservation planning task.

Results from Boreal and temperate forests in Europe demonstrate the importance of forest habitat quality models and the development of a spatial framework to integrate information of protected and private forest land. The present analysis provides a large ecological scale approach, where certain natural disturbances can be considered. Forest management operations, can be targeted for the maintenance of the quality of the protected areas. The work opens questions regarding the needs for a comprehensive adaptive forest management.

Adaptive management seems to be the only options that will contribute to balance multiple objectives under changing environmental conditions to improve forest management in a wide range of territories.