Miao Yu

Entropic Unbalanced Optimal Transport : Application to Full-Waveform Inversion and Numerical Illustration

Keywords: Transport-optimal entropique non-équilibré, divergence de Sinkhorn, inversion de formes d'onde, Unbalanced entropic optimal transport, Sinkhorn divergence, full-waveform inversion

Les méthodes de tomographie sismique visent à inférer les propriétés physique et reconstruire le "modèle", i.e. les structures de l'intérieur de la Terre, à partir des ondes mécaniques -radiées par des sources naturelles ou anthropogéniques -enregistrées par des récepteurs en surface sous la forme de sismogrammes. Les méthodes d'inversion de formes d'onde ont été activement développées dans les contextes académiques et industriels et sont devenus des outils puissants pour améliorer l'estimation des propriétés physiques et des structures d'objets géologiques depuis les échelles globales jusqu'aux échelles locales de la géophysique d'exploration. L'inversion de formes d'onde est formulée comme un problème d'optimisation non linéaire, associé à un système d'équations aux dérivées partielles. Il est classiquement résolu par des méthodes d'optimisation locale via la minimisation itérative d'une fonction coût qui mesure la diérence entre les sismogrammes observés et synthétiques, et utilisent des méthodes d'état adjoint. Les méthodes d'état adjoint permettent le calcul des dérivées de la fonction coût par rapport aux paramètres du modèle en combinant le champ d'onde direct et un champ d'onde adjoint gouverné par un système d'équations adjointes et des conditions adjointes complémentaires. Les méthodes d'inversion de forme d'onde, qui inversent simultanément les courtes et grandes longueurs d'onde, sourent malheureusement en pratique de dicultés qui restreignent leur utilisation pratique. Leurs capacités se détériorent du fait du décit en basse fréquence des observations et d'un bon modèle initial. Des limitations qui sont associées à la nature mal posée du problème inverse qui peut facilement être piégé dans un minimum local.

Une direction proposée, an de réduire la dépendance vis à vis du modèle initial, est de remplacer la fonction classique, basée sur une distance de type moindres carrés, par de nouvelles fonctions coûts, pouvant impliquer une transformation non linéaire du signal, an de promouvoir la convexité et élargir le bassin d'attraction du minimum global.

La théorie du Transport Optimal (OT) a récemment été utilisée dans le cadre des problèmes inverses et de l'apprentissage automatique. Le transport optimal généralise les propriétés de la distance euclidienne au carré à l'espace des distributions de probabilité. La valeur optimale (au carré) du transport lui-même dénit une distance appelée distance 2-Wasserstein. Cette quantité est à nouveau convexe mais maintenant sur l'ensemble des distributions de probabilité.

Le Transport Optimal est déjà utilisé en FWI, cette thèse en fait partie. L'approche OT est encore largement ouverte sur trois fronts : les formes d'onde sismiques ne sont ni positives ni de masse totale normalisée. La convexité par rapport au modèle n'est pas garantie et nalement le calcul réel de la distance OT est coûteuse.

Dans ce travail, nous utilisons et combinons -d'un point de vue académique -deux extensions récentes d'OT dans le contexte FWI. D'abord la distance OT "non-équilibrée", qui dénit rigoureusement une distance sur l'ensemble des mesures de Radon positives contournant ainsi le problème de normalisation des données (mais pas le problème de positivité). Puis le cadre du Transport Optimal entropique et en particulier la variante simple et facile à calculer appelée divergence de Sinkhorn fournissant une bonne approximation de la distance 2-Wasserstein. La divergence de Sinkhorn peut être naturellement étendue au transport "non-équilibré".

Nous utilisons ces outils pour construire et mettre en oeuvre une fonction coût OT "non-équilibrée". Nous discutons de son utilisation dans le contexte FWI au travers d'un certain nombre d'exemples académiques et de problèmes de référence classiques.

Entropic Unbalanced Optimal Transport :

Application to Full-Waveform Inversion and Numerical Illustration par Miao Yu

Résumé étendu en français

La tomographie sismique vise à inférer un modèle des propriétés physiques et des structures de l'intérieur de la Terre à partir des ondes mécaniques rayonnées par des sources sismiques naturelles ou articielles, et enregistrées à la surface par des récepteurs sous la forme de sismogrammes, également appelés traces sismiques. Il s'agit d'un problème important en géophysique avec des applications sociétales et industrielles telles que la découverte et l'exploitation de ressources énergétiques, l'investigation de sites pour le génie civil, les gisements minéraux et l'approvisionnement en eau souterraine, l'aléa sismique et l'atténuation des risques.

C'est également problème dicile et intrinsèquement mal posé. Les interactions entre les ondes sismiques et le milieu géologique dépendent de la bande de fréquence, de la durée et de l'énergie des sources ; des échelles d'hétérogénéité du milieu ; du trajet de propagation des ondes entre la source et les récepteurs. L'extraction d'information à partir des sismogrammes enregistrés en surface n'est donc pas une simple tâche en particulier du fait que ces interactions sont intégrées le long du trajet de propagation. Le développement de nouveaux systèmes d'acquisition et l'utilisation de nouvelles sources (naturelles et/ou articielles) ont favorisé de nouvelles méthodes d'imagerie sismique haute résolution exploitant une quantité et une diversité de données de plus en plus importantes.

La problème d'inversion de forme d'onde (FWI) est formulé comme un problème d'optimisation non linéaire associé à une fonction coût mesurant les diérences entre données observées et prédites, et qui est contraint par un système d'équations aux dérivées partielles. Celles-ci modélisent (sous certaines approximation et paramétrisation physiques) la propagation du champ d'onde et dénissent implicitement une application de l'espace des paramètres du modèle dans l'espace des données. En pratique, le problème d'optimisation est résolu de manière itérative dans le cadre de méthodes d'optimisation locale de type gradient.

Au cours des dernières décennies, les méthodes FWI ont été activement étudiées et ont tirées parti : (1) du développement de nouvelles méthodes de modélisation physique et numérique de la propagation 3D du champ d'ondes -de l'acoustique au viscoélastique -dans des milieux hétérogènes et complexes ; [START_REF] Gangbo | The geometry of optimal transportation[END_REF] de l'augmentation rapide de la puissance et des capacités du calcul haute performance ; (3) de la quantité croissantes des données large bande.

Les méthodes FWI inversent simultanément, e.g. sans séparation d'échelle, les courtes et grandes longueurs d'onde du modèle. Cependant ces méthodes sourent dans la pratique de dicultés qui restreignent leur utilisation en production. Le problème d'optimisation peut être piégé dans un minimum local, éventuellement éloigné du minimum global, du fait de : la non linéarité du problème d'optimisation qui le rend dépendant de la connaissance a-priori du modèle initial ; la bande de fréquence limitée des observations induite par la nature des sources et la géométrie du système d'acquisition, et en particulier un decit en basses fréquences ; la présence de bruit dans les observations dont l'origine est multiple (environnemental, instrumental) et relativement important dans les basses fréquences ; la complexité de la physique des ondes sismiques et de la faible sensibilité de la onction coût vis à vis des composants du modèle dont les nombres d'onde sont plus petits que les longueurs d'ondes sismiques résolues par l'opérateur de propagation.

La thèse présente un rapide aperçu des méthodes d'imagerie sismique et des fondements de la méthode FWI, ainsi que de la méthode de l'état adjoint qui permet dans le cadre des méthodes d'optimisation locale un calcul ecace du gradient de la fonction coût via la corrélation point à point du champ d'onde incident prédit par les équations d'état du système et du champ d'onde adjoint régi par un système d'équations adjointes. La méthode est illustrée dans le cadre de de l'approximation acoustique qui est celui de cette thèse. Les diérentes limitations associées à la méthode FWI sont discutées avec une rapide revue des diérentes approches développées pour atténuer ces limitations. 4.3 (top) The input densities p (blue) and q (curve) as functions of x and y. (center) Evolution of the couplings l at iteration l of the Sinkhorn algorithm. (bottom) solution ?

List of Figures

" of problem (4.4) for several values of ". 5.1 Illustration of the effect of for the transformation P defined by (5.4), applied to the function f (t)=t, 8t 2 [ 1,1], which changes sign at t =0 .A s increases, P (f ) is getting smoother and smoother in the vicinity of t =0 , but increasing artificial mass is getting created that lead to possible bias in optimal transport. ............................. (5,10,30,35Hz) . Bottom figures: the negative (bottom-left figure )   and the positive (bottom-right figure) parts of the two wavelets. ................. 5.9 Comparison between the L 2 (solid blue curve), the balanced (normalised) S ✏, (solid purple curve), and the unbalanced OT ✏, (solid red curve) and S ✏, (solid yellow curve) misfit functions as a function of the time shift between Ricker and Ormsby wavelets for decreasing ✏ and . . La thèse fournit un rapide aperçu des fondements de ma méthode FWI, et de l'état adjoint, avec une illustration dans le cade de l'approximation acoustique scalaire, ainsi qu'un revue des différentes méthodes proposées pour prévenir les problèmes associés à cette méthodes, dans l'espace des données et/ou des images.

Cette thèse est une contribution à une nouvelle approche qui se développe activement aujourd'hui pour la formulation de nouvelles fonctions coût des les problèmes inverses (et l'apprentissage statistique) basés sur les distances de transport optimal (OT). Les distances OT transforment la mesure intensive locale des moindres carrés en une mesure globale, sensible a tous les défauts de correspondance en espace, temps et amplitude. Les distances basées sur l'OT sont de nature lagrangienne et en particulier convexes par rapport aux translations et aux dilatations. Cette propriété en fait de bons candidats pour résoudre le problème du "cycle skipping".

A ce jour, plusieurs familles de fonctions coût basées sur l'OT ont été proposées : Dans cette thèse, nous avons développé une troisième famille de fonctions coût, continues et différentiables, basées sut l'OT issues en combinant deux extensions récentes:

• La variante OT « déséquilibrée » introduite récemment et qui définit rigoureusement une distance sur l'espace des mesures de Radon positives contournant ainsi le problème de la normalisation des données.

Cette classe de distance OT a été déjà fait l'objet d'une application FWI par Li et collaborateurs en 2021.

• La méthode de régularisation entropique de OT qui permet une implémentation efficace à l'aide de l'algorithme de Sinkhorn. L'OT entropique ne définit pas réellement une distance mais fournit une approximation du second ordre de la distance 2-Wasserstein sous la forme dite de divergence • Nous supposons que la source sismique est ponctuelle et connue, alors qu'en pratique la source doit être déterminée ou calibrée, généralement à partir d'une arrivée directe.

• Nous supposons que les observables sont des échantillons ponctuels de la solution de l'équation des ondes acoustiques, plutôt que des versions filtrées (à la fois dans l'espace et dans le temps) des perturbations de pression mesurées par les récepteurs.

• Nous ne traitons pas la problématique des bruits d'acquisition dans les observations (ex : dysfonctionnement des détecteurs, bruit sismique ambiant, diffusion incohérente des structures qui ne sont pas

Chapter 1

Introduction

Seismic tomography aims at inferring physical properties and reconstructing quantitatively structures of the Earth interior from the mechanical waves that are radiated by natural and man-made seismic sources and recorded at the surface at receivers in the form of seismograms also called seismic traces. This is a fundamental problem in geophysics with applications to societal concerns such as the discovery and exploitation of energy resources, site investigation for civil engineering, mineral deposits and underground water supplies, seismic hazard and risk mitigation.

With respect to other physics-based osculating methodse.g. gravimetry, geomagnetism, and electromagnetism -seismic waves have large source-dependent penetration depths and provide high-resolution reconstruction due to their source-dependent short wavelengths. They are today recorded with increasing accuracy by extended arrays with a large number of receivers -seismometers, accelerometers, hydrophones, geophones, fibre optic cables -at the surface (land and sea) or near-surface (boreholes) of the Earth.

In seismology, passive source tomography that makes use of natural seismic sourcese.g. earthquakes radiating very energetic seismic waves -have been actively developed in the last decades, and provide significant insights of the deep Earth interior from global to regional scales. The main limitation being the distribution and the frequency content of the natural sources, together with the available data coverage at these scales. In exploration geophysics, the development of controlled man-made sources -e.g. explosive and vibro-seis sources on land, air-gun sources at sea -and of long and very long offset acquisition systems has foster the development of high-resolution seismic imaging methods exploiting increasingly large amount of data. This led to higher and higher resolution subsurface models at local scales, from dozens of meters to dozens of kilometers, depending on the source energy and the acquisition system geometry.

Interactions between the waves and the geological medium depend on: the limited frequency-band, the duration and radiated energy of the seismic sources; the characteristic heterogeneity scales of the medium;

and of the source-receiver wave propagation path and distance. Extracting information from the recorded seismograms is not a trivial task as these interactions get integrated over the propagation distance. Usable data in seismic tomography have long been restricted to secondary observables such as travel times, phase speeds or waveforms -a small portion of the full recorded seismograms -by the limited capability of the physics and numerical wave propagation modelling.

Most discoveries related to the structure of the Earth interior -such as the Earth inner core, asthenosphere and the major seismic discontinuities -have long been based on ray theory describing short-wavelength energy propagation and travel-time attributes. In the ray theoretical framework, the arrival time of the seismic phases are sensitive to the wave speeds along the ray path connecting source and receiver. Travel time tomography [START_REF] Aki | Determination of the three-dimensional seismic structure of the lithosphere[END_REF][START_REF] Luo | Wave-equation travel time inversion[END_REF][START_REF] Pratt | Combining wave-equation imaging with travel time tomography to form high-resolution images from crosshole data[END_REF][START_REF] Schuster | Wavepath eikonal travel time inversion: Theory[END_REF][START_REF] Trampert | Global phase velocity maps of love and raileigh waves between 40 and 150 seconds[END_REF][START_REF] Ekström | Measurements and global models of surface wave propagation[END_REF][START_REF] Nemeth | Dynamic smoothing in crosswell traveltime tomography[END_REF][START_REF] Zelt | Modelling strategies and model assessment for wide-angle seismic travel time data[END_REF][START_REF] Boschi | Whole earth tomography from delay times of p, pcp, and pkp phases: Latera heterogeneities in the outer core or radial anisotropy in the mantle?[END_REF]Raw linson and Sambridge , 2003;[START_REF] Romanowicz | Global mantle tomography: progress status in the past 10 years[END_REF] is a fast and cost effective tool.

It is however limited to well-isolated body-wave phases or surface waves, for which the fundamental and higher modes are well separated. High-frequency approximation does not take into account the source finite-frequency band, and is only applicable to smooth media, in which the characteristic scales of the inhomogeneities are much larger than the dominant wavelength. As such the velocity structure derived from travel-time tomography is only sub-optimal, and ray-based method can yield distorted results and fail at caustics [START_REF] Čeverný | Seismic, ray theory[END_REF].

