Maya De Buhan
email: maya.de-buhan@parisdescartes.fr

Marie Kray
email: marie.kray@unibas.ch

A new approach to solve the inverse scattering problem for waves: combining the TRAC and the Adaptive Inversion methods

The aim of this paper is to propose a new method to solve the inverse scattering problem. This method works directly in the time-dependent domain, using the wave equation and proceeds in two steps. The first step is the Time-Reversed Absorbing Condition (TRAC) method to reconstruct and regularize the signal and to reduce the computational domain. The second step is the Adaptive Inversion method to solve the inverse problem from the TRAC data, by using basis and mesh adaptation. This strategy allows to recover the position, the shape and the properties of the scatterer in a precise and robust manner.

Introduction

In this paper, we propose a new method to solve the following inverse problem: we aim at reconstructing the location, the shape and the wave propagation speed of an unknown inclusion, surrounded by a medium whose properties are known, from boundary measurements.

Literature about inverse problems

The literature on inverse problems for waves is quite huge and the list that follows is not exhaustive. Let us first mention a large literature on the classically named inverse scattering problem. The goal is here to determine only the position and the shape of the scatterer from measurements of the far field pattern. We refer to the surveys by Colton, Coyle, and Monk [START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF] and Cakoni and Colton [START_REF] Cakoni | Qualitative methods in inverse scattering theory: an introduction[END_REF] on the inverse scattering theory.

Concerning numerical methods, let us now mention iterative methods like the least-square method, which are summarized by Tarantola [START_REF] Tarantola | Inverse Problem Theory And Methods For Model Parameter Estimation[END_REF], or like level-set methods as proposed by Burger and Osher [START_REF] Burger | A survey on level set methods for inverse problems and optimal design[END_REF]. Another theoretical point of view is given by Ammari and Kang [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF] about least-square algorithms for detection of inclusions. In parallel, they present the linear sampling method with the same application. Sampling methods are alternative to reconstruct sources or detect inclusions. The more famous is the linear sampling method introduced by Colton and Kirsch [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF]. It has similarities with the MUltiple SIgnal Classification algorithm (MUSIC) developed by Devaney [START_REF] Devaney | Super-resolution processing of multi-static data using time reversal and MUSIC[END_REF], and discussions between the linear sampling method and the MUSIC algorithm can be found in [START_REF] Cheney | The linear sampling method and the MUSIC algorithm[END_REF][START_REF] Kirsch | The MUSIC-algorithm and the factorization method in inverse scattering theory for inhomogeneous media[END_REF]. As sampling methods, we can cite the probe method [START_REF] Ikehata | Reconstruction of an obstacle from the scattering amplitude at a fixed frequency[END_REF], the singularsource method [START_REF] Ben Hassen | The singular sources method for 3D inverse acoustic obstacle scattering problems[END_REF] and the factorization method [START_REF] Lechleiter | The factorization method is independent of transmission eigenvalues[END_REF]. A survey about sampling methods is also proposed by Potthast [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF]. By the same author, we can mention another method of reconstruction called the point-source method described in [START_REF] Potthast | Point Sources and Multipoles in Inverse Scattering Theory[END_REF]. Time reversal techniques have also applications to inverse scattering problems. For example, Prada and Fink [START_REF] Prada | Eigenmodes of the time reversal operator: a solution to selective focusing in multiple-target media[END_REF][START_REF] Mordant | Highly resolved detection and selective focusing in a waveguide using the D.O.R.T. method[END_REF] developed the D.O.R.T. method, based on the decomposition of the time reversal operator to recover the location of the inclusion. A theoretical study about time reversal techniques in non-homogeneous media can be found in [START_REF] Fouque | Wave propagation and time reversal in randomly layered media[END_REF], then application to imaging are developed by Borcea et al. [START_REF] Borcea | Imaging and time reversal in random media[END_REF]. A lot of these methods are working in the time-harmonic domain and use Fourier transform to treat the signal. For applications in the time-dependent domain, we refer to Lines and Chandler-Wilde [START_REF] Lines | A time domain point source method for inverse scattering by rough surfaces[END_REF], who work on the point-source method, to H. Haddar et al. [START_REF] Haddar | An improved time domain linear sampling method for Robin and Neumann obstacles[END_REF] for the linear sampling method, and to Cassier and Hazard [START_REF] Cassier | Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: mathematical justification of the foldy-lax model[END_REF] for time reversal techniques.

Note that some of these methods do not need to know the type of inclusions to be found, soundsoft or sound-hard mostly but in general, they do not try to recover the properties of the scatterer. In our case, we are interested in reconstructing the shape but also the wave propagation speed inside the inclusion. We also refer to the work by Serranho [START_REF] Serranho | A hybrid method for inverse scattering for shape and impedance[END_REF] for shape and impedance reconstruction and by Ben Hassem et al. [START_REF] Ben Hassen | Application of the linear sampling method to identify cracks with impedance boundary conditions[END_REF] and Cakoni et al. [START_REF] Cakoni | The determination of boundary coefficients from far field measurements[END_REF] for impedance reconstruction on boundary conditions. As for the inverse scattering problem, the first issue of concern is that of uniqueness, i.e. does the far field pattern (respectively the Dirichlet-to-Neumann map in a bounded domain) uniquely determine the properties of the inclusion? This question was solved by Sylvester and Uhlmann in their fundamental work [START_REF] Sylvester | A global uniqueness theorem for an inverse boundary value problem[END_REF] in the time-harmonic case. Concerning stability results in the time-dependent domain, we can mention the results by Imanuvilov and Yamamoto [START_REF] Imanuvilov | Global lipschitz stability in an inverse hyperbolic problem by interior observations[END_REF] for a single internal measurement and by Klibanov and Yamamoto [START_REF] Klibanov | Lipschitz stability of an inverse problem for an acoustic equation[END_REF] for a boundary measurements. These results are based on Carleman estimates. A lot of results about coefficient inverse problems using Carleman estimates can be found also in Klibanov and Timonov [START_REF] Klibanov | Carleman Estimates for Coefficient Inverse Problems and Numerical Applications[END_REF] and references therein. From a numerical point of view, all existing methods for determining the wave propagation speed without linearizing the problem are based on nonlinear optimization methods. We mention in particular the numerous works by Beilina and Klibanov [START_REF] Beilina | Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems[END_REF]. These methods give good reconstruction results but suffer from extremely time consuming computations. In this article, we will show that in our case, we are able to drastically decrease this computational time and make the reconstruction conceivable in practice.

