Section: New Results

Inverse Problems

Reconstruction of an elastic scatterer immersed in a homogeneous fluid

Participants : Hélène Barucq, Rabia Djellouli, Élodie Estécahandy.

The determination of the shape of an obstacle from its effects on known acoustic or electromagnetic waves is an important problem in many technologies such as sonar, radar, geophysical exploration, medical imaging and nondestructive testing. This inverse obstacle problem (IOP) is difficult to solve, especially from a numerical viewpoint, because it is ill-posed and nonlinear. Its investigation requires as a prerequisite the fundamental understanding of the theory for the associated direct scattering problem, and the mastery of the corresponding numerical solution methods.

In this work, we are interested in retrieving the shape of an elastic obstacle from the knowledge of some scattered far-field patterns, and assuming certain characteristics of the surface of the obstacle. The corresponding direct elasto-acoustic scattering problem consists in the scattering of time-harmonic acoustic waves by an elastic obstacle Ωs embedded in a homogeneous medium Ωf, that can be formulated as follows:

Δ p + ( ω 2 / c f 2 ) p = 0 in Ω f · σ ( u ) + ω 2 ρ s u = 0 in Ω s ω 2 ρ f u · n = p / n + e i ( ω / c f ) x · d / n on Γ σ ( u ) n = - p n - e i ( ω / c f ) x · d n on Γ lim r + r p / r - i ( ω / c f ) p = 0 (1)

where p is the fluid pressure in Ωf whereas u is the displacement field in Ωs, and σ(u) represents the stress tensor of the elastic material.

This boundary value problem has been investigated mathematically and results pertaining to the existence, uniqueness and regularity can be found in [92] and the references therein, among others. We have obtained a new result proving the well-posedness of the problem when the fluid-solid interface is only lipschitzian. This has been published in the Journal of Mathematical Analysis and Applications [20] . We then propose a solution methodology based on a regularized Newton-type method for solving the IOP. The proposed method is an extension of the regularized Newton algorithm developed for solving the case where only the Helmholtz equation is involved, that is the acoustic case by impenetrable scatterers [86] . The direct elasto-acoustic scattering problem defines an operator F:Γp which maps the boundary Γ of the scatterer Ωs onto the far-field pattern p. Hence, given one or several measured far-field patterns p˜(x^), corresponding to one or several given directions d and wavenumbers k, one can formulate IOP as follows:

Find a shape Γ such that F ( Γ ) ( x ^ ) = p ˜ ( x ^ ) ; x ^ S 1 .

At each Newton iteration, we solve the forward problem using a finite element solver based on discontinuous Galerkin approximations, and equipped with high-order absorbing boundary conditions. We have first characterized the Fréchet derivatives of the scattered field and the characterization has been published in the Journal of Inverse and Ill-posed problems [18] . It is worth noting that they are solutions to the same boundary value problem as the direct problem with other transmission conditions. This work has been the object of several talks [63] , [50] , [36] . Elodie Estécahandy has defended her PhD thesis [14] in September 2013 and two papers will be submitted soon.

hp-adaptive inversion of magnetotelluric measurements

Participants : Hélène Barucq, Julen Alvarez Aramberri, David Pardo.

The magnetotelluric (MT) method is a passive electromagnetic exploration technique. It makes use of natural electric fields which propagate permanently into the Earth. Electric fields induce magnetic waves which can be detected at the surface to produce a map of the subsurface from the determination of the resistivity distribution. Magnetotelluric method is based on the mathematical relation between the magnetic and telluric variations which involve the electric resistivity of the subsurface. It is particularly relevant for the detection of metallic ores and for the study of geothermal sites. It is also used for oil and gas exploration because it provides information on sedimentary basins. It performs well on depth scales varying between a few tens of meters to hundred of kilometers, following the pioneering works of Tikhonov and Cagniard . Magnetotelluric measurements are governed by polarized Maxwell's equations in such a way that Helmholtz equations have to be solved. The geological mapping is constructed from the solution of an Inverse problem which requires computing the Impedance and/or the Resistivity distributions. In this work, we assimilate Earth with a horizontally layered model with possible 2D heterogeneities. Both the size of the direct problem and the required computational times may be excessively large. Indeed, on the one hand, the model of the source requires defining a horizontally sufficiently large thick plate to avoid undesirable effects that could take place around the edges. On the other hand, the inversion of MT measurements typically requires the computation of an accurate solution at the receivers located at different positions. Since traditional hp-goal oriented techniques [98] , [97] provide an accurate solution in one single point, we use a multi-goal-oriented algorithm [99] to obtain accurate solutions at all receivers. To get accurate quantities at several positions, it is necessary to increase the size of the mesh. This induces high computational costs in particular because the solution of the inverse problem is based on reiterated solutions of the direct problem. To decrease the computational costs required to perform the inversion, we propose an adaptive multi-dimensional inversion algorithm, which consists in increasing step by step the dimension in which the direct problem and the inversion are solved. At first step, we compute the 1D primary field with a semi-analytical solution and we invert the 1D problem. After that, we introduce the 2D heterogeneities. Regarding the direct problem, we compute the secondary field, thereby, drastically reducing the size of the computational domain for this problem. Then, we perform the inversion using the solution to the 1D Inverse Problem as a regularization term, increasing the robustness of the inversion algorithm.

Reverse Time Migration with Elastic Wave Equations

Participants : Hélène Barucq, Henri Calandra, Julien Diaz, Jérôme Luquel.

Even if RTM has enjoyed the tremendous progresses of scientific computing, its performances can still be improved, in particular when applied to strong heterogeneous media. In this case, images have been mainly obtained by using direct arrivals of acoustic waves and the transition to elastic waves including multiples is not obvious, essentially because elastic waves equations are still more computationally consuming. We have thus chosen to consider high-order Discontinuous Galerkin Methods which are known to be well-adapted to provide accurate solutions based upon parallel computing. Now, one of the main drawback of RTM is the need of storing a huge quantity of information which is prohibitive when using elastic waves. For that purpose, we apply the Griewank algorithm  [88] following Symes' ideas  [101] for the acoustic RTM. The idea is to find a compromise between the number of wave equations to solve and the number of numerical waves that we have to store. This is the so-called Optimal Checkpointing. By reducing the occupancy of the memory, RTM should be efficient even when using elastic waves. The next step is the derivation of accurate imaging conditions, which could take advantage of all the information contained in the elastic wavefield. For acoustic media, Claerbout  [83] proposed an imaging condition which is widely used and turns out to be sufficient to accurately reproduce interfaces. But Claerbout conditions do not take wave conversions into account and, since P-wave and S-wave interact with each other, it might be relevant to use an imaging condition including these interactions. This has been done successfully by J. Tromp and C. Morency  [102] for seismology applications based upon the inversion of the global Earth. Their approach is based upon the adjoint state and it involves sensitivity kernels which are defined from the propagated and the back-propagated fields. Now, it has been shown in  [93] that full wave form inversions using these sensitivity kernels may be polluted by numerical artifacts. One solution is to use a linear combination of the sensitivity kernels to delete artifacts. In this work, we propose then a new imaging condition which construction is inspired from  [93] with some approximations required to keep admissible computational costs. We illustrate the properties of the new imaging condition on industrial benchmarks like the Marmousi model. In particular, we compare the new imaging condition with other imaging conditions by using as criteria the quality of the image and the computational costs required by the RTM. The results will be presented at the 2014 ECCOMAS conference in Barcelona.