La thèse est une longue aventure, semblable à la descente d'un long fleuve, des sources perdues dans les montagnes jusqu'à son immense embouchure. C'est un voyage duquel on ressort grandi, une aventure qui nous forme, nous façonne et modifie profondément notre perception du monde. Ce fleuve, aux eaux parfois paisibles, mais aussi tumultueuses par moment, requiert toute notre force mais aussi celle de notre entourage. Je n'aurais pas facilement navigué jusqu'à l'embouchure que représente la soutenance sans l'aide et le soutien de nombreuses personnes que je tiens à remercier ici.

Tout d'abord, je voudrais chaleureusement remercier Romain et Ludovic, mes deux directeurs de thèse. Vous qui m'avez découvert et suivi depuis le master, qui m'avez fait découvrir le monde de la recherche en profondeur durant mon stage de M2 et m'avez donné le goût de continuer à approfondir un petit peu le savoir humain au travers de cette thèse. Vous avez été d'un grand soutien, toujours présents, même lors des coup durs, et toujours prêts à trouver des solutions. Travailler avec vous a été une immense chance, et faire partie du groupe SEIS-COPE fut une expérience forte que peu de doctorants ont la chance de vivre. Ce sentiment d'appartenir à une équipe, de pouvoir collaborer avec d'autres personnes activement, de partir ensemble en conférence à l'autre bout du monde. Ce sont des moments que je n'oublierai pas.

C'est donc tout naturellement que je voudrais remercier les collègues de SEISCOPE avec qui j'ai vécu de formidables aventures. Merci Philippe Le Bouteiller pour les incroyables moments passés ensemble et pour notre amour des pots ! Un grand merci au reste de l'équipe pour tous les bons moments passés ensemble et les souvenirs partagés ! Merci aux membres de mon Jury, Hervé Chauris, René-Edouard Plessix, Philippe Roux, Jean-Pierre Vilotte et Jérémie Messud. Pour une grande partie d'entre vous, j'avais déjà eu le plaisir de vous rencontrer lors de conférences ou du meeting annuel de SEISCOPE. C'était un grand honneur de vous avoir dans mon jury et l'intérêt que vous avez porté à mon travail me touche particulièrement. Je vous remercie aussi pour la très enrichissante discussion que nous avons pu avoir ensemble lors de ma soutenance.

Merci aux collègues d'ISTerre, avec une pensée toute particulière pour Jacques Brives, toujours là pour discuter, rigoler, et avec qui j'ai passé de supers moments. Merci pour ton soutien. Merci aussi aux amis de licence, Jonathan Sekkat et William Moraglia. Vous avez tous les deux vécu des moments difficiles, et j'espère que le futur sera radieux pour vous. Vous avez été de magnifiques compagnons de route durant ces années universitaires et je ne vous oublierai pas.

Merci à Guillaume et Anne, rencontrés lors de la licence à Nice, et qui sont devenus au fil de temps des amis très proches. Merci pour les bons moments passés ensemble, d'être venus en Norvège me voir lors de mon séjour chez Equinor. Tout cela restera d'incroyables souvenirs pour la vie, et il y en aura encore plein d'autres à écrire j'espère ! Un merci tout particulier pour Grégoire Guillet, rencontré lors du master. Tu as été incroyable, et les moments passés ensemble, à créer les pires playlists de mauvais goût sur Spotify, à réussir à me motiver pour aller skier (et par la même occasion avoir les pires courbatures de ma vie !), mais aussi à me proposer un masque de ski pour couper des oignons puis de passer à deux doigts de mettre le feu au chalet. Même si on n'est souvent pas rassuré , je te souhaite le meilleur pour le futur et tu resteras toujours un ami exceptionnel.

iii Enfin, et j'ai envie de dire évidemment , un immense merci à Julien Thurin. On s'est rencontré en L1 à Nice, mais c'est à partir de la L2 et de notre parcours géophysique que nous avons commencé à être si proches. On s'est suivi durant toutes ces années, de la licence, du parcours géophysique, du M1/M2, jusqu'à la thèse. C'est complètement incroyable d'avoir fini par faire tous les deux une thèse dans la même équipe avec les mêmes directeurs de thèse, après toutes ces années partagées ! Du voyage au Bénin ensemble pour une mission de terrain à planter des géophones par 40 degrés dans du sol dur comme du béton, aux vacances à la montagne ensemble avec Greg, et tous les congrès qu'on a fait ensemble durant la thèse. Ton coté hipster incroyable, et ta capacité à être instantanément à 400% sur un projet me fascineront toujours ! Tu as vécu, comme moi, des moments pas faciles durant tes années de thèse et tu as aussi réussi à affronter tout cela et revenir toujours plus fort. Tu es maintenant marié et papa d'un magnifique petit Arthur qui est tellement mignon. J'espère que ce bonheur ne cessera jamais. Toi qui es devenu au fil du temps un de mes amis les plus proches, je te souhaite le meilleur dans ta vie et j'espère qu'on gardera toujours le contact ! Un grand merci à toute ma famille, et particulièrement à ma soeur, qui ne sait toujours pas si je suis plus scientifique que cuisinier ! C'était un immense plaisir d'avoir pu vous avoir tous ensemble présent dans l'amphi après une année 2020 plus que chaotique et d'avoir réussi à partager un petit mais magnifiquement réussi pot de thèse ensemble ! Vous m'avez toujours soutenu durant ces années de fac, du premier jour de licence jusqu'à ce 14 Juin 2021 que vous, comme moi, n'êtes pas prêts d'oublier.

Un grand merci à mon papa, lui aussi docteur, mais en biologie. Merci d'avoir toujours aiguiser ma curiosité scientifique, des premiers microscopes et loupe binoculaire que tu m'avais rapportés alors que j'étais encore au collège, mais aussi de cette incroyable visite d'ATLAS au LHC qui m'avait tellement fasciné et encore plus motivé dans mon amour pour les sciences. Merci pour tes conseils et de m'avoir aiguillé vers la géologie alors que j'étais parti pour m'inscrire en biologie pour marcher dans tes pas. C'était un grand conseil, et sans toi je ne pense pas que j'aurais réussi comme ça. Mes années de fac en géophysique furent une grande révélation et j'ai vraiment adoré cette liberté et ce sentiment d'être si proche de la science. C'est cette motivation qui m'a permis d'aller jusqu'au bout et je ne sais pas si j'aurais si bien réussi dans une autre branche ! J'imagine ta fierté aujourd'hui de savoir ton fiston docteur comme toi.

Un merci tout particulier à ma maman, qui souvent était ma bouée de secours lorsque les rapides se faisaient trop intenses. Tu as toujours été derrière moi, prête à m'aider à n'importe quel moment. Tu as toujours été là, durant tout ma scolarité, à me suivre, me motiver, alors que je n'étais pas forcement le meilleur des élèves au lycée. Tu m'as toujours poussé de l'avant, donné la confiance qu'il me manquait quand je n'étais pas sûr de moi. Merci pour ton soutien indéfectible durant mes années de fac. Je sais à quel point je t'ai rendu fière à chaque fois que je réussissais mes années avec brio. Merci aussi pour tes heures de relecture méticuleuse de textes scientifiques relativement hermétiques pour les non-initiés, sacrée fin carrière à la meilleure des profs d'anglais ! Je ne te remercierai jamais assez pour tout ça. Je sais à quel point tu es fière de me savoir aujourd'hui docteur et ce diplôme sera toujours intimement lié à toi ! Et évidemment, un immense et incommensurable merci à Sophie, ma moitié qui m'a supporté (dans tous les sens du terme) pendant toute ces années ! On a traversé de nombreuses épreuves ensemble, et bien que, si par moment le fleuve semblait nous avoir coulés, nous avons réussi à ne jamais nous noyer et à toujours revenir à flots. Je ne pense pas que j'aurais réussi à arriver jusque-là sans toi, et pour cela je te serai à toujours éternellement reconnaissant.

J'espère que le futur sera radieux, et qu'on trouvera notre équilibre et que tu réussiras enfin à t'épanouir complètement, que ça soit en devenant une botaniste aguerrie collectionnant les graines ou tout autre chose qui te passionne. C'est à moi maintenant de te soutenir après tout ce que tu m'as donné.

And now, a part in English.

I would like to thank good family friends from the US, Tony and Judy Schumacher, who were always kind to me and who even woke up at 3 am to watch me live on YouTube during my defense. You were in my heart and still are, and I hope to come soon to the US to visit you.

I would also thank all the SEISCOPE team, Hugo Sanchez, Marco Salcedo, Wei Zhou, Yang Li, Paul Wellington, and all the others. You were some fantastic teammates. I would like to give some special thanks to Pengliang Yang. It was a pleasure working with you a bit on the code. I know life hasn't been that easy with you in the past years, but I hope you are now finding a nice balance. You deserve it! Another special thanks to Yubing Li, for all the good moments spent together, speaking of science but also of culture, society, and many other subjects. Your open mind was a delight, and your motivation was terrific, always ready to grab the best food and take a 30min Uber to get some good craft beer! Thanks, Andrzej Górszczyk, for the excellent discussion together, for the beers, and for your spot on scientific feedback. You are someone amazing, and I wish you all the best for the future! A special thanks to Nishant Kamath for all the time spend together, from barbecues to discussions on good craft beers. Thanks for the scientific works we did together working on the Valhall data! You also are one of a kind, and I hope the future will be bright for you and that we will keep in touch! Finally, an extraordinary thanks to an extraordinary friend. Marwan Irnaka (first of his name!), I met you during my Ph.D. as you Joined SEISCOPE for your Ph.D. Since then, we got closer and closer. You quickly join the gang in our office, and from then we became close friends. Thanks for your support in the good but also hard times, thanks for all the chit-chat in the office, thanks for the fantastic memories during conferences with you. You became one of my best friends, and I will never forget you. You are one of the nicest people I know (and maybe on Earth?), with consistently good intentions. Even if I tried to give you some French course on "how to bitch" you were always too kind for this. You are amazing and unique, and I hope the future will be bright for you. Now that you are back in Indonesia, I definitely need to come and visit you there (with maybe one master code cracked?). v Abstract Full-waveform inversion is a seismic imaging method known to yield high-resolution results that go towards directly interpretable reconstructed models. Unfortunately, it suffers from a flaw: the need for accurate initial models to converge toward meaningful reconstruction of subsurface parameters. This limitation is due to the non-convexity of the least-squares distance conventionally used as the distance measurement function, which translates into the cycle-skipping issue that has been documented since the introduction of full-waveform inversion itself. Thus, finding solutions to improve full-waveform inversion robustness to cycle-skipping has been the subject of a large number of studies. From all the propositions made to improve this method, one of the most documented is replacing the classical least-squares norm with alternative misfit functions.

We first propose a review of multiple propositions of alternative misfit functions. We explain which principles these methods are based on and illustrate how they are designed to better handle kinematic mismatch than the least-squares misfit function. A set of carefully designed synthetic benchmarking tests is then introduced to assess the behavior of a selection of alternative misfit functions containing two optimal transport-based misfit functions, a wiener filter based, an instantaneous enveloped based, and finally, a normalized integration based misfit function. On these several canonical synthetic tests, each formulation is pushed to its limits, allowing us to establish the pros and cons of each formulation.

This work can also be seen as an attempt to promote a more systematic cross-comparison of alternative misfit functions on fair benchmarking setups, as we observed that alternative misfit functions are often proposed on carefully designed synthetic setups that maximize their benefits. This observation can be related, in our view, to the discrepancy between the number of alternative formulations made and assessed on synthetic cases in the literature and actual case studies performed with them on field data.

Finally, after selecting what appears as the most promising candidates -a graph-space optimal transport-based misfit function -we compare it to the least-squares distance in a case study based on the three-dimensional ocean bottom cable data from the Valhall field. The motivation behind using this data set comes from the wide variety of studies already performed with it. This allows us to work in a relatively controlled framework while facing field data difficulties. The comparison is first performed starting from a reflection traveltime tomography initial model used in previous studies. Then, a second comparison is performed starting from a crude, linearly varying in-depth one-dimensional initial velocity model. This last setup illustrates the robustness to cycle-skipping obtained through a change of the misfit function inside the full-waveform inversion formalism. The encouraging and meaningful results obtained from this case study demonstrate that cycle-skipping robustness of full-waveform inversion can be drastically improved using the proposed graph-space optimal transport-based misfit function.

Résumé

L'inversion de forme d'onde complète est une méthode d'imagerie sismique permettant d'obtenir des résultats à haute résolution pouvant même parfois être directement interprétés. Malheureusement, cette méthode souffre d'un défaut : la nécessité d'avoir un modèle initial suffisamment précis pour assurer la convergence de l'inversion vers une reconstruction réaliste du sous-sol. Cette limitation est due à la non-convexité de la norme euclidienne qui est traditionnellement utilisée comme fonction de mesure de distance. Cette non-convexité se traduit par le problème du saut ou ambiguité de phase, qui a été documenté depuis l'introduction même de l'inversion de forme d'onde complète. C'est donc naturellement qu'un grand nombre d'études ont été réalisées afin de trouver des solutions pour améliorer la robustesse au problème du saut de phase. Parmi toutes les propositions faites pour améliorer cette méthode, une des plus documentées est le remplacement de la fonction de mesure de distance par d'autres fonctions alternatives.

Nous proposons d'abord une analyse de différentes alternatives pour remplacer la fonction de mesure de distance traditionnellement utilisée. Les principes fondamentaux sur lesquels sont basées ces alternatives sont introduits, en association avec une analyse de leurs avantages en terme de robustesse au saut de phase. Une sélection de certaines de ces fonctions coût alternatives est ensuite testée sur différents tests synthétiques soigneusement conçus pour évaluer leurs comportements, poussant chacune des différentes formulations dans ses retranchements. Les avantages et inconvénients de chacune des propositions sont établis grâce à l'analyse effectuée sur ces différents tests. Une évaluation finale est effectuée afin de définir les formulations les plus prometteuses pour une application sur données réelles.

Ce travail sur les fonctions coût alternatives s'inscrit dans la continuité d'une observation que nous avons réalisée : la disparité entre le grand nombre de fonctions coût alternatives proposées dans la litterature et le nombre reduit d'applications sur données réelles avec de telles alternatives. En effet, trop souvent les fonctions coût alternatives sont proposées sur des tests synthétiques soigneusement conçus pour mettre en avant leurs bénéfices de façon "fictive". Ce travail prône une comparaison juste et reproductible des fonctions coût alternatives.

Finalement, après avoir sélectionné la fonction coût qui nous paraissait la plus prometteusele transport optimal du graphe -nous la comparons à la norme euclidienne dans une application sur le jeu de données réelles de fond de mer du champ de Valhall. Le choix de ce jeu de données est motivé par l'abondante littérature disponible, avec de nombreuses inversions de forme d'onde complète déjà réalisées. Cela nous permet d'avoir plus de contrôle dans notre analyse. La comparaison entre les deux fonctions coût est d'abord effectuée en partant du modèle initial provenant d'une tomographie en réflexion, qui est traditionnellement utilisé avec ces données. Puis, une deuxième comparaison est effectuée à partir d'un modèle initial mono-dimensionnel variant linéairement avec la profondeur. Ce deuxième cas sert à illustrer l'amélioration de la robustesse au saut de phase obtenue grâce au transport optimal du graphe. Les résultats présentés dans cette étude sont très encourageants et démontrent bien que l'inversion de forme d'onde complète peut être rendue plus résistante au saut de phase grâce à un changement de fonction coût.

General Introduction

Studying the Earth's interior is a complex task, as only the surface is directly observable. For a long time, the only information geologists could gather came from surface observations, where only geological outcrop (visible exposure of bedrock or ancient surface deposits) can be directly observed. Using these observations and some fundamental principles of geology (uniformitarianism, original horizontality, superposition, cross-cutting relationships and Walther's law), reconstruction of geological structures has been made possible, such as mountain chain formation. Adding to these surface observations, the first few meters to kilometers could be sparsely explored through drilling. This sparse information can help constrain the reconstruction of subsurface structures made from surface observations. Still, complex areas, or locations with no outcrop, stayed out of reach, and understanding of such locations remains poor.

To gain knowledge where the eyes cannot see, we rely on one significant branch of physics: geophysics -the physics of the Earth -a vast domain of study that includes the physics of space, atmosphere, oceans, and the interior of the planet. Here, we narrow down our focus on the Earth's interior, to what the layperson thinks when he eared "geophysics". Figure 1: (a) Some of the possible ray paths for seismic waves passing through the Earth (from Fowler et al., 1990). (b) Jeffreys-Bullen (J-B) travel-times curves for an earthquake focus at the surface. Travel-times are represented as a function of epicentral distance. Each line represents a given travel path through the Earth's interior (after [START_REF] Bullen | An introduction to the theory of seismology[END_REF].

The use of seismic waves or mechanical waves presents multiple advantages. First, they are accurately observable and recorded at the surface of the Earth. Second, they can be generated from the largest and most powerful events on Earth: earthquakes, and third, they contain one of the shortest wavelengths of "geophyiscal waves" [START_REF] Lay | Modern global seismology[END_REF], apart from the electromagnetic waves, which are shorter but also present way shorter penetration depth. This last point could potentially lead to high-resolution reconstruction.

INTRODUCTION

Going high resolution with full-waveform inversion

Since seismic waves traveling through the Earth have relatively small wavelengths, the idea to extract the maximum resolution out of them is tempting. In the '80s, a high-resolution seismic imaging method with the potential resolution of up to half the propagated wavelength was proposed by [START_REF] Lailly | The seismic problem as a sequence of before-stack migrations[END_REF] and Tarantola (1984a). This method, referred to as full-waveform inversion (FWI), has been proposed as an attempt to "bridge the gap" between high wavenumber migration imaging principle [START_REF] Claerbout | Towards a unified theory of reflector mapping[END_REF] and low wavenumber approach of traveltime tomography [START_REF] Aki | Three-dimensional seismic velocity anomalies in the crust and upper-mantle under the usgs, california seismic array[END_REF]. is not dense enough. Despite this, excellent results are obtained using FWI at global scale, as shown in Figure 5 with shear-wave-speed reconstruction results on the last GLAD-M25 global FWI model. This model displays several notable features such as hotspots and subduction slabs. In recent years, strong development of passive seismic has been made [START_REF] Stehly | A study of the seismic noise from its longrange correlation properties[END_REF][START_REF] Stehly | Tomography of the Alpine region from observations of seimic ambient noise[END_REF][START_REF] Lu | High-resolution surface wave tomography of the european crust and uppermost mantle from ambient seismic noise[END_REF][START_REF] Lu | Imaging Alpine crust using ambient noise wave-equation tomography[END_REF]. These methods are based on ambient noise correlations which use recordings of noise generated by wide spread multiple sources (such as ocean waves, wind, human activity, micro-seismicity) to retrieve Green's functions between receivers couples. These Green's functions can be translated into "virtual sources" that can be used for FWI.

All of these examples confirm that FWI is now a mature seismic technic that is able, when it can achieve convergence, to yield breathtaking results that could be interpreted directly or results in drastically improved migrated images. [START_REF] Lei | Global adjoint tomography-model glad-m25[END_REF]. Remarkable hotspots are highlighted with white arrow.

INTRODUCTION

Non convexity of full waveform inversion

Still, FWI suffers from a critical flaw in its conventional formulation: the non-convexity of the least-squares norm used as a distance measurement function. This limitation has been documented since the beginning of the FWI formalism and results in the so-called cycle-skipping issue. In a simple way, cycle-skipping occurs when calculated data are shifted by more than half a period compared to the observed data. This results in a condition to apply FWI: start from a good enough initial model that predicts the data within half a period. If this requirement is satisfied, FWI should converge toward the correct solution. If not, FWI suffers from cycleskipping and is led toward a local minimum.

Practically, solutions have been proposed to mitigate the cycle-skipping issue, and the vast majority of FWI applications rely on three preliminary steps to ensure successful results at exploration scales. The first step is to pick travel time from the data. This is a labor-intensive task requiring human supervision, which is very time-consuming. From these picked travel times, the second step is to build an initial model using traveltime tomography. Ideally, the model should be as high-resolution as possible, for example, using stereotomography methods [START_REF] Billette | Velocity macro-model estimation from seismic reflection data by stereotomography[END_REF]. The last step is to build a multi-scale workflow for FWI, which follows a frequency-continuation principle that results in inverting first for the lowest frequencies available and then progressively introducing higher frequencies.

The requirement for initial models generated by tomography clearly illustrates the famous "gap" in resolution between tomography and migration. A possible solution would be to make travel-time tomography methods high-resolution. Unfortunately, getting higher wavenumber content from traveltime tomography is not an easy task, and the methods suffer from intrinsic limitation. Ideally, the solution would be to drag FWI toward low wavenumber content or even entirely remove the need for initial models generated with traveltime tomography and only relies on FWI. This motivates the idea to reformulate FWI to make it more flexible, and "automatic", thus making it less dependent on the initial model quality. By changing its formulation and reducing its non-convexity, the range of applications tackled with FWI can be increased while reducing the uncertainties associated with complex workflows.

In an attempt to make FWI more flexible, two main families can be identified. The first one relies on extension strategies. Here the goal is to artificially extend the research space to mitigate the non-convexity, then progressively constrain the artificial degrees of freedom to converge toward a non-relaxed formulation of FWI. The second one relies on replacing the non-convex least-squares-based misfit function with alternative misfit functions that present increased convexity. In the second family, many propositions have been introduced by the FWI community.

In this thesis, the focus is made on the second family. We first want to systematically compare alternative misfit functions coming from already published methods or newly proposed solutions. It is crucial to keep in mind that replacing the least-squares norm is not an easy task as it often translates into a loss of resolution power, increased complexity, and increased computational cost. This has made the least-squares norm the most used misfit function even if it suffers from a fatal flaw. To test if alternative misfit functions can improve FWI robustness, we have selected a panel of five methods representing different ways of reformulating the problem. One of the scientific contributions of this thesis is the benchmarking of alternative INTRODUCTION misfit functions in a systematic way, allowing us to exhibit their pros and cons.

Finally, we perform in-depth testing on one of the most promising solutions, a graph-space optimal transport-based misfit function, on a field dataset starting from a 1D starting model that is not based on prior inversion results. The second scientific contribution of this thesis showed that robust FWI by misfit function replacement is possible and provides encouraging results and that reformulating the original FWI problem is the way forward.

Outline of the manusctipt

Chapter 1 introduces FWI and the history behind it. An overview of the formulation of FWI is performed with details on both the forward problem and the inverse problem. Details are given on how to solve the inverse problem of FWI, with an overview of the different approaches available and the introduction of the limitation of the most commonly used method.

Chapter 2 comes from a published paper in Geophysics journal [START_REF] Pladys | On cycle-skipping and misfit functions modification for full-wave inversion: comparison of five recent approaches[END_REF]. It focuses on a subset of alternative misfit functions and presents several synthetic benchmarks that are carefully designed to assess the capacity of a given misfit function. An in-depth analysis of the pros and cons of each formulation is performed. The conclusion is drawn on the most promising methods.

Chapter 3 comes from a submitted paper in Geophysics journal [START_REF] Pladys | Robust fwi with graph space optimal transport: application to 3d obc valhall data[END_REF]. We take one of the best candidates of the previously tested misfit function, graph-space optimal transport, and use it on a field data application. This application is performed in two parts. First, using the initial model already used in the literature to assess graph-space optimal transport capacity to perform similarly to the conventional L 2 norm in a controlled environment. Second, starting from a crude 1D initial model to introduce strong cycle-skipping and assess the capacity of graph-space optimal transport in one of the most challenging cases.

Chapter 4 presents some explorations I performed during my years of research, with some ideas of new alternative misfit functions and some workarounds developed for practical application.

Chapter 5 finally gives conclusion and perspectives of this manuscript.

Chapter 1

Full-waveform inversion

A rapid overview of full-waveform inversion development

A rapid overview of full-waveform inversion development

In 1985, Claerbout's analysis on seismic imaging methods shown there was a "gap" between migration imaging principle and travel-time tomographic imaging principle (Figure 1.1). The first family, relying on the migration principle, provides high-wavenumber content of the subsurface (which translates into small sharp details) but requires a good model of the Earth's interior to provide accurate results. The second family, travel-time tomography, focuses on reconstructing a low-wavenumber model containing the overall kinematic of the Earth's interior that is not easily interpretable. As an attempt to bridge the "gap" between these two families, full-waveform inversion (FWI) was first formulated by [START_REF] Lailly | The seismic problem as a sequence of before-stack migrations[END_REF] and [START_REF] Tarantola | Linearized inversion of seismic reflection data[END_REF].

FWI is a seismic imaging method aiming at reconstructing high-resolution models (up to half the shortest wavelength) of the mechanical properties of the subsurface [START_REF] Devaney | Geophysical diffraction tomography[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain, part II: Fault delineation in sediments using crosshole data[END_REF][START_REF] Plessix | Full waveform inversion of a deep water ocean bottom seismometer dataset[END_REF]. FWI relies on a fundamental principle: wavefields traveling through the Earth's interior gather pieces of information along their path on medium physical parameters (such as P-wave velocity, S-wave velocity, density, attenuation or anisotropy parameters). Therefore, the idea behind FWI is simple: match calculated seismic data (computed by solving the wave equation numerically) to observed seismic data by iteratively updating a model of the Earth's subsurface. This type of problem is called an "inverse problem" as it boils down to estimate the parameters of partial derivative equations instead of estimating its solution. This problem can also be seen as a data-fitting procedure.

As mentioned, the main interest of FWI is its high-resolution potential that could help filling the "gap" in seismic imaging methods. We can analyze the wavenumber content of the gradient of the misfit function used to update the model, for a given couple of source/receiver of a fixed diffraction point in a homogeneous medium, in the plane wave approximation (Figure 1.2):

k = 2ω c cos θ 2 , (1.1)
where k is the wavenumber modulus, ω the angular frequency of the monochromatic plane waves, c the velocity, and θ the angle between the source and receiver rays to the diffraction point. The dependency between the wavenumber content and the illumination angle is made clear, with maximum wavenumber for zero degrees angle, meaning a reflection with zero offsets between source and receiver rays (from [START_REF] Wu | Diffraction tomography and multisource holography applied to seismic imaging[END_REF]. This illustrates that FWI resolution is constrained by three main factors. First, the most straightforward one is the frequency content of the data used. Second, the wavenumber content at depth is lower as velocity is generally higher (which translates into larger wavelengths). Third, structures at depth are often constrained by small scattering angles. This is why it is important to have various propagation regimes and various arrivals types originating at the same structure to broaden the range of scattering angles (transmission, reflections and wide-angle reflections). Ideally, with omnidirectional illumination of a structure, the resolution can theoretically be half wavelength. Practically, optimal illumination conditions are never reached. It also exhibits the need for broad frequency data with good quality low frequency to easily start inversion and avoid cycle-skipping and high frequency to obtain high-resolution reconstruction.

Interestingly, when the FWI problem was formulated, the community missed two main ingredients to ensure successful results. The first one is that it was introduced in the early age of computational science. That led to serious difficulty solving the wave equation even in limited simple 2D cases as computational power was fairly limited. The second one was the Figure 1.1: Shematic illustration of the accuracy and resolution of tomography, full-waveform inversion and migration methods. The solid black line corresponds to the resolution gap between tomography and migration technics (as proposed by [START_REF] Claerbout | Imaging the Earth's interior[END_REF]. This gap is supposed to be filled by full-waveform inversion (figure from [START_REF] Yao | A review on reflection-waveform inversion[END_REF]. availability of appropriate seismic data, and a better understanding of the type of data needed to sufficiently constrain the problem (seismic data were mainly short-offset seismic reflection surveys). That explained why FWI was difficult to apply [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: numerical results[END_REF][START_REF] Cary | Automatic 1-D waveform inversion of marine seismic refraction data[END_REF][START_REF] Crase | Robust elastic nonlinear waveform inversion: application to real data[END_REF][START_REF] Jin | Two-dimensional asymptotic iterative elastic inversion[END_REF][START_REF] Lambaré | Iterative asymptotic inversion in the acoustic approximation[END_REF]. These limitations (inappropriate data and computational power) slowed the development of the method, but, in the '90s, Pratt's work brought the essential bricks that unleashed FWI understanding: the formulation of a 2D frequency-domain solver that considerably lighten the computational cost [START_REF] Pratt | Frequency-domain elastic modeling by finite differences: a tool for crosshole seismic imaging[END_REF][START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. Part I: acoustic wave-equation method[END_REF][START_REF] Pratt | Combining wave-equation imaging with traveltime tomography to form high-resolution images from crosshole data[END_REF][START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by wavefield inversion[END_REF][START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain, part I: theory and verification in a physical scale model[END_REF]. By seeing FWI in terms of wavenumber covering and acquisition offset, he showed the potential of FWI for crosshole acquisition with transmitted waves. This work was essential to lay the foundation for the first field data applications in the years 2000 when surface acquisition presenting sufficiently large offset coverage started to show up and finally unlock FWI applications with meaningful updates at more than a few hundred meters.

Formulation of FWI

Historically, the method was born in the context of exploration (related to the oil & gas industry). Since then, FWI has been continuously developed and improved, from 2D acoustic mono-parameter [START_REF] Ravaut | Seismic imaging of a thrust belt by traveltime and frequency-domain waveform inversions[END_REF][START_REF] Ravaut | Multi-scale imaging of complex structures from multi-fold wide-aperture seismic data by frequency-domain full-wavefield inversions: application to a thrust belt[END_REF][START_REF] Sourbier | FWT2D: A massively parallel program for frequency-domain full-waveform tomography of wide-aperture seismic data-Part 2: Numerical examples and scalability analysis[END_REF], 3D visco-acoustic mono or multi-parameters [START_REF] Plessix | Three-dimensional frequency-domain full-waveform inversion with an iterative solver[END_REF][START_REF] Malinowski | Full-waveform modeling and inversion of wide-aperture data from the Polish basin[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 full-waveform inversion[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 approximation[END_REF][START_REF] Górszczyk | Toward a robust workflow for deep crustal imaging by FWI of OBS data: The eastern nankai trough revisited[END_REF][START_REF] Górszczyk | Graph-space optimal transport concept for time-domain full-waveform inversion of ocean-bottom seismometer data: Nankai trough velocity structure reconstructed from a 1d model[END_REF][START_REF] Kamath | Multiparameter full-waveform inversion of 3D OBC data from the valhall field[END_REF][START_REF] Pladys | Graph space optimal transport based FWI: 3D OBC valhall case study[END_REF] to 3D visco-elastic inversion [START_REF] Zhu | Seismic structure of the European upper mantle based on adjoint tomography[END_REF][START_REF] Raknes | Three-dimensional elastic full waveform inversion using seismic data from the sleipner area[END_REF]Trinh et al., 2019b,c,a) Meanwhile, rapidly after the introduction of FWI at exploration scale, global to regional scale FWI from the academic community emerged, with groups from Caltech/Princeton [START_REF] Tromp | Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels[END_REF][START_REF] Liu | Efficient modeling of q for 3d numerical simulation of wave propagation[END_REF][START_REF] Tape | Finite-frequency tomography using adjoint methods-methodology and examples using membrane surface waves[END_REF][START_REF] Bedle | S velocity variations beneath north america[END_REF][START_REF] Tape | Seismic tomography of the southern California crust based on spectral-element and adjoint methods[END_REF][START_REF] Bozdag | Global adjoint tomography: first-generation model[END_REF][START_REF] Tromp | Seismic wavefield imaging of earth's interior across scales[END_REF], Berkeley [START_REF] Panning | Importance of crustal corrections in the development of a new global model of radial anisotropy[END_REF][START_REF] French | Broad plumes rooted at the base of the earth's mantle beneath major hotspots[END_REF] and Munich/Utrecht/ETH [START_REF] Fichtner | Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods[END_REF][START_REF] Fichtner | Hessian kernels of seismic data functionals based upon adjoint techniques[END_REF][START_REF] Fichtner | Multiscale full waveform inversion[END_REF].

More recently, smaller-scale approaches have emerged and have made drastic improvements in the past years [START_REF] Köhn | 2d full waveform inversion applied to a strongly-dispersive love wave field dataset[END_REF][START_REF] Wittkamp | Individual and joint 2-d elastic full-waveform inversion of rayleigh and love waves[END_REF][START_REF] Smith | Tunnel detection at yuma proving ground, arizona, usa-part 2: 3d full-waveform inversion experimentstunnel detection at ypg-part 2: 3d fwi[END_REF][START_REF] Irnaka | 3D elastic full waveform inversion for subsurface characterization: Study of a shallow seismic multicomponent field data[END_REF]. This remains a challenging topic as many difficulties are faced at this scale, such as propagation of surface waves, poorly consolidated area, strong heterogeneities of mechanical properties in the first tens to hundreds of meters.

Finally, the interest of FWI is now going beyond the scope of seismic waves, with successful application on ground penetrating radar (GPR) [START_REF] Klotzsche | Full-waveform inversion of cross-hole ground-penetrating radar data to characterize a gravel aquifer close to the Thur River, Switzerland[END_REF][START_REF] Minet | Soil surface water content estimation by full-waveform GPR signal inversion in the presence of thin layers[END_REF][START_REF] Klotzsche | Crosshole GPR fullwaveform inversion of waveguides acting as preferential flow paths within aquifer systems[END_REF][START_REF] Klotzsche | 3-D characterization of high-permeability zones in a gravel aquifer using 2-D crosshole GPR full-waveform inversion and waveguide detection[END_REF][START_REF] Lavoué | Two-dimensional permittivity and conductivity imaging by full waveform inversion of multioffset GPR data: a frequency-domain quasi-Newton approach[END_REF][START_REF] Klotzsche | Review of crosshole groundpenetrating radar full-waveform inversion of experimental data: Recent developments, challenges, and pitfalls[END_REF][START_REF] Giannakis | A machine learning-based fast-forward solver for ground penetrating radar with application to full-waveform inversion[END_REF], or even in the medical imaging field that starts to use FWI [START_REF] Guasch | Full-waveform inversion imaging of the human brain[END_REF][START_REF] Marty | Acoustoelastic full-waveform inversion for transcranial ultrasound computed tomography[END_REF]. This proves that now FWI is a mature method used in a wide variety of scales, targeting multiple imaging problems.

In the following, I introduce the physical and mathematical concepts related to FWI before developing on the actual limitations of the conventional FWI formulation.

Formulation of FWI

Full-waveform inversion is formulated as a minimization of a function F measuring the distance between observed (recorded) data d obs and synthetics data d cal [m] calculated in a discrete physical model of the subsurface parameters m. These data, d obs and d cal [m], are extracted from the wavefield at receivers' locations. FWI aims to retrieve the model m from the observed data, categorizing FWI into the inverse problem family. By taking into account the complete signal information (all the phases and all the amplitudes), this method is supposed to provide a high-resolution reconstruction of the subsurface model. Several strategies can be used to solve eq. 1.2, but first, it is required to detail the forward problem A[m] which is compulsory to obtain d cal [m]

The forward problem

One way of seeing the seismic waves is as optic rays based on a high-frequency approximation [START_REF] Červený | Seismic Ray Theory[END_REF][START_REF] Chapman | Fundamentals of seismic waves propagation[END_REF][START_REF] Virieux | Theory and observations -body waves: ray methods and finite frequency effects[END_REF]. This would give access to the traveltime, which leads to raypath connecting a source-receiver couple, that can be solved by different methods [START_REF] Zelt | Seismic traveltime inversion for 2-D crustal velocity structure[END_REF][START_REF] Bishop | Tomographic determination of velocity and depth in laterally varying media[END_REF][START_REF] Le Meur | Seismic tomography of the gulf of Corinth: a comparison of methods[END_REF][START_REF] Popovici | Three dimensional traveltimes using the fast marching method[END_REF][START_REF] Lelièvre | Computing first-arrival seismic traveltimes on unstructured 3-D tetrahedral grids using the fast marching method[END_REF][START_REF] Vidale | Finite-difference calculation of traveltimes in three dimensions[END_REF][START_REF] Bretaudeau | First-arrival delayed tomography using 1st and 2nd order adjoint-state method[END_REF][START_REF] Bouteiller | A discontinuous Galerkin fast-sweeping Eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media[END_REF].

The second way, which is of interest here, is modeling the full seismic wavefield, referred to as forward modeling. To do so, we consider the wave equation to approximate as close as possible the real physics of seismic wave propagation in a given medium.

I will detail the elastodynamic equations in the following parts and derive the acoustic and the anisotropic visco-acoustic modeling. These two modeling choices are motivated by their usage in the FWI application performed in this thesis. The first one -acoustic -is the most "basic" modeling we can think of. It only accounts for the pressure waves. The second one introduces anisotropic effects and attenuation, while still only modeling for the pressure waves. This approximation of the physics is important to keep in mind for field data application as the Earth is elastic, and contains more complex anisotropy than the VTI approximation used here. This, by definition, will never allow obtaining a perfect fit of the modeling to the field data as we will miss a part of the physics. This could also lead to modification of physical properties such as P-waves velocity to take into account other physical effects that is not taken into account in the modeling, leading to the creation of non-physical updates. This issue is commonly referred to as cross-talk [START_REF] Operto | A guided tour of multiparameter full waveform inversion for multicomponent data: from theory to practice. The Leading Edge[END_REF], but that will not be the focus of this work.

The elastodynamic equations

In the following, we follow Einstein convention (summation over repeated indices) for these equations. Generally, to describe the elastodynamic system, we rely on the following Newton's law

ρ∂ t v i = ∂ j σ ij , (1.5)
and on Cauchy's generalization to the Hooke's law in tensorial form

σ ij = c ijkl ε kl , (1.6)
where ρ is density, σ is the stress tensor, c ijkl is the stiffness tensor, ε is the strain tensor with i, j, k, l ∈ [x, y, z]. Strains are related to displacement vector u = (u x , u y , u z ) through the expression

ε ij = 1 2 (∂ i u j + ∂ j u i ) , (1.7) 1.2 Formulation of FWI whose first time derivative is εij = 1 2 (∂ j ui + ∂ i uj ) = 1 2 (∂ j v i + ∂ i v j ) . (1.8)
Simplification of the stiffness tensor can be obtained due to the symmetry as well as the unique energy definition (Carcione, 2001, p. 55):

c ijkl = c jikl = c ijlk = c klij , (1.9)
which reduces the components to only 21.

We introduce the particle velocity vector as v = (v x , v y , v z ) T and the stress vector as σ = (σ xx , σ yy , σ zz , σ yz , σ xz , σ xy ) T . We can now define the first time derivative for the Hooke's law as the following matrix expression:

∂t         σxx σyy σzz σyz σxz σxy         σ =         c11 c12 c13 c14 c15 c16 • c22 c23 c24 c25 c26 • • c33 c34 c35 c36 • • • c44 c45 c46 • SY M • • c55 c56 • • • • • c66         C         εxx εyy εzz 2 εyz 2 εxz 2 εxy         ε =         c11 c12 c13 c14 c15 c16 • c22 c23 c24 c25 c26 • • c33 c34 c35 c36 • • • c44 c45 c46 • SY M • • c55 c56 • • • • • c66         C         ∂x 0 0 0 ∂y 0 0 0 ∂z 0 ∂z ∂y ∂z 0 ∂x ∂y ∂x 0         D   vx vy vz   v , (1.10)
where we have introduced a differential operator D and the symmetric stiffness matrix C using Voigt indexing, (11) → 1, (22) → 2, (33) → 3, (23) = (32) → 4, (13) = (31) → 5, (12) = (21) → 6.

Anisotropic VTI acoustic wave approximation

We can now simplify the stiffness matrix C depending on which type of anisotropy we want. In the case of orthorhombic anisotropy, the stiffness matrix C reduces to

C =         c 11 c 12 c 13 0 0 0 • c 22 c 23 0 0 0 • • c 33 0 0 0 • • • c 44 0 0 • SY M • • c 55 0 • • • • • c 66        
.

(1.11)

The stiffness matrix C can be further simplified with a special case of orthorhombic anisotropy known as VTI anisotropy. The simplification is obtained by assuming a symmetric (x, y)-plane and setting

σ xx = σ yy , σ xz = σ yz , ε xx = ε yy , ε xz = ε yz , yielding c 11 = c 22 , c 13 = c 23 , c 44 = c 55 .
The resulting stiffness matrix becomes [START_REF] Duveneck | Acoustic VTI wave equations and their application for anisotropic reverse-time migration[END_REF]) Finally, to obtain the acoustic wave approximation in the presence of VTI anisotropy, we force shear wave velocity on the symmetry axis to be zero in the elastic case in the sense that c 55 = c 66 = 0, c 12 = c 11 . Therefore only first 3 rows of C remains:

C =         c 11 c 12 c 13 0 0 0 • c 11 c 13 0 0 0 • • c 33 0 0 0 • • • c 55 0 0 • SY M • • c 55 0 • • • • • c 66         , ( 1
C =         c 11 c 11 c 13 0 0 0 • c 11 c 13 0 0 0 • • c 33 0 0 0 • • • 0 0 0 • SY M • • 0 0 • • • • • 0         . (1.13)
To be rigorous, we can note that anisotropy does not exist in a pure acoustic regime. While forcing shear velocity to be zero, it does not prevent the propagation of fictitious shear waves [START_REF] Alkhalifah | Acoustic approximations for processing in transversely isotropic media[END_REF].

The resulting VTI acoustic system is written as

                     ρ∂ t v x = ∂ x σ xx ρ∂ t v y = ∂ y σ yy ρ∂ t v z = ∂ z σ zz ∂ t σ xx = c 11 (∂ x v x + ∂ y v y ) + c 13 ∂ z v z ∂ t σ yy = c 11 (∂ x v x + ∂ y v y ) + c 13 ∂ z v z ∂ t σ zz = c 13 (∂ x v x + ∂ y v y ) + c 33 ∂ z v z .
(1.14)

Let us denote σ xx = σ yy := g, and σ zz := q. The pressure is the averaged trace of diagonal stress matrix:

p = - 1 3 i σ ii = - 1 3 (2g + q) .
(1.15)

Finally, we can introduce Thomsen's anisotropy parameters δ and as [START_REF] Duveneck | Stable P-wave modeling for reverse-time migration in tilted TI media[END_REF][START_REF] Zhang | A stable TTI reverse time migration and its implementation[END_REF]:

     c 11 = ρV 2 p (1 + 2 ) = ρV 2 h c 13 = c 33 √ 1 + 2δ = ρV 2 p √ 1 + 2δ c 33 = ρV 2 p = κ , (1.16)
where κ is the bulk modulus while V h = V p √ 1 + 2 denotes the horizontal velocity. Allowing for the external sources, the VTI acoustic system 1.14 becomes

               ρ∂ t v x = ∂ x g + f vx ρ∂ t v y = ∂ y g + f vy ρ∂ t v z = ∂ z q + f vz ∂ t g = κ[(1 + 2 )(∂ x v x + ∂ y v y ) + √ 1 + 2δ∂ z v z ] + f σ xx/yy ∂ t q = κ[ √ 1 + 2δ(∂ x v x + ∂ y v y ) + ∂ z v z ] + f σzz .
(1.17)

Formulation of FWI

We can further simplify the system by setting and δ to zero in eq 1.17. This yields the isotropic acoustic wave equation: (1.18) where we have considered the pressure as defined in eq 1.15.

           ρ∂ t v x = -∂ x p + f vx ρ∂ t v y = -∂ y p + f vy ρ∂ t v z = -∂ z p + f vz ∂ t p = -κ(∂ x v x + ∂ y v y + ∂ z v z ) + f p ,

Anisotropic VTI viscoacoustic wave approximation

We considered before the propagation of waves in purely elastic solid. However, as waves propagate through the Earth, the wavefront loses energy (through thermal loss, for instance). This is called attenuation, and it is an important feature that should be taken into account. This effect (reduction of amplitude over distance) has been observed and confirmed from a wide range of experimental tests and field observations. Depending on the target, accounting for attenuation can be crucial to obtain good results of FWI when applied to field data; see [START_REF] Kamath | Multiparameter full-waveform inversion of 3D OBC data from the valhall field[END_REF] for instance.

In a first attempt to account for attenuation, constant-Q models were developed to approximate seismic attenuation in seismic exploration and seismology. It is mathematically simple to implement Q in the frequency domain. However, it is more complex in the time domain.

In order to simulate attenuation for time-domain modeling, the generalized Maxwell body (GMB), which is equivalent to the superposition of several standard linear solid (SLS) mechanisms, is used to approximate the constant Q over a specified frequency range.

The SLS mechanisms correspond to attenuation processes (e.g. interstitial atom relaxation, grain boundary relaxation, thermoelasticity, the diffusional motion of dislocations and point defects) that affect the relaxation time following the wave propagation. The general SLS rheology explains these processes very well. Some of them can be modeled with one mechanism, and others using several relaxation mechanisms.

A practical problem in superposing SLS mechanisms is determining the appropriate number of mechanisms, providing the desired constant-Q behavior, and saving computational costs for solving wave equations, particularly in 3D modeling. In general, three SLS mechanisms are considered to be accurate enough for 3D simulations in geophysical prospecting and global seismology study. Less there are SLS mechanisms, and fewer memory-variable equations have to be solved, resulting in time-saving. Now, let us introduce, for a 1D case, the rheology of generalized Maxwell body (GMB) [START_REF] Emmerich | Incorporation of attenuation into time-domain computation of seismic wavefield[END_REF][START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF]:

∂ t σ = M u ( ε -L =1 Y ξ ) ∂ t ξ + ω ξ = ω ε, = 1, 2, . . . , L , (1.19)
where M u is the unrelaxed modulus, ξ is the memory variable independent of the model properties to delineate the historical energy dissipation associated with -th circular frequency ω and the dimensionless anelastic coeffcients Y . From eq 1.19, we may easily find that moving

FULL-WAVEFORM INVERSION

from elastics to viscoelastics is simply subtracting the weighted memory variables ξ with the strain term ˙ in lossless elastic wave equation. By analogy with the VTI acoustic wave equation, the VTI viscoacoustic system in 3D reads

                                               ∂ t            σ xx σ yy σ zz σ xy σ xz σ yz            σ =            c 11 c 11 c 13 0 0 0 • c 11 c 13 0 0 0 • • c 33 0 0 0 • • • 0 0 0 • SY M • • 0 0 • • • • • 0            C            εxx -L =1 Y ξ xx εyy -L =1 Y ξ yy εzz -L =1 Y ξ zz 2 εxy -L =1 Y ξ xy 2 εxz -L =1 Y ξ xz 2 εyz -L =1 Y ξ yz            ε-L =1 Y ξ ∂ t    ξ xx ξ yy ξ zz    ξ +ω    ξ xx ξ yy ξ zz    ξ = ω    εxx εyy εzz    ε = 1, 2, . . . , L , (1.20)
where Y := Y κ the shear modulus µ (and the corresponding anelastic coefficients) are forced to be zero in the elastic case in order to perform the acoustic simplification. Allowing for the Newton's law in eq 1.5, the VTI viscoacoustic wave equation is now defined as

                                   ρ∂ t v x = ∂ x σ xx ρ∂ t v y = ∂ y σ yy ρ∂ t v z = ∂ z σ zz ∂ t σ xx = c 11 (∂ x v x + ∂ y v y ) + c 13 ∂ z v z -L =1 Y [c 11 (ξ xx + ξ yy ) + c 13 ξ zz ] ∂ t σ yy = c 11 (∂ x v x + ∂ y v y ) + c 13 ∂ z v z -L =1 Y [c 11 (ξ xx + ξ yy ) + c 13 ξ zz ] ∂ t σ zz = c 13 (∂ x v x + ∂ y v y ) + c 33 ∂ z v z -L =1 Y [c 13 (ξ xx + ξ yy ) + c 33 ξ zz ] ∂ t ξ xx + ω ξ xx = ω ∂ x v x , = 1, 2, . . . , L ∂ t ξ yy + ω ξ yy = ω ∂ y v y , = 1, 2, . . . , L ∂ t ξ zz + ω ξ zz = ω ∂ z v z , = 1, 2, . . . , L .
(1.21) Similar to VTI acoustic case, based upon the notation in eq 1.15, the VTI viscoacoustic wave equation including the external sources can then be expressed as

                         ρ∂ t v x = ∂ x g + f vx ρ∂ t v y = ∂ y g + f vy ρ∂ t v z = ∂ z q + f vz ∂ t g = c 11 (∂ x v x + ∂ y v y ) + c 13 ∂ z v z -L =1 Y [c 11 (ξ xx + ξ yy ) + c 13 ξ zz ] + f g ∂ t q = c 13 (∂ x v x + ∂ y v y ) + c 33 ∂ z v z -L =1 Y [c 13 (ξ xx + ξ yy ) + c 33 ξ zz ] + f q ∂ t (ξ xx + ξ yy ) + ω (ξ xx + ξ yy ) = ω (∂ x v x + ∂ y v y ), = 1, 2, . . . , L ∂ t ξ zz + ω ξ zz = ω ∂ z v z , = 1, 2, . . . , L .
(1.22) Note that the memory variables in the symmetric plane (x, y), namely ξ xx and ξ yy , are summed together as one quantity to reduce the volume of memory allocation in practical implementa-1.2 Formulation of FWI tion. With the notation ξ xx + ξ yy := ξ g , ξ zz := ξ q , system 1.22 becomes

                         ρ∂ t v x = ∂ x g + f vx ρ∂ t v y = ∂ y g + f vy ρ∂ t v z = ∂ z q + f vz ∂ t g = c 11 (∂ x v x + ∂ y v y ) + c 13 ∂ z v z -L =1 Y [c 11 ξ g + c 13 ξ q ] + f g ∂ t q = c 13 (∂ x v x + ∂ y v y ) + c 33 ∂ z v z -L =1 Y [c 13 ξ g + c 33 ξ q ] + f q ∂ t ξ g = -ω ξ g + ω (∂ x v x + ∂ y v y ), = 1, 2, . . . , L ∂ t ξ q = -ω ξ q + ω ∂ z v z , = 1, 2, . . . , L .
(1.23)

Discretization and absorbing boundary coniditions

The wave equation is a linear partial differential equation and can be described and solved with various schemes [START_REF] Kelly | Synthetic seismograms -a finitedifference approach[END_REF][START_REF] Marfurt | Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[END_REF][START_REF] Virieux | SH wave propagation in heterogeneous media: Velocity-stress finite difference method[END_REF][START_REF] Dablain | The application of high order differencing for the scalar wave equation[END_REF][START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF][START_REF] Brossier | Parsimonious finite-volume frequency-domain method for 2-D P-SV-wave modelling[END_REF]. The choice of the schemes and the choice of discretization will influence the numerical resolution in which the wavefield is resolved.

One of the most popular methods used to solve the wave equation relies on finite-difference schemes. The main advantages of this method are the good scalability, effectiveness, and simplicity of implementation.

As this manuscript focuses on more robust FWI, and as the applications will be performed on marine field data, the excellent scalability of finite-difference approaches will be needed to tackle large 3D field data. Therefore, this work relies on using a code based on finitedifference, based on staggered-grid methods [START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite difference method[END_REF][START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF]. The discretization is performed with fourth-order in space and second-order in time. To localize sources and receivers that do not directly fall on the finite-difference grid, we rely on Hicks interpolation [START_REF] Hicks | Arbitrary source and receiver positioning in finite-difference schemes using Kaiser windowed sinc functions[END_REF]. This allows recasting the position on the grid using windowed sinc functions.

Finally, artificial reflections generated at the boundary are mitigated by modifying the equations at a continuous level. In this work, we rely on absorbing boundary conditions, more precisely, perfectly matched layers [START_REF] Bérenger | A perfectly matched layer for absorption of electromagnetic waves[END_REF] for 2D acoustic FWI. For 3D visco-acoustic, we rely on "sponge" absorbing boundaries as PML would be challenging to implement in this case [START_REF] Cerjan | A nonreflecting boundary condition for discrete acoustic and elastic wave equations[END_REF].

Solution to the FWI problem

Solving the minimization problem of FWI can be done in two ways. As developed later on in this manuscript, the misfit function conventionally used for FWI -the least-squares norm -is non-convex and presents local minima. Therefore, one of the first ideas that would come to the mind would be relying on global optimization to solve the FWI problem. This would ensure convergence toward the global minimum. The second idea that can be used is to rely on local optimization. Details of the two different approaches will be introduced in the following, with their advantages and disadvantages.

FULL-WAVEFORM INVERSION

Global solution to the FWI problem

This solution was explored in the early years of FWI formulation [START_REF] Sen | Non-linear one-dimensional seismic waveform inversion using simulated annealing[END_REF][START_REF] Stoffa | Nonlinear multiparameter optimization using genetic algorithms: Inversion of plane-wave seismograms[END_REF] but also more recently [START_REF] Tran | Two-dimensional inversion of full waveforms using simulated annealing[END_REF][START_REF] Datta | Estimating starting models for full waveform inversion using a global optimization method[END_REF][START_REF] Aleardi | 1d elastic full-waveform inversion and uncertainty estimation by means of a hybrid genetic algorithm-gibbs sampler approach[END_REF][START_REF] Mazzotti | Two-grid genetic algorithm full-waveform inversion[END_REF][START_REF] Sajeva | Genetic algorithm full-waveform inversion: uncertainty estimation and validation of the results[END_REF].

Global optimization method aims at finding the global minima or maxima of a continous function (here a misfit function). The function can possibly be non-convex.

f [m] : Ω ⊂ R n → R , (1.24)
with Ω is the (not necessarily convex) solution space, discretized over n parameters.

By handling potentially non-convex misfit function, global optimization methods can naturally deal with local minima in Ω. Therefore these methods naturally allow tackling the ill-posed inverse problem that FWI represents. Global solutions aim to explore the solution space Ω to evaluate all possible minima (or maxima) and find the global one. This could be illustrated in the example shown in Figure 1.3, global optimization method would allow exploring the solution space Ω entirely and find the possible solutions (here, the red dots) while avoiding the local minimum represented by the blue dot.

The main disadvantage of global search approaches is the difficulty faced when tackling problems with a large number of degrees of freedom. Recent applications using Hamiltonian Markov-Chain Monte-Carlo (H-MCMC) tackle problems with up to thousands of degrees of freedom (Fichtner et al., 2018a,b;[START_REF] Gebraad | Bayesian elastic full-waveform inversion using hamiltonian monte carlo[END_REF][START_REF] Gebraad | Bayesian elastic full-waveform inversion using hamiltonian monte carlo[END_REF]. The issue is that even small 2D synthetic studies correspond to tens to hundreds of thousands of degrees of freedom. This is known as the "curse of dimensionality". Indeed, the computational cost of this method is directly linked with the solution space Ω, which is directly related to the number of degrees of freedom representing the model. Realistically, and considering the trend of moving toward high-frequency 3D FWI application, global optimization seems too computationally intensive to be used. Despite this issue, the method remains of interest for FWI as it could be used for uncertainty estimation [START_REF] Martin | A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion[END_REF][START_REF] Bardsley | Randomize-then-optimize: A method for sampling from posterior distributions in nonlinear inverse problems[END_REF][START_REF] Biswas | 2d full-waveform inversion and uncertainty estimation using the reversible jump hamiltonian monte carlo[END_REF][START_REF] Sajeva | Genetic algorithm full-waveform inversion: uncertainty estimation and validation of the results[END_REF][START_REF] Thurin | Uncertainties estimation in Full Waveform Inversion using Ensemble methods[END_REF] and low-frequency inversion.

Local solution to the FWI problem

For problem space larger than 10 1 to 10 3 , the solution used to solve the FWI inverse problem relies on local optimization techniques. This method aims at finding the closest minimum for a given starting point, using local information provided by the misfit function (such as the gradient and curvature). This highlights a crucial parameter: the choice of initial model m 0 used to start the inversion. This initial model m 0 defines a subspace A ⊂ Ω, in which the misfit function is minimized. This drastically reduces the number of potential solutions but also implies that only the closest minimum is reached. Therefore a good enough initial model is needed to allow convergence of the scheme toward the global minimum if the misfit function used is not globally convex.

The formalism consists in iteratively minimizing the discrepency between observed and synthetic data by building the following sequence of models where m 0 the initial model, α is the step length, ∆m is the model update (which corresponds to the descent direction) and k the index of the current iteration.

m k+1 = m k + α k ∆m k , (1.25)
At each iteration k, both the step length α and the direction of descent ∆m are updated to provide a better data fit

f [m k+1 ] < f [m k ] ,
(1.26)

The steplength α should satisfy the Wolfe criterion [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF] which are

1. Sufficient decrease f [m + α∆m] ≤ f [m] + c 1 α∇f [m] T ∆m , (1.27)
where ∇f [m] the gradient. Practicaly c 1 is set to 10 -4 in practice.

Curvature condition

∇f [m + α∆m] T ∆m ≥ c 2 ∇f [m] T ∆m , (1.28)
with c 2 = 0.9 in practice.

Satisfying the Wolfe conditions ensures the convergence towards the nearest local minimum [START_REF] Nocedal | Numerical Optimization[END_REF].

To build the model update ∆m k , we usually rely on an approximate solution of the Newton equation

H[m]∆m = -∇f [m] , (1.29)
where H[m] is the Hessian operator, the matrix of second-order derivatives of the function f [m]. As the computational cost of inverting for the Hessian operator is too high, we rely on an approximation P : This steepest-descent algorithm is known to be slow and exhibits convergence difficulties even on elementary problems.

P [m] H[m] -1 . ( 1
To avoid this issue, we rely on quasi-Newton methods that build an approximation of the inverse Hessian for P . The most famous of this method is called l-BFGS and build the approximation through the stored gradient and model values [START_REF] Nocedal | Updating Quasi-Newton Matrices With Limited Storage[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF]. This method will be the one used in this thesis.

Gradient computation

From the description of the quasi-Newton scheme, we see that in practice the critical quantity to compute is the gradient of f [m]. In this section, we explain how we can compute it. We define the number of sources (shots) as N s (with s a subscript referring to shot number) and the number of receivers as N r . We assume that the number of receivers is the same for each shot.

First, we consider the case of the least-squares norm: For convenience, we now consider a simple case with N s = 1. For a perturbation dm we have

min m f [m] = Ns s=1 F (d cal,s [m], d obs,s ) , (1.32
f [m + dm] = 1 2 d cal (m + dm) -d obs 2 1 2 d cal (m) -d obs + J(m)dm + o( dm 2 ) 2 , (1.37)
where

J(m) = ∂d cal ∂m =     ∂(d cal ) 1 ∂m 1 ∂(d cal ) 1 ∂m 2 . . . ∂(d cal ) 2 ∂m 1 ∂(d cal ) 2 ∂m 2 . . . . . .     (1.38) 1.2 Formulation of FWI is the Jacobian matrix. It can be used to re-write f [m + dm] as f [m + dm] = 1 2 d cal (m) -d obs 2 + d cal -d obs , J(m)dm + o( dm 2 ) 1 2 d cal (m) -d obs 2 + J(m) T (d cal -d obs ) , dm + o( dm 2 ) . (1.39) Therefore f (m + dm) -f (m) = J(m) T (d cal -d obs ), dm + o( dm 2 ) , (1.40) and ∇f (m) = J(m) T (d cal -d obs ) . (1.41)
This corresponds to the "direct" method to obtain the gradient. Unfortunately, this method relies on a Jacobian matrix which is computationally too intensive to compute or too large to be stored (I/O storage is too slow).

Therefore, we rely on the adjoint-state approach [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF]. This method relies on using the Lagrangian functional. The main goal of the adjoint state method is to compute the gradient of a functional, f [m], where f depends on u(m), without explicitly computing the Jacobian matrix J.

We introduce the Lagrangian operator in the case of N s = 1 with an arbitrary distance measurement function F :

L(m, u, d cal , λ, µ) = F (d cal , d obs ) + A(m)u -b, λ + d cal -Ru, µ .
(1.42)

Let ū(m) be the solution of the forward problem for a given m (eq 1.36). Similarly, if d cal is extracted from ū(m), we denote it by dcal (m) (eq 1.35). We thus have

L(m, ū(m), dcal (m), λ, µ) = F ( dcal (m), d obs ) = f [m] . (1.43)
Therefore, by deriving the previous equation with respect to the model parameter m we obtain an expression of the gradient of f [m] as:

∂L(m, ū(m), dcal (m), λ, µ) ∂m = ∇f [m] . (1.44)
Using the chain rule on the Lagrangian, we obtain

∂A(m) ∂m ū(m), λ + ∂L(m, ū(m), dcal (m), λ, µ) ∂u ∂ ū(m) ∂m + ∂L(m, ū(m), dcal (m), λ, µ) ∂d cal ∂ dcal (m) ∂m = ∇f [m] .
(1.45) We can identify in ∂ ū(m) ∂m and ∂ dcal (m) ∂m operators similar to the Jacobian one. As we do not want to compute such operators (it was to avoid them that we introduce the Lagrangian), we define λ and µ such as

∂L(m, ū(m), dcal (m), λ, µ) ∂u = 0 , (1.46) ∂L(m, ū(m), dcal (m), λ, µ) ∂d cal = 0 . (1.47) FULL-WAVEFORM INVERSION
To obtain λ and µ that satisfies the two latter equations, we express the derivatives of the Lagrangian with respect to u and d cal using first order Taylor devlopments. We have (1.48) and similarly ∂L(m, ū(m), dcal (m), λ, µ) ∂d cal = ∂F ∂d cal -µ .

∂L(m, ū(m), dcal (m), λ, µ) ∂u = A(m) T λ -R T µ ,
(1.49)

This provides an equivalent to the system 1.46 as:

A(m) T λ -R T µ = 0 (1.50) ∂F ∂d cal -µ = 0 , (1.51)
which can be rewritten as

A(m) T λ = R T µ (1.52) µ = ∂F ∂d cal , (1.53)
with the index T representing transpose. From the two latter equations, we can now extract conditions to cancel the terms in the gradient expression 1.45 which are:

1. µ is the adjoint source, which is just the derivative with respect to d cal of the distance measurement function F 2. λ is the solution of an adjoint equation assocaited with the operator A(m) T with a source term equal to R T µ

This finally allows to have the expresion of the gradient relying only on the correlation part as

∂A ∂m u, λ = ∇f [m] ,
(1.54) with λ the adjoint field

A(m) T λ = -R T ∂F ∂d cal .
(1.55)

A very interesting property of this method is that a change in the misfit function only translates into a modification of the source term used for the adjoint field back-propagated from the receivers.

The gradient computation requires solving two wave propagation problems: one incident field and one adjoint field. This translates into a reduction of the computational cost compared to forming or storing the Jacobian matrix. Since the introduction of the FWI formalism in the '80s, the least-squares norm -referred to as the L 2 norm -has been traditionally used as the "default" norm to measure the distances between the observed data and the calculated data (Tarantola, 1984a[START_REF] Tarantola | Theoretical background for the inversion of seismic waveforms including elasticity and attenuation[END_REF][START_REF] Nolet | Seismic tomography with applications in global seismology and exploration geophysics[END_REF].

Let us remind the formulation of the L 2 norm already introduced in eq 1.33 for an isolated traces such as N r = N s = 1

F L 2 (d cal , d obs ) = 1 2 T 0 |d cal (t) -d obs (t)| 2 dt . (1.56)
As demonstrated earlier, FWI formalism is based on the adjoint theory, requiring an adjoint source associated with the misfit function. For the L 2 norm, the adjoint source S is defined as

S L 2 (t) = ∂F L 2 ∂d cal (t) (1.57) = ∂ ∂d cal (t) 0.5 × (d cal (t) -d obs (t)) 2 .
(1.58) By using the chain rule, we obtain the adjoint source

S L 2 (t) = 0.5 × 2 × (d cal (t) -d obs (t)) × ∂(d cal (t) -d obs (t)) ∂d cal (t) (1.59) = d cal (t) -d obs (t) .
(1.60)

In the case of the L 2 norm, the obtained adjoint source is simply the difference between d cal and d obs . This is obviously straightforward to implement, it is not affected by any stability condition and associated computational cost is close to zero.

The L 2 misfit function presents several advantages, such as noise robustness and "usage" of the complete signal information. This explains why it is still widely used, and it is still the "default" choice for field data application. The main limitation of the L 2 norm is its lack of convexity for time-delays, leading to cycle-skipping being one of its most substantial weaknesses.

The cycle-skipping limitation

The L 2 distance measure is inherently affected by the phase ambiguity, known as cycle-skipping, which appears when the synthetic data are shifted by more than half a period compared to the observed ones (Figure 1.4). Unfortunately, finding perfectly convex misfit functions is challenging, resulting in FWI being strongly non-linear when based on a non-convex misfit function (such as the least-square norm conventionally used). This leads to the common issue of converging toward a non geologically informative local minimum, an issue that has been FULL-WAVEFORM INVERSION faced since the beginning of FWI applications [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: numerical results[END_REF], and which translates into an erroneous reconstruction of the velocity model [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF]. This nonconvexity with respect to time-shifts link with Claerbout's "gap" between seismic imaging methods. Indeed, if the misfit function could be made perfectly convex, it would allow starting from a low-wavenumber initial model and iteratively build high-wavenumber updates, thereby bridging the gap.

The non-convexity of the L 2 misfit function with respect to time-shifts can be illustrated easily in a straightforward and schematic example based on time-shifted Ricker wavelets. Let us introduce a reference signal composed of one Ricker wavelet in the center, seen as the observed data. Then, a second signal is created from a shifted in time version of the reference signal (with a time-shift going from -1.5 to 1.5 s). The distance between the reference signal and the calculated signal is then obtained using the L 2 norm. Results are presented in Figure 1.5.

In this example, it can be observed that from time-shift superiors to ±0.12 s, the L 2 norm is not convex and presents local minimum and a flat part for time-shift superior to ±0.4 s. This non-convex behavior of the L 2 norm for time shift superiors to half the signal period corresponds to the so-called cycle-skipping issue. As FWI relies on local optimization schemes to be solved, converging toward a local minimum implies an erroneous reconstructed V P model, with no possibilities to improve it. Practically, the origin of time shifts between observed and calculated data can be easily tracked down to erroneous low wavenumber content in the model, which generally relates to the quality of the initial model. What could be referred to as a "good" initial model will predict the data within less than half a period, allowing the L 2 -based FWI scheme to converge toward the global minimum. A degradation in the initial model's quality with erroneous low-wavenumber background velocities will likely translate into cycle-skipped arrivals, leading to artifacts during the reconstruction of the V P model. This illustrates the need for a good initial model already

n n + 1 n -1 n Time (s ) n + 1 n n -1 T /2 T /2 n + 1 n -1
Figure 1.4: Schematic example of the cycle-skipping/phase-ambiguity issue on sinusoidal signals. As soon as the initial shift is larger than half a period, the signal's fit using a least-squares distance is performed up to one or several phase shifts. One may try to fit the n + 1 dashed wriggle of the top signal with the n continuous wriggle of the middle signal moving in the wrong direction. The bottom dashed signal predicts the n wriggle in less than a half-period leading to a correct updating direction (from [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF]. inside the global valley when starting the FWI process.

Non-convexity in FWI

The link between the initial model "quality" and correct FWI reconstruction can be easily illustrated, as shown in a canonical test based on a spherical inclusion with a cross-hole acquisition (sources on the left, receivers on the right). The target is a simple square model of 1 by 1 km, of homogeneous V P = 1300 m.s -1 , with a spherical inclusion in the center of 100 m radius with V P = 1700 m.s -1 (Figure 1.6 (a)). Two initial models are used, both being completely homogenous V P , with one at the correct background velocity of the true model (1300 m.s -1 ), while the second one is at 1700 m.s -1 (respectively (b) and (c) in Figure 1.6). The data are generated with acoustic modeling using a Ricker wavelet with a central frequency of 3 Hz. The difference between d obs and the two d cal (one for each initial model) is presented in Figure 1.7. It clearly shows that the 2nd initial model, containing a faster V P , makes the signal arrive "too soon" compared to the observed data, creating out of phases arrivals by more than half a period. The results obtained with FWI starting from each initial model are presented in Figure 1.8. With no surprise, the reconstruction is correct for initial model 1, while results from initial model 2 present strong high-velocity artifacts, with no reconstruction of the spherical inclusion and no correct background velocity. This exhibits the limitation of L 2 based FWI for inaccurate initial model.

FULL-WAVEFORM INVERSION

Classical approaches

The ideas developed here are solutions to successfully apply L 2 -based FWI, relying on data and initial models available. This is often the combination of the folllowing elements.

Starting from a good initial model

The first condition to avoid convergence toward a local minimum is to start FWI from an initial model m 0 , which is inside the global minimum basin of attraction. Traditionally, FWI relies on using a good enough initial model that contains a correct approximation of the low wavenumber content of the model. Usually, such models come from reflection tomography or stereotomography [START_REF] Lambaré | Stereotomography[END_REF]. Unfortunately, these methods are based on travel-time or reflected events picking, which is a time-consuming task that requires extensive human expertise. Moreover, picking can introduce errors and uncertainties which will be propagated through FWI. This again links to the famous "gap" as we attempt to bring higher wavenumber content from tomography methods.

Data hierarchy

Complementary to the requirement of a good initial model, it is often required to select "smartly" the data for inversion, following an approach known as data hierarchy [START_REF] Bunks | Multiscale seismic waveform inversion[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] Shipp | Two-dimensional full wavefield inversion of wide-aperture marine seismic streamer data[END_REF][START_REF] Wang | Reflection seismic waveform tomography[END_REF][START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. As low-frequency data are less subject to cycle-skipping (as illustrated in Figure 1.9), the idea is Figure 1.10: Illustration of frequency-continuation approach used in FWI; here, on a shallow near-surface target of 30 m (from [START_REF] Irnaka | 3D elastic full waveform inversion for subsurface characterization: Study of a shallow seismic multicomponent field data[END_REF]. FWI is started on a 3 to 25 Hz frequency band, then higher frequency are progressively introduced with an upper limit of 45 Hz then finally 65 Hz. A clear improvement in resolution (smaller wavenumber content) is observed for higher frequency data.

to start the inversion using the lowest available frequencies while restricting higher frequency. Then, higher frequencies are progressively introduced following a multi-scale approach as shown in Figure 1.10 ( [START_REF] Sirgue | Efficient waveform inversion and imaging : a strategy for selecting temporal frequencies[END_REF].

Another level of data hierarchy consists of restricting the temporal/offset selection of the data to reduce the number of propagated wavelengths interpreted. By doing so, the shallow part of the reconstructed models can be recovered first, while later on, deeper updates are progressively introduced [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] Sears | Elastic full waveform inversion of multicomponent ocean-bottom cable seismic data: Application to alba field, u. k. north sea[END_REF][START_REF] Kamei | On acoustic waveform tomography of wide-angle obs data -strategies for preconditioning and inversion[END_REF][START_REF] Górszczyk | Toward a robust workflow for deep crustal imaging by FWI of OBS data: The eastern nankai trough revisited[END_REF]. We use this strategy for the field-data application in this thesis, and an illustration can be found in Figure 3.11. Practically, these two levels of data hierarchy are often used together to obtained successful FWI applications on field data, no matter what other methods are used to tackle the non-linearity of FWI.

FULL-WAVEFORM INVERSION

Reformulating the FWI problem

Nonetheless, even if the conditions details previously should allow successful FWI applications, they are not easy to gather or can even be insufficient to tackle some specific cases (such as sub-salt imaging). Therefore, finding other solutions to mitigate the non-linearity of FWI has been the motivation for a large number of studies in the past decades. Two main lines of investigation can be identified, which both boil down to reformulating the FWI problem.

Extension strategies

One way of seeing the non-linearity of FWI is that FWI is based on a reduced-space approach used to solve the PDE-constrained optimization problem, which is non-linear. Therefore, to mitigate this non-linearity, adding artificial degrees of freedom to the problem has been proposed. This corresponds to the first line of investigation and regroups methods often referred to as "extension strategies". The philosophy behind this is to artificially introduce supplementary degrees of freedom to the FWI problem in the hope to match data in early iterations (and therefore reduce cycle-skipping sensitivity). To ensure convergence, these artificial degrees of freedom are constrained to converge to zero to recover a correct and physical velocity reconstruction.

These methods have been historically developed with model-based extension, known as migration velocity analysis (MVA) methods [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF]. MVA is based on scale separation assumption: subsurface parameters can be decomposed into two parts, a smooth macro-velocity background model and a sharp high-wavenumber reflectivity model. Using this scale separation, MVA introduces artificial degrees of freedom in the reflectivity model. The goal is to update the macro-velocity model to focus the energy of the extended reflectivity model at zero in the artificial dimension. Mathematical analysis shows that in a transmission regime, the new extended problem of FWI asymptotically converges toward a travel-time tomography problem known to be convex. Unfortunately, MVA suffers from two main drawbacks. The first one relates to the high computational cost of building extended reflectivity hypercubes. The second one corresponds to a more fundamental difficulty and is related to the difficulty faced by MVA to handle properly complex data with multi-arrival and multiple reflections [START_REF] Cocher | Seismic Iterative Migration Velocity Analysis: two strategies to update the velocity model[END_REF].

More recently, as an attempt to solves the issue affecting MVA, a new class of extension strategies has emerged, with a different way of introducing artificial degrees of freedom, for example, at the source level with source extension strategy. These methods are referred to as matched source waveform inversion (MSWI) (Huang et al., 2018a,b). These technics have shown promising results in 2D synthetic cases using frequency-domain FWI. Still, their implementation for time-domain FWI is under development [START_REF] Wang | Full-waveform inversion with the reconstructed wavefield method[END_REF][START_REF] Aghamiry | Accurate and efficient wavefield reconstruction in the time domain[END_REF]. The reason behind the preference for frequency-domain FWI for MSWI boils down to the fact that MSWI requires the solution of a square wave propagation problem that is easily obtained through the harmonic wave equation based on a direct solver, which is much more difficult to obtain for time-domain time-stepping algorithms. This represents a disincentive to 3D field data application as time-domain FWI remains the methods of choice for this kind of targets (frequency-domain FWI is for now limited to moderate-size targets as there is a lack of scalability of the direct solvers used to solve harmonic equations, see [START_REF] Li | On the comparison of MUMPS and STRUMPACK for 3D frequency-domain elastic wave modeling[END_REF] for a recent review on the status of direct solvers to solve large scale harmonic wave equation problems).

Non-convexity in FWI

MSWI methods can be seen as an equivalent of a class of methods that have been previously introduced to FWI: wavefield reconstruction inversion (WRI) method (van Leeuwen andHerrmann, 2013, 2016;[START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF]. A parallel can also be drawn with deconvolution approaches, with one of the most documented ones being adaptive waveform inversion (AWI) by [START_REF] Warner | Adaptive waveform inversion: Theory[END_REF], which can be recast as an MSWI technique. This is interesting as AWI can also be seen as an alternative misfit function (which will be the second line of investigation for reformulation of the FWI problem), indicating that the distinction between the two families is not as straightforward as one could think.

Finally, we can note that another type of extension strategy can be made with the extension being made at the receiver location [START_REF] Métivier | Extension strategy for time-domain full waveform inversion using receiver relocalization[END_REF]. This new formulation bypasses the needs for square wave propagation operators and is therefore directly applicable to time-domain FWI with reasonable computational cost. This method has, however, not yet been tested on field data so far.

Alternative misfit functions

Another way of seeing the non-linearity of FWI relates to how the data is interpreted. As low-wavenumber content of the initial model mainly generates kinematic shifts in the data (time-shifts), the non-convexity of the least-squares (L 2 ) misfit function with respect to timeshifts partially explains the non-linearity of FWI [START_REF] Jannane | Wavelengths of Earth structures that can be resolved from seismic reflection data[END_REF]. This non-convexity with respect to time-shifts of the L 2 norm is the cause of the well-documented cycle-skipping issue.

Hence, the second line of investigation to mitigate the non-linearity of FWI is to reformulate FWI by replacing the L 2 misfit function used to compare the observed and calculated data with a more convex alternative. This idea is at the origin of a large variety of more-convex alternative misfit functions that have been proposed in the literature. This will be the subject of the analysis performed in the next chapter. This chapter directly comes from the published paper [START_REF] Pladys | On cycle-skipping and misfit functions modification for full-wave inversion: comparison of five recent approaches[END_REF].

Introduction abstract

Full waveform inversion, a high-resolution seismic imaging method, is known to require sufficiently accurate initial models to converge toward meaningful estimations of the subsurface mechanical properties. This limitation is due to the non-convexity of the least-squares distance with respect to kinematic mismatch. We propose a comparison of five misfit functions promoted recently to mitigate this issue: adaptive waveform inversion, instantaneous envelope, normalized integration, and two methods based on optimal transport. We explain which principles these methods are based on and illustrate how they are designed to better handle kinematic mismatch than a least-squares misfit function. By doing so, we can exhibit specific limitations of these methods in canonical cases. We further assess the interest of these five approaches for application to field data based on a synthetic Marmousi case study. We illustrate how adaptive waveform inversion and the two methods based on optimal transport possess interesting properties, making them appealing strategies applicable to field data. Another outcome is the definition of generic tools to compare misfit functions for full-waveform inversion.

Introduction

Full waveform inversion (FWI) is a high-resolution seismic imaging method dedicated to reconstructing the mechanical properties of the subsurface [START_REF] Devaney | Geophysical diffraction tomography[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain, part II: Fault delineation in sediments using crosshole data[END_REF][START_REF] Plessix | Full waveform inversion of a deep water ocean bottom seismometer dataset[END_REF][START_REF] Raknes | Three-dimensional elastic full waveform inversion using seismic data from the sleipner area[END_REF][START_REF] Górszczyk | Toward a robust workflow for deep crustal imaging by FWI of OBS data: The eastern nankai trough revisited[END_REF]. It is formulated as an iterative process based on minimizing a function measuring the misfit between observed and calculated data over a space of model parameters describing the subsurface. The resolution improvement FWI can procure, compared with standard tomography methods, is used to significantly improve depth-migration images or even produce directly interpretable quantitative estimates of the subsurface mechanical properties [START_REF] Shen | High-resolution full-waveform inversion for structural imaging in exploration[END_REF]. FWI is applied at multiple scales, from global and regional scales in seismology to exploration scale for the oil & gas industry, and even, more recently, at near-surface scale for geotechnical applications. A thorough review of FWI and its applications can be found in [START_REF] Virieux | An introduction to Full Waveform Inversion[END_REF].

FWI suffers from a significant shortcoming in its classical formulation: the non-convexity of the least-squares (L 2 ) misfit function on which it is conventionally based. This non-convexity of the misfit function is an issue because the iterative process on which is based FWI is a local optimization algorithm. Standard size for realistic applications makes global optimization strategies beyond modern high-performance computing platforms current and predictable capabilities. Therefore, if the initial model used is too far away from the global minimum, FWI converges toward a potentially non geologically informative local minimum. This constraint leads to the need for an accurate enough initial model to ensure convergence toward the global minimum of the misfit function.

In a physical sense, the non-convexity of the L 2 misfit function is associated with a phenomenon known as cycle-skipping. It appears when the calculated data are shifted (in time) from more than half a period (corresponding to the signal dominant frequency) compared to the observed data. If the time-shift between observed and calculated data is larger than half a period, the minimization of the L 2 norm between the two signals will "skip" a phase and align the two signals on the closest phase (hence the name, cycle-skipping). This ambiguity translates into an erroneous reconstruction of the velocity model [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF].

ON CYCLE-SKIPPING AND MISFIT FUNCTIONS MODIFICATION FOR FULL-WAVE INVERSION: COMPARISON OF FIVE RECENT APPROACHES

This limitation of FWI has been documented since its origin [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: numerical results[END_REF]. To address this limitation in practical cases, the workflow generally relies on data hierarchy [START_REF] Bunks | Multiscale seismic waveform inversion[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] Shipp | Two-dimensional full wavefield inversion of wide-aperture marine seismic streamer data[END_REF][START_REF] Wang | Reflection seismic waveform tomography[END_REF][START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. The historical approach consists in interpreting first the lowest frequency available (around 2 to 4 Hz for seismic exploration targets), then progressively introducing higher frequency data, following a multi-scale approach [START_REF] Sirgue | Efficient waveform inversion and imaging : a strategy for selecting temporal frequencies[END_REF]. The lowest frequencies are, by definition, less subject to cycle-skipping. The second level of data hierarchy can then be defined by playing on temporal and/or offset selection of the data. The idea is to reduce the number of propagated wavelengths that are interpreted simultaneously, hence reducing the risk of cycleskipping. In practice, this second level corresponds to first reconstructing the near-surface and progressively introducing deeper updates referred to as layer stripping approach.

Successful practical applications at the exploration scale often rely on the conjunction of these approaches as well as the design of an accurate initial starting model, obtained, for instance, through reflection tomography or stereotomography [START_REF] Lambaré | Stereotomography[END_REF]. Nonetheless, the conditions detailed previously to obtain a satisfactory FWI result are not always gathered. For instance, low-frequency data around 2 to 4 Hz are not always available or of sufficient quality. Moreover, obtaining low-frequency can increase the cost of acquisition, or can sometimes not be physically possible, or can even compromise the quality of the high frequency needed to obtain a very high resolution. Accurate initial model building can also be a time-consuming and challenging task requiring strong human expertise as it generally relies on tomography methods based on travel-time or reflected event picking. It also relies on prior information coming from geology or well logs; all of these require human expertise. This makes FWI less robust and reduces its range in terms of applications.

Mitigating the sensitivity to initial model quality has been the motivation for a large number of studies in the past decades. Two main lines of investigations can be identified, both leading to the reformulation of the conventional least-squares FWI problem.

Considering the first line, we regroup methods that can be cast under the frame of "extension strategies". It is not our purpose to give an extensive overview of these methods here, but we try to sketch their main ingredients. The philosophy of extension strategies consists in introducing supplementary degrees of freedom to the FWI problem, which can match the data in the early iterations of the FWI process to avoid cycle-skipping. Relaxing iteratively the use of these artificial degrees of freedom should lead to a correct subsurface model estimation.

Historically, these methods derive from migration velocity analysis (MVA) [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF]. MVA relies on the scale separation assumption. The subsurface parameters to recover are decomposed as a smooth macro-velocity model and a high wavenumber content reflectivity model. Artificial degrees of freedom are introduced at the reflectivity level by introducing an extra dimension on offset, subsurface offset, or time-lag. The MVA problem is formulated as the iterative update of the macro-velocity model to focus the energy of the "extended" reflectivity model at zero in the artificial dimension. These methods have benefited from in-depth mathematical research work, leading to a clear understanding of their foundations, thanks to the theory of pseudo-differential operators. However, their application to field data is still limited, mainly because of two issues. First, the repeated construction of high-dimensional reflectivity cubes is computationally demanding. Second, the macro-velocity model construction through MVA is complicated as soon as complex data with multi-pathing and multiple reflections are considered.

Introduction

More recently, another class of extension strategies has emerged. As opposed to model space extension, the artificial degrees of freedom are introduced at the source level, following a source extension strategy (Huang et al., 2018a;[START_REF] Van Leeuwen | Mitigating local minima in full-waveform inversion by expanding the search space[END_REF]. These methods have shown interesting promises in 2D synthetic case studies. However, their application to 3D field data seems still limited, mainly because of the difficulty of applying these methods in the time-domain. Current solutions either rely on relatively crude approximations [START_REF] Wang | Full-waveform inversion with the reconstructed wavefield method[END_REF] or on a sophisticated iterative solution, which increases the computational cost of the approach significantly [START_REF] Aghamiry | Accurate and efficient wavefield reconstruction in the time domain[END_REF].

The second investigation line relies on reformulating the FWI problem using an alternative measure of the distance between observed and calculated data, namely a different misfit function. A large variety of approaches have been proposed on this framework. The first proposed along this line is to use cross-correlation measurements [START_REF] Luo | Wave-equation traveltime inversion[END_REF], a strategy later revisited by van [START_REF] Van Leeuwen | A correlation-based misfit criterion for waveequation traveltime tomography[END_REF]. The idea behind this is that cross-correlation should give access to the time-shifts between synthetic and observed traces. A misfit function based on the minimization of these time-shifts, resembling a tomography misfit function, should thus be less prone to cycle-skipping. The original approach of [START_REF] Luo | Wave-equation traveltime inversion[END_REF] was labeled as "wave equation tomography" strategy.

However, when seismic traces contain multiple seismic events, the cross-correlation measurement might fail to give a correct estimation of a potential time-shift. This is why deconvolution based approaches have been later promoted, first by [START_REF] Luo | A deconvolution-based objective function for wave-equation inversion[END_REF], then improved by [START_REF] Warner | Adaptive waveform inversion: Theory[END_REF]. The latter approach has been labeled as "adaptive waveform inversion" (AWI) and is based on a normalized deconvolution of the synthetic and observed seismic traces. It has shown very interesting properties both on synthetic and field data. The deconvolution of the traces yields a Wiener filter, which is then normalized and serves as an input for the misfit function. The misfit function penalizes the energy of the filter away from a bandpass Dirac filter, which would have been obtained in the correct subsurface model. Note that AWI shares some similarities with the extended source approach and can indeed be recast in the frame of these methods (Huang et al., 2018a). This indicates that the separation between extended methods and misfit function reformulation methods is not as watertight as one could think. Nevertheless, it is useful to draw a landscape of the investigations around the cycle-skipping issue in FWI.

Another family of misfit function modifications relies on transforming the signal itself prior to comparison through a least-squares distance. Extracting the instantaneous phase and envelope [START_REF] Fichtner | Theoretical background for continental-and global-scale full-waveform inversion in the time-frequency domain[END_REF][START_REF] Bozdag | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF] has been successfully used in seismology. The goal of the instantaneous phase is to avoid amplitude prediction issues, as earthquake source and receiver calibration are significant challenges in seismology. The use of the envelope to mitigate the cycle-skipping issue has also been developed in the framework of seismic exploration [START_REF] Wu | Seismic envelope inversion and modulation signal model[END_REF]). An interesting alternative consists of using a normalized integration of the signal, namely the cumulative distribution of the traces. This approach has been promoted by [START_REF] Donno | Estimating the background velocity model with the normalized integration method[END_REF].

Finally, optimal transport distances have also been promoted to derive alternative misfit functions for FWI. The motivation is to benefit from the convexity of the optimal transport distance with respect to translation and dilation, which provides a misfit function convex with respect to time-shifts, this being a good proxy for convexity with respect to seismic velocities [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF][START_REF] Métivier | Optimal transport for mitigating cycle skipping in full waveform inversion: a graph space transform approach[END_REF]. The main difficulty in applying optimal transport in the framework of FWI is that the optimal transport theory is developed to compare

ON CYCLE-SKIPPING AND MISFIT FUNCTIONS MODIFICATION FOR FULL-WAVE INVERSION: COMPARISON OF FIVE RECENT APPROACHES

probability distributions, therefore positive functions with the same total integral. Seismic data do not fulfill this assumption.

To overcome this difficulty, different options have been promoted. For instance, one can rely on a prior transformation of the signal, such as extraction of positive and negative parts, squaring the data, affine scaling, exponential transform, softmax transform [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF][START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport norm[END_REF][START_REF] Yang | Application of optimal transport and the quadratic Wasserstein metric to full-waveform inversion[END_REF][START_REF] Yang | Analysis of optimal transport and related misfit functions in full-waveform inversion[END_REF]. This has been shown effective in some synthetic cases. However, relevant seismic information might be lost in the process of these transformations.

One solution is to rely on a specific optimal transport distance, which can be extended to comparing non-positive data. This is the Kantorovich-Rubinstein optimal transport (KROT) approach, which has been promoted in Métivier et al. (2016c,a,b), and which has been successfully applied to 3D synthetic elastic data [START_REF] He | Land seismic multi-parameter full waveform inversion in elastic VTI media by simultaneously interpreting body waves and surface waves with an optimal transport based objective function[END_REF] as well as to field data [START_REF] Poncet | Fwi with optimal transport: a 3D implementation and an application on a field dataset[END_REF][START_REF] Messud | Multidimensional optimal transport for 3d FWI: Demonstration on field data[END_REF][START_REF] Sedova | Acoustic land full waveform inversion on a broadband land dataset: the impact of optimal transport[END_REF]. One interest of this approach is its ability to account for lateral coherency in 2D or 3D shot gathers. One shortcoming is that, even if the valley of attraction is wider, compared with the L 2 approach, the convexity property of the optimal transport distance with respect to time-shifts is lost.

Another option has been promoted more recently. Considering each discrete seismic traces as point clouds and computing the optimal transport distance between synthetic and observed points clouds provide a new distance measurement. This specific optimal transport problem can be cast as a linear assignment problem, for which efficient solvers exist, for point clouds containing a few hundred to thousands of points, a situation we encounter for realistic scale exploration case studies [START_REF] Métivier | Optimal transport for mitigating cycle skipping in full waveform inversion: a graph space transform approach[END_REF][START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF]. The benefit of this graph-space optimal transport (GSOT) strategy is its ability to recover the convexity with respect to time-shifts. Compared with the KROT approach, GSOT is a trace-by-trace strategy that does not make it possible to account for lateral coherency. GSOT has been successfully applied to 3D synthetic and field data (He et al., 2019a;[START_REF] Pladys | Assessment of optimal transport based FWI: 3d OBC valhall case study[END_REF][START_REF] Li | Joint FWI for imaging deep structures: A graphspace OT approach[END_REF][START_REF] Górszczyk | Mitigating the nonlinearity of the crustal scale full waveform inversion through the graph space optimal transport misfit function[END_REF].

As can be seen, numerous investigations motivated by the inherent ill-posedness of the FWI problem have been lead in parallel. To our knowledge, no cross-comparison has been proposed so far, which is undoubtedly a lack. The first motivation of this study is to start developing tools that could be used to benchmark different FWI strategies. However, beyond a simple comparison of FWI strategies, we would like to highlight specific characteristics that an ideal misfit function should satisfy to render the FWI problem less ill-posed. Cycle-skipping is certainly an issue, but we also show that other criteria than robustness with respect to cycle-skipping should be considered, such as:

• sensitivity to the signal polarity;

• applicability in the framework of complex/multi-arrival data;

• number of tuning parameters and sensitivity to these parameters;

• sensitivity to wrong amplitude prediction and inaccurate wavelet estimation.

To illustrate these properties, we select a series of synthetic case studies of increasing complexity, from time-shifted Ricker traces to a realistic Marmousi II case study (not in inverse crime settings). We restrict our attention to five misfit functions, which have been promoted recently and have shown promising results: adaptive waveform inversion (AWI), instantaneous envelope 2.2 General FWI framework and misfit function formulation (IE), normalized integration method (NIM), KROT, and GSOT. We consider extended space strategies out of the scope of this study to keep it reasonably simple, and also because, as stated before, we consider that alternative misfit strategies have shown more promising results than extended space strategies so far in terms of practical applications. The tests that we develop here could, however, be used to benchmark extended space strategies also.

General FWI framework and misfit function formulation

The comparison between misfit functions is made simple by the FWI formalism (reviewed in the following section), more precisely by the adjoint state strategy used to compute the gradient at each iteration of the minimization loop. However, let us recall the main result: a modification of the misfit function results only in modifying the adjoint source. Therefore, implementing different misfit functions in the same FWI code can be done directly by isolating misfit function evaluation and adjoint source computation in different subroutines.

General framework

The FWI problem can be written as The solution of the minimization problem 2.1 is computed through local optimization following the iteration

min m f [m] = F (d cal [m], d obs ) , ( 2
m k+1 = m k + α k ∆m k (2.4)
starting from an initial guess m 0 . In eq. 2.4, α k is the steplength, which should satisfy the Wolfe criterion [START_REF] Nocedal | Numerical Optimization[END_REF], and ∆m k is the descent direction, given by

∆m k = -P [m k ]∇f [m k ] , (2.5)
where ∇f (m k ) is the gradient of the misfit function f [m] and P [m k ] a preconditioner approximating the inverse Hessian operator

P [m k ] H[m k ] -1 , H[m k ] = ∇ 2 f [m k ]. (2.6)
Following the adjoint state strategy [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF], the gradient is given by

∇f [m] = ∂A ∂m u, λ , (2.7)

ON CYCLE-SKIPPING AND MISFIT FUNCTIONS MODIFICATION FOR FULL-WAVE INVERSION: COMPARISON OF FIVE RECENT APPROACHES

where (., .) is the Euclidean scalar product in the wavefield space, and λ is the adjoint field, solution of the adjoint equation

A(m) T λ = s , (2.8)
where s is the generic adjoint source, given by

s = -R T ∂F ∂d cal .
(2.9)

Note that in the case of the L 2 norm, we recover immediately that

s = -R T (Ru[m] -d obs ) , (2.10)
i.e. the adjoint source is equal to the residual (difference between observed and calculated data).

Next, we review the formulas for the five misfit functions selected in this study, as well as their corresponding adjoint sources. For convenience, we will introduce the distance measurement function associated with each strategy for a single source/receiver couple, except for the KROT strategy. The calculated and observed data will be denoted by d cal (t) and d obs (t) unless stated otherwise. Except for KROT, the final misfit function is built as a sum over each source/receiver couple of this distance measurement function, and by linearity, the resulting adjoint source is also obtained by summation.

Adaptive waveform inversion

We give here the AWI formalism. We have

F AW I (d cal , d obs ) = T 0 |P(τ )w(τ )| 2 dτ T 0 |w(τ )| 2 dτ , (2.11)
where w(t) is the Wiener filter which either transforms the calculated d cal (t) into the observed data d obs (t) (forward AWI) or the opposite way around (reverse AWI). Both implementations are discussed in [START_REF] Warner | Adaptive waveform inversion: Theory[END_REF]. Also, the computation of w(t) can be implemented either in the time-domain or the frequency-domain. In both cases, a water level ε is required to stabilize the deconvolution operation.

The role of the function P(τ ) is to penalize energy at non-zero time lag. There are several possibilities to define this penalty function. Here we focus only on a Gaussian formulation defined as

P(τ ) = e -τ 2 /σ 2 , (2.12)
where σ is a tuning parameter controlling the width of the Gaussian function away from 0 time-lag. This σ tuning parameter is defined in seconds and corresponds to the maximum expected time-shift between the observed and calculated data.

In the case of a frequency-domain reverse AWI implementation, the adjoint source for a single-trace reads In eq. 2.15, ζ is a user-defined damping ratio, ranging from 10 -2 to 10 -5 in our experiment.

∂F AW I ∂d cal = dτ (P(τ ) -2F (d cal , d obs )) w(τ )p(t + τ ) dτ w 2 (τ ) , ( 2
A large ζ will help when trying to tackle large time-shift, with a "smoothing/regularizing" effect. Large ζ is also required if there is noise on the data. A smaller ζ will help preserve small features present in the signal. In terms of computational cost, the overhead associated with the computation of the Wiener filter is negligible, and the AWI strategy can be easily implemented.

Instantaneous envelope

The separation of the phase and envelope information of the signal relies on the use of the analytical function defined as follows. For a given time signal d(t), the analytical signal d(t) is defined as

d(t) = d(t) + iH[d(t)] , (2.16)
where H is the Hilbert function which can be defined in the time domain as

H[d(t)] = 1 π P +∞ -∞ d(τ ) t -τ dτ ,
(2.17)

where P stands for the Cauchy principal value. Practically, we do not use the time formulation of the Hilbert function, but rather a frequency domain formulation that gives us the analytical signal in a three-step approach [START_REF] Marple | Computing the discrete-time" analytic" signal via fft[END_REF]):

• Compute the Fourier transform of d(t) using an FFT

• Change the negative frequency to zeros

• Compute the inverse Fourier transform

This directly gives us access to the analytical signal and, by extension, to the Hilbert transform by taking its imaginary part

H[d(t)] = I[ d(t)] .
(2.18)

The analytical signal allows to separate the signal as the combinaison of the instantaneous phase φ(t) and the instantaneous envelope E(t): t) .

d(t) = E(t)e iφ(
(2.19)

Thus, the intantaneous enveloppe E(t) can be simply defined as:

E(t) = R[ d(t)] 2 + I[ d(t)] 2 .
(2.20)

ON CYCLE-SKIPPING AND MISFIT FUNCTIONS MODIFICATION FOR FULL-WAVE INVERSION: COMPARISON OF FIVE RECENT APPROACHES

We can define a new distance-measurement function using instantaneous envelope as

F IE (d cal , d obs ) = 1 2 T 0 |E cal (t) -E obs (t)| 2 dt , (2.21)
where E cal and E obs are instantaneous envelopes of the calculated and observed data respectively. Following [START_REF] Bozdag | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF], the adjoint source is defined as:

∂F IE ∂d cal = (E cal (t) -E obs )d cal (t) E cal (t) + ε -H (E cal (t) -E obs )H(d cal (t)) E cal (t) + ε , (2.22)
with ε a water level defined as ε = (max

t E obs (t))ζ.
(2.23)

Contrary to AWI, in the following experiments, ζ is fixed and taken at ζ = 10 -5 for IE. We have verified that the results with IE are not sensitive to this choice.

The instantaneous envelope misfit formulation is straightforward to implement thanks to the algorithm from [START_REF] Marple | Computing the discrete-time" analytic" signal via fft[END_REF]. No tuning parameter is required, and the computation cost overhead is negligible. 2013) consider the least-squares difference between the cumulative distributions Q obs and Q cal . For a given time signal d(t), its normalized cumulative distribution Q(t) is defined by

Normalized Integration Method

Q(t) = t 0 d(τ ) 2 dτ T 0 d(τ ) 2 dτ .
(2.24)

The NIM misfit function thus relies on the distance measurement

F N IM (d cal , d obs ) = 1 2 T 0 |Q cal (τ ) -Q obs (τ )| 2 dτ , (2.25)
where Q cal (t) and Q obs (t) are the cumulative distributions associated with d cal (t) and d obs (t) respectively.

The corresponding adjoint source is

∂F N IM ∂d cal (t) = 2d cal (t) T 0 Q cal (t) T t (Q cal (τ ) -Q obs (τ ))dτ - T 0 Q cal (τ )(Q cal (τ ) -Q obs (τ ))dτ .
(2.26)

The NIM implementation is straightforward and does not require any tuning parameters.

2.2 General FWI framework and misfit function formulation

Kantorovich-Rubinstein optimal transport distance

In the frame of the KROT approach, we consider the data as a function of both time and receiver position, such that we denote the calculated and observed data as d cal (x r , t) and d obs (x r , t) respectively.

The KROT is based on a particular instance of optimal transport distance, namely the 1-Wasserstein distance. It can be applied to non-positive data, provided mass conservation is satisfied i.e. (2.27)

For a given shot in seismic data, this corresponds to the summation over each trace of the mean value in time of the trace. We consider this mean value is equal to 0 (this is the zerofrequency noise, which is usually removed from the data prior to inversion). Therefore the mass conservation assumption is satisfied for seismic data.

On this basis, the KROT distance can be written as

F KROT (d cal , d obs ) = max ϕ∈Lip 1 xr T 0 ϕ(x r , t) d cal (x r , t) -d obs (x r , t) dx r dt , (2.28)
where Lip 1 is the set of 1-Lipschitz functions for the 1 distance

Lip 1 = ϕ(x r , t), |ϕ(x r , t) -ϕ(x r , t )| < |x r -x r | + |t -t | .
(2.29)

The adjoint source is then given by

∂F KROT ∂d cal = ϕ(x r , t), (2.30)
where

ϕ(x r , t) = arg max ϕ∈Lip 1 xr T 0 ϕ(x r , t) (d cal (x r , t) -d obs (x r , t)) dx r dt . (2.31)
Compared with previous misfit functions, the final misfit is obtained here by summation over shot gather, and not a summation over source/receiver couples (not a trace-by-trace approach).

From the above equations, we see that the computation of the KROT misfit function and its corresponding adjoint source requires solving a constrained maximization problem per shot gather. Details on how to solve this problem are given in [START_REF] Métivier | An optimal transport approach for seismic tomography: Application to 3D full waveform inversion[END_REF]. The proximal splitting algorithm ADMM is used [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] and the resulting algorithm has complexity in O(N log N ), where N = N r × N t with N r the number of receivers and N t the number of time samples. Compared with the previous misfit functions, the computational cost overhead is non-negligible. Tuning parameters will be associated with a prior scaling of the data to make its maximum amplitude close to 1, and the number of iterations required to solve the constrained maximization problem.

ON CYCLE-SKIPPING AND MISFIT FUNCTIONS MODIFICATION FOR FULL-WAVE INVERSION: COMPARISON OF FIVE RECENT APPROACHES

Graph space optimal transport distance

Let (t i , d(t i ), i = 1, . . . , N ) be the discrete graph of the time function d(t). This discrete graph is a point cloud containing N points. The GSOT distance measurement is formulated as

F GSOT (d cal , d obs ) = min σ∈S(N ) N i=1 c iσ(i) , (2.32)
where c ij is the L 2 distance between the points of the discrete graph of d cal and d obs , namely

c ij = |t i -t j | 2 + η 2 |d cal (t i ) -d obs (t j )| 2 , (2.33)
and S(N ) is the ensemble of permutations of (1 . . . N ). The function F GSOT corresponds to the 2-Wasserstein distance between the discrete graph of the calculated trace d cal (t) and the observed trace d obs (t).

The scaling parameter η in eq. 2.33 controls the convexity of the misfit function f GSOT with respect to time-shifts. In practice, we define it as

η = τ A , (2.34)
where τ is a user-defined parameter corresponding to the maximum expected time-shift between observed and calculated data in the initial model, and A is the maximum amplitude discrepancy between observed and calculated data.

The adjoint source of the misfit function f GSOT [m] is computed from ∂f GSOT ∂ cal using the adjoint-state strategy. We prove in [START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF] the following equality: denoting σ * the minimizer in eq. 2.32, we have

∂f GSOT ∂ cal = 2 d cal -d σ * obs , (2.35) where d σ * obs (t i ) = d obs (t σ * (i) ) .
(2.36)

In this sense, the GSOT approach can be viewed as a generalization of the L 2 distance: the adjoint source is equal to the difference between calculated and observed data at time samples connected by the optimal assignment σ * . Similarly, as the KROT approach, the solution of the problem 2.32 provides the information to compute both the misfit function and the adjoint source.

The numerical algorithm used to solve the linear assignment problem 2.32 is the auction algorithm [START_REF] Bertsekas | The auction algorithm for the transportation problem[END_REF]. For problems involving less than 1000 points, the auction algorithm is very efficient. In seismic exploration, Nyquist sampling yields traces containing a number of points within this order of magnitude. Consequently, we have designed an efficient numerical strategy, yielding lower computational overhead than the KROT approach. On 3D field data application, we observe 15 to 20% computation time increase for gradient computation on the lowest frequency bands compared with classical L 2 . This computational cost overhead decreases when the frequency band increases as the total complexity of the GSOT problem is O(ω 3 ), while the complexity of the wave propagation solver is in O(ω 4 ). For more details, the reader can refer to [START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF].

A simple convexity analysis based on time-shifted Ricker wavelets

Compared with previous approaches, the computational cost overhead is comparable with AWI, IE, and NIM while being lower than KROT. In terms of implementation, as for KROT, the solution of the assignment problem requires specific solvers, which makes the GSOT implementation less trivial than for AWI, IE, or NIM. In terms of tuning parameters, the more important parameter is the parameter τ , which controls the convexity of GSOT misfit function with respect to time-shifts.

A simple convexity analysis based on time-shifted Ricker wavelets

We start by investigating the convexity of the proposed misfit functions with respect to timeshifts. We fix a reference signal composed of one Ricker wavelet in the center, seen as the observed data. The calculated data is the same Ricker wavelet, shifted in time with a timeshift going from -1.5 s to 1.5 s. We compute the distance between the reference signal and the calculated signal using the five selected misfit functions, depending on the input time-shift.

Results are presented in Figure 2.1.

The results obtained here with alternative misfit functions might not reflect the performance of the algorithms with total accuracy, both in terms of computational efficiency and inversion results. Algorithms might not have been implemented in the most optimal way or in the way the original authors intended. Subtle choices of tuning parameters might improve the inversion results in some cases. However, the primary purpose of this comparison is to seek to understand how the data is interpreted within each of these strategies and how this affects the inversion results in each case. We intend to provide the reader with sufficient material to infer the main properties and philosophy behind the compared methods.

Let us first analyze the results obtained with L 2 waveform misfit, the reference for FWI. As expected, L 2 misfit displays a narrow basin of attraction, with local minima and a flat part for time-shift superior to 0.4 s. The local minimum appears when the time-shift is larger than 0.12 s, which corresponds to half the Ricker wavelet period. This validates that the L 2 misfit function presents low robustness for shifted-patterns, leading to cycle-skipping when signals are shifted by more than half a period. In such cases, L 2 misfit function does not guarantee convergence toward the global minimum.

We can now compare the selected alternative misfit functions to the L 2 misfit. From the obtained results, we can define two groups. The first one contains GSOT, AWI, and NIM, characterized by a large basin of attraction. The second group contains IE and KROT, characterized by a "slightly" larger basin of attraction than L 2 , but not as wide as the first group members.

Understanding why the first group members exhibit the convexity property is essential. Starting with GSOT, if the input parameters τ is correctly set to the maximum expected timeshift of 1.5 s, the convexity to shifted-patterns is expected as there is a direct link between the τ parameters and the width of the basin of attraction as shown in [START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF].

The same convexity property is observed with AWI. With σ set to 1.5 s, the results are satisfying with a large basin of attraction. Similarly as the τ parameter from GSOT, σ directly controls the convexity to shifted-patterns. Note that we use ζ = 10 -5 in this analysis as we predict signal with machine precision. INVERSION: COMPARISON OF FIVE RECENT APPROACHES Finally, to understand the robustness of the NIM approach, we display in Figure 2.2 the quantities Q obs and Q cal (for three time-shifts, -1.5 s, -0.1 s and in-phase). This makes visible the drastic modification of the signal shape induced by NIM. The NIM cost function boils down to be the area under the curve delimited by Q obs -Q cal . We see clearly that this area increases with time-shifts, illustrating the convexity to shifted patterns observed with NIM.

Moving to the second group, to understand why the IE misfit only slightly increases the width of the valley of attraction compared with L 2 , we display in Figure 2.3 the quantities E obs and E cal . Here we can observe the increase of temporal support of the signal induced by the envelope. This "broader" temporal support of the instantaneous envelope directly translates into the increase of the width of the valley of attraction as IE relies on a L 2 norm between E obs and E cal .

Finally, we present in Figure 2.4 the function ϕ(t) solution of the maximization problem defined in eq. 2.31, which defines the KROT distance, together with the residuals d obs (t)d cal (t). We can observe that when Ricker wavelets start to overlap at -0.3 s, we obtain a convexity that classical L 2 cannot achieve. This can be understood by looking at the function

ϕ(t) [d obs (t) -d cal (t)].
The area below the curve defined by this function corresponds to the KROT misfit function. This area remains constant as long as the two signals do not overlap and monotonically decrease as soon as the two signals overlap, reaching 0 at 0 time-shift.

On a second test, presented in Figure 2.5, we introduce a second Ricker wavelet that remains in phase. This test aims at validating the robustness to cycle-skipping when multiple arrivals are considered. From the results obtained, we observe that all misfit functions behave similarly as on the previous test except for AWI. In this case, the shape of the misfit function seems affected by oscillations near 0 time-shift, reducing the effective convexity to the one of classical L 2 formulation. This seems to be related to one of the potential issues of deconvolution based misfit function: the sensitivity to cross-talks between multiple events. To analyze this sensitivity of AWI to multi-arrivals, we display the Wiener filters together with the penalty function and the combination of both (Figure 2.6). In test B (where one wavelet is always in-phase), the Wiener filter presents a strong peak at 0 time-lag due to the in-phase arrivals. Because of finite frequency effect, it is not a Dirac delta function but a bandpass Dirac delta function. The oscillations of the bandpass delta function combine in a destructive/constructive manner when the two time-shift peaks (one for each Ricker wavelet) get closer to each other. These interferences are at the origin of the local minima observed.

A simple convexity analysis based on time-shifted Ricker wavelets

FWI tests on two canonical examples

This section attempts to assess the pros and cons of the selected alternative misfit functions on two schematic FWI tests, focusing on a different aspect of the information contained in a dataset. The first test focuses only on transmission with a cross-hole acquisition. The second test focuses mainly on reflection information. These two tests can be seen as a way of assessing if the proposed misfit function can improve the FWI robustness (cycle-skipping in transmission in the first test) while preserving the ability to correctly interpret reflection information (reflector positioning and imaging in the second test)

Both tests are performed in 2D using our 2D/3D time-domain acoustic modeling and inversion code in inverse crime settings (observed and calculated data are computed on the same grid, without noise introduced in the data). Besides, we use a constant density model and invert only for the P-wave velocity model. In both cases, the l-BFGS algorithm is used to minimize the misfit function, with FWI stopping criterion being a line search failure. The source wavelet is a Ricker wavelet with a central frequency f ref = 3 Hz. The gradient is smoothed using a Gaussian filter with horizontal and vertical correlation lengths equal to 0.3 times the local wavelength

λ loc (x, z) = 0.3 v P (x, z) f ref .
(2.37)

FWI Test 1: transmission configuration

Case study presentation

This first case study focuses on transmitted energy. The exact model is defined as a square of 1000 m sides with homogeneous V P = 1300 m.s -1 containing a spherical inclusion of 100 m radius in the center with V P = 1700 m.s -1 (Figure 2.7). The acquisition mimics a crosshole setting, with 96 sources on the left side of the model and 256 receivers on the right side.

The spacing is 10 m between sources and 3.8 m between receivers. The boundaries are all set to absorbing layers [START_REF] Bérenger | A perfectly matched layer for absorption of electromagnetic waves[END_REF] to avoid reflections and only focus on transmitted events. The relatively strong contrast between the background and the anomaly generates an identifiable diffraction pattern in the data. In this experiment, no preconditioning is applied to the gradient. The lower and upper V P bound constraints are respectively set to 1000 and 2500 m.s -1 .

We introduce three starting homogeneous models (Figure 2.7). The first is at the true model background velocity (1300 m.s -1 ). The second is at V P = 1700 m.s -1 , setting a challenging FWI problem as the starting model is as fast as the inclusion. The third case is even more challenging, with a starting homogeneous V P model at 1900 m.s -1 . FWI results are presented in Figure 2.8 with reconstructed V P at the final iteration. Figure 2.9 presents traces for a single source-receiver couple representing the shortest path through the spherical inclusion (straight horizontal path at 500 m depth). Traces are extracted from data generated in the true model, initial model, and final reconstructed model for all misfit functions.

FWI tests on two canonical examples

Results from initial model 1

We start the analysis with the "reference" initial model. As shown in Figure 2.9, this model does not generate cycle-skipping (arrivals in the true model are less than half a period away from the arrivals in the initial model). The objective is to retrieve the high-velocity spherical inclusion in the center of the model. As expected, the L 2 misfit function produces a correct result: the inclusion is retrieved correctly, and the final data are in phase with the true data. The vertical resolution is higher than the horizontal resolution as expected from the cross-hole configuration. This has a lateral smoothing effect on the reconstructed anomaly, which explains why its peak amplitude (around 1500 m.s -1 ) is lower than the amplitude of the true anomaly. The five selected misfit functions produce equivalently good results in this configuration. In all cases, the spherical anomaly is reconstructed with a similar resolution, and the data fit is equivalent. In terms of parameter settings, we choose here τ = 0.2 s for GSOT and σ = 0.2 s for AWI, a choice motivated by the absence of cycle-skipping. For AWI, we use ζ = 10 -2 in this transmission test (for all three models) to maximize the kinematic effects of AWI that work better when ζ is relatively high, which acts as a regularization effect.

Results from initial model 2

As can be observed in Figure 2.9, the second initial model generates clear cycle-skipping in the data. In this case, we expect the L 2 misfit function to fail in reconstructing the anomaly. Indeed, the L 2 fails to converge and reaches the boundary set for the inversion. The final synthetic trace does not match the observed trace. It is interesting to observe that four of the five selected misfit functions succeed in reconstructing the background and the anomaly and produce final synthetic traces in phase with the observed trace in this already quite challenging test. The only alternative misfit function that fails is KROT, which could be expected from the previous section (weak increase of robustness to cycle-skipping). AWI, IE, NIM and GSOT show that the increase in convexity procured by these formulations is enough here to make convergence achievable. The data-fit obtained with these methods is good in this case. In terms of tuning parameters, τ and σ are increased to 0.35 s for GSOT and AWI, according to the time-shift between the reference and the initial traces in the initial model.

Results from initial model 3

Finally, the initial model 3 generates an even more substantial cycle-skipping effect than model 2 (Figure 2.9). L 2 and KROT still fail to converge to the correct model, as it was already the case starting from model 2. IE starts to exhibit diagonal cycle-skipping artifacts associated with the longest source/receiver paths in this more challenging setting. This is expected from the time-shift convexity analysis performed before: IE robustness to cycle-skipping is limited. AWI also starts to exhibit artifacts close from the acquisition, while central anomaly is correctly reconstructed (with σ = 0.6 s). NIM and GSOT (with τ increased to 0.6 s) achieve a relatively satisfactory reconstruction of the background and anomaly, similar to the results obtained from the previous background models. 2.4.2 FWI Test 2: reflection configuration

ON CYCLE-SKIPPING AND MISFIT FUNCTIONS MODIFICATION FOR FULL-WAVE INVERSION: COMPARISON OF FIVE RECENT APPROACHES

Case study presentation

This second case study focuses on reflected energy. We consider two different true models, composed of a homogeneous background at 1500 m.s -1 and a velocity layer 100 m thick at 300 m depth (Figure 2.10). In the first case, the velocity of the layer is set to 1600 m.s -1 , while in the second case, the velocity of the layer is set to 1400 m.s -1 . The starting model is homogeneous at the correct background velocity of 1500 m.s -1 (Figure 2.10). The surface acquisition comprises 96 sources and 512 receivers located close to the surface at 42 m depth. The spacing is 20 m between sources and 3.8 m between receivers. We implement PML absorbing conditions on the bottom and lateral sides of the medium to mimic a medium of infinite extension in these directions and a free surface condition on the top of the model. A simple linear in-depth preconditioner is also applied to compensate for geometrical spreading effects and accelerate the convergence. The lower and upper V P boundaries for the inversion are respectively set to 1200 m.s -1 and 1800 m.s -1 .

This test analyzes how the reflected data is interpreted by FWI depending on the choice of misfit function. The difference between the two exact models is only the sign of the velocity change at the layer level: in one case, velocity increases; in the second case, it decreases. This induces a change of polarity of the reflected wave, as clearly visible in Figure 2.11. We want to identify how the different misfit functions are sensitive to this change of polarity. INVERSION: COMPARISON OF FIVE RECENT APPROACHES

Results analysis

The reconstructed models are presented in Figure 2.12. Traces from the observed and synthetic data in the initial and final models for zero offset couple (source and receiver at the same position) located in the middle of the acquisition are presented in Figure 2.13. For visualization purposes, we cropped over the reflection after the first arrival.

The L 2 results are coherent with the expectation, with a correct reconstruction of the layer in both cases. The L 2 norm is sensitive to amplitude variation and polarity and is expected to interpret reflected events correctly. As the background velocity is known, there is no cycleskipping in the initial model for the two target models. The data fit in both cases is perfect.

Together with L 2 misfit function, results obtained with IE, KROT, and GSOT are equivalently correct. This is expected from KROT and GSOT, which should behave similarly as L 2 when cycle-skipping does not occur. GSOT relies on τ = 0.2 s in this experiment. This is somehow more surprising from IE, as one could think that the polarity of reflected events might be lost in the envelope extraction process. However, this is not true in this case but might be due to the inverse crime settings we are using. There is indeed a subtle change in the envelope of the observed data between model 1 (positive layer) and model 2 (negative layer) (Figure 2.14), which is enough to guide the inversion in the right direction. However, this is probably possible only because the first arrival is correctly predicted. Small inaccuracies in predicting the first arrival might be enough to impede a correct reconstruction using IE.

The results obtained with NIM are less satisfactory. In particular, the reconstruction of the negative layer is altered by strong positive artifacts beneath the layer. In the opposite case, negative artifacts also pollute the reconstruction of the positive layer, although the strength of these artifacts seems weaker. Analyzing the data fit shows that NIM has difficulties reproducing the reflection pattern in both cases (Figure 2.13). Spurious oscillations appear, which can be associated with the positive artifacts observed on the model reconstruction. This lack of sensitivity to the polarity is somehow expected. From NIM formulation, Q obs should be more or less the same independently of model 1 or model 2 being used. This is illustrated in Figure 2.14, where Q obs (t) is presented for both models (positive and negative layers). We can observe that the difference between the two true models leads to a very marginal modification of Q obs compared to Q cal . This is likely the explanation of the difficulties faced by NIM in interpreting the reflected waves correctly.

Finally, the results obtained with AWI are incorrect for both the negative and positive V P anomaly. From observing the data, we can see that the direct waves exhibit a clear dominance in amplitude over the reflected events. Therefore, we expect the Wiener filter to be dominated by the direct waves and only show a small imprint of the reflected waves that are of small amplitude (≈ 1% of peak amplitude). This is illustrated by the Wiener filter shown in Figure 2.15 for both positive and negative layer models. They indeed present a main event around 0 lag, corresponding to the in-phase direct wave. Around 0.3 s, the imprint of reflected events is very weak but still visible in the Wiener filters. This motivates us to use a large σ = 1 s to maximize the information coming from the small reflected waves, together with a small ζ = 10 -5 . However, these settings do not make it possible to obtain satisfactory results with AWI. To have a deeper understanding of why AWI fails to reconstruct a proper V P model in this case, we perform a sensitivity analysis of the misfit function with respect to the value of V P in the layer. We compute the AWI misfit value (and L 2 misfit value for reference) between d obs and d cal (V P ). Here d obs corresponds to a shot gather in the center of the acquisition generated in the true model. d cal (V P ) corresponds to data generated in different models similar to the true model, with as only varying parameter the layer velocity (ranging from ±100 m.s -1 around the layer velocity of the true model). The results of this analysis are presented in Figure 2.16. The L 2 results are coherent with the expectation: the misfit function is convex with respect to the variation of the layer velocity and presents a minimum when the velocity of the layer used to generated d cal (V P ) is similar to the one of the true model used to generate d obs , so respectively 1400 and 1600 m.s -1 . For AWI, we observed that the minimum is not aligned with the correct velocity (1440 m.s -1 in the first case, 1610 m.s -1 in the second). This exhibits the loss of sensitivity of AWI in this case, explaining the failure of convergence of the FWI. Only reducing the σ below 0.04 s would make AWI behaves more like L 2 and converge to a result similar to the one of NIM. We do not think such a parameterization is interesting as it prevents the advantages introduced by AWI, which is improved convexity, and still introduces artifacts in the reconstructed model and computational overhead compared to L 2 . This second FWI test is a good illustration of the potential limitation that an alternative misfit function mainly focused on resolving time-shift could introduce. Here AWI and NIM have difficulties providing satisfactory results when the main arrival is correctly predicted. The point to point approaches that classical L 2 procures is here the "reference", making possible to fit the small perturbation properly following the main arrivals. AWI and NIM being more kinematic oriented, it is not surprising that this setup is challenging for such formulations. We design a synthetic case study using the Marmousi II P-wave velocity model (Figure 2.17) [START_REF] Martin | Marmousi2: An elastic upgrade for Marmousi[END_REF] to continue our analysis on a more realistic FWI configuration. We use a fixed spread surface acquisition model with 128 sources and 169 receivers. The source spacing is 132 m, and the receiver spacing is 100 m. The data is generated with a 4 Hz centered Ricker wavelet high-pass filtered to remove energy below 2 Hz (wavelet is visible in Figure 2.24). The recording time is set to 7 s. PML absorbing layers are used on the bottom and lateral sides of the model to mimic a medium of infinite extension in these directions, while a free surface condition is applied on top.

In the first case, referred to as "inverse crime inversion", we model the data in the constant density acoustic approximation and use the same grid for modeling and inversion to remain in the inverse crime settings. The mesh spacing is 25 m in this case.

In the second case, referred to as "more realistic inversion", we use the variable density Marmousi II model and a refined 10 m grid to generate the data. White noise bandpassed between 2 Hz and 10 Hz is added to the data to reach a signal to noise ratio of 15%. The inversion is done on a 25 m grid, using a density model derived from the initial V P model through Gardner's law [START_REF] Gardner | Formation velocity and density-the diagnostic basics for stratigraphic traps[END_REF] (Figure 2.18). A wavelet estimation is done before inversion. Performing the inversion in this more realistic framework, away from the usual inverse crime settings, makes it possible to assess the effects of incorrect amplitude prediction on the different misfit functions and better judge their usability toward field data applications.

The optimization is performed using the l-BFGS algorithm. The regularization of the gradient is defined as 0.3 of the local wavelength. Pseudo-hessian preconditioning is used [START_REF] Choi | Frequency-Domain Elastic Full Waveform Inversion Using the New Pseudo-Hessian Matrix: Experience Of Elastic Marmousi 2 Synthetic Data[END_REF]Yang et al., 2018a). The lower and upper V P boundaries for the inversion are respectively set to 1000 m.s -1 and 5200 m.s -1 . Inversion is performed without any frequency continuation approaches or other multi-scale strategies for all the misfit functions considered (including L 2 ).

Inverse crime inversion

Case study description

We rely here on two starting models. The first one, called S500 (Figure 2.17 c), is derived from the true V P model using a Gaussian smoothing with a correlation length of 500 m. This starting model preserves the long wavenumber content of the true model. The second one is a linearly increasing vertical 1D (Figure 2.17 d) model, ramping from 1500 m.s -1 at seabed to 4500 m.s -1 at depth. This initial model does not contain long-wavelength structures inherited from the true model. The data-fit obtained through this initial model (Figure 2.19 b) is affected by cycle-skipping. In comparison, the data-fit obtained with S500 is globally better (Figure 2.19 a), with more in-phase arrivals (especially on 0 to 3 km offset diving waves). In this model, the L 2 results give, at first order, a good reconstruction of the Marmousi model. However, we can observe on the left part that the horizontal layers are not correctly

ON CYCLE-SKIPPING AND MISFIT FUNCTIONS MODIFICATION FOR FULL-WAVE INVERSION: COMPARISON OF FIVE RECENT APPROACHES

reconstructed and present an up-shift (0 < x < 3 km), associated with a low-velocity anomaly on the shallow left part of the model (around x = 1 and z = 0.8 km). This corresponds to the part where strong reflections are generated. Because the background velocity is incorrectly predicted in the early iterations, the arrivals corresponding to these reflections are cycle-skipped. As we illustrated earlier, L 2 misfit function being unable to tackle cycle-skipping effects, FWI cannot update the medium correctly to fit these arrivals.

AWI provides a clear improvement over classical L 2 : the horizontal layers are correctly positioned on the left part. The central part at depth (8 < x < 13 km, (z > 2 km) is improved compared to L 2 , with better contrast and more lateral coherency in the layers structure. Moreover, the low-velocity anomaly on the shallow left part of the model is removed. These results are obtained with σ = 0.25 s and ζ = 10 -5 . IE also improves the reconstructed model. The relatively small improvement in cycleskipping robustness introduced by the envelope is enough to mitigate the artifacts on the left part of the model (0 < x < 3 km) and flattens the layers compared to the L 2 result. One drawback is the slight degradation in the reconstruction of the central part at depth (9 < x < 13 km, z > 2 km). Still, such a simple formulation is enough to improve the FWI workflow over the classical L 2 in this case.

The case of NIM misfit is interesting. Here, we can see that it fails to converge, producing an erroneous reconstructed model. This illustrates the limitation of NIM when applied to more realistic cases where the data contains multiple arrivals, multiple phases, and potentially mixed phases. Integrating all these pieces of information into a single observable (the cumulative distribution) does not make it possible to reconstruct the subsurface velocity. As we can see in the data-fit, NIM can also not fit the vast majority of the signal.

The KROT misfit function, as AWI and IE, can prevent the appearance of the left side artifacts observed with the L 2 reconstruction (0 < x < 3 km). As for IE, the relatively small improvement in terms of attraction valley width provided by KROT is sufficient to improve the results significantly. Besides, KROT can account for the lateral coherency of the data, which might also help stabilize the inversion.

Finally, GSOT also produces a significant improvement over the L 2 result, with almost no artifacts on the reconstructed model. We use τ = 0.25 s, similarly to AWI parameterization. This illustrates that GSOT, as AWI, while being able to significantly enlarge the valley of attraction of the misfit function on simple convexity cases, can also be used in a more realistic framework that mixes transmitted and reflected energy with relatively complex multi-arrival data.

Regarding the data-fit, excepted for NIM, all the misfit functions can provide a good data-fit in this case. Starting from this initial model, strong artifacts appear on the L 2 results. We observe long-wavelength low-velocity anomalies on both left and right parts of the models typical of cycle-skipping induced artifacts. The data-fit analysis confirms this observation: only early 2.5 Marmousi case study: toward a more realistic case study arrivals in the near offset are correctly fitted. Diving waves arriving at larger offsets on the left (5 s and -7.5 km offset) and right parts (4 s and 6 km offset) of the gather are cycle-skipped.

Results starting from 1D initial model

Without any surprise, NIM cannot provide a meaningful estimate of the V P model, as it is already the case starting from the S500 initial model.

As the 1D initial model generates large time-shifts, and since both IE and KROT are only marginally improving cycle-skipping robustness, it is not surprising to observe artifacts on the associated reconstructed V P models. IE results present strong artifacts, mainly on the left part of the model (0 < x < 7 km), while the shallow right part (11 < x < 16 km, x <= 2 km) presents an improved reconstruction compared to L 2 . This is confirmed by the data-fit, where all the arrivals on the right part (offset between 1 and 8 km) are correctly predicted, whereas data-fit on the left part (offset between -8 and 0 km) is degradated compared to L 2 data-fit. Conversely, KROT provides a more accurate reconstruction in the left part of the model. The strong low-velocity anomalies observed in the left part (0 < x < 7 km) of the L 2 and IE reconstructions are reduced and appear only in the deep part of the model (z > 2 km). This is consistent with the data-fit, where we see that using KROT, the long offset diving waves for negative offsets of the shot gather are correctly fitted.

AWI manages to provide a clear improvement over classical L 2 , mainly on the center part of the model (5 < x < 14 km), while some artifacts on both sides of the reconstructed model are still present. These parts are more difficult to reconstruct as they are illuminated by waves traveling along the longest paths of the medium, increasing cycle-skipping risk. We set σ = 0.6 s to try to capture as large as possible time-shifts. To get the best results possible, we used ζ = 10 -2 to obtain the helping smoothing effect required to tackle large time-shifts introduced by this initial model. The data-fit obtained with AWI is good, with only some out-of-phase arrivals for late diving waves (around -8 to -6 km offset).

Using GSOT, the reconstructed V P model also presents a clear improvement over classical L 2 . We use τ = 0.6 s in this case. At the first order, most of the artifacts are removed. Still, some artifacts are present close to the edges (0 < x < 1 and 15 < x < 17 km), which is expected from the lack of illumination in these parts. Some other artifacts are visible in the center part of the structure at depth (10 < x < 12 km and 2 < z < 3.5 km). The data-fit appears to be good, with no out of phase arrivals for all offsets.

Error reduction analysis

Besides this qualitative analysis of the results, we can provide quantitative comparisons by analyzing the data error and model error evolution along with iterations. We use here the following relative L 1 model error definition

Err(V P ) = 100 M M i=1 |V P,i -V true P,i | V true P,i (2.38)
where V true P is the true model, M the number of points in the model and i denotes one pixel of the grid used to describe the models at the discrete level.

For the second experiment only (1D initial model), we present the evolution of • the convergence rate (misfit error with respect to the iterations); INVERSION: COMPARISON OF FIVE RECENT APPROACHES

• the L 2 convergence rate (L 2 error with respect to the iterations);

• the model convergence rate (model error with respect to the iterations);

• the model vs. data convergence rate (model error with respect to the misfit error).

For model error, we truncate the model by 1 km on the left and right sides and 625 m at depth to remove the model areas that are not well illuminated.

For alternative misfit function definition, the L 2 -based convergence rate is interesting as moving away from L 2 local minima should be made visible by an increase of the L 2 error with respect to the iterations. Also, the fourth item is interesting, as, ideally, we look for a monotonic decrease of the model error with respect to the misfit error. Besides, to improve the readability, we have excluded from these figures the results corresponding to NIM. The method does not produce reliable results in both cases.

The error reduction analysis is shown in Figure 2.23. First, we observe that KROT and AWI present a relatively slow convergence rate on the cost evolution, while IE and GSOT have a faster convergence rate. L 2 convergence is in between. KROT follows more or less the same as the L 2 misfit function. This is somehow expected, as the valley of attraction of KROT is expected to be similar to the one of the L 2 misfit function. Note, however, that in the early iterations, KROT displays a small increase of the L 2 error, which clearly states that the two misfit functions follow a different minimization path. IE, AWI, and GSOT display another trend: the L 2 error is increased in the first iterations before being strongly decreased in a second stage. The substantial decrease of the L 2 error appears the latest for AWI (after 100 iterations) and the earliest for GSOT (after 30 iterations). GSOT achieves the smallest L 2 misfit, followed by AWI and KROT. The model convergence rate classifies the misfit functions into two groups: one that does not reduce model error compared to the starting point, with L 2 , IE, and KROT; and a second group that decreases the model error with AWI and GSOT. In the second group, only GSOT provides a constant decrease with respect to the iterations, while AWI start to increases the model error until 100 iterations, followed by a decrease. The final reduction of model error obtained with KROT and IE are smaller than the one attained by the L 2 , still, this does not explicitly compared to better interpretable results overall. AWI and GSOT obtain the best reduction of model error. Finally, looking at the model vs. data convergence, only GSOT provides a quasi-monotonic decrease. For all the others, the model error starts by increasing with the reduction of the misfit.

A more realistic inversion

Case study description

Similar to the previous inverse crime inversion, we perform FWI starting from two different initial models. The first one is derived from the true Marmousi model using a lighter Gaussian smoothing, referred to as the S250 model with a correlation length of 250 m (Figure 2.17 b). The second one is the S500 model already used in the inverse crime settings (Figure 2.17 c). We did not re-use the 1D initial model as it proves to be too difficult for any alternative misfit functions to provide convergence in this more realistic case. Not predicting the data to machine precision generates a more challenging benchmark. The wavelet used for FWI is obtained through a source estimation in the initial model based on short offset (100 m) only to decouple the influence of the initial V P model as much as possible. The obtained inverted wavelet is presented in Figure 2.24. We can observe that the inverted wavelet is close to the true wavelet, but some noticeable amplitude and waveform differences are visible. These amplitude effects are induced by the use of a "true" density model for the data generation compared to Gardner's one (Figure 2.18) used for wavelet estimation and to the white noise added to the data.

The data-fits for these two initial models are presented in (a) of Figures 2.26 and 2.27. As expected, the data-fit is better using S250 initial model, while the data-fit generated with S500 initial model displays more out of phase arrivals.

We compare the results obtained using L 2 , AWI, IE, KROT, and GSOT misfit functions. We do not include NIM results here, as we have already shown how the method fails to produce meaningful results in the previous inverse crime settings. A maximum of 500 FWI iterations is performed for both initial models. INVERSION: COMPARISON OF FIVE RECENT APPROACHES

Results starting from S250 initial model

Starting from the S250 initial model (Figure 2.25 left column), the main expected difference with the previous "inverse crime" setup is an inaccurate amplitude prediction (which would be the case if considering field data). Data-fits are presented in Figure 2.27. Classical L 2 can provide an acceptable result. Good reconstruction in the well-illuminated area is achieved, with no visible artifacts in the center part and only a small low-velocity artifact visible at x = 2 km z = 0.8 km and a high-velocity artifact at x = 16 km z = 1 km. The data-fit obtained with L 2 is quite satisfying with most of the arrivals in phase.

This time, IE results are clearly degraded compared to the classical L 2 one. The reconstructed V P model is tainted with high wavenumber oscillation and strong artifacts. This is an indication that the IE approach is sensitive to a correct amplitude prediction. This validates the interest of a more realistic framework, making us able to detect this kind of limitation. The data-fit presents many out-of-phase arrivals, coherent with the small artifacts present everywhere in the reconstructed V P .

Again, AWI improves over the L 2 results. We use a relatively small σ = 0.2 s here as the maximum time-shifts expected are relatively small with this good initial model. We used ζ = 10 -2 as noise requires a relatively large amount of damping, moreover as illustrated before, a larger damping value helps when facing challenging FWI setups. The deep center part is improved with a more coherent deep-layer structure. The left (x = 2 km and z = 0.8 km) and right (x = 16 km z = 1 km) side artifacts present in L 2 results are also partially mitigated. Surprisingly, the data-fit obtained with AWI is poor for large offset arrivals (from -8 km to -3 km and 3 to 8 km). This degradation of the data-fit is slightly counter-intuitive and does not correlate with the improvement of the reconstructed V P model observed.

Finally, KROT and GSOT reconstructed models both present similar improvement compared to the L 2 one. We can observe an increase in terms of high wavenumber content. Interestingly, the deep center part (9 < x < 13 km, z > 2 km), which is the main target of interest of the Marmousi model (an anticlinal structure) is more resolved using KROT and GSOT compared to L 2 . For GSOT, we use τ = 0.2 s in this case. The data-fit obtained with both methods is good, with almost all arrivals in phase. Only some first arrivals between -4 to -2 km offset are still not well explained. The GSOT data-fit appears to be slightly better INVERSION: COMPARISON OF FIVE RECENT APPROACHES than the KROT one.

Results starting from S500 initial model

Starting from the S500 model, reconstructed V P results are presented in Figure 2.25 right column, while data-fit are presented in Figure 2.27.

Here, the classical L 2 fails to reconstruct a meaningful V P model. Many artifacts are present on the model that may come in part from cycle-skipping. This would prevent any interpretation of the reconstructed model. The data-fit present out-of-phase arrivals, even if the majority would appear to be in-phase. This again illustrates potential convergence toward a local minimum that makes possible to fit the data with non-meaningful V P updates.

With no surprise, IE fails to reconstruct a meaningful V P estimate. The reconstructed V P model suffers from many artifacts. The data-fit is clearly degraded compared to L 2 , which is 2.5 Marmousi case study: toward a more realistic case study AWI reconstructed model produces here an improvement over L 2 or IE, with the central part and right part of the Marmousi model more or less retrieved. However, significant artifacts are present in the left part of the model (1 < x < 6 km) associated with an erroneous reconstruction of the central part at depth (9 < x < 13 km, z > 2 km). Here we increase σ to 0.4 s, and keep ζ = 10 -2 . The data-fit is degraded with out of phase arrivals for offsets between -8 to -3 km as well as between 2 to 8 km.

KROT produces satisfactory results here. This is interesting as KROT only marginally improves cycle-skipping robustness. Here, it manages to perform well in this complexified case. This is a good indication that the difficulties induced in this more realistic inversion are not only cycle-skipping but also amplitude mismatch (due to density) and noise. As KROT introduces lateral coherency and has a regularizing effect on noise, it is not surprising to observe a better behavior in this case. The data-fit obtained with KROT is good with almost all arrivals in phase, except for some transmitted waves from -3 to -1 km offset and some long offset arrivals INVERSION: COMPARISON OF FIVE RECENT APPROACHES around -8 to -7 km.

Finally, GSOT provides a good reconstructed V P model. The central part (9 < x < 13 km, z > 2 km) is well reconstructed. The layers show more lateral coherency compared with KROT. Furthermore, left side artifacts are reduced compared to KROT. Again and similarly to AWI, τ is increased to 0.4 s to account for the larger time-shifts introduced by the degraded initial model. The data-fit is also good, with improvement over the KROT for the long offset arrivals around -8 to -7 km.

Error reduction analysis

A similar analysis for the different misfit functions is presented for this inverse crime inversion of Marmousi. The model error is calculated in a similar zone as in the previous experiment.

Starting from the S500 initial model (Figure 2.28), we observe that KROT and AWI present again a relatively slow convergence rate (AWI being the slower), while L 2 , IE and GSOT have a faster convergence rate. The L 2 data-error is again interesting, with GSOT and KROT performing the most substantial reduction of L 2 data error (with an initial jump to pass a L 2 local minimum for GSOT at the first iteration). While IE increases the cost drastically for the first two iterations, it then fails to reduce the data error. We can note that KROT is not following L 2 misfit function behavior anymore compared to the inverse crime Marmousi case. Regarding AWI, we can observe that it starts to increase the L 2 data error until 20 iterations, then rapidly reduce for 10 iterations, to finish with a constant increase afterward. This time, the model error displays a strong increase for L 2 and IE misfit functions, which is coherent with the artifacts present in the reconstructed V P models. AWI is also increasing the model error as it is also affected by artifacts, but less drastically than L 2 and IE, which is visible on the reconstructed V P model. KROT and GSOT manage to decrease the model error continuously. Looking at the model vs. data convergence, only KROT and GSOT present monotonic behavior, while L 2 , IE, and AWI are increasing the model error.

Computational cost

The computational overhead induced by the alternative misfit function selected in this review varies from +2 to +30% compared to L 2 misfit. These values are coherent with the values documented in the literature. The key feature here is that even a +30% computational overhead is not a blocking feature and is affordable with modern computing facilities. For us, the key feature is the "physical" performance of the misfit function that translates into an improvement of FWI robustness.

Discussion

Discussion

Among the five misfit functions compared here, namely NIM, IE, AWI, KROT, GSOT, three of them show a significant improvement in convexity with respect to a time-shift: NIM, AWI, and GSOT. However, when applied to a realistic case (Marmousi), NIM fails to produce a meaningful V P estimate. Conversely, for AWI, while difficulties are identified on schematic examples, including multiple arrivals (a situation known to be problematic for correlation and deconvolution approaches), satisfactory results are obtained when applied to the Marmousi case, both within and without the inverse crime settings. GSOT also appears as an interesting strategy, providing satisfactory results in all the tests performed here.

Interestingly, while IE and KROT show less robustness to strong cycle-skipping, the small increase in the valley of attraction they provide is sufficient to enhance the velocity reconstruction in the Marmousi test in inverse crime settings. However, IE fails when it comes to non-inverse crime settings, that is when noise corrupts the data, and amplitude prediction cannot be guaranteed anymore. On the contrary, KROT reveals relatively robust to these settings, probably benefiting from its ability to account for the lateral continuity of events in shot-gather representation and for the robustness of optimal transport based distances with respect to the presence of noise [START_REF] Engquist | Optimal transport for seismic full waveform inversion[END_REF].

From the experiment performed in this article, KROT, AWI, and GSOT appear as an interesting alternative to the least-squares distance from the perspective of field data application. In cases where no strong cycle-skipping is expected, KROT should perform well, and this is supported by several field data applications already performed on exploration data [START_REF] Messud | Multidimensional optimal transport for 3d FWI: Demonstration on field data[END_REF][START_REF] Sedova | Acoustic land full waveform inversion on a broadband land dataset: the impact of optimal transport[END_REF][START_REF] Carotti | Optimal transport full waveform inversion -applications[END_REF]. The computational cost of KROT is relatively higher than that of AWI and GSOT; however, its ability to account for the lateral coherency of the data in shot-gather panels makes it an appealing strategy. For 3D data cubes, cutting it into 2D slices and summing over the slices is a good compromise. To deal with larger kinematics inaccuracy, AWI and GSOT should be preferred options. AWI has already been successfully applied to field data [START_REF] Warner | Robust adaptive waveform inversion[END_REF][START_REF] Ravaut | Imaging beneath a gas cloud in the north sea without conventional tomography[END_REF][START_REF] Debens | Full-bandwidth adaptive waveform inversion at the reservoir[END_REF][START_REF] Roth | Improving seismic image with high resolution velocity model from awi starting with 1d initial modelcase study barents sea[END_REF][START_REF] Guasch | Adaptive waveform inversion: Practice[END_REF][START_REF] Warner | Extending the velocity resolution of waveform inversion below the diving waves using awi[END_REF]. We, however, show here that it could suffer from some limitations in the case of complex data containing multiple arrivals. GSOT has been mostly applied to synthetic data by now (He et al., 2019a;Provenzano et al., 2020a). Nevertheless, field data applications are ongoing [START_REF] Pladys | Graph space optimal transport based FWI: 3D OBC valhall case study[END_REF][START_REF] Górszczyk | Relaxing the initial model constraint for crustal-scale full-waveform inversion with graph space optimal transport misfit function[END_REF].

Conclusion

This article is dedicated to comparing misfit function reformulation for FWI, which aims at mitigating cycle-skipping. The first result drawn is that the link between cycle-skipping and the non-convexity with respect to time-shifts of the least-squares distance is evident from the different tests we provide. However, when no such cycle-skipping occurs (sufficiently accurate initial model), least-squares FWI performs well, even for complex data including multiple phases, mixed phases, noise, and when amplitude prediction cannot be performed accurately (as is the case for field data). Therefore, efficient reformulation of the FWI misfit function should not rely only on a better convexity to time-shifts to replace the least-squares norm advantageously but should also exhibit robustness with respect to these settings, which are always met on field data applications.

Chapter 3

Robust FWI with graph space optimal transport: application to 3D OBC Valhall data

Introduction abstract

Improving full-waveform inversion to make it more robust to cycle-skipping has been the subject of a large number of studies. From the several families of approaches developed, one of the most documented consists in modifying the least-squares distance defining the discrepancy between observed and calculated data. From all the propositions made to improve and replace the least-squares distance, only a few of them have been applied to field data. One of the methods proposed recently, the graph space optimal transport distance, presents appealing properties for field data applications. This study compares it with the least-squares distance in an analysis performed on the three-dimensional ocean bottom cable data from the Valhall field. This data has already been at the heart of several full-waveform inversion studies, making it a good candidate to evaluate the properties of this new misfit function. We first perform this comparison starting the inversion from the reflection traveltime tomography model used in previous studies. We then perform a second comparison from a crude, linearly varying in-depth one-dimensional velocity model. Starting from this model, least-squares-based full-waveform inversion fails to provide a meaningful estimate of the pressure-wave velocity model due to cycle skipping. We illustrate how the graph-space optimal transport-based full-waveform inversion mitigates this issue. A meaningful estimate of the pressure-wave velocity model is obtained in the zone sampled by both diving and reflected waves, down to almost two kilometers depth. To our knowledge, this is the first application of a graph space optimal transport-based fullwaveform inversion to three-dimensional field data.

Introduction

Full waveform inversion (FWI) is a seismic imaging method aiming at reconstructing highresolution models (up to half the shortest wavelength) of the mechanical properties of the subsurface [START_REF] Devaney | Geophysical diffraction tomography[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain, part II: Fault delineation in sediments using crosshole data[END_REF][START_REF] Plessix | Full waveform inversion of a deep water ocean bottom seismometer dataset[END_REF][START_REF] Raknes | Three-dimensional elastic full waveform inversion using seismic data from the sleipner area[END_REF][START_REF] Górszczyk | Toward a robust workflow for deep crustal imaging by FWI of OBS data: The eastern nankai trough revisited[END_REF]. The method is an iterative process based on minimizing a misfit function between observed and calculated data over a space of model parameters describing the subsurface. The improvement of resolution that FWI provides over standard tomography methods makes it possible to significantly improve depth-migration images or directly produce interpretable images of the subsurface physical properties. This method is used at multiple scales; from global and regional scales [START_REF] Fichtner | Full waveform tomography for radially anisotropic structure: New insights into present and past states of the Australasian upper mantle[END_REF][START_REF] Tape | Seismic tomography of the southern California crust based on spectral-element and adjoint methods[END_REF][START_REF] Bozdag | Global adjoint tomography: first-generation model[END_REF] to seismic exploration targets for oil & gas industry [START_REF] Plessix | Full waveform inversion of a deep water ocean bottom seismometer dataset[END_REF][START_REF] Stopin | Multiparameter waveform inversion of a large wide-azimuth low-frequency land data set in Oman[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 approximation[END_REF] and even near-surface scale [START_REF] Bretaudeau | 2D elastic full waveform imaging of the near surface: Application to synthetic and a physical modelling data sets[END_REF][START_REF] Groos | The role of attenuation in 2D full-waveform inversion of shallow-seismic body and Rayleigh waves[END_REF][START_REF] Schäfer | 2D full waveform inversion of recorded shallow seismic Rayleigh waves on a significantly 2D structure[END_REF][START_REF] Irnaka | Towards 3d 9c elastic full waveform inversion of shallow seismic wavefields -case study ettlingen line[END_REF]. A thorough review of FWI and its applications is given in [START_REF] Virieux | An introduction to Full Waveform Inversion[END_REF].

As powerful as this method is, it suffers from a significant shortcoming in its classical formulation: the non-convexity with respect to time-shifts of the least-squares (L 2 ) misfit function used to calculate the distance between observed and synthetic data. This non-convexity of the misfit function is an issue as the iterative process used in FWI is based on local-optimization algorithms. This leads to the so-called cycle-skipping issue. This limitation of FWI in its classical formulation has been documented since FWI has been introduced [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: numerical results[END_REF], and it has been of great interest to overcome it.

A large amount of studies has been published, proposing different approaches. One of OBC VALHALL DATA the historical approaches to overcome this limitation in practical cases is to rely on a data hierarchy workflow [START_REF] Bunks | Multiscale seismic waveform inversion[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] Shipp | Two-dimensional full wavefield inversion of wide-aperture marine seismic streamer data[END_REF][START_REF] Wang | Reflection seismic waveform tomography[END_REF][START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. This approach consists in interpreting first the lowest frequency available, generally 2 -4 Hz for seismic exploration targets, then progressively introducing higher frequency data following a multi-scale approach [START_REF] Sirgue | Efficient waveform inversion and imaging : a strategy for selecting temporal frequencies[END_REF]. A second level of data-hierarchy can be defined by modifying the temporal and/or offset selection of the data used during inversion. The idea is to reduce the number of propagated wavelengths that are interpreted simultaneously. Current industrial applications generally rely on these two levels of data hierarchy, combined with a robust starting model, obtained, for instance, through reflection traveltime tomography or stereotomography [START_REF] Lambaré | Stereotomography[END_REF].

Nonetheless, the conditions to apply this workflow are not always satisfied. For instance, low-frequency data around 2 -4 Hz are not always available or of sufficient quality. Moreover, obtaining low frequency increases the cost of the acquisition and can also compromise the quality of the high frequency needed to obtain high-resolution model reconstructions. Accurate initial model building can also be a time-consuming and challenging task requiring strong human expertise as it generally requires accurate traveltime and/or reflected event picking. Besides, prior information coming from geology and well logs are often needed. These constraints make FWI less robust and reduce its potential range in terms of applications.

Several methods have been introduced to improve robustness to cycle-skipping. The first group can be named as "extension strategies" and relies on introducing supplementary degrees of freedom to the FWI problem [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF][START_REF] Symes | Algorithmic aspects of extended waveform inversion[END_REF]Huang et al., 2018a;[START_REF] Van Leeuwen | Mitigating local minima in full-waveform inversion by expanding the search space[END_REF][START_REF] Wang | Full-waveform inversion with the reconstructed wavefield method[END_REF][START_REF] Aghamiry | Accurate and efficient wavefield reconstruction in the time domain[END_REF], which can be used to artificially match the data at early iterations of the FWI process, avoiding cycle-skipping.

The second group consists in reformulating the FWI problem using an alternative measure of the distance between the observed and calculated data, namely a different misfit function. Numerous approaches have been proposed, such as cross-correlation [START_REF] Luo | Wave-equation traveltime inversion[END_REF][START_REF] Van Leeuwen | A correlation-based misfit criterion for waveequation traveltime tomography[END_REF] and deconvolution based misfit function [START_REF] Luo | A deconvolution-based objective function for wave-equation inversion[END_REF][START_REF] Warner | Adaptive waveform inversion: Theory[END_REF], or by modifying the signal itself, making the L 2 norm between this new observable more convex with, for instance, instantaneous envelope [START_REF] Fichtner | Theoretical background for continental-and global-scale full-waveform inversion in the time-frequency domain[END_REF][START_REF] Bozdag | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF]. We want to keep in mind that replacing the L 2 norm is not an easy task, as, despite its simplicity, the L 2 misfit presents excellent and interesting properties. First, it is robust to Gaussian noise. Second, it presents an excellent resolution power, translating into high-resolution reconstruction that FWI is well known for. Third, it is straightforward to implement and the computational cost of the misfit function evaluation is negligible compared to most of the proposed alternative misfit functions. These advantages have made the L 2 misfit the "state of the art" for FWI at exploration scales and could explain why L 2 is still widely used even if many alternative misfit functions have been proposed to mitigate the cycle-skipping issue.

Indeed, there is a discrepancy between the many propositions for alternative misfit functions compared to the number of actual field data applications. We think that this discrepancy could be explained by the -often not deeply discussed -intrinsic limitations of these alternative formulations. For instance, cross-correlation-based misfit functions might have difficulties handling complex data, including missing arrivals or multiples. Deconvolution-based strategies require a penalization/weighting function, which can be difficult to set. Such settings are often case-dependent, making FWI less of an automated process. Instantaneous envelope intrinsically modifies the signal shape and discards information coming from the phase (which 3.1 Introduction is essential to interpret the polarity of reflected events in the data correctly).

A new class of misfit functions based on optimal transport (OT) has been introduced recently. The motivation is to benefit from the convexity of the optimal transport distance with respect to translation and dilation, which provides a misfit function convex with respect to time and amplitude shifts, a good proxy towards convexity with respect to velocities perturbations [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF][START_REF] Métivier | Optimal transport for mitigating cycle skipping in full waveform inversion: a graph space transform approach[END_REF]. Another important motivation to use OT as a misfit function is the ability to take into account the coherency of the seismic signal in an adequate space, be it a common shot or receivers gather. However, OT can only be applied to positive quantities and cannot be directly applied to seismic traces. To circumvent this difficulty, three main strategies have been developed.

The first one proposes to bring back the problem to the comparison of positive quantities by modifying the signal before solving the OT problem [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF][START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport norm[END_REF][START_REF] Yang | Application of optimal transport and the quadratic Wasserstein metric to full-waveform inversion[END_REF][START_REF] Yang | Analysis of optimal transport and related misfit functions in full-waveform inversion[END_REF]. A nonlinear transform is applied to the data in a trace-by-trace framework to transform each of them as probability measures. However, modifying the signal and altering the polarity information might be detrimental to a stable and satisfactory reconstruction of the subsurface mechanical properties.

The second one relies on the dual formulation of a specific instance of optimal transport distance, namely the 1-Wasserstein distance (Métivier et al., 2016a,b,c). This formulation can be naturally extended to the comparison of signed data. It benefits from its ability to be applied directly to 2D and 3D data, taking into account the coherency of the seismogram in the receiver and/or sources direction. However, even if the attraction valley to the global minimum is enlarged compared with the least-squares approach, the application of the 1-Wasserstein distance to signed data loses the convexity with respect to time-shift, which was the original motivation to use OT in the framework of FWI (see [START_REF] Métivier | Optimal transport for mitigating cycle skipping in full waveform inversion: a graph space transform approach[END_REF] for a review on different OT formulations). This strategy has been successfully applied to several field datasets [START_REF] Poncet | Fwi with optimal transport: a 3D implementation and an application on a field dataset[END_REF][START_REF] Messud | Multidimensional optimal transport for 3d FWI: Demonstration on field data[END_REF][START_REF] Sedova | Acoustic land full waveform inversion on a broadband land dataset: the impact of optimal transport[END_REF][START_REF] Carotti | Optimal transport full waveform inversion -applications[END_REF][START_REF] Hermant | Imaging complex fault structures on-shore oman using optimal transport full waveform inversion[END_REF].

Finally, the third one considers each discrete seismic trace as point clouds and computes the optimal transport distance between point clouds associated with synthetic and observed traces. This method is called the graph space optimal transport (GSOT) and presents the main characteristic of preserving the convexity with respect to time shifts [START_REF] Métivier | Optimal transport for mitigating cycle skipping in full waveform inversion: a graph space transform approach[END_REF][START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF]. GSOT has already been successfully applied to 3D synthetic and field data (He et al., 2019a;[START_REF] Pladys | Assessment of optimal transport based FWI: 3d OBC valhall case study[END_REF][START_REF] Li | Joint FWI for imaging deep structures: A graphspace OT approach[END_REF][START_REF] Górszczyk | Mitigating the nonlinearity of the crustal scale full waveform inversion through the graph space optimal transport misfit function[END_REF] and is a promising candidate to tackle cycle-skipping on field datasets.

Therefore, this study focuses on applying the GSOT strategy against the classical L 2 misfit on a 3D OBC dataset from the North Sea, the Valhall field data. This dataset has been one of the first used to make proof of concept of the resolution power that FWI can bring on field data as shown in [START_REF] Sirgue | Full waveform inversion: the next leap forward in imaging at Valhall[END_REF]. Since then, this dataset has been used several times for FWI application [START_REF] Prieux | On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study[END_REF][START_REF] Gholami | Which parametrization is suitable for acoustic VTI full waveform inversion? -Part 2: application to Valhall[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Operto | On the role of density and attenuation in 3D multiparameter visco-acoustic VTI frequency-domain FWI: an OBC case study from the North Sea[END_REF][START_REF] Kamath | Multiparameter full-waveform inversion of 3D OBC data from the valhall field[END_REF]. This dataset can been seen as a "calibrated reference" for testing FWI formulations, such as frequency-domain FWI in [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 approximation[END_REF] or in time-domain with attenuation in [START_REF] Kamath | Multiparameter full-waveform inversion of 3D OBC data from the valhall field[END_REF]. Here, we are using this dataset to compare the GSOT misfit function to the conventional L 2 norm through a time-domain 3D visco-acoustic VTI FWI.

To make this comparison, we first consider a canonical case where the initial model is the OBC VALHALL DATA same as the one used in the aforementioned studies, which is derived from reflection tomography. This initial model ensures FWI converges toward a plausible estimation of the subsurface using conventional L 2 misfit function.

Then, we introduce a very "crude" initial model in which calculated data are shifted by more than one cycle compared to observed data, which is the typical case scenario for cycle-skipping. We show how L 2 -based FWI fails from this "crude" starting model, whereas GSOT-based FWI manages to correctly interpret data to provide a plausible reconstruction of the subsurface (down to almost 2 km depth, which corresponds to the zone of the data sampled by both diving and reflected waves). This constitutes one of the first applications of the graph-space optimal transport misfit function to 3D field data at the exploration scale.

The paper is organized as follows. In the next section, we present the modeling and inversion algorithm used for our FWI application. Then we detail the global methodology of GSOT for FWI. This is followed by a presentation of the Valhall field application, from the geological situation to the initial model and dataset presentation. We detail our FWI workflow and analyze the results in two cases: from the reflection traveltime tomography initial model and then from a 1D initial model. Results are then discussed, followed by conclusion and perspectives, which are given in a final Section.

Methodology

Modeling

This study is performed in the frame of 3D time-domain FWI. We rely on the anisotropic viscoacoustic time-domain modeling and inversion algorithm developed by Yang et al. (2018a), based on the following partial differential equations:

                   ρ∂tvx = ∂xg ρ∂tvy = ∂yg ρ∂tvz = ∂zq ∂tg = c11(∂xvx + ∂yvy) + c13∂zvz -L =1 Y [c11ξ g + c13ξ q ] ∂tq = c13(∂xvx + ∂yvy) + c33∂zvz -L =1 Y [c13ξ g + c33ξ q ] ∂tξ g
= -ω ξ g + ω (∂xvx + ∂yvy), = 1, 2, . . . , L ∂tξ q = -ω ξ q + ω ∂zvz, = 1, 2, . . . , L .

(3.1)

In system 3.1, c 11 , c 13 , and c 33 are the stiffness tensor coefficients, ρ is the density, v x , v y , v z are the horizontal and vertical displacement velocities respectively, while g and q are related to the normal stress components σ xx , σ yy and σ zz through

g = σ xx = σ yy q = σ zz . (3.2)
This simplification is due to the VTI approximation. Similarly, the memory variables ξ q and ξ g are related to the memory variables ξ xx , ξ yy and ξ zz associated with the normal stress components through

ξ g = ξ xx + ξ yy ξ q = ξ zz . (3.3)

Methodology

These memory variables are used to model the viscosity of the medium following the generalized Maxwell body theory. Each represents one relaxation mechanism. We use three relaxation mechanisms to approximate a constant attenuation within the considered frequency band (L=3). The variables Y are therefore calibrated depending on the target quality factor representing the attenuation in the considered media. This calibration is done through the solution of a least-squares problem. The details of this calibration can be found in Yang et al. (2016a) for instance.

In the VTI approximation, the stiffness tensor coefficients are related to the vertical P-wave velocity, the density, and the Thomsen anisotropy parameters and δ through

c 11 = ρV 2 P (1 + 2 ) c 33 = ρV 2 P c 13 = ρV 2 P (1 + 2δ) .
(3.4)

The discretization of this system of partial differential equations is performed using a fourthorder in space and second-order in time staggered grid finite-difference method [START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite difference method[END_REF][START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF].

A flat free surface condition is applied on top of the model to represent the water/air interface. Sponge layers [START_REF] Cerjan | A nonreflecting boundary condition for discrete acoustic and elastic wave equations[END_REF] are applied on the other faces of the model to mimic a medium of infinite extensions in these directions. This numerical method is preferred over the use of perfectly matched layers (PML, [START_REF] Bérenger | A perfectly matched layer for absorption of electromagnetic waves[END_REF] both because of the anisotropy, which might cause PML instabilities, as well as the presence of attenuation and relaxation mechanisms, which complexifies the PML implementation.

Finally, windowed sync interpolation is used to simulate source and receivers off-grid points accurately [START_REF] Hicks | Arbitrary source and receiver positioning in finite-difference schemes using Kaiser windowed sinc functions[END_REF].

Inversion

General formulation

FWI is an iterative process which relies on the minimization of a misfit function. Classically, it is the L 2 misfit function defined as follows

f L 2 [m] = s r h d cal [m](x r , t; x s ), d obs (x r , t; x s ) , (3.5)
where

h(d 1 (t), d 2 (t)) = 1 2 T t=0 |d 1 (t) -d 2 (t)| 2 dt . (3.6)
The observed and synthetic trace calculated in model m, associated with source x s and receiver x r , are denoted by d obs (x r , t; x s ) and d cal [m](x r , t; x s ) respectively.

To solve this minimization problem, we rely on a local optimization scheme. We use a preconditioned quasi-Newton l-BFGS algorithm [START_REF] Nocedal | Updating Quasi-Newton Matrices With Limited Storage[END_REF], implemented in the SEISCOPE optimization toolbox (Métivier and Brossier, 2016).

We compute the gradient following the adjoint state method [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF]. This method makes it possible to easily replace the L 2 norm as a change of the misfit function only translates ROBUST FWI WITH GRAPH SPACE OPTIMAL TRANSPORT: APPLICATION TO 3D OBC VALHALL DATA to a modification of its associated adjoint-source to obtain the gradient. See Yang et al. (2018a) for a review. Since the method relies on reverse time propagation of the wavefield, which is numerically unstable with attenuation, our code relies on the checkpoint-assisted reverse forward simulation (CARFS) strategy proposed by [START_REF] Yang | Wavefield reconstruction in attenuating media: A checkpointing-assisted reverse-forward simulation method[END_REF] to provide a stable and yet efficient implementation for large scale problems with attenuation.

The preconditioner we use is the wavefield preconditioner presented in [START_REF] Kamath | Multiparameter full-waveform inversion of 3D OBC data from the valhall field[END_REF] An anisotropic non-stationary Gaussian smoothing is also applied to the gradient. Inversions are performed for P-wave velocities (V P ) only.

Optimal transport for FWI

Optimal transport (OT) distances are derived from the OT theory introduced by the French mathematician Gaspard Monge more than two centuries ago [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF]. More precisely, they rely on the OT relaxation proposed by [START_REF] Kantorovich | On the transfer of masses[END_REF]. The distances, also called Wasserstein distances, have an intrinsic property of particular interest for the definition of inverse problems: they are convex with respect to translation and dilation of the compared quantities. This convexity with respect to translation has been, in particular, the motivation to introduce it in the framework of FWI to obtain a distance measurement convex with respect to time-shifts [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF]. However, OT distances are defined for comparing probability distributions, which are by definition positive and normalized. So this new distance cannot be directly applied to seismic data, which is oscillatory (a generalization of OT to signed distribution is still an open question from a mathematical point of view, see [START_REF] Ambrosio | Gradient flow of the Chapman Rubinstein Schatzman model for signed vortices[END_REF] and [START_REF] Mainini | A description of transport cost for signed measures[END_REF] for instance.

This study focuses on a recent proposition made to apply OT to seismic data: the graph space optimal transport (GSOT) strategy, proposed in [START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF].

This formulation of OT distance should preserve the signal unmodified while also preserving the convexity to shifted patterns. It relies on the idea of comparing the discrete graph of the data rather than the data itself: each 1D trace in time becomes a point cloud of Dirac delta functions (of amplitude 1) in a 2D space made of the time dimension, and a new amplitude dimension. This transformation does not affect the signal shape but makes it possible to deal with positive mass (the Dirac delta functions). Thus the OT distance can be applied directly.

The corresponding misfit function is formulated as

f GSOT [m] = s r h d cal [m](x r , t; x s ), d obs (x r , t; x s ) , (3.7)
where this time

h(d 1 , d 2 ) = min σ∈S(Nt) Nt i=1 c iσ(i) (d 1 , d 2 ) . (3.8)
S(N t ) denotes the ensemble of permutations of {1, . . . , N t }, and c ij the L 2 distance between the discrete points of the graph (t i , d 1 (t i )) and (t j , d 2 (t j )):

c ij (d 1 , d 2 ) = |t i -t j | 2 + η 2 |d 1 (t i ) -d 2 (t j )| 2 .
(3.9)

The function h corresponds to the 2-Wasserstein distance between the discrete graph of the "calculated" trace d 1 (t) and the "observed" trace d 2 (t).

Field data presentation and full-waveform inversion workflow

The scaling parameter η controls the convexity of the misfit function f GSOT with respect to time shifts. In practice, we define it as

η = τ A , (3.10)
where τ is a user-defined parameter corresponding to the maximum expected time shift between observed and calculated data in the initial model, and A is the maximum peak amplitude difference between observed and calculated data. This ensures the convexity of the GSOT distance for time up to approximately τ .

The final cost function we use for the purpose of FWI application with N s shots containing N r receivers is defined as: (3.11) where w s,r is a trace-by-trace weighting factor, typically used to restore the AVO trend in the data. This trend is removed from the trace-by-trace GSOT approach, as the amplitude of each trace is treated separately through the normalization factor η.

min m f GSOT [m] = Ns s=1 Nr r=1 w s,r h d s,r cal [m], d s,r obs ,
The adjoint source of the misfit function f GSOT [m] is computed from ∂h ∂ cal using the adjointstate strategy [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF]. The following result is proved in [START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF]. Denoting σ * the minimizer in eq. 3.7, we have

∂h ∂ cal = 2 d cal -d σ * obs , (3.12)
where d σ * obs (t i ) = d obs (t σ * (i) ) .

(3.13)

The GSOT approach can thus be viewed as a generalization of the L 2 distance: The adjoint source is equal to the difference between calculated and observed data at time samples connected by the optimal assignment σ * . The solution of the problem eq. 3.7 provides the information to compute both the misfit function and the adjoint source.

To solve eq. 3.7 efficiently, we use the auction algorithm [START_REF] Bertsekas | The auction algorithm for the transportation problem[END_REF], dedicated to the solution of linear sum assignment problems such as eq. 3.7. Despite a relatively high computational complexity in O(N 3 t ), it is quite efficient for small instances of such problems. Resampling the data close to the Nyquist frequency (at the exploration scale, under acoustic approximation, traces are often around few hundreds of time steps after resampling) yields such small-scale problems making the GSOT feasible for realistic scale FWI applications, as seen in the application presented here. nearly constant water depth of 70 m. Valhall reservoir lies along the Lindesnes Ridge, which trends NNW [START_REF] Munns | The Valhall field: a geological overview[END_REF][START_REF] Leonard | Valhall Field in Geology of Norwegian Oil and Gas Fields[END_REF]. The field has been discovered in 1975 and is used since then for oil production. It is characterized as an anticlinal in chalk in the Upper Cretaceous Hod and Tor formations, which form the reservoir at a depth of approximately 2400 m. Trapped gas in Tertiary shale is present above the reservoir [START_REF] Sirgue | Full waveform inversion: the next leap forward in imaging at Valhall[END_REF][START_REF] Prieux | On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[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 approximation[END_REF]. The Tertiary overburden is relatively simple and free of complex structure [START_REF] Hall | Fractured reservoir characterization using P-wave AVOA analysis of 3D OBC data[END_REF]. Chalk compaction resulting from pressure depletion and water weakening has led to seabed subsidence (Field: VALHALL -Norwegianpetroleum.no https://www.norskpetroleum.no/en/facts/field/valhall).

Field data presentation and full-waveform inversion workflow

Dataset presentation

The seismic data are 4-components acquired by ocean-bottom cables (OBC), with wide aperture/azimuth acquisition. The covered zone is a surface of 145 km 2 . Twelve receiver cables are deployed on the seabed, containing 2048 receivers with an inline spacing of 50 m and a cable spacing of 300 m. A total of 50824 shots are available, located 5 m below the sea surface. The layout of this 3D acquisition is presented in Figure 3.2. In this study, we use only the pressure component of one acquisition performed in 2011 as part of the Valhall Life of Field Seismic (LoFS) project [START_REF] Barkved | Life of Field Seismic Implementation -Another "first at valhal[END_REF]. From the raw data provided in SEGY format, only a simple pre-processing is applied. As our FWI code relies on source parallelization, source-receiver reciprocity is applied to process the hydrophone as explosive sources and the shots as hydrophones, hence sensibly reducing the impact on computer resources. The data is then de-spiked before a quality control over the complete dataset: the energy (RMS) of each gather is calculated to manually remove faulty gathers (the one with a large variation of RMS amplitude). The last step is to create frequency bands for the inversion using a minimum-phase band-pass filter in two distinct bands: 2.5-5 Hz and 2.5 -7 Hz (referred to as band 1 and band 2 in the following). This goes in hand with time-decimation; from ∆t = 4 ms to ∆t = 8 ms on the first band, and from ∆t = 4 ms to ∆t = 5 ms on the second band (Fig. 3.3).

Initial models

An initial V P model was made accessible to us thanks to AkerBP. It has been obtained through reflection traveltime tomography. It is referred to as TOMO initial and is presented in Figure 3.4. This initial model has been used in several publications using this Valhall dataset and has proven its capacity to give satisfactory results with L 2 -based FWI as it predicts the arrival within half a period of the considered starting frequencies (2.5 -5 Hz), avoiding cycle-skipping issue [START_REF] Prieux | On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study[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 approximation[END_REF][START_REF] Operto | On the role of density and attenuation in 3D multiparameter visco-acoustic VTI frequency-domain FWI: an OBC case study from the North Sea[END_REF]. The associated density model is derived from V P TOMO using Gardner's law [START_REF] Gardner | Formation velocity and density-the diagnostic basics for stratigraphic traps[END_REF], defined as 3.3 Field data presentation and full-waveform inversion workflow ρ = 309.6 * V 0.25 P . This relation is a fair average for brine-saturated rock (excluding evaporites), which is coherent with the expected geology of the Valhall field. AkerBP also provided us the anisotropy model and δ , and while their structure is not complicated, they are of significant influence in the modeling. The anisotropic parameter η define by [START_REF] Alkhalifah | Acoustic approximations for processing in transversely isotropic media[END_REF] as

η = -δ 1 -2δ (3.14)
is shown in Figure 3.6. We can see that maximum anisotropy reach values around 15% near the reservoir. Finally, to introduce attenuation, a simple two-layer Q P model, with 1000 in the water column and 200 in the sediments, is used [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 approximation[END_REF]. While it is a simple model, it has been proven to be of great importance to explain the data. Using a more complex Q P model (for example derived from V P ) only results in marginal improvements. Moreover, it introduces complexity and uncertainties that we prefer not to deal with [START_REF] Kamath | Multiparameter full-waveform inversion of 3D OBC data from the valhall field[END_REF].

As the TOMO initial model is good enough to match the data within half a period of the 2.5 -5 Hz frequency-band used to start FWI (Fig. 3.7), it does not represent any challenge regarding the cycle-skipping issue. This is why we introduce a new "crude" V P starting model called 1D initial presented in Figure 3.5. It is a purely 1D vertical starting model, based on a linearly increasing profile with one main interface around 2400 m depth. This 1D starting model generates strong cycle-skipping, as made visible in Figure 3.8. As traditional L 2 -based FWI cannot tackle this cycle-skipping in the data, it should be a good candidate to benchmark the capability of GSOT FWI. Associated with this V P 1D model, the density model is derived using the preceding Gardner's law. The anisotropy model and δ , as the attenuation model Q P are kept similar to those used in the TOMO setup.

Full-waveform inversion workflow

To obtain the final FWI reconstructed V P model, we rely on several elements that compose our complete workflow. Each part of it plays a critical role in obtaining the best results possible, from wavelet estimation, data selection, to post-processing the FWI model. The generic workflow that we use in this article is presented in Figure 3.9. We detail the different parts composing the workflow in the following sections.

Source wavelet estimation

The first step before running FWI is to obtain a proper source wavelet that is used during the inversion to generate synthetic data. This is a crucial step as an incorrect wavelet estimation could induce artifacts into the reconstructed model.

The wavelet is estimated by solving a linear deconvolution problem in the frequency domain, following the methodology described in [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part I: theory and verification in a physical scale model[END_REF]. The wavelet inversion is performed on a single random subset of 240 shot-gathers. The sample of 240 wavelets (one per shot-gather) is then averaged to produce a single wavelet, assuming that all hydrophones have the same coupling response.

To minimize the cross-talk between the V P model and the wavelet, we rely on a carefully designed data weighting strategy, focusing on short-offset only as presented in Figure 3.10. The data weighting consists in using 400 m of offset at full amplitude, then ramping down to zero OBC VALHALL DATA at 1200 m offset. A tail mute is applied to remove the Schölte waves. Before being used for inversion, the wavelet is then manually checked and tapered to ensure its causality. With this methodology, the wavelet is only estimated at the beginning of each frequency band and kept fixed during FWI steps. The wavelet for the first frequency band and its associated spectrum is presented in Figure 3.10. We can see that no oscillations are present after 2.7 s, and wavelets generated from TOMO or 1D initial models are similar, which validates that this careful data selection mitigates the potential leakage in the wavelet estimation, which could come from the 3.3 Field data presentation and full-waveform inversion workflow P-wave velocity model. Finally, the wavelet spectrum is coherent with the data (2.5 to 5 Hz).

V P -init 1D (e) 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 4 1 5 Y (km) Depth (km) (d) Depth (km) < d < e (c)

Random shot subsampling and frequency continuation

The FWI workflow relies on a frequency continuation approach. Inversion is performed first on the frequencies ranging between 2.5 -5 Hz, then on a second band between 2.5 -7 Hz.

Respectively, the model grid spacing is set to 70 m and 50 m, ensuring at least five grid points for the smallest wavelength. These parameters are the same as in [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 approximation[END_REF] and OBC VALHALL DATA As previously mentioned, the dataset contains 2048 shot-gathers. Since our FWI code relies on source-parallelization, the complete dataset could not be reasonably fit in a single inversion on the HPC facility we have access to. Therefore, we rely on a source-subsampling strategy, similar to the one described in [START_REF] Warner | Anisotropic 3D full-waveform inversion[END_REF]. This makes it possible to divide the dataset into batches of pseudo-randomly selected shot-gathers (120 in our case). The pseudorandom selection implies that the previously used shot-gather could not be selected in the next subsample until all of the available ones are used once. The subsample of shot-gathers is changed every time the memory limit for l-BFGS is reached (3 iterations is the maximum memory of l-BFGS in this study). The source-subsampling strategy makes it possible to fit the FWI problem on relatively small HPC facilities while mitigating the acquisition foot-print in the reconstructed model.

Hierarchical data weighting strategy

The first data weighting strategy tried on the dataset for FWI is simple: only remove the Schölte waves and muted the trace near-zero offset in a radius of 350 m (see the approach in [START_REF] Kamath | Multiparameter full-waveform inversion of 3D OBC data from the valhall field[END_REF]. While this direct and straightforward approach (using all the data directly) could be justified with TOMO initial models as it is not supposed to generate cycleskipping, this approach could not be applied to tackle a crude initial model as the 1D one, even with the GSOT misfit function. Because of the strong cycle-skipping generated with the OBC VALHALL DATA 1D initial model, we need to rely on a more careful data weighting strategy to maximize the capability of the GSOT. This leads to introduce a six-step data weighting strategy, presented in Figure 3.11. The three first steps only focus on diving waves with a strict time windowing while restricting to offset to the first 4 km (later referred to as DW SO for Diving Wave -Short Offsets), 8 km offset (later referred to as DW MO for Diving Wave -Medium Offsets), and full offset (later referred to as DW FO for Diving Wave -Full Offsets). Then, the three next steps release the time windowing progressively, starting with 8 km offset (later referred to as RT1 MO for Release Time 1 -Medium Offsets), then on full offset (later referred to as RT1 FO for Release Time 1 -Full Offsets); to finally finish with full offset and a complete release of time windowing (later referred to as RT2 FO for Release Time 2 -Full Offsets). Please also note that for each data weighting, the Schölte waves are additionally muted, as we cannot model them in the acoustic approximation.

Field data presentation and full-waveform inversion workflow

Model post-processing

The last step of the FWI workflow is a post-processing applied to the updated V P model at the end of each FWI step. It consists in cutting the model using a stencil (based on the shape of the acquisition) and extrapolating the value outside the stencil with a nearest-neighbors algorithm. This is performed for each depth slice. The stencil shaped is adapted for each depth slice to consider the maximum illumination the acquisition can provide at depth. This aims to remove the area on the edges of the model that are never updated during FWI due to the lack of illumination and remove artificial low-velocity zones created on the border of the well-illuminated zone.

To be consistent and perform a fair comparison between our two starting models TOMO and 1D , the complete workflow detailed before is applied to both starting models. By doing so, we validate that our workflow can tackle the dataset properly and provide satisfactory results. Then, changing only the misfit function from L 2 to GSOT, we reconduct the complete inversion process to check how GSOT compares to L 2 in this controled environment. ROBUST FWI WITH GRAPH SPACE OPTIMAL TRANSPORT: APPLICATION TO 3D OBC VALHALL DATA

Full-waveform inversion results

Starting from the reflection tomography model

The first FWI results presented in this article are based on the TOMO initial model. The final reconstructed V P at 7 Hz is shown in Figure 3.12 using the L 2 misfit function, and in Figure 3.13 using the GSOT misfit function. Starting from the TOMO initial model, and as we do not expect cycle skipping, we keep the τ parameter of GSOT to 0.2 s. This value is low enough to always allow for fast convergence (similar to L 2 ) while preserving the potential improvement that GSOT can provide. Using a smaller value of τ would make results similar to L 2 while making τ larger would degrade the capacity of convergence of FWI due to a flatter attraction basin which is not needed in this case. The results are almost identical between both formulations on this setup. Some minor differences can, however, be observed. On the shallow depth slice (a) at 200 m, a small reduction of the acquisition imprint is observable using the GSOT misfit function. On the depth slice (c) at 1 km, a slight improvement of contrast between the low velocity anomaly and the sediment background is also visible. Finally, on the vertical slices (d-g), more lateral coherency in the geological structures is visible using the GSOT misfit function. The difference globally remains marginal, but this similarity between the L 2 and GSOT misfit functions results is satisfying in itself. Indeed, alternative misfit functions generally bring some drawbacks, such as loss of resolution power, which is not the case here. When performing QC on the V P profile extracted from the reconstructed model versus the sonic log (Fig. 3.22), we can see that the GSOT result is almost perfectly following L 2 result, which is consistent with the observation made directly on the V P model. A slight improvement of the well velocity fit at depth for GSOT misfit function can be observed for the well log 1 (which passes through the center of the target). Generally, FWI results improve the fit to the log compared to the initial model, showing that our FWI workflow provides robust and reliable results.

The data fit presented in Figure 3.15 shows that GSOT can improve over the L 2 FWI data fit (which is already very satisfying). On the rec A cable A 2D CRG (through the low velocity anomaly), we can observe some data fit improvement with GSOT, for example, at -6.5 km offset and 7 s. On the rec B cable B 2D CRG, in an area with relatively mild variations in velocities and well away from the "gas cloud", we can observe that, globally, the data fit is almost the same for L 2 and GSOT, with all phases correctly explained. This is not surprising as this CRG focuses on the part of the model with sedimentary geology only, which the initial model better explains. However, we can still see the advantages of GSOT this time with more arrivals present for late time, as visible in the -7 to -4 km offset, from 5 to 7 s. This "reference" result validates two essential aspects of this study. This first one is that our FWI workflow is adapted to the dataset and can provide satisfactory results when using L 2 -based FWI when a good enough initial model is used as the TOMO model. This point is of first importance as it validates that our FWI workflow is consistent with the literature on this dataset. The second conclusion is that the GSOT misfit function can tackle field data problems without significant issues. Modifying the misfit function and keeping all other parameters similar, GSOT can provide comparable results to those provided by the L 2 FWI. We can even see slight improvement both in terms of the P-wave velocity model estimation and data fit with the GSOT misfit function. These encouraging results make it possible to push the analysis forward: can GSOT help tackle the cycle-skipping issue and improve FWI

Starting from the 1D initial model

To validate the capability of GSOT to tackle large cycle-skipping, we use the 1D initial model that we introduced previously. This 1D initial model generates cycle-skipping, even on the mid-offset diving waves (mainly at -4 and -8 km offset), as clearly illustrated in Figure 3.8.

First, we compare the data fit obtained with L 2 and GSOT at the 2nd step of our FWI workflow as presented in Figure 3.16. For this first two steps of the FWI workflow, τ is set to 0.35 s, which is enough to handle the cycle-skipping generated in the data. First-order observation may indicate a similar data fit with GSOT and L 2 , but a more in-depth analysis shows several differences. At offset higher than 5 km, we observe a degradation with outof-phase arrivals for L 2 on the rec A cable A CRG, while GSOT results present a significant improvement on this part. Then, late arrivals events are better explained with GSOT. On the second rec B cable B CRG, we observe a more continuous reconstruction of the first events while also reducing out of phase one. Late arrivals are also better reconstructed with GSOT.

Full-waveform inversion results

On the second CRG, the data fit obtained with GSOT is already quite good for such an early stage of the inversion. In Figure 3.17 is presented the same CRG, but without applying the data weighting used at this early stage (DW MO data weighting), but instead the final relaxed data weighting (RT2 FO). Here, we can see that we are starting to predict data at larger offset when using GSOT compared to L 2 .

After the first two steps of FWI (DW SO and DW MO), we obtain the reconstructed model presented in Figures 3.18 L 2 reaches its limit and is likely affected by strong cycle-skipping. This is why we stop the L 2 inversion at this stage of our workflow: pushing forward the inversion by introducing more data does not help. Conversely, GSOT-based FWI can provide promising and meaningful V P updates, with the recovery of correct background velocities at depth and even key features such as the definition of strong low velocity anomalies (slices (c) (d) (f)).

As GSOT results are encouraging after only two passes on the first frequency band, we apply our complete workflow and perform the complete inversion similar to our reference inversion, with six passes on 5 Hz data and 7 Hz data. Regarding the τ parameter for GSOT, as said earlier, the first two steps on the first frequency band (2.5 to 5 Hz) were performed using τ = 0.35 s. The next four steps on the first band used a reduced τ = 0.2 s as it is enough to tackle the shift present in the data (which already illustrates the improvement that GSOT achieves in the first step of the inversion). For the second frequency band (2.5 to 7 Hz), only the first step of the workflow uses a slightly relaxed τ of 0.25 s, while the remaining five steps use τ = 0.2 s. Again, τ is relaxed at the begining to mitigate cycle-skipping, while a smaller τ is used afterward to preserve resolution power and speed up convergence.

We present the final data fit at 7 Hz for GSOT in Figure 3.20. We can see that a relatively good data fit is obtained, with most of the arrivals correctly explained. Still, some late arrivals are out of phase, mainly for large offsets (larger than 6 km) We obtain the final reconstructed V P presented in Figure 3.21. The results are promising and show a clear improvement in resolution compared to the early results at 5 Hz, with the main target structures retrieved above 2 km depth. The shape of the low velocity anomaly is correctly retrieved (slice (c) (d) (f)). Lateral resolution is very significantly improved, allowing the definition of narrow low-velocity (150 m wide) anomalies not resolved in the L 2 inversion (slice (e)).

When performing a QC using the sonic log presented in Figure 3.22, comparing the final GSOT results with the early L 2 results (that could not be pushed further with meaningful updates) clearly illustrates that GSOT-based FWI performs meaningful updates of the model that follows the trends of the V P sonic log. Conversely, the L 2 -based FWI updates are not correcting for the wrong initial model, and updates go in the wrong direction at several depths (for example, 1200 m and 1750 m). A comparison between the sonic log and GSOT-based FWI results for the two starting models is shown in Figure 3.23. For the two logs outside the center target area (Log 2 & 3), a good agreement of the FWI results is observed, with reconstructed V P models following the same trend. For Log 1 in the center of the target area, we can expect that reconstructed V P results will be incorrect below approximately ≈ 2 km due to the poor illumination. Above 1.4 km depth, a good agreement can be observed between the reconstructed models obtained using the two starting models. OBC VALHALL DATA Figure 3.17: Same as Figure 3.16 but with final relaxed data weighting (RT2 FO) for display. This exhibits the improvement of datafit obtained in area which are not yet inverted. Here, black arrows point to improve fit and coherency of the diving wave with GSOT, whereas yellow arrows point to improved data-fit of reflected events with GSOT.

The results are still not perfect, and the main issue is the presence of a low-velocity update around 500 m depth, as made visible on slice (b). This low-velocity update (not present in the reference FWI results starting from TOMO initial model) introduces a down vertical shift of layer under this perturbation. This is why slice (c) is extracted 100 m under (1.1 km instead of 1 km depth). This vertical shift does not affect the shape of the low velocity anomaly but only its depth. One possible way to avoid this artifact would be to modify the early stage of the inversion, for example, with different data-selections or modifications of the initial model for a slightly better one (for example, based on sonic log information). While this would have probably improved the results, we decided to keep a crude 1D initial model to stay as generic as possible and assess the capacity of GSOT in a setup without prior information.

Computational costs

Computational cost analysis is performed for one gradient estimation for both L 2 and GSOT misfit functions. Computation is performed on Haswell E5-2690V3@2.6 GHz Intel nodes containing 24 cores and 64 Gb of memory per node. We use 4 cores per source as our finite difference scheme uses an OpenMP parallelization.

The computational costs on the first frequency band (2.5 -5 Hz) and second frequency band (2.5 -7 Hz) are summerized in the table. The gradient column corresponds to the time spent to calculate the gradient (containing all required wave modelings), while the misfit column isolates the time spent computing the misfit function. The total time is the sum of gradient and misfit, and the ratio compares L 2 (put at 100%) with GSOT. This analysis shows that while the computational complexity of the solution of the gradient estimation scales to O(ω 4 )), the computation complexity of the GSOT computation is in O(ω 3 ), as noted in [START_REF] Métivier | A graph space optimal transport distance as a generalization of L p distances: application to a seismic imaging inverse problem[END_REF].

The overhead cost induced by GSOT is therefore reduced on the higher frequencies, which are the most expensive ones computationally speaking (for example, with a maximum frequency of 5 Hz, the first frequency band is relatively fast to compute, making a 20% overhead acceptable). This is one key feature that makes the GSOT misfit function able to tackle field data applications as higher frequency drastically increases the computational cost, and GSOT overhead will become smaller. OBC VALHALL DATA V P at 7 Hz -init 1D -GSOT

Discussion

Discussion

Facing challenging field data applications with alternative misfit functions is not widely documented in the literature. One of the only other alternative misfit functions that has been applied successfully to field data is adaptive waveform inversion (AWI) [START_REF] Warner | Robust adaptive waveform inversion[END_REF][START_REF] Ravaut | Imaging beneath a gas cloud in the north sea without conventional tomography[END_REF][START_REF] Debens | Full-bandwidth adaptive waveform inversion at the reservoir[END_REF][START_REF] Roth | Improving seismic image with high resolution velocity model from awi starting with 1d initial modelcase study barents sea[END_REF][START_REF] Guasch | Adaptive waveform inversion: Practice[END_REF][START_REF] Warner | Extending the velocity resolution of waveform inversion below the diving waves using awi[END_REF] or Kantorovich-Rubinstein optimal transport (KROT) [START_REF] Poncet | Fwi with optimal transport: a 3D implementation and an application on a field dataset[END_REF][START_REF] Messud | Multidimensional optimal transport for 3d FWI: Demonstration on field data[END_REF][START_REF] Sedova | Acoustic land full waveform inversion on a broadband land dataset: the impact of optimal transport[END_REF].

Regarding the final results obtained with the GSOT misfit function, we first validate that it can improve over L 2 -based FWI in a controled environment (starting from a good enough V P TOMO ) with improved datafit and improved lateral coherency and reduced acquisition imprint on the model side. The computational overhead induced by GSOT stays limited thanks to the computational complexity of GSOT being one order of magnitudes smaller than the complexity of the gradient estimation for a given maximum frequency. This behavior translates into a good scaling property when facing high-frequency data, which is now one trend for field data FWI applications. When tackling a difficult 1D initial model, which generates cycle-skipping, we show that GSOT-based FWI provides good model updates and good datafit where classical L 2 -based FWI fails. The control of the convexity provided by the GSOT misfit function is performed through the choice of the τ parameter. This parameter is directly linked to the observed time-shift between calculated and observed data, making it easy to tune, depending on the initial fit of the data and the expected maximum time shift. We adapt it from 0.2 s to 0.35 s in the initial stage of the workflow when we switch from the initial tomography model to the initial 1D model.

We also emphasize that results under 2 km suffer from a substantial lack of illumination and are therefore limited to the resolution of our initial model, explaining why no meaningful updates are present for depth superior to 2 km. This limitation in terms of depth reconstruction (under 2 km) is not surprising as the FWI alone is not expected to present enough illumination. In this case, it would require the use of reflected wave inversion (RWI) or joint full-waveform inversion (JFWI). Some preliminary but encouraging results were obtained combining GSOT and JFWI [START_REF] Provenzano | Joint FWI of diving and reflected waves using a graph space optimal transport distance: Synthetic tests on limitedoffset surface seismic data[END_REF]: GSOT adds the convexity necessary to predict the reflected data, enabling robust velocity updates of the model at depth.

Conclusion

In this work, we focus on the application of a new misfit function: graph space optimal transport. This formulation shows a clear improvement over L 2 in a controled environment, unleashing the potential of FWI to perform meaningful updates when starting from a crude, 1D initial model. This clear improvement in cycle-skipping robustness, combined with a simple setup (only one physical parameter to define) and a reasonable computational overhead, illustrates that GSOT is a good candidate to improve FWI robustness and therefore make FWI more accessible and easily applicable.

Chapter 4

Exploring new ideas and overcoming practical limitations

Overcoming practical issue for 3D KROT applications

This chapter presents some of the ideas and workaround that I have developed during my years of research. This could be seen as the "graveyard" of multiples ideas, that could seem interesting but were never pushed further for multiples reasons. Some ideas were simply not suitable, some were interesting but present problematic flaws, and better ones simply replaced some.

I think it is important to remind that science is not a straight path. It often requires going through multiple paths to have a better understanding of the possibilities. Unfortunately, the scientific community often focuses only on working solutions but not so much on what is not working. This leads researchers worldwide to repeatedly try the same "doomed-from-thebeginning" idea that nobody ever published about.

In a small and humble attempt to improve this, I present some of what I could refer to as "interesting failures", that still brings meaningful insights on the road ahead for better and more robust FWI by possibly removing some dead-ends for future research in this field.

Overcoming practical issue for 3D KROT applications

As illustrated in Figure 4.1, the Kantorovich-Rubinstein optimal transport (KROT) does not translate into drastic convexity improvement. The main benefit comes from the fact that it can be defined in 1D, 2D, or 3D, allowing us to consider the lateral coherency of shot gathers. By taking lateral coherence information of the data into account, KROT allows having a beneficial effect on the convergence of the FWI schemes. It also helps "regularize" the results, with a denoising effect that is welcome for field data application. This has been illustrated in recent publications [START_REF] Messud | Multidimensional optimal transport for 3d FWI: Demonstration on field data[END_REF][START_REF] Sedova | Acoustic land full waveform inversion on a broadband land dataset: the impact of optimal transport[END_REF][START_REF] Carotti | Optimal transport full waveform inversion -applications[END_REF].

For 3D field data, keeping computational costs low is essential, and it is essential to minimize the potential overhead introduced by the KROT formulation. To mitigate the computational overhead of 3D KROT, we rely on a new idea: extract lines of receivers (Figure 4.2) inside the 3D seismic cube [START_REF] Poncet | Fwi with optimal transport: a 3D implementation and an application on a field dataset[END_REF]. This approach creates 2D panels on which 2D KROT can be applied. The 3D residual will be reconstructed using all 2D residuals, and the cost function will be the summation of every panel cost.

The algorithm for this new optimal transport work as follows:

• Choose the direction of receivers line (in X or Y) and number of receivers lines to extract in the 3D cube (these parameters should make sense with the acquisition)

• A width of X or Y values for a line is calculated depending on the size of the acquisition

• A line is extracted from the 3D cube by selecting all receivers inside a given range of x rec or y rec values

• 2D KROT is applied on the panel. The cost function is increased by the value of the cost for this 2D panel, outputted 2D residual for the line is remapped to the 3D residual (each traces of 2D residual is put back where it was in the complete 3D residual)

• Code move to next line until all lines are done This yields a more flexible KROT implementation for 3D data, which presents a reasonable computational overhead (increase around 25% max compared to L 2 ) while preserving the ability of KROT to take advantage of lateral coherency of the data compared to 1D KROT.

The only potential issue is that multi-2D implementation assumes that the extracted 2D line is straight. This can be circumvented by performing receivers re-interpolation on a regular cartesian grid directly on the data, which allows performing a regular line extraction afterward (Figure 4.3).

Exploring new ideas of alternative misfit functions

In this section, I present some ideas that I have explored. These ideas have been explored mainly on an early "prototyping" stage, focusing on increasing convexity for time-shifts on simple 1D shifted patterns tests. These ideas have never passed this early stage and have never been implemented into FWI for different reasons that I will detail.

Instantaneous phase and unwrapped phase misfit functions

First, in this section, I re-introduce the concept behind instantaneous phase and envelope (as proposed by [START_REF] Bozdag | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF]. One way to separate the signal phase and envelope information is to rely on the analytical signal based on the Hilbert transform.

Hilbert transform and analytical signal

The analytical signal relies on the utilization of the Hilbert function H which can be defined in the time domain as

H[d(t)] = 1 π P +∞ -∞ d(τ ) t -τ dτ , (4.1)
where P stands for the Cauchy principal value.

In practical application, we do not rely on this time formulation to get the Hilbert transform but rather on a frequency formulation proposed by [START_REF] Marple | Computing the discrete-time" analytic" signal via fft[END_REF]. This method relies on a simple three-step approach:

1. Compute the Fourier transform of d(t) using an efficient FFT solver 2. Change the negative frequency to zero 3. Compute the inverse FFT This simple algorithm gives direct access to the analytical signal in a computationally efficient manner. From this analytical signal, we can simply retrieve the Hilbert transform as the imaginary part of the analytical signal d(t)

H[d(t)] = I[ d(t)] . (4.2)
This analytical signal can be written regarding the instantaneous phase φ(t) and the instantaneous envelope E(t) as

d(t) = E(t)e iφ(t) , (4.3)
where the instantaneous envelope is defined as (4.4) and the instantaneous phase as

E(t) = R[ d(t)] 2 + I[ d(t)] 2 ,
φ(t) = atan I[ d(t)] R[ d(t)] + ε , (4.5)
4.2 Exploring new ideas of alternative misfit functions or alternatively using the atan2 that has been introduced in Fortran and was originally intended to return a correct and unambiguous value for the angle in cartesian coordinates

φ(t) = atan2 I[ d(t)], R[ d(t)] . (4.6)
The atan2 formulation of instantaneous phase is prefered as it is uncondtionally stable and does not require a water level ε.

This decomposition between the instantaneous phase and instantaneous envelope can be used to define alternative misfit functions as detailed in the following. In chapter 2, we have already seen that instantaneous envelope could provide a potential improvement over the L 2 norm. In the following, I want to give more details on the instantaneous phase and some other ways to obtain "envelope-like" observables that could be used for alternative misfit functions.

Instantaneous phase misfit

Using the instantaneous phase defined on equation 4.6, an instantaneous phase misfit has been proposed by [START_REF] Bozdag | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF]. The motivation behind the utilization of the phase information is to reduce the sensitivity to amplitude variation, which mainly impacts global seismology. At the crustal scale, these amplitude problems are less present.

A representation of the instantaneous phase of a signal is shown on Figure 4.4. It exhibits that the instantaneous phase is oscillatory and could not increase cycle-skipping robustness as the width of the temporal support carrying information is now smaller than the original signal. This instantaneous phase misfit function is defined as

F IP = 1 2 T 0 [φ cal (t) -φ obs (t)] 2 dt . (4.7)
This misfit function is not designed to improve the convexity for time delays (Figure 4.5), but rather improve robustness to amplitude miss-match. It is sensitive to polarity information, and its resolution is prone to be good. One of the issues of removing amplitude information is that the weighting between noise and the arrivals will be the same, leading to an increased sensitivity to noise. From the instantaneous phase can be derived the unwrapped phase, which can potentially lead to an interesting formulation. The unwrapping of the phase is done by a simple algorithm: whenever the jump between consecutive angles is greater than or equal to π, the algorithm shifts the angles by adding multiples of ±2π until the jump is less than π.

Using an unwrapped phase, we can obtain a drastic improvement of convexity for time delays, which mimic NIM misfit function (Figure 4.6). Unfortunately, unwrapping of phase is an arduous task that is incredibly sensitive to noise (Figure 4.7), and is still challenging today in optics or electronic fields [START_REF] Wang | One-step robust deep learning phase unwrapping[END_REF]. It would be impossible to extract a clean and stable unwrapped phase of a field data signal. Therefore all the attempts to use the unwrapped phase were not further explored. This idea remains an excellent example of tempting ideas that could work well on synthetic traces but are doomed to failure when analyzed a bit closer. It represents the large amount of "dead-ends" that I have faced when trying to find alternative misfit functions.

Spline based envelope misfit function

The Hilbert transform is one way to extract the "true" envelope of the signal. By "true" it directly relates to the structure of a signal, which is composed of phase and amplitude, which translates into the fact that the instantaneous envelope can be recombined with the instantaneous phase to recreate the original signal. However, another possibility is to extract an "envelope" from the signal in a more physical sense, with no direct link to the phase.

To create this new envelope, we first need to extract all the maxima (and minima) of the signal; then, spline interpolation is performed between these maxima points (or minima point) to create a spline envelope approximation of the signal as proposed by [START_REF] Han | Spline envelope full-waveform inversion[END_REF]. Thus, we can define two separate envelopes: the upper part (from the maxima) and the lower part (from the minima). This is the first drastic difference to the unique envelope created by the Hilbert transform. These lower and upper envelopes can be merely average [START_REF] Hu | Multistep full-waveform inversion based on waveform-mode decomposition[END_REF] to create an "average" spline envelope (ASE). The original idea is performed using cubic spline interpolant. From the extensive test performed, a better version was obtained using a Piecewise Cubic Hermite Interpolating Polynomial (PCHIP).

The processing of the signal is performed in two steps. The first one (Algorithm 1) computes the first and second derivative of the signal, then a second part (Algorithm 2) extracts the top and bottom peaks. From the peaks, the PHCIP algorithm can be applied to get the top and bottom envelope.

An example of a spline envelope and average spline envelope is presented on Figure 4.8. In this case, on synthetic data, the average envelope obtained from the top and bottom envelope could make an excellent observable as the width of the temporal support has been increased while preserving phase and amplitude information.

The convexity improvement obtained with the ASE misfit is slightly improved compared to the instantaneous envelope, as shown in Figure 4.9. The advantages of ASE over IE would not only be related to slightly increased convexity but also to the fact that it integrates some polarity (phase) information that IE discards. This could potentially translate into improved results for practical application where phase information is crucial.

The effect of noise on the spline envelope is interesting. Test of data without noise vs.

Algorithm 1 calculate the first and second derivative of the signal and extract binary version of the first derivative 1: sn = s + small white noise 2: for i = 2 to N t -1 do 3:

ds(i) = (sn(i + 1) -sn(i -1))/(2 * dt) 4:
dds(i) = (sn(i + 1) -2 * sn(i) + sn(i -1))/(dt * dt) 5: end for 6: for i = 1 to N t do if dds(i) < 0 and sign(i) = sign(i + 1) then 5:

p top = p top + 1 6: t top (p top ) = t(i) 7: env top (p top ) = s(i) 8:
else if dds(i) > 0 and sign(i) = sign(i + 1) then The top peaks (blue circle) and bottom peaks (red circle) are extracted using algorithms 1 and 2. From these minima and maxima, the top envelope (blue line) and bottom envelope (red line) are extracted using a piecewise cubic Hermite interpolating polynomial. The mean of these two envelopes, corresponding to the average spline envelope (purple line), presents a nice shape that could be seen as a low-frequency version of the original signal.

two levels of noise is presented on Figure 4.10. It can be observed that for noiseless data (Figure 4.10a) or data with a small amount of noise (Figure 4.10b), the ASE present wider temporal support, which will translate into increased cycle-skipping robustness. We can even observe that the ASE common shot gather is slightly cleaner at zero offset for the signal with a small amount of noise. However, when strong noise is present (Figure 4.10c), the resulting signal is a denoized version, without the increased temporal support width. While this denoising effect is impressive on its own, it does not translate into improvement for cycle-skipping robustness.

The conclusion on ASE would be that the algorithms 1 and 2 present two main issues:

• It is not a differentiable operation, meaning that any attempt to derive an adjoint source from any misfit based on this operation will be unsuccessful.

• It relies on the addition of subtle white noise to stabilize the results. For synthetic data, there will be no issue. For real data with strong noise, the peaking will result in too many minima and maxima. When interpolated, this will only create a "denoise" version of the signal and not an envelope. Peak extraction is not a differentiable operation: this is a dramatic drawback and prevents us from applying this misfit function to FWI when using the adjoint state method to obtain gradient. Still, ASE properties seem interesting, and from toy experiment on traces shows it could be promising. The strong denoising effect could also be of interest for FWI as a pre-processing step to clean out data strongly contaminated by noise.

Diffusion based misfit function

Another idea explored for a new observable to improve the L 2 norm misfit function is a Diffusion-based observable. The simple idea is to take each trace and set it as the initial condition for a 1D diffusion problem solved by finite differences. By doing so, we hope to extract a "smooth" version of the observable with a wider temporal support. Since it is a diffusion problem, we can finely control and adjust the "diffusion time", which allows for precise control on how much the signal is diffused and therefore smooth. At first, we can diffuse for a long time (until we obtained a smooth enough observable to mitigate the cycle-skipping issue), then during the non-linear iteration of the FWI process, slowly reduce the diffusion time to make the observable less diffused and closer to the original signal. Ideally, in the end, it will be a non-diffusion misfit that will be similar to the classical L 2 waveform misfit. This ability to smoothly control the widening of temporal support translate into control of the cycle-skipping robustness, which is very appealing.

Let us remind the formulation of the diffusion equation as

∂u ∂t + κ ∂ 2 u ∂x 2 = 0 , (4.8)
where κ is the diffusion parameters.

In our case, since we do not care about the physical meaning of the diffusion we performed, we can set κ = 1. Since we set d(t) as an initial space condition, we can rewrite the diffusion equation as

∂d ∂τ + ∂ 2 d ∂t 2 = 0 . (4.9)
The new quantity D(t, τ ) is computed as the diffusion of the original signal d(t) for a given time τ (the diffusion time). If τ = 0, we do not diffuse and D(t, τ = 0) = d(t). The more τ increases, the more D(t, τ ) will be a diffused version of d(t).

Instead of considering the least-squares misfit (L 2 norm) between d obs and d cal , we consider the least-squares difference between D obs and D cal . An illustration of the effect of diffusion on a signal is visible on Figure 4.11. This clearly illustrates the increased temporal support for the diffused version compared to the original signal. Also, it shows the ability to control "how much" increased temporal support is possible by modifying the diffusion "time". This potentially translates into increased convexity of the misfit function with respect to time-shift as shown in Figure 4.12.

F dif f usion (τ ) = 1 2 T 0 [D cal (d cal , t, τ ) -D obs (d obs , t, τ )] 2 dt . ( 4
This behavior is further illustrated on Figure 4.13, with two Ricker wavelets, one unfiltered and one high passed above 2 Hz. For the unfiltered Ricker wavelet, we can see a clear improvement in the temporal support width, which is the feature we want to observe. However, for the highpass filtered version, it can be observed that the signal shape is not preserved, while the associated spectrum is still good and coherent. This raises concern with the suppression of zero-frequency information. To push the analysis deeper, noise can be added as presented on Figure 4.14. Here, it clearly illustrates that adding white noise to the data starts to interfere with the diffused version and altered the signal. To validate this hypothesis, we observed the effect of diffusion on a field-data trace (Figure 4.15). It can be observed that the spectrum of the diffused data is not clean and does not represent a pure shift down in frequency. This behavior can be linked to the fact that diffusion is expected to behave like a linear low-pass filter, meaning it amplifies information from the pre-existing low frequencies in the data. When looking at a high-passed Ricker wavelet or a field data (which contains slight noise and no zero frequency information), it can be observed that there is no added benefit of using diffusion to create a new observable.

This idea illustrates the potential false "good idea" which can be appealing at first while presenting no advantages to field data application. It is a good reminder that playing with synthetic data and especially Ricker wavelet (which contains close to zero frequency information) can introduce a skew conclusion. This kind of artificial boosting of a Ricker wavelet's lowfrequency content has been observed some times on some ideas presented during conferences and should raise a warning.

Exploring new ideas of alternative misfit functions

Beat-tone misfit function

Finally, another explored idea is to generate a beat tone using a cosine function on the data.

To do so, we multiply our signal d(t) in time with a cosine at a given frequency f c :

B(t) = d(t) * cos(2πf c t) (4.11)
If we have a signal of a given central frequency f , the beat generated should have a high and a low-frequency component such as f -f c and f + f c . We then have to remove the highfrequency component by applying a low pass filter, for example keeping only the frequencies under the central frequency f of the original signal. This would theoretically allow shifting the frequencies of the signal toward the low frequency. This behavior is illustrated on Figure 4.16 using a simple Ricker wavelet at 10 Hz. In this case, introducing a carrier at f c = 2 Hz results in a split of the frequency content, with a part around 5 Hz on the low-frequency side and one around 15 Hz for the high-frequency side. After applying a low-pass filter under 7 Hz we obtain a low-frequency version of our original signal, which looks appealing.

Again, this should theoretically translate into increased convexity of the misfit function as illustrated in Figure 4.17.

To further illustrate the effect, we apply beat tone with f c = 4 Hz on a trace from field data from Valhall as presented on Figure 4.18. The signal is then low-passed with a maximum cutoff frequency of 5 Hz. The original signal central frequency is around 6 Hz, and after applying the beat tone, we again have a split of the spectrum with peaks around 2.5 and 10 Hz. After the low-pass filter, we obtain a low-frequency version of the original signal. However, one of the first issues would be the sensitivity to the carrier frequency f c as illustrated on Figure 4.19. Here, with f c = 12 Hz, no beat tone effect is present. There is no "automatic and easy" way of finding the proper carrier frequency to optimize the beat tone effect from the experiment. Moreover, a second issue is also present and corresponds to the very unpredictable behavior of beat tone as it depends on the cosine phase. This effect is illustrated in Figure 4.20. A slight shift in time could lead to drastic modification of the beat tone version. This unstable behavior will effectively lead to a basin of attraction as large as the classical L 2 misfit. We can observe here that as the carrier is always the same, a possible shift in time of the wavelet will dramatically change the signal shape.

Investigation published during my thesis

Investigation published during my thesis

Here, I want to present a list of the investigations I have performed and published during my thesis that contains meaningful insights on robust FWI and alternative misfit functions. Some of the ideas are presented in-depth in this thesis, but some were put aside for the clarity of the manuscript. I would refer interested readers to my work for broader pictures on robust FWI and alternative misfit functions. • Improving cycle-skipping robustness with GSOT-based FWI on Valhall: Pladys, A., Brossier, R., and Métivier, L. (2020). Graph space optimal transport based FWI: 3D OBC valhall case study. In SEG Technical Program Expanded Abstracts 2020.

• Published paper on alternative misfit functions for FWI: Pladys, A., Brossier, R., [START_REF] Pladys | On cycle-skipping and misfit functions modification for full-wave inversion: comparison of five recent approaches[END_REF]. On cycle-skipping and misfit functions modification for full-wave inversion: comparison of five recent approaches. Geophysics, in press.

• Submitted paper on 3D OBC field data application of GSOT-based FWI: Pladys, A., Brossier, R., [START_REF] Pladys | Robust fwi with graph space optimal transport: application to 3d obc valhall data[END_REF]. Robust fwi with graph space optimal transport: application to 3d obc valhall data. Geophysics, submitted.

I also wanted to express my regrets for not sharing my work performed on a 3D OBC North Sea dataset in collaboration with Equinor. Unfortunately, the ownership of the dataset was transferred during my thesis to another company, and all the authorizations were lost. This was my first time getting around field data, and the application focused on 1D KROT analysis, then on a multi-2D KROT analysis.

Conclusion and perspectives Conclusion

This thesis manuscript aims at comparing alternative misfit functions for full-waveform inversion. In the beginning, a general introduction of the history of seismic imagining and, more generally, imaging in geophysics is performed. This introduces the different methods developed along the year and introduces the method at the core of this manuscript: full-waveform inversion (FWI). The pros and cons of this method are developed and make a transition toward the first Chapter. In Chapter 1, a generic introduction to the FWI formalism is made, from the forward modeling part of the problem to the inversion part with the multiple formalisms that can be used to solve inverse problems. The non-convexity issue is introduced, followed by a presentation of the cycle-skipping issue. An introduction to possible solutions to mitigate FWI non-convexity is then performed. In Chapter 2, alternative misfit functions are benchmarked on several synthetic tests, from simple 1D examples to more complex 2D FWI cases to assess the pros and cons of each proposition. In Chapter 3, an application of the most promising method (graph space optimal transport) is applied to a 3D field dataset. Finally, in Chapter 4, some ideas for other alternative misfit functions and some workaround performed for practical applications are presented.

Throughout this work, several important points have been discussed, and I want to reemphasize them as a conclusion:

• Is robustness to cycle-skipping the only important feature of an alternative misfit function?

As introduced in Chapter 1, the primary and central issue of FWI comes from the nonconvexity of the misfit function conventionally used: the least-squares norm (L 2 norm). This L 2 norm is at the origin of the cycle-skipping issue. Therefore, when the idea came of replacing the L 2 norm with alternative misfit functions, it seemed obvious -at firstto focus on reducing the non-convexity of the misfit function. This led to the proposition of many more convex (with respect to time-shifts) misfit functions. This should have theoretically solved the cycle-skipping issue and improve FWI robustness. Nonetheless, the number of field data applications using alternative misfit functions remains small compared to the number of proposed solutions.

This discrepancy is, for me, explained by the fact the L 2 norm is pretty good for all other aspects: it preserves signal shape (no alteration in amplitude/phase information), it exhibits strong noise robustness, it -when not affected by cycle skipping -delivers very high-resolution results that FWI is known for and finally, the computational cost CONCLUSIONS AND PROSPECTS is negligible. These advantages have made the L 2 norm the "way to go" for FWI. By looking at cycle-skipping robustness alone while not preserving the other important key features that a misfit function should exhibit, propositions of alternative misfit functions can reveal themselves not adapted to realistic FWI usage (such as field data), even if they reduce cycle-skipping sensitivity "on paper".

• How to assess the real world capacity of an alternative misfit function? Assessing misfit functions in the most generic and reproducible way is important. The literature misses much cross-comparison between several propositions. As an attempt to address this issue, I proposed several synthetic cases (presented in Chapter 2). From 1D toy example that only focuses on shifted patterns, to simple 2D FWI tests in transmission and reflection, to finally more realistic 2D cases with Marmousi, with a shared parameterization between misfit functions, cross-comparison is made easier and understanding potential drawbacks of a given formulation is made simpler. This allowed us to restrict the potential good down to what is for us the best candidate and further test it on field data in Chapter 3.

If an alternative misfit function is directly tested on a more realistic case, it can be easily tweaked to provide artificially good results, or it can fail and be difficult to identify why. I want to promote a more systematic approach to validate alternative misfit functions capacity and hope this work is a step in this direction.

• Are alternative misfit function some kind of magic button that could make FWI automatic? After extended tests, using what is for me the best misfit function at my disposal, from simple 2D FWI setup to 3D field data application, it never behaves as a "magic button". Even if, on simple synthetic cases, the increase of robustness that graph space optimal transport (GSOT) provides can simplify the workflow or allows to start from really crude initial models, extended tests have shown that it still requires the use of appropriate workflows when trying to tackle more complex situations. This is not a drawback, as in such cases, L 2 -based FWI simply fails to provide results even when using all the possible approaches. For me, GSOT is a clear improvement to FWI robustness, which constantly (for the tests I performed) improves results over L 2 , but at a given computational overhead. FWI remains a complex imaging process that is far from automatic, and providing tools to increase its robustness is a step in the good direction but not a magic leap toward a completely robust tool.

Perspectives

This work can be separated into two parts. The first one is a presentation of several misfit functions with an attempt at cross-benchmarking them. The analysis can be pushed forward on multiple aspects:

• Non-exhaustive selection of misfit functions. Indeed, the selection of misfit functions that are benchmarked together is non-exhaustive. Moreover, we compared some of the best challengers to some simpler propositions, which where never pushed further on. This was made as an attempt to exhibit pros and cons and to create more systematic ways of comparing misfit functions for FWI. By focusing on time-domain FWI formulation and only on misfit functions, we restricted ourselves from possible good candidates who rely on other kinds of reformulation of the FWI problem.

Other misfit functions could be considered and benchmarked against our selection using the same framework.

• Exploring and benchmarking other types of solutions.

In this thesis, my focus has never been on extension strategies. Nevertheless, it could be of interest to compare them with alternative misfit functions directly. Making a direct comparison of alternative misfit function to extension strategies could bring new insight into robust FWI.

The second part of this work consists of a field data application comparing GSOT to L 2 norm. Here, several possibilities should be considered to improve the results and push the analysis forward:

• Considering other misfit functions for field data.

From the analysis performed on several misfit functions, we decided to focus only on the most promising candidates at our disposal. Moreover, being the solution developed in the team, we had the greatest confidence with it, which is important when tackling complex and large datasets that require substantial computational resources. For the amount of time at my disposal, pushing the analysis on the best candidate make sense. Future works, comparing other good alternative misfit functions, such as AWI, would bring more insight into their behaviors when facing field data.

• Improving results using complementary methods.

As it was shown, results obtained on the Valhall dataset when starting from the 1D model were impressive. Still, some imperfections in the shallow layer reconstruction are visible compared to the results obtained from the reflection tomography initial model. These kinds of issues could be tackled with other methods such as RWI or JFWI [START_REF] Bibliography Zhou | Full Waveform Inversion of Early Arrivals and Reflections for Velocity Model Building and Case Study with Gas Cloud Effect[END_REF][START_REF] Zhou | Joint full waveform inversion of early arrivals and reflections: a real obc case study with gas cloud[END_REF][START_REF] Brossier | From rwi to jfwi: including diving waves in reflection-based velocity model building[END_REF][START_REF] Zhou | Velocity model building by waveform inversion of early arrivals and reflections: a 2d ocean-bottom-cable study with gas cloud effects[END_REF][START_REF] Li | Joint FWI for imaging deep structures: A graphspace OT approach[END_REF][START_REF] Provenzano | Joint FWI of diving and reflected waves using a graph space optimal transport distance: Synthetic tests on limitedoffset surface seismic data[END_REF][START_REF] Yao | A review on reflection-waveform inversion[END_REF]. Also, considering a more complex physic for the modeling could help explain the data better. This would provide an exciting insight into the capacity of alternative misfit functions for more complex workflow.

• Pushing inversion further.

The last line of investigation would be to push the inversion further. The choice of stopping at 7 Hz was mainly motivated by computational resources at our disposal. This link with the previous point and pushing the FWI workflow further will bring more knowledge on alternative misfit function behavior on field data.

I see this work as a stepping stone to encourage more rigorous and systematic analysis and cross-comparison of alternative misfit functions for FWI. Pushing FWI forward by increasing its robustness can only make this method more appealing, while unlocking targets that are still out of reach for now because of lacking good initial model. We show that field data applications using one of the last and most promising methods (GSOT) is feasible and yield very encouraging results and hope that FWI will continue moving forward with more interesting applications in the years to come.

List of Figures

Figure 2 :

 2 Figure 2: Schematic of an onshore and offshore setup in seismic exploration surveys (Energy Information Australia).

Figure 3 :

 3 Figure 3: Refraction/reflection FWI results at exploration scale. (a) correspond to the initial model used to start FWI workflow. (b) final 25 Hz FWI results. From Dickinson et al. (2017).

INTRODUCTIONFigure 4 :

 4 Figure 4: Small scale FWI results at 65 Hz. Here the target is the Ettlingen trench near Karlsruhe, Germany. The low-velocity anomaly in the center corresponds to fill-in materials used to level the ancient trench. From Irnaka, 2021.

Figure 5 :

 5 Figure 5: Example of global FWI for vertically polarized shear-wave-speed perturbations in GLAD-M25 at 250 km depth[START_REF] Lei | Global adjoint tomography-model glad-m25[END_REF]. Remarkable hotspots are highlighted with white arrow.

  Contents 1.1 A rapid overview of full-waveform inversion development . . . . . 1.2 Formulation of FWI . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 The forward problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Solution to the FWI problem . . . . . . . . . . . . . . . . . . . . . . . 1.3 Non-convexity in FWI . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Non-convexity and cycle skipping limitations . . . . . . . . . . . . . . 1.3.2 Classical approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Reformulating the FWI problem . . . . . . . . . . . . . . . . . . . . .

Figure 1

 1 Figure 1.2: Couple of source/receiver of a fixed diffraction point in a homogeneous medium.

  Mathematically, FWI can be formulated as finding the minimum of a misfit function f defined as min m f (m) = F (d cal [m], d obs ) , (1.2) where F is a generic function measuring the distance between d obs and d cal [m]. Under general notation, d cal [m] is obtained through the extraction of the values of wavefield at the receivers location such that d cal [m] = Ru[m] , (1.3) FULL-WAVEFORM INVERSION where R is an extraction operator and u[m] is the solution of the wave propagation problem A[m]u = b , (1.4) with A[m] a generic wave propagation operator (from acoustic to visco-elastic).

  .12) FULL-WAVEFORM INVERSION which now involves only 5 independent coefficients due to the relation c 12 = c 11 -2c 66 .

Figure 1

 1 Figure 1.3: Representation of misfit function f [m] : Ω ⊂ R 2 → R presenting three local minima, two of which are global and denoted by a red dot and a third one denoted by a blue dot. From Thurin (2020).

1. 3

 3 Non-convexity in FWI 1.3 Non-convexity in FWI 1.3.1 Non-convexity and cycle skipping limitations 1.3.1.1 The least-squares norm misfit function

Figure 1 . 5 :Figure 1

 151 Figure 1.5: Simple time-shifted Ricker wavelet (4 Hz) test for L 2 misfit function. (a) presentation of the setup, with reference signal (seen as d obs ) in solid black, and shifted signal in dashed black (seen as d cal ). (b) Value of the L 2 norm with respect to time-shift.

Figure 1

 1 Figure 1.7: FWI schematic example for cycle skipping: traces for a source-receiver couple at z = 500 m for d obs (solid black), d cal in the initial model (dashed black) and d cal in final reconstructed V P model (solid red). (a) correspond to initial model 1, (b) to initial model 2.

FWIFigure 1

 1 Figure 1.8: FWI schematic example for cycle skipping: reconstructed V P models using FWI starting from (a) homogeneous intial model 1, (b) homogeneous initial model 2.

Figure 1

 1 Figure 1.9: Same 1D setup as Figure 1.5, but this time comparing the influence of frequency on the convexity. (a) 0.5 Hz Ricker wavelet, (b) 2 Hz Ricker wavelet, (c) 4 Hz Ricker wavelet. (d) L 2 misfit function values with respect to time-shifts for the three different frequency cases. The width of the basin of attraction of the misfit function is enlarged for lower frequencies.

  .1) where the subsurface parameters are denoted by m, d obs is the observed data, d cal [m] is the synthetic data, and F is a generic function measuring the misfit between d obs and d cal . Under general notation, d cal [m] is obtained through the extraction of the values of wavefield at the receivers location such that d cal [m] = Ru[m] , (2.2) where R is an extraction operator and u[m] is the solution of the wave propagation problem A[m]u = b , (2.3) with A[m] a generic wave propagation operator (from acoustic to visco-elastic).

  Donno et al. (

d

  cal (x r , t)dx r dt = xr T 0 d obs (x r , t)dx r dt .

Figure 2

 2 Figure 2.1: Comparison of several misfit functions in a simple 1D case for one shifted arrival. The arrival is set to be a Ricker wavelet with a central frequency of 4 Hz. (a) represents the signal used for the test (with only one arrival at the center). The fixed reference signal is displayed in continuous black. The shifted signal is displayed in dotted black (here for +1.5 s). (b) represents the normalized misfit function values with respect to the time-shift (from -1.5 s to 1.5 s).

Figure 2

 2 Figure 2.2: (a) quantities Q obs and Q cal for three time-shifts (-1.5 s in green, -0.1 s in blue and "in phase" in dashed red). (b) the area under the curve for Q obs -Q cal for the three time-shifts.

Figure 2 Figure 2 Figure 2

 222 Figure 2.3: E obs and E cal for three time-shifts (-0.5 s, -0.2 s and in phase).

Figure 2

 2 Figure 2.6: AWI analysis with two setups: test A and test B. (a) test A both wavelets shift, (b) test B only the left wavelet is shifted, the right one being always in phase (similarly to Figure 2.5). (c) shape of AWI misfit function with respect to time-shift in both cases. The Wiener filters presented under are shown for a time-shift of -0.5 s (black cross on the misfit). (d) Wiener filters (w(t)) and the Gaussian penalty function P(t). (e) the Wiener filters multiplied by the penalty function.

  Figure 2.7: FWI Test 1: (a) true model, (b) initial model 1 with V P = 1300 m.s -1 , (c) initial model 2 with V P = 1700 m.s -1 and (d) initial model 3 with V P = 1900 m.s -1 .

Figure 2

 2 Figure 2.9: FWI Test 1: Extracted traces along the shortest path (horizontal straight line at 500 m depth passing through the spherical inclusion). Observed data are in solid black, synthetic data in the initial model in dashed black and final reconstructed synthetic data in solid red. First column corresponds to initial model 1, second column to initial model 2 and third column to initial model 3. Each line corresponds to reconstructed V P model from FWI using respectively L 2 (a-c), AWI (d-f), IE (g-i), NIM (j-l), KROT (m-o) and GSOT (p-r) misfit functions.

Figure 2 Figure 2 Figure 2

 222 Figure 2.10: FWI Test 2: (a) Homogeneous initial model (V P = 1500 m.s -1 ). (b) True model with a negative V P anomaly, (c) true model with a positive V P anomaly.

Figure 2 .Figure 2 Figure 2 Figure 2

 2222 Figure2.13: FWI Test 2: Extracted traces for a center shot at zero offset. Observed data are in solid black, synthetic data in the initial model in dashed black, and final reconstructed synthetic data in solid red. The left column corresponds to a negative V P anomaly while the right column to positive V P . (a,b) shows the complete traces. Each subfigures under are cropped on the magenta box to emphasize the polarity reversal introduced by the layer. They correspond to traces calculated in the reconstructed V P model using respectively L 2 (c,d), AWI (e,f), IE (g,h), NIM (i,j),KROT (k,l) and GSOT (m,n) misfits functions.

Figure 2

 2 Figure 2.17: Test 3.a: (a) True Marmousi II model. (b) S250 initial model, (c) S500 initial model and (d) 1D initial model.

  Figure 2.18: Test 3.a: Density model obtained from V P using Gardner's law for: (a) True Marmousi II model, (b) S250 initial model and (c) S500 initial model.

  Reconstructed models from the 1D initial model are presented in the right column of Figure 2.20. The associated data-fit obtained after FWI are presented in Figure 2.22.

Figure 2 .

 2 Figure 2.19: Test 3.a: Common shot gather (CSG) for field data overlapped by synthetic data in (a) S500 initial model, and (b) 1D initial model. Field data in black and white, overlapped by red to blue synthetic data with transparency. Red and white visible mean out of phase, black and blue mean in phase.

Figure 2

 2 Figure 2.20: Test 3.a: Inverse crime FWI final reconstructed V P model for Marmousi. Left column corresponds to S500 initial model, right column to 1D initial model. The lines respectively correspond to the final reconstructed V P model using L 2 (a,b), AWI (c,d), IE (e,f), NIM (g,h), KROT (i,j) and GSOT (k,l).

Figure 2 .

 2 Figure 2.21: Test 3.a: Overlapped common shot gathers for synthetic data in the final reconstructed V P model starting from S500 initial model vs field data. Each subfigure corresponds to misfit function, with L 2 (a), AWI (b), IE (c), NIM (d), KROT (e), and GSOT (f).

Figure 2 .

 2 Figure 2.22: Test 3.a: Same as Figure 2.21 but starting from 1D initial model.

Figure 2

 2 Figure 2.23: Test 3.a: Costs evolution in the inverse crime Marmousi for 1D initial model. (a) evolution of cost functions over iterations, (b) true L 2 cost evolution over iterations, (c) model error reduction over iterations and finaly (d) model error vs. the data error reduction.

2. 5

 5 Figure 2.24: Test 3.b: Comparison of true wavelet used to generate the observed data (for both test 3.a and 3.b) versus the inverted wavelet used for inversion in more realistic test 3.b. Under the associated amplitude spectrum of both true wavelet and estimated wavelets.

Figure 2 Figure 2

 22 Figure 2.25: Test 3.b: More realistic FWI final reconstructed V P model for Marmousi. Left column corresponds to S250 initial model, right column to S500 initial model. The lines respectively correspond to the final reconstructed V P model using L 2 (a,b), AWI (c,d), IE (e,f), KROT (g,h) and GSOT (i,j).

Contents 3 . 1

 31 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Field data presentation and full-waveform inversion workflow . . . 3.3.1 Geological situation, dataset and initial models . . . . . . . . . . . . . 3.3.2 Full-waveform inversion workflow . . . . . . . . . . . . . . . . . . . . . 3.4 Full-waveform inversion results . . . . . . . . . . . . . . . . . . . . . 3.4.1 Starting from the reflection tomography model . . . . . . . . . . . . . 3.4.2 Starting from the 1D initial model . . . . . . . . . . . . . . . . . . . . 3.4.3 Computational costs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . This chapter directly comes from the submitted paper Pladys et al. (2021a).

  Figure 3.1: Location of the Valhall field on the North Sea (from Thurin, 2020)

Figure 3

 3 Figure 3.2: Layout of the Valhall acquisition overlapped on an horizontal V P slice at 1 km obtained by GSOT-based FWI (from this study). Location of sources (gray dots) and receivers (blue diamonds). Two receivers positions (A and B) are located with black stars. Cables A (x = 2950 m), B (x = 5530 m) and C (x = 3080 m) are identifed. Black dots denote the position of the three V P well sonic-logs.

Figure 3

 3 Figure 3.3: 2D common-receiver gathers extracted for receiver A along cable A for: raw data (top), band 1 data 2.5 -5 Hz (middle) and band 2 data 2.5 -7 Hz (bottom). White dashed arrows point on the Schölte waves which are muted for the inversion. Blue and yellow arrows respectively point on the reflexion from the low velocity zone and the reflexion on the top of the reservoir. Black arrows point on the diving waves.

Figure 3

 3 Figure 3.4: Slices of the initial model V P TOMO . (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km.

Figure 3

 3 Figure 3.5: Same as Figure 3.4 for initial model V P 1D .

Figure 3

 3 Figure 3.6: Same as Figure 3.4 for anisotropic parameter η.

3. 3

 3 Figure 3.7: 2D common-receiver gathers at 5 Hz. Synthetic data generated into the initial model TOMO are displayed in a blue/white/red color scale, field data are overlapped in grayscale with transparency. The best result is achieved when black and blue are the only colors visible. Red and white are shown when data are not in phase. (a) receiver A along cable A (through the low velocity anomaly). (b) receiver B along cable B.

Figure 3 . 8 :

 38 Figure 3.8: Same as Figure 3.7 into the 1D initial model.

Figure 3

 3 Figure3.9: Complete FWI workflow used on the OBC Valhall dataset. At the core of the process lies the FWI iterations. Then several loops are nested one into another, from internal FWI iterations, source-subsampling, data selection, to finally the outer one of frequency continuation.

Figure 3 Figure 3

 33 Figure 3.12: Slices of the 7 Hz FWI reconstructed V P using L 2 misfit function starting from TOMO initial model. (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km.

Figure 3

 3 Figure3.14: Comparison of V P profiles extracted from the TOMO initial model (dashed red), 7 Hz FWI models using L 2 (solid yellow) and 7 Hz FWI model using GSOT (solid purple) with sonic log (solid black). Left subfigure corresponds to the Log 1 at the center of the target. Middle subfigure to the Log 2, and right subfigure to Log 3 (which is far away from the target).

  Figure 3.15: 2D common-receiver gathers at 7 Hz starting from the TOMO initial model. Synthetic data (blue/white/red color scale) generated into the final reconstructed V P using (a,b) L 2 misfit function, (c,d) GSOT misfit function. (a,c) receiver A along cable A (through the low velocity anomaly). (b,d) receiver B along cable B. Field data are overlapped in grayscale with transparency. Black arrows point to improvement obtained with GSOT.

Figure 3 .

 3 Figure 3.16: 2D common-receiver gathers at 5 Hz starting from the 1D initial model with data weighting apply (DW MO). Synthetic data (blue/white/red color scale) generated into the final reconstructed V P at 2nd workflow step using: (a,b) L 2 misfit function, (c,d) GSOT misfit function. (a,c) receiver A along cable A (through the low velocity anomaly). (b,d) receiver B along cable B. Field data are overlapped in grayscale with transparency. Black arrows point to area where GSOT improves the datafit.

Figure 3

 3 Figure 3.18: Slices of the 5 Hz FWI reconstructed V P using L 2 misfit starting from 1D initial model. (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1.1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km.

Figure 3

 3 Figure 3.19: Same as Figure 3.18 using GSOT. Here, results are consistant until 1.5 to 2 km depth compared to L 2 -based FWI. Characteristic structures of the Valhall field are recovered.

Figure 3 .

 3 Figure 3.20: 2D common-receiver gathers at 7 Hz starting from the 1D initial model. Synthetic data (blue/white/red color scale) generated into the final reconstructed V P using GSOT. (a) receiver A along cable A (through the low velocity anomaly). (b) receiver B along cable B. Field data are overlapped in grayscale with transparency.

Figure 3

 3 Figure 3.21: Slices of the 7 Hz FWI reconstructed V P using GSOT misfit starting from 1D initial model. (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1.1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km.

Contents 4 . 1

 41 Overcoming practical issue for 3D KROT applications . . . . . . . 4.2 Exploring new ideas of alternative misfit functions . . . . . . . . . 4.2.1 Instantaneous phase and unwrapped phase misfit functions . . . . . . 4.2.2 Spline based envelope misfit function . . . . . . . . . . . . . . . . . . . 4.2.3 Diffusion based misfit function . . . . . . . . . . . . . . . . . . . . . . 4.2.4 Beat-tone misfit function . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Investigation published during my thesis . . . . . . . . . . . . . . .

Figure 4

 4 Figure 4.1: Simple time-shifted Ricker wavelet test for crosscorrelation misfit function. (a) presentation of the setup, with reference signal (seen as d obs ) in solid black, and shifted signal in dashed black (seen as d cal ). (b) Value of the KROT misfit function with respect to timeshift.

Figure 4 . 3 :

 43 Figure 4.3: Same as Figure 4.2, but for the interpolated dataset. Here the extracted 2D lines are now perfectly straight, with regular receivers spacing.

Figure 4 Figure 4 . 5 :

 445 Figure 4.4: Representation of a simple signal d(t) (black line) and the associated instantaneous phase φ(t) (red line).

Figure 4 Figure 4 . 7 :

 447 Figure 4.6: Simple time-shifted Ricker wavelet test for the unwrapped phase misfit function. (a) presentation of the setup, with reference signal (seen as d obs ) in solid black, and shifted signal in dashed black (seen as d cal ). (b) Value of the unwrapped phase misfit function with respect to time-shift.

  12: end for EXPLORING NEW IDEAS AND OVERCOMING PRACTICAL LIMITATIONS Algorithm 2 extract top and bottom peaks 1: p top = 0 2: p bot = 0 3: for i = 2 to N t -1 do 4:

Figure 4 . 8 :

 48 Figure 4.8: Spline envelopes and average spline envelopes for a given seismic trace (solid black).The top peaks (blue circle) and bottom peaks (red circle) are extracted using algorithms 1 and 2. From these minima and maxima, the top envelope (blue line) and bottom envelope (red line) are extracted using a piecewise cubic Hermite interpolating polynomial. The mean of these two envelopes, corresponding to the average spline envelope (purple line), presents a nice shape that could be seen as a low-frequency version of the original signal.

Figure 4

 4 Figure 4.9: Simple time-shifted Ricker wavelet test for ASE misfit function. (a) presentation of the setup, with reference signal (seen as d obs ) in solid black, and shifted signal in dashed black (seen as d cal ). (b) Value of the ASE and IE misfit functions with respect to time-shift.

  Figure 4.10: Extraction of average spline envelope based using PCHIP for three levels of noise (a: without noise, b: small amount of noise, c: strong noise).

  .10) where D cal = D(d cal , t, τ ) and D obs = D(d obs , t, τ ) are modified versions of the wavefield d cal (t) and d obs (t) based on the diffusion equation for a given diffusion time τ .

Figure 4

 4 Figure 4.11: Test of diffusion on a simple signal (solid red) for two diffusion values. The less diffused version (solid purple) and the two times more diffused version (solid blue).

Figure 4 Figure 4

 44 Figure 4.12: Simple time-shifted Ricker wavelet test for crosscorrelation misfit function. (a) presentation of the setup, with reference signal (seen as d obs ) in solid black, and shifted signal in dashed black (seen as d cal ). (b) Value of the diffusion-based misfit function with respect to time-shift. Here the diffusion time τ is set to 3000.

Figure 4 Figure 4

 44 Figure 4.14: Same setup as Figure 4.13 but with noise added on the signal.

Figure 4

 4 Figure 4.16: Beat tone illustration on a simple unfiltered Ricker wavelet and associated spectrum. First line corresponds to the raw Ricker wavelet, the second line to the beat tone version and the third line to a low-passed version of the beat tone version to remove high-frequency content.

Figure 4 .Figure 4 Figure 4 Figure 4

 4444 Figure 4.17: Simple time-shifted Ricker wavelet test for crosscorrelation misfit function. (a) presentation of the setup, with reference signal (seen as d obs ) in solid black, and shifted signal in dashed black (seen as d cal ). (b) Value of the beat-tone misfit function with respect to time-shift. Here f c = 0.5 Hz and the low-pass filter is set to 1 Hz.

  Some of the possible ray paths for seismic waves passing through the Earth(from Fowler et al., 1990). (b) Jeffreys-Bullen (J-B) travel-times curves for an earthquake focus at the surface. Travel-times are represented as a function of epicentral distance. Each line represents a given travel path through the Earth's interior (after[START_REF] Bullen | An introduction to the theory of seismology[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . 2 Schematic of an onshore and offshore setup in seismic exploration surveys (Energy Information Australia). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Refraction/reflection FWI results at exploration scale. (a) correspond to the initial model used to start FWI workflow. (b) final 25 Hz FWI results. From Dickinson et al. (2017). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Small scale FWI results at 65 Hz. Here the target is the Ettlingen trench near Karlsruhe, Germany. The low-velocity anomaly in the center corresponds to fill-in materials used to level the ancient trench. From Irnaka, 2021. . . . . . .5Example of global FWI for vertically polarized shear-wave-speed perturbations in GLAD-M25 at 250 km depth[START_REF] Lei | Global adjoint tomography-model glad-m25[END_REF]. Remarkable hotspots are highlighted with white arrow. . . . . . . . . . . . . . . . . . . . . . . . . . . . .1.1 Shematic illustration of the accuracy and resolution of tomography, fullwaveform inversion and migration methods. The solid black line corresponds to the resolution gap between tomography and migration technics (as proposed by[START_REF] Claerbout | Imaging the Earth's interior[END_REF]. This gap is supposed to be filled by full-waveform inversion (figure from[START_REF] Yao | A review on reflection-waveform inversion[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Couple of source/receiver of a fixed diffraction point in a homogeneous medium. 1.3 Representation of misfit function f [m] : Ω ⊂ R 2 → R presenting three local minima, two of which are global and denoted by a red dot and a third one denoted by a blue dot. From Thurin (2020). . . . . . . . . . . . . . . . . . . . . 1.4 Schematic example of the cycle-skipping/phase-ambiguity issue on sinusoidal signals.As soon as the initial shift is larger than half a period, the signal's fit using a least-squares distance is performed up to one or several phase shifts. One may try to fit the n + 1 dashed wriggle of the top signal with the n continuous wriggle of the middle signal moving in the wrong direction. The bottom dashed signal predicts the n wriggle in less than a half-period leading to a correct updating direction (from[START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF]. . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 1.5 Simple time-shifted Ricker wavelet (4 Hz) test for L 2 misfit function. (a) presentation of the setup, with reference signal (seen as d obs ) in solid black, and shifted signal in dashed black (seen as d cal ). (b) Value of the L 2 norm with respect to time-shift. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 FWI schematic example for cycle skipping: (a) true model with spherical inclusion, (b) homogeneous initial model 1, (c) homogeneous initial model 2. . . . . 1.7 FWI schematic example for cycle skipping: traces for a source-receiver couple at z = 500 m for d obs (solid black), d cal in the initial model (dashed black) and d cal in final reconstructed V P model (solid red). (a) correspond to initial model 1, (b) to initial model 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 FWI schematic example for cycle skipping: reconstructed V P models using FWI starting from (a) homogeneous intial model 1, (b) homogeneous initial model 2. 1.9 Same 1D setup as Figure 1.5, but this time comparing the influence of frequency on the convexity. (a) 0.5 Hz Ricker wavelet, (b) 2 Hz Ricker wavelet, (c) 4 Hz Ricker wavelet. (d) L 2 misfit function values with respect to time-shifts for the three different frequency cases. The width of the basin of attraction of the misfit function is enlarged for lower frequencies. . . . . . . . . . . . . . . . . . . . . . 1.10 Illustration of frequency-continuation approach used in FWI; here, on a shallow near-surface target of 30 m (from Irnaka, 2021). FWI is started on a 3 to 25 Hz frequency band, then higher frequency are progressively introduced with an upper limit of 45 Hz then finally 65 Hz. A clear improvement in resolution (smaller wavenumber content) is observed for higher frequency data. . . . . . . 2.1 Comparison of several misfit functions in a simple 1D case for one shifted arrival. The arrival is set to be a Ricker wavelet with a central frequency of 4 Hz. (a) represents the signal used for the test (with only one arrival at the center). The fixed reference signal is displayed in continuous black. The shifted signal is displayed in dotted black (here for +1.5 s). (b) represents the normalized misfit function values with respect to the time-shift (from -1.5 s to 1.5 s). . . . . . . 2.2 (a) quantities Q obs and Q cal for three time-shifts (-1.5 s in green, -0.1 s in blue and "in phase" in dashed red). (b) the area under the curve for Q obs -Q cal for the three time-shifts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 E obs and E cal for three time-shifts (-0.5 s, -0.2 s and in phase). . . . . . . . . 2.4 Detail for ϕ(t) from KROT. (a) the setup with d obs (in black) and d cal for two time-shifts (-0.3 s in red and -0.15 s in blue). (b) shape of L 2 and KROT misfit function with respect to time-shifts, red and blue cross represent the positions of the two time-shifts selected. (c) and (d) respectively display ϕ(t) and d obs -d cal for the two time-shifts of -0.3 s and -0.15 s. (e) the area under the curve for ϕ(t)(d obs -d cal ) quantity for the two time-shifts. This last quantity is used to get the misfit function value after time integration. . . . . . . . . . . . . . . . . 2.5 Same as Figure 2.1 but with two Ricker wavelets with one shifted (left) and one in phase (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 2.15 FWI Test 2: (a) Extracted traces d obs , (b) Wiener filters w(t) and P, and (c) w(t) × P for center shot at zero offset at first iteration. Data associated to negative V P anomaly are in solid blue, and positive V P anomaly in solid red. . 2.16 FWI Test 2: Misfit value for L 2 and AWI with respect to the layer velocity for the two inversion cases: (a) negative and (b) positive V P anomaly. Misfit is calculated between d obs (data in the true medium) and a d cal (V P ) generated in a medium with a correct background V P = 1500 m.s -1 and a as only varying parameter the layer V P , ranging from ±100 m.s -1 around the original velocity of layer in d obs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Test 3.a: (a) True Marmousi II model. (b) S250 initial model, (c) S500 initial model and (d) 1D initial model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Test 3.a: Density model obtained from V P using Gardner's law for: (a) True Marmousi II model, (b) S250 initial model and (c) S500 initial model. . . . . . 2.19 Test 3.a: Common shot gather (CSG) for field data overlapped by synthetic data in (a) S500 initial model, and (b) 1D initial model. Field data in black and white, overlapped by red to blue synthetic data with transparency. Red and white visible mean out of phase, black and blue mean in phase. . . . . . . . . . 2.20 Test 3.a: Inverse crime FWI final reconstructed V P model for Marmousi. Left column corresponds to S500 initial model, right column to 1D initial model. The lines respectively correspond to the final reconstructed V P model using L 2 (a,b), AWI (c,d), IE (e,f), NIM (g,h), KROT (i,j) and GSOT (k,l). . . . . . . . 2.21 Test 3.a: Overlapped common shot gathers for synthetic data in the final reconstructed V P model starting from S500 initial model vs field data. Each subfigure corresponds to misfit function, with L 2 (a), AWI (b), IE (c), NIM (d), KROT (e), and GSOT (f). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.22 Test 3.a: Same as Figure 2.21 but starting from 1D initial model. . . . . . . . . 2.23 Test 3.a: Costs evolution in the inverse crime Marmousi for 1D initial model. (a) evolution of cost functions over iterations, (b) true L 2 cost evolution over iterations, (c) model error reduction over iterations and finaly (d) model error vs. the data error reduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.24 Test 3.b: Comparison of true wavelet used to generate the observed data (for both test 3.a and 3.b) versus the inverted wavelet used for inversion in more realistic test 3.b. Under the associated amplitude spectrum of both true wavelet and estimated wavelets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.25 Test 3.b: More realistic FWI final reconstructed V P model for Marmousi. Left column corresponds to S250 initial model, right column to S500 initial model. The lines respectively correspond to the final reconstructed V P model using L 2 (a,b), AWI (c,d), IE (e,f), KROT (g,h) and GSOT (i,j). . . . . . . . . . . . . . 2.26 Test 3.b: Overlapped common shot gathers for synthetic data in the final reconstructed V P model starting from S250 initial model vs field data. (a) corresponds to the data-fit in the S250 initial model. Then, each subfigure corresponds to misift function, with L 2 (b), AWI (c), IE (d), KROT (e), and GSOT (f). . . . . 2.27 Test 3.b: Same as Figure 2.26 but starting from S500 initial model. . . . . . . LIST OF FIGURES 2.28 Test 3.b: Costs evolution in the more realistic Marmousi for S500 initial model. (a) evolution of cost functions over iterations, (b) true L 2 cost evolution over iterations, (c) model error reduction over iterations and finaly (d) model error vs. the data error reduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Location of the Valhall field on the North Sea (from Thurin, 2020) . . . . . . . 3.2 Layout of the Valhall acquisition overlapped on an horizontal V P slice at 1 km obtained by GSOT-based FWI (from this study). Location of sources (gray dots) and receivers (blue diamonds). Two receivers positions (A and B) are located with black stars. Cables A (x = 2950 m), B (x = 5530 m) and C (x = 3080 m) are identifed. Black dots denote the position of the three V P well sonic-logs. . . 3.3 2D common-receiver gathers extracted for receiver A along cable A for: raw data (top), band 1 data 2.5 -5 Hz (middle) and band 2 data 2.5 -7 Hz (bottom). White dashed arrows point on the Schölte waves which are muted for the inversion. Blue and yellow arrows respectively point on the reflexion from the low velocity zone and the reflexion on the top of the reservoir. Black arrows point on the diving waves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Slices of the initial model V P TOMO . (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Same as Figure 3.4 for initial model V P 1D . . . . . . . . . . . . . . . . . . . . . 3.6 Same as Figure 3.4 for anisotropic parameter η. . . . . . . . . . . . . . . . . . . 3.7 2D common-receiver gathers at 5 Hz. Synthetic data generated into the initial model TOMO are displayed in a blue/white/red color scale, field data are overlapped in grayscale with transparency. The best result is achieved when black and blue are the only colors visible. Red and white are shown when data are not in phase. (a) receiver A along cable A (through the low velocity anomaly). (b) receiver B along cable B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Same as Figure 3.7 into the 1D initial model. . . . . . . . . . . . . . . . . . . . 3.9 Complete FWI workflow used on the OBC Valhall dataset. At the core of the process lies the FWI iterations. Then several loops are nested one into another, from internal FWI iterations, source-subsampling, data selection, to finally the outer one of frequency continuation. . . . . . . . . . . . . . . . . . . . . . . . . 3.10 On top, weighted data for source inversion displayed on a 5 Hz 2D commonreceiver gather (receiver A cable A). On the bottom, the estimated wavelet (left) and associated spectrum (right). . . . . . . . . . . . . . . . . . . . . . . . 3.11 2D common-receiver gathers extracted for receiver A along cable A with different data weighting applied on them. From top to bottom: first break & short offset (DW SO), first break & medium offset (DW MO), first break & full offset (DW FO), first time release & medium offset (RT1 MO), first time release & full offset (RT1 FO), second time release & full offset (RT2 FO). . . . . . . . . . . . . . . LIST OF FIGURES 3.12 Slices of the 7 Hz FWI reconstructed V P using L 2 misfit function starting from TOMO initial model. (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km. 3.13 Same as Figure 3.12 using GSOT misfit function. . . . . . . . . . . . . . . . . . 3.14 Comparison of V P profiles extracted from the TOMO initial model (dashed red), 7 Hz FWI models using L 2 (solid yellow) and 7 Hz FWI model using GSOT (solid purple) with sonic log (solid black). Left subfigure corresponds to the Log 1 at the center of the target. Middle subfigure to the Log 2, and right subfigure to Log 3 (which is far away from the target). . . . . . . . . . . . . . . . . . . . . . 3.15 2D common-receiver gathers at 7 Hz starting from the TOMO initial model. Synthetic data (blue/white/red color scale) generated into the final reconstructed V P using (a,b) L 2 misfit function, (c,d) GSOT misfit function. (a,c) receiver A along cable A (through the low velocity anomaly). (b,d) receiver B along cable B. Field data are overlapped in grayscale with transparency. Black arrows point to improvement obtained with GSOT. . . . . . . . . . . . . . . . . . . . . . . . 3.16 2D common-receiver gathers at 5 Hz starting from the 1D initial model with data weighting apply (DW MO). Synthetic data (blue/white/red color scale) generated into the final reconstructed V P at 2nd workflow step using: (a,b) L 2 misfit function, (c,d) GSOT misfit function. (a,c) receiver A along cable A (through the low velocity anomaly). (b,d) receiver B along cable B. Field data are overlapped in grayscale with transparency. Black arrows point to area where GSOT improves the datafit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Same as Figure 3.16 but with final relaxed data weighting (RT2 FO) for display. This exhibits the improvement of datafit obtained in area which are not yet inverted. Here, black arrows point to improve fit and coherency of the diving wave with GSOT, whereas yellow arrows point to improved data-fit of reflected events with GSOT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Slices of the 5 Hz FWI reconstructed V P using L 2 misfit starting from 1D initial model. (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1.1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km. . . . . . . . . . . . . 3.19 Same as Figure 3.18 using GSOT. Here, results are consistant until 1.5 to 2 km depth compared to L 2 -based FWI. Characteristic structures of the Valhall field are recovered. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 2D common-receiver gathers at 7 Hz starting from the 1D initial model. Synthetic data (blue/white/red color scale) generated into the final reconstructed V P using GSOT. (a) receiver A along cable A (through the low velocity anomaly). (b) receiver B along cable B. Field data are overlapped in grayscale with transparency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Slices of the 7 Hz FWI reconstructed V P using GSOT misfit starting from 1D initial model. (a-c) Horizontal slices at (a) 200 m depth, (b) 500 m depth and (c) 1.1 km depth. (d-e) Inline vertical slices for (d) x = 2.95 km and (e) x = 3.95 km. (f-g) Cross-line vertical slices at (f) y = 9 km and (g) y = 6 km. . . . . . . . . .

CONTENTS 3

 3 Robust FWI with graph space optimal transport: application to 3D OBC Valhall data 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Field data presentation and full-waveform inversion workflow . . . . . . . . . . 3.3.1 Geological situation, dataset and initial models . . . . . . . . . . . . . . 3.3.2 Full-waveform inversion workflow . . . . . . . . . . . . . . . . . . . . . . 3.4 Full-waveform inversion results . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Starting from the reflection tomography model . . . . . . . . . . . . . . 3.4.2 Starting from the 1D initial model . . . . . . . . . . . . . . . . . . . . . 3.4.3 Computational costs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Exploring new ideas and overcoming practical limitations 4.1 Overcoming practical issue for 3D KROT applications . . . . . . . . . . . . . . 4.2 Exploring new ideas of alternative misfit functions . . . . . . . . . . . . . . . . Investigation published during my thesis . . . . . . . . . . . . . . . . . . . . . .

	CONTENTS
	Conclusion References
	x	ix

4.2.1 Instantaneous phase and unwrapped phase misfit functions . . . . . . . 4.2.2 Spline based envelope misfit function . . . . . . . . . . . . . . . . . . . . 4.2.3 Diffusion based misfit function . . . . . . . . . . . . . . . . . . . . . . . 4.2.4 Beat-tone misfit function . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 and perspectives Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

	Chapter 2
	On cycle-skipping and misfit
	functions modification for full-wave
	inversion: comparison of five recent
	approaches

Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 General FWI framework and misfit function formulation . . . . . 35 2.2.1 General framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.2.2 Adaptive waveform inversion . . . . . . . . . . . . . . . . . . . . . . . 36 2.2.3 Instantaneous envelope . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.2.4 Normalized Integration Method . . . . . . . . . . . . . . . . . . . . . . 38 2.2.5 Kantorovich-Rubinstein optimal transport distance . . . . . . . . . . . 39 2.2.6 Graph space optimal transport distance . . . . . . . . . . . . . . . . . 40 2.3 A simple convexity analysis based on time-shifted Ricker wavelets 41 2.4 FWI tests on two canonical examples . . . . . . . . . . . . . . . . . 46 2.4.1 FWI Test 1: transmission configuration . . . . . . . . . . . . . . . . . 46 2.4.2 FWI Test 2: reflection configuration . . . . . . . . . . . . . . . . . . . 49 2.5 Marmousi case study: toward a more realistic case study . . . . . 58 2.5.1 Common framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 2.5.2 Inverse crime inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 2.5.3 A more realistic inversion . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.6

•

  First, my Master's thesis, focusing on improving FWI robustness with alternative misfit functions:[START_REF] Pladys | Sensibility of several misfit functions for Full Waveform Inversion: theoretical analysis, implementation and comparative study on multiple realistic cases[END_REF]. Sensibility of several misfit functions for Full Waveform Inversion: theoretical analysis, implementation and comparative study on multiple realistic cases. Master's thesis, Univ. Grenoble Alpes.• Some work on key properties that alternative misfit functions should satisfy: Pladys, A.,Brossier, R., and Métivier, L. (2017). FWI alternative misfit functions: what properties should they satisfy. In 79th EAGE Conference and Exhibition 2017, Paris, page Tu P1 01.

• Case study for

GSOT-based FWI on Valhall: Pladys, A., Brossier, R., Irnaka, M., Kamath, N., and Métivier, L. (2019)

. Assessment of optimal transport based FWI: 3d OBC valhall case study. In SEG Technical Program Expanded Abstracts 2019, pages 1295-1299.

Remerciements Acknowledgments

This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.f r), sponsored by AKERBP, CGG, CHEVRON, EQUINOR, EXXON-MOBIL, JGI, SHELL, SINOPEC, SISPROBE and TOTAL. This study was granted access to the HPC resources of the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07 13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche, and the HPC resources of CINES/IDRIS/TGCC under the allocation 046091 made by GENCI.

I thank AKERBP ASA and their partner Pandion Energy for providing the dataset and permission to present this work. Special thanks to Einar Jarle Kjos, Ross Milne and Espen Birger Raknes for their help and support on the Valhall data. I thank Denes Vigh and Xin Cheng from Schlumberger for their support and warm welcome, for their in-depth and interesting scientific discussion, and for providing me with the opportunities to come to their office to present my work.

I thank Equinor for the opportunities to work on 3D OBC real data applications in their office for two summer internships. I thank Lyubov Skopintseva and Joachim Mispel for their help and support. I sincerely thank Céline Ravaut for the opportunities you provided to me, for the help you provided when I started working for the first time on field data, and for the interesting scientific discussion.

Contents

Field data presentation and full-waveform inversion workflow