Efforts to overcome the limitations of ray theory led to the development of finite-frequency travel-time tomography, a semi-analytical extension of the ray-based inverse problem accounting for spatially extended, finite-frequency 3D wave-sensitivity kernels, together with the integration of different data sets [START_REF] Čeverný | Fresnel volume ray tracing[END_REF][START_REF] Masters | A shear-velocity model of the [and discussion[END_REF][START_REF] Friederich | Propagation of seismic shear and surface waves in laterally heterogeneous mantle by multiple forward-scattering[END_REF][START_REF] Friederich | The s-velocity structure of the east asian mantle for inversion of shear and surface waveforms[END_REF][START_REF] Ritsema | Complex shear wave velocity structure imaged beneath africa and iceland[END_REF][START_REF] Dahlen | Fréchet kernels of finite-frequency traveltimes -i. theory[END_REF][START_REF] Romanowicz | Anomalous splitting of free oscillations: A reevaluation of possible interpretations[END_REF][START_REF] Gu | Models of the mantle shear velocity and discontinuities in the pattern of lateral heterogeneities[END_REF][START_REF] Montelli | Finite-frequency tomography reveals a variety of plumes in the mantle[END_REF]Yoshizawa andKennett , 2004, 2005;[START_REF] Boschi | Global multiresolution models of surface wave propagation: comparing equivalently regularized born and ray theoritical solutions[END_REF][START_REF] Sigloch | Two-stage subduction history under north america inferred from multiple-frequency tomography[END_REF]. In particular asymptotic finite-frequency kernels allowed to exploit weighted, lowamplitude energy wave packets in the seismograms [START_REF] Li | Global mantle shear velocity model developed using nonlinear asymptotic coupling theory[END_REF][START_REF] Mégnin | The three-dimensional shear velocity structure from the inversion of body, surface and higher-mode waveforms[END_REF][START_REF] Gung | Q tomography of the upper mantle using three-component long-period waveforms[END_REF] leading to improved global tomographic models.

In exploration seismology, more controlled setup and multifold acquisition systems have enabled accurate, passive seismic imaging methods making use of both travel times and reflectivity energy information. However limited offsets of the seismic reflection surveys and limited frequency band of the source make seismic imaging poorly sensitive to intermediate wavelengths [START_REF] Jannane | Wavelengths of earth structures that can be resolved from seismic reflection data[END_REF].

As such, two broad regimes of wave interactions were mainly considered: a transmission regime, in which wave interactions are dominated by the slow variations of the physical properties, and the seismic phases and their propagation direction are slightly perturbed with a dominant forward scattering; a reflectivity regime, in which wave interactions are dominated by the fast variations of the physical properties, and the seismic phases and their propagation directions are strongly perturbed with a significant back scattering. However, it is worth to note that other types of interactions such as multiple scattering may also contribute to the recorded traces.

As one phase is mixing different interaction regimes during their propagation, the recorded seismic traces contain both kinematic information and reflectivity information. Following the early breakthroughs by [START_REF] Claerbout | Toward a unified theory of reflector mapping[END_REF][START_REF] Claerbout | Fundamentals of geophysical data processing with applications to petroleum prospecting[END_REF], two-step seismic imaging workflows orchestrating a tomography mode and a migration mode were actively developed, assuming that reflectivity result from small perturbation in velocity under the Born approximation, and that data -up to a first-order perturbation -linearly depend on the short Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

wavelengths of the velocity model.

In tomography mode, long wavelength components of the velocity model are reconstructed using the kinematic information of the wave propagation. In migration mode, short wavelengths images of the subsurface structures, which appear as a discontinuity for a given seismic signal frequency, are reconstructed after kinematic corrections by amplitude summation and back-projection using different types of migration methods [START_REF] Claerbout | Downward continuation of moveout-corrected seismograms[END_REF][START_REF] Gazdag | Wave equation migration with the phase shift methode[END_REF][START_REF] Stolt | Migration by fourier transform[END_REF][START_REF] Baysal | Reverse time migration[END_REF][START_REF] Nemeth | Least-square migration of incomplete reflection data[END_REF][START_REF] Yilmaz | Seismic data analysis[END_REF][START_REF] Biondi | Angle-domain common-image gathers for migration velocity analysis by wavefield continuation imaging[END_REF] formulated in the time or the frequency domain.

The notion of full wave form inversion (FWI), which is the context of this thesis, was first introduced to the geophysical community in the early '80s together with the elegant and physically insightful adjoint-state method [START_REF] Bamberger | Une application de la théorie du contrôle à un problème inverse sismique[END_REF][START_REF] Bamberger | Inversion of normal incidence seismograms[END_REF]Lail ly , 1983;Tarantola , 1984aTarantola , ,b, 1987[START_REF] Tarantola | A theoretical background for the inversion of seismic waveforms including elasticity and attenuation[END_REF]. It was later extended to the frequency domain [START_REF] Pratt | Inverse inverse theory applied to multi-source cross-hole tomography. part 1: acoustic wave-equation method[END_REF][START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. part 2: Elastic waveequation method[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain, part i: Theory and verification in a physical scale model[END_REF][START_REF] Sirgue | Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF].

FWI goes one-step further by inverting simultaneously the short and long wavelengths with the potential to bridge the gap between the transmission and the reflectivity regimes, overcoming the limitations of the conventional sequential seismic-imaging workflow involving a tomography step followed by a migration step.

FWI is formulated as a nonlinear PDE-based optimisation problem that considers the full wavefield information to estimate the Earth interior properties. The adjoint-state methods allows the computation of the derivative of an objective function -measuring the difference between the observations and the predicted observables -with respect to the model parameters by combining the predicted incident wavefield governed by the wave equations, together with appropriate initial and boundary equations, and the adjoint wavefield governed by a set of adjoint equations, together with adjoint subsidiary terminal and boundary conditions.

In the last decades, FWI has been actively studied in both academia and industry, taking advantage of:

(1) the development of new physical and numerical modelling methods for 3D full wavefield propagationfrom acoustic to viscoelastic -in heterogeneous, complex media; (2) the rapidly increasing capability and capacity of high-performance computing;

(3) the increasing amount and coverage of broadband seismic data.

FWI has proven to be a powerful tool that dramatically improved the capability to estimate physical properties and structures of various geological targets from global to regional scales [START_REF] Komatitsch | The spectral-element method, beowulf computing, and global seismology[END_REF][START_REF] Tromp | Seismic tomography, adjoint methods, time reversal and bananadoughnut kernels[END_REF]Fichtner et al., 2006a,b;[START_REF] Tape | Finite-frequency tomography using adjoint methods -methodology and examples using membrane surface waves[END_REF][START_REF] Liu | Finite-frequency sensitivity kernels for global seismic wave propagation based upon adjoint methods[END_REF][START_REF] Fichtner | Theoretical background for continental -and global -scale full waveform inversion in the time-frequency domain[END_REF][START_REF] Fichtner | Ful l seismic waveform model ling and inversion[END_REF][START_REF] Peter | Forward and adjoint simulations of seismiv wave propagation on fully unstructured hexahedral meshes[END_REF][START_REF] Fichtner | Multiscale full waveform inversion[END_REF][START_REF] Zhu | Seismic structure of the european upper mantle based on adjoint tomography[END_REF][START_REF] Komatitsch | Anelastic sensitivity kernels with parsimonius storage for adjoint tomography and full waveform inversion[END_REF][START_REF] Tromp | Seismic wavefield imaging of the earth'sinterior across scales[END_REF] and local scales in exploration geophysics [START_REF] Pratt | Two-dimensional velocity models from wideangle seismic data wevefield inversion[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain -part 2: Faukt delineation in sediments using crosshole data[END_REF][START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF][START_REF] Brossier | Seismic imaging of complex onshore structures by 2d elastic frequency-domain full-waveform inversion[END_REF][START_REF] Plessix | Full waveform inversion of a deep water ocean bottom seismometer dataset[END_REF][START_REF] Sirgue | Full waveform inversion: the next leap forward in imaging at valhall[END_REF][START_REF] Plessix | Full waveform inversion and distance separated simultaneous sweeping: a study with a land seismic data set[END_REF][START_REF] Warner | Anisotropic 3d fullwaveform inversion[END_REF][START_REF] Stopin | Multiparameter waveform inversion of a large wide-azimuth low-frequency land data set in oman[END_REF][START_REF] Vigh | Elastic full-waveform inversion application using multicomponent measurements of seismic data collection[END_REF][START_REF] Operto | Efficient 3-d frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to valhall in the visco-acoustic vertical transverse isotropic approximationb[END_REF][START_REF] Shen | Full waveform inversion: The next leap forward in subsalt imagin[END_REF][START_REF] Borisov | Application of 2d full waveform inversion on exploration land data[END_REF].

In practice, however, solving FWI problems using local, nonlinear optimisation methods is facing challenges that preclude routine use.

The physics of seismic waves is complex and the objective functional is poorly sensitive to wavenumber components of the model that are shorter than seismic wavelengths and only through information that is homogenised by the wave propagation operator [START_REF] Cupillard | Non-periodic homogeneization of 3-d elastic media for the seismic wave equation[END_REF][START_REF] Capdeville | Elastic full-waveform inversion based on the homogeneization method: theoritical framework and 2-d numerical illustrations[END_REF][START_REF] Hedjazian | Multiscale seismic imaging with inverse homogenization[END_REF]. Accurate and computationally efficient models of wave propagation -from Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 pure acoustic waves to anisotropic viscoelasticity -together with well-defined physics-based and data-based parametrisation of the model need to be carefully considered in the FWI, specially in the case of multiparameter estimation and multi-mode modelling that remain very challenging [START_REF] Operto | A guided tour of multiparameter full-waveform inversion with multi-component data: From theoty to practice[END_REF]. Despite the rapidly increasing computational capabilities, numerical models can only process limited frequencies in 3D. While increasing the level of parametrisation allows more realistic models, the inversion problem becomes more ill-posed and even non-unique when the model space dimension increases.

The accuracy deteriorates from the lack of low frequencies, observation noise and poor starting models [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: Numerical results[END_REF][START_REF] Mora | Nonlinear two-dimensional elastic inversion of multioffset seismic data[END_REF][START_REF] Luo | Wave-equation travel time inversion[END_REF][START_REF] Fichtner | Theoretical background for continental -and global -scale full waveform inversion in the time-frequency domain[END_REF][START_REF] Virieux | An introduction to full waveform inversion[END_REF], unless some prior information is employed [START_REF] Bunks | Multiscale seismic inversion[END_REF]. In exploration geophysics, while ultra long-offsets allow to recover partially low frequencies, the seismic traces contain a full-variety of integrated, frequencydependent information associated with different resolution and interaction regimes during the propagation.

This makes the inversion problem highly nonlinear with multimodal objective function. When using local gradient-based optimisation methods, a poor initial model can cause the nonlinear optimisation problem to be trapped easily in a local minimum e.g. [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. Global optimisation methods -such as Monte Carlo [START_REF] Sambridge | Inverse theory, monte carlo method[END_REF][START_REF] Ray | Frequency domain full waveform elastic inversion of marine seismic data from the alba field using bayesian trans-dimensional algorithm[END_REF][START_REF] Biswas | 2d full-waveform inversion and uncertainty estimation using the reversible jump hamiltonian monte carlo[END_REF][START_REF] Ray | Low frequency full waveform seismic inversion within a tree based bayesian framework[END_REF], and Hamiltonian Monte Carlo (Fichtner et al., 2019;Fichtner and Zunino, 2019;[START_REF] Gebraad | Bayesian elastic full-waveform inversion using hamiltonian monte carlo[END_REF] sampling methods, simulated annealing [START_REF] Datta | Estmating a starting model for full-waveform inversion using global optimization method[END_REF], particle swarm optimisation (Chen and Wang, 2017) -could avoid theoretically this problem and quantify the uncertainties. However, owing to the large dimension of the model space, these methods remain currently intractable and far from reaching the computational efficiency needed to handle realistic large-scale seismic inversion problems in exploration geophysics.

Seismic data are oscillatory signals by nature, as the zero-frequency part is not recorded [START_REF] Chauris | Seismic imaging: a practical approach[END_REF]. Since the optimisation problem is highly nonlinear with multimodal objective functions a poor initial kinematic model may cause the solution to be trapped in local minima. In the standard form of FWI, the objective function, measuring the misfit between the predicted observables and the observations, is based on the least-squares distance, which is oscillatory and non convex in directions associated with wavenumber components of the model that are longer than a half the dominant wavelength. For large kinematic errors, local optimisation algorithms will match incorrect phases between the predicted and the observed signals, making the method prone to the so-called cycle-skipping effect. The inversion converges in a local minimum, possibly far from the global minimum, due to kinematic errors introduced by the poor starting model.

Different strategies have been developed to mitigate the non convexity of the FWI problem and to reduce the dependency on the initial model:

• Improve the kinematic accuracy of the initial model using for example reflection tomography [START_REF] Woodward | A decade of tomography[END_REF] or stereo-tomography methods [START_REF] Lambaré | Stereotomography[END_REF], so that the initial model is not too far from the correct model and the kinematic errors are reduced. Following a similar approach, it has also been proposed to enhance the tomographic components at the early iterations and gradually reduce its weights towards convergence [START_REF] Tang | Tomographically enhanced full wavefield inversion[END_REF], or to further separate the gradient components based on the scattering angle at the imaging point [START_REF] Alkhalifah | Scattering-angle based filtering of the waveform inversion gradients[END_REF]. Although these methods do enhance Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

the low wavenumber components of the FWI gradient, they reduce the flexibility of the FWI workflow and can be rather time consuming.

• Mitigate multi-modalities and enlarge the basin of attraction of the objective function using a data hierarchical approach: seeding it with low frequency data only, and slowly enlarge the data bandwidth as the iterations of the gradient-based methods progress [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part i: Theory and verification in a physical scale model[END_REF][START_REF] Shipp | Two-dimensional full wavefield inversion of wide-aperture marine seismic streamer data[END_REF][START_REF] Brossier | Seismic imaging of complex onshore structures by 2d elastic frequency-domain full-waveform inversion[END_REF][START_REF] Wang | Reflection seismic waveform tomography[END_REF]. Extension to the time domain involves different temporal and/or offset selection of the observations. The idea is to reduce the number of propagated wavelengths to be recovered simultaneously. Among the many factors that affect in practice the success of this approach, the lowest starting frequency is an essential one. Broadband data with high signal-to-noise ratio may not be always available or of sufficient quality preventing reliable model reconstruction. Moreover, smart selection of the data is required for the different inversion stages, which remains often rather empirical. Following similar approach, low frequency data extrapolation based on the phases and amplitudes in the observed band has also been proposed, fitting smooth non-oscillatory functions to represent and extrapolate wave physics to the unrecorded frequency band [START_REF] Li | Phase and amplitude tracking for seismic-event separation[END_REF].

Other solutions based on a reformulation of the FWI have been extensively studied in the past decades to reduce the non-linearity of the continuous PDE-based optimisation problem and to increase the convexity of the inversion in higher-dimensions. This has led to different extension strategies that aim to restore the dimensional balance between the model and the data space by considering artificial, extended models through the introduction of redundant coordinates in the model space. The common feature in all model extension strategies is their tendency to promote the same sort of long-wavelength velocity updates as does travel time tomography, i.e. to extract kinematic information from the data.

Enforcing a consistency condition that link the unphysical, extended model to the original physical one is widely viewed as a reasonable way of retrieving the background model in which waves propagate. There are at least two different ways of enforcing the consistency relation:

• A first strategy in the image (reflectivity) domain was originally developed within the framework of wave equation migration analysis [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF] and relies on a scale separation assumption. The idea is to introduce additional degrees of freedom at the reflectivity level extending the model space along subsurface offset, plane-wave ray-parameter, or time-lag [START_REF] Shen | Differential semblance velocity analysis by wave-equation migration[END_REF][START_REF] Shen | Wave equation migration velocity analysis by differential semblance optimization[END_REF][START_REF] Sava | Time-shift imaging condition in seismic migration[END_REF][START_REF] Sun | Waveform inversion via nonlinear differential semblance optimization[END_REF][START_REF] Biondi | Simultaneous inversion of full data bandwidth by tomographic fullwaveform inversion[END_REF]. The extended modelling operator is now surjective under very general conditions in 3D, and as such an extended model can always be found to match any data, regardless of wether the background kinematic model is correct. The construction delocalises scattering events and the extended image (reflectivity) should take on near-zero values because the wave fields are not supposed to interact constructively. The consistency condition is then derived by penalising the additional parameters in order to focus the energy of the extended reflectivity image [START_REF] Chauris | Two-dimensional velocity macro model estimation from seismic reflection data by local differential semblance optimization: applications to synthetic and real data sets[END_REF][START_REF] Biondi | Angle-domain common-image gathers for migration velocity analysis by wavefield continuation imaging[END_REF][START_REF] Sava | Time-shift imaging condition in seismic migration[END_REF] the additional computational cost to the already computationally demanding FWI, associated to the construction of the extended reflectivity images. Moreover since only primary reflections are used to generate the images, considerable preprocessing -such as removing the multiples, muting the direct arrivals, and diving waves -is required.

• A second strategy is the source-receiver extension strategy introducing the seismic source indexes as additional coordinates in the model parameter space [START_REF] Van Leeuwen | Mitigating local minima in full-waveform inversion by expanding the seach space[END_REF][START_REF] Warner | Adaptive waveform inversion: Theory[END_REF][START_REF] Huang | Full-waveform inversion via soure-receiver extension[END_REF][START_REF] Huang | Source-independent extended waveform inversion based on space-time source extension: Frequency-domain extension[END_REF], initially including trace-based Wiener filters to compare the observed and modelled data and later on the whole wavefield so that the reconstructed wavefield fits the data by design. The imaging condition involves incident and adjoint fields computed in the extended model independently for each source. The consistency condition again relies on the concept of differential semblance optimisation [START_REF] Symes | Velocity inversion by differential semblance optimization[END_REF], which is a good though not fully compelling choice of regularisation for model-extended waveform inversion. The applicability of these methods is by no mean feat when multiple arrivals are present in the transmission data [START_REF] Plessix | Automatic cross-well tomography by semblance and differential semblance optimization: Theory and gradient computation[END_REF], as well as their extension to the time domain [START_REF] Aghamiry | Accurate and efficient data-assimilated wavefield reconstruction in the time domain[END_REF]. Current solutions either rely on crude approximations [START_REF] Wang | Full waveform inversion with the reconstructed wavefield method[END_REF] or on advanced iterative solutions, which are computationally demanding [START_REF] Aghamiry | Accurate and efficient data-assimilated wavefield reconstruction in the time domain[END_REF] Apart from finding an adequate initial kinematic model, another line of investigation is to reformulate the FWI using alternative ways to measure the difference between predicted observables and observations in an attempt to mitigate the non convexity of the objective function with respect to the model parameters and to enlarge the basin of attraction of the global minimum. Different families of misfit function with enhanced convexity have been proposed in the past decades:

• A first family involves transforming the signal before measuring the difference between the predicted observable and observations in the sense of the least-squares distance. Examples of such transformation include: reconstruction of the envelope and the unwrapped phase of the signal [START_REF] Wu | Seismic envelope inversion and modulation signal model[END_REF][START_REF] Luo | Seismic envelope inversion: reduction of local minima and noise resistance[END_REF], extraction of the instantaneous phase and enveloppe of the signal by the Hilbert transform [START_REF] Fichtner | Theoretical background for continental -and global -scale full waveform inversion in the time-frequency domain[END_REF][START_REF] Bozdaǧ | Misfit functions for waveform inversion baseds on instantaneous phase and envelope measurements[END_REF][START_REF] Alkhalifah | Taming waveform inversion non-linearity through phase unwrapping of the model and objective function[END_REF][START_REF] Luo | Time-domain inversion using instantaneous phase information with damping[END_REF], normalised integration of the signal leading to trace cumulative distribution [START_REF] Liu | The normalised integration method -and alternative to full waveform inversion?[END_REF][START_REF] Donno | Estimating the background velocity model with the normalised integration method[END_REF].

• Another family involves more global measures in replacement of the least-squares distance. Correlationbased objective functions, promoting the minimisation of travel-time shifts, were initially proposed [START_REF] Luo | Wave-equation travel time inversion[END_REF][START_REF] Dahlen | Fréchet kernels of finite-frequency traveltimes -i. theory[END_REF][START_REF] Zhao | Three-dimensional fréchet differential kernels for seismicdelay time[END_REF]van Leeuwen andMulder , 2008, 2010;[START_REF] Choi | Application of multi-source waveform inversion to marine streamer data using the global correlation norm[END_REF], leading to the so-called wave equation tomography strategy. However, crosscorrelation-based measures tend to fail in the case of multiple energetic seismic arrivals. Matching filter-based objective functions [START_REF] Luo | A decovolution-based objective function for wave-equation inversion[END_REF][START_REF] Warner | Adaptive waveform inversion: Theory[END_REF][START_REF] Huang | Full-waveform inversion via soure-receiver extension[END_REF][START_REF] Huang | Source-independent extended waveform inversion based on space-time source extension: Frequency-domain extension[END_REF][START_REF] Zhu | Building good starting model for full-waveform inversion using adaptive matching filter misfit[END_REF] were then developed to circumvent these limitations, making use of a normalised deconvolution, i.e. a Wiener filter, of the predicted and observed seismic traces , leading to the so-called adaptive wave inversion method. The resulting objective function penalises the energy of Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

the filter away from a band-pass Dirac filter, which would have resulted from a correct model. However, these objective functions can still suffer from cycle-skipping and from a loss of resolution in the case of complex multi-arrival signals.

Matching-filter-based objective functions can be seen as an attempt to regularise the ill-conditioned inverse problem by controlling the slowly decaying energy associated to multiple arrivals and by removing effects of the small eigenvalues -associated to wavenumber components of the model that are shorter than seismic wavelengths -of the linearised forward wave operator.

Finally, a new trend has recently emerged in inverse problems and machine learning promoting optimal transport (OT) distances to measure the difference between the predicted observables and the observations, [START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport[END_REF]Yang and Engquist , 2018) to be compared using this metric. This can alter the phase and amplitude information, and prone non-convexity with respect to time shifts. Computing the 2-Wasserstein metric using the Monge-Ampère formulation either with optimised finite-difference, linear solvers [START_REF] Benamou | Numerical solution of the Optimal Transportation problem using the Monge-Ampère equation[END_REF][START_REF] Froese | A numerical method for the elliptic monge-amperè equation with transport boundary conditions[END_REF] or semi-discrete strategy [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF] add significant computational cost to the FWI and in practice application of the 2-Wasserstein distance has often been restricted 1D-optimal transport, i.e. trace-by-trace comparison, for which an analytical solution exists.

defining
• A second family, proposed by Métivier and collaborators (Métivier et al., 2016a,b;Métivier et al., 2016), resorts to the 1-Wasserstein metric, a particular instance of the Kantorowitch-Rubinstein metric, associated to the dual formulation of the optimal transport problem. This OT metric is now defined for signed measures and can be lifted to handle unbalanced optimal transport, avoiding ad-hoc signal transformation and normalisation. It leads to a convex optimisation problem under linear constraints that can be solved using proximal splitting methods [START_REF] Combettes | Proximal splitting methods in signal processing, in Fixed-point algorithms for inverse problems in science and engineering[END_REF], reducing significantly the additional computational compared to the 2-Wasserstein metric. The convexity of the objective function is however not guaranteed with regard to time and amplitude transformations.

• A third family recently introduced by Métivier et al [START_REF] Métivier | Optimal transport for mitigating cycle skipping in full-waveform inversion: A graph-space transform approach[END_REF][START_REF] Metivier | A graph space optimal transport distance as a generalisation of the l p distances: application to a seismic imaging inverse problem[END_REF] This new objective function shows promising convexity properties with respect to time and amplitude transformation, and in 2D the distance can be computed with moderate additional cost as a linear assignment problem. However extension to higher domain dimensions is still under development and remains to be evaluated, and it might increase significantly the computational cost.

To date, OT-based misfit functions have mostly been illustrated using 2D academic and benchmark cases, and only few of them have been practically demonstrated for field data [START_REF] Poncet | Fwi with optimal transport: a 3d implementation and an application on a field dataset[END_REF][START_REF] Messud | Multidomensional optimal transport for 3d fwi: Demonstration of field data[END_REF][START_REF] Sedova | Acoustic land full waveform inversion on a broadband land dataset: The impact of optimal transport[END_REF]Górszczyk et al., 2020Górszczyk et al., , 2021;;[START_REF] Pladys | On cycle-skipping function modification for fullwaveform inversion: comparison of five recent approaches[END_REF].

In this thesis, we developed a third family of OT-based, continuous and differentiable objective functional, derived from combining two recent extensions of OT:

• an "unbalanced" OT variant introduced in (Chizat et al., 2018a;[START_REF] Chizar | Unbalanced optimal transport: Dynamic and kantorovich formulations[END_REF][START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF]Kondratyev et al., 2016a), which rigorously defines a distance on the set of positive Radon measures thus bypassing the data normalisation issue. This class of OT distance was proposed but not investigated in details in [START_REF] Li | Application of an unbalanced optimal transport distance and a mixed l1/wasserstein distance to full waveform inversion[END_REF].

• an entropic regularisation of the OT framework and its efficient implementation using the Sinkhorn algorithm Cuturi (2013) [START_REF] Peyré | Computational optimal transport[END_REF]. Entropic OT does not actually define a distance but provides a second-order approximation of the 2-Wasserstein distance, and the so-called Sinkhorn divergence variant [START_REF] Genevay | Learning generative models with sinkhorn divergences[END_REF][START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF] leads to a very efficient computational implementation. The Sinkhorn divergence can be naturally extended to unbalanced OT [START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF].

We provide an overview of the theoretical framework of these two recent OT extensions leading to smooth and differentiable objective functions in the context of FWI. We numerically illustrate, in the 2D acoustic case, the use of these objective functions in the context of FWI with academic examples and the Marmousi benchmark data sets.

The numerical illustrations make however a number of physical assumptions for sake of simplicity that need to be kept in mind:

• We use the classical acoustic approximation of elastic waves and assume constant density in the 2D illustrations ; in this case, the sound speed velocity is the only model parameter. In the acoustic approximation mode conversion is not considered and the medium is isotropic. In practice, when dealing with heterogeneous media such an approximation can be challenged Cance and Capdeville (2015)

• We assume that the seismic source is ponctual and known, whereas in practice the source needs to be determined or calibrated, usually from a direct arrival.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

• We assume that the observables are point samples of the solution of the acoustic wave equation, rather than filtered versions (both in space and time) of pressure disturbances as measured by the receivers.

• We do not deal with the issue of acquisition noise in the observations (e.g. malfunctioning detectors, ambient seismic noise, incoherent scattering from structures that are not to be imaged).

As such the 2D numerical illustrations presented in this work has to be considered as academic illustrations that need to be further investigated and evaluated in more realistic contexts.

Outline of the manuscript

Chapter 2 (Seismic Background) presents the Full-waveform Inversion (FWI) problem of reconstructing a background model and the classical adjoint state method used to compute the gradient of the misfit function, and provides a discussion the so-called "cycle skipping" problem in the case of the classical least-squares misfit function.

Chapter 3 (OT Background) a general presentation of the Optimal Transport (OT) theory and of the main contributions using OT misfit functions in the FWI context is provided.

Chapter 4 (Entropic regularization of OT and its generalizations) presents recent branches of OT. First its entropic penalization and companion Sinhkhorn algorithm . Then its "unbalanced" extension allowing to define a bona fide distance between positive distributions with total mass. Finally how the two extension can be combined and how to correct the bias introduced by the entropic penalization (Sinkhorn Divergence).

Chapter 5 (The unbalanced optimal transport misfit function) summarizes the construction of the unbalanced misfit function from the model onward. We also illustrate its properties on two simple 1D and 2D parametric models.

A general setting

In the context of full-waveform inversion, the state of the dynamical system is determined by the wavefield

u(x,t):Ω ⇥ [0,T] ! R nu for some n u .W ea s s u m et h a tu belongs to a subset of W = L 2 ([0,T], [L 2 (Ω)] nu ),
The physical or computational domain Ω ⇢ R no can be open or bounded with boundary Γ.

The evolution of the dynamical system is governed, under some physical approximation, by a state equations, i.e. a system of partial differential equations, represented by the differential operator L,t h e behavior of which depends on model parameters m:

L (u(x,t), m(x)) f (x,t)=0, 8x 2 Ω, 8t 2 [0,T], (2.1)
where the seismic source f :

Ω ⇥ [0,T] ! R np is in W,a n dm : Ω !M⇢R nm the model parameters that
define a vector-valued field, with M the admissible model parameters space. The solution of (2.1) defines a physical realisation of the system u(m)=u(x,t; m) for fixed model parameters m.

The differential operator L governing the wavefield propagation, i.e. from acoustic to viscoelastic, is generally a linear operator:

L(u, m) ⌘L[m] u, (2.2)
where L [m] is the wave-propagation operator associated to the physical approximation.

When L [m] is of order Q in times, the state equations (2.2) must be complemented with Q 1 initial conditions at t =0(we assume a null initial state to simplify)

@ q u @t q (x,t = 0) = 0; 8x 2 Ω, 8q =0,..,Q 1, (2.3)
and possibly with appropriate boundary conditions along Γ when the physical domain is bounded.

In the data space D, observations d obs are seismograms recorded at n r receiverse.g. seismometers, hydrophones, geophones -for a source f :

d obs : {d obs (x r ,t),r=1,...,n r } ; t 2 [0,T] , (2.4)
where x r denotes the spatial position vector of the receiver r.

The predicted observables are defined accordingly as

d cal (m):{d cal (x r ,t; m),r=1,...,n r } = Ru(x,t; m); t 2 [0,T] , (2.5)
where R : W!Dis an extraction operator mapping the physical realisation u(m) onto the receiver positions

x r R : u(x,t; m) !{u(x 1 ,t; m),...,u(x nr ,t; m)} . (2.6)
Note that the definition of R supposes that u(x,t; m) is continuous at points x = x r ; this holds generally when the seismic source is regular enough and when the medium m(x) is constant around the receivers

x = x r
FWI is classically formulated as a nonlinear optimisation problem associated to an objective functional J : M!R attached to the observations and that acts on the physical observables:

J (m)=h (d cal (m), m; d obs )=h (Ru(m), m; d obs ) , (2.7)
where the functional h(•), which may depend explicitly on m, measures the difference between the observations and the predicted observables associated to the physical realisation u(m)=u(x,t; m).

2.2

The optimisation problem and its solution.

The optimisation problem can be stated as finding the optimal model parameters m ? such that J (m ? ) is the global minimum of J for a given set of observations: such as steepest descent, conjugate gradient or Newton-like methods can be formulated as:

m ⇤ = argmin m2M J (m) where J (m)=h (d cal (m), m; d obs ) . ( 2
m k+1 = m k + ↵ k g k with g k • r m J (m k ) < 0, (2.9)
where J (m k+1 ) < J (m k ),a n dg k is the local descent direction, with r m (•) the gradient operator with respect to the model parameters, known as the Fréchet derivative.

The step length ↵ k 2 R + , must be chosen by a line search process so that J (m k ↵ k g k ) is minimal, and must satisfies the strong Wolfe conditions [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF][START_REF] Wolfe | Convergence conditions for ascent methods. ii: some cirrections[END_REF], which are sufficient conditions to assure the convergence:

J (m k + ↵ k g k ) < J (m k )+c 1 ↵ k r m (m k ) • g k (2.10) |r m J (m k + ↵ k g k )| <c 2 |r m J (m k ) • g k | , (2.11) with 0 <c 1 <c 2 < 1.
The choice of the step length ↵ k is therefore a balance between stability and speed.

The steepest descent algorithm corresponds to g k = r m J (m k ) and the conjugate gradient algorithm

to g k = r m J (m k )/ k g k 1 , with k 0; whereas for Newton methods g k = H 1 J (m k )r m J (m k )
where the functional Hessian H J is defined as

H J (m k )= r m @J @m (m k ), (2.12)
and positive definite Hessian assures J (m k+1 ) < J (m k ).

Gradient descent typically converges slowly, which is a significant impediment for large-scale problems, whereas the Newton method, which takes into account the curvature of the objective function, converges faster than gradient descent methods in the neighborhood of a local minimum when the Hessian of J is close to being positive semi-definite.

Newton method is in general much more complicated to set up than gradient descent methods and much more computationally demanding as the Hessian of J is a large matrix, costly to store and to invert.

Practical alternatives are the quasi-Newton methods which attempt to partially invert the Hessian of J .

m k+1 = m k ↵ k Q k r m J (m k ), (2.13)
where Q k ⇡H 1 J , is symmetric and positive definite. Among the quasi-Newton methods, the l-BFGS [START_REF] Nocedal | Updating quasi-newton matrices with limited storage[END_REF][START_REF] Liu | On the limited memory bfgs method for large scale optimisation[END_REF], a limited memory implementation of the Broyden-Fletcher-Goldfarb-Shanno quasi-Newton method, leads to an efficient iterative algorithm for computing recursively Q k involving only the functional derivative of J at iteration k and at the l previous iterations.

The objective functional (2.7) can be relatively insensitive to wavenumber components of the model that are shorter than the seismic wavelengths. The eigenvalues of the Hessian of J corresponding to these components are nearly zero and the optimisation problem (2.8) becomes locally ill-posed.

Institut de Physique du Globe de Paris

Miao YU, Ph.D. Thesis, 2021

The eigenvectors associated to these eigenvalues are directions in which J (m) have very small curvature in the vicinity of m ? , the curvature being twice the second directional derivative in the eigen-directions. Small perturbations of the data or of the model m induce modifications of J that may result in large movements of its global minimum in problematic directions in the vicinity of the null-space of the Hessian of J .

All gradient-based minimisation algorithms critically rely on the ability to efficiently compute the directional derivative of the objective function with respect to the model parameters.

Adjoint-state methods and gradient computation

Adjoint-state methods [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF][START_REF] Chavent | Identification of functional parameters in partial differential equations[END_REF][START_REF] Chavent | History matching by use of optimal theory[END_REF][START_REF] Bamberger | Inversion of normal incidence seismograms[END_REF][START_REF] Tarantola | A theoretical background for the inversion of seismic waveforms including elasticity and attenuation[END_REF]Lail ly , 1983;[START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF][START_REF] Fichtner | Ful l seismic waveform model ling and inversion[END_REF] allow to compute the directional derivatives of J with optimal efficiency.

The adjoint wave equations can be derived from the wave equations, and the properties of the adjoint wave field, solution of the adjoint wave equations, are determined by the adjoint source, which is completely specified by the misfit function (2.7).

Generalisation of the adjoint method allows to compute the Hessian of J , which plays a fundamental role in Newton-like methods of non-linear minimisation.

The variation of the wave field u, a functional of m, with respect to m in the direction of m is given by the functional derivative, known as the Gâteau derivative:

u = r m u m, (2.14)
where r m (•) is the gradient operator with respect to the model parameters, e.g. the Fréchet derivative.

The corresponding functional derivative of J is given by:

J = r m J m = r m h(Ru(m), m; d obs ) m + ⌦ r u h(Ru(m), m; d obs ), u ↵ W , (2.15)
where h•, •i W denotes the inner product in W,

hu 1 , u 2 i W = Z T 0 ✓ Z Ω u 1 (x,t)u 2 (x,t) dx ◆ dt, (2.16)
and the chain rule is used here together with the definition of u.

The main problem with (2.15) is that explicitly computing the full kernel of the operator r m u can be highly inefficient. The adjoint-state method is a very good way to eliminate r m u m so that J can be computed in more favorable complexity.

In order to achieve this, the state equation L [m] u f = 0 can be linearised, assuming a sufficiently smooth seismic source, leading to:

 r m L [m] m u + L [m] u = 0 (2.17)
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 which shows that u can be eliminated by composition with L [m] on the left. The main idea of the adjointstate method is that a copy of L [m] u can materialise in (2.15) provided that r u h is seen as the adjoint of L [m] applied to some field λ,

[L [m]] † λ = r u h(Ru(m), m; d obs )= R T r d cal h(d cal (m), m; d obs ), (2.18)
with λ naturally called the adjoint field and (2.18) the adjoint-state equation. Then,

J = r m h(Ru(m), m; d obs ) m ⌦ (L [m]) † λ, u ↵ W = r m h(Ru(m), m; d obs ) m ⌦ λ, L [m] u ↵ W
Taking into account (2.17), the first variation of J , now reads

J = r m h(Ru(m), m; d obs ) m + ⌦ λ, [r m L[m] m]u ↵ W (2.19)
which is often much easier to compute than (2.15)

Physics approximation: the acoustic wave

Physical approximation of the dynamical system is an important component in full-waveform inversion. In this thesis, the behaviour of the dynamical system is described by the 2D acoustic wave equation.

Linearised isotropic, inviscid fluid flow, in the physical domain Ω ⇢ R 2 , leads to

@u @t = 1 ⇢ ? (x) rp, @p @t =  ? (x)r • u (2.20)
where the model parameters m 2Dare ⇢ ? : Ω ! R + , the background mass density, and  ? : Ω ! R,t h e background bulk modulus with  ? = ⇢ ? [c ? ] 2 with c ? the sound speed.

The physical state of the dynamical system is w =(u,p) 2W, with u :

Ω ⇥ R + ! R 2 the displacement perturbation and p : Ω ⇥ R + ! R the pressure perturbation.
The behaviour of the system is described by the linear wave operator

L(m)w(x,t)=  @ @t + S(m) w(x,t) with S(m, x)= 0 @ 0 1 ⇢ ? (x) r  ? (x)r • 0 1 A (2.21)
which defines a hyperbolic system of partial differential equations, and S is anti-self-adjoint, i.e. S † = S, with respect to the inner product

hw, wi = 1 2 Z Ω ✓ ⇢ ? u • ũ + 1  ? p p◆ dx (2.22)
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 in which the factor 1/2 is chosen to be consistent with the physics convention for the energy.

The acoustic wave operator can be rewritten in term of the pressure perturbation only

[L p (m) p](x,t)=  1  ? (x,t) @ 2 @t 2 r• ✓ 1 ⇢ ? (x) r ◆ p(x,t) (2.23)
leading to the scalar state equation

[L p (m) p](x,t) f p (x,t)= 1  ? (x) @ 2 p(x,t) @t 2 r• ✓ 1 ⇢ ? (x) rp(x,t) ◆ f p (x,t)=0 (2.24)
where f p : Ω ⇥ R + ! R is the pressure perturbation associated to a seismic source, generally considered as a point source: f p (x,t)=s(t) (x x s ) where x s 2 Ω is the position of the source and s(t) the source time function often modelled by a causal Ricker wavelet

s(t)= ⇥ 1 2⇡ 2 f 2 0 (t t 0 ) 2 ⇤ e ⇡ 2 f 2 0 (t t 0 ) 2 (2.25)
where the time shift t 0 assures the causality, i.e. s(t)=0 , for t<0,a n df 0 is the characteristic frequency of the wavelet.

The state equation ( 2.24) needs to be completed by the initial conditions

@p @t (x, 0) = ṗ0 (x)=0 and p(x, 0) = p 0 (x) = 0; 8x 2 Ω. (2.26)
When the physical domain is limited by an upper free surface Γ 1 , the associated boundary condition is

p| Γ 1 = p(x,t)=0, 8x 2 Γ 1 , 8t 2 [0,T] . (2.27)
When the physical domain is unbounded, the computational domain must be truncated by artificial interfaces Γ 2 together with an appropriate absorbing boundary condition (ABCs). The simplest form of ABC is the zero-order absorbing boundary condition [START_REF] Clayton | Absorbing boundary consitions for acoustic and elastic wave equations[END_REF][START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF][START_REF] Nataf | Absorbing boundary conditions and perfectly matched layers in wave propagation problems, in Direct and Inverse Problems in Wave Propagation and Applications[END_REF] derived from the approximation of the one-way wave equation,

@p @t (x,t)+c ? (x)rp • n =0, 8x 2 Γ 2 , c ? (x)= s  ? (x) ⇢ ? (x) , (2.28)
where n the unit outward normal associated to Γ 2 .

Designing good absorbing boundary conditions is a somewhat difficult problem that has a long history.

The currently most popular solution to this problem is to slightly expand the computational domain through an absorbing perfectly matched layer (PML) [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. The PML shows great superiority in the context of acoustic wave simulation over the classical ABC mentioned above [START_REF] Arbanel | A mathematical analysis of the pml method[END_REF][START_REF] Diaz | A time domain analysis of pml models in acoustics[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Gao | Comparison of artificial absorbing boundaries for acoustic wave equation modelling[END_REF]. In particular, convolutional implementations of complex frequency shifter perfectly matched layer (Roden andGedney , 2000; Festa andVilotte , 2005; Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 [START_REF] Komatitsch | An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave equation in poroelastic media[END_REF] have become widely used in acoustic wave simulation [START_REF] Pasalic | Convolutional perfectly matched layer for isotropic and anisotropic acoustic wave equations[END_REF][START_REF] Gao | Comparison of artificial absorbing boundaries for acoustic wave equation modelling[END_REF]. Here, to simplify the presentation, we adopt the first order absorbing condition (2.28) This leads to the forward modelling problem

8 > > > > > > > > > > < > > > > > > > > > > : 1  ? (x) @ 2 p(x,t) @t 2 r• ✓ 1 ⇢ ? (x) rp(x,t) ◆ f p (x,t) = 0; 8x 2 Ω, 8t 2 [0,T] @p @t (x, 0) = ṗ0 (x)=0 and p(x, 0) = p 0 (x) = 0; 8x 2 Ω p(x,t) = 0; 8x 2 Γ 1 , 8t 2 [0,T] @p @t (x,t)+c ? (x)rp(x,t) • n(x)=0, 8x 2 Γ 2 .
(2.29)

In order to specify the adjoint acoustic wave operator

L † p let consider 8 2W † ; ⌦ L p (m)p, ↵ W = ⌧ 1  ? (x) @ 2 p(x,t) @t 2 r• ✓ 1 ⇢ ? (x) rp(x,t) ◆ , W = Z T 0 Z Ω 1  ? @ 2 p @t 2 r• ✓ 1 ⇢ ? rp ◆ dx dt.
(2.30)

After two integration by part:

Z T 0 Z Ω ✓ 1  ? (x) @ 2 p @t 2 ◆ dx dt = Z T 0 Z Ω ✓ 1  ? @ 2 @t 2 ◆ pdx dt + Z T 0 ✓ 1  ? @p @t ◆ t=T t=0 dx Z Ω ✓ 1  ? (x) @ @t ◆ p t=T t=0
dx.

(2.31)

The initial conditions (2.26) imply that the first two last terms in the right hand side can be eliminated when imposing the terminal conditions for (x,T)=0 and @ @t (x,T) = 0; 8x 2 Ω.

(2.32)

Similarly applying the divergence theorem and the identity

r • ✓ 1 ⇢ ? rp ◆ r• ✓ 1 ⇢ ? pr ◆ = r • ✓ 1 ⇢ ? rp ◆ pr • ✓ 1 ⇢ ? r ◆ , (2.33)
together with the boundary condition (2.27) yields to where ds is a line element and n is the unit outward-pointing normal on the boundaries.

Z T 0 Z Ω r • ✓ 1 ⇢ ? rp ◆ dx dt = Z T 0 Z Ω r • ✓ 1 ⇢ ? r ◆ pdx dt + Z T 0 Z Γ 1 ✓ 1 ⇢ ? rp ◆ • n ds dt + Z T 0 Z Γ 2 ✓ 1 ⇢ ? rp ◆ • n ds dt Z T 0 Z Γ 2 ✓ 1 ⇢ ? pr ◆ • n ds dt,
Making use of the absorbing boundary condition (2.28)

Z T 0 Z Γ 2 ✓ 1 ⇢ ? rp ◆ • n ds dt = Z T 0 Z Γ 2 ✓ 1 ⇢ ? ◆ 1 c @p @t dsdt = Z T 0 Z Γ 2 ✓ 1 ⇢ ? p ◆ 1 c @ @t Z Γ2 1 ⇢ ? 1 c pds t=T t=0 , (2.35)
where the last term can be eliminated when imposing the initial conditions (2.26)o np and the terminal condition (2.32)o n .

Gathering the previous results leads

⌧ 1  ? @ 2 p @t 2 r• ✓ 1 ⇢ ? rp ◆ , W = ⌧ 1  ? @ 2 @t 2 r• ✓ 1 ⇢ ? r ◆ ,p W + Z T 0 Z Γ 2 1 ⇢ ? ✓ 1 c @ @t r • n ◆ pdsdt Z T 0 Z Γ 1 ✓ 1 ⇢ ? rp ◆ • n ds dt.
(2.36)

The two last terms in (2.36) can be eliminated by imposing to the adjoint field the boundary condition

8 < : (x,t) = 0; 8x 2 Γ 1 , 8t 2 [0,T] @ @t cr • n = 0; 8x 2 Γ 2 , 8t 2 [0,T] .
(2.37)

Finally, the equations governing the adjoint wavefield is

8 > > > > > > > > > > < > > > > > > > > > > : 1  ? @ 2 @t 2 r• ✓ 1 ⇢ ? r ◆ = R T r d cal h (d cal = Rp, m; d obs ); 8x 2 Ω, 8t 2 [0,T] @ @t (x,T)=0and (x,T) = 0; 8x 2 Ω | Γ 1 = (x,t) = 0; 8x 2 Γ 1 , 8t 2 [0,T] @ (x,t) @t c ? (x) r (x,t) • n(x)=0, 8x 2 Γ 2 , 8t 2 [0,T] .
(2.38)

The adjoint state is solution of a backward problem where initial conditions (2.26) are replaced by the terminal condition (2.32). Homogeneous free surface boundary condition should be the same for the forward and adjoint wave field, whereas absorbing boundary conditions involving time derivatives needs to be properly time-reversed.

Introducing the change of variable t ! T t, the adjoint problem can be rewritten in the same form as the forward problem with the auxiliary unknown: ˜ = (x,T t) and the time-reverse source R T r d cal h(d cal ; d obs )(T t). The result needs then to be time-reverse at each point x. This is a nice property that enable to solve the forward and the adjoint problem with the solver.

Note however that is not in general the physical wave field run backward in time because of the limited Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 sampling of the receivers, and instead is introduced purely out of computational convenience.

It is interesting to notice that changing the way to measure the distance between observed and modelled data (i.e. modifying the definition of h) only amounts to modify the source term in the adjoint equation.

A simplification: the constant density model

To go further, we will assume the following hypothesis 1. The density is assumed to be constant. More precisely, we assume that

⇢ ? (x)=⇢ 0 , ? (x)=⇢ 0 [c ? (x)] 2 .
In this situation, the only parameter is c ? (x), the velocity law in Ω.

2. The admissible parameter set M is a convex subset of a finite dimensional vectorial space. Practically, we choose M scalar non negative functions m (x), x 2 Ω such that P m m (x)=1 , x 2 Ω and we define (0 <v <v + are some given parameters)

M = ( c ? (x)= M X m=1 v m m (x),v  v m  v + , 8m =1,..,M ) . ( 2
h(d cal , m, d obs )= h(d cal , d obs ) In this situation, (2.19)a m o u n t st o J = ⌦ , [r c ? L p (c ? ) c ? ]p ↵ W . (2.40)
With these simplifications, the pressure p(x,t) is the solution of

8 > > > > > > > > > > < > > > > > > > > > > : 1 [c ? (x)] 2 @ 2 p(x,t) @t 2 ∆p(x,t)=⇢ 0 f p (x,t); 8x 2 Ω, 8t 2 [0,T] @p @t (x, 0) = 0 and p(x, 0) = 0; 8x 2 Ω p(x,t) = 0; 8x 2 Γ 1 , 8t 2 [0,T] @p(x,t) @t + c ? (x) rp(x,t) • n(x)=0, 8x 2 Γ 2 , 8t 2 [0,T] . (2.41) It is possible to show that when f p (x,t) is in W and c ? is in M as given in (2.39) then p is in C 1 ([0,T],L 2 (Ω))\
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

C 0 ([0,T],H 1 (Ω))
. Now, we define the source term for the adjoint equation as

S ad (x,t)) = R T r d cal h (d cal = Rp, d obs )(x,t), (2.42)
or, more explicitely,

S ad (x,t)= X r (x x r )r d cal h (d cal = Rp, d obs )(x r ,t). (2.43)
The adjoint state is solution of

8 > > > > > > > > > > < > > > > > > > > > > : 1 [c ? (x)] 2 @ 2 (x,t) @t 2 ∆ (x,t)=S ad (x,t) 8x 2 Ω, 8t 2 [0,T] @ @t (x,T)=0and (x,T) = 0; 8x 2 Ω (x,t) = 0; 8x 2 Γ 1 , 8t 2 [0,T] @ (x,t) @t c ? (x) r (x,t) • n(x)=0, 8x 2 Γ 2 , 8t 2 [0,T] .
(2.44)

Now, let us define the M -vector v =( v 1 ,v 2 ,...,v M ) T and p(x,t; v), (x,t; v) as p(x, t),a n d (x, t),t h e solutions of (2.41)-( 2.44)w i t hc ? (x)= P m v m m (x) 2M; the functional is now a simple function of vector v and is given by

J(v)= h (Rp(x,t; v); d obs )= h ((p(x 1 ,t; v),p(x 2 ,t; v),...,p(x R ,t; v)); d obs ) ,
and we have, according to (2.40),

J(v)= Z T 0 Z Ω (x,t)( [ r c ? L p (c ? ) c ? ]p(x,t)) dxdt. Since L p (c ? )= 1 [c ? ] 2 @ 2 @t 2 ∆+ Γ 2 ✓ 1 c ? @ @t n •r ◆ , we have [r c ? L p (c ? ) c ? ]= 2 c ? [c ? ] 3 @ 2 @t 2 Γ 2 ✓ c ? [c ? ] 2 @ @t ◆ , or, when c ? (x)= P m v m m (x), [r c ? L p (c ? ) c ? ]= M X m=1 ✓ 2 m v m [c ? ] 3 @ 2 @t 2 + Γ 2 m v m [c ? ] 2 @ @t ◆ , so that J(v)= X m  Z T 0 ✓Z Ω (x,t) 2 m (x) [c ? (x)] 3 @ 2 p(x,t) @t 2 dx + Z Γ 2 (x,t) m (x) [c ? (x)] 2 @p(x,t) @t dx ◆ dt v m ,
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

and finally, the gradient of J at point v is given by

@ J @v m (v)= Z T 0 ✓Z Ω (x,t) 2 m (x) [c ? (x)] 3 @ 2 p(x, t) @t 2 dx + Z Γ 2 (x,t) m (x) [c ? (x)] 2 @p(x, t) @t dx ◆ dt.
(2.45)

Discretisation of the problem

The discretisation of the wave equation in heterogeneous medium is a classical topic and has been investigated

for a long time. Here, we adopt the simplest method but others are possible. The starting point of the problem is to rewrite the wave problems (2.41)a n d( 2.44) in a weak form, i.e.

Find p, in C 1 ([0,T],H 1 (Ω)) \ C 2 ([0,T],L 2 (Ω)) such that for all p t , t in H 1 (Ω), 8 > > > > > > > > > > < > > > > > > > > > > : m ✓ d 2 p dt 2 ,p t ; c ? ◆ + b ✓ dp dt ,p t ; c ? ◆ + k(p, p t )=S p (p t ,t) m ✓ d 2 h dt 2 , t ; c ? ◆ b ✓ d h dt , t ; c ? ◆ + k( h , t )=S ( t ,t) p(•,t = 0) = 0, dp dt (•,t = 0) = 0, (•,t = T )=0, d dt (•,t = T )=0.
(2.46)

In this formulation, appears three bilinear forms and two linear forms :

• The mass bilinear form

m(u, u t ; c ? )= Z Ω u(x) u t (x) [c ⇤ (x)] 2 dx.
• The stiffness bilinear form

k(u, u t )= Z Ω ru(x) • ru t (x) dx.
• The boundary mass bilinear form

b(u, u t ; c ? )= Z Γ 2 u | Γ 2 (x) u t | Γ 2 (x) c ⇤ (x) dx.
• The two linear forms associated to the second terms

S p (p t ,t)=⇢ 0 Z Ω f p (x,t) p t (x) dx,S ( t ,t)= Z Ω S ad (x,t) t (x) dx.
This weak formulation is obtained by multiplying each wave equation by a test function then integrating the result over Ω and finally by using the Green's Theorem as well as the boundary conditions.

Remark: When the sources are ponctual, this formulation makes no sense since H 1 (Ω) contains functions which are not continuous on Ω. This problem can be avoided by regularizing the Dirac mass, i.e. replacing (x x o ) by some smooth positive function r (x x o ) supported by a small neighborrow of x 0 such that

R Ω r (x x o )=1.
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Semi discretisation in space

We assume that Ω is a square in the plane, Γ 1 is its top horizontal line and Γ 2 is composed of the three other sides of the square boundary. We mesh Ω with a regular grid of step h. The set of the nodes of this grid is denoted by Ñ and the associated node is x n ; in order to distinguish between interior nodes and nodes located on the boundaries we set

Ñ = N N Γ 1 , N = N 0 N Γ 2 .
The grid is decomposed in squared cells of side h; the set of all cells is denoted by C: Ω= L C2C C;a t each node n 2N we associate C[n] the set of cells whose one vertex is node n. This set is composed of four cells for n 2N 0 , two or one cell(s) for n 2N Γ 2 ; we define N h n (x), n 2N the only continuous function whose restriction on each cell is linear and such that

N h n (x n )=1, N h n (x n 0 )=0,n6 = n 0 2N 0 ; its easy to see that the support of N h n (x) is [ C2C[n] C.
The semi discretization of problems (2.41)a n d( 2.44) consists in looking for approximate solutions in V h (Ω) ⇢ H 1 (Ω), the vector space spanned by the family of functions

N h n (x),n2N p(x,t) ' p h (x,t)= X n2N p n (t)N h n (x), (x,t) ' h (x,t)= X n2N n (t)N h n (x);
p h and h are defined as the solution of

Find p h , h in C 2 ([0,T],V h (Ω)) such that for all N h n , n 2N 8 > > > > > > > > > > < > > > > > > > > > > : m h ✓ d 2 p h dt 2 ,N h n ; c ? ◆ + b h ✓ dp h dt ,N h n ; c ? ◆ + k(p h ,N h n )=S p (N h n ,t) m h ✓ d 2 h dt 2 ,N h n ; c ? ◆ b h ✓ d h dt ,N h n ; c ? ◆ + k( h ,N h n )=S ad (N h n ,t) p h (•,t = 0) = 0, dp h dt (•,t = 0) = 0, h (•,t = T )=0, d h dt (•,t = T )=0.
To obtain this semi-discrete problem, we made the substitutions p ! p h and ! h and restricted the tests functions to element of V h (Ω). The linear and bilinear forms have been approximated (the substrict h )b y using a quadrature rule on each cell; more precisely, we assume that the restriction of c ? is constant on each cell

c ? (x)= X C2C v c 1 C (x).
We split the integrals according to

Z Ω u(x) N h n (x) [c ? (x)] 2 dx = X C2C(n) Z C u(x) N h n (x) [c ? (x)] 2 dx = X C2C(n) 1 v 2 c Z C u(x) N h n (x) dx. Z Γ 2 u | Γ 2 (x) N h n (x) [c ? (x)] dx = X C2C(n) Z Γ 2 \@C u | Γ 2 (x) N h n (x) [c ? (x)] dx = X C2C(n) 1 v c Z Γ 2 \@C u | Γ 2 (x) N h n (x) dx,
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 then we use the quadature rule (V(C) is the set of the four vertexes of cell C and V(@C) is the set of the two

vertexes of Σ C = @C \ Γ 2 when Σ C is not empty) Z C u(x) N h n (x) dx ' h 2 4 X xm2V(C) u(x m ) N h n (x m )= h 2 4 u(x n ), Z Γ 2 \@C u | Γ 2 (x) N h n (x) dx ' h 2 X xm2V(@C) u | Γ 2 (x m ) N h n (x m )= h 2 u | Γ 2 (x n ).
With this process, we obtain

m h (u h ,N h n ; c ? )= ✓ h 2 4 X C2C(n) 1 v 2 c ◆ u h (x n ),b h (u h ,N h n ; c ? )= ✓ h 2 X C2C(n) 1 v c ◆ [u h ] | Γ 2 (x n ).
This approximation method is known as Mass Lumping; it is essential to obtain in fine a time explicit scheme. Now, we define the two diagonal mass matrices ( m n is the Kronecker symbol)

M(v) n,m = m n ( h 2 4 X C2C(n) 1 v 2 c ), B(v) n,m = 8 > < > : m n ✓ h 2 P C2C(n) 1 vc ◆ , if x n 2N Γ 2 0 if x n 2N 0 (2.47)
and the stiffness matrix

K n,m = Z Ω rN h n (x) • rN h m (x) dx,
which corresponds for n or m in N 0 (interior nodes) to a nine points stencil for h 2 times the Laplacian. Note that the coefficients of matrix K are bounded by a pure constant, we deduce that there exists a pure non negative number  such that

(KU, U )   (U, U ).
(2.48)

These definitions being made, the semi discrete problem can be put in the following matricial form

Find P (t)=(p n (t)) n2N and Λ(t)=( n (t)) n2N in C 2 ([0,T], R #N ) such that 8 > > > > > > > > > < > > > > > > > > > : M(v) d 2 P (t) dt 2 + B(v) dP (t) dt + KP (t)= Sp (t) M(v) d 2 Λ(t) dt 2 B(v) dΛ(t) dt + KΛ(t)= Sad (t) P (t = 0) = 0, dP dt (t = 0) = 0, Λ(t = T )=0, dΛ dt (t = T )=0.

Discretisation in time

We pick a time step ∆t; our aim is to compute approximations of vectors P (t) and Λ(t) at discrete instant

k∆t P k ' P (k∆t), Λ k ' Λ(k∆t),k =0,...,K, K∆t = T.
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

For that, we introduce the finite difference operators

[@ 2 ∆t 2 U ] k = U k+1 2U k + U k 1 ∆t 2 , [@ 2∆t U ] k = U k+1 U k 1 2∆t ,
and define the totally discretised problem Find P k and Λ k in R #N ,k=0,...,K such that for all k>0,

8 > > > > > > < > > > > > > : M(v)[@ ∆t 2 P ] k + B(v)[@ 2∆t P ] k + KP k = Sk p = Sp (t = k∆t) M(v)[@ ∆t 2 Λ] k B(v)[@ 2∆t Λ] k + KΛ k = Sad (t = k∆t) P 0 =0,P 1 =0, Λ K =0, Λ K 1 =0.
Note that computing P k+1 (resp. Λ k 1 ) from P k ,P k 1 (resp. Λ k+1 , Λ k ) involves the inversion of matrix v); this matrix being diagonal (thanks to Mass Lumping) the scheme is explicit.

M(v)+ ∆t 2 B(

Stability and CFL condition.

The stability of this scheme relies on a conservation of a pseudo energy. For simplicity, we drop the source term and consider instead initial conditions. The model problem is now Find U k in R #N ,k=0,...,K such that for all k>0, 1) .

8 > > > < > > > : M(v)[@ ∆t 2 U ] k + B(v)[@ 2∆t U ] k + KU k =0 1 2 (U 1 + U 0 )=u (0) , U 1 + U 0 ∆t = u (
We take the scalar product of the equation in

U k by [@ 2∆t U ] k ,w eo b t a i n ⇣ M(v)[@ ∆t 2 U ] k , [@ 2∆t U ] k ⌘ + ⇣ KU k , [@ 2∆t U ] k ⌘ + ⇣ B(v)[@ 2∆t U ] k , [@ 2∆t U ] k ⌘ =0. (2.49)
Let us define

V k+ 1 2 =[@ ∆t U ] k+ 1 2 = U k+1 U k ∆t , [ ∆t U ] k+ 1 2 = U k+1 + U k 2 , it is straightforward to show that ∆t ⇣ M(v)[@ ∆t 2 U ] k , [@ 2∆t U ] k ⌘ = 1 2 ⇣ M(v) ⇣ [@ ∆t U ] k+ 1 2 [@ ∆t U ] k 1 2 ⌘ , ⇣ [@ ∆t U ] k+ 1 2 +[@ ∆t U ] k 1 2 ⌘⌘ = 1 2 ⇣ M(v)V k+ 1 2 ,V k+ 1 2 ⌘ 1 2 ⇣ M(v)V k 1 2 ,V k 1 2 ⌘ ,
the last line being obtained thanks to the symmetry of matrix. M(v).I nt h es a m ew a y ,w eh a v e ∆t

⇣ KU k , [@ 2∆t U ] k ⌘ = 1 2 ⇣ KU k ,U k+1 ⌘ 1 2 ⇣ KU k ,U k 1 ⌘ .
Institut

U k a = ⇣ B(v)[@ 2∆t U ] k , [@ 2∆t U ] k ⌘ 0 E k+ 1 2 = 1 2 ⇣ M(v)V k 1 2 ,V k 1 2 ⌘ + 1 2 ⇣ KU k ,U k+1 ⌘ ,
we get from (2.49),

E k+ 1 2 E k 1 2 +∆t U k a =0,k 1,
and by a simple summation

E k+ 1 2 + k X q=1 ∆t U q a = E 1 2 .
(2.50)

At this point, the problem is that the potential energy 1 2 KU k ,U k+1 might be a negative term since vector U is not taken at the same instant. However, we can rewrite

U k+1 =[ ∆t U ] k+ 1 2 + ∆t 2 [@ ∆t U ] k+ 1 2 ,U k =[ ∆t U ] k+ 1 2 ∆t 2 [@ ∆t U ] k+ 1 2 ,
and so (due to the symmetry of Matrix K),

⇣ KU k ,U k+1 ⌘ = ⇣ K[ ∆t U ] k+ 1 2 , [ ∆t U ] k+ 1 2 ⌘ ∆t 2 4 ⇣ K[@ ∆t U ] k+ 1 2 , [@ ∆t U ] k+ 1 2 ⌘ .
The energy term can be rewritten

E k+ 1 2 = 1 2 ⇣ M(v)V k 1 2 ,V k 1 2 ⌘ + 1 2 ⇣ K ∆t U k+ 1 2 , ∆t U k+ 1 2 ⌘ , M(v)=M(v) ∆t 2 4 K.
Our claim is that if this matrix is definite positive, i.e. ∆t 2 4 K < M(v), then our explicit scheme is stable. We can get a crude approximation of this condition using (2.48):

∆t 2 4 (KU, U )   ∆t 2 v 2 + 4h 2 h 2 v 2 + (U, U )   4 ∆t 2 v 2 + h 2 (M(v)U, U ),
with v + the maximum of the components of v. Finally, there exists a pure constant ↵ 2   4 such that when

C = ↵ 2 ∆t 2 v 2 + h 2 < 1( CFL condition), (2.51)
then E k+ 1 2 is always positive and consequently (see (2.50)a n du s eU q a > 0)

0 E k+ 1 2 E 1 2 .
It is easy to see that this ap r i o r iestimation implies that for each node n

⇣ U k+1 n U k n ⌘ 2  4v 2 + ∆t 2 h 2 E 1 2  4 ↵ 2 E 1 2 ,
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

and U k n can only grow linearly with k : the scheme is stable.

The estimates for the discretized problem with a source term is more involved and we fail to derive a simple estimate involving some norm of the source term in this case.

Accuracy: the dispersion effect

Analysis of the precision of such schemes are usually done for infinite grids (i.e. Ω = R 2 ) and constant velocity. The scheme is simply

h 2 c 2 [@ 2 ∆t 2 U ] k + KU k =0.
The scheme being invariant by time and spatial translations, there exists particular solutions of the form (discrete plane wave solution)

U k = U (0) ,e i!k∆t , with 1 h 2 KU ? = 4sin 2 (! ∆t 2 ) c 2 ∆t 2 U ? .
U ? is an eigenvector of matrix K. Let x n = ihx + jhẑ, we look for eigenvectors in the form

U ? n = e i kh (!)•xn = e i(k h x (!)ih+k y h (!)jh) .
Plugging this expression, we obtain the dispersion relation of the scheme

4sin 2 (! ∆t 2 ) c 2 ∆t 2 = 1 h 2 K(k h x (!)h, k h y (!)h); (2.52)
K is the symbol of matrix K. We skip its exact definition, but only retain that due to both the centered nature and the consistency of the scheme, we have

1 h 2 K(k h x (!)h, k h y (!) 2 )=[k h x (!) 2 + k h y (!))] ⇣ 1+O((k h x (!)h) 2 )+O((k h y (!)h) 2 ) ⌘ .
Now the Taylor development of sin x x gives immediately

! 2 c 2 =[k h x (!) 2 + k h y (!)) 2 ] ⇣ 1+O((k h x (!)h) 2 )+O((k h y (!)h) 2 )+O((!∆t) 2 ) ⌘ ,
and the relation of dispersion

! 2 c 2 (k 2
x +k 2 y ) =1is asymptotically recovered: the scheme is consistent. However, there is a second order error in space and time and this error is the cause of what is called the Pollution effect; as the wave propagates in the medium, the errors accumulate and this, all the more since the duration of the wave increases. It is usually mentioned that 10 points per wavelengh is enough to get a "precise" result, let us say ✏ =0.1% for some h. As the matter of fact, this result has been obtained in the seventies, when the computers only allows for a propagation over roughly 6 wavelengh. As the errors are second order in time and space, a propagation along 4 ⇥ 6 wavelengh would induce an error of 2✏. To get the same error we need to divide the space step h (as well as the time step for a constant CFL) by p 2.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

The analysis of the error of the numerical schemes for wave equations in heterogeneous medium has not been studied a lot. Usually, a thumb rule consists of choosing F c , a cut frequency of the seismic source and take n points per wavelenghts at this frequency and for the slowest velocity in the medium. The formula is

h = v F c n . (2.53)
The choice of n depends, as we said above, on the size of the domain of propagation and on the required accuracy.

Note that for our inverse problem, we do not have a precise knowledge of the values of the velocities involved. The simplest way to control the dispersion effect is to find their values into a given and not too large interval.

The discrete functional and its gradient

At this point we have a process to compute an approximation of p h (x,t) at nodes x n of grid and time k∆t.

Our discrete observed data set D h,∆t are seismograms recorded at n r receivers located on the grid

d obs : {d obs (x nr ,k∆t),n r 2N r ⇢N,k 2 K r ⇢{1,...,K} .} , (2.54)
The discrete predicted observables are defined as

d cal (v): n P (v) k nr = R h,∆t P (v); n r 2N r ,k2 K r o , (2.55)
where

R h,∆t :[R #N ] K !D h,∆t is the discrete extraction operator. R h,∆t : P =(P k n ) k=1,..,K,n2N ! (P k n ) k2Kr,n2Nr (2.56)
The discrete functional is now

J (v)= hh,∆t (d cal = R h,∆t P (v); d obs ) .
To compute the first variation of the discrete functional J , we move parameter from v to v + v,thesolution

P = P (v)
of the direct problem move to P + P and to the first order we have

J (v)= X k2Kr X nr2Nr r cal h(x nr ,k∆t) P k nr = K X k=1 ⇣ R T h,∆t r d cal h(R h,∆t P (v); d obs ), P k ⌘ . (2.57)
Now, to the first order we have

8 > > < > > : M(v)[@ ∆t 2 P ] k + B(v)[@ 2∆t P ] k + K P k = M(v)[@ ∆t 2 P ] k B(v)[@ 2∆t P ] k P 1 = P 0 =0.
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

We take the scalar product of each of these equations by Λ k and sum over k =1,...,K. Thanks to time discrete integrations by parts and using the null initial and final conditions, we get

Λ K =Λ K 1 = P 1 = P 0 =0 8 > > > > > > > > > > > > > > < > > > > > > > > > > > > > > : K X k=1 ⇣ M(v)[@ ∆t 2 P ] k , Λ k ⌘ = K X k=1 ⇣ M(v)[@ ∆t 2 Λ] k , P k ⌘ K X k=1 ⇣ B(v)[@ 2∆t P ] k , Λ k ⌘ = K X k=1 ⇣ B(v)[@ 2∆t Λ] k , P k ⌘ K X k=1 ⇣ K P k , Λ k ⌘ = K X k=1 ⇣ KΛ k , P k ⌘ .
Adding the three equations and using the equations satisfied by Λ k and P k , we get

K X k=1 ⇣ M(v)[@ ∆t 2 P ] k , Λ k ⌘ K X k=1 ⇣ M(v)[@ ∆t 2 P ] k , Λ k ⌘ = K X k=1 ⇣ S k ad , P k ⌘ . (2.58)
Returning back to (2.57), we get

J h,∆t = K X k=1 ⇣ M(v)[@ ∆t 2 P ] k , Λ k ⌘ + ⇣ B(v)[@ 2∆t P ] k , Λ k ⌘ , (2.59)
as soon as (S k ad ) n = r cal h(x nr ,k∆t) for k 2 K, n 2N r , and 0 elsewhere.

From this we deduce the expression of the gradient of J h,∆t ,

@J h,∆t @v c = K X k=1 ✓ @M(v) @v c [@ ∆t 2 P ] k , Λ k ◆ + ✓ @B(v) @v c [@ 2∆t P ] k , Λ k ◆ . ( 2

.60)

What we have done here is to compute the exact gradient of the discrete functional. It is nowadays well known that discretizing directly the gradient of the continuous functional given in (2.45) by some other consistent scheme leads to problems of coherency when using a solver.

From the very expressions of both diagonal matrices M and B (see (2.47)), it is seen that the two scalar products in (2.60) involve only a finite number of terms at each time step k∆t. The remaining problem is that P k is computed forward in time while Λ k is computed backward in time. Storing all the P k 's is memory consuming (solving the scheme only requires to store two successive instants). To overcome this the P k 's are first computed for k =1,..,K, then we restart from the two last instants P K and P K 1 and recompute backward in time the P k and Λ k , k = K, K 1..., 1; in order to avoid instability problems, it is nevertheless necessary to store all the P k n at node n 2N Γ 2 (i.e. on the absorbing boundary) since the boundary condition is only absorbant in the forward direction.

Institut de Physique du Globe de Paris

Miao YU, Ph.D. Thesis, 2021

Algorithms and Computational framework

We have now all the ingredients for designing the algorithm for our Full Waveform Inversion method. A simple FWI workflow can be summarised by the flow chart of Figure 2.1.

For a given seismic source and an initial model: In PySIT, the 2D acoustic wave forward and adjoint problems are solved using standard second-order finite differences in space and time on a uniform mesh in space together with a second-order explicit integration scheme in time with constant time step. Complex shifted frequency (CFS) PMLs are implemented following [START_REF] Sim | Nonreflecting boundary conditions for the time dependent wave equation[END_REF][START_REF] Grote | Local nonreflecting boundary condition for time-dependent multiple scattering[END_REF]. Thus, there is slight differences from what we have presented in the previous subsections (first order absorbing conditions instead of PML, a different molecule for the discrete Laplacian,..) but all the ideas remain the same: we have a CFL condition, see (2.51), and the dispersion of the waves must be controlled using a large enough number of points per shortest wavelenghts (i.e. for the maximum frequency of the signal and the minimum velocity).

Least-squares misfit function and the cycle-skipping problem

Since the introduction of full-waveform inversion, the most commonly used objective function measures the L 2 norm difference between the synthetic observables and the observations.

As previously seen, in the data space D, observations d obs are common shot-gather of seismograms, often called seismic traces, recorded at n r receivers for a source s:

d s obs : {d s obs (x r ,t),r=1,...,n r } ,t2 [0,T] (2.61)
where x r denotes the spatial position vector of the receiver r

The predicted observables are defined accordingly for a physical realisation satisfying the state equation L[m]u f s , together with the appropriate initial and boundary conditions

d s cal (m):{d s cal (x r ,t; m),r=1,...,n r } = Ru s (x,t; m); t 2 [0,T] (2.62)
where R : W!Ddefined in (2.6) and maps the physical realisation onto the receiver positions x r .

The least-squares based objective function is defined as

J (m)= h(d cal (m); d obs )= 1 2 ns X s=1 Z T 0 |d s cal (m) d obs | 2 dt = 1 2 ns X s=1 Z T 0 |Ru s (m) d obs | 2 dt (2.63)
where n s is the number of seismic source used in the seismic survey. The monotonicity of the optimal transport with respect to time-and-space shifts and dilations, together with improved robustness with respect to noise, has attracted the interest of the image processing and machine learnings. Today, there is a growing interest of the FWI community to extend the optimal transport in the context of seismic oscillatory signals for developing alternative continuous and differentiable objective functions. In the next section, a brief overview of the theoretical optimal transport framework of the optimal transport and of some of the current limitations in the context of FWI is provided.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Chapter 3

A brief overview of the Optimal Transport.

This chapter provides a very brief overview of optimal transport introducing the main notions needed in this thesis, comprehensive presentation of optimal Transport can be found for example in [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF] and [START_REF] Peyré | Computational optimal transport[END_REF], and the references therein.

Optimal transport (OT), is an old problem originally formulated by Monge (1781)asfindingtheoptimal way of moving soil from two places through a volume preserving transport in order to make embarkments.

Assuming constant density and no mass creation or loss during the transport, the problem is to find a transport map that minimises the mechanical work , i.e. the mass µ(x) times the moving distance. The later being defined commonly as a L p transport distance with p 1.

Monge's problem

Let consider two probability measures: a "source" µ and a "target" ⌫ that have a compact support X 2 R d and Y 2 R d , respectively, and the density defined with respect to the Lebesgue measure.

Mathematically, the problem is to find a transport map T from X to Y such that MP(µ, ⌫):

T = arg inf T 2M Z X kx T (x)k p dµ(x), (3.1)
where the displacement cost, i.e. the ground cost, is here defined as the L p distance, while value function in (3.1) is called the transport cost.

Let denote M := {T : X ! Y |T # µ = ⌫} the set of all measure preserving maps between µ and ⌫ with T # µ the pushforward of the measure µ, defined as:

[T # µ](B)=µ(T 1 (B)), For any measurable subset B 2 Y
The Monge formulation (3.1) defines a non-convex optimisation problem with non linear constraints. However when, for example, µ and ⌫ have densities, then the optimal transport map T ? exist a give a natural interpolation between the two measures. In particular for the L p distance, the map T t (x)=(1 t)x + tT ? (x) describes the path of particle x and furthermore the measure of µ pushed forward by T t is the geodesic, i.e.

shortest path for the OT induced metric, between µ and ⌫.

If the map is smooth and the measures absolutely continuous, the constraint can be written as an equation for the Jacobian of T :

µ(x)=⌫ (T (x)) |det rT (x)| (3.2)
3.2 Monge Ampère and semi-discrete formulation Since Brenier's theorem [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] it is known that the quadratic cost OT problem, i.e. p =2, has nice properties: under very general assumptions on the densities µ and ⌫, there exists a unique transport map, which is characterised as the gradient of a convex potential '(x): T (x)=r'(x).

Using this result in (3.2) yields the Monge-Ampère form of Monge problem :

µ(x)=⌫ (r'(x)) det r 2 '(x) r'(X) ⇢ Y (3.3)
where ' is a convex function.

When ⌫(x): =µ(x ⌧ ), the convex potential is given analytically by '(x)= 1 2 kx 2 k ⌧ • x + C (this is unique up to a constant). Plugging this in (3.1) the transport cost is k⌧ k 2 , i.e quadratic with respect to the shift ⌧ . We will illustrate this property in chapter 5.

Several numerical methods have been proposed to discretise and solve this problem:

• A first class of methods [START_REF] Benamou | Numerical solution of the Optimal Transportation problem using the Monge-Ampère equation[END_REF][START_REF] Benamou | Minimal convex extensions and finite difference discretisation of the quadratic Monge-Kantorovich problem[END_REF] is based on the discretisation of X using cartesian grids and optimised finite difference stencils. The measure µ is approximated as a discrete measure prescribed on the grid but ⌫ requires a given analytic function or an oracle based on interpolation of ⌫, evaluated at r'(x) in Y . This is the approach followed in particular by Engquist and collaborators [START_REF] Engquist | Application of the wasserstein metric to seismic signals[END_REF][START_REF] Engquist | Optimal transport for seismic full waveform inversion[END_REF]Yang et al., 2018) in the context of full-waveform inversion.

• A second class of methods [START_REF] Mérigot | Xhapter 2 -optimal transport: discretization and algorithms[END_REF][START_REF] Lévy | Notions of optimal transport theory and how to implement them on a computer[END_REF] relies on a semidiscrete approach, in which (3.3) is solved from a dual point of view. The measure ⌫ is assumed to be an arbitrary empirical measure, i.e. a weighted point cloud, and a continuous oracle is still needed but only for µ.

Both approaches lead to a non-linear set of equations that can be solved with a quasi-newton solver in O(N log N ) operations, N being the number of points of the discrete measure, but with a problem-dependent constant that may be large.

Kantorovich Formulation

When both probability measures µ and ⌫ have discrete support, e.g. as in numerical approximations, rearranging the mass using a one-to-one map T may become impossible, see Figure 3 Kantorovich (1942) suggested to optimise instead a relaxed functional with respect to the transport plans:

KP(µ, ⌫)= inf 2Π(µ,⌫) Z X⇥Y kx yk p d (x, y). (3.4)
The set of admissible transport plans

Π(µ, ⌫)={ 2P(X ⇥ Y ) | P X # = µ, P Y # = ⌫} (3.5)
is the set of probability measures on the product space X ⇥ Y , also called "couplings", with prescribed marginal laws µ and ⌫,a n dP X : X ⇥ Y ! X and P Y : X ⇥ Y ! Y stand for the projections onto each space.

Transp ort plans still transp ort the µ-mass repartition to the ⌫-mass repartition but the mass at location

x can be split and sent to different locations y, see blueFigure 3.1 for illustration. The cost is again the total work, i.e. mass times ground cost.

For the L p ground cost, p>1, the Monge and the Kantorovich formulations are equivalent. They both define the same distance on the set of probability measures: the Wasserstein distance denoted W p p (µ, ⌫),

W p p (µ, ⌫)= inf 2Π(µ,⌫) Z X⇥Y kx yk p d (x, y)= inf T 2M Z X kx T (x)k p µ(x) dx (3.6)
and, when µ has continuous densities, the support of any optimal transport plan ? is contained on the graph of a function T ? . In particular this implies ? (A, B)=µ ({x : x 2 A, T ? (x) 2 B}) and furthermore that the optimal plane defines a mapping between µ and ⌫.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

The dual problem associated to (3.6) is formulated as

W p p (µ, ⌫)= sup 2C(X), 2C(Y ) s.t. (x)+ (y) kx yk p , 8(x,y)2X⇥Y. Z X (x) dµ(x)+ Z Y (y) d⌫(y). (3.7)
When freezing ⌫ and considering the Wasserstein distance as a function of µ, the Wasserstein distance is the supremum of affine functions µ parameterised by ( , ), and as such is convex with respect to µ. The dual formulation (3.7) leads to a convex, linear optimisation problem with linear constraints.

Under the primal form (3.6), the optimal transport is defined for probability measures, while the dual form (3.7) admits a notion of first variation over the space of general mesures, provided they have the same total mass, see Santambrogio (2015) section 7.2 for details.

Discrete Kantorovich problem

When X and Y are discrete, i.e. X = {x i ,i =1 ,...,N},Y = {y j ,j =1 ,...,N} where the x i and y j are in R d , µ and ⌫ can be written as the sums of Dirac masses:

µ(x)= N X i=1 µ i x i and ⌫(y)= N X j=1 ⌫ j y j .
The weights µ i and ⌫ j have to satisfy the positivity and the total mass balance constraints:

8 > > < > > : µ i > 0,⌫ j > 0, X i µ i =1, X j ⌫ j =1. (3.8)
The set of admissible plans becomes:

Π(µ, ⌫)= 8 < : 2S N,N : N X j=1 ij = µ i , N X i=1 ij = ⌫ j , 8i, j 9
= ; with S N,N = 8 < : ij 2 R N ⇥N + : N X i=1 N X j=1 ij =1 9 =
;

The marginal constraints are given by summing over lines and columns the matrix .

Finally the optimisation problem becomes:

min 2Π(µ,⌫) N X i=1 N X j=1 c ij ij ,c ij = |x i y j | p (3.9)
where c ij is the ground cost between x i and y j .

Institut The dual formulation associated to (3.7) yields

OT(µ, ⌫) = max u2R N ,v2R N u i +v j c ij X i u i µ i + X j v j ⌫ j (3.11)
In the context of this thesis, the discretisation (x i ,y j ) is assumed to be static, and defined by the receiver positions and by the time sampling of the acquisition, and (µ, ⌫) are vectors in R N with coefficients (µ i )s

and (⌫ j )s.

The convexity of the map µ 2 R N 7 ! OT(µ, ⌫) again follows as the supremum of affine functions, and Danskin's theorem provides the gradient with respect to µ

@ @µ i OT(µ, ⌫)=u ? i , @ @⌫ j OT(µ, ⌫)=v ? j , 8 i, j (3.12)
where (u ? ,v ? ) are the maximisers of (3.11).

The W 1 distance

The case p =1is special. The Kantorovitch metric is well defined but the Monge problem has no unique solution [START_REF] Kantorovich | On a functional space and certain extremum problems[END_REF][START_REF] Villani | Optimal Transport, Grundlehren der mathematischen Wissenschaffen[END_REF].

Assuming X = Y , the dual formulation, a particular instance of the Kantorowitch-Rubinstein theorem, e.g. Santambrogio (2015) section 3.2.1, yields

W 1 (µ, ⌫)=s u p Lip( )1 Z X (x) d(µ ⌫). (3.13)
where Lip( ) denotes the minimal Lipschitz constant for .

The 1-Lipschitz functions , for the ground cost associated to the `1 distance on R d , satisfy

: X ! R, | (x) (y)||x y| , 8(x, y) 2 X ⇥ X with |x y| = d X i=1 |x i y i | , 8(x, y) 2 R d
The Kantorovitch problem can be generalised when the total mass between µ and ⌫ is not conserved, leading Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

to the generalised W1 distance

W 1 (µ, ⌫)= s u p Lip( )1s.t. k k 1 C Z X (x) d(µ ⌫). (3.14)
which is a particular instance of the KR norm [START_REF] Bogachev | Measure thory[END_REF] defined in the space of Radon measures.

This can also be seen as a generalisation of the L 1 norm (Lel lmann et al., 2014).

The generalised W1 metric, also referred to as the KR distance, was first introduced in the context of fullwaveform inversion by Métivier and collaborators (Métivier et al., 2016a,b;Métivier et al., 2016).

One advantage of the KR norm in the context of FWI is that there is no need to normalise the data to be positive and mass balanced. The KR norm has no direct connection with optimal transport when µ and ⌫ are no more required to be probability measures.

The problem (3.14) leads to a convex optimisation problem under linear constraints that can be solved efficiently using proximal splitting methods [START_REF] Combettes | Proximal splitting methods in signal processing, in Fixed-point algorithms for inverse problems in science and engineering[END_REF], reducing significantly the computational cost compared to the 2-Wasserstein metric.

Optimal Transport distances in the context of full-waveform inversion

Classical optimal transport defines distances W p p , p>1 on the set of probability measures, assuming µ>0,⌫>0, and

Z X dµ = Z Y d⌫.
where µ and ⌫ have for support X 2 R d and Y 2 R d , respectively.

In the context of full-waveform inversion, the observables d ob and d cal are oscillatory signals in time, recorded at receiver positions over a time window. As such the positivity assumption breaks down and the total mass conservation needs to be carefully checked in practice.

Apart from data normalisation, there exists a significant amount of work from the mathematics community to generalise optimal transport to what is called the unbalanced optimal transport, either by penalising the marginal constraints [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transport problem[END_REF]Piccoli andRossi, 2014, 2016;[START_REF] Lombardi | Eulerian models and algorithms for unbalanced optimal transport[END_REF][START_REF] Gangbo | Unormalised optimal transport[END_REF] or by rigorously defining an "unbalanced" distance on the set of positive Radon measures (Chizat et al., 2018b;[START_REF] Chizar | Unbalanced optimal transport: Dynamic and kantorovich formulations[END_REF]Kondratyev et al., 2016b;[START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF]. The later will be considered in connexion with the entropic approximation of optimal transport and the Sinkhorn divergence formulation.

Optimal transport and oscillatory signals

In the mathematics community [START_REF] Ambrosio | Gradient flow of the chapman-rubinstein-schatzman model for signed vortices[END_REF][START_REF] Mainini | A description of transport cost for signed measures[END_REF] different approach have been explored

to extend optimal transport to signed measure. One approach leads to split the predicted d cal and the This amounts to compute the optimal transport distance between the positive and negative parts of d obs d cal leading to the objective function h(d cal ,d obs = W 2 2 ( dcal , dobs ). This type of transformation suffers however from a convexity loss with respect to translation and from a loss of sensibility as translation increases [START_REF] Métivier | Optimal transport for mitigating cycle skipping in full-waveform inversion: A graph-space transform approach[END_REF].

Another approach [START_REF] Engquist | Application of the wasserstein metric to seismic signals[END_REF][START_REF] Engquist | Optimal transport for seismic full waveform inversion[END_REF] This strategy assures convexity, but provides no guarantee that the positive part, respectively the negative part, of the predicted signal has the same total mass than the positive part, respectively the negative part, of the observed signal. Moreover such a transformation is not differentiable, which is problematic when computing the Fréchet derivatives of the objective function in the context of local gradient-based optimisation methods.

Signal transformation into probability distribution

Another approach to achieve data positivity and mass balance is to directly transform the seismic data into probability densities by linear or nonlinear scaling functions [START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport[END_REF][START_REF] Yang | Optimal transport for seismic inverse problems[END_REF]Engquist , 2017, 2018;[START_REF] Engquist | Optimal transport based seismic inversion: Beyond cycle skipping[END_REF]. The objective function is now defined as h(d cal ,d obs )=W 2 2 ( dcal , dobs ) A first type of transformation is the affine transformation

dcal = d cal + b hd cal + bi , dobs = d obs + b hd obs + bi ,b > 0 (3.17)
where b is chosen such that dcal and dobs are both positive.

This kind of transformation is straightforward and ensures the total mass conservation. The convexity of the objective with respect to translation is however affected under such transformation as the constant b leads to the creation of artificial mass in the optimal transport.

A This is a particular instance of the KR norm defined on the space of signed measures [START_REF] Bogachev | Measure thory[END_REF]Lel lmann et al., 2014) that is lifted by adding a simple bound constraint on the dual potential to deal with unbalanced mass transport. As such there is no need to normalised the data to be positive and mass balanced.

The resulting objective function is differentiable and allows for simultaneous comparison of 2D and 3D shot gathers following the numerical strategy of Métivier et al. (2016a); Métivier et al. (2016). This strategy was shown in [START_REF] Métivier | Optimal transport for mitigating cycle skipping in full-waveform inversion: A graph-space transform approach[END_REF] to be equivalent to the one proposed by [START_REF] Mainini | A description of transport cost for signed measures[END_REF] in the case of the 1-Wasserstein distance.

However the KR norm has no direct connection with optimal transport once we no longer require d obs and d cal to be probability measures [START_REF] Vershik | Long history of the monge-kantorovich transportation problem[END_REF], and the convexity is not guaranteed with respect to large translation. This research topic has been very active since for two reasons: the method comes with a simple computational method called the Sinkhorn algorithm for which there are now countless variants and acceleration techniques; the class of problems it can applied to is very flexible and includes in particular the recent extension of OT to positive radon measure called "unbalanced" optimal transport (Chizat et al., 2018a;[START_REF] Chizar | Unbalanced optimal transport: Dynamic and kantorovich formulations[END_REF][START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF]. See section 4.2. We discuss its use in connection with positive/negative mass splitting in chapter 5.

Two problems arise when using entropic OT as a proxy for the classical OT distance or unbalanced OT distance. First, the numerical stability and cost of Sinkhorn algorithm deteriorates when " goes to 0.

Second, the entropic OT cost is not a mathematical distance on the set of probability measure anymore. A simple variant called Sinkhorn divergence [START_REF] Genevay | Learning generative models with sinkhorn divergences[END_REF][START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF]

Entropic regularisation and Sinkhorn algorithm

Entropic regularisation

This section is presented in the discrete Kantorovich setting, section 3.2.2. The entropic regularisation of (3.4)-(3.5) has been introduced in the OT context by [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], the Shannon entropy is defined as Ent( )

def. = X i,j ij (log( ij ) 1), 8(i, j) 2 [1,n] ⇥ [1,m] (4.1)
The Entropic regularised problem becomes

OT " (µ, ⌫)= min 2Π(µ,⌫) h , ci " Ent( ), (4.2)
where ">0 is a parameter, and c is the ground cost, with the associated pairwise cost matrix c ij = c(x i , x j ) evaluated on the support of µ and ⌫. We will assume (as is the case in FWI) that µ and ⌫ have the same support, and will restrict to the quadratic ground cost (i.e. p =2)

c ij = kx i x j k 2 (4.3)
Later, in the context of FWI, x =(t, x) will be a point in time t and offset x.

The problem (4.2)i sa n"-strongly convex program with a unique optimal solution ? " . The solution ?

" converges to the solution with maximal entropy within the set of all optimal solutions of the Kantorovitch problem, i.e. OT " "!0 ! OT.

Problem (4.2) can be reformulated as a projection onto the transport plans of the Gibbs kernel associated to the ground cost with respect to the Kullback-Leibler divergence

OT " (µ, ⌫)= min 2Π(µ,⌫) h , ci + " X i,j ij (log( ij ) 1), =m i n 2Π(µ,⌫) " X i,j KL ij | K " ij . (4.4)
where K " ij def.

=e x p { c ij /"}, is the Gibbs kernel associated to the ground cost, and the Kullback-Leibler divergence is defined as

KL(s | t) def. = 8 > < > : slog( s t ) s for(s, t) > 0, +1 for t =0 (4.5)
The Kullback-Leibler divergence or relative entropy takes two arguments s and t, and is a particular instance of statistical divergences that has nice analytical, and computational properties. As illustrated in Figure 4.1, it is strictly convex with its minimum at s = t, e.g. t =1in this case.

A key insight is that, as " increases, the entropic part of the cost function favours mass spreading and smooths out small scale transport of the non-regularised OT. The optimal coupling becomes less and less sparse. While this has the effect of both accelerating computational algorithms and promoting faster convergence, it introduces a penalisation bias in the transport that tends to destroy the distance properties.

In particular, OT " (µ, µ) > 0 and µ is no more necessarily the minimum of ⌫ 7 ! OT " (µ, ⌫) as it may become cheaper to diffuse from µ than not moving. This is illustrated by the 1D example in Figure 4.3.

Remark on the reference measure In order to simplify the presentation of this chapter, we used the penalization Ent( ij )=KL( ij |1). This is the relative entropy with respect to the uniform measure 1 (Lebesgue). A standard practice is to use KL( ij |µ i ⌫ j ) which yields a smaller value for the the penalization Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 and also a more accurate (in the sense of the non entropic problem) support of the entropic plan. All formulations and algorithms are easily extended to this more general case.

inf ij sup u i ,v j L( ij ,u i ,v j ):= ij c ij + " X i,j ij (log( ij ) 1) u i ( X j ij µ i ) v j ( X i ij ⌫ j ) , (4.6)
The Lagrangian L is convex in ,a n dt h efi r s tv a r i a t i o ni n ij yields:

c ij + " log( ij ) u i v j =0, 8(i, j) 2 [1,n] ⇥ [1,m] (4.7)
where {u i ,i =1,...,n} and {v j ,j =1,...,m} are the Lagrangian multipliers associated to the n-constraints Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 P j ij = µ i and the m-constraints

P i ij = ⌫ j , respectively.
The unique optimal coupling (plan) for entropic regularisation has the specific form

ij =exp ⇢ u i + v j c ij " = a i K " ij b j , 8(i, j) 2 [1,n] ⇥ [1,m] (4.8)
where

K " ij =exp n c ij " o ; a i =exp n u i " o ; b j =exp n v j " o . (4.9)
with a and b being non negative vectors.

Using the marginal constraints in Π(µ, ⌫), it leads to the nonlinear system

a i 0 @ m X j=1 K " ij b j 1 A = µ i , and b j n X i=1 K " ji a i ! = ⌫ j (4.10)
which involves two unknowns a i and b j

a i = µ i P m j=1 K " ij b j and b j = ⌫ j P n i=1 K " ji a i (4.11) or in compact form ( K " =(K " ij ) (i,j)2[1,n]⇥[1,m] .i sn o wam a t r i x )
a (K " b)=µ, and b K "T a = ν (4.12)

where corresponds to the entry side multiplications of vectors, and the unknowns a and b are a = µ ↵ (K " b), and b = ν ↵ (K "T a) (4.13) where ↵ corresponds to the entry side divisions of vectors.

That problem is known in the numerical analysis community as the matrix scaling problem, see Nemirovski and Rothblum (1999) and references therein. An intuitive way is to solve it iteratively. The Sinkhorn algorithm [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF] constructs iteratively the solution

(l+1) = a (l+1) K " b (l) , (l+1) ij = a (l+1) i K " ij b (l+1) j . (4.14)
where a (l+1) = µ ↵ (Kb (l) ) and b The convergence error at iteration l can be monitored using the residual in the Kullback-Leibler sense :

(l+1) = ⌫ ↵ (K T a (l+1) ), a (l+1) i = µ i P m j=1 K " ij b (l) j and b (l+1) j = ⌫ j P n i=1 K " ji b (l+1) i
err (l) KL = n X i=1 KL 0 @ a (l+1) i m X j=1 K " ij b (l) j µ i 1 A + m X j=1 KL b (l)+1 j n X i=1 a (l+1) i K " ji ⌫ j ! . (4.16)
Alternatively, the residual of the Kantorovich potentials in sup norm between successive iterations can also be used It is also known, e.g. Léonard (2014) proposition ( 13), that when " ! 0, the unique minimiser ? " of OT " (µ, ⌫) converges to the maximal entropy plan among the possible optimal transport plans of OT(µ, ⌫).

err (l) 1 = "
An example borrowed from [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF], and reproduced Figure 4.3, illustrates the behaviour of the optimal transport plan when " decreases. The optimal entropic plan ?

" converges to the unregularized optimal transport plan solution of (3.10) and concentrates on the graph of the optimal transport map. The middle line shows how the mass is progressively shifted away from the diagonal during the Sinkhorn iterations for a fixed ".

Numerical stability

When implementing and using the Sinkhorn algorithm it is important to keep in mind that the Kernel K " scales like exp ⌧ d /" where ⌧ is the scale of the ground transport and d the dimension of the physical space, Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 which depends on the data µ and ⌫. The numerical stability therefore depends on these parameters, and the convergence of the Sinkhorn algorithm deteriorates when " ! 0.

The Sinkhorn algorithm often fail to terminate as soon as some of the elements of the kernel K " become too negligible for the machine precision. This can then result in a matrix product K " b or K "T a with very small entries with respect to the machine precision, resulting in a division by 0 numerically in the Sinkhorn update (4.15). Such issues can be partly resolved by carrying out computations in the log domain or using a technique called "-scaling, e.g. [START_REF] Schmitzer | Stabilized Sparse Scaling Algorithms for Entropy Regularized Transport Problems[END_REF].

Cost

Denoting by N the number of points used in the discretisation of the marginal densities,i.e. {a 1 ,...,a N } and {b 1 ,...,b N } a naive Sinkhorn algorithm involves a O n iter ⇥ (n ⇥ m) 2 complexity. where n iter is the number of Sinkhorn iterations to converge, which would be prohibitive in practical FWI applications. For small " the sparsity of the Kernel can be used to reduce the cost, e.g. [START_REF] Schmitzer | Stabilized Sparse Scaling Algorithms for Entropy Regularized Transport Problems[END_REF].

Speed-up for separable kernels

As discussed in Chizat (2017)a n di nPeyré and Cuturi (2019) Remark 4.15, an important particular case, relevant in FWI applications, for which the complexity of each Sinkhorn iteration can be significantly reduced is when each index i and j considered in the ground cost-matrix c ij can be described as a d-uple taken in Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

the cartesian product of d finite sets [1,n 1 ] ,..., [1,n d

]. i =(i k ) d k=1 ,j =(j k ) d k=1 2 [1,n 1 ] ⇥ ...⇥ [1,n d ]
When the ground cost is additive along the sub-indices, namely there exists d matrices C 1 ,...,C d , each of respective size n 1 ⇥ n 1 ,...,n d ⇥ n d such that

c ij = d X k=1 c k i k j k ,
then the kernel appearing in the Sinkhorn iterations has a separable multiplicative structure

K ij =Π d k=1 K k i k j k
leading to fast and exact matrix-vector multiplication, speeding-up the Sinkhorn algorithm.

Let illustrate the Kernel separability in the context of FWI. For 2D horizontal shot-gathers, the predicted d cal (x i ,t j )=d cal ij and the observed d obs (x j ,t j )=d obs ij are discrete realisations at the i =1 ,...,n r receiver spatial positions and the j =1,...,n t time samples, which define on a 2D grid in the space and time domain.

As such d cal and d obs are instantiated as n r ⇥ n t matrices, and as a consequence the Lagrange multiplier u and v. For the quadratic cost, c iljm = |x i x l | 2 + |t j t m | 2 , the associated kernel matrix yields

K " ilmj =exp |x i x l | 2 + |t m t j | 2 " ! = K x il K t mj (4.18)
where K x il =e x p

⇣ |x i x l | 2 /" ⌘ and K t mj =e x p ⇣ |t m t j | 2 /"
⌘ are n r ⇥ n r and n t ⇥ n t convolution matrices defined on the receiver and time finite sets, respectively. In this case a multiplication by K " can be carried out more efficiently.

Such a separable multiplicative structure allows for a fast and exact multiplication by K when applying each 1-D convolution matrix along "slice" of the 2D shot-gather grid. Let rewrite a and b as A and B, respectively, to emphasise the fact that the multipliers are reshaped as n r ⇥ n t matrices. Then computing KB , which would naively require (n r ⇥ n t ) 2 operations with a naive implementation, can be obtained by applying two 1-D convolutions separately, as

⇣ K x K t B T ⌘ T = K x BK t
to recover a n r ⇥ n t matrix in n as a tensor of suitable size. The relative gain is even higher at n 2

r 1 n t + n 2 r 2 n t +(n r 1 + n r 2 ) n 2 t operations instead of n 2 r 1 n 2 r 2 n 2 t .
In order to get a simplified view of the tensorization gain, let us assume n t ' n r := N Acceleration by the Successive Over-Relaxation (SOR)

Since the Sinkhorn algorithm is a fixed-point algorithm, standard linear or even non linear relaxation schemes to enhance the conditioning of the fixed-point mapping near the solution and improve the linear convergence rate. A Successive Over-Relaxation (SOR) algorithm, with an acceleration parameter ✓ 2 [1,[START_REF] Gangbo | The geometry of optimal transportation[END_REF] is presented in Chizat (2017), which yields.

8 > > > > > > < > > > > > > : a (l+1) i =(a (l) i ) 1 ✓ µ i P m j=1 K " ij b (l) j ! ✓ b (l+1) j =(b (l) j ) 1 ✓ ⌫ j P n i=1 K " ji b (l+1) i ! ✓ (4.19)
Numerical experiments in [START_REF] Chizat | Transport optimal de mesures positives: modèles, méthodes numériques, applications[END_REF] illustrate that the convergence rate can be improved by up to orders of magnitude when using SOR, and a detailed convergence analysis is provided, which explicit the best choice of the acceleration parameter ✓. All the experiments in thesis have been performed for ✓ =1.4.

Computation of the Gradient

From (4.6), convex duality yields OT " (µ i ,⌫ j )= sup

(u i ,v j ) = X i u i µ i + X j v j ⌫ j " X ij a i K " (x i ,y j ) b j . (4.20)
As in the non entropic case the dual formulation shows that the map µ 2 R N 7 ! OT " (µ, ⌫) is the supremum of linear functional, hence convex. The gradient is easily deduced as in (3.12),

8(i, j) 2 [1,N] ⇥ [1,N] @ @µ i OT " (µ, ⌫)=u ? i = " log a ? i ; @ @⌫ j OT " (µ, ⌫)=v ? j = " log b ? j (4.21)
where (u ? ,v ? ) are the maximisers of the entropic dual problem (4.20). .

Unbalanced Optimal Transport

In order to deal with seismic unbalanced data we want to drop the hard marginal constraints in Π(µ, ⌫). It allows mass to be transported and also created or destroyed, but at a cost. The "good" additional cost turns out to be the Kullback-Leibler divergence between the data and the marginal of the plan. The Sinkhorn algorithm is very easily adapted to this distance that eliminates the total mass balance constraint.

Recently

Penalisation of the marginal constraints

The hard marginal constraints are replaced by the following soft penalisations:

OT (µ, ⌫)=min N X i=1 N X j=1 c ij ij + 1 N X i=1 KL 0 @ N X j=1 ij µ i 1 A + 2 N X j=1 KL N X i=1 ij ⌫ j ! (4.22)
where KL is the Kullback-Leibler divergence defined in (4.5). The 1 > 0 and 2 > 0 are parameters quantifying the cost of unbalanced mass relaxation. If µ and ⌫ are balanced, the classical optimal transport can be recovered by letting 1,2 ! +1.

A nice interpretation of (4.22) is that, depending on the value of , available mass on both side will be either transported or just created or destroyed locally if it is cheaper through the KL penalization of the marginal constraint.

In the case of the quadratic ground cost (4.3), the value function OT is called the Gaussian Hellinger-Kantorovich(GHK) distance, [START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF].

Entropic regularised unbalanced OT problem

Entropic regularisation can also be applied to (4.22)a si n( 4.4), which yields:

OT ", (µ, ⌫)= min 2Π(µ,⌫) " KL ⇣ K " ⌘ + 1 N X i=1 KL 0 @ N X j=1 ij µ i 1 A + 2 N X j=1 KL N X i=1 ij ⌫ j ! . (4.23)
The dual formulation of the entropic regularised unbalanced optimal transport problem in discrete form follows again from convex duality OT ", (µ i ,⌫ j )= sup As explained in Chizat et al. (2018b) the optimal plans are again in the form (4.8) and can be computed iteratively using the modified Sinkhorn iterations :

(u i ,v j ) 1 X i ✓ exp ⇢ u i 1 1 ◆ µ i 2 X j ✓ exp ⇢ v j 2 1 ◆ ⌫ j " X ij a i K " (x i ,y j ) b j .
8 > > > > > > > > > > < > > > > > > > > > > : a l+1 i = µ i X j (K " (x i ,y j )b l j 1 1 + " , 8i, b l+1 j = ⌫ j X i K T " (x i ,y j ) a l+1 i 2 2 + " , 8j. (4.25)
As before a i =e x p{u i /"} and b j =e x p{v j /"} and the Sinkhorn algorithm obeys to the same convergence properties. In practice we use the error (4.17) to monitor the convergence.

The dual formulation (4.24) shows again that the map µ i 7 ! OT ", (µ i ,⌫ j ) is concave. Applying Danskin's theorem once more the differentials of OT ", (µ i ,⌫ j ) with respect to µ i and ⌫ i can be directly obtained as the differentials with respect to µ i and ⌫ i of the right hand side of (4.24) for the optimal (u ? i ,v ? i )s:

@ @µ i OT ", (µ, ⌫)= 1 ✓ exp ⇢ u ? i 1 1 ◆ = 1 ⇣ a "/ 1 i 1 ⌘ @ @⌫ j OT ", (µ, ⌫)= 2 ✓ exp ⇢ v ? j 1 1 ◆ = 2 ⇣ b "/ 2 i 1 ⌘ .
(4.26)

Sinkhorn divergence

Sinkhorn divergence for balanced OT problems

As previously seen in Section 4.1, the entropic regularised optimal transport problem (4.4) can be efficiently solved with the Sinkhorn algorithm as long as the regularisation parameter " is large enough. The larger the " is, the faster the convergence of the iterative Sinkhorn algorithm. The entropic regularisation however introduces a bias increasing with ".

In the continuous setting asymptotic results by [START_REF] Conforti | A formula for the time derivative of the entropic cost and applications[END_REF]a n dPal (2019)o nO T "

when " ! 0, reminding that OT " ! OT when " ! 0, offer a better understanding of the entropic bias

OT " (µ, ⌫) OT(µ, ⌫) ' 2 " p 2⇡"+ " [KL (µ | ⌫)+KL (⌫ | ⌫)] + O(" 2 ). (4.27)
This result holds for measures (µ, ⌫) which have smooth densities with respect to the Lebesgue measure, still denoted here as µ, ⌫ by abuse of notation.

The definition of the KL divergence in (4.5)i sa l s oa na b u s eo fn o t a t i o n ,

KL (µ | ⌫) def. = 8 > > < > > : Z µ(x) ✓ log ⇢ µ(x) ⌫(x) 1 ◆ dx if µ absolutely continuous w.r.t ⌫.
+1 otherwise (4.28)

It is worth to note that the densities µ and ⌫ are decoupled in the first-order term. A simple idea to remove Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

these terms is to use the symmetric entropic OT between the marginals and to remove them, keeping in mind that OT(µ, µ)=0

S " (µ, ⌫)=OT " (µ, ⌫) 1 2 {OT " (µ, µ)+OT " (⌫, ⌫)} (4.29)
which is, formally at least, a second order approximation of OT(µ, ⌫) with respect to ".

This formulation has been studied under the name "Sinkhorn divergence" in [START_REF] Genevay | Learning generative models with sinkhorn divergences[END_REF][START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF]. It does not define a distance but has nice properties. We immediately see that S " (µ, µ)=0for

any ">0, a property that was lost for OT " . It is also proven in [START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF]; [START_REF] Feydy | Interpolating between optimal transport and mmd using sinkhorn divergences[END_REF] that it is symmetric in µ and ⌫, and remains positive and convex with respect to µ and ⌫. It also metrises the weak convergence of measures.

The Sinkhorn divergence (4.29) can be computed as well as its gradient just by applying the Sinkhorn algorithm three times. For small ", The Sinkhorn divergence S " for the quadratic ground cost is therefore a better approximation of OT 0 = W 2 2 than OT " , inheriting the good properties of Entropic Optimal Transport: strict convexity, smoothness and the Sinkhorn computational algorithm.

Sinkhorn divergence for unbalanced OT problem

The Sinkhorn divergence has a natural extension to the unbalanced case [START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF], yielding S ", (µ, ⌫)=OT ", (µ, ⌫)

1 2 {OT ", (µ, µ)+OT ", (⌫, ⌫)} + " 2 (m(µ) m(⌫)) 2 . (4.30)
where OT ", (µ, ⌫) is the entropic regularised unbalanced optimal transport problem defined in (4.23), while m(µ)= P i µ i and m(⌫)= P j ⌫ j are the total mass of µ and ⌫, respectively, which can be different. It is proved in [START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF]thatS ", (µ, ⌫) is again positive, definite, and convex. A comparison between OT ", (µ, ⌫) with S ", (µ, ⌫) is illustrated using 1D wavelets in the next chapter.

Computing the unbalanced Sinkhorn divergence and its gradient

In order to solve problem (4.30), let first look at OT ", (µ, ⌫). Following section 4.2.2,O T ", (µ, ⌫) can be solved as a fixed-point problem with the Sinkhorn iterations 4.25. For the symmetric case OT ", (µ, µ),t h e Sinkhorn algorithm simplifies [START_REF] Séjourné | Sinkhorn divergences for unbalanced optimal transport[END_REF].

The dual problem becomes a concave maximisation problem, symmetric with respect to a unique unknown vector u i :

OT ", (µ i ,µ j )=sup (u i 2 X i ✓ exp ⇢ u i 1 1 ◆ µ i " X ij a i K " (x i ,y j ) U j . (4.31)
where a i =exp( u i /"). At iteration l +1, the optimal vector a l+1 i can be computed as

a l+1 i = µ i P j K " (x i ,x j )a l j ! + " (4.32)
Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

As in (4.26), the first variation of OT ", (µ, µ) with respect to µ can be directly obtained from the dual formula (4.31):

@ @µ i OT ", (µ, µ)= 2 ⇣ a "/ i 1 ⌘ . (4.33)
The Sinkhorn divergence S ", can be computed using the Sinkhorn iterations (4.25)tosolv ethestandard problem OT ", (µ, ⌫) and using the Sinkhorn iterations (4.32) to solve the symmetric problem OT ", (µ, µ)

and OT ", (⌫, ⌫).

The gradient of the unbalanced Sinkhorn divergence S ", (µ i ,⌫ j ) with respect to µ i is finally given by:

@ @µ i S ", (µ, ⌫)= @ @µ i OT ", (µ, ⌫) 1 2 @ @µ i OT ", (µ, µ)+" (m(µ) m(⌫)) (4.34)
where the first term is given by ( 4.26) and the second term by (4.33).

Metric and domain scaling

A 2D illustration of metric scaling

In the context of FWI, the data are just pixelized images to be compared. Using a local misfit function (L 2

for instance) the relative size of the domain in time and offset does not impact the minimization. This is not the case when using a non local OT misfit.

When comparing seismograms, a natural question is how to measure the displacement of mass in time and offset space and how to choose the dimensions of along these two axis. This is equivalent to use a normalized domain [0, 1] 2 and scaling the ground cost (4.3)( x =(t, x)) We introduce a rescaled cost

c ijml = cos 2 (✓ M ) |x i x m | 2 +sin 2 (✓ M ) |t j t l | 2 (4.35)
depending on an angle ✓ M 2]0,⇡/2[. We are going to discuss how the translation of specific data depends on ✓ M .

Let the observed data

µ(t, x)= {cos(✓ D ) x+sin(✓ D ) t=0}
represent an idealized plane wave event. A single straight line carrying a weight of 1 and making an angle ✓ D with the offset axis. The other part of the data (the simulation):

⌫ ⌧ (t, x)=µ(t ⌧ sin(✓ D ),x ⌧ cos(✓ D ))
is a simple translation depending on a parameter ⌧ in the direction (cos(✓ D ), sin(✓ D )). This emulates the misfit associated with an erroneous constant background velocity and as explained in the introduction will cause cycle skipping when using L 2 misfit.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Assuming for simplicity that we are in a periodic box, the classical (non entropic) optimal transport map is a pure parallel translation of the support of µ. This is a consequence of Brenier's theorem which can be extended to ground costs formed as convex functions the vector y x such as (4.35) [START_REF] Gangbo | The geometry of optimal transportation[END_REF]. The translation may however not occur in the normal direction ✓ D + ⇡ 2 because the shortest path now depends on the distance induced by (4.35)a n dt h ea n g l e✓ M .

The OT misfit can therefore be reduced to computing the shortest path squared from the origin to the translated line. It is given by solving the convex program:

d 2 ⌧ = argmin (x,t),s.t. cos(✓ D ) x+sin(✓ D ) t=⌧ cos(✓ M ) 2 2 |x| 2 + sin(✓ M ) 2 2 |t| 2 = C ⌧ 2 2 (4.36)
where

C = 1 cos(✓ D ) 2 cos(✓ M ) 2 + sin(✓ D ) 2 sin(✓ M ) 2 .
(4.37) Real or synthetic complex seismic data may be interpreted a collection of local event possibly misplaced.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

The discussion above raises the interesting but still untouched question of building a data dependent ground cost function where (↵, ) depends on (t, x). One could even consider more general metrics such as Finsler metric for example Benamou et al. (2018).

On the (t, x) range of Entropic OT plans

In the entropic version of OT, the plan (4.8) is a diagonal scaling of the interaction Kernel which "governs" mass displacement. In its separable version (4.18) it becomes

K ijml = e |x i x m | 2 2 " e |t j t l | 2 2 " . (4.38)
Let us discuss the choice of the relative size of the domain in time and offset and their discretization. Let (T,N t ) and (X, N x ) be the maximum time, respectively offset, and number of points in time, respectively offset, used in the computation of the misfit. That gives a computational grid size (dt, dx)=( T Nt , X Nx ). Setting k x dx = x i x m and k t dt = t j t l we rewrite (4.38):

K " ijml = e (k x dx) 2 2 " e (k t d t ) 2 2 " = e (k x X) 2 2 "N 2 x e (k t T ) 2 2 "N 2 t (4.39)
Freezing " and assuming (N t ,N x ) are given either by the data (number of lines) or though time sampling to limit the computational cost. The range of "horizontal" (in offset) or respectively "vertical" (in time) mass displacement is therefore governed by (k t T,k x X). Setting the finite precision 0 at exp 36 the Kernel above will vanish when

k t T =6 p 2 "N t or k x X =6 p 2 "N x .
Entropic Optimal Transport therefore allows to choose a smooth transport window of grid points in time and offset. It is possible to use (T,X) to control the width of this window. Let (k max t ,k max x ) be a prescribed maximum allowed displacement in terms of grid points and set .40) (also note that (T,X) can be interpreted as penalization of time/offset displacement in the ground cost as in the previous section ).

T =6 p 2 "N t k max t or X =6 p 2 "N x k max x . ( 4

Trace by trace misfit with entropic OT

We recall that the system solved by the Sinkhorn algorithm (4.11) is equivalent (in tensorized form (4.38

)) to the marginal constraints: Setting k max x =1in particular corresponds to 1D "trace by trace" misfit. The offset part of the Kernel K x im will vanish except when k x =0, i.e. x i = x m ,i n (4.38) where it is just 1. The 2D marginal constraints (4.41) will therefore reduce for all offsets x i to the 1D marginal constraints : The unbalanced optimal transport misfit function

8 > > > < > > > : X ml A ij K x im K t jl B ml = µ ij 8 i, j X ij A ij K x im K t jl B ml = ⌫ ml 8m, l (
8 > > > < > > > : X l A ij K t jl B il = µ ij 8 j X j A ij K t jl B il = ⌫ il 8 l
The general abstract form of the misfit to minimize (see (2.8)) is

J (m)=h (d cal (m), d obs ) (5.1)
where m is the model we strive to reconstruct and d obs is some given observation. The map m 7 ! d cal (m) is the forward model and in this thesis we used a 2D constant-density acoustic finite difference wave modeling.

We also limit ourselves to synthetic data sets that is d obs := d obs (m true ), the observation is constructed with the same forward mode ling using a true model which is known a priori.

We recall that we have implemented on top of PySIT (https://pysit.org/) is an "an open source toolbox for seismic inversion and seismic imaging developed by Russell J. Hewett and Laurent Demanet at MIT".

It was designed to allow for fast prototyping of the h component of the misfit function (5.1). Developers can implement their own misfit on top of the available forward models (in our case 2D constant-density acoustic finite difference wave modeling) and optimisation method (we have been using the Limited-memory BFGS method).

The rest of this chapter explains how we constructed the h function using the entropic OT costs and shows its behavior on simple parametric models.

The misfit function

In an attempt to preserve phase information, we follow the strategy of comparing separately the positive and negative parts (f + = max(f, 0) and f = max( f, 0)) of the observed and predicted seismic signal as proposed in [START_REF] Engquist | Application of the wasserstein metric to seismic signals[END_REF][START_REF] Engquist | Optimal transport for seismic full waveform inversion[END_REF]. This is to take care of the positivity constraint on the data. As explained in chapter 4 instead of normalizing and using the classic W (5.3) This will not be used in this chapter but in order to ensure misfit function differentiability, we will split the signal, here denoted f ,i n t oas m o o t hp o s i t i v ef + = P (f ) and smooth negative f = P ( f ) part, where P is given by

P (f )= 1 2 ⇣ f + p f 2 +( f 0 ) 2 ⌘ (5.4)
The parameter is a small parameter that controls the smoothness of P and f 0 carries the dimension of the signal (typically f 0 = kd obs k 1 .). This function is plotted in in Figure 5.1, for a simple function f (t)=t, f 0 =1and different values of . In summary the misfit function, used for FWI in the last chapter, is h S ", (d cal , d obs )=S ", (d cal ) + , (d obs ) + + S ", (d cal ) , (d obs ) .

( 5.5) and for completeness the adjoint-source to be used in the adjoint wave propagation (2.42) is given by

@ @d cal h (d cal ; d obs )= @ @d cal S ", (d cal ) + , (d obs ) + @ @d cal P (d cal ) @ @d cal S ", (d cal ) , (d obs ) @ @d cal P ( d cal ) .
( 5.6) where the differentials of the Sinkhorn divergence (4.30) with respect to d cal is given by (4.34).

Figure 5.2 shows that the smoothness of adjoint source, given by the pixel wise multiplication of the S ", gradient and the P derivative depends on the parameter .

1/2D parametric examples

In this section we skip the smoothing related to the sign splitting of the data. We will use misfits (5.2)a n d ( 5.3) and they will be referred to as OT ", and S ", . We will omit the when the data is balanced and there is no need for unbalanced OT.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Translation/shift of a 1D signal

We are considering two copies of the same Ricker wavelet, with 20Hz peak frequency, denoted f (t),p l a y i n g the role of the observed signal d obs (t),a n dg(t), playing the role of the predicted signal d cal (t). The two wavelets are shifted in time such that g(t; s)=f (t s) see Figure 5.3. We will study the effect on the misfit of the variations in s the time shift.

We recall the classical result (see section 3.2): the squared Wasserstein distance W 2 2 is quadratic, hence convex with respect to the translation of identical data. Sign splitting will preserve the mass balance and the translation and as such will not change this property. This was already verified numerically in [START_REF] Engquist | Application of the wasserstein metric to seismic signals[END_REF]. We here compare the behaviour of the L 2 , OT " ,a n dS " misfit functions, Figure 5.4.I n this example, the data is balanced so we do not need the unbalanced version of the transport distance. For such a simple example there is not much difference between the entropic OT " and the Sinkhorn divergence but the debiasing effect discussed in section 4.3.1 is noticeable. We are very close to a perfect s 7 ! 2 s 2 (the 2 factor comes from sign splitting). Sinkhorn divergence is known to be a good proxy for non entropic OT.

In all experiments " = 10 2 and = 10. In the next experiment, Figure 5.6, we consider the misfit behaviour as a function of the time shift between the two Ricker wavelets when a Gaussian noise is added to one of the wavelet (linked for instance to the acquisition of d obs ), see Figure 5.5.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 We compare the L 2 and the unbalanced OT ", and S ", misfit functions for an increasing level of Gaussian noise (1, 5, 10, 20 %). The data is now unbalanced because of the noise. We also consider the balanced S " after normalisation of the data (as in [START_REF] Engquist | Application of the wasserstein metric to seismic signals[END_REF][START_REF] Engquist | Optimal transport for seismic full waveform inversion[END_REF]).

As expected L 2 global minima remains centred but more local minima appear due to the noise. The balanced (normalised) S " remains convex as expected but the global minimum is shifted up because of the noise. It is indeed built to be positive convex and vanish when comparing two identical measures, this never is the case with added noise. A constant is added for all the mass creation. The position of the minimum depicted by a black cross is biased because the noisy and normalized signal does not satisfy the initial Ricker symmetry anymore.

We recall that unbalanced OT remains a distance on radon measures, the marginal constraint on the data is penalised using the Kullback-Leibler divergence KL see (4.22)( w eh a v eu s e d 1 = 2 = ). Unbalanced Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Sinkhorn Divergence S ✏, is also covex positive and corrects the unbalanced OT ✏, by removing the diagonal terms and therefore achieves a lower miminum. Again it cannot be zero as we are comparing different measures. The unbalanced version seems however to do a better job as the minimum is closer to the zero shift.

Finally, we remark that the noise is akin to add mass everywhere therefore impacting the transport.

Because we use the quadratic ground cost (4.3), the transport will favor moving the nearby noise mass when possible instead of fetching the shifted Ricker. This effect increases with the noise level, reduces the range of the transport and the modulus of convexity of the misfit. In figure 5.7 and or a fixed noise (5%) we explore the behavior of the different misfit for different values of ✏ and .I nt h efi r s tl i n e✏ is decreased, as expected OT ✏, and S ✏, (the second order approximation in ✏ of OT converge). We also notice that normalizing the data induce a significant bias to the minimum. The second line, when the marginal constraint are locally relaxed, the relative diminution of the mass transport leads to degeneracy of the convexity property with respect to translation (the amount of mass beeing transported decrease).

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 We perform the same experiment, Figure 5.9, when instead of adding noise to the recorded wavelet, we change its frequency content, see Figure 5.8. We roughly have the same observation as in figure 5.7.W e notice that the normalisation approach does much better, probably because the support of the wavelet is localised.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Linear 2D model

A simple two-parameter model, borrowed from Métivier et al. (2016a), allows the 2D representation of the L 2 ,O T ", and S ", misfit functions. We use sign splitting for the positivity for the transport misfits. Let consider a 2D-acoustic wave approximation model, with constant density, parameterised by the background acoustic velocity v 0 and a depth gradient ↵,s u c ht h a t

v = v 0 + ↵z (5.7)
The physical domain is a 2D rectangular domain 17 km long and 3.5 km deep. A single seismic source is located at x =8 .45 km at a depth of 50 m. The radiated acoustic waves are recorded at 168 receivers located at the same depth position that the source and regularly deployed with a spacing of 100 m for

x 2 [0.15, 16, 85] km. The source wavelet is a Ricker wavelet centred at 5 Hz.

The observed signals d obs (x r ,t) correspond to a numerical physical realisation, using second-order finite difference stencils on a regular finite difference grid and a second-order leap-from time integration, for a model parameterised with v ? 0 =2km/s and ↵ ? =0.67 s 1 . The L 2 ,O T ", and S ", misfit functions are compared for predicted signals d cal (x r ,t) corresponding to numerical physical realisations when exploring the parameter space v 0 2 [1.75, 2.25] km/s and ↵ 2 [0.49, 0.91] s 1 with ∆v 0 =0.0125 km/s and ∆↵ =0.015 s 1 , leading to a discretised model space of 29 ⇥ 41 points. For the OT misfit functions, the entropic parameter is chosen as " =0.01.

The L 2 misfit function, in Figure 5.10, exhibits multiple local minima, and convergence to the global minimum with local gradient-descent methods requires to start with a good enough initial model, otherwise the optimisation problem will be trapped in a local minimum.

Institut de Physique du Globe de Paris

Miao YU, Ph.D. Thesis, 2021 The S ", misfit function, in Figure 5.12, exhibits nice convexity and sensitivity properties, together with a large basin of attraction toward the global minimum. It is worth to remind that while the mapping d cal 7 ! S ", (d obs ,d cal ) is strictly convex, here this is the mapping (v 0 ,↵) 7 ! S ", (d obs ,d cal (v 0 ,↵)) that is actually plotted.

Institut de Physique du Globe de Paris

Miao YU, Ph.D. Thesis, 2021 2D FWI numerical illustrations

Preliminaries

This section provides numerical illustrations for the use of the misfit function, see (5.6), based on the Sinkhorn divergence formulation of unbalanced, entropic optimal transport , see tol = 10 09 Sinkhorn convergence error (4.17). With the choice of parameters above and an optimal transport grid size considered of the order 10 2 in space and 128 in time, the number of iterations to achieve this precision is always of the order of 10 2 iterations and never exceeded 10 3 iterations. Transport window in space (in number of points) (see section 4.4.2).

Velocity anomaly inclusion

This is the classical circular velocity anomaly inclusion inclusion model in a homogeneous medium, and is inspired from [START_REF] Pladys | On cycle-skipping function modification for fullwaveform inversion: comparison of five recent approaches[END_REF]. The physical model is defined on a 1000 ⇥ 1000 ms q u a r ed o m a i n ,a n di s parametrised by a homogeneous background acoustic wave velocity set to 1300m/s, and a circular inclusion of 100 m radius centred in space, with a homogeneous wave velocity 1700m/s (⇡ 30% anomaly).

The computational domain includes additional perfectly matched layers along each of the boundaries.

The constant-density, acoustic wave equation is solved in space on a regular grid of 256 ⇥ 256 points using a sixth-order, centred finite difference scheme and a second-order leap-frog integration scheme in time with constant time step satisfying the CFL condition, as implemented in PySIT.

The FWI inversion problem is solved using the l-BFGS method, as implemented in PySIT, and the number of FWI iterations is initially fixed to 25 iterations. Two homogeneous initial models are considered: initial model-A with acoustic wave velocity set to 1300 m/s; and initial model-B with acoustic wave velocity set to 1700 m/s.

Both the predicted d cal and the observed d obs data, i.e. the latter corresponding to the "true" model, are physical realisations obtained with the same wave propagation solver, and without considering additional noise to the data, i.e. the inversion crime setting.

Unless otherwise specified, the S ", -misfit function is constructed as discussed in Chapter 5. The observed and predicted signals are decomposed into positive and negative parts using the P -transforn in (5.4), and transported separately using the Sinkhorn divergence formulation of unbalanced optimal transport.

The numerical experiments focus on transmission and reflection configurations.

Transmission configuration

The acquisition system is composed of 48 sources evenly-spaced below the upper boundary of the physical As expected, the resolution for both initial models is higher in the horizontal direction that in the vertical direction due to the acquisition geometry, with a vertical smoothing of the reconstructed velocity anomaly.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

The shape of the inclusion is not fully recovered possibly due to limited offset and deficit in long-wavelength components of the transmitted waves and lack of reflected wave information.

Transmission (initial B) 25 iterations 100 iterations convergence convergence The S ", -misfit function provides reasonably accurate results, and the observed and predicted seismic traces are in phase after 25 iterations, as illustrated in Figure 6.5, even though the exact background velocity is not fully recovered.

As expected, the resolution is higher in the horizontal direction that in the vertical direction due to the acquisition geometry, leading to vertical smoothing of the reconstructed velocity anomaly.

The peak-amplitude of the velocity anomaly is not fully recovered after 25 iterations, see Figure 6.4 and With both the L 2 and S ", based inversion, the shape of the inclusion is not fully reconstructed, due limited offset and relative deficit in long-wavelength components. The resolution in the horizontal direction is higher than in the vertical one, leading to vertical smoothing.

The L 2 and S ", adjoint-sources are different, as shown in Figure 6.10, leading to different sensitivity kernels.

L 2 Gradient S ", Gradient The S ", adjoint source is more complicated with the signature of a mass-spreading over the whole domain. This is primarily due to the penalisation of the marginals in the unbalanced Sinkhorn divergence formulation, and more specifically to the gradient of Kullback-Leibler divergence, and possibly to the separate 2D-transport of the positive and negative parts of the seismic signals.

The initial and final predicted data, in the time-receiver domain, are shown in Figure 6.11 together with the "true" data, i.e. the observed data. A diffraction pattern can be observed as well as a resolution deficit in the long wavelengths components, which precludes accurate reconstruction of the inclusion shape. FWI results, using different implementations of the S ", -misfit function are shown in Figure 6.12 and S ", (2D -P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P )S ", (1D-Time-P ) S ", (2D-P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P ) SD (1D-Time-P ) Figure 6.13: Velo city profile along a vertical section of the physical domain for different FWI iterations and different implementations of the S ", -misfit function.

Reflection configuration

The The peak-amplitude of the velocity anomaly is reasonably well retrieved, see Figure 6.17 The initial and final predicted data are shown in the time-receiver domain Figure 6.18. Reflection phase associated to the lower boundary of the inclusion, and the weaker multiple reflections behind, are poorly resolved. The background velocity also is not exactly retrieved.

The S ", based inversion exhibits a monotonous convergence rate and the misfit keeps continuing to decrease after 25 iterations, see Figure 6.16.

Initial model B

Reconstructed models after 25 iterations with the L 2 and the S ", based inversion are shown in Figure 6.19.

The L 2 -misfit function now does not suffer from cycle-skipping, and both the L2 and S ", based inversion reconstruct meaningful models. The observed and predicted data are in phase, see Figure 6.9 for the S ",misfit function.

The peak-amplitude of the velocity anomaly is not yet fully recovered after 25 iterations with the L 2misfit function, whereas with the S ", -misfit function a relatively more accurate peak-amplitude is recovered, see Figures 6.19 L 2 Gradient S ", Gradient S ", (2D -P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P )S ", (1D-Time-P ) S ", (2D-P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P ) SD (1D-Time-P ) 6.3 Layered reflection Models

One-layer reflection model

The one-layer model is inspired from [START_REF] Pladys | On cycle-skipping function modification for fullwaveform inversion: comparison of five recent approaches[END_REF].

The physical model is defined on a rectangular domain of 2000 ⇥ 1000 m dimension, and is parametrised by a homogeneous background acoustic wave velocity set to1500 m/s and by a 100 m thick layer at 300 m depth with a homogeneous wave velocity of 1600 m/s (postive anomaly) or 1400 m/s (negative anomaly).

The surface acquisition is the same for both models with 48 evenly-spaced seismic sources and 251 evenly- Reconstructed models, after 50 iterations, with the L 2 and S ", based inversion are shown in Figure 6.27.

The L 2 and S ", based inversion provide similar meaningful reconstructions.

Results for the L 2 based inversion were expected since the L 2 -misfit function is sensitive to amplitude variation and polarity. Results for the S ", -based inversion are less intuitive since the S ", -misfit function is based on signal decomposition into positive and negative parts and does not guarantees the preservation of phase and polarity information. The data are in phase for the S ", based inversion, after 50 iterations, see Figure 6.28.

The initial and final predicted data for both L 2 and S ", based inversion are shown in the time-receiver domain Figure 6.32.

The peak-amplitude of the velocity anomaly is correctly retrieved by the L 2 based inversion after 50

iterations, but only partially retrieved by the S ", based inversion, see Figure 6.29. The S ", based inversion recovers only half of the anomaly, with relatively more vertical smoothing than in the L 2 case, see also Figure This could possibly result from the entropic regularisation, which tends to spread mass and smooth discontinuity. The S ", based inversion shows slower convergence rate compared to the L 2 based inversion, see Figure 6.30, but the misfit distance keep continuing to decrease after 50 iterations.

The L 2 and S ", adjoint sources and optimisation gradients, at the first iteration, are shown in Figure 6.31. The S ", adjoint source shows the signature of a mass-spreading over the whole domain as previously discussed, and the gradient is smoother than the L 2 one, even though both capture similar features. Results obtained with the offset-by-offset S ", , i.e. 1-D unbalanced transport in space, are less satisfactory.

In particular the reconstruction of the positive anomaly layer is altered by strong horizontal smoothing artefacts.

S ", (2D -P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P )S ", (1D-Time-P ) Reconstructed models, after 50 iterations, with the L 2 and S ", based inversion are shown in Figure 6.27.

Both inversion provide similar meaningful reconstructions.

For the L 2 based inversion, the recovered peak-amplitude of the velocity anomaly slightly exceeds the "true" one, see Figure 6.35, whereas for the S ", based inversion, the peak-amplitude is only partially recovered after 50 iterations.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

S ", (2D -P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P )S ", (1D-Time-P ) The L 2 and S ", adjoint sources and optimisation gradients, at the first iteration, are shown in Figure 6.37. The gradient associated to S ", appears smoother compared to the L 2 one, even though both capture similar features.

The convergence rate of the S ", based inversion is similar to the L 2 one, and the misfit distance keep continuing to decrease after 50 iterations see Figure 6.36.

Inversion results for different implementations of the S ✏, -misfit function are shown in Figure 6.38 and The results obtained with the offset-by-offset S ", -misfit function, i.e. 1-D unbalanced transport in space, are more meaningfull than in the case of a positive velocity anomaly, see section 6.3.2.

S ", (2D-P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P ) SD (1D-Time-P ) Figure 6.39: Velo city profile along a vertical section of the physical domain for different FWI iterations and different implementations of the S ", -misfit function.

Three-layer reflection model

Model configuration

The three-layer model is inspired from [START_REF] Yang | Optimal transport for seismic inverse problems[END_REF].

The physical model is defined on a rectangular domain of 15 ⇥ 6 km dimension, with three homogeneous layers: a first layer [0, 1] km, with a homogeneous wave velocity set to 1 km/s; a second layer [1,[START_REF] Gangbo | The geometry of optimal transportation[END_REF] km with a homogeneous wave velocity set to 2 km/s; and a third layer Both inversions converge slowly, see Figure 6.43, with a relatively better convergence rate, after 500 iterations for S ", .

The first arrivals in the observed and predicted data are in phase but smaller-amplitude reflected phases behind are not properly recovered for larger offset, see Figure 6.44. A closer look at the reconstructed velocity models during the L2 and S ", based inversion iterations, see It is not intuitive that the S ", based inversion can recover velocity information in the third layer while no waves through this region return to the receivers. Velocity information of the third layer is printed in the amplitude of the reflected phases associated to the third-layer interface and in the head waves, the later being more easily retrieved with the 2D S ", -misfit function, even for limited offset.

As the S ", based inversion gradually recovers deep velocity information, the reconstructed velocity model in the shallower region exhibits higher-amplitude, zero-mean oscillations around the "true" velocity model, compared to the L 2 based inversion. The amplitude of these oscillations slowly decrease as inversion iterate.

Both the L2 and S ", based inversion overshoot the peak-amplitude of the velocity contrast at the third-layer interface.

The initial and final predicted data are shown in the time-receiver domain Figure 6.46 for both the L 2 and S ", based inversion. The L 2 and S ", adjoint sources and optimisation gradients, at the first iteration, are shown in Figure 6.45.

The initial predicted data do not contain the reflected phases at the third-layer interface, which has a strong signature in the L 2 and the S ", adjoint sources beside kinematic errors on reflected phases associated to the second-layer interface and the diving waves.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 The S ", contains additional signature of a mass-spreading over the whole domain. This can be associated to the penalisation of the marginals in the Sinkhorn divergence formulation of unbalanced entropic optimal transport, and more specifically of the local gradient of the Kullback-Leibler. The level of mass is linked to how mass-unbalanced are the predicted and observed data, and as such to the initial model. For strongly unbalanced observed and predicted data the transport cost relatively decreases. The gradient in the S ", based inversion is smoother at the third-layer interface, and with with less subsurface artefacts compared to L 2 based inversion.

Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

S ", (2D -P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P )S ", (1D-Time-P ) The velocity contrast at the third-layer interface is not correctly retrieved. This is not actually surprising as reflected waves polarity are lost in the squared-transformation.

Interestingly, results obtained with the offset-by-offset S ", -misfit function, i.e. The amplitude of the zero-mean oscillations in the subsurface region seems better mitigated than with 2D transport.

S ", (2D-P )S ", (2D-Squared Amp.)

S ", (1D-Offset-P ) SD (1D-Time-P )

Model configuration

This model is derived from the synthetic Marmousi I, originally generated at the Institut Français du Petrole (IFP) using a 2D acoustic wave solver [START_REF] Versteeg | The marmousi experience: velocity model determination on a synthetic complex data set[END_REF].

The model is defined on a 9000 ⇥ 3000 m domain and contains many reflectors, steep dips, and strong velocity variations in bot lateral and vertical dimensions, with a minimum velocity of 1500 m/s and a maximum velocity of 5500 m/s.

The physical model here is defined on a mini-square domain extracted from the Marmousi I model, with origin x = 4680 ma n d3000 ⇥ 3000 m dimension. The constant density, acoustic wave velocity model is shown in Figure 6.49

The surface acquisition is composed of 48 evenly-spaced seismic sources and 150 evenly spaced receivers, both at 500 m depth below the upper interface. The sources are point sources and the source-time function is a Ricker wavelet with central frequency f c = 10Hz, and the recording time is set to 3 s.

The computational domain includes additional perfectly matched layers along each side of the physical domain. The constant-density acoustic wave equation is solved in space on a regular grid with 150 ⇥ 150 points, i.e. dx = dz = 20 m, using a sixth-order finite difference scheme and in time using a second-order leap-frog integration scheme with constant time step satisfying the CFL condition, as implemented in PySIT.

The inversion problem is solved using the l-BFGS gradient-based method, as implemented in PySIT, with a fixed number of FWI iterations set to 200 iterations.

The initial model, Figure 6.49, is derived from the true velocity model using a smooth low-pass filter with a frequency of 1/300, which preserves part of the long-wavelength components of the true model. Reconstructed models with the L 2 and S ", based inversion are presented in Figure 6.50.

Both the L 2 and S ", based inversion recover similar velocity models in the shallow region (above 1. The initial and final predicted data are shown in the time-receiver domain Figure 6.51 for both the L 2 and S ", based inversion.The L 2 and S ", adjoint sources and optimisation gradients, at the first iteration, are shown in Figure 6.52.

The predicted data at the first iteration miss, as expected, the reflected phases and are only reasonably in phase at short offset for the diving waves. This is reflected into the L 2 adjoint source.

The S ", adjoint source is more complicated with signature of mass-spreading over the whole domain as a result of mass-unbalanced predicted and observed data induced by the smooth initial model.

The gradient in the L 2 and S ", based inversion are different and smoother for Sl ", .S i g n a t u r eo ft h e limited offset is visible in both cases.

The convergence rate of the S ", based inversion, over 200 iterations, is relatively faster than the L 2 one, rapidly reducing the misfit distance by one order of magnitude after 20 iterations. The misfit distance keeps continuing slowly decreaseing after 200 iterations. The computational domain includes additional perfectly matched layers along each side of the physical domain. The constant-density acoustic wave equation is solved in space on a regular grid of 176 ⇥ 150 points, with dx = dz = 20m, using a sixth-order finite difference scheme and in time using a second-order leap-frog integration scheme with constant time step satisfying the CFL condition, as implemented in PySIT.

The inversion problem is solved using the l-BFGS gradient-based method, as implemented in PySIT, with a fixed number of FWI iterations set to 250 iterations.

Two initial models are considered: a S-500 initial model, which is derived from the "true" velocity model applying a smooth low-pass filter with a frequency set to 1/500; a 1D-GRAD initial model, which is a homogeneous model with a linearly depth-increasing velocity from 1500 m/s at the bottom of the water layer to 3500 m/s at depth. While the S-500 initial model conserves some of the long-wavelength content of the true model, the 1D-GRAD initial does not contain anymore signature of the long wavelength components of the true model. With both initial models it is a challenging task for FWI to recover large low-wavenumber discrepancies between the initial and "true" model.

Both the predicted d cal and the observed d obs data, i.e. the latter corresponding to the "true" model, are physical realisations obtained with the same wave propagation solver, and without considering additional noise to the data, i.e. the inversion crime setting.

Unless otherwise specified, the S ", -misfit function is formulated as in Chapter 5. The observed and predicted signals are each decomposed into positive and negative parts, using the P -transform defined in (5.4), and transported separately using the Sinkhorn divergence formulation for unbalanced optimal transport.

Results

Reconstructed physical models with the L 2 and S ", based inversion are shown in The observed data and the predicted data at the last iteration for both the L 2 and S ", based inversion are reasonably in phase in the shallow region (above 1500 m).

In the predicted data at the first iteration reflected phases, as expected, are missing whereas diving phases are reasonably in phase only at short offset. This reflects into the L 2 and the S ", adjoint sources.

The S ", adjoint source is smooth and capture information of both diving and reflection phases. The source captures also signature of a mass-spreading over the whole domain as a result of the unbalanced-mass between predicted and observed data. At first iteration, the he L 2 and the S ", based inversion gradients look different with different sensitivity regions at depth.

The S ", based inversion convergence rate appears faster than the L 2 -based inversion one and the misfit distance keeps continuously decreasing after 250 iterations.

1D-GRAD initial model

Reconstructed models with the L 2 and S ", based inversion are not satisfactory after 250 iterations., even if some long-wavelength components start emerging in a more satisfactory way for S ", . It is an extremely The predicted data at last iteration does not fit well the observed data for both the L 2 and S ", based inversion. Long wavelength error components appears for L 2 possibly as a result of cycle-skipping issue.This appears less an issue for S ", .

The L 2 and S ", adjoint source are again different and does not capture the same space and time information in the early and late times of the source. The gradient of the S ", exhibit large lobe-sided, long-wavelength updates in the shallow part (above 800 m) and relatively less sensitivity in the deeper region. The gradient of the L 2 based inversion exhibits more complicated and stronger updates in the region between 800-2000 m depth.

The convergence rate, over the first 250 iterations, appears faster with the S ", based inversion compared to the L 2 based inversion, and the misfit keeps continuing decreasing. Chapter 7

Conclusion

This thesis aims at investigating application of entropic unbalanced optimal transport in the context of fullwaveform seismology. We provide a brief overview of full-waveform inversion in seismology and of the main issues encountered in practice and that are associated to the inherent ill-posedness and the non-convexity of the inversion problem, together with a brief review of different research directions that have been attempting

to resolve these issues.

The premise of this study is that transport based distances can be a reasonably effective tool, partly because of their Lagrangian nature, for designing objective functions alternative to classical L p based misfit functions that improve convexity and enlarge the basin of attraction of the full-waveform inversion problem. This has been already promoted by a number of studies investigating transport distances, such as 2-Wasserstein distance, earth mover's distance and Graph optimal distance. This Thesis can be seen as a contribution to this new trend.

More specifically in Chapter 4, we present the Sinkhorn divergence formulation of entropic unbalanced optimal transport, correcting the bias introduced by the entropic penalisation, that benefits from the easy The numerical illustrations provided in this Thesis make however a number of assumptions, for sake of simplicity and computing cost, that need to be kept in mind when reaching conclusions. First they all assume constant density, acoustic wave propagation (only one model parameter) and as such mode conversion and anisotropy are not considered. Second, all observables are point samples solution of the acoustic wave equation, using the same solver and without additional noise such as acquisition noise (full control set-up, so called inversion crime setting). Third seismic sources are assumed to be known and modelled as point sources with a Ricker wavelet source-time function.

As such, the numerical illustrations provided in this Thesis have to be considered as academic illustrations that need to be further investigated and evaluated in more realistic contexts before drawing any practical conclusion in term of field data applications.

At the end of this document, we cannot unfortunately provide a definitive conclusion on the relevance of using Entropic Unbalanced OT in FWI even just in the context of academic models. We here give where we stand and the various ideas arising from this work.

Based on chapter 6 we can say, as expected, that S ✏, provides similar results as the 2-Wasserstein distance promoted by the Engquist group under the use of signal sign splitting. It is yet unclear whether the introduction of unbalanced OT, avoiding signal normalisation, improves the inversion or not.

Introducing unbalanced (to avoid normalisation) and entropic (for its numerical performance and flexibility) OT is appealing but introduces two new parameters ✏ and which have so far simply been set by trial and visual assessment on the numerical tests. We do not have a clear understanding of the impact of these parameters on the convexity of the misfit and how they can be optimised. One of the conclusion is therefore that this should be first investigated in more details in the 2D parametric example (section 5.2).

The potential use of these parameters in simulated annealing strategies could also be investigated.

On the positive side, our Entropic misfit function is easy to implement and the Sinkhorn divergence formulation really allows to set ✏ such that the number of Sinkhorn iterations to convergence remains acceptable. In all our experiments, it remained on the order of the discretisation of the offset line and in practice, one OT gradient step never exceeded 10 times one L 2 gradient step and was more like 5 times on average. This is encouraging and indicates that the extension of this method to 3D FWI is possible thanks to the tensorisation technique (section 4.1.2).

Another interesting research direction raised by this study is the choice of the ground metric, the metric used to measure the distance (in time ⇥ offset) traveled by the mass (here the positive or the negative part of the signal). A first discussion just on scaling the dimensions of the common shot gather is given in section 4.4 with two unrelated conclusions: first, scaling the domain will affect the modulus of convexity for simple translating signals, and second when using entropic OT the scale of the domain will constrain the transport range and can indeed be used to control it. In the OT theory and also numerical practice, it is legit to replace The strategy developed in this thesis remains based on the awkward interpretation of seismic wavefield as probability densities and we pay the heavy price of poorly understood non-linear signal transformations even with the unbalanced OT (see above) normalisation, even though it avoids signal normalisation that often dampen the features and reduces the intensity range. In contrast with the L p distances that only considers differences in intensity, OT distances are transport distances and as such less sensitive to high frequency perturbations, the transport being on the order of the wavelength of the perturbation.

The representation of the data as a graph in the (R amplitude ⇥ R time ⇥ R d 1 offset ) space proposed by the Metivier group is more natural and gives good results. There is a numerical price to pay though and it has remain until today limited to trace-by-trace comparisons. The application of Entropic OT maybe of interest for its numerical efficiency and also the range controlling property mentioned above.

Finally let us remark that the "graph" strategy proposed above can be pushed further, e.g. a discretised seismogram line (a collection of ordered points) can be interpreted as a Dirac measure in R d 1 offset ⇥ R nt time ⇥ R amplitude and the OT distance can be defined on probability distributions over this space. This strongly reduces the number of samples of the data distributions but the computational burden now resides on the computation and storage of the ground cost. Interestingly, an efficient software exists (developed by J. Feydy https://www.kernel-operations.io/geomloss/), based on an online strategy to evaluate the ground cost and on a sophisticated GPU implementation, to compute Sinkhorn divergences in high dimensional spaces.

It may be interesting to test if OT distances based on deformations in these higher dimensional spaces can be useful.

Institut de Physique du Globe de Paris

Miao YU, Ph.D. Thesis, 2021

2. 1

 1 Flow chart of a simple full-waveform inversion workflow. ..................... 2.2 Simple example illustrating cycle-skipping artefact. The solid black line signal in the middle represents the observed sinusoidal signal with period T . The upper dashed line signal represents a predicted signal, which is simply the observed signal shifted by a time delay greater than, T/2. In this case the local optimisation problem will update the model so that the n +1-th phase of the predicted signal matches the n-th phase of the observed signal. The lower dashed line signal represents a predicted signal, which is the observed signal shifted by a time delay of less than T/2. In that case the local optimisation problem will rightly update the model such that the n-th phases of the observed and predicted signals match. From Virieux and Operto (2009). .......................................... 3.1 Kantorovich's transport plan between two discrete distributions of masses defined on spaces X and Y , allowing to distribute mass from the one location to multiple locations. ....... 4.1 KL function. s 7 ! KL(s|1) + 1. ................................... 4.2 Influence of the regularization parameter " on the convergence rate for (4.16)o fS i n k h o r n ' s algorithm (source: Peyré and Cuturi (2019)) ...........................

  (source: Benamou et al. (2015)). The coupling are plotted in grey level in x ⇥ y space. ...................................... 4.4 Colormap of C (4.37) as a function of ✓ M (the Metric Angle) and ✓ D (the Propagation Angle).

5. 2

 2 The top row are a three layer synthetic model recording (d obs )a n dd a t a( d cal ) for one of the intermediate models in the simulation. The bottom line show the S ", gradient for the positive part and the P derivative of the signal needed to smooth the gradient. ............ v 5.3 Top figure: Ricker wavelets f (t) (solid blue line) and g(t) (solid red line), the latter being a shifted version for a translation s =0 .25. Bottom Figures: the negative (bottom-left figure) and the positive (bottom-right figure) parts of the two Ricker wavelets. Also their negative and positive parts. .......................................... 5.4 Misfit with respect to the shift s. .................................. 5.5 Top figure: Ricker wavelets g(t) (solid blue line) and f (t) (solid red line), with an added Gaussian noise. Bottom figures: the negative (bottom-left figure) and the positive (bottomright figure) parts of the two Ricker wavelets. ........................... 5.6 Comparison between the L 2 (solid blue curve), the balanced (normalised) S ✏, (solid purple curve), and the unbalanced OT ✏, (solid red curve) and S ✏, (solid yellow curve) misfit functions as a function of the time shift between the two Ricker wavelets for increasing Gaussian noise added to one of the wavelet. From top-left to bottom-right N (0, 0.01), N (0, 0.05), N (0, 0.1), N (0, 0.2):1 , 5 , 1 0a n d2 0%. ..................................... 5.7 Comparison between the L 2 (solid blue curve), the balanced (normalised) S ✏, (solid purple curve), and the unbalanced OT ✏, (solid red curve) and S ✏, (solid yellow curve) misfit functions as a function of the time shift between the two Ricker wavelets (one with added 5% noise) for decreasing ✏ and . .......................................... 5.8 Top figure: Ricker wavelet g(t) (solid blue line) and Ormsby Wavelet f (t) (solid red line) with characteristic frequencies

5 .

 5 figure); FWI results with the L 2 -misfit function (bottom-left figure) and the S ", -misfit function (bottom-right figure). ...................................... 6.5 Seismic traces with the S ", -misfit function: initial (solid-green line) and predicted data (solidred line) after 25 FWI iterations, together with observed data (solid-blue line). ........ 6.6 Velocity profile along a vertical section of the physical domain at a surface position x = 600 m as reconstructed by the S ", based inversion. .......................... 6.7 Wavefield in the time-receiver domain for the common shot gather associated to the centred source: "true" wavefield (upper-right figure); initial wavefield (upper-left figure); wavefield at the final FWI iteration (lower-right figure); difference between the "true" wavefield and the wavefield at the final FWI iteration (lower-left figure). ...................... 6.8 Final reconstructed models: initial model (upper-left figure) and targeted model (upper-right figure); FWI results with the L 2 -misfit function (bottom-left figure) and the S ", -misfit function (bottom-right figure). ...................................... 6.9 Seismic traces: initial (solid-green line) and predicted data (solid-red line) after 25 FWI iterations, together with observed data (solid-blue line). ....................... 6.10 Ajoint-source in the time-receiver domain at the first FWI iteration. .............. 6.11 Wavefield in the time-receiver domain for the common shot gather associated to the centred source: "true" wavefield (upper-right figure); initial wavefield (upper-left figure); wavefield at the final FWI iteration (lower-right figure); difference between the "true" wavefield and the wavefield at the final FWI iteration (lower-left figure). ...................... 6.12 Reconstructed models for different implementation strategies of the S ", -misfit function: 2-D transport after decomposition of the signals into poistive and negative parts (upper-left figure); Squared-transform of the signals before 2-D transport (upper-left figure); 1-D transport in space (offset-by-offset) after decomposition of the signals into positive and negative parts (lower-left figure); 1-D transport in tile (trace-by-trace) (in time) after decomposition of the signals into positive and negative parts (lower-right figure). The later two implementations make use of the separability of the Gibbs matrix in the Sinkhorn algorithm and of the appropriate scaling of the transport domain. ............................ 6.13 Velocity profile along a vertical section of the physical domain for different FWI iterations and different implementations of the S ", -misfit function. ....................... 6.14 Final reconstructed models: initial model (upper-left figure) and targeted model (upper-right figure); FWI results with the L 2 -misfit function (bottom-left figure) and the S ", -misfit function (bottom-right figure). ...................................... 6.15 Seismic traces: initial (solid-green line) and predicted data (solid-red line) after 25 FWI iterations, together with observed data (solid-blue line). ....................... 6.16 Convergence rate for the S ", based inverstuin ...........................

  islinked to a graph-space optimal transport distance, following original developments in transport L p distance Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 for signal analysis by Thorpe et al. (2017). Discrete graphs of the predicted and observed signals are compared rather than the signals themselves. The time-signal is represented, after time discretisation, by a point cloud in the 2D time-amplitude domain, which ensure the positivity of the compared data.

Chapter 6 ( 2

 62 2D FWI numerical illustrations) illustrates numerically the use in the FWI context of the misfit function presented in chapter 4 on synthetic data, through simple canonical models, and discusses the results. Chapter 7 (Conclusion) gives partial conclusions supported by this study. We also describe a possible follow up work program. Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 Chapter Seismic full-waveform inversion: an overview This section provides a brief overview of the full-waveform inversion and of the adjoint-state formulation, following closely the reviews of Tarantola (1988); Plessix (2006); Virieux and Operto (2009); Fichtner (2011); Virieux et al. (2014); Chauris (2021)

. 8 )

 8 This is not a trivial task as J depends on m through d cal and the relation between d cal and m is given by the state equation (2.1)-(2.2)a n d( 2.6) and is clearly non linear. Global optimisation methods, which imply exploring the admissible model space M, still remain unfeasible today for most seismology applications owing to the large number of model parameters. In practice, the minimisation proceeds iteratively toward the global minimum m ? , starting from a plausible initial model m 0 , using local optimisation approach. Gradient-based methods Nocedal and Wright (2006), Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

  (2.34) Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

  .39) For example, given the partition Ω = L m C m , function m (x) can be 1 on C m and 0 in Ω\C m . Other choices are of course possible. 3. Function h(d cal , m, d obs ) does not depends on m = c ? :

Figure 2 . 1 :

 21 Figure 2.1: Flow chart of a simple full-waveform inversion workflow.

  Both seismic observations and predicted observables are oscillatory signals by nature. As such the leastsquares based misfit function (2.63) is oscillatory and non convex in directions associated with wavenumber Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 components of the model that are longer than seismic wavelengths, i.e. longer than a half-period of the dominant inverted wavelet. For large kinematic errors, local optimisation algorithms will match incorrect phases between the predicted and observed signals, making the method prone to the so-called cycle-skipping effect, i.e. the inversion converges in a local minimum possibly far from the global minimum due to poor starting model. This can be illustrated in Figure2.2 with a very simple example where both the observed and predicted signals are sinusoidal signals of period T .Apart from finding a good initial kinematic model, reformulating the FWI with alternative ways to measure the difference between predicted observable and observations to mitigate the non convexity with respect to time and space shifts has long been recognised in the FWI community as a proxy for the convexity of the FWI problem with respect to the model wave velocities[START_REF] Jannane | Wavelengths of earth structures that can be resolved from seismic reflection data[END_REF].

Figure 2 . 2 :

 22 Figure 2.2: Simple example illustrating cycle-skipping artefact. The solid black line signal in the middle represents the observed sinusoidal signal with period T . The upper dashed line signal represents a predicted signal, which is simply the observed signal shifted by a time delay greater than, T/2. In this case the local optimisation problem will update the model so that the n +1-th phase of the predicted signal matches the n-th phase of the observed signal. The lower dashed line signal represents a predicted signal, which is the observed signal shifted by a time delay of less than T/2. In that case the local optimisation problem will rightly update the model such that the n-th phases of the observed and predicted signals match. From Virieux and Operto (2009).

. 1 .Figure 3 . 1 :

 131 Figure 3.1: Kantorovich's transport plan between two discrete distributions of masses defined on spaces X and Y , allowing to distribute mass from the one location to multiple locations.

  observed d obs signals into a positive and a negative part and to recombine them as dcal = d + cal + d obs dobs = d + obs + d cal (3.15) Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

  with regard to the positivity assumption is to split the predicted d cal and the observed d obs signals into a positive and a negative part, and to scale them by their respective total mass. The objective function h(d cal ,d obs ) is defined as the sum of the associated Wasserstein distances h(d cal ,d obs )=W 2 d + = max {d, 0} and d = max { d, 0}.

A 4

 4 more recent strategy introduced by Métivier et al. (2018); Metivier et al. (2019) is linked to a graphspace optimal transport distance, following the original developments of transport L p distance for signal analysis by Thorpe et al. (2017). After discretisation in time, the predicted d cal and the observed d obs signals are mapped into the 2D time-amplitude graph space as point clouds that are compared rather than the signals themselves. The graph-space optimal transport distance can be efficiently computed as a linear assignment problem.This signal graph representation ensures positivity. This new objective function shows promising convexity properties with respect to translation and amplitude variation. This has been demonstrated for field data applications (Górszczyk et al., 2020, 2021; Pladys Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 et al., 2021). However extension to 2D and 3D source gathers remain to be evaluated as it might increase significantly the computational cost. Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 Chapter Entropic regularization of OT and its generalizations A new class of numerical methods based on "entropic regularisation" has been introduced for OT computations in Cuturi (2013) (see also the book Peyré and Cuturi (2019)). It penalises OT(µ, ⌫) with the Shannon entropy of the transport plan. A small parameter, often denoted as ", controls this penalisation.

Figure 4 . 1 :

 41 Figure 4.1: KL function. s 7 ! KL(s|1) + 1.

( 4

 4 .15) initialised with an arbitrary vector b 0 = m , where m is the m-vector of ones.. Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

  converges linearly, e.g. see[START_REF] Peyré | Computational optimal transport[END_REF]Remark4.13,buttheasymptotic convergence rate deteriorates when " ! 0. In the case of quadratic optimal transport (4.3)t h ea l g o r i t h m behaves as O(1 "). An example borrowed from[START_REF] Peyré | Computational optimal transport[END_REF], and reproduced in Figure4.2, illustrates this effect.

Figure 4 . 2 :

 42 Figure 4.2: Influence of the regularization parameter " on the convergence rate for (4.16)o fS i n k h o r n ' sa l g o r i t h m (source: Peyré and Cuturi (2019))

Figure 4 .

 4 Figure 4.3: (top) The input densities p (blue) and q (curve) as functions of x and y.( c e n t e r ) E v o l u t i o n o f t h e couplings l at iteration l of the Sinkhorn algorithm. (bottom) solution ? " of problem (4.4)f o rs e v e r a lv a l u e so f". (source: Benamou et al. (2015)). The coupling are plotted in grey level in x ⇥ y space.

  , a new metric on positive Radon measures has been independently introduced in Chizat et al. (2018b), Kondratyev et al. (2016b)a n dLiero et al. (2018). Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Figure 4 .

 4 Figure 4.4 shows the value of this parameter for different regimes (✓ D ,✓ M ).

Figure 4 . 4 :

 44 Figure 4.4: Colormap of C (4.37)a saf u n c t i o no f✓ M (the Metric Angle) and ✓ D (the Propagation Angle).

  4.41) Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

  solve the 1D trace by trace OT " problem. This remark naturally extends to unbalanced (soft marginal constraints) and Sinkhorn divergence variants. Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021Chapter 5

Figure 5 . 1 :

 51 Figure 5.1: Illustration of the effect of for the transformation P defined by (5.4), applied to the function f (t)= t, 8t 2 [ 1, 1], which changes sign at t =0.A s increases, P (f ) is getting smoother and smoother in the vicinity of t =0,b u ti n c r e a s i n ga r t i fi c i a lm a s si sg e t t i n gc r e a t e dt h a tl e a dt op o s s i b l eb i a si no p t i m a lt r a n s p o r t .

Figure 5 . 2 :

 52 Figure 5.2: The top row are a three layer synthetic model recording (d obs )anddata(d cal )foroneofthein termediate models in the simulation. The bottom line show the S ", gradient for the positive part and the P derivative of the signal needed to smooth the gradient.

Figure 5 . 3 :

 53 Figure 5.3: Top figure: Ricker wavelets f (t) (solid blue line) and g(t) (solid red line), the latter being a shifted version for a translation s =0 .25.B o t t o mF i g u r e s :t h en e g a t i v e( b o t t o m -l e f tfi g u r e )a n dt h ep o s i t i v e( b o t t o m -r i g h t figure) parts of the two Ricker wavelets. Also their negative and positive parts.

Figure 5 . 4 :

 54 Figure 5.4: Misfit with respect to the shift s.

Figure 5 . 5 :

 55 Figure 5.5: Top figure: Ricker wavelets g(t) (solid blue line) and f (t) (solid red line), with an added Gaussian noise. Bottom figures: the negative (bottom-left figure) and the positive (bottom-right figure) parts of the two Ricker wavelets.

Figure 5 . 6 :

 56 Figure 5.6: Comparison between the L 2 (solid blue curve), the balanced (normalised) S ✏, (solid purple curve), and the unbalanced OT ✏, (solid red curve) and S ✏, (solid yellow curve) misfit functions as a function of the time shift between the two Ricker wavelets for increasing Gaussian noise added to one of the wavelet. From top-left to bottom-right N (0, 0.01), N (0, 0.05), N (0, 0.1), N (0, 0.2):1 , 5 , 1 0a n d2 0%.

Figure 5 . 7 :

 57 Figure 5.7: Comparison between the L 2 (solid blue curve), the balanced (normalised) S ✏, (solid purple curve), and the unbalanced OT ✏, (solid red curve) and S ✏, (solid yellow curve) misfit functions as a function of the time shift between the two Ricker wavelets (one with added 5% noise) for decreasing ✏ and .

Figure 5 . 8 :

 58 Figure 5.8: Top figure: Ricker wavelet g(t) (solid blue line) and Ormsby Wavelet f (t) (solid red line) with characteristic frequencies (5,10,30,35Hz) . Bottom figures: the negative (bottom-left figure) and the positive (bottom-right figure) parts of the two wavelets.

Figure 5 . 9 :

 59 Figure 5.9: Comparison between the L 2 (solid blue curve), the balanced (normalised) S ✏, (solid purple curve), and the unbalanced OT ✏, (solid red curve) and S ✏, (solid yellow curve) misfit functions as a function of the time shift between Ricker and Ormsby wavelets for decreasing ✏ and .

Figure 5 .

 5 Figure 5.10: The L 2 misfit function as a function of the model parameters v 0 and ↵. The red star indicates the location of the global minimum, while in the right figure the intersection of the two dashed lines indicate the reference model v ? 0 ,↵ ? in the model space.

Figure 5 .

 5 Figure 5.11: The OT ", misfit function as a function of the model parameters v 0 and ↵. The red star indicates the location of the global minimum. The blue star indicates the position of the reference model v ? 0 ,↵ ? in the left the figure, while in the right figure this position is indicated as the intersection of the two dashed lines.

Figure 5 .

 5 Figure 5.12: Misfit function based on the Sinkhorn divergence (S ", ), with respect to v 0 and ↵.T h er e ds t a ri st h e location of the global minimum, while the position of the reference model v ? 0 ,↵ ? is indicated as the intersection of the two dashed lines.

  (4.30), through canonical 2D models. The first model, i.e. the Camembert model, focuses on transmission and reflection configurations. The second model, i.e. layered model, focuses on simple reflection configurations. The third model, i.e. the Marmousi model, focus on a more realistic configuration. All the numerical examples are performed using 2D constant-density acoustic wave modelling and local gradient-based optimisation as implemented in the PySIT platform. This is reviewed in section 2. The seimic source is assumed to be a point source with a Ricker wavelet for time function. For all the numerical experiments, the main parameters values associated to the S ", -misfit function and the iterative Sinkhorn algorithm are summarised below Parameter (used value) Use " = 10 2 Entropic regularisation, see section 4.1.1,2 = 10 1 KL unbalanced penalisation of the marginal constraints on d cal and d obs , see section 4.2.2. = 10 3 P -transform into positive/negative part of d cal and d obs , see section 5.1.

  n t = 128 Number of discretisation points in time used for the OT computation. This implies a sampling of the signals which discretization vary according to the model and n r . n r problem-dependent Number of discretisation points in offset for the OT computation. It also is by default the number of receivers used in the acquisition configuration and the FD discretization of the forward model. k max t = [25, 1] Transport window in time (in number of points) (see section 4.4.2), k max t =1 corresponds to the trace by trace misfit.

  Figure 6.1: Numerical inversion after 25 iterations

Figure 6 . 2 :

 62 Figure 6.2: Inversion convergence

Figure 6 .

 6 Figure 6.3: S ", inversion: transmission (initial model B)

Figure 6 . 6 ,Figure 6 . 4 :

 6664 Figure 6.6, as well as the exact background velocity. The shape of the inclusion is not fully recovered also possibly due to limited offset and relative deficit in long-wavelength components. The initial and final predicted data, are shown in the time-receiver domain Figure 6.7. A diffraction pattern can be observed due to the relatively strong velocity anomaly in the inclusion. The improved convexity of the S ", -misfit function avoids the FWI problem to be trapped in local minima and enable the convergence toward the global minimum.

Figure 6 . 5 :

 65 Figure 6.5: Seismic traces with the S ", -misfit function: initial (solid-green line) and predicted data (solid-red line) after 25 FWI iterations, together with observed data (solid-blue line).

Figure 6 . 6 :Figure 6 . 7 :Figure 6 . 8 :

 666768 Figure 6.6: Velo city profile along a vertical section of the physical domain at a surface p osition x = 600 ma s reconstructed by the S ", based inversion.

Figure 6 . 9 :

 69 Figure 6.9: Seismic traces: initial (solid-green line) and predicted data (solid-red line) after 25 FWI iterations, together with observed data (solid-blue line).

Figure 6 .

 6 Figure 6.10: Ajoint-source in the time-receiver domain at the first FWI iteration.

Figure 6 .

 6 Figure 6.11: Wavefield in the time-receiver domain for the common shot gather asso ciated to the centred source: "true" wavefield (upper-right figure); initial wavefield (upper-left figure); wavefield at the final FWI iteration (lowerright figure); difference between the "true" wavefield and the wavefield at the final FWI iteration (lower-left figure).

Figure 6 .

 6 Figure 6.13 for comparison. A first implementation makes use of the squared-transformation, without normalisation, of the observed and predicted signals. With this signal transformation, the 2D S ", -misfit function leads to reasonably correct model reconstruction, and relatively improved estimation of the peak-amplitude of the velocity anomaly compared to the 2D S ", -misfit function based on the decomposition into positive and negative parts of the signals. Another implementation makes use of the separability of the Sinkhorn kernel matrix and of appropriate scaling of the transport domain, see section 4.4.2. The S ", -misfit can be reduced to either 1-D unbalanced transport in time (trace-by-trace) or 1-D unbalanced transport in space (offset-by-offset).

Figure 6 .

 6 Figure 6.12: Reconstructed models for different implementation strategies of the S ", -misfit function: 2-D transport after decomposition of the signals into poistive and negative parts (upper-left figure); Squared-transform of the signals before 2-D transport (upper-left figure); 1-D transport in space (offset-by-offset) after decomposition of the signals into positive and negative parts (lower-left figure); 1-D transport in tile (trace-by-trace) (in time) after decomposition of the signals into positive and negative parts (lower-right figure). The later two implementations make use of the separability of the Gibbs matrix in the Sinkhorn algorithm and of the appropriate scaling of the transport domain.

6 Figure 6 . 14 :

 6614 Figure 6.14: Final reconstructed models: initial model (upper-left figure) and targeted model (upper-right figure); FWI results with the L 2 -m i s fi tf u n c t i o n( b o t t o m -l e f tfi g u r e )a n dt h eS ", -misfit function (bottom-right figure).

Figure 6 .

 6 Figure 6.15: Seismic traces: initial (solid-green line) and predicted data (solid-red line) after 25 FWI iterations, together with observed data (solid-blue line).

  . The upper interface and horizontal extension of the inclusion is meaningful, whereas the lower inclusion boundary is not properly resolved and appears as "melted" because reflections from the lower inclusion boundary have not been sufficiently recorded due to the limited offset. Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

Figure 6 .Figure 6 .

 66 Figure 6.16: Convergence rate for the S ", based inverstuin

Figure 6 . 19 :

 619 Figure 6.19: Final reconstructed models: initial model (upper-left figure) and targeted model (upper-right figure); FWI results with the L 2 -m i s fi tf u n c t i o n( b o t t o m -l e f tfi g u r e )a n dt h eS ", -misfit function (bottom-right figure).

Figure 6 .

 6 Figure 6.20: Seismic traces: initial (solid-green line) and predicted data (solid-red line) after 25 FWI iterations, together with observed data (solid-blue line).

Figure 6 .Figure 6 .

 66 Figure 6.22: Ajoint-source in the time-receiver domain at the first FWI iteration.

Figure 6 .

 6 Figure 6.24: Reconstructed models for different implementation strategies of the S ", -misfit function: 2-D transport after decomposition of the signals into poistive and negative parts (upper-left figure); Squared-transform of the signals before 2-D transport (upper-left figure); 1-D transport in space (offset-by-offset) after decomposition of the signals into positive and negative parts (lower-left figure); 1-D transport in tile (trace-by-trace) (in time) after decomposition of the signals into positive and negative parts (lower-right figure). The later two implementations make use of the separability of the Gibbs matrix in the Sinkhorn algorithm and of the appropriate scaling of the transport domain.

Figure 6 .

 6 Figure 6.25: Velo city profile along a vertical section of the physical domain for different FWI iterations and different implementations of the S ", -misfit function.

Figure 6 . 26 :

 626 Figure 6.26: Convergence rate of the L 2 (left) and S ", (right) based inversion

  Figure 6.27: Final reconstructed models: initial model (upper-left figure) and targeted model (upper-right figure); FWI results with the L 2 -m i s fi tf u n c t i o n( b o t t o m -l e f tfi g u r e )a n dt h eS ", -misfit function (bottom-right figure).

Figure 6 .

 6 Figure 6.28: Velocity profiles reconstructed with the L 2 and the S ", based inversion at receivers 67 and 84: initial (solid-green line) and predicted data (solid-red line) after 50 iterations, together with the velocity profile in the "True" model (solid-blue line).

Figure 6 .

 6 Figure 6.30: Convergence rate with the L 2 (left) and the S ", (right) based inversion.

6. 27 .

 27 Institut de Physique du Globe de ParisMiao YU, Ph.D. Thesis, 2021

Figure 6 . 31 :Figure 6 .

 6316 Figure 6.31: Adjoint source and gradient at first iteration for the L 2 and S ", misfit functions.

Figure 6 .Figure 6 . 34 :

 6634 Figure 6.33: Reconstructed models for different implementation strategies of the S ", -misfit function: 2-D transport after decomposition of the signals into poistive and negative parts (upper-left figure); Squared-transform of the signals before 2-D transport (upper-left figure); 1-D transport in space (offset-by-offset) after decomposition of the signals into positive and negative parts (lower-left figure); 1-D transport in tile (trace-by-trace) (in time) after decomposition of the signals into positive and negative parts (lower-right figure). The later two implementations make use of the separability of the Gibbs matrix in the Sinkhorn algorithm and of the appropriate scaling of the transport domain.

Figure 6 . 35 :Figure 6 . 36 :Figure 6 . 37 :

 635636637 Figure 6.35: Velo city profiles in the physical domain at surface p osition x = 1000 m for different iterations with the L 2 and the S ", based inversion.

Figure 6 .

 6 Figure 6.38: Reconstructed models for different implementation strategies of the S ", -misfit function: 2-D transport after decomposition of the signals into poistive and negative parts (upper-left figure); Squared-transform of the signals before 2-D transport (upper-left figure); 1-D transport in space (offset-by-offset) after decomposition of the signals into positive and negative parts (lower-left figure); 1-D transport in tile (trace-by-trace) (in time) after decomposition of the signals into positive and negative parts (lower-right figure). The later two implementations make use of the separability of the Gibbs matrix in the Sinkhorn algorithm and of the appropriate scaling of the transport domain.

Figure 6

 6 Figure 6.39. The 2-D S ", misfit function based on squared-amplitude signal transformation, leads to a meaningful reconstruction similar to the L 2 results in Figure 6.34. The peak-amplitude of the velocity anomaly is only partially reconstructed after 50 iterations, and similar to the one obtained with the 2-D S ", -misfit function based on signal decomposition into positive and negative parts, see Fig 6.39. The results obtained with the trace-by-trace S ", -misfit function, i.e. 1-D unbalanced transport in time, are very similar to the L 2 results, Figure 6.27. The peak-amplitude of the velocity anomaly is well retrieved Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021

[ 2 , 1 ]Figure 6 . 41 :

 21641 Figure 6.41: Final reconstructed models: initial model (upper-left figure) and targeted model (upper-right figure); FWI results with the L 2 -m i s fi tf u n c t i o n( b o t t o m -l e f tfi g u r e )a n dt h eS ", -misfit function (bottom-right figure).

Figure 6 . 42 :

 642 Figure 6.42: Velo city profiles in the physical domain at surface p osition x = 7500 mf o rd i ff e r e n ti t e r a t i o n so ft h e L 2 and the S ", based inversion

Figure 6 . 43 :

 643 Figure 6.43: Convergence rate for the L 2 (left) and the S ", (right) based inversion

Figure 6 .

 6 Figure 6.44: Seismic traces at receivers 84 and 101: initial (solid-green line) and predicted data (solid-red line) after 500 FWI iterations with the S ", -misfit function; observed data (solid-blue line)

Figure 6 .

 6 Figure 6.42,sho wsthattheS ", based inversion recovers gradually part of the velocity information below the upper third-layer interface. Both the L 2 and S ", based inversion seem to converge toward a local minimum with more artificial deep reflectors, see Figure 6.41, which could be the signature of an ill-posed inversion problem rather than a cycle-skipping problem.

Figure 6 .

 6 Figure 6.46: Wavefield in the time-receiver domain for a common shot gather asso ciated to the centred source: "true" wavefield (upper-right figure); initial wavefield (upper-left figure); inverted wavefield and difference between the "true" wavefield and inverted wavefield for L 2 (lower-left figures) and S ", (lower-right figures).

Figure 6 .

 6 Figure 6.47: Reconstructed models for different implementation strategies of the S ", -misfit function: 2-D transport after decomposition of the signals into poistive and negative parts (upper-left figure); Squared-transform of the signals before 2-D transport (upper-left figure); 1-D transport in space (offset-by-offset) after decomposition of the signals into positive and negative parts (lower-left figure); 1-D transport in tile (trace-by-trace) (in time) after decomposition of the signals into positive and negative parts (lower-right figure). The later two implementations make use of the separability of the Gibbs matrix in the Sinkhorn algorithm and of the appropriate scaling of the transport domain.

  1-D unbalanced transport in space, are meaningful and similar to the 2-D S ", results based on the signal decomposition into positive and negative parts. The inversion gradually recovers also the velocity below the third-layer interface. In this configuration, i.e. homogeneous, horizontally-layered medium, 1-D unbalanced transport in space retrieves velocity information associated to head waves and reflection waves at the third-layer interface. The peakamplitude of the velocity anomaly exceeds slightly the one retrived by the 2D S ", based inversion. Results obtained with the trace-by-trace S ", misfit function, i.e. 1-D unbalanced transport in time, are Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 also meaningful in this case, and the inversion gradually recovers the velocity below the third-layer interface.

Figure 6 .

 6 Figure 6.48: Velo city profile along a vertical section of the physical domain for different FWI iterations and different implementations of the S ", -misfit function.

Figure 6 .Figure 6 . 50 :Figure 6 .

 66506 Figure 6.49: Mini-square Marmousi I: initial (upper-left figure) and true (upper-right figure) models.

  5 km), with accurately resolved fine layers and normal faults, but have more trouble resolving a high resolution and accurate velocity model in the deep region because reflections from the dipping reflectors and the anticline structure have not been sufficiently recorded due to the limited offset. Pseudo velocity logs at two surface locations reconstructed with the L 2 and the S ", based inversion are compared in Figure 6.54. Velocity models in both inversion recover reasonably well the true velocity model above 1500 m, whereas quality of the inverted models degrades with depth. Both inversion fail to update at depth (below 2000 m) the low wavenumber structure of the velocity model and place reflectors at wrong Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 positions.

Figure 6 . 52 :Figure 6 .

 6526 Figure 6.52: Adjoint source and gradient at first iteration for the L 2 and S ", misfit functions

Figure 6 . 54 :

 654 Figure 6.54: Pseudo velocity logs at two surface locations x = 6480 m (first row) and x = 7380 m (second row) for the L 2 (left) and S ", based inversion

Figure 6 Figure 6 .Figure 6

 666 Figure 6.55: Fracture zone Marmousi I I mo del: results for L 2 (second row) and S ", (third row) misfit functions, with the S-500 (second column) and 1D-GRAD (third column) initial models.

Figure 6 .

 6 Figure 6.56: S-500 initial model: pseudo velocity logs at two surface locations x = 10730 m (first row) and x = 11430 m (second row) for the L 2 (left) and S ", based inversion

Figure 6 .

 6 Figure 6.57: S-500 initial model: wavefield in the time-receiver domain for a common shot gather associated to a source et the middle of the domain: observed wavefield (upper-right figure); initial wavefield (upper-left figure); final wavefield for the inverted model and difference between the observed and final wavefield for L 2 (lower-left figures) and S ", (lower-right figures).

Figure 6 .

 6 Figure 6.59: S-500 initial model: FWI convergence with the L 2 (left) and the S ", (right) misfit functions

Figure 6 .

 6 Figure 6.60: Reconstructed models for different implementation strategies of the S ", -misfit function: 2-D transport after decomposition of the signals into poistive and negative parts (upper-left figure); Squared-transform of the signals before 2-D transport (upper-left figure); 1-D transport in space (offset-by-offset) after decomposition of the signals into positive and negative parts (lower-left figure); 1-D transport in tile (trace-by-trace) (in time) after decomposition of the signals into positive and negative parts (lower-right figure). The later two implementations make use of the separability of the Gibbs matrix in the Sinkhorn algorithm and of the appropriate scaling of the transport domain.

  Figure 6.63.

Figure 6 .Figure 6 .

 66 Figure 6.61: 1D-GRAD initial model: FWI convergence with the L 2 (left) and the S ", (right) misfit functions.

Figure 6 .

 6 Figure 6.64: 1D-GRAD initial model: wavefield in the time-receiver domain for a common shot gather associated to asourceetthemiddleofthedomain: observ edw a v efield(upper-righ tfigure);initialw a v efield(upper-leftfigure);final wavefield for the inverted model and difference between the observed and final wavefield for L 2 (lower-left figures) and S ", (lower-right figures).

  to use and cost effective Sinkhorn algorithm associated to the entropic penalisation, and specifically discuss its use in the context of full-waveform inversion. The S ", distance still interprets signals as positive measures, requiring ad-hoc signal transformation into positive and negative parts (discussed in chapter 5) but unbalanced transport avoids ad-hoc normalisation methods that can often dampen features and reduce the intensity range of a signal. This new misfit function is assessed in the most generic and reproducible way through several synthetic cases (presented in chapter 5 and 6) from 1D toy models (focussing on shifted patterns) to 2D FWI numerical illustrations on simple canonical configurations in transmission and reflection, and on more realistic heterogeneous configurations extracted from Marmousi models. Comparison between misfit functions based on L 2 and different formulations of S "," , associated to different signal transformations and ground cost metrics linked to transport domain scaling, is made simple by the adjoint state formalism of the FWI (presented in Chapter 2). For reproducibility, all misfit functions and numerical tests developed in this Thesis have been implemented within the Python Seismic Imaging Toolbox (PySIT), an open-source platform developed by Russel J. Hewett and Laurent Demanet in the Imaging and Computing group of the Department of Mathematics at MIT.

  the quadratic euclidean metric by a Riemanian or even a Finsler metric. The open question is wether such a metric depending on the model could have better convexity properties and can this new non-linearity in Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021the inverse problem be tackled theoretically and numerically.

•

  Une première famille, proposée à l'origine par Engquist et ses collaborateurs est liée à la métrique 2-Wasserstein. Elle est limitée à des mesures de probabilité. Les signaux sismiques nécessitent alors une Une deuxième famille, proposée par Métivier et collaborateurs fait appel à la métrique 1-Wasserstein, instance particulière de la métrique de Kantorowitch-Rubinstein, associée à la formulation duale du problème de transport optimal. Cette métrique OT est désormais définie pour des mesures signées et peut être utilisée pour gérer un transport optimal déséquilibré, évitant la transformation et la normalisation ad hoc du signal. Cela conduit à un problème d'optimisation convexe sous contraintes linéaires qui peut être résolu en utilisant des méthodes d'optimisation convexe non lisses qui permettent de réduire significativement le coût de calcul par rapport à la métrique 2-Wasserstein. La convexité de la fonction objectif n'est cependant pas garantie vis-à-vis des transformations en temps et en amplitude.• Une troisième famille récemment introduite par Métivier et collaborateurs est liée à une distance de transport optimale dans l'espace étendu du graphe des données. Faisant ainsi suite aux développements

transformation et une normalisation ad-hoc pour être comparés à l'aide de cette métrique. Cela peut altérer les informations de phase et d'amplitude, et favoriser la non-convexité par rapport aux décalages temporels. Le calcul de la métrique 2-Wasserstein en utilisant la formulation de Monge-Ampère soit avec des différences finies optimisées, des solveurs linéaires, ou une stratégie semi-discrète ajoute un coût de calcul significatif à le FWI. Dans la pratique l'application de la distance 2-Wasserstein a souvent Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 été restreinte au transport 1D-optimal, c'est-à-dire à la comparaison trace par trace moins couteuse. • originaux de la distance de transport Lp pour l'analyse du signal par Thorpe et collaborateurs. Des graphes discrets des signaux prédits et observés sont comparés plutôt que les signaux eux-mêmes. Le signal temporel est représenté, après discrètisation temporelle, par une mesure empirique dans le domaine 2D temps-amplitude, ce qui assure la positivité des données comparées. Cette nouvelle fonction coût montre des propriétés de convexité prometteuses en ce qui concerne les transformations en temps et en l'amplitude. En 2D, la distance peut être calculée avec un coût supplémentaire modéré en tant qu'un problème d'affectation linéaire. Cependant, l'extension à des dimensions de domaine plus élevées reste à développer et à évaluer, et cela pourrait augmenter considérablement le coût de calcul.

  Nous utilisons l'approximation acoustique classique des ondes élastiques et supposons une densité constante dans les illustrations 2D ; dans ce cas, la vitesse du son est le seul paramètre du modèle.Dans l'approximation acoustique, la conversion de mode n'est pas prise en compte et le milieu est isotrope. En pratique, lorsqu'il s'agit de milieux hétérogènes, une telle approximation peut être remise en cause[START_REF] Cance | Validity of the acoustic approximation for elastic waves in heterogeneous media[END_REF]

Sinkhorn, qui réduit significativement le coût de calcul. La divergence de Sinkhorn peut être étendue à un OT déséquilibré. Cette approche nécessite une separation du signal en partie positive et partie negative ainsi qu'un filtrage. Nous donnons dans ce document un aperçu du cadre théorique de ces deux extensions récentes de l'OT conduisant à des fonctions objectifs lisses et différentiables dans le contexte de FWI, dans le contexte général de la théorie du transport optimal, ainsi qu'une revue des différentes approches existantes. Ces extensions Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 récentes sont en particulier analysées dans le contexte d'une application ) la méthode de FWI au travers d'études paramétriques 1D et 2D. Leur implémentation, exploitant le.s propriétés de tensorisation pour des espaces discrets cartésiens sont présentées ainsi qu'une évaluation de la complexité de l'algorithme de Sinkohrn.

Nous illustrons numériquement l'utilisation de ces fonctions coût dans le cadre de FWI au travers de cas académiques et d'un example de référence (Marmousi-I). Les illustrations numériques font cependant un certain nombre d'hypothèses physiques par souci de simplicité qui doivent être gardées à l'esprit :

•

  [START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport[END_REF] Yang and Engquist , 2018; Yang et al., 2018;[START_REF] Sun | The application of an optimal transport to a preconditionned data matching function for robust waveform inversion[END_REF]

	new objective or fidelity functions. Those newly proposed methods transform the local, sample-by-
	sample least-squares measure to a global one, trace by trace, or even source gather by source gather. OT-based
	distances are appealing because of their Lagrangian nature and convexity with respect to translations and
	dilations of a prescribed probability distribution. To date several families of OT-based functions have been
	proposed:
	• A first family, proposed originally by Engquist and collaborators (Engquist and Froese, 2014; Engquist
	et al., 2016; Engquist and Yang, 2019; Engquist et al., 2020; Engquist and Yang, 2021) is linked to the 2-Wasserstein
	metric associated to balanced optimal transport, which shows nice convex properties with respect to
	time shifts and rotations. The classical 2-Wasserstein distance however only applies to probability
	measures. Raw oscillatory seismic signals require ad-hoc transformation and normalisation

  •i denote the L 2 scalar product on R N ⇥N ,( 3.9) can be written in short form, This is a linear optimisation problem with N 2 unknowns and 2 ⇥ N linear constraints. The best known linear programming methods have cubic complexity, see remark 3.3 in[START_REF] Peyré | Computational optimal transport[END_REF], and cannot be practically used for realistic seismic applications.

	OT(µ, ⌫)= min 2Π(µ,⌫)	h , ci.	(3.10)

de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 Let h•,

  To mitigate this effect, it has been proposed to add in (3.18)a n d( 3.19) the same transformation applied d cal and d obs , leading to symmetric transformations.Other type of transformations originally proposed in the context of FWI least-squares objective function by[START_REF] Liu | The normalised integration method -and alternative to full waveform inversion?[END_REF][START_REF] Donno | Estimating the background velocity model with the normalised integration method[END_REF] are of relevance for optimal transport. They consider the normalised cumulative distributions associated to positive transforms of d cal and d obs ,

	Q cal (x r ,t)=	R t 0 dcalc (x r ,t) dt R T 0 dcalc (x r ,t) dt	,Q obs (x r ,t)=	R t 0 dobs (x r ,t) dt R T 0 dobs (x r ,t) dt	(3.20)
	where		8
			> > > >	|d| ,		[absolute value of the signal] ,
			>
	d =	<	d 2 ,		[square of the signal] ,	(3.21)
			>
			>
			> > > :	E(d),	[enveloppe of the signal]
	second class of transformation is the exponential transformation
		dcal =	exp(bd cal ) hexp(bd cal )i	, dobs =	exp(bd obs ) hexp(bd obs )i	,b > 0	(3.18)
	or the more stable Softplus version of the scaling,
	dcal =	log(exp(b dcal ) + 1) hlog(exp(b dcal ) + 1)i	, dcal =	log(exp(b dcal ) + 1) hlog(exp(b dcal ) + 1)i	,b > 0	(3.19)

Under such transformations the influence of the high-positive values and of the high-negative values are not Institut de Physique du Globe de Paris Miao YU, Ph.D. Thesis, 2021 the same.

Other strategies for full-waveform inversion

Seemingly uncorrelated with the previously mentioned strategies, Métivier et al

(Métivier et al., 2016a,b; Métivier et al., 2016)

developed an approach based on a generalised dual W 1 distance (3.14), see Section 3.2.2.

  1/d where N is the total number of points in time ⇥ offset. Then the cost of one Sinkhorn iteration is O(N 1+1/d ) in place of O(N 2 ). Of course one still has to iterate Sinkhorn until convergence.

  2 2 quadratic Wasserstein distance, we can use the unbalanced OT modification and its various entropic approximations. A first choice is to use (4.23) h OT ", (d cal , d obs )=OT ", d + cal , d + obs ) + OT ", d cal , d obs (5.2) but we will mostly focus on its "debiased" Sinkhorn divergence version (4.30) h S ", (d cal , d obs )=S ", d + cal , d + obs + S ", d cal , d obs .

Institut de Physique du Globe de ParisMiao YU, Ph.D. Thesis, 2021

Institut de Physique du Globe de Paris

Miao YU, Ph.D. Thesis, 2021

in time with constant time step satisfying the CFL condition as implemented in PySIT.

The inversion problem is solved using the l-BFGS gradient-based method, as implemented in PySIT, with a fixed number of iterations set to 500 iterations. The initial model is a model where the deepest third layer is unknown, see Figure 6.40.

Both the predicted d cal and the observed d obs data, i.e. the latter corresponding to the "true" model, are physical realisations obtained with the same wave propagation solver, and without considering additional noise to the data, i.e. the inversion crime setting.

Unless otherwise specified, the S ", -misfit function is formulated as in Chapter 5. The observed and predicted signals are decomposed into positive and negative parts, with the P -transform defined in (5.4), and transported separately using the Sinkhorn divergence formulation for unbalanced optimal transport. In this model configuration, there is no back-scattered information from the interior of the third layer returning to the receivers.