EN FR
EN FR


Section: New Results

Understanding the interior of the Earth and the Sun by solving inverse problems

Time-Domain Full Waveform Inversion based on high order discontinuous numerical schemes

Participants : Hélène Barucq, Julien Diaz, Pierre Jacquet.

Full Waveform Inversion (FWI) allows retrieving the physical parameters (e.g. the velocity, the density) from an iterative procedure underlying a global optimization technique. The recovering of the medium corresponds to the minimum of a cost function quantifying the difference between experimental and numerical data. In this study we have considered the adjoint state method to compute the gradient of this cost function.

The adjoint state can be both defined as the adjoint of the continuous equation or the discrete problem. This choice is still under study and complementary results has been presented at WAVES 2019 conference in Vienna [38].

The FWI has been largely developed for time-harmonic wave problems essentially because of computational time which is clearly below the one of corresponding time-dependent problems. However, the memory cost in large 3D domain is overflowing the computer capabilities, which motivates us to develop a FWI algorithm in the time-domain. To fully exploit the information from the seismic traces, while preserving the computational cost, it is important to use an accurate and flexible discretization. For that purpose we study several time schemes such as Runge-Kutta 2/4 or Adam-Bashforth 3 and regarding the space discretization, we employ Discontinuous Galerkin (DG) elements which are well-known not only for their h and p adaptivities but also for their massively parallel computation properties.

In the work-flow of DIP a Reverse Time Migration (RTM) code has been developed in collaboration with Total using their Galerkin Discontinuous acoustic time domain solver. Then this code served as a prototype of the time domain FWI code called utFWI (Unstructured Time-Domain Full Waveform Inversion) (https://bil.inria.fr/fr/software/view/3740/tab). Thanks to this code, 2D acoustic multi-scale reconstructions has been performed. Several optimizers such as gradient descent, non linear conjugate gradient and limited BFGS have also been developed.

This work is a collaboration with Henri Calandra (TOTAL). The time domain FWI results has been presented at Total conference MATHIAS 2019 in Paris [37] and also during the Fall Meeting 2019 AGU in San Francisco [47].

Box-Tomography imaging in the deep mantle

Participant : Yder Masson.

Box Tomography is a seismic imaging method (Masson and Romanowicz, 2017) that allows the imaging of localized geological structures buried at arbitrary depth inside the Earth, where neither seismic sources nor receivers are necessarily present. The big advantage of box-tomography over standard tomographic methods is that the numerical modelling (i.e. the ray tracing in travel time tomography and the wave propagation in waveform tomography or full waveform inversion) is completely confined within the small box-region imaged. Thus, box tomography is a lot more efficient than global tomography (i.e. where we invert for the velocity in the larger volume that encompasses all the sources and receivers), for imaging localized objects in the deep Earth. Following a successful, yet partial, application of box tomography to the imaging of the North American continent (i.e. Clouzet et al, 2018), together with Barbara Romanowicz and Sevan Adourian at the Berkeley Seismological Laboratory, we finished implementing the necessary tools for imaging localized structure in the Earth's lower mantle. The following tasks have been completed:

  • Modify the global wave propagation solver Specfem_3D_globe in order to compute Green's functions in our current reference Earth model (SEMUCB).

  • Modify the local wave propagation solver RegSEM_globe so that the wavefield can be recorded and stored at the surface of the modeling domain.

  • Implement real-time compression for an improved management of computed data.

Preliminary results have been presented at the American geophysical union fall meeting 2019 [24]. In the near Future, the latest implementation of Box-Tomography will be deployed to investigate the deep structure under the Yellowstone hotspot down to 20 seconds period.

Wave-propagation modeling using the distributional finite difference method (DFD).

Participant : Yder Masson.

In the last decade, the Spectral Element Method (SEM) has become a popular alternative to the Finite Difference method (FD) for modeling wave propagation in heterogeneous geological media. Though this can be debated, SEM is often considered to be more accurate and flexible than FD. This is because SEM has exponential convergence, it allows to accurately model material discontinuities, and complex structures can be meshed using multiple elements. In the mean time, FD is often thought to be simpler and more computationally efficient, in particular because it relies on structured meshed that are well adapted to computational architectures. This motivated us to develop a numerical scheme called the Distributional Finite Difference method (DFD), which combines the efficiency and the relative simplicity of the finite difference method together with an accuracy that compares to that of the finite/spectral element method. Similarly to SEM, the DFD method divides the computational domain in multiple elements but their size can be arbitrarily large. Within each element, the computational operations needed to model wave propagation closely resemble that of FD which makes the method very efficient, in particular when large elements are employed. Further, large elements may be combined with smaller ones to accurately mesh certain regions of space having complex geometry and material discontinuities, thus ensuring higher flexibility. An exploratory code allowing to model 2D and 3D wave propagation in complex media has been developed and demonstrated the interest of the proposed scheme. We presented numerical examples showing the accuracy and the interest of the DFD method for modeling wave propagation through the Earth at the EGU general assembly, 2019 and at the AGU fall meeting 2019 [41].

Time-harmonic seismic inverse problem

Participants : Hélène Barucq, Florian Faucher.

We study the inverse problem associated with the propagation of time-harmonic wave. In the seismic context, the available measurements correspond with partial reflection data, obtained from one side illumination (only the Earth surface is available). The inverse problem aims at recovering the subsurface Earth medium parameters and we employ the Full Waveform Inversion (FWI) method, which employs an iterative minimization algorithm of the difference between the measurement and simulation.

In particular, we investigate the use of new misfit functionals, based upon the reciprocity-gap. The use of such functional is only possible when specific measurements are available, and relates to Green's identity. The feature of the cost function is to allow a separation between the observational and numerical sources. In fact, the numerical sources do not have to coincide with the observational ones, offering new possibilities to create adapted computational acquisitions, and possibilities to reduce the numerical burden.

This work is a collaboration with Giovanni Alessandrini (Università di Trieste), Maarten V. de Hoop (Rice University), Romina Gaburro (University of Limerick) and Eva Sincich (Università di Trieste).

This work has given rise to a publication [11] and a preprint [53]. It has also been presented in several conferences [32], [35], [33], [34].

Convergence analysis for the seismic inverse problem

Participants : Hélène Barucq, Florian Faucher.

The determination of background velocity by Full Waveform Inversion (FWI) is known to be hampered by the local minima of the data misfit caused by the phase shifts associated with background perturbations. Attraction basins for the underlying optimization problems can be computed around any nominal velocity model and guarantee that the misfit functional has only one (global) minimum. The attraction basins are further associated with tolerable error levels which represent the maximal allowed distance between the (observed) data and the simulations (i.e., the acceptable noise level). The estimates are defined a priori, and only require the computation of (possibly many) first and second order directional derivatives of the (model to synthetic) forward map. The geometry of the search direction and the frequency influence the size of the attraction basins, and complex frequency can be used to enlarge the basins. The size of the attraction basins for the perturbation of background velocities in the classical FWI (global model parametrization) and the data space reflectivity (MBTT) reformulation are compared: the MBTT reformulation increases substantially the size of the attraction basins. Practically, this reformulation compensates for the lack of low frequency data.

This work is a collaboration with Guy Chavent (Inria Rocquencourt) and Henri Calandra (TOTAL). The results have been published in Inverse Problems  [12] and in the Research Report [43].

Eigenvector representation for the seismic inverse problem

We study the seismic inverse problem for the recovery of subsurface properties in acoustic media. In order to reduce the ill-posedness of the problem, the heterogeneous wave speed parameter is represented using a limited number of coefficients associated with a basis of eigenvectors of a diffusion equation, following the regularization by discretization approach. We compare several choices for the diffusion coefficient in the partial differential equations, which are extracted from the field of image processing. We first investigate their efficiency for image decomposition (accuracy of the representation with respect to the number of variables and denoising). Next, we implement the method in the quantitative reconstruction procedure for seismic imaging, following the Full Waveform Inversion method.

This work is a collaboration with Otmar Scherzer (University of Vienna) and the results have been documented in [49].

Outgoing solutions for the scalar wave equation in Helioseismology

Participants : Hélène Barucq, Florian Faucher, Ha Pham.

We study the construction and uniqueness of physical solutions for the time-harmonic scalar wave equation arising in helioseismology. The definition of outgoing solutions to the equation in consideration or their construction and uniqueness has not been discussed before in the context of helieoseismology. In our work, we use the Liouville transform to conjugate the original equation to a potential scattering problem for Schrödinger operator, with the new problem containing a Coulomb-type potential. Under assumptions (in terms of density and background sound speed) generalizing ideal atmospheric behavior, we obtain existence and uniqueness of variational solutions.

Solutions obtained in this manner are characterized uniquely by a Sommerfeld-type radiation condition at a new wavenumber. The appearance of this wavenumber is only clear after applying the Liouville transform. Another advantage of the conjugated form is that it makes appear the Whittaker special functions, when ideal atmospheric behavior is extended to the whole domain. This allows for the explicit construction of the outgoing Green kernel and the exact Dirichlet-to-Neumann map and hence reference solutions and radiation boundary condition.

This work has given rise to a report of 135 pages, [45]. Some part have been extracted for a publication which has recently been accepted in ESAIM: M2AN.

Consequently to this work, ongoing research includes the fast construction of the Green's kernel, which is possible thanks to the family of special functions obtained from the previous analytical study, [36]. As part of the ANTS associate team, applications to helioseismology is also ongoing. This work is a collaboration with Damien Fournier and Laurent Gizon (Max Planck Institute at Göttingen). The “Ants workshop on computational helioseismology” has also been organized in this context.

Modeling the propagation of acoustic wave in the Sun

Participants : Hélène Barucq, Juliette Chabassier, Marc Duruflé, Nathan Rouxelin.

We study time-harmonic propagation of acoustic waves in the Sun in the presence of gravity forces.

Galbrun’s equation, a Lagrangian description of acoustic wave propagation, is usually used in helioseismology. However, the discretization of this equation with high-order discontinuous Galerkin methods leads to poor numerical results for solar-like background flow and geometries. As better numerical results were obtained by using another model, the Linearized Euler’s Equations, we investigate the equivalence between those two models.If compatible boundary conditions are chosen, it should be possible to reconstruct the solution of Galbrun’s equation by solving the Linearized Euler’s Equations and then a vectorial transport equation. It turns out that finding those boundary conditions is quite difficult and not always possible.

We also have constructed a reduced model for acoustic wave propagation in the presence of gravity. Under some additional assumptions on the background flow and for high frequencies, the Linearized Euler’s Equations can be reduced to a scalar equation on the pressure perturbation. This equation is well-posed in a usual variational framework and it will provide a useful reference solution to validate our numerical methods. It also seems that a similar process could be used in more realistic cases.

Equivalent boundary conditions for acoustic media with exponential densities. Application to the atmosphere in helioseismology

Participants : Juliette Chabassier, Marc Duruflé, Victor Péron.

We present equivalent boundary conditions and asymptotic models for the solution of a transmission problem set in a domain which represents the sun and its atmosphere. This problem models the propagation of an acoustic wave in time-harmonic regime. The specific non-standard feature of this problem lies in the presence of a small parameter δ which represents the inverse rate of the exponential decay of the density in the atmosphere. This problem is well suited for the notion of equivalent conditions and the effect of the atmosphere on the sun is as a first approximation local. This approach leads to solve only equations set in the sun. We derive rigorously equivalent conditions up to the fourth order of approximation with respect to δ for the exact solution u. The construction of equivalent conditions is based on a multiscale expansion in power series of δ for u. Numerical simulations illustrate the theoretical results. Finally we measure the boundary layer phenomenon by introducing a characteristic length that turns out to depend on the mean curvature of the interface between the subdomains. This work has been published in Applied Mathematics and Computation [15].