New approach

Our strategy combines two methods recently developed by the authors:

1. the Time-Reversed Absorbing Condition method (TRAC) first introduced in [START_REF] Assous | Time Reversed Absorbing Condition: Application to inverse problem[END_REF]:

It combines time reversal techniques and absorbing boundary conditions. Since the seminal paper by Fink et al. [START_REF] Fink | Imaging through inhomogeneous media using time reversal mirrors[END_REF], time reversal is a subject of very active research. The main idea is to take advantage of the reversibility of wave propagation phenomena such as it occurs in acoustics, elasticity or electromagnetism in a non dissipative unknown medium to backpropagate signals to the sources that emitted them. Number of industrial applications have already been developed: touchscreen, medical imaging, non destructive testing and underwater communications. However, the resolution of these applications is restricted by the diffraction limit. The motivation of the TRAC method comes from the experiment of Larmat et al. [START_REF] Larmat | Time-reversal imaging of seismic sources and application to the great Sumatra earthquake[END_REF].

Our idea is to overcome the diffraction limit by removing a subdomain and imposing a special boundary condition on the artificial border. Our previous study shows us that this condition must be a time reversed absorbing boundary condition. Absorbing boundary conditions (ABCs) or non-reflecting boundary conditions (NRBCs) have been developed since the seminal paper by Engquist and Majda [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] for planar boundaries, afterwards for spherical boundaries by Bayliss and Turkel [START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF], for the wave-like equation. Further works are also proposed by Grote and Keller in [START_REF] Grote | Exact nonreflecting boundary conditions for the time dependent wave equation[END_REF] for wave-like equation, in [START_REF] Grote | Nonreflecting boundary conditions for Maxwell's equations[END_REF] for Maxwell's equations and in [START_REF] Grote | Nonreflecting boundary conditions for elastodynamic scattering[END_REF] for elastic waves.

The TRAC method has at least two applications: the reduction of the size of the computational domain and the determination, from boundary measurements, of the location and volume of an unknown inclusion. The method does not rely on any a priori knowledge of the physical properties of the inclusion. The second application has already been well developed in [START_REF] Assous | Time Reversed Absorbing Condition: Application to inverse problem[END_REF][START_REF] Assous | Time Reversed Absorbing Condition in the Partial Aperture Case[END_REF].

In this paper we focus on the first application. We use the TRAC method to reconstruct the signal in a truncated domain that encloses the inclusion. This enables to reduce the size of the computational domain where we solve the inverse problem, now from virtual internal measurements. Moreover the TRAC method is fairly insensitive with respect to the noise on the recorded data and then allows us to work with regularized data without using any filtering method.

2. the Adaptive Inversion (AI) method initially proposed for the viscoelasticity equation in de Buhan and Osses [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF]: This method solves the inverse problem by minimizing a non-quadratic functional that classically depends on the gap between the experimental measurements and the simulated data. For the minimization, we use a Newton-type algorithm. The originality of this method comes from the parametrization of the problem. Instead of looking for the value of the unknown parameter at each node of the mesh, it projects the parameter into a basis composed by eigenvectors of the Laplacian operator. Then, the AI method uses an iterative process to adapt the mesh (see P. Frey [START_REF] Frey | Le maillage facile[END_REF]) and the basis of eigenfunctions to the previous approximation. This process allows to concentrate the vertices of the mesh and the variations of the eigenfunctions in the region of the inclusion interface. In this way, the interface can be more precisely captured and the reconstruction is improved without increasing the computational time. This method has been successfully applied in [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF] to recover a parameter in the integro-differential hyperbolic system of the viscoelasticity from internal measurements. Previously, thanks to a Carleman estimate, we showed that this reconstruction was possible from only one measurement in de Buhan and Osses [START_REF] De Buhan | A stability result in parameter estimation of the 3D viscoelasticity system[END_REF].

The novelty of our work is threefold. Firstly, we present a new study on the regularizing power of the TRAC method. Secondly, we adapt the Adaptive Inversion method to the case of the wave equation and we propose a new anisotropic version of the iterative process. Finally, we present a large variety of numerical examples to illustrate the efficiency of the combination of both methods. In particular, our strategy allows (a) to reduce the computational cost, (b) to stabilize the inverse problem and (c) to improve the precision of the results.

Specifications and assumptions

Let us now describe in details the framework of our study. We hope this will help the reader to replace without any ambiguity this work in the large literature on inverse problems. The context of our method is the following: we want to recover the position, the shape and the wave propagation speed of a penetrable inclusion. In particular, the surrounded medium is known and for the sake of simplicity it is supposed homogeneous. The wave propagation speed inside the inclusion can be non constant but is superior to the one of the surrounded media. We refer to Baudouin et al [START_REF] Baudouin | A global Carleman estimate in a transmission wave equation and application to a one-measurement inverse problem[END_REF] for this necessary condition to guarantee the stability of the reconstruction in the presence of an interface.

Our method works directly in the time-dependent domain. The physical equation is the acoustic wave equation: our study consists on a singlewave imaging using ultrasound waves only, unlike the work of Fink and Tanter [START_REF] Fink | Multiwave imaging and super resolution[END_REF] who propose also multiwave imaging by combining ultrasound waves and shear waves. Initially, the domain is at rest and we highlight the inclusion by a given incident wave with compact support in time. We deliberately use only one incident wave, coming from a single source, like in the point-source method [START_REF] Lines | A time domain point source method for inverse scattering by rough surfaces[END_REF] and unlike other methods, for example the linear sampling method [START_REF] Haddar | An improved time domain linear sampling method for Robin and Neumann obstacles[END_REF] that is using several sources. This restrictive framework obviously makes the problem more difficult but is of interest for various applications. For example, in biomedical applications we aim at reducing the irradiation of the patients and in a military context, it is a way to remain furtive. The total field is then recorded on the source-receivers array (SRA) that encloses the interest area. This SRA is assumed to be continuous and of full aperture. This assumption is reasonable by considering application in medical imaging and will be relaxed in a future study. Note that we do not record any information about the normal derivative of the total field on the SRA. For the sake of simplicity, we work in a two-dimensional domain, but we hope to extend our work to three-dimensional cases. Until now, we only used synthetic data with multiplicative noise, see [START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF]. For this study to be meaningful, we are expecting collaborations with experimental teams that could be interested by the type of results presented here.

The remainder of this paper is divided in three parts. We recall in a first section the principle of the TRAC method. Then we present the Adaptive Inversion method for the wave equation in a second section. Finally, we propose some numerical results obtained thanks to the combination of both methods that illustrate the quality of the reconstruction.

First step: the TRAC method

The TRAC method was first introduced in [START_REF] Assous | Time Reversed Absorbing Condition: Application to inverse problem[END_REF] in the general case with full aperture of the sourcereceivers array (SRA), then in [START_REF] Assous | Time Reversed Absorbing Condition in the Partial Aperture Case[END_REF] in the time-dependent domain with partial aperture of the SRA and discrete receivers. Here we consider the TRAC method in the case of the wave equation with a full aperture SRA. We first recall the principle of the method in this particular case. Then we illustrate the TRAC method and propose a numerical analysis of the error of reconstruction. Finally we present the application to the reduction of the size of the computational domain, first step before solving the inverse problem.

Principle of the TRAC method

The context of the forward experiment is the following: we consider an incident wave u I impinging on an inclusion D characterized by different physical properties from the homogeneous surrounding medium, c(x) = c 0 in R 2 \ D. The total field u T can be decomposed into the incident and scattered fields, so u T = u I + u S . For simplicity, we consider the problem in two space dimensions and we assume that the total field satisfies the linear wave equation:

∂ 2 u T ∂t 2 -∇ • (c 2 (x)∇u T ) = 0 in (0, ∞) × R 2 , (1)
together with homogeneous initial conditions. The scattered field u S has a finite energy at any time.

Let Γ R be a curve that defines a bounded domain Ω and encloses the inclusion D, see Figure 1.

After a time T f , the total field u T is negligible in Ω. The scattered field u S is recorded on Γ R on the time interval [0, T f ]. Let u S R (t, x) := u S (T f -t, x) denote the scattered time reversed field which also satisfies [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF].

Our aim with the TRAC method is to reconstruct the time reversed field u S R from the measurements on Γ R . For this purpose, we derive a boundary value problem whose solution is u S R in Ω. Yet we know neither the physical properties nor the exact location of the inclusion D, but only the physical properties of the surrounding medium. Therefore, we introduce a subdomain B which encloses the inclusion D. Then, we have to determine a boundary condition for u S R on the boundary ∂B such that the solution to this problem will coincide with u S R restricted to Ω \ B. Let us choose B as a ball of radius r. On ∂B, we approximate the radiation condition satisfied by u S with the first-order Bayliss-Turkel (BT 1 ) absorbing boundary condition:

ABC(u S ) := ∂u S ∂t + c 0 ∂u S ∂r + c 0 u S 2r = 0 on (0, T f ) × ∂B. (2)
Next we time reverse this relation using

u S R (t, •) = u S (T f -t, •).
In doing so, we note that ∂/∂r = -∂/∂n on ∂B where n is the outward normal to Ω \ B. Hence, by multiplying the relation above by -1, we obtain:

TRAC(u S R ) := ∂u S R ∂t + c 0 ∂u S R ∂n -c 0 u S R 2r = 0 on (0, T f ) × ∂B. (3)
Hence the time reversed problem for the scattered field, analogous to (1), reads:

               ∂ 2 u S R ∂t 2 -c 0 ∆u S R = 0 in (0, T f ) × Ω \ B, TRAC(u S R ) = 0 on (0, T f ) × ∂B, u S R (t, •) = u S (T f -t, •) on (0, T f ) × Γ R , homogeneous initial conditions. (4)
Note the anti-absorbing term (-c 0 u S R /2r) in the TRAC condition. In [START_REF] Assous | Time reversed absorbing conditions[END_REF] we proved an energy estimate in three space dimensions to state the stability of the time-reversed problem with a timereversed absorbing condition.

For the following applications, we need to reconstruct the total field. It is necessary in order to solve the inverse problem with the proposed method and more logically because we record the total field in practice. Hence the boundary condition on ∂B has to be changed:

TRAC(u S R ) = TRAC(u T R -u I R ) = TRAC(u T R ) -TRAC(u I R ).
Then the time reversed problem for the total field reads:

               ∂ 2 u T R ∂t 2 -c 0 ∆u T R = 0 in (0, T f ) × Ω \ B, TRAC(u T R ) = TRAC(u I R ) on (0, T f ) × ∂B, u T R (t, •) = u T (T f -t, •) on (0, T f ) × Γ R ,
homogeneous initial conditions.

( 5) Note that to solve this problem, we need to know the incident wave in a neighborhood of ∂B. In practice we assume to know the source which emits and the medium without inclusion, then we compute the incident wave numerically.

Ω \ B Γ R := SRA B D
We have assumed so far, for the sake of simplicity, that the surface B is a ball. Since we are finding an approximate location of the inclusion, this is usually sufficient. For an elongated body, a ball can be replaced by an ellipse or spheroidal surface. Absorbing boundary conditions for these cases have been developed in [START_REF] Medvinsky | Local absorbing boundary conditions for elliptical shaped boundaries[END_REF][START_REF] Medvinsky | On surface radiation conditions for an ellipse[END_REF] and for more general surfaces see, for example [START_REF] Antoine | Bayliss-turkel like radiation conditions on surfaces of arbitrary shape[END_REF]. As shown above, a first-order TRAC method simply reverses the sign of the non-differentiated term of the corresponding first-order absorbing boundary condition. Thus, a first-order TRAC for a general subdomain in two space dimensions is given by

TRAC(u S R ) := ∂u S R ∂t + c 0 ∂u S R ∂n -c 0 κ u S R 2 = 0 on (0, T f ) × ∂B, ( 6
)
where κ is the curvature of the boundary ∂B.

Numerical error analysis of the reconstruction with the TRAC method

On Figure 2, we display the results of the TRAC method for a reflecting fish-shaped inclusion. We also proposed other time reversal simulations. Figure 2 is composed of five subfigures, each of them contains four snapshots of the wave propagation from initial time t = 0 on the left column (corresponding to final time t = T f of the forward problem) to final time t = T f on the right column (corresponding to initial time t = 0 of the forward problem).

The first subfigure displays the perfect time reversal. It is simply the forward scattering field generated by the reflecting fish for decreasing times. Note that we know the forward scattering field for this case because we simulate the synthetic data. We use it to illustrate the TRAC method in an academical way. In practice, we are not assumed to know the forward scattering field in the whole domain, only on the recording boundary. The second subfigure illustrates the TRAC method with a subdomain B which correctly satisfies the assumption: B encloses inclusion D. This numerical time reversal experiment was performed with discrete receivers, spaced of a quarter of the central wavelength, and with 20% of noise on the recording data. In this case, we observe that the TRAC method enables us to reconstruct the scattered field in the truncated computational domain Ω \ B. The third subfigure shows the classical time reversal. We indeed perform time reversal, but without removing the artificial subdomain B necessary in the TRAC method. Moreover we do not assum to know neither the location nor the properties of the inclusion, then the idea of classical time reversal is to solve the boundary value problem with c(x) = c 0 in the whole domain Ω. Because of this error on the medium properties, we are not able to reconstruct correctly the scattered field. On the fourth subfigure, we display a case of the TRAC method where the assumption of B enclosing D is not satisfied. The tail of the fish stays indeed in the truncated computational domain Ω \ B. These pictures show the weight of this assumption because we are not able to reconstruct the scattered field anymore. Thanks to this observation, we developed an identification technique of the presence of an inclusion, or of the absence of an inclusion in the truncated computational domain, see [START_REF] Assous | Time Reversed Absorbing Condition: Application to inverse problem[END_REF]. Finally the fifth subfigure is simply an illustration of the choice of the boundary condition imposed on ∂B. In this case, we correctly enclose inclusion D with the subdomain B. If the condition is not TRAC, but for example a homogeneous Neumann condition, the numerical time reversal experimental does not give the correct scattered field.

To quantify the error due to the TRAC method, we propose to compute the relative L 2 -error for nine different tests and for different levels of noise on the recording data. This error is computed in a two-wavelength thick ring of the truncated computational domain Ω\B containing the border ∂B. The results are displayed in Table 1. This table is composed of five lines and eleven columns. The first line gives us the number of the considered test, when the other four lines give the error for each test depending on a level of noise, respectively 0%, 5%, 10% and 20%. The last column gives a mean value of the error depending on the level of noise on the recording data. The results are really satisfactory. We indeed observe that even if the error is already of 5% or 10% without noise, it is stable for increasing level of noise on the recording data. This error is certainly due to the approximation of the TRAC method which adds mesh interpolation errors. On average, when the noise level is 10%, we obtain a scattered field with only 5.66% of error. The TRAC method is fairly insensitive with respect to the noise on the recording data and regularizes the solution at least in a neighborhood of the artificial boundary ∂B.

Reduction of the size of the computational domain

An application of the TRAC method is the reduction of the size of the computational domain. The principle is to reconstruct the field in a part of the computational domain and define a new surface where the receivers seem virtually to be located. For example, with the TRAC method we are able to move virtually the receiver border Γ R to the closer boundary ∂B. Then we solve the inverse problem in the truncated and smaller computational domain B with measurements from ∂B. This is reminiscent of the redatuming method that uses a paraxial approximation of the wave equation to propagate the signal in the domain and move virtually the SRA, see the seminal work of Berryhill [START_REF] Berryhill | Wave-equation datuming[END_REF], and more recently [START_REF] Liu | Image enhancement with wave-equation redatuming: application to GPR data collected at public transportation sites[END_REF][START_REF] Mulder | Rigorous redatuming[END_REF]. In our method, we are able to do the same but without the paraxial approximation. Then the reconstruction must be more accurate because does not require the accuracy of the chosen paraxial approximation.

Moreover the TRAC method regularizes the reconstructed field with respect to the noise on the recorded data, by assuming that the noise is due to measurement errors, see Table 1. Hence we do not need any filter method to smooth our resulting signal before performing the inverse problem. The location of subdomain B can be chosen thanks to discrimination criteria, developed in [START_REF] Assous | Time Reversed Absorbing Condition: Application to inverse problem[END_REF] and [START_REF] Assous | Time Reversed Absorbing Condition in the Partial Aperture Case[END_REF], which tell us when the trial subdomain B encloses entirely the unknown inclusion D. The idea here is not to get a subdomain that is very close to the inclusion. We need simply a subdomain that encloses the inclusion and that enables us to reconstruct the signal.

The first step of our technique is to perform the TRAC method to reconstruct the scattered field, or the total field similarly. Then we use the regularized data on the neighborhood of ∂B, denoted by ω, to solve an inverse problem in a smaller computational domain, see Figure 3 and next section. Notice that we do not know the normal derivative on the SRA but only measurements of the field. To perform the TRAC method, we need also to know the incident wave, i.e. the source, but it is not needed then to solve the inverse problem.

Second step : the Adaptive Inversion method

The Adaptive Inversion method has been initially proposed for the viscoelasticity equation in [START_REF] De Buhan | Logarithmic stability in determination of a 3D viscoelastic coefficient and a numerical example[END_REF]. The originality of this method comes from the parametrization of the problem. Instead of looking for the value of the unknown coefficient at each node of the mesh, it projects the coefficient into a basis composed by the eigenvectors of the Laplacian operator. Then, it uses an iterative process to adapt the mesh and the basis and improve the reconstruction. This method takes as an entry the value of the observed solution in an annulus enclosing the unknown domain. And this is exactly what we obtained in the previous section using the TRAC method.

Inverse problem

Let B ∈ R 2 be the domain obtained at the first step by the TRAC method and let B another open domain strictly containing B. We introduce ω = B \ B, as shown on Figure 3. We denote by u T R the solution reconstructed by the TRAC method in ω. The latter is the observation taken as entry for the inverse problem.

We consider the wave equation in B with Dirichlet boundary condition:

                 ∂ 2 v ∂t 2 -∇ • (c 2 (x)∇v) = 0 in (T 1 , T 2 ) × B, v(T 1 , •) = 0 in B, ∂v ∂t (T 1 , •) = 0 in B, v = u T R on (T 1 , T 2 ) × ∂ B, (7)
with c 2 (x) = c 2 0 + p(x)χ B (x), for all x ∈ B, where χ B is the characteristic function of B. The time T 1 (respectively T 2 ) is the first time (resp. last time) when the solution u T R in ω in non zero. We have 0 < T 1 < T 2 < T f . The inverse problem that we consider can be written as follows: Find p such that the solution v of system (7) associated to p coincide with u T R in (T 1 , T 2 ) × ω. To solve the inverse problem, we look for the minimizer of the non-quadratic functional This minimization of the gap between the observations and the solution v of system (7) associated to p is a quite classical way for solving an inverse problem. And because the inverse problem is ill-posed in the L 2 ((T 1 , T 2 ) × ω)-norm considered for J, we add a regularizing term to stabilize the method. The choice of the regularization parameter α is crucial. If α is too small, the reconstruction will be unstable but if α is too large, the method will converge to another function than the p we are looking for. A compromise between stability and precision have to be found. In practice, we choose α such that both terms in J have the same order of magnitude at the initialization of the algorithm. We solve the minimization problem by a Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. To do that, we have to compute:

J(p) = 1 2 T 2 T 1 ω |v(t, x) -u T R (t, x)| 2 dxdt + α 2 B |∇p(x)| 2 dx. (8)
∇J(p, δp) = T 2 T 1 ω (v(t, x) -u T R (t, x))δv(x, t)dxdt + α B ∇p(x) • ∇δp(x)dx,
where δv is the solution of the following linearized problem:

∂ 2 δv ∂t 2 -∇ • (c 2 (x)∇δv) = ∇ • (δp(x)χ B ∇v) in (T 1 , T 2 ) × B,
with zero initial and boundary conditions. Now, if we introduce v * the solution of the adjoint problem:

                 ∂ 2 v * ∂t 2 -∇ • (c 2 (x)∇v * ) = v -u T R 0 in (T 1 , T 2 ) × ω, in (T 1 , T 2 ) × B, v * (T 2 , •) = 0 in B, ∂v * ∂t (T 2 , •) = 0 in B, v * = 0 on (T 1 , T 2 ) × ∂ B,
we can simplify the expression of ∇J as follows:

∇J(p, δp) = - B δp(x) T 2 T 1 ∇v(x, t) • ∇v * (x, t)dt dx + α B ∇p(x) • ∇δp(x)dx.
In this form, for each p, we can compute the term

T 2
T 1 ∇v(x, t) • ∇v * (x, t)dt once and then obtain the gradient for each direction δp by a simple computation of integral.

Parametrization

The unknown coefficient p is not a constant but a function of x in B. In order to recover it, we have first to project the function into a discrete space. A natural choice could be to use a P 1 Lagrange finite elements basis leading thus to recover the values of p at each node of the computational mesh. Unfortunately, this leads to very large minimization spaces if we want to represent correctly the variations of p.

In our method, we propose to use another L 2 basis and to look for the coefficient p in the space spanned by the K first eigenfunctions of the Laplacian operator, that is

p(x) = K i=1 p i ϕ i (x),
where

-∆ϕ i = σ i ϕ i in B, ϕ i = 0 on ∂B. (9)
The unknown are the p i for i = 1 to K. This allows to decouple the size of the mesh and the dimension of the minimization space (size K).

On Figure 4, we plot some eigenfunctions of the Laplacian operator defined in [START_REF] Beilina | Adaptivity with relaxation for ill-posed problems and global convergence for a coefficient inverse problems[END_REF]. The figure is composed of two lines and three columns. The top line shows a two dimensional view of the chosen eigenfunctions, and the bottom line shows a three dimensional view of the same eigenfunctions. From left to right, we display eigenfunctions n o 0, 1 and 49. On Figure 5, we compare on a simple example both parametrizations: with finite elements basis functions and with Laplacian eigenfunctions. The Figure is composed of three lines and three columns. The top line displays the exact coefficient. The middle line considers the case of the finite elements basis whereas the bottom line shows the case of the Laplacian eigenfunctions. Both cases are presented for K = 50 basis functions. Then we proposed, for both cases, the mesh on the left, a two dimensional view in the middle and a three dimensional view on the right for the coefficient decomposition. We observe that with the same number of functions, here K = 50, the Laplacian basis allows a better representation of the coefficient p.

Mesh and basis adaptation

Then, we propose an adaptive method to improve the accuracy of the reconstruction. After computing a first solution p (0) in the initial mesh using the method described in the previous sections, we use this solution to adapt the mesh and to find another basis that better represents the coefficient. The mesh adaptation used here is classical. It is based on the Hessian of the solution p (0) and it enables to concentrate the nodes of the mesh where the solution varies to decrease the approximation error without increasing the computational time. See Figure 7(a) for an example. We could also envisage an anisotropic mesh adaptation to better capture the interface of D, see Frey and Georges [START_REF] Frey | Le maillage facile[END_REF].

The basis adaptation considered is the following: we now look for the unknown coefficient p in the space spanned by the K first eigenfunctions of an elliptic operator, that is where the matrix function A is chosen depending on p (0) . More precisely, we propose the following form for A A = 1 |∇p (0) | q Id, with q ∈ N. [START_REF] Ben Hassen | The singular sources method for 3D inverse acoustic obstacle scattering problems[END_REF] Remark : To avoid problems when |∇p 0 | becomes close to zero, we replace it numerically by

p(x) = K i=1 p i ψ i (x), with -∇ • (A(x)∇ψ i ) = σ i ψ i in B, ψ i = 0 on ∂B, (10)
|∇p 0 | ←-max(|∇p 0 |, ε), (12)
where ε > 0 is a small non null parameter.

On Figure 6, in the two first lines, we compare the basis obtained with A for q = 1 and q = 2. In a line, we plot successively the eigenfunctions n o 0, 1 and 49 of the operator defined in [START_REF] Ben Hassen | Application of the linear sampling method to identify cracks with impedance boundary conditions[END_REF] and the projection of the exact coefficient p proposed in Figure 5(a) in the basis for K = 50. We also compute the relative L 2 -norm in each case. We observe that these choices of matrices A allow to concentrate the variations of the eigenfunctions in the regions where coefficient p (0) varies, that is near the interface of the inclusion D. When the power q of the gradient increases, the localization is improved. And we observe that in both cases the relative L 2 -norm is better than the one for the parametrizations of Figure 5.

In the last line of Figure 6, we present another choice for A which is anisotropic. The idea is to take into account the orientation of the interface of D and to accord a preference to variations of the basis functions in the direction of the gradient of p (0) . To do that, as presented in Figure 8, we take a new orthonormed system (x, X 1 , X 2 ) whose first axis is locally oriented by the gradient of p (0) and we define the transformation matrix P as follows:

P = 1 |∇p (0) | 2     ∂p (0) ∂x 1 - ∂p (0) ∂x 2 ∂p (0) ∂x 2 ∂p (0) ∂x 1     .
In the new system, we choose to give more weigh to the direction of the gradient, by setting:

A = 1 |∇p (0) | 2 P CP -1 , with C =   ε 2 |∇p (0) | 2 0 0 1   .
This final choice for A generates a really better basis of representation for p as proved by the corresponding relative L 2 -norm computed in Figure 6.

Finally, the Adaptive Inversion method proposes a numerical process with four steps: at each step we use the solution obtained at the previous step to adapt the mesh and to built the matrix A.

In fact, we use successively the three choices of A, with q = 1, then with q = 2 and finally with q = 2 and anisotropy. On Figure 7, we show the result of this process in the approximation of the exact coefficient p proposed in Figure 5(a). In the first line are the corresponding meshes which are adapted from the previous solution. In the second and third lines, we present the projection of p in the basis. The last line gives the relative L 2 -error between the projected and the exact coefficients.

We are aware that, until now, the method is essentially heuristic and will need a complete mathematical study to understand what is going on and to justify the strategy. Nevertheless, the

(a) A = |∇p (0) | -1 Id, Relative L 2 -error = 1.34% (b) A = |∇p (0) | -2 Id, Relative L 2 -error = 1.05% (c) A = |∇p (0) | -2 P CP -1 , Relative L 2 -error = 0.75%

2.90%

1.34% 0.64% 0.15% results in the inverse problem are very encouraging as shown in Section 4. Let us mention that for the same problem of coefficient reconstruction, Beilina and Klibanov [START_REF] Beilina | Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems[END_REF] propose a numerical algorithm and prove its global convergence. They also prove that mesh adaptation, which possesses the relaxation property, ensures the accuracy improvement, see Beilina et al. [START_REF] Beilina | Adaptivity with relaxation for ill-posed problems and global convergence for a coefficient inverse problems[END_REF].

(d) Relative L 2 -error

Combination of the TRAC and the Adaptive Inversion methods

Numerical data

To reconstruct the shape and the properties of the inclusion, we first generate the synthetic data, then proceed the two steps. To avoid the inverse crime, we use neither the same meshes nor the same numerical schemes to simulate the synthetic data and to solve the inverse problem.

Step 0: synthetic data In this paper, we work with synthetic data. The context of the forward problem is displayed on Figure 9: a punctual source highlights the inclusion D, then the resulting total field is recorded on the SRA.

To create our synthetic data, we use a punctual source in space, with compact support in time. It is chosen to be a Ricker function, see [START_REF] Lines | A time domain point source method for inverse scattering by rough surfaces[END_REF]:

f (t) =    1 -2π 2 (ν 0 t -1) 2 e -π 2 (ν 0 t-1) 2 , if t ∈ (0, T s ) and x = x s , 0, otherwise, (13)
where ν 0 is the central frequency of the signal, T s is the emitting time of the source, equal to two periods, and x s is the location of the source outside the domain delimited by the SRA. The efficient spectrum of this signal, i.e frequencies which contribute with more than 20% of the signal, is [2ν 0 /7, 2ν 0 ], see Figure 10. Now we define two characteristic wavelengths: the minimal wavelength λ min which corresponds to the maximal frequency of the efficient spectrum, and the central wavelength λ which corresponds to the central frequency of the Ricker function. All lengths will be given with respect to both wavelengths as a reference. The computational domain is a ball of radius 10λ, the punctual source is located on the border of this domain and the SRA is a concentric ball of radius 4λ. Inside the domain delimited by the SRA, we have an inclusion of a characteristic size varying between λ/2 and 2λ. The propagation speed c in the medium reads:

c 2 (x) = c 2 0 = 1 in R 2 \ D, c 2 D = 3 in D. (14)
To discretize the forward problem, we use a P 1 -finite elements method in space, with a mesh conformed to the shape of the inclusion, and an explicit leap-frop scheme in time. The characteristic size of a cell of the mesh is h = λ min /10 and the CFL is taken equal to √ 2/8. All simulations of the forward problem are executed with the software FreeFem++ [START_REF] Hecht | FreeFem++. Numerical Mathematics and Scientific Computation[END_REF].

On the recorded data, we add a Gaussian noise:

u T (t j , x i ) := (1 + Coeff × randn(t j , x i )) × u T (t j , x i ) , (15)
where x i denotes the position of receiver i, t j denotes the j th step of time, randn satisfies a centered reduced normal law and Coeff is the level of noise. Note that the model of noise, that we chose, is a multiplicative noise, not an additive one. Our choice is motivated by [START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF]. It allows to perturbate the recorded signal in order to avoid the inverse crime and to model the recording measurements from the receivers. It is not an ambient noise.

Step 1: TRAC method

To perform the TRAC method, we remove a ball of radius 1.4λ which is concentric to the SRA in this case (not necessary, see [START_REF] Assous | Time Reversed Absorbing Condition in the Partial Aperture Case[END_REF]). We work with a new mesh, independent on the previous one for the forward problem. The radius of the new computational domain is 5λ, containing the SRA of radius 4λ. Just note that in this study the SRA is continuous and with full aperture. In Ω \ B, the propagation speed is c(x) = c 0 .

To discretize the TRAC problem, we use a P 1 -finite elements method in space, with a mesh conformed to subdomain B and an explicit leap-frop scheme in time. The characteristic size of a cell of the mesh is still h := λ min /10 = 0.01 and the CFL = √ 2/8. All simulations of the TRAC method are also executed with the software FreeFem++ [START_REF] Hecht | FreeFem++. Numerical Mathematics and Scientific Computation[END_REF]. Once the TRAC method is performed, the first user gives the reconstructed signal and the mesh on Figure 11(b) to the second user. The former does not give any indication about the shape and the properties of the inclusion, that he used to create the synthetic data. In particular, the second user does not need to know the incident wave to solve the inverse problem.

Step 2: Adaptive Inversion method In order to solve the inverse problem, we discretize the equations (7) in space using P 1 Lagrange Finite Elements on the mesh of B shown on Figure 11(c). On Figure 11, we illustrate the reduction of the domain obtained thanks to the TRAC method. From left to right, we present the entire mesh 11(a) including the source of the incident wave, the SRA and the inclusion D, the mesh 11(b) on which the TRAC method is performed, the mesh 11(c) used to solve [START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF] and the mesh 11(d) of the domain B where we recover p. Without the TRAC method, at each iteration of the minimization process of the inverse problem, we would have to solve the forward problem in the whole domain, which is of radius 10λ. We then reduce the degree of liberty of our computational domain from 138678 nodes to 2448 nodes.

In time, we discretize the equations ( 7) using a θ-scheme with θ = 0.5 (implicit centered scheme). Because the chosen scheme is implicit, we do not have to satisfy any stability condition and can reduce the number of time steps by taking only one step over three. Moreover, we also reduce the total time T 2 -T 1 < T f of simulation because we reduced the size of the domain. In fact, we only have to propagate the incident wave from the virtual SRA on ∂ B. Finally, we only have 250 time steps and not 2263 as for the forward problem. At the end, we have divided the computational time by at least 500. This is the reason why we can afford the use of an adaptive method with several steps (four here). Moreover, we always initialize a step with the solution obtained at the previous step, such that the first step might be costly but the following steps converge rapidly.

To solve the inverse problem, we use at each step a basis of 100 functions. The parameter ε in [START_REF] Berryhill | Wave-equation datuming[END_REF] is taken equal to 10 -5 . To solve the eigenvalue problem [START_REF] Ben Hassen | Application of the linear sampling method to identify cracks with impedance boundary conditions[END_REF] on the mesh 11(d), we use the toolbox included in FreeFem++ with the standard parameters. In the same way, we use the BFGS algorithm and the mesh adaptation already implemented in the software FreeFem++. The convergence error of the BFGS algorithm is chosen as 5 • 10 -4 . For the mesh adaptation, we pay a special attention to have the more refined mesh but without increasing its number of vertices.

Numerical results

In this section, we present our numerical results obtained by combining the TRAC and the AI methods. We performed six tests with different shapes of the inclusion D. For each test however, we take the propagation speed described in [START_REF] Burger | A survey on level set methods for inverse problems and optimal design[END_REF].

As a first illustration, we propose the case where the inclusion is a full disk of radius λ. On Figure 12, we show the exact unknown coefficient and the recovered coefficient in the whole domain Ω. This is to show that at this scale, both results are very similar. On Figure 13, we display the results of the successive iterations of the AI methods doing a zoom on B, where we actually solved the inverse problem. This test is performed with the exact data, i.e. by assuming we know the total field in the whole domain. We did not execute the TRAC method as a first step here. In practice, we interpolate the forward total field in the correct mesh to solve then the inverse problem by the AI method. On the top line of Figure 13, we show the two-dimensional view and on the bottom the corresponding three-dimensional view. Figure 13(a) is the result of the 0 th step with the Laplace operator eigenfunctions basis [START_REF] Beilina | Adaptivity with relaxation for ill-posed problems and global convergence for a coefficient inverse problems[END_REF]. Then Figure 13(b) illustrates the first step, using a matrix A with exponent q = 1, in [START_REF] Ben Hassen | The singular sources method for 3D inverse acoustic obstacle scattering problems[END_REF]. In the same spirit, Figures 13(c) and 13(d) are the second and third steps respectively with q = 2 for the matrix A and then with q = 2 and anisotropy for the matrix A. Moreover, at each step, we choose the Tikhonov regularization term α such that both terms in the functional J introduced in (8) have the same value at the initialization of the step. This implies that the value of α is decreasing at each step. The results presented on Figure 13 are really significative. We indeed observe the progress of the method and the improvement of the reconstruction of the shape and the properties of the inclusion. Especially, we are able to find a good approximation of the shape from the zeroth step, but the coefficient is openly oscillating. After each iteration, the solution is smoother and smoother, but without loosing the gradient that marks the shape of the inclusion.

For the test of the λ-radius disk, we also performed the combination of the TRAC method and the AI method. The results are displayed on Figure 15. This time, we use the signal reconstruction and the reduction of the size of the computational domain given by the TRAC method. On this Figure, the top line is used to show the 2D-view, the middle line shows the corresponding 3D-view and on the bottom line, we present the relative L 2 error in the coefficient, computed as follows :

error = B |p -p ex | 2 dx B |p ex | 2 dx
, where p ex is the exact coefficient. Figure 15(a) is actually the exact coefficient p ex , given as a reference. Then Figure 15(b) is the result obtained with the exact data, i.e. the TRAC method is not used as a first step. We assume that we know the forward total field in Ω \ B. We propose this test to emphasize the role of the TRAC method. We want here to show that the results obtained with the TRAC data are as satisfying as with the exact ones. Figures 15(c) and 15(d) illustrate the shape and property reconstruction by the AI method after the TRAC method respectively without noise on the recorded data and with 20% level of noise. Once more, the results are satisfactory. The shape is correctly reconstructed, if we take into account that the AI method allows to find a smooth function associated to our piecewise constant coefficient. On the pictures, we have superimpose the shape of the exact coefficient to the recovered one in order to better compare, but obviously this shape was an unknown of the reconstruction. Note that the scale is the same for all pictures, and thus the reconstructed coefficient has the desired intensity: outside the inclusion c 2 = 1 and inside the inclusion c 2 = 3.

You can also enjoy results for a pentagon on Figure 16, a holed ellipse on Figure 17, two disjoint ellipses on Figure 18, a fish on Figure 19 and a star on Figure 20. All those Figures are composed as Figure 15 described above. The third test corresponds to a case where c is not constant inside the inclusion and we could also have a non piecewise constant function. Once more, all the tests are quite satisfactory. Notice however that we are limited by the size of the inclusions which are here of the same order as the wavelength (Rayleight limitation). Non convex shapes and strong angles (as for the fish and the star) are harder to reconstruct. To get a better resolution in these cases, we can execute the AI method with more than K = 100 eigenfunctions. The more eigenfunctions you take to perform the AI method, the better you can recover small details.

Finally, on Figure 14 we display the adapted mesh for each test, at the end of the global process. One can clearly observe that the mesh is refined close to the boundary of the inclusion and very coarse elsewhere. Moreover those meshes keep the same number of nodes as the original one, when we start the process.

Comments and conclusion

Here we give some numbers to have an idea of the cost of our strategy in the case of the pentagon. To create once the synthetic data with a 138678 nodes mesh and 2263 time steps, the CPU time is of 1200 s. We perform then the TRAC method with CPU time equal to 2503 s and the AI method with a CPU time of 9588 s for the first step, 4094 s for the second step, 1868 s for the third step and 1169 s for the last step. This gives us finally a total CPU time equal to almost 6 hours for the whole process. Remember from §4.1 that thanks to the TRAC method, we have been able to reduce the computational time of the forward problem by more than 500. Without using the TRAC method, one might have to solve the inverse problem in the whole computational domain of the forward problem. And in the AI method, we need to solve several times the forward problem (almost 200 times in the zeroth step). Despite this huge reduction, the CPU time remains high. We could reduce it more by taking coarser meshes, less eigenfunctions in the basis and ask for less convergence in the minimization process, but leading then to less precise results.

Moreover let us precise that the choice of subdomain B to perform the TRAC method does not require too many trials. Thanks to the classical time reversal, we can easily have an approximation of the location of inclusion D. As B does not need to be too close to the inclusion, we may find a satisfying subdomain B in one shot. This can be validated thanks to criteria developed with the TRAC method. Thus this first step will not increase the computational cost.

Finally, we are optimistic in improving our method to overcome the Rayleigh limitation and also to reconstruct non-homogeneous inclusions.

Figure 1 :

 1 Figure 1: Computational domain to perform the TRAC method, delimited by the source-receivers array ΓR and the subdomain B, such that B encloses inclusion D.

  (a) Perfect time reversed scattered field as reference (b) TRAC method with B enclosing inclusion D (c) Classical time reversal: computational domain Ω and B = ∅ (d) TRAC method with B not enclosing D: the tail stays in the computational domain Ω \ B.(e) On the artificial boundary ∂B, homogeneous Neuman conditions

Figure 2 :

 2 Figure 2: Illustration of the TRAC method for a reflecting fish-shaped inclusion: snapshots of the wave propagation from initial time t = 0 (corresponding to final time t = T f of the forward problem) on the left to final time t = T f (corresponding to initial time t = 0 of the forward problem) on the right.

Figure 3 :

 3 Figure 3: New virtual source receivers array (SRA) on ∂B: new data for the inverse problem in the coloured area ω, new computational domain B = ω ∪ B.

Figure 4 :Figure 5 :

 45 Figure 4: Eigenfunctions n o 0, 1 and 49 of the Laplacian operator.

Figure 6 :

 6 Figure 6: Comparison between the different choices for the matrix A: (left) eigenfunctions n o 0, 1 and 49 of the operator (10), (right) projection of p into the corresponding basis.

Figure 7 :Figure 8 :

 78 Figure 7: Illustration of the Adaptive Inversion process.

Figure 9 :

 9 Figure 9: Geometry of the forward problem: location of the source with respect to the SRA and the inclusion.

  Spectrum of the source for a central frequency ν0 = 5.

Figure 10 :

 10 Figure 10: (a) Source and (b) its spectrum: central frequency (dashed line), efficient spectrum (between dotted lines).

Figure 11 :

 11 Figure 11: Illustration of the reduction of the computational domain: (a) mesh for the forward problem, conformed to the inclusion D (here a pentagone), radius of the computational domain R = 10λ, 138678 vertices; (b) mesh for the TRAC method, conformed with the SRA and the subdomain B, radius R = 5λ, 31383 vertices; (c) mesh for the inverse problem in B, conformed with the domain B, radius R = 2λ, 2448 vertices; (d) mesh to recover the coefficient p in B, radius R = 1.4λ, 3573 vertices.

Figure 12 :Figure 13 :

 1213 Figure 12: Result of the Adaptive Inversion method in the entire domain Ω: (a) the exact coefficient, (b) the recovered coefficient at final step.

Figure 14 :

 14 Figure 14: Meshes built with FreeFem++ after adaptation by the AI method.

Figure 15 :Figure 16 :Figure 17 :

 151617 Figure 15: Shape and properties reconstruction of a penetrable disk by using both TRAC and AI methods: (a) Propagation speed profile inside and outside the inclusion. (b) Result obtained with exact data. (c) Result obtained with 0%-noisy TRAC data. (d) Result obtained with 20%-noisy TRAC data.

Figure 18 :

 18 Figure 18: Shape and properties reconstruction of two penetrable ellipses by using both TRAC and AI methods: (a) Propagation speed profile inside and outside the inclusion. (b) Result obtained with exact data. (c) Result obtained with 0%-noisy TRAC data. (d) Result obtained with 20%-noisy TRAC data.

Figure 19 :

 19 Figure 19: Shape and properties reconstruction of a penetrable fish by using both TRAC and AI methods: (a) Propagation speed profile inside and outside the inclusion. (b) Result obtained with exact data. (c) Result obtained with 0%-noisy TRAC data. (d) Result obtained with 20%-noisy TRAC data.

Figure 20 :

 20 Figure 20: Shape and properties reconstruction of a penetrable star by using both TRAC and AI methods: (a) Propagation speed profile inside and outside the inclusion. (b) Result obtained with exact data. (c) Result obtained with 0%-noisy TRAC data. (d) Result obtained with 20%-noisy TRAC data.

Table 1 :

 1 Relative L 2 -error for nine different tests and for different levels of noise on the recording data.

	1	1b	2	2b	3	3b	4	5	6	Mean
	0% 6.28% 10.31% 3.14% 3.98% 1.93% 2.01% 3.92% 7.37% 5.18% 4.90%
	5% 6.39% 10.38% 3.37% 4.17% 2.31% 2.43% 4.15% 7.50% 5.31% 5.11%
	10% 6.72% 10.63% 4.03% 4.74% 3.21% 3.32% 4.68% 7.84% 5.78% 5.66%
	20% 8.04% 11.47% 5.98% 6.42% 5.32% 5.41% 6.28% 8.87% 7.21% 7.22%

Acknowledgements

The authors want to thank warmly Frédéric Nataf and Franck Assous for their support and advices throughout this work.