Keywords: Waveform inversion, wavefield decomposition, seismic imaging

2-13 Principle of cycle-skipping problem [Virieux and Operto, 2009]. If the phase mismatch is less than one half of the period, the local optimization methods could adjust the correct phase. In the contrary case, the local optimization methods adjust these two signals with one period shift and the inversion falls into a local minimum. [Virieux and Operto, 2009] . . . . . .

I would like to thank all the people who have contributed to this thesis work.

First of all, I would like to thank my advisers Hervé Chauris and Daniela Donno. I thank them for their availability, patience, support and encouragement given along my stay in the geophysics group. I thank them also for having trusted in my ability to develop this research even though I had not much experience in geophysics before the thesis. I especially acknowledgement their patience with me and explain to me even the basic knowledge and formulations. This really helps me to start the thesis in the irst year. The discussions with them have always been helpful for me. Thanks for their availability even on weekends to answer my questions or to help me revising reports and presentations. Their passion and persistence for the research also encourage me.

I would like to thank the members of my jury. Thanks to Paul Sava and Gilles Lambaré for having accepted to review my thesis in detail. Thanks to Stéphane Operto and Romain Brossier for examining my thesis.

I thank Total for inancially and technically supporting my thesis. I especially acknowledge Henri Calandra from Total who has brought me to the area of geophysics and initiated this PhD project. I thank François Audebert from Total for providing the real data set used for the real application in this thesis and for his help to apply for the computational resources from Total. I thank Florent Ventimiglia from Total for kindly helping me with troublesome computer and account problems.

Many thanks to the members in the geophysics team of the Geosciences Center at MINES ParisTech. I especially thank Véronique Lachasse for her help to facilitate the administrative afairs since my irst day at the oice. Thanks to the Alexandrine Gesret and Mark Noble for their availability and advices.

I would like to thank all my oicemates during the three and half years. Nidhal Belayouni, Sébastien Penz, Carlos Pérez, Elise Vi Nu Bas, Charles-Antoine Lameloise, Emmanuel Cocher, Sven Schilke, Yves-Marie Batany, Yubing Li and Hao Jiang for their friendship and support. Thanks Elise for ofering me facilities when I arrived in Fontainebleau.

I also thank Dominique Vassiliadis from the Geosciences Center and Evelyne Peyrillou from the delegation to help me solving administrative matters regarding the visa, contract, residence cart, etc.

Finally, I would like to thank my family for always respecting and supporting my professional and life choices. My research would not have been possible without their support. Before continuing with the text of this research, I would like to especially thank my boyfriend Haisheng Wang who voluntarily came to France to accompany me and has always been at my side and supporting me during these years.

3

Abstract

The knowledge of the Earth subsurface poses economic, environmental, human and scientiic issues. Seismic imaging is a procedure to image the Earth subsurface from the data observed at the surface. In the context of hydrocarbon exploration, seismic imaging techniques are widely used to characterise the irst few kilometres of the Earth's interior. Full Waveform Inversion (FWI) is one of the eicient seismic imaging method. Recent advances in high performance computer make FWI feasible for large applications. In theory, FWI could reconstruct a high-resolution subsurface image provided that low frequency and wide angle/aperture/azimuth data are available. FWI is a data-itting procedure and is resolved as an optimization problem. Depending on the frequency content of the data, the objective function of FWI may be highly nonlinear and has many local minima. If a data set mainly contains relections, this problem particularly prevents the gradient-based methods from recovering the long wavelengths of the velocity model.

The model of the subsurface could be separated into two parts by scale separation. One part is the short-wavelength part which contains the singularities of the model and the other part is the long-wavelength part which is a smooth version of the model. In this thesis, I propose a variant of FWI based on the scale separation of these two parts to mitigate the nonlinearity of the problem.

In the irst section, methodologies of the conventional FWI and the new proposed method are presented. The new method is a Relection-based Waveform Inversion (RWI) method . It consists of decomposing the gradient of FWI into a short-wavelength part and a long-wavelength part and the inversion is performed in an alternating fashion between these two parts. The gradient decomposition is achieved by decomposing the waveields into their one-way components. Diferent waveield decomposition methods are also presented.

In the second section, we implement the FWI and the new method to several case studies. For numerical modeling, we use a inite-diference approach to resolve the acoustic wave equation with constant density in the time domain. The model update is based on the L-bfgs algorithm and the waveield is decomposed using the 2D FFT-based method in the -domain. These case studies show the diiculties associated with FWI to recover the long-wavelength part of the velocity model when low frequency and large-ofset data are absent, and the initial model is far from the true one. The new method shows its robustness in this case especially for constructing the long-wavelength model. Keywords: Waveform inversion, waveield decomposition, seismic imaging. 5 Résumé La connaissance du sous-sol de la Terre pose les enjeux économiques, environnementaux, humains et scientiiques. L'imagerie sismique est une procédure à imager le sous-sol à partir des données observées à la surface. Dans le cadre de l'exploration des hydrocarbures, les méthodes d'imagerie sismique sont beaucoup utilisées pour caractériser les premiers kilomètres de l'intérieur de la Terre. Full Waveform Inversion (FWI) est l'une des méthodes d'imagerie sismique eicace. Les progrès récents en superordinateur rendent FWI possible pour les grandes applications. En théorie, FWI pourrait reconstruire une image du sous-sol à haute résolution, à condition que les données de basse fréquences et grandangle/ouverture/azimut sont disponibles. FWI est une procédure d'ajuster des données, et est résolu comme un problème d'optimisation. En fonction du contenu en fréquence des données, la fonction objectif de FWI peut être fortement non linéaire et présente de nombreux minima locaux. Pour des données qui contiennent principalement des rélexions, ce problème empêche notamment les méthodes basées sur le gradient de retrouver les longues longueurs d'onde du modèle de vitesse.

Le modèle du sous-sol peut être séparé en deux parties par separation d'échelle. Une partie est la partie courte longueur d'onde qui contient les singularités du modèle. l'autre partie est la partie longue longueur d'onde qui est une version lisse du modèle. Dans cette thèse, nous proposons une variante de FWI basée sur la séparation d'échelle entre les deux parties pour atténuer la non-linéarité du problème.

Dans la première section, les méthodologies de la FWI classique et la nouvelle méthode proposée sont présentés. La nouvelle méthode est une méthode de l'inversion des formes d'ondes basées sur les rélexions. Il consiste à décomposer le gradient de FWI en une partie de courte longueur d'onde et une partie de longue longueur d'onde, et l'inversion est efectuée d'une manière alternée entre ces deux parties. La décomposition du gradient est obtenu par la décomposition des champs d'ondes en leurs composants unidirectionnels. Diférentes méthodes de décomposition des champs d'ondes sont également présentées.

Dans la deuxième section, nous appliquons la FWI et la nouvelle méthode à plusieurs études de cas. Pour la modélisation numérique, nous utilisons une approche diférence inie pour résoudre l'équation des ondes acoustiques avec une densité constante dans le domaine temporel. La mise à jour du modèle est basé sur l'algorithme L-BFGS et les champs d'ondes sont décomposés en utilisant la méthode basée sur FFT 2D dans le domaine -. Ces études montrent les diicultés liées à FWI pour récupérer la partie longue longueur d'onde du modèle de vitesse lorsque les données de basse fréquence et grands ofsets sont absentes, et le modèle initial est loin du vrai modèle. La nouvelle méthode présente sa robustesse dans ce cas en particulier pour la construction du modèle de longue longueur d'onde.

List of igures

1-1 Wavenumber at a scattering point [START_REF] Huang | Resolution limits for wave equation imaging[END_REF]. . . . . . . . 1-2 Separation of scales between the frequency content of the velocity model and of the seismic data [Claerbout, 1982] [START_REF] Huang | Resolution limits for wave equation imaging[END_REF]. . [2][3][4][5][6] Wavepath for the irst-order free-surface multiple [START_REF] Huang | Resolution limits for wave equation imaging[END_REF] [3][4][5][6][7][8][9][10][11][12] Gradients in equation (3.11). Correlation of (left) 0 and , (right) and 0 . To be compared with Figure 3-8 (last two panels). . . . . . . . . . [3][4][5][6][7][8][9][10][11][12][13] One snapshot of the waveield propagated in the model in Figure 2 [Chauris et al., 2006]. . . . . . . . . . . . . . . 3-25 A direction example of curvelet [Chauris et al., 2006]. . . . . . . . . . . . . The direction (left) and scale (right) ilter in the curvelet domain [START_REF] Chauris | Seismic demigration/migration in the curvelet domain[END_REF] [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] The diference between the FWI model (Figure 5678) and the initial model (Figure 5 L'imagerie sismique permet d'imager l'intérieur de la Terre, pour analyser les propriétés physiques du sous-sol. L'imagerie sismique pose des questions économiques, humaines, environnementales et scientiiques. Dans la géophysique d'exploration, cette technique est principalement utilisée pour la détection d'hydrocarbures. Avec le développement des ordinateurs de haute performance, l'imagerie sismique 3D est de plus en plus sophistiquée [Biondi, 2006], pour obtenir de meilleures résolutions. Le résultat de l'imagerie sismique 3D est similaire à une imagerie médicale. Mais l'imagerie sismique se prolonge plusieurs kilomètres ou plus dans la Terre. Et surtout, les sources et récepteurs ne sont qu'à la surface. Les méthodes d'imagerie sismique continuent à s'améliorer, et les ordinateurs les plus avancés maintenant permettent aux scientiiques de traiter les données en quelques jours plutôt qu'en mois, et accélèrent la découverte et la production inale du pétrole et du gaz.

Dans l'imagerie sismique, la propriété physique nous voulons imager est généralement le modèle de vitesse de la propagation des ondes, en utilisant les données sismiques enregistrées à la surface pendant l'acquisition terrestre ou marine. Les données enregistrées apportent des informations sur le sous-sol. La loi physique reliant le modèle de vitesse et les données est décrite par l'équation de propagation des ondes. Cette relation est habituellement non linéaire, et l'imagerie sismique est considérée comme un problème d'optimisation non linéaire.

Le modèle de vitesse que nous cherchons à caractériser pourrait être séparé en deux parties. Une première partie est la partie des grandes longueurs d'onde qui peut être considérée comme une version lisse du modèle, et qui est désignée comme le macro-modèle. L'autre partie est la partie des courtes longueurs d'onde qui contient toutes les singularités du modèle, et est appelée le modèle de rélectivité ou l'image migrée. La procédure standard dans l'imagerie sismique consiste à reconstruire les grandes longueurs d'onde en premier, suivie par la reconstruction des courtes longueurs d'onde. La qualité de la mise à jour des courtes longueurs d'onde dépend de l'exactitude du macro modèle.

Avec le développement des ordinateurs de haute performance, l'imagerie sismique 3D et 4D s'est développée. Elle fournit des images du sous-sol de haute résolution. Néanmoins, l'imagerie sismique présente plusieurs limitations communes. Par exemple, la résolution de l'imagerie sismique est limitée par la bande de fréquence de la source, de l'atténuation et de la géométrie d'acquisition. A cause de la bande limitée de la source et de l'acquisition inie, souvent il est diicile d'obtenir une résolution suisante. De plus, les basses fréquences dans la source et les données de grands ofsets et grandes ouvertures sont importantes pour la reconstruction des grandes longueurs d'ondes du modèle. Cependant, les données de grands ofsets et grandes ouvertures présentent des déis de déploiement et inanciers.

L'inversion des forme d'ondes (en anglais, Full Wavefrom Inversion, FWI) est l'une des techniques principales de l'imagerie sismique. Nous pouvons citer [Fichtner, 2010] et [Virieux and Operto, 2009] pour une revue tutorielle de cette méthode. FWI déinit un problème inverse non linéaire dans le domaine des données [Tarantola, 1984, Lailly, 1983], cherchant à minimiser les moindres carrés des diférences entre les données observées et simulées. Le principe de FWI est d'utiliser tous les types d'ondes (ondes directes, rélexions, réfractions, et multiples) pour résoudre les diférents paramètres du modèle (la vitesse, la densité, l'atténuation, . . . ). Les informations du champ d'ondes complet (phase et l'amplitude) sont tous pris en considération en même temps, sans avoir besoin d'introduire explicitement le temps de trajet. FWI est capable de traiter des modèles complexes et de fournir des images de haute résolution. Dans l'industrie, FWI est principalement utilisée pour l'exploration pétrolière. Toutefois, elle est également appliquée à d'autres domaines, tels que l'électromagnétisme.

Les méthodes d'optimisation utilisées pour FWI pourraient être divisées en méthodes d'optimisation globales et méthodes d'optimisation locales. Les méthodes d'optimisation globales, telles que Monte-Carlo [Jin and Madariaga, 1994], les algorithmes génétique [Jin and Madariaga, 1993], permettent potentiellement d'explorer l'ensemble de l'espace pour trouver le minimum global de la fonction objective. Les méthodes d'optimisation locales, telles que les méthodes basées sur le gradient [Tarantola, 1984], commencent à partir d'un modèle initial, et calculent la direction de la mise à jour du modèle en utilisant le gradient de la fonction objective. Le gradient peut être calculé en utilisant la méthode de l'étatadjoint [Plessix, 2006], puis le modèle est mis à jour et des itérations sont répétées jusqu'à ce que la critère de convergence est remplie. Bien que les méthodes globales soient plus robustes pour résoudre les problèmes non linéaires, elles demandent un coût de calcul élevé, et la taille de l'espace des modèles et la taille des données à traiter sont au-delà de la capacité de calcul actuelle pour les applications réelles. Généralement, une méthode d'optimisation locale est utilisée pour la FWI et le problème d'optimisation est résolu de manière itérative pour compenser la non-linéarité [Tarantola, 1984].

FWI a été initialement appliquée en temps [Tarantola, 1984, Lailly, 1983]. La première application 2D de FWI en géophysique d'exploration a été mise en oeuvre par [START_REF] Gauthier | Twodimensional nonlinear inversion of seismic waveforms[END_REF]]. [Tarantola, 1986] and [Mora, 1987] ont appliqué FWI en temps au cas élastique. FWI appliquée sur des données réelles est mis en oeuvre par [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF]. FWI en fréquence a été proposée depuis les années 90 [Pratt andWorthington, 1990, Pratt, 1999]. FWI multi-échelle a été étudiée dans le domaine temporel [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF] et dans le domaine fréquentiel [START_REF] Sirgue | Eicient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF] en vue d'atténuer la non-linéarité du problème inverse. Cette approche consiste à d'abord inverser les données de plus basses fréquences et d'introduire progressivement des fréquences élevées dans l'inversion. Avec le développement de l'ordinateur haute performance, FWI est entré dans une nouvelle ère depuis 2009, lorsque les premières applications 3D de FWI ont été obtenus avec des données réelles [Plessix, 2009, Sirgue et al., 2010].

Un problème inverse est considéré comme le pendant du problème direct qui relie les paramètres de modèle avec les données que nous observons. Le problème direct de FWI simule la propagation des ondes dans un modèle donné. Les données réelles sont décrites plus précisément par la modélisation élastique plutôt que la modélisation acoustique, où la vitesse des ondes P, la vitesse des ondes S, la densité et l'atténuation sont tous considérées. Toutefois, il est encore une pratique courante d'appliquer la FWI acoustique sur des données réelles en considérant la propagation des ondes acoustiques. L'équation des ondes acoustiques pourrait être résolue dans le domaine temporel [Tarantola, 1984, Mora, 1987, Mora, 1989] ou dans le domaine fréquentiel [Pratt andWorthington, 1990, Sirgue andPratt, 2004]. La résolution numérique de l'équation diférentielle partielle pourrait être réalisée par des méthodes numériques, telles que la méthode des diférences inies [Virieux, 1986, Moczo et al., 2004], éléments inis [Marfurt, 1984], etc. Lorsque la discrétisation spatiale est efectuée sur une grille régulière, la méthode des diférences inies est généralement choisie car elle est rapide et facile à mettre en oeuvre. FWI est capable de potentiellement imager des structures complexes. Cependant, elle présente plusieurs limitations [Virieux and Operto, 2009]. Pour atténuer la nonlinéarité de FWI, de nombreuses formulations alternatives sont proposées, qui sont principalement liées aux choix de la représentation des données, aux choix de la norme des résidus des données, et aux choix du critère d'optimisation.

Lorsqu'on traite les données de sismique rélexion, l'utilisation des méthodes d'inversion basées sur les rélexion (en anglais Relection Waveform Inversion, RWI) [Mora, 1989, Chavent et al., 1994, Xu et al., 2012, Biondi et al., 2012, Tang and Lee, 2013, Wang et al., 2013, Brossier et al., 2015, Zhou et al., 2015[START_REF] Wang | Inversion of seismic refraction and relection data for building long-wavelength velocity models[END_REF] semblent prometteuses. Ces approches reposent sur la séparation explicite du modèle en un modèle de grandes longueurs d'onde et un modèle de rélectivité. Le modèle de rélectivité est obtenu en premier par la migration et ensuite se comporte comme les sources primaires en profondeur, ce qui permet de produire les transmissions entre les rélecteurs et la surface. Ces transmissions sont utiles pour la mise à jour du modèle de grandes longueurs d'ondes.

Des données de rélexions contiennent des informations de temps de trajet ainsi que l'amplitude, donc l'inversion de la forme d'onde des rélexions devrait en principe être capable de reconstruire à la fois les grandes longueurs d'ondes et les courtes longueurs d'ondes [Snieder et al., 1989, Hicks andPratt, 2001]. Toutefois, en pratique, la qualité de l'inversion est dégradée à cause des données à bande limitée et ofsets limités. Dans le domaine temporel, lorsqu'une grande gamme de fréquences est utilisée en même temps, ces deux composantes sont mélangées pendant la FWI. Ce couplage est un problème majeur pour la FWI [START_REF] Snieder | Retrieving both the impedance contrast and background velocity: A global strategy for the seismic relection problem[END_REF], puisque la mise à jour des grandes longueurs d'ondes a généralement une amplitude plus faible que la mise à jour des courtes longueurs d'onde. Par conséquent, lorsque ces deux composants sont mélangées ensemble, la mise à jour du modèle est dominée par les courtes longueurs d'ondes. Il est donc naturel d'inverser séparément ces deux parties.

Dans cette thèse, nous proposons une méthode d'inversion en deux étapes obtenue par la décomposition de la formule de gradient de FWI en une partie de grandes longueurs d'ondes et une partie de courtes longueurs d'ondes [START_REF] Wang | Taking advantage of waveield decomposition in full waveform inversion[END_REF].

Dans le chapitre 2, nous passons en revue les notation basiques de la FWI, y compris la fonction objective et la résolution du problème direct en utilisant le schéma de diférences inies. Le calcul du gradient de la fonction objective et les méthodes d'optimisation locales sont également détaillées. La résolution du gradient est analysée pour étudier la résolution de FWI.

Dans le chapitre 3, nous présentons la méthode d'inversion proposée, et montrons comment les deux étapes d'inversion séparés, soit une pour la mise à jour des courtes longueurs d'ondes et l'autre pour la mise à jour des grandes longueurs d'ondes, sont formulés. Nous illustrons cette approche en utilisant un modèle de couche 1D. Diférentes méthodes de décomposition des champs d'ondes sont également introduites, l'une dans le domaine de Fourier, l'une en utilisant le vecteur de Poynting et la troisième dans le domaine des curvelets.

Dans le chapitre Seismic imaging allows to image the Earth's interior, to analyze the physical properties of the subsurface. It poses economic, human, environmental and scientiic issues. In exploration geophysics, seismic imaging is mainly used for the detection of hydrocarbons. With the development of the high performance computers, sophisticated 3D seismic imaging is investigated [Biondi, 2006], yielding high-resolution images of the subsurface. The result of 3D seismic imaging is similar to an X-ray scan or medical imaging that extends several kilometres or more into the Earth. 4D seismic imaging [START_REF] Zhang | Understanding reservoir architectures and steam-chamber growth at Christina Lake, Alberta, by using 4D seismic and crosswell seismic imaging[END_REF] allows imaging potential luid migration paths within the reservoir, by inverting for time-lapse parameters. Methods for seismic imaging continues to improve, and more advanced computers now enable scientists to accelerate the data processing, speeding up the discovery and inal production of oil and gas.

In seismic imaging, the physical property we want to image is usually the model of the wave propagation velocity, using the seismic data recorded during land or marine acquisitions. Seismic acquisitions use controlled seismic sources (explosive and vibroseis for land acquisitions or air guns for marine acquisitions). Then, seismic waves propagate in the subsurface and relected or refracted waves are recorded by the seismic sensors placed at the surface (geophones for land acquisitions or hydrophones for marine acquisitions). The recorded data bring information of the subsurface. The physical law relating the velocity model and the data is described by the wave propagation equation. This relationship is usually nonlinear, and seismic imaging is considered as a nonlinear optimization problem.

1.1.2

Review of classic methods/scale separation

The wave-propagation velocity model that we would like to characterize could be divided into two parts. One part is the long-wavelength part which can be considered as a smooth version of the model, and is referred to as the macro or background model. The other part is the short-wavelength part which contains all the singularities of the model, and is referred to as relectivity model or migrated image. The standard procedure in seismic imaging consists of recovering the long wavelengths irst, followed by the reconstruction of the short wavelengths, as the quality of the short-wavelength update depends on the accuracy of the background model. The background model could be recovered either by traveltime tomography [START_REF] Bishop | Tomographic determination of velocity and depth in laterally varying media[END_REF], Pratt and Chapman, 1992, Billette and Lambaré, 1998, Woodward et al., 2008] or by Migration Velocity Analysis (MVA) [START_REF] Symes | Velocity inversion by diferential semblance optimization[END_REF], Chauris and Noble, 2001, Sava and Biondi, 2004]. Traveltime tomography uses the traveltimes to estimate the wave propagation velocity. In industry, the ray-based tomography approaches involved migration, residual move-out picking, demigration and linear or non-linear velocity updates. Nonlinear slope tomography uses a local focusing criterion without using a pre-deined shape of the relectors or of the RMO curves [START_REF] Guillaume | Kinematic invariants: an eicient and lexible approach for velocity model building[END_REF].Migration Velocity Analysis (MVA) is performed in the image domain and aims at building a relationship between migrated image perturbation and model perturbation. This approach a priori does not require picking.

The short wavelengths provide the ine structure of the subsurface model, allowing to localize relectors in depth. The short wavelengths can be obtained by migration techniques, assuming the background velocity is correct. There are two major categories of migration methods: ray-based methods [Beylkin, 1985, Bleistein, 1987], which are based on the high frequency assumption, and wave-equation based methods [START_REF] Baysal | Reverse time migration[END_REF], Whitmore et al., 1983]. [START_REF] Etgen | An overview of depth imaging in exploration geophysics[END_REF] gives a comparison of diferent migration methods. The principle of migration is formulated by [Claerbout, 1971], and it consists of propagating the source signal and the recorded data into to the medium and cross correlate these two waveields. The zero-lag cross correlation gives the locations of relectors. There exist other imaging conditions, such as deconvolution-based imaging condition [START_REF] Valenciano | 2-D deconvolution imaging condition for shot-proile migration[END_REF], source/receiver-normalized imaging condition [START_REF] Kaelin | Imaging condition for reverse time migration[END_REF], extended imaging condition [START_REF] Sava | Time-shift imaging condition in seismic migration[END_REF]. [START_REF] Chattopadhyay | Imaging conditions for prestack reverse-time migration[END_REF] and [START_REF] Sava | Overview and classiication of waveield seismic imaging methods[END_REF] give a summary of the imaging conditions. The classical correlation-based migration is qualitative, as it only provides a relectivity image. Alternately, quantitative migration [START_REF] Lambaré | Iterative asymptotic inversion in the acoustic approximation[END_REF], Jin et al., 1992, Lameloise et al., 2015, Symes, 2015] allows imaging the values of the physical parameters. Recent developments have shown that for migration-based velocity analysis, quantitative migration is preferable as it provides a more accurate migration image.

Diiculties in seismic imaging

With the development of advanced computers, 3D and 4D seismic imaging are investigated, yielding more accurate subsurface images. Nevertheless, seismic imaging presents several common limitations:

• Insuicient resolution. The quality of seismic imaging resolution is limited by the source wavelet bandwidth, the wave propagation attenuation and the acquisition geometry. In the framework of the single-scattering approximation [Devaney, 1982, Miller et al., 1987], the wavenumber at a point of the model depends on the local wavelength and the scattering angle (Figure 1-1).

According to Figure 1-1, high resolution can be obtained with small relection angle and small wavelength, and small wavelength corresponds to small velocity and high frequency. However, during the wave propagation, waves sufer from attenuation efect. Apart from the geometric spreading, intrinsic attenuation is also presented due to the nature of some material. High frequency waves are particularly attenuated during the propagation, making diicult to image the deep part of the Earth.

Besides, as we will show in Chapter 2, numerical resolution of the wave equation imposes the stability condition on the spatial sampling, which depends on the maximum frequency of the source wavelet. The use of a high frequency source requires iner spatial and temporal sampling, making the problem more computational extensive.

• Acquisition limits. The quality of the migrated image depends closely on the acquisition geometry and acquisition devices. [Mora, 1989] and [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF] show that the long wavelengths of the model could be retrieved through the use of long-ofset (diving waves, refractions) and transmission data. If the data lacks of low frequencies, then long-ofset and wide-aperture/azimuth acquisitions [Sirgue Figure 1-1: Wavenumber at a scattering point [START_REF] Huang | Resolution limits for wave equation imaging[END_REF]]. et al., 2010, Shipp and Singh, 2002] are necessary in order to image deep parts of the model [Sirgue, 2006]. However, long-ofset and wide-aperture/azimuth acquisitions usually present deployment and inancial challenges.

• Requirement of good velocity model. [Versteeg, 1993] analysed the sensitivity of the depth migrated image to the velocity model and found that the accuracy of the depth image closely depends on the accuracy of the long wavelengths of the model, and, if an incorrect model is used, the migration image is unfocused and relectors are mispositioned. In practice, the inversion of the background velocity model still remains a challenge.

Full waveform inversion

We now concentrate on a very popular method used for imaging the Earth. We refer to [Virieux and Operto, 2009] and [Fichtner, 2010] for recent reviews.

History

Full waveform inversion (FWI) is one technique for seismic imaging which develops rapidly. FWI deines a nonlinear inverse problem in the data-domain [Tarantola, 1984, Lailly, 1983], seeking to minimize the least-squares diferences between observed and simulated data.

The principle of FWI is to use all types of waves (direct waves, relections, refractions, and multiples) to resolve diferent model parameters (velocity, density, attenuation). The information of the complete waveield: traveltime, phase and amplitude are all taken into consideration at once. FWI is capable of dealing with complex models and of delivering high resolution images with a resolution of half the minimal wavelength [START_REF] Sirgue | Eicient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF]. In industry, FWI is mainly used for the oil exploration; however, it could also be applied to other domains, such as electromagnetics and medical imaging.

The optimization methods used to solve the FWI nonlinear problem could be divided into global optimization methods and local optimization methods. Global optimization methods, such as Monte-Carlo [Jin and Madariaga, 1994], genetic [Jin and Madariaga, 1993] and simulated annealing algorithms [START_REF] Sen | Nonlinear one-dimensional seismic waveform inversion using simulated annealing[END_REF], allow to potentially explore the whole model space to ind the global minimum of the objective function. Local optimization methods, such as gradient-based methods [Tarantola, 1984], start from an initial model, and calculate the model update direction by computing the gradient of the objective function. The gradient could be calculated using the adjoint-state method [Plessix, 2006] and then the model is updated and iterations will be repeated until the convergence criterion is met. Although global methods are more robust to handle nonlinear problems, they require a high computational cost, and the size of the model space and the size of data to process is beyond current computation capacity for real data applications. Usually a local optimization method is used and the optimization problem is resolved iteratively to compensate the nonlinearity [Tarantola, 1984].

FWI was originally applied in the time domain [Tarantola, 1984, Lailly, 1983]. The irst 2D application of FWI in exploration geophysics was implemented by [START_REF] Gauthier | Twodimensional nonlinear inversion of seismic waveforms[END_REF]. [Tarantola, 1986] and [Mora, 1987] applied the time-domain FWI to the elastic case. FWI applied on real data set is implemented by [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF]. The frequencydomain FWI was proposed since the '90s [Pratt andWorthington, 1990, Pratt, 1999]. The frequency-domain FWI is equivalent to the time-domain FWI when all frequencies are considered simultaneously [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF]]. From a series of discrete frequencies, FWI in the frequency domain aims at retrieving model perturbations [Sirgue, 2006]. Rules have been designed to select frequencies [Sirgue, 2006]. Multiscale FWI was investigated in the time domain [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF] and in the frequency domain [START_REF] Sirgue | Eicient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF] in order to mitigate the nonlinearity of the inverse problem. The method consists of starting by inverting the lower frequency data and of progressively introducing high frequencies in the inversion. Moreover, long-ofset data could be used to retrieve irst the long wavelengths of the model [Shipp andSingh, 2002, Sirgue, 2006].

With the development of high performance computing, FWI entered in a new era since 2009, when the irst 3D FWI applications were successfully obtained with real data [Plessix, 2009, Sirgue et al., 2010]. The redundancy of the 3D acquisitions, together with the availability of long-ofset and wide-azimuth data allow reducing the non linearity of the problem.

Physics of wave propagation

An inverse problem is considered as the "inverse" of the forward problem which relates the model parameters to the data we observe. The forward problem in FWI simulates the wave propagation in a given model. Among the recorded seismic waves, we can distinguish between P waves and S waves. P waves, also called pressure waves or primary waves, are characterized by particle motion in the same direction as wave propagation. S waves, also called shear waves or secondary waves, are characterized by particle motion in the plane perpendicular to the direction of wave propagation. P waves and S waves are all body waves but surface waves could also be present in real data. Real data are more accurately described by elastic rather than acoustic modeling, where P-velocity, S-velocity, density and attenuation are all considered. However it is still common practice to apply acoustic FWI on real data. By considering the acoustic wave propagation, S waves are neglected and only the P waves are modeled. The acoustic wave equation could be resolved in the time domain [Tarantola, 1984, Mora, 1987, Mora, 1989] or in the frequency domain [Pratt andWorthington, 1990, Sirgue andPratt, 2004]. The numerical resolution of the partial diferential wave equation could be achieved by numerical methods, such as inite diference [Virieux, 1986, Moczo et al., 2004], inite element [Marfurt, 1984], etc. When the spatial discretization is done using a regular grid, the inite-diference method is usually chosen because it is fast and easy to implement.

1.2.3

Current limitations of FWI FWI is capable of potentially imaging complex structures. However, it presents several limitations [Virieux and Operto, 2009]:

• Local minima. The objective function of FWI is highly nonlinear with respect to the model parameters and has many local minima [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF]. These local minima prevent gradient-based techniques from inding the global minimum if the initial model is far from the global solution. This is known as the cycle-skipping problem. A phase mismatch of less than one half of the wavelength is required between the exact model and the initial model, to avoid falling into a local minimum.

• Lack of low frequency data. [Claerbout, 1982] and [START_REF] Jannane | Wavelengths of earth structures that can be resolved from seismic relection data[END_REF]] irst demonstrate that relectivity model is linearly derived from the relected waves in the data, while the large-scale velocity model does not linearly depend on the data. It seems that the data are not sensitive to middle range scale wavelengths. To avoid the non-linearity, multiscale FWI is needed to reconstruct the model, starting from low frequencies and progressively adding higher frequencies. When the low frequencies are missed in the data, which is often the case with real data, the FWI usually fail to retrieve the long-wavelength part of the model, and only acts as a migration mode.

Figure 1-2: Separation of scales between the frequency content of the velocity model and of the seismic data [Claerbout, 1982].

(velocity, density, attenuation). Since the '90s, FWI has been mainly used for imaging the P-wave velocity model by considering the acoustic wave propagation. This approximation is quite valid when FWI is mainly based on direct arrivals (diving waves, transmissions), as these waves are less sensitive to density perturbations. But for relected data, which are sensitive to density perturbations, the P-wave velocity model does not always allow the reliable interpretation of the complex structure of the subsurface, as the density model, the S-wave velocity model and the attenuation factor all have a great inluence on the amplitude of relection data. [START_REF] Mulder | Exploring some issues in acoustic full waveform inversion[END_REF] study the efect of acoustic inversion of elastic data through a marine data set and conclude that variable-density acoustic inversion of marine data can have some value if the subsurface is not too complex and the target is not too deep.

For complex models, multi-parameter inversion should be considered. However simultaneous inversion of several parameters is still challenging. The main obstacle is that diferent parameters have coupling efects and diferent orders of magnitude, making the inversion poorly conditioned and more nonlinear. In [START_REF] Operto | A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice[END_REF], Prieux et al., 2013, Malinowski et al., 2007], the authors propose guidelines to choose suitable strategies for multi-parameter FWI.

Alternative formulations

To mitigate the non-linearity of FWI, many alternatives formulations are proposed, which are mainly related with:

• Choice of the representations of the data, such as logarithm of the complex value of the data [START_REF] Shin | Waveform inversion using a logarithmic waveield[END_REF], separation of data amplitude and phase [Shin andMin, 2006, Bozdağ et al., 2011], envelope of the data [START_REF] Chi | Full waveform inversion method using envelope objective function without low frequency data[END_REF], Wu et al., 2014, Bozdağ et al., 2011], normalized integration of the data [START_REF] Donno | Estimating the background velocity model with the normalized integration method[END_REF], energy lux of the data [Causse, 2002]. Among these possible alternatives, the Laplace-domain FWI [START_REF] Shin | Waveform inversion in the Laplace domain[END_REF] has been shown to be efective to build a smooth velocity model. The authors show that by transforming the wave in the Laplace domain, the objective function of the ℓ 2 norm of the logarithmic waveield appears to be more convex and artiicial frequencies smaller than the frequencies in the source wavelet are created.

• Choice of data residual norm. The ℓ 2 norm is based on the assumption of Gaussian distribution of data uncertainties [Tarantola, 2005]. It may be not valid for all cases. Moreover, it is sensitive to large errors [Tarantola, 2005]. ℓ 1 norm [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF], which is not based on Gaussian assumption, has been proven to be more robust in the presence of noise in the data [START_REF] Brossier | Which data residual norm for robust elastic frequency-domain full waveform inversion[END_REF]. Besides, as ℓ 1 norm ignores the amplitude in the residual, the inversion is less sensitive to the large data error [Virieux and Operto, 2009]. We can also cite the Huber norm [START_REF] Guitton | [END_REF]Symes, 2003, Ha et al., 2009], the sech norm and the Cauchy criterion norm [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF], Amundsen, 1991]. These alternative norms could be considered as intermediate between ℓ 1 norm and ℓ 2 norm.

• Choice of minimization criterion. The conventional criterion measures the leastsquares data misit which could be quite nonlinear with respect to model. Other criterion are proposed to mitigate the nonlinearity of the problem, such as data diferential based criterion [START_REF] Chauris | Investigating the diferential waveform inversion[END_REF] or data correlation based criterion [van Leeuwen and Mulder, 2008]. [van Leeuwen and Mulder, 2008] show that compared to the least-squares functional, most of the weighted norms of the correlation function have a large basin of attraction and respond smoothly to a perturbation of the true velocity model.

Relection-based waveform inversion 1.3.1 Principles

When dealing with relection data, the use of Relection-based Waveform Inversion (RWI) [Mora, 1989, Chavent et al., 1994, Xu et al., 2012, Biondi et al., 2012, Tang and Lee, 2013, Wang et al., 2013, Brossier et al., 2015, Zhou et al., 2015[START_REF] Wang | Inversion of seismic refraction and relection data for building long-wavelength velocity models[END_REF] appears promising.

The method relies on the explicit separation of the model into a background model and a relectivity model. The relectivity model is inverted irst by migration and then serves as primary sources in depth, which allows computing the transmission wavepaths from the relectors to the surface. These transmission wavepaths are useful for the update of the background model. A new relectivity model must be generated by migration according to the background velocity model update at each nonlinear iteration and the two-step worklow is iterated until a ixed convergence criterion is met [START_REF] Brossier | Velocity model building from seismic relection data by full-waveform inversion[END_REF]. Migration Based Traveltime Tomography (MBTT) [START_REF] Chavent | Automatic determination of velocities via migration-based traveltime waveform inversion: A synthetic data example[END_REF] is based on a combination of migration and modelling. It combines the least-squares migration and the multiscale FWI to mitigate the nonlinearity. [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF] use the waveield decomposition to separate the migration part and the tomography part in the gradient. Their idea is to mix these two parts together and set diferent weights to enhance the tomography part. [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF]] determine, at each iteration, the weights by solving an optimization problem, which is not a trivial task. The advantage is that, it is easier to compute the misit function at each iteration and then to control the convergence. [START_REF] Brossier | Velocity model building from seismic relection data by full-waveform inversion[END_REF] propose a relection-based full waveform inversion. They use a correlation-based misit function instead of the classic least-squares data misit function to avoid cycle-skipping. Besides, the inversion is performed in the pseudo-time domain instead of depth domain to avoid recomputing the relectivity at each iteration of the inversion. [START_REF] Zhou | Combining diving and relected waves for velocity model building by waveform inversion[END_REF] propose to introduce diving waves along with relections to improve the reconstruction of the shallow parts of the model, which in turn improves the imaging of the deeper parts. [START_REF] Alkhalifah | The natural combination of full and image-based waveform inversion[END_REF] In recent years, RWI has shown its robustness compared to conventional FWI. However, this method is not mature yet, and several limitations could be listed. Firstly, as RWI is based on the relection data, other waves, such as diving waves, refractions, which are useful for the determination of the long-wavelength model, are neglected. [START_REF] Zhou | Combining diving and relected waves for velocity model building by waveform inversion[END_REF] propose to combine the relections with the diving waves to enhance the long-wavelength components of the model. [START_REF] Wang | Inversion of seismic refraction and relection data for building long-wavelength velocity models[END_REF] propose to combine the relections with refractions. In [START_REF] Wang | Taking advantage of waveield decomposition in full waveform inversion[END_REF] and [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF] the proposed methods do not use relections explicitly and could a priori incorporate diving waves and refractions automatically. Further investigations on this point are needed. Secondly, during the inversion, at each nonlinear iteration a new relectivity must be recomputed, causing a slow convergence rate. [START_REF] Brossier | Velocity model building from seismic relection data by full-waveform inversion[END_REF] propose to perform the migration in the pseudo-time domain as the zero-ofset traveltime is preserved during the inversion, to avoid the re-computation of the relectivity model at each iteration.

Thirdly, the nonlinear step is sensitive to the relectivity model obtained in the migration step. Thus building a true-amplitude relectivity is important. The relectivity is usually obtained by migration methods. However, as the migration uses the adjoint of forward-modeling operator instead of the inverse operator, a single iteration does not necessarily preserve the true amplitude of migration images. Alternatively Least-Squares Migration (LSM) can be used. LSM is an iterative migration method resolving the linearized inversion problem. It has been shown to improve amplitude information and to better focus migrated images. Several authors show that the LSM provides more reliable results than conventional migration images [Clapp, 2005, Valenciano, 2008, Zhang et al., 2005].

1.4

Ob jective and outline of the thesis Seismic relection data contain traveltime information as well as relection-amplitude information, therefore waveform inversion of relections should in principle be able to recover both the long wavelengths and the short wavelengths of the model [Snieder et al., 1989, Hicks andPratt, 2001]. However, in practice, the quality of the inversion is degraded due to the band-limited and ofset-limited data. In the time domain, when a large frequency range is used at once, these two components are coupled during FWI. This coupling is a major problem for FWI [START_REF] Snieder | Retrieving both the impedance contrast and background velocity: A global strategy for the seismic relection problem[END_REF], as the long-wavelength update usually has a smaller amplitude than the short-wavelength update. Therefore, when these two components are mixed together, the velocity update is mainly driven by the shortwavelength update. It is therefore natural to try to invert separately the long-wavelength and the short-wavelength components of the velocity model. Following [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] and [START_REF] Zhou | Acoustic wave-equation traveltime and waveform inversion of crosshole seismic data[END_REF], we propose a two-step inversion worklow achieved by decomposing the gradient formula of FWI into the long wavelength part and short wavelength part [START_REF] Wang | Taking advantage of waveield decomposition in full waveform inversion[END_REF]. From the literature, the FWI gradient is computed as the cross-correlation of the forward propagated source waveield and the back-propagated residual waveield [Tarantola, 1984, Lailly, 1983, Plessix, 2006]. With two-way modeling, both down-going and up-going components are present in the propagating waveields. After decomposition of the forward and back-propagated waveields into their down-and up-going components, the correlation of the two initial waveields actually gives four terms. The back-scattered correlations provide the positions of the relectors (short-wavelength components), while the forward-scattered correlations give information along the propagation paths (long-wavelength components). For Reverse-Time Migration (RTM) imaging, the forward-scattered correlations are usually removed because they create long-wavelength artifacts in the inal migrated image [Liu et al., 2011, Yoon andMarfurt, 2006]. But for velocity analysis, these correlations help in updating the long wavelength components of the velocity model [START_REF] Diaz | Understanding the reverse time migration backscattering: noise or signal?[END_REF], Wang et al., 2013, Brossier et al., 2015].

In Chapter 2, I review the basics of FWI, including the objective function, and the resolution of the direct problem using the inite-diference scheme. The calculation of the gradient of the objective function and the available local optimization methods are also detailed.

In Chapter 3, I outline the proposed inversion procedure, and show how the two separate inversion steps, either for short-wavelength or long-wavelength update, are formulated. I illustrate the method using a 1D layer model. Diferent waveield decomposition techniques are also introduced, one in the Fourier domain, one using the Poynting vector and an alternative in the curvelet domain.

In Chapter 4, I apply FWI, multiscale FWI and the proposed inversion method on two synthetic models and compare the results. I also analyze the inluence of key parameters during the inversion.

In Chapter 5, I apply FWI and the proposed inversion method on the Brunei real dataset, and present the preprocessings on the data and the diiculties we meet with the acoustic assumption.

Publications:

• Wang, F., Chauris, H., Donno, D., and Un problème inverse vise à reconstruire le modèle du sous-sol à partir d'un ensemble de données observées. L'imagerie sismique déinit un problème inverse en géophysique. En entrée, les données sismiques sont enregistrées aux positions des récepteurs. Les solutions que nous recherchons sont les propriétés physiques du sous-sol. Le problème directe consiste à simuler les données sismiques dans un modèle connu. L'inversion des formes d'ondes est une procédure pour ajuster les données. Le principe de FWI est d'utiliser tous les types d'ondes (ondes directes, rélexions, réfractions, multiples, . . . ) pour résoudre les diférents paramètres du modèle (la vitesse, la densité, l'atténuation). Le champ d'onde n'a pas besoin d'être décomposé en temps de trajet et amplitudes par exemple.

FWI a été proposée par [Lailly, 1983] et [Tarantola, 1984]. Cette méthode est largement utilisée dans la sismologie globale et dans l'exploration pétrolière. Parmi les très grand nombre de références, nous pouvons citer [Fichtner, 2010] et [Virieux and Operto, 2009] pour une revue tutorielle de cette méthode.

Dans ce chapitre, nous introduisons les formulations classiques de FWI, y compris la fonction objective, la résolution du problème direct en utilisant le schéma des diférences inies. Le calcul du gradient de la fonction objective et les méthodes d'optimisations locales sont également détaillées. L'analyse de résolution du gradient est mise en oeuvre pour étudier l'inluence de la coniguration d'inversion sur la résolution du gradient.

La fonction objective de la FWI vise à minimiser l'écart entre les données enregistrées et les données calculées [Tarantola, 1984]. FWI est généralement considérée comme un problème d'optimisation. Vue l'échelle du modèle, souvent les méthodes d'optimisation locales sont utilisées.

Le problème direct de FWI consiste à modéliser la propagation des données sismiques dans un modèle donné. Les données observées sont obtenues par la résolution de l'équation des ondes. Dans le cadre de cette thèse, nous nous limitons à l'équation des ondes acoustiques. L'équation des ondes acoustiques peut être résolue dans le domaine temporel [Tarantola, 1984, Mora, 1987, Mora, 1989] ou fréquentiel [Pratt andWorthington, 1990, Sirgue andPratt, 2004]. Elle peut être résolue analytiquement ou numériquement. Les méthodes numériques, par exemple, diférences inies [Virieux, 1986, Crase et al., 1990, Moczo et al., 2004], éléments inis [Marfurt, 1984, Choi et al., 2008], sont généralement utilisées. Dans cette thèse, nous utilisons la méthode des diférences inies, d'ordre 4 en temps et d'ordre 8 en espace. La méthode des diférences inies est un moyen naturel pour résoudre des équations aux dérivées partielles car elle est rapide et facile à mettre en oeuvre. Elle ne repose pas sur l'hypothèse haute fréquence. Elle est donc plus générale que les méthodes basées sur les rais et fournit une implémentation relativement eicace par rapport aux méthodes des éléments inis. La méthode des diférences inies estime directement chaque terme diférentiel en utilisant le développement de Taylor sur une grille régulière. En plus, nous utilisons les conditions absorbantes aux bords du modèle pour simuler un milieu inini.

Pour le problème d'inversion de données réelles, la signature de la source est généralement inconnue et doit être estimée en même temps que les paramètres du modèle. Dans cette thèse, nous utilisons la méthode proposée par [Pratt, 1999] pour estimer la source.

La mise à jour du modèle par FWI est obtenue avec le gradient et le Hessien de la fonction objective. Le gradient est la dérivée de premier ordre de la fonction objective par rapport aux paramètres du modèle. Le gradient relie la perturbation des données à la perturbation du modèle. La mise à jour du modèle est basée sur le gradient. En fait, la direction inverse du gradient donne la plus grande pente de la fonction objective. Le gradient peut être calculé eicacement avec la méthode état adjoint [Plessix, 2006], et est obtenu comme la corrélation entre le champ direct et le champ rétropropagé des résidus. Dans ce sens, le gradient de FWI est similaire au principe d'imagerie proposée par [Claerbout, 1971].

Le Hessien est la dérivée du deuxième ordre de la fonction objective par rapport au modèle. L'inverse du Hessien permet de corriger le facteur d'atténuation géométrique et illumination inégale, et de déconvoluer le gradient causé par la bande limitée des données [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF], Brossier et al., 2009].

La résolution de l'inversion dépend directement de la résolution du gradient. La résolution du gradient est fonction des diférents types d'ondes et des dispositifs de l'acquisition. Concrètement, les ondes directes, les ondes plongeantes, les ondes transmises et les multiples donnent une résolution plus basse que les ondes de première rélexion [START_REF] Huang | Resolution limits for wave equation imaging[END_REF], donc fournissent un gradient plus basse fréquence. Pour l'acquisition, l'inversion avec une source plus basses fréquences [Claerbout, 1985, Jannane et al., 1989] et des données de plus grands ofsets [Mora, 1989, Pratt et al., 1996] fournissent un modèle plus basse fréquence. D'où viennent les stratégies d'inversion multi-échelles [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF], Sirgue, 2003] et la stratégie de multi-étapes qui commence par les grands ofsets pour d'abord retrouver les grandes longueurs d'ondes du modèle.

Selon les méthodes utilisées pour calculer l'inverse du Hessien, les méthodes de mise à jour du modèle de vitesse peuvent être divisées en plusieurs catégories. Les méthodes basées sur le gradient supposent que le Hessien est un diagonal ne dépendant que d'un scalaire. Elles donnent une solution simple. La méthode du gradient conjugué combine la direction actuelle et la direction du gradient précédente pour accélérer la convergence. La méthode de Newton calcule le Hessien complet, tandis que le Gauss-Newton ne considère que le premier terme du Hessien. La méthode de quasi-Newton [Nocedal, 1980] repose sur l'estimation du Hessien, permettant de proiter de l'information contenue dans le Hessien aux itérations préc dentes, sans payer le prix d'une vraie itération du Newton complet. L-BFGS est une variante de BFGS. Elle garde en mémoire seulement un nombre limité d'itérations, généralement entre 3 à 15 itérations, étant moins exigeante en mémoire que BFGS. L'algorithme L-BFGS sera utilisé dans les applications présentées dans les chapitres 4 et 5 de cette thèse.

Le succès des méthodes d'optimisation locales dépend de l'exactitude du modèle initial. Comme les fonctions objectives peuvent présenter des minima locaux, le modèle initial doit être situé proche du minimum global pour assurer la convergence . Si le déphasage entre les données calculées et les données observées ne dépasse pas la moitié de la période du signal, les méthodes d'optimisation ajusteront ces deux données sans ambiguïté de phase. Sinon, les méthodes d'optimisation ajusteront ces deux données avec un déphasage d'une ou plusieurs périodes, provoquant la convergence vers un minimum local. Ceci est le phénomène de "cycle-skipping". Le modèle initial pour la FWI peut être obtenu par tomographie du temps de trajet ou par l'analyse de vitesse par migration. Pour atténuer la non linéarité, la méthode multi-échelle est souvent utilisée.

Introduction

An inverse problem aims at reconstructing the underlying model from a set of observed data. Seismic imaging deines an inverse problem in geophysics. The input is the recorded seismic data at receivers positions, and the solution we search for are the physical properties of the subsurface. The forward problem consists of simulating the seismic data in a known model. Full Waveform Inversion (FWI) is an inverse problem and is deined as a data-itting procedure. The principle of FWI is to use all types of waves (direct waves, relections, refractions, and multiples) to resolve diferent model parameters (velocity, density, attenuation). The information of the complete waveield: travel time, phase and amplitude are all taken into consideration. Full waveform inversion was pioneered by [Lailly, 1983] and [Tarantola, 1984]. This method is widely used in global seismology and in oil exploration. Among the extremely large number of references, we refer to [Fichtner, 2010] and [Virieux and Operto, 2009] for a tutorial review of the method. In this chapter, I review the classic Full Waveform Inversion (FWI) formulation, including the objective function, the resolution of the direct problem using the inite-diference scheme. The derivation of the gradient of the objective function and the local optimization methods are also detailed. The resolution analysis of the gradient is carried out to study the inluence of inversion setup on the gradient resolution.

2.2

Ob jective function FWI is stated as a non-linear optimization problem. [Tarantola, 2005] demonstrates that, by assuming Gaussian distribution of the data uncertainties, the maximum of the probability density function of the model is achieved by minimizing the least-squares functional of data residual. Based on this probabilistic study, the objective function of FWI is formulated, aiming at minimizing the diferences between simulated and observed data with the ℓ 2 norm. The objective function () in the time domain is therefore written as [Tarantola, 1984]:

() = 1 2  s,r  T 0  cal ()(s, r, ) -obs (s, r, )  2 , (2.1)
where is the model parameter, s = ( x , z ) and r = ( x , z ) are the shot and receiver positions, is the recorded time length, is the time variable, cal are the simulated data, and obs are the observed data, both recorded at the surface. The least-squares diferences are summed over all the sources s and all the receivers r.

To simplify, equation (2.1) can be re-written in vectorial notation as:

() = 1 2 ∆ t • ∆, (2.2)
with ∆ = calobs denoting the data residual, and (.) t denoting the transpose operator.

In the time domain, the seismic data are represented by the seismograms of dimension t x r , where t is the number of recorded time samples and r is the number of seismic traces within a shot gather. The objective of the inversion is to minimize the functional () in equation ( 2.2) by modifying the model . In order to ease the extraction of information from the data and to improve the algorithm convergence, a weight is usually applied to the data. The weighted least-squares objective function is deined as:

() = 1 2 ∆ t t ∆.
(2.3)

Commonly used weights are data selection functions that, for example, give emphasis to long/short ofsets, select events around the irst arrival, etc...

Forward problem 2.3.1 Wave equation

The forward problem in FWI consists of modeling the seismic data propagation in a given model . The simulated data cal is obtained by solving the wave equation. The physical law relating the model parameters and the data is the visco-elastic wave equation. For the scope of this thesis, we limit ourselves to the simplest case, by considering the constant density acoustic approximation of the wave propagation, which reads:

1 2 2 2 -∆ = , (2.4)
where is the velocity, is the waveield, and is the source term. at the surface provides cal . In the acoustic approximation, is generally the pressure. We discuss in Chapter 4 the efect of variable density.

Numerical resolution

The acoustic wave equation can be resolved in the time domain [Tarantola, 1984, Mora, 1987, Mora, 1989] or in the frequency domain [Pratt andWorthington, 1990, Sirgue andPratt, 2004]. It can be resolved analytically or numerically. In the time domain, the numerical resolution of the partial diferential wave equation could be achieved by numerical methods, such as inite diference [Virieux, 1986, Crase et al., 1990, Moczo et al., 2004], inite element [Marfurt, 1984, Choi et al., 2008], discontinuous Galerkin [Cockburn et al., 2000, Dumbser and[START_REF] Dumbser | [END_REF]. [Virieux et al., 2009] and [Plessix, 2007, Plessix, 2009] present some analyses of the complexity of diferent numerical modeling methods. The inite-diference method is a natural way to solve partial diferential equations. It is not based on high frequency assumption, thus it is more general than ray-based methods [START_REF] Cerveny | Ray method in seismology[END_REF], and it provides a relatively eicient implementation compared to inite-element approaches. The method directly estimates each diferential term using the Taylor expansion on a regular grid. Suppose that the 2D model dimension is model x model , the spatial samplings are ∆ and ∆, the recording time is , the time step is and the number of time step is t , then the 2 nd order Taylor expansion of a waveield at the th time step provides the estimation of the 2 nd order temporal derivative of as follows:

2 ( ij , n ) 2 ≃ n+1 ij -2 n ij + n-1 ij ∆ 2 , (2.5)
where , are the spatial coordinates, and n is the temporal coordinate. Similarly, we get the 2 nd order spatial derivative of :

∆( ij , n ) ≃ ( n i+1j + n ij+1 ) -4 n ij + ( n i-1j + n ij-1 ) ∆ℎ 2 ,
(2.6) by setting ∆ℎ = ∆ = ∆.

Inserting equation ( 2.6) and equation ( 2.5) into equation (2.4) gives the explicit resolution of at time step + 1 and spatial position (, ) using values of at previous iterations.

n+1 ij = 2 n ij -n-1 ij + 2 ∆ 2 ∆ℎ 2 ( n i+1j + n ij+1 + n i-1j + n ij-1 -4 n ij ) + 2 ∆ 2 .
(2.7)

In this work we use the scheme with 4 th order in time and 8 th order in space for the wave propagation modeling. It is based on the same principle as the above formulations. This wave propagation modeling operator is auto-adjoint as explained later.

In this conventional inite-diference implementation, 5 coeicients (including the central position and the 4 neighboring positions on the standard Cartesian coordinate system) are considered for 2D modeling and 7 coeicients (including the central position and the 6 neighboring positions on the standard Cartesian coordinate system) are considered for 3D modeling. [START_REF] Operto | 3D inite-diference frequency-domain modeling of viscoacoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF] propose to use not only the neighboring positions on the standard Cartesian coordinate system, but also those on the 45 ∘ rotated Cartesian coordinate system. For this implementation, 9 coeicients are considered for 2D modeling and 27 coeicients are considered for 3D modeling. [Virieux, 1986] irst introduced the inite-diference scheme on a staggered grid for modeling the seismic rupture problem. Note that for constant density, velocity and pressure are deined at the same position while for variable density, the density is known at half-integer positions, and the velocity is known at integer positions. This staggered grid scheme resolves the problems with instabilities in models with high velocity contrasts and with grid dispersion in media with high Poisson ratio, at least for the conventional grid inite diference scheme.

The spatial and temporal sampling should be chosen carefully to avoid numerical dispersion and to ensure the stability of the solution. The spatial sampling ∆ℎ is deined such that

∆ℎ ≤ min λ = min max 1 λ , (2.8)
where min corresponds to the smallest wavelength, which depends on the maximum frequency max and the minimum velocity min . λ denotes the number of grid points per wavelength. The number of points per wavelength depends on the spatial and temporal order in the Taylor expansion, usually 5.5 points per wavelength is chosen for order 8 in space and order 4 in time.

The temporal sampling ∆ is deined by the Courant-Friedrichs-Lewy (CFL) condition to ensure the stability of the forward propagation:

∆ max  1 (∆) 2 + 1 (∆) 2 ≤ 1, (2.9)
where max is the maximum velocity, and ∆ and ∆ are the spatial discretizations.

The spatial discretization of the inite diference method is achieved on a regular grid, making it simple to implement. Nevertheless, the utilization of regular grid also gives rise to certain limitations as the spatial sampling is ixed and not adapted to local physical properties, thus the method is not as lexible as numerical methods using irregular grids, such as inite elements and discontinuous Galerkin. [Marfurt, 1984] compares the accuracy of inite-diference and inite-element for wave propagation modeling and observes that inite-element scheme better models curvi-linear interfaces and thin beds than inite-diference scheme. [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF] propose the extension of the initediference method with loating interfaces on regular grids. However, for the scope of this work, we consider the classical inite-diference implementation.

Boundary conditions

Equation (2.4) deines the wave propagation in an ininite domain. However, because of the limited acquisition extension and of the limited computational capacity, the models we consider are usually limited in a inite domain. If we do not deine any boundary condition at the borders of the model, waves will be artiicially relected from the borders of the model. To attenuate these parasite waves, spatial boundary conditions must be applied.

Absorbing boundary conditions are widely used to simulate an ininite medium. Ideally, absorbing boundary conditions should let pass the waveields coming out of the model but prevent waveields coming back into the model. Some types of absorbing boundary conditions rely on the modiication of local attenuation properties by adding additional layers around the model [START_REF] Cerjan | A nonrelecting boundary condition for discrete acoustic and elastic wave equations[END_REF]. The additional layers work as sponges. However, these sponges layers do not always work perfectly. Waves impinging with normal incidence are optimally absorbed by the sponges, however parasite relections would appear for other incident angles [START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. [Berenger, 1994] introduced a highly efective absorbing boundary condition based on the Perfectly Matched Layer (PML) for electromagnetics. In the continuous limit, it is proven that a PML interface between a regular medium and the ictitious perfectly matched medium completely absorbs incident waves from the regular medium, regardless of their incidence angle.

To be realistic, free surface boundary condition is also required, which describes the free surface of the Earth, thus simulating relections from the Earth-air interface. However, in this thesis, I do not deine the free surface, and I consider absorbing boundary layers all around the model.

Source signature estimation

For a real data inversion problem, the source signature is usually unknown and has to be estimated along with the model parameters. [Tarantola, 1984] estimates the source using the adjoint-state method along with the other model parameters. The source wavelet is resolved iteratively using the gradient and the Hessian of the objective function with respect to the source wavelet.

[ Pratt, 1999] propose a source signature estimation method for each discrete frequency during the iterative inversion procedure, by resolving a linear inverse problem as the relationship between the waveield and the source is linear (refer to equation (2.4)). The solution of the linear inverse problem is given by:

= t * t * ,
(2.10) where = cal is the simulated waveield in the actual velocity model and = obs is the observed waveield. The source signature could be estimated during the inversion procedure once the forward waveield is modeled. In principle, the source should be estimated at each iteration of the inversion process. For the synthetic examples shown in Chapter 4, the source is deined as a given Ricker wavelet, while for the real data test in Chapter 5, the source wavelet is estimated according to equation (2.10).

We have validated the implementation of the method described in equation ( 2.10) on a 2D synthetic example. A Ricker wavelet with central frequency of 12 Hz is used for modeling the observed data. The initial guess of the source wavelet is a Dirac wavelet with frequency bandwidth between 0 Hz and 40 Hz. For each frequency between 0 Hz and 40 Hz, we estimate the source wavelet using equation (2.10). Figure 2-2 displays the initial Ricker wavelet and the estimated wavelet in the frequency domain and in the time domain. We note that the reconstructed source wavelet is very close to the original one. The velocity is correct at the surface in this case.

2.5

Linearization of the inverse problem FWI is considered as a non-linear optimization problem. Usually a local optimization method is used. If we assume that the data is locally approximately linear with the model and that the initial model 0 is close to the true one, then the objective function is locally approximately quadratic in the neighbourhood of the starting model 0 , and we search for the minimum of the objective function in the neighbourhood of 0 , with

= 0 + ∆.
(2.11) Suppose the model perturbation vector ∆ is relatively small compared to model vector , then the second-order Taylor expansion gives:

() = ( 0 +∆) = ( 0 )+∇ m ( 0 )∆+ 1 2 ∆ t (∇ m ) ( 0 )∆+(∆ 3 ). (2.12)
We denote the Hessian matrix as

= (∇ m ) ,
(2.13) and we express the derivative of with respect to using equation ( 2.12) with the locally linear assumption:

∇ m () = ∇ m ( 0 ) + ( 0 )∆.
(2.14)

The minimum of the objective function is reached at if the gradient at this location is zero. By setting ∇ m () to zero we have:

∆ = --1 ∇ m ( 0 ).
(2.15)

The irst term -1 is the pseudo-inverse of the Hessian matrix. The second derivative term is the gradient vector of the objective function. By substituting equation ( 2.15) into equation (2.11), we get the Newton solution that allows updating the initial model [Tarantola, 1984, Virieux andOperto, 2009].

If the objective function is quadratic, leading to the linear inverse problem case, the convergence of the problem could be achieved in one iteration using equation (2.15). If the objective function is not exactly quadratic, the Newton method will not converge in one iteration, and the problem must be resolved iteratively to compensate the problem non-linearity. In this case, the model update is given by:

k+1 = k --1 ( k )∇ m ( k ).
(2.16)

() = ∇ m () =  cal () )  t ( cal () -obs ) = () t ∆,
(2.17)

where is the Fréchet derivative matrix. If we consider the linear direct problem, represents the Born modeling operator, with ∆ = ∆ [Tarantola, 2005].

As shown in equation (2.15), the resolution of local optimization problems is based on the gradient. In fact, the opposite direction of the gradient gives the largest slope of the objective function. [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF] gives a physical interpretation of the gradient vector by illustrating the partial derivative waveield using a point difraction model. They conclude that the partial derivative waveield represents forward propagation using virtual sources at the difracting locations. The gradient is the correlation of the partial derivative waveield and the residual vector. This correlation tests the similarity between these two waveields and allows getting the difracting points which give rise to the data residuals. These difracting points are initially not present in the model. [Sirgue, 2003] gives an interpretation of the gradient vector by performing the Singular Value Decomposition (SVD) on the Fréchet derivative matrix. He expresses the gradient vector in the model eigenvector basis and concludes that the gradient vector is the sum of the projected components of true model perturbation along the model eigenvectors, multiplied by the model eigenvalues. Thus, the gradient direction is the direction of the true model perturbation in a stretched model space.

Instead of calculating the Fréchet derivative , which is not trivial, the gradient can be eiciently calculated with the adjoint-state method [Plessix, 2006], by using the augmented functional and Lagrange multipliers. With the adjoint-state method, and with the fact that the modeling operator is auto-adjoint, the computation of the gradient of the objective function is simpliied and is given by:

= ∇ m () = - 2 3  s  T 0 (s, ) 2 (s, ) 2 = ⊗ .
(2.18)

From equation ( 2.18), we can observe that the construction of the gradient requires two waveields, the forward propagated waveield , resolving the direct wave equation (2.4), and the back-propagated waveield , satisfying the adjoint wave equation:

1 2 2 2 -∆ = cal -obs , (2.19)
with the data residual calobs as the source term. This equation is back propagated with ( max = 0). The gradient is computed as the cross-correlation of and . This formula is similar to the migration imaging principle of [Claerbout, 1971]. The gradient obtained after the irst iteration of the local optimization looks like a migrated image obtained by Reverse-Time Migration (RTM).

Hessian and preconditioning

The Hessian is a matrix of dimensions model x model . It is the second-order derivative of the objective function with respect to the model parameter :

() = 2 () 2 = t () () + () ∆.
(2.20)

From equation ( 2.15) we observe that the Hessian relates the gradient with the true model perturbation. The gradient only does not give a correct model update direction as the computed gradient may sufer from geometrical spreading and characteristics from the source and acquisition setup. The inverse of Hessian in equation ( 2.15) helps to correct the geometrical spreading, uneven illumination and deconvolves the gradient [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF], Brossier et al., 2009].

Referring to equation ( 2.20), we can explicitly separate the Hessian into two terms. Of the two terms, the irst term is straightforward to compute if the Fréchet derivatives are available, while the second term is much more diicult to compute. [START_REF] Métivier | Full waveform inversion and the truncated newton method[END_REF] propose a matrix-free method to compute the full Hessian term using 2D adjoint-state method and inner linear optimization. If the full inverse of Hessian is computed, it yields the Full-Newton method. If only the irst term is considered, it is the approximate Hessian and it yields the Gauss-Newton method [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF].

The diagonal terms of the approximate Hessian of the Gauss-Newton solution contain the zero-lag autocorrelations and therefore represent the square of the amplitude of the partial-derivative waveield. Applying the inverse of these diagonal terms to the gradient partly compensates the geometric spreading efect. In the framework of surface seismic experiments, this scaling provides a good balance between shallow and deep perturbations. The of-diagonal terms of the approximate Hessian are computed by correlation between partial-derivative waveields associated with diferent model parameters.

For 1D media, the approximate Hessian is a band-diagonal matrix, and the numerical bandwidth decreases as the frequency increases. For a inite range of frequencies, however, the Hessian is no longer diagonal and not even diagonally dominant [Chavent andPlessix, 1999, Pratt et al., 1998]. The of-diagonal elements of the approximate Hessian account for the limited-bandwidth efects related to the source term. Applying its inverse to the gradient can be interpreted as a deconvolution of the gradient from these limitedbandwidth efects. This deconvolution of gradient enables to focus and sharpen the FWI images [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF].

The explicit resolution of the inverse of the approximate Hessian matrix is diicult, as the matrix size is too large to be directly used in practical applications. Sometimes, it is even impossible, as the Hessian may be not positive and requires regularization. Some authors propose to consider only the diagonal terms of the approximated Hessian [Shin et al., 2001, Chavent and[START_REF] Chavent | [END_REF]. But as the gradient of the least-squares functional multiplied by a diagonal matrix cannot perfectly deconvolve the band-limited efect, iterations are needed to gradually construct the inverse of the approximate Hessian [START_REF] Lambaré | Iterative asymptotic inversion in the acoustic approximation[END_REF], Nemeth et al., 1999]. Some approximations of the diagonal of the approximate Hessian are also proposed to construct amplitude-preserving migration images [Shin et al., 2001, Plessix andMulder, 2004]. [Virieux and Operto, 2009] illustrates the scaling efect performed by the diagonal elements of the approximate Hessian and the deconvolution efect performed by the ofdiagonal elements of the approximate Hessian. They show that the inversion image is sharper when using the approximate Hessian instead of setting it as a scalar.

It is also possible to precondition the gradient without using the Hessian, but by exploiting some a priori information. For example, weighting the data residuals in time or in depth linearly is an appropriate way to compensate the geometrical spreading.

2.7

Resolution analysis

Gradient formulation

The gradient of the objective function of FWI is calculated as the cross-correlation of the forward waveield and back-propagated waveield . The gradient formulation (equation (2.18)) can also be written in the frequency domain using Green's function as [START_REF] Sirgue | Eicient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF]]:

= -2  s  r ( * 0 (x, s) * 0 (x, r)∆(r, s)), (2.21)
where 0 (x, s) and 0 (x, r) are the Green's functions for the source and receiver waveield, and ∆ is the data residual at the receiver positions.

As the velocity update of FWI is based on the gradient, the resolution analysis of the gradient is necessary to study the inluence of diferent parameters on the spectrum of the gradient, and to ind the reason why FWI fails to recover the long wavelengths of the model. [START_REF] Huang | Resolution limits for wave equation imaging[END_REF] derive formulas of resolution limits for diferent types of waves, including diving waves, primary relections, multiple relections, and difractions. The resolution analysis was initially realized to study seismic migration. However we can beneit from it for the spectral analysis of waveform inversion. Here, I give a summary of the analysis and the conclusions. Readers could refer to [START_REF] Huang | Resolution limits for wave equation imaging[END_REF] for a detailed derivation.

Resolution analysis of diferent waves

Resolution analysis of a scattering point

In the framework of the inverse-scattering and difraction tomography [Devaney, 1982, Miller et al., 1987], the spatial frequency vector k sampled in the image at a difraction point is related to the wavelength and the half-scattering angle . We denote the directions of the the source and the receiver plane wave by d s and d r . With the assumption of the far ield, and by replacing the Green's function with the plane-wave approximation, for a single frequency we have:

0 (x, s) ≈ ( 0 d s • x), 0 (x, r) ≈ ( 0 d r • x).
(2.22) (2.23)

The coverage of the wavenumber k of the gradient at the scattering point (Figure 23) is deined according to:

k = 0 (d s + d r ) = 2 cos n = 4 cos n, (2.24)
where is the velocity and is the angular frequency. According to equation (2.24), the resolution is optimal for zero-ofset data, and is equal to λ 2 . For FWI, low frequency (small ) and long ofset (large ) data provide long-wavelength velocity update, while high frequency (large ) and short ofset (small ) data provide short-wavelength velocity update. The resolution analysis of other types of waves is also based on this singlescattering resolution analysis.

Resolution analysis of primary relections

For primary relections, in the far-ield approximation, the horizontal resolution limit ∆ and the vertical resolution limit ∆ (Figure 2-4e) are deined according to:

 ∆ = Zλ 4X , ∆ = λ 4 , (2.25)
where is the relector depth, is the source-receiver ofset, and is the wavelength.

Resolution analysis of diving waves

For diving waves, the vertical limitation ∆ is equal to half width of the irst Fresnel zone at the point middle-way between the source and the receiver (Figure 2-4a). It depends on the travel length and the wavelength :

2∆ = √ , (2.26)

Resolution analysis of relection-related transmission

The relection-related transmission provides the rabbit-ear wavepaths shown in Figure 2-4b and 2-4c. The resolution analysis is similar to that of the diving waves. It is calculated using the mirror-relection scheme (Figure 2345) as if the receiver is located at the mirror position (0, 2d) in a homogeneous model. The resolution limit 2∆ perpendicular to the ray at the midpoint is equal to the width of the irst Fresnel zone, given by:

2∆ = √ , (2.27)
where is the travel length and is the wavelength.

Resolution analysis of difraction-related transmission

The difraction-related transmission is similar to the relection-based transmission, except that the difraction point is considered as the virtual receiver and the secondary source. The propagation distance is halved compared to the relection transmission, leading to a narrower wavepath with the same travel length and wavelength . The resolution limit is given by (Figure 2-4d):

2∆ =  /2.
(2.28)

Resolution analysis of multiple-related transmission

The multiple-related transmission simulates the repeated wavepath between the source and the relector and between internal relectors. The resolution limit is computed according Figure 2-6: Wavepath for the irst-order free-surface multiple [START_REF] Huang | Resolution limits for wave equation imaging[END_REF].

to the mirror-relection scheme (Figure 2-6 and 2-7):

2∆ = √ , (2.29)
where is the travel length and is the wavelength. From the previous analysis, we conclude that the resolution limit of diferent waves mainly depends on the scattering angle, the wavelength, and the wave propagation distance. These parameters are function of the seismic acquisition setup. For example, long ofsets are favorable for recording wide-angle and long-propagation distance data.

Resolution analysis and acquisition setup

In this section, we show the inluence of the source bandwidth and the data ofset range on the wavenumber components of the FWI gradient through a 1D layer model example.

Resolution analysis and source bandwidth

We perform conventional FWI in the time domain for the 1D layer model. The initial model is homogeneous at 2 km/s. We analyze four FWI tests with diferent source wavelets. These four source wavelets are all Ricker wavelets, but the central frequencies for the four tests are respectively 4 Hz, 8 Hz, 12 Hz and 16 Hz (Figure 23456789). Here, low frequencies are present in these source wavelets. The inversion is performed in the time domain with all frequencies mixed together. Figure 2-10 shows the gradients calculated at the irst iteration of FWI for each test and the wavenumbers of the vertical slice taken in the middle of each gradient model. We observe that the spectrum of the gradient is closely related to the bandwidth of the source wavelet and the frequency bandwidth in the gradient increases with the increase of the frequency bandwidth in the source wavelet. This result corresponds to equation 2.24 and to the conclusion that for the same acquisition geometry, lower frequency data provide longer-wavelength velocity update, while higher frequency data provide shorter-wavelength velocity update. In order to retrieve the long wavelengths of the model, low-frequency data are favorable [Sirgue, 2003]. And when low frequencies are missing from the data, which is often the case in real data case, the FWI usually fails to retrieve the long-wavelength part of the model [Claerbout, 1985, Jannane et al., 1989].

Resolution analysis and ofset [Mora, 1989] and [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF] show that FWI acts in two modes, a migration mode, related to the short wavelengths of the model which could be obtained by iterative migration, as well as a tomography mode, related to the long wavelengths of the model which could be retrieved with the beneit of long-ofset data.

We compute the gradient of FWI using the 1D layer model of Figure 2-8, and using data with diferent ofset ranges. The central frequency of the Ricket wavelet is 12 Hz and the initial model is homogeneous. The ofset ranges are respectively [0-0.17] km, From the analysis in sections 2.7.2 and 2.7.3, we conclude that in order to retrieve the long-wavelength components of the model, low frequency and long-ofset data are important [Sirgue, 2006]. Based on this idea, ofset-windowing [Mora, 1989, Pratt et al., 1996] strategies are proposed for FWI. Moreover, relection-related transmission [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF], Tang and Lee, 2013, Wang et al., 2013, Brossier et al., 2015], and multiple-related transmission [START_REF] Staal | Velocity estimation using internal multiples[END_REF] can also be taken into consideration for the reconstruction of the long-wavelength model.

2.8

Velocity model update methods

Supposing that the initial model 0 is situated in the basin of attraction of the global minimum of the objective function, the local optimization methods should be able to converge to this global minimum in a limited number of iterations. When added to the initial velocity, the velocity perturbations calculated using equation (2.11) lead to an updated velocity model, which is used as the starting model for the next iteration.

Depending on the methods used to calculate the Hessian matrix, the velocity model update methods could be divided into several categories.

Steepest descent (or gradient descent) methods

The gradient descent method assumes that the Hessian is a scalar, giving a simple solution for the model perturbation vector:

∆ = -.
(2.30)

This method does not require the estimation of the inverse of the Hessian and only the descent step length needs to be determined. The role of the step length can also be considered as converting the unit of the gradient vector to model unit [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF]]. Generally, is determined by line-search methods [Nocedal, 1980, Gauthier et al., 1986, Tarantola, 2005] by minimizing the objective function with respect to . Alternatively, can be estimated by quadratic or cubic interpolation [START_REF] Vigh | Developing earth models with full waveform inversion[END_REF]. The idea is to determine a pair ( 1 , 2 ) such that

 ( k + 1 ∆) < ( k ) ( k + 2 ∆) > ( k + 1 ∆) .
(2.31)

The quadratic interpolation passing through these three points

( k , k + 1 ∆, k + 2 ∆
) provides the minimum of the parabola in the interval [0, 2 ], giving the optimal descent step opt .

Knowing that the gradient represents the direction of the steepest descent of the objective function, we can always ind an that reduces the objective function in the opposite direction of the gradient.

Conjugate-gradient method

The convergence rate of the gradient descent method is generally quite slow as the successive directions estimated by the gradient are orthogonal if the step length is optimally chosen. The convergence rate can be improved by the conjugate-gradient method, which combines the current gradient direction with the previous gradient direction to accelerate the convergence rate [Mora, 1987, Crase et al., 1990, Hu et al., 2011].

The current gradient direction k uses the present gradient direction k and previous conjugate-gradient direction k-1 :

k = k + k k-1 .
(2.32) There are various ways to determine . One well-known choice of is expressed as [START_REF] Polak | Note sur la convergence de méthodes de directions conjuguées[END_REF]

k = ( k -k-1 ) t k ( k-1 ) t k-1 .
(2.33)

Newton and Gauss-Newton methods

The Newton method computes the full Hessian, while the Gauss-Newton only considers the irst term of the Hessian in equation 2.20 [Pratt et al., 1998, Virieux andOperto, 2009]. These two methods have a faster convergence rate than the previously described gradient-based methods but they require the calculation of extra direct problems for the calculation of the inverse of the Hessian and are considered quite time-consuming for real data inverse problems. [START_REF] Métivier | Full waveform inversion and the truncated newton method[END_REF] propose a matrix-free method to compute the full Newton term using 2D adjoint-state method and inner linear optimization.

2.8.4

Quasi-Newton method

The quasi-Newton method relies on the estimation of the Hessian, allowing to take advantage of the information contained in the Hessian without paying the price of a full Newton iteration. The BFGS (named after Broyden, Fletcher, Goldfarb and Shanno) algorithm [Nocedal, 1980] estimates the inverse of the Hessian from gradients and objective functions of previous iterations and increases the precision of the estimation along with iterations.

The limited-memory version (L-BFGS) of the algorithm [Nocedal, 1980] is a variant of BFGS. It keeps in memory only a limited number of iterations, typically between 3 to 15 iterations, being less demanding in memory than BFGS. The L-BFGS algorithm will be used in the applications presented in Chapter 4 and Chapter 5 of this thesis. [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF]] review diferent optimization methods and compare their quality and convergence rate. [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF] show that the gradient-based methods are less robust than the Newton-based methods, and may fail to converge no matter how many iterations. [START_REF] Métivier | Toward Gauss-Newton and exact Newton optimization for full waveform inversion[END_REF], Castellanos et al., 2015] compare the FWI results using Newton, Gauss-Newton and quasi-Newton methods. Figure 2-12 shows an example of the optimization scheme for diferent optimization methods.

2.9

Initial model and non-linearities

The success of the local optimization methods depends on the accuracy of the initial model.

As the objective functions may present some local minima, the initial model should be located in the neighbourhood of the global minimum to ensure the success of convergence.

If the phase shift between the calculated data and the observed data does not exceed half period of the signal, the optimization method will adjust these two data without phase ambiguity. Otherwise, the optimization method will adjust these two data with a phase shift of one or several periods (Figure 2-13), causing the convergence to a local minimum. This is the so-called cycle-skipping phenomenon.

The initial model for FWI can be obtained by traveltime tomography [START_REF] Bishop | Tomographic determination of velocity and depth in laterally varying media[END_REF], Pratt and Chapman, 1992, Billette and Lambaré, 1998, Woodward et al., 2008] or Migration Velocity Analysis (MVA) [START_REF] Symes | Velocity inversion by diferential semblance optimization[END_REF], Chauris and Noble, 2001, Sava and Biondi, 2004]. Traveltime tomography uses the traveltimes to estimate the wave propagation velocity. Generally, RMO events are picked after prestack depth migration [START_REF] Guillaume | Kinematic invariants: an eicient and lexible approach for velocity model building[END_REF]. Migration Velocity Analysis (MVA) is performed in the image domain and aims at building a relationship between migrated image perturbation and model perturbation. This approach has the advantage of not requiring picking.

In waveform inversion, strategies of using diferent data at diferent stages of inversion are proposed to mitigate the non-linearity of the problem. For the same acquisition geometry, lower frequency data provide longer-wavelength velocity update, while higher frequency data provide shorter-wavelength velocity update [Sirgue, 2003]. Based on this idea, [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF] perform the inversion in the time domain using scale decomposed preconditioned with the diagonal term of the inverse Hessian [Sirgue, 2003].

source and data, starting from the lowest available frequencies and gradually increasing the frequency bandwidth. This point is also illustrated in [START_REF] Mulder | Exploring some issues in acoustic full waveform inversion[END_REF]] by showing that the basin of attraction of the objective function becomes narrower for increasing frequencies. Similarly to [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF], the multiscale inversion can be implemented in the frequency domain using increasing frequencies [Sirgue, 2003, Mulder andPlessix, 2008]. Guidelines to properly choose successive frequencies or frequency ranges are also provided [Sirgue andPratt, 2004, Brossier et al., 2009]. The elastic inversion using the multiscale scheme is carried out by [START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. [Virieux and Operto, 2009]. If the phase mismatch is less than one half of the period, the local optimization methods could adjust the correct phase. In the contrary case, the local optimization methods adjust these two signals with one period shift and the inversion falls into a local minimum. [Virieux and Operto, 2009] For each frequency, we perform 25 iterations of FWI. The waveield is initially computed in the time domain and then the chosen discrete frequency is extracted. Compared to the conventional FWI, the multiscale FWI (Figure 2-18 and Figure 2345678910111213141516171819) better recovers the long wavelength of the anomaly.

Another strategy to avoid local minima is data-windowing [Sirgue, 2003, Brenders andPratt, 2007]. It consists of inverting the early arrivals irst as they provide lower wave numbers in the velocity gradient [START_REF] Sears | Elastic full waveform inversion of multi-component OBC seismic data[END_REF]. They mainly update the longwavelength components of the shallow part of the model [Sirgue, 2003]. Time windowing of early arrivals is achieved by damping the amplitudes of later arrivals in the data residuals.

The strategy of selecting data with diferent ofsets is also proposed. [Mora, 1989] and [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF] observe that the short wavelengths of the model could be obtained by iterative migration, and the long wavelength of the model could be retrieved with Based on this idea, [START_REF] Shipp | Two dimensional full waveield inversion of wide-aperture marine seismic streamer data[END_REF] propose that the long-ofset data should be used in the early stages of inversion, in order to recover the long-wavelength components of the model. At later stages, short-ofset data are inverted to provide the short-wavelength components. This strategy of data windowing from long ofsets to short ofsets is eicient, provided that there is no cycle skipping in the long-ofset data. However, as long-ofset data correspond to long propagation distance, they introduce increasing non-linearities in an inaccurate model, which makes the application of FWI to long-ofset data still challenging [Sirgue, 2006]. The opposite strategy of windowing the data from short-ofset to long-ofset in order to mitigate the risk of cycle skipping has been also proposed [Sirgue, 2003].

Conclusion

I have reviewed diferent aspects of classical FWI, including the objective function, the methods for the numerical modeling of waveield propagation, the computation of the gradient and the Hessian of the objective function, and the velocity model update methods. I have also presented a resolution analysis of the gradient and show that long-ofset data and source with low frequencies help to recover the long-wavelength part of the model. Finally, I have reviewed the diiculties that FWI meet and some methods to mitigate the nonlinearity of FWI.

Chapter 3

Waveform inversion based on waveield decomposition La fonction objective de l'inversion des formes d'ondes peut être fortement non linéaire et présente beaucoup de minima locaux. Si les données contiennent principalement des rélexions, ce problème empêche notamment les méthodes classiques basées sur le gradient de récupérer les grandes longueurs d'ondes du modèle [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF]. Avec l'analyse de la résolution du gradient présentée dans le chapitre 2, nous présentons ici une variation de l'approche FWI en proitant de la transmission liée à la rélexion pour construire les grandes longueurs d'ondes du modèle. La méthode est basée sur la séparation du gradient en une partie de courtes longueurs Cette méthode d'inversion basée sur la décomposition est similaire á d'autres méthodes d'inversion basées sur la rélexion. Nous comparons en particulier avec la méthode proposée par [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF]. La comparaison des gradients calculés dans un modèle 1D montre la similarité entre ces deux méthodes.

Les méthodes de décomposition des champs d'ondes sont également présentées dans ce chapitre. Nous discutons de trois méthodes de décomposition : l'une dans le domaine de Fourier, l'une utilisant le vecteur de Poynting et la troisième dans le domaine des curvelets. La méthode par la transformé de Fourier est la plus populaire, et a été largement util-isée pour supprimer des artefacts de haute fréquence dans les images migrées, mais avec cette méthode les basses fréquences sont iltrées dans les résultat inal. La méthode de Poynting est la plus rapide car elle pourrait être calculée au cours de la propagation des ondes, mais elle soufre du problème d'artefacts numériques fortes (instabilités). La méthode de décomposition par les curvelets donne le meilleur résultat car il peut prendre en compte les diférentes directions en même temps, avec le coût de calcul le plus élevé des trois. Dans cette méthode, le champ d'onde est décomposé en champs d'onde localement plans. La corrélation permet de savoir dans quelle direction l'onde se propage, sans avoir à faire d'hypothèse sur la direction de propagation. Ce n'est pas le cas avec la méthode de Fourier pour laquelle il faut savoir si les ondes se propagent essentiellement horizontalement ou verticalement. La méthodologie a été développée pour corréler les champs directement dans le domaine des curvelets. La nouvelle méthode d'inversion et les méthodes de décomposition sont illustrées par un modèle de couche 1D. L'application de la nouvelle méthode d'inversion aux modèles synthétiques et aux données réelles est abordée dans les chapitres 4 et 5.

Introduction

The objective function of the Full Waveform Inversion (FWI) may be highly nonlinear and has many local minima. If a data set mainly contains relections, this problem particularly prevents the classical gradient-based methods from recovering the long wavelengths of the model [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF]. In this chapter, we present a variation of FWI by separating the gradient of FWI into a migration part and a tomographic part. We call this new method Decomposition-based Waveform Inversion (DWI). The scale separation in the gradient is based on the waveield decomposition into their one-way components. The waveield decomposition methods are also presented in this chapter. The new inversion method and the decomposition methods are illustrated through a 1D layer model. The application of the new inversion method on synthetic models and real data set will be addressed in Chapter 4 and Chapter 5.

Methodology

Based on the resolution analysis of the gradient presented in Chapter 2, we propose a new inversion method taking advantage of the relection-related transmission to construct the long wavelengths of the model. The method is based on the gradient separation into a short-wavelength part and a long-wavelength part. The scale separation in the gradient is achieved by the waveield decomposition. Given a non-smooth model , both forward waveield and back-propagated waveield contain up-and down-going waveields, as shown in Figure 3-1 (top panel). They can therefore be separated into:

 = + + -, = + + -, (3.1)
where + is the direct forward waveield, -is the relected forward waveield, -is the direct backward waveield and + is the relected backward waveield. By substituting the decomposed waveields in equation ( 3.1) into the gradient formula (equation (2.18)), the gradient can be decomposed into four terms:

= ⊗ = + ⊗ -+ -⊗ + + + ⊗ + + -⊗ -, (3.2)
where ⊗ indicates the correlation between two waveields. For the irst two terms, the waveields + and -, as well as -and + , propagate in the same directions (both downward or upward) and coincide in time at the positions of relectors (Figures 3-1 second and third panel). Thus, their correlations + ⊗ -and -⊗ + provide the shortwavelength update of the model (relectivity) and correspond to the migration imaging principle proposed by [Claerbout, 1971]. The other two terms + ⊗ + and -⊗ -are illustrated in Figure 3-1 (fourth and ifth panel). As the back-propagated waveield is reversed in time, waveields + and + , as well as -and -propagate in the opposite direction and coincide in time along the whole propagation path that connects the source to the relector, and the relector to the receiver. Thus, their correlations provide the long wavelength update of the model (background velocity) [START_REF] Zhou | Acoustic wave-equation traveltime and waveform inversion of crosshole seismic data[END_REF].

Suppose that the velocity model consists of a background model 0 and a perturba- tion model : = 0 + . If all frequencies are considered at the same time, the gradient of the short-wavelength update is given by

r = + ⊗ -+ -⊗ + , (3.3)
while the long-wavelength update is given by

m 0 = + ⊗ + + -⊗ -.
(3.4)

Inversion strategy

Note that the conventional FWI gradient formula (equation (3.2)) sums the long-wavelength term m 0 and short wavelength term r with equal weights. For relection-dominant data, the long-wavelength term m 0 is usually much weaker than the short wavelength term r , of a factor 2 , where is a typical relectivity value. [Mora, 1987] and [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF] propose to enhance the long-wavelength term m 0 in the gradient by applying diferent weights to r and m 0 . This is not trivial as the weighting factor is diicult to determine and may vary with models and iterations. [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF] propose to determine the weight at each iteration by solving a linear optimization problem which involves the calculation of the inverse of the Hessian. Instead of setting diferent weights, we propose to separately invert the short and long wavelengths of the model. Let start with a smooth initial velocity model 0 . First we retrieve the short wavelengths of the model using equation (3.2) instead of equation (3.3) to reduce the computing cost, as it avoids to decompose the waveields into diferent parts. This is feasible because the gradient amplitudes are dominated by the short wavelengths especially in the irst stages of the inversion. It is important to note that in this step, we use zero or short ofset data and we perform iterative inversion. The short-wavelength model obtained in this irst step contains relectors which are needed for the successive step.

Note that the short-wavelength update step behaves like a least-squares migration, but it is not exactly the same as the separation of the long wavelength and the short wavelength is not explicit in this step, thus the inversion is not linearized. Instead, we perform the iterative inversion on zero-ofset data with all frequencies simultaneously. The use of short-ofset data is crucial because the long-wavelength components of the model are not updated during this step, as no moveout information is available. The iterative inversion helps to correct the relector amplitude as well as to achieve better focusing of the relectors [START_REF] Plessix | Frequency-domain inite-diference amplitude -preserving migration[END_REF], which allows for a more accurate long wavelength update. However, the relectors are mispositioned in depth at this stage as the used background velocity model is probably erroneous.

After the iterative inversion, the forward waveield and backward waveield found at the end of the irst step are now decomposed into their up-and down-going components. We perform the long-wavelength update using equation (3.4). In this step, longer ofset data is used to get the moveout information. A smoothing ilter is applied to the gradient, in order to remove the residual short-wavelength components, that are caused by acquisition limitations. We carry out one iteration of inversion in this step and the velocity is updated using the smoothed gradient.

Once we have updated the long-wavelength model, we use this smooth model to update once more the short-wavelength model using the iterative inversion, as the positions of relectors are not consistent with the actual velocity model any more and should be updated accordingly. At each global iteration, the relectivity is set to zero before being estimated. Then we continue alternating between these two steps. The key point of the proposed inversion method is this alternating scheme. The comparison of the conventional FWI and the proposed method is shown in Figure 3-2.

The ofset limitation and the application of iterative inversion in the irst step and the smoothing of the gradient in the second step are the key issues to ensure the decoupling between the long wavelength and short wavelength updates of the velocity model, as we will describe in more details in the following subsections. A 1D-layer model will be used to better illustrate the method.

Choice of ofset range

In the irst step of the inversion, we use zero-ofset or short-ofset data in order to preserve the vertical traveltime. In this step, we do not want to introduce the long-wavelength update. As the zero-ofset traveltime is preserved, there is no moveout information in the data residual.

Iterative inversion

In the irst step of inversion, we perform the iterative inversion to derive an optimal relectivity model. As the second step of inversion is sensitive to the relectivity model obtained in the irst step, building a true-amplitude relectivity is important. Otherwise, residuals at short ofsets will be introduced in the long-wavelength update and degrade the inversion result. An example describing this issue will be shown in Chapter 4. The relectivity is usually obtained by migration methods. However, as the migration uses the adjoint of forward-modeling operator instead of the inverse operator, they could not preserve the true amplitude of migration images. For Relection-based Waveform Inversion (RWI), usually Least-Squares Migration (LSM) is used in the irst step of inversion. LSM is an iterative migration method resolving the linearized inversion problem. It has been shown to improve amplitude information and to focus migrated images. The sensitivity of the update of the macro model to the relectivity model is also analyzed in [START_REF] Lameloise | Improving the gradient of the image-domain objective function using quantitative migration for a more robust migration velocity analysis[END_REF], Symes, 2015] by showing that the global minimum of the objective function of the Diferential Semblance Optimization (DSO) does not necessary coincide with the exact model due to the artifacts present in the relectivity model, and that the inversion result can be improved by performing the iterative migration or the quantitative migration preceding the update of the long-wavelength part of the model. In our method, as the scale separation is not explicit, instead of performing the LSM, we perform the nonlinear inversion in the irst step to maintain the compatibility between the two steps of inversion.

3.3.3

Gradient smoothing

In the second step of the inversion, we try to retrieve the long-wavelength components of the model. In order to remove the residual short-wavelength components in the gradient, we apply a smoothing operator to the gradient. After 2D Fourier transform of the gradient, a 2D low-pass elliptic ilter is applied to the transformed gradient. The elliptic ilter is deined by the model dimension and the maximum preserved wavenumber. A taper zone is applied to prevent artifacts caused by discontinuity. After iltering the gradient in the wavenumber domain, an inverse Fourier transform is performed to recover the iltered gradient in the space domain.

1D layer model case

Before presenting the waveield decomposition scheme, let use the 1D layered model in Figure 2-8 to study the characteristics of the proposed method. A Ricker wavelet with central frequency of 12 Hz and temporal sampling of 1.5 ms is used as the source function.

The initial velocity model is displayed in Figure 3-3. This is the relectivity model and the model contains mainly short-wavelength components. Here, we artiicially use the exact relectivity to analyse the diferent waveields. Later, the relectivity model is obtained through iterative inversion. Starting from this model we calculate the forward waveield and backpropagated waveield for the single sources located at 0.85 km at the surface. Usually, all the available source positions are used, but in this case we use a single source to more easily illustrate the behavior of the proposed algorithm. Then, we decompose and into their up-and down-going parts to get four correlation terms:

+ ⊗ -, -⊗ + , + ⊗ + and -⊗ -. The decomposition is achieved using the 2D Fourier-transform based method, which will be described in the next section. For one source and one receiver setup, the four correlation terms are shown in Figure 3-4 and Figure 3-5. In these igures, the source is ixed at 0.85 km and the receivers are at three 3-8 on the right. We can observe that the spectrum of ⊗ contains both high wavenumbers and low wavenumbers. Instead, for the decomposed case, the spectra of + ⊗ -and -⊗ + mainly contain the high-wavenumber part of ⊗ , while + ⊗ + and -⊗ - mainly contain the low-wavenumber part of ⊗ .

3.4

Comparison with similar methods

The decomposition-based FWI has many common features with other Relection-based Waveform Inversion (RWI) methods. Here, we particularly compare our method to the one proposed by [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] as an example. The method proposed by [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] is based on a linearization (Born approximation) and on the explicit separation of the velocity model into a smooth background model 0 and a relectivity model . The forward problem is described by the wave equation:

() = ()(x -s), (3.5)
where () = -2 -∆ is the forward propagator. is the waveield, composed of 0 and , with = 0 + . 0 is the waveield propagating in the smooth model 0 and is the waveield generated by the model perturbation . The objective function is deined as: where obs is the observed waveield generated by the true model perturbation. At irst, we deduce by using the Born approximation for the waveields and 0 :

() = 1 2 ‖ -obs ‖ 2 , ( 3
 () = ()(x -s), ( 0 ) 0 = ()(x -s).
(3.7)

By subtracting these two equations, we get the formula to compute :

( 0 ) = 2 0 .
(3.8)

The adjoint-state method [Plessix, 2006] is used to compute the gradient of the objective function in equation (3.6) with respect to 0 and . Two adjoint variables 0 and are introduced and the augmented function  is:

 = 1 2 ‖ -obs ‖ 2 +    [( 0 ) 0 -()(x -s)]dsdxd +    0 [( 0 ) -2 0 ]dsdxd.
(3.9)

We compute the gradient of  with respect to 0 and respectively:

 =    (-2 0 (s, x, ) 0 (s, x, ))dsdxd, (3.10)
 0 =    [- 2
0 (s, x, )(s, x, ) -2 (s, x, ) 0 (s, x, )]dsdxd. (3.11) 0 and are computed by the two back-propagation wave equations:

 ∂  J ∂δP = * 0 + ( -obs ) = 0, ∂  J ∂P 0 = * -2 0 = 0.
(3.12)

where * denotes the adjoint operator.

From equation (3.11) we observe that four waveields are involved in the gradient calculation for the background model 0 . They are respectively the direct source waveield 0 and the direct receiver waveield 0 , both propagating from the surface to the subsurface relector, as well as the perturbed source waveield and perturbed receiver waveield , both propagating from the relector back to the surface. The relectivity serves as a secondary source in depth that emits up-going waves towards the surface.

The decomposition-based gradient in equation (3.4), related to the long-wavelength update, has a similar structure compared to the gradient in equation (3.11). We irst show the comparison using the model in Figure 3-3. We compute the original forward waveield and the back-propagated waveield (Figure 3-9), then the decomposed waveields as described in equation (3.4) (Figure 3-10). The four components in the long-wavelength gradient (equation (3.11)) are shown in Figure 3-11. The gradient calculated with all the sources and receivers, using equation (3.11) is displayed in Figure 3-12.

If we compare the correlation of + and + , with the correlation of 0 and , we can observe that with the time, they both move along the wavepath from the surface to the relector. The same conclusion could be drawn when we compare the correlation of and -, with the correlation of and 0 , while they both move along the wavepath from the relector to the surface. The gradients in Figure 3-12 and in Figure 3-8 (last two panels) are also similar. In this sense, we could consider decomposition-based FWI as an analogue of the method proposed by [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF]. The mayor diference resides in how to calculate the four components in the gradient. For the method proposed by [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] they are obtained by resolving the linearized wave equation, and for the method proposed in this thesis, they are obtained by decomposing the full waveields after resolving the nonlinear wave equation.

However, as the method of [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] is based on the linearized modeling, it assumes that the velocity perturbations are small and the data contain only primaries. While similarly to classical FWI, the proposed decomposition-based FWI method resolves the nonlinear wave equation; therefore refractions, diving waves and multiples could potentially be included in the dataset. [Symes, 2008] suggests that MVA is a solution method for the linearized waveform inversion problem.

Besides the method of [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF], there are in the literature other similar methods, also based on separating diferent components in the gradient. [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF] also propose the gradient separation by the waveiled decomposition. The diference is that they recombines the two parts in the gradient with diferent weights. [START_REF] Brossier | Velocity model building from seismic relection data by full-waveform inversion[END_REF], Zhou et al., 2015] propose to get the tomographic part by subtracting the migration part that is obtained in the smooth model from the full gradient. [START_REF] Zhou | Combining diving and relected waves for velocity model building by waveform inversion[END_REF] propose to combine the relections with the diving waves to enhance the long-wavelength components of the model. [START_REF] Wang | Inversion of seismic refraction and relection data for building long-wavelength velocity models[END_REF] propose to combine the refractions with relections. For relection data, the long-wavelength part is updated using the tomographic part in the gradient by including the migration part within the density model. [Alkhalifah, 2015] propose to do the separation in the angle domain by iltering the gradient in the wavenumber domain according to the scattering angles and then to use the gradient corresponding to big scattering angles for long-wavelength update and the gradient corresponding to small scattering angles for short-wavelength update. Migration Based Traveltime Tomography (MBTT) [START_REF] Chavent | Automatic determination of velocities via migration-based traveltime waveform inversion: A synthetic data example[END_REF] combines the least-squares migration and the multiscale FWI to mitigate the nonlinearity.

Waveield decomposition methods

The gradient decomposition into four terms is possible with the beneit of the waveield decomposition. The objective of the waveied decomposition is to obtain the one-way component of the waveield with respect to the vertical direction (donwgoing/upgoing) and the horizontal direction (leftgoing/rightgoing). The waveield decomposition is mainly used in the literature for:

1. Providing a new imaging condition. The conventional zero-lag cross-correlation imaging condition of RTM is subject to strong migration artifacts [START_REF] Liu | An efective imaging condition for reverse-time migration using waveield decomposition[END_REF].

The decomposition of the full waveields to their one-way components and the application of the imaging condition to the appropriate combinations of the waveield components could remove the undesired noise.

2. Retrieving the long wavelength components of the velocity model. The gradient of FWI could be separated into a migration part and a tomography part, after the waveield decomposition. Thus, the waveield decomposition is used for recovering the long-wavelength part of the model [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF], Wang et al., 2013, Wang et al., 2015a].

Here, we give a review of diferent methods to produce partial waveields or to decompose the full waveields. In general, the waveield decomposition methods could be divided into one-way based methods and two-way based methods.

Methods based on one-way wave equation

To produce one-way waveields, the straightest way is to use the one-way approximation of the full wave equation. The one-way wave equation is easy to implement, but has dip limitations. In [START_REF] Zhang | Imaging complex salt bodies with turning-wave one-way wave equation[END_REF], the authors propose a modiied version of one-way wave equation to decompose waveields. With this method, turning waves can be properly imaged and the imaging capability of one-way wave equation is greatly improved.

Methods based on two-way wave equation

Methods based on two-way wave equation has no dip limitation, and are able to deal with multiples and turning waves. They produce a better resolution and are more robust with respect to strong velocity variations. There are two major categories to obtain the one-way component of waveields using the two-way wave equation. The irst one is by suppressing internal relections, and the second one is by decomposing waveields.

Suppression of internal relections

These methods aim at modifying the impedance of the media to suppress internal relections:

• Impedance matching. [Whitmore et al., 1983] and [Baysal et al., 1984] proposed a method to suppress undesired relections by forcing the acoustic impedance, = (, , )(, , ), to be a constant, and deriving an appropriate impedance-matched counterpart to the full wave equation.

• Smoothing of slowness ield [START_REF] Loewenthal | Suppressing the unwanted relections of the full wave equation[END_REF]. This method is based on the smoothing of the slowness ield throughout the media by an operator which is longer than the wavelengths associated with the acoustic wave. The reason of using slowness instead of velocity is that the integral of the original slowness curve and the integral of the smoothed slowness curve are the same, while the integral of the smoothed velocity curve difers from the original one and causes a shift in the arrival time. • Add of a damping term [START_REF] Fletcher | Suppressing artifacts in prestack reverse time migration[END_REF]. This method applies a directional damping term to the wave equation in areas of the velocity model where unwanted relections occur. A linear derivative operator is used. This linear derivative operator radiates in one direction and damps in the opposite direction.

Decomposition of waveields

These methods compute the full waveield and then decompose it according to the propagation direction. We analyze here three methods for achieving the decomposition of waveields, namely the 2D Fourier transform, the Poynting-vector method and the curvelet transform.

• 2D Fourier Transform (FT) decomposition.

[ [START_REF] Liu | An efective imaging condition for reverse-time migration using waveield decomposition[END_REF] applied waveield decomposition in the F-K domain to suppress the low-frequency artifacts in RTM images. This method was irst applied to vertical seismic proiles [START_REF] Hu | Wave-ield transformations of vertical seismic proiles[END_REF]. The waveield is decomposed into its upgoing and downgoing components by using the 2D Fourier transform:

 z+ (, z ) =   (, z ) for z ≥ 0 0 for z < 0 ,  z-(, z ) =  0 for z ≥ 0  (, z ) for z < 0 , (3.13)
where  (, z ) is the 2D Fourier transform of (, ). z+ (, ) and z-(, ) are the decomposed down-and up-going waveields. They are the inverse Fourier transform of  z+ (, z ) and  z-(, z ) respectively.

The horizontal decomposition is obtained similarly, by replacing by . 1. It can not handle vertical directions and horizontal directions at the same time.

 x+ (, x ) =   (, x ) for x ≥ 0 0 for x < 0 ,  x-(, x ) =  0 for x ≥ 0  (, x ) for x < 0 . ( 3
Usually, vertical decomposition is performed for models with sources and receivers at the surface. Therefore, in this case, artifacts will be generated by waveield propagating laterally. When a point source is used, direct waves propagate horizontally along the surface. That is why artefacts tend to appear near the surface.

2. As we separate data in thedomain, we need to apply a smooth cut in the domain (Figure 3-15) to reduce the artefacts due to discontinuities. Thus information near the zero axes in the F-K domain is partially iltered and this information correspond to the low frequency component of the waveield. Compared to Figure 3-13, the two images in Figure 3-17 look like a low-cut iltered version. However, applying the smooth cut has also the advantage of reducing the artefacts near the surface as presented in the previous paragraph (1).

• Poynting-vector decomposition.

Another way to decompose a waveield is by analyzing the instantaneous wave propagation directions of the waveields using the Poynting vector [START_REF] Yoon | Reverse-time migration using the Poynting vector[END_REF]. The Poynting vector measures the energy low and computes the wave propagation direction. It is deined as:

= -∇ . (3.15)
The principle of this method is to evaluate the sign of the product ∇ ∂P ∂t . This idea could be interpreted using the expression under the high frequency approximation of the waveield: (3.16) where is the traveltime between the source and the image point. Is satisies the Eikonal equation. (s, x) is the amplitude term, solution of the transport equation. () is the signature of the source. The term acting at the source depends on the dimension (1D, 2D, 3D). For the derivation of the spatial gradient, both and depends on x. But the smooth variation of can be neglected in front of the oscillatory kernel, yielding

≃ (s, x) (k) ( -(s, x)),
∇ ≃ -(s, x) (k+1) ( -(s, x))∇ (s, x), (3.17) and = (s, x) (k+1) ( -(s, x)).
(3.18)

Then the product of ∇ and ∂P ∂t yields

-∇ = 2 (s, x)[ (k+1) ( -(s, x))] 2 ∇. (3.19)
The sign of this product only depends on ∇ , which indicates the direction of waveield propagation with time. Therefore, the waveield is decomposed in the vertical direction as follows: We show an example of the implementation of this method using the same waveield as in the FT-based method.

z+ (, ) =  (
1. We calculate the partial derivative of the waveield with respect to and , and we get ∂P ∂z and ∂P ∂t . This can be done during the waveield propagation, thus requires no extra computation. Two snapshots of ∂P ∂z and ∂P ∂t are shown in Figure 3-19. 2. We then correlate these two waveields ∂P ∂z and ∂P ∂t . A snapshot of the product ∂P ∂z ∂P ∂t is shown in Figure 3-20. The original waveield is separated into downgoing and up-going parts according to the sign of the product at each location. According to equation (3.20), at each location, if the product is positive, then the original waveield at this location corresponds to the down-going part, otherwise, it corresponds to the up-going part. The decomposition result is shown in Figure 3-21. We notice that the down-going and up-going waveields are not completely separated. Parasite residuals are present, especially in the middle of the waveield. If we set the positive values in Figure 3-20 to 1 and negative values to 0, we get the image in Figure 3-22 on the left. In this image, points in black correspond to the locations of down-going waveield, while points in white correspond to the up-going waveield. This image betters shows the locations of artefacts. These artefacts are like "salt and pepper" artefacts, and they are due to the numerical approximation of the partial derivative. For the spatial positions whose partial derivatives are close to zero, the sign may be incorrect, thus yielding a wrong estimation of the mask. The points corresponding to near-zero partial derivatives are a) points at the intersections of downgoing and upgoing waveield; b) points propagating perpendicularly to the decomposition direction, for example, horizontally propagating waveield when the vertical decomposition is performed, which are usually near the surface; and c) points where is maximal, which are along the central part of the wavelet of the waveield. The artifacts decrease with the increasing of the spatial discretization, because the derivative value is better estimated. 3. We apply a median ilter to suppress the salt and pepper artefact of the mask.

The result is shown in Figure 3-22 on the right. The decomposed waveields using the median-iltered mask are shown in Figure 3-23. These two images present less artifacts than the initial results in Figure 3-21.

Figure 3-25: A direction example of curvelet [Chauris et al., 2006].

• Curvelet decomposition.

The main limitation of the FT-and Poynting-based methods is that they can not deal with the vertical and the horizontal directions at the same time. To circumvent this problem, we have designed a new method to decompose waveields using curvelets in order to attenuate the artifacts that FT-or Poynting-based decomposition methods cause.

Curvelets form a basis well-adapted for representing smooth images. Curvelets are multi-directional [Candes et al., 2006, Chauris and[START_REF] Chauris | [END_REF] and can decompose waveields in several directions at the same time. A curvelet is characterized by three parameters: scale, direction and position [START_REF] Candes | Fast discrete curvelet transforms[END_REF]. A direction and scale example of curvelets are shown in Figure 3-24 and 3-25. Curvelets form a tight frame, such that a function can be decomposed as a sum of curvelets. Besides, [START_REF] Candes | Curvelets: A surprisingly efective nonadaptive representation for objects with edges[END_REF] show that the curvelet representation of wave propagators is optimally sparse. For a 2D image with smooth discontinuities (, ), the reconstructed image n (, ), obtained with the irst most important coeicients, approximates the original image by:

‖ -n ‖ 2 L 2 ≈ (log) 3 2 .
(3.22)

The optimal convergence rate for all the decomposition methods is 1 n 2 . Therefore, compared to the convergence rate for Fourier transform ( 1 √ n ) and for the wavelet transform ( 1n ), the curvelet representation is more eicient. The curvelet is widely used in geophysics for data denoising [START_REF] Hennenfent | Seismic denoising with nonuniformly sampled curvelets[END_REF]] and data compression [Ma, 2011]. [START_REF] Chauris | Seismic demigration/migration in the curvelet domain[END_REF] apply the demigration operator in the curvelet domain to reduce computation cost.

From a practical point of view, we implement the curvelet transform as [START_REF] Chauris | Seismic demigration/migration in the curvelet domain[END_REF] did. First, we compute the 2D Fourier transform of the original image. Then, for all scales and all directions, we ilter the transformed image by the predeined curvelet functions. Finally, we apply an inverse 2D Fourier transform. This provides the curvelet coeicients. The inverse curvelet transform consists of applying a forward 2D Fourier transform for each scale and direction, iltering once more, adding the contributions for all scales and all directions, and performing an inverse Fourier transform to obtain the original image. In this thesis, we consider all scales at the same time, and the numbers of directions are 64 for all the scales. The number of directions could also be reduced by two at every two scales to satisfy the parabolic scaling [START_REF] Chauris | Seismic demigration/migration in the curvelet domain[END_REF]. To deine curvelet ilters, we use polar coordinates instead of Cartesian coordinates. This naturally leads to consider rotations instead of shears [START_REF] Candes | The curvelet representation of wave propagators is optimally sparse[END_REF] and avoids special considerations for handling edge efects. The design of the ilters is based on the approach proposed by [START_REF] Simoncelli | Shiftable multiscale transforms[END_REF].

Each ilter d,l (, ) is deined in the polar coordinate system (, ) and constructed as the product of two ilters, the direction ilter d () and the scale ilter l () (Figure 3-26). Suppose we have 2 N +1 diferent directions. The directional part d () is deined as follows:

d () =  2 √ 5 3 [ -P i 4 ( -1)] for ∈ [-P i 2N + P i 4N ( -1), P i 2N + P i 4N ( -1)] 0 otherwise . (3.23)
The factor 2 √ 5 is a normalization factor, and one could check that the sum of all squared ilters equal to 1. This ensures a perfect reconstruction scheme. The radial ilters l () are deined in a similar way, with a change of variable from to log, as described in [START_REF] Simoncelli | Shiftable multiscale transforms[END_REF].

Let 1 (x) and 2 (x) be two successive snapshots, for diferent 2D spatial positions x = (, ), as shown in Figure 3-27. The objective is to decompose the waveields in two diferent parts: a part for which the waveield expands from the source, and the remaining contracting part. More precisely, 1. We irst decompose the two images 1 (x) and 2 (x) in the curvelet domain for diferent directions. Let now consider two curvelets 1 d (x) and 2 d (x), for the selected direction .

2. We decompose the curvelet domain in sub-domains. For each sub-domain Ω i , the reconstructed images in the space domain for direction are 1 (x) and 2 (x) , obtained by

1 (x) =  Ω i dy 1 d (y) * d (x -y), (3.24) 2 (x) =  Ω i dy 2 d (y) * d (x -y), (3.25)
where d is the ilter deining the inverse curvelet transform for direction . This is simply a shifted version of the ilter, with a weight equal to the curvelet coeicient. The star * indicates that we consider the complex conjugate because the d coeicients and possibly the ilters are indeed complex.

3. We compute a correlation function d () that computes the correlation between 1 and 2 . The local correlation is evaluated along the direction, and characterized by the normal d as follows 3.26) where is a scalar value, denoting the shift along the direction. The meaning of d is the following: if the maximum of d is obtained for a positive value, then the waveield is considered as in expansion, otherwise in contraction.

d () =  dx 1 (x) * 2 (x + n d ), (
4. We repeat steps (1) to (3) for each direction and each sub-domain, and we sum all the contributions. The inal result is shown in Figure 3-29.

In the following, we will show that instead of reconstructing images from curvelet coeicients and then correlating the two images in the space domain, as done in equation (3.26), we could also achieve this goal in the curvelet domain directly. The correlation function only depends on the curvelet coeicients.

After step (1), in the curvelet domain, let select the coeicients belonging to a particular spatial zone Ω i . This is the zone where we want to evaluate the behavior of the waveields. For the evaluation of the local correlation, we consider coeicients in Ω i . For the correlation, the summation is performed over the full space. We replace 1 and 2 in equation (3.26) by their expression in equation (3.24) and 3.25 and interchange the order of integration, leading to

d () =  Ω i dy 1 d (y)  Ω i dz 2* d (z)  dx * d (x -y) d (x + n d -z).
(3.27)

In practice, the d ilters are very oscillatory in the spatial domain. The main contribution is thus obtained for z = y + n d . We thus have

d () ∝  Ω i dy 1 d (y)  Ω i dz 2* d (z)(z -y + n d ) (3.28) =  Ω i dy 1 d (y) 2* d (y + n d ).
(3.29)

We observe that the original correlation in the space domain (equation (3.26)) is equivalent to the correlation in the curvelet domain (equation (3.29)). Compared to the direct implementation in the space domain, we do not have to reconstruct the images for all directions nor for all sub-domains. All the computations are performed with curvelet coeicients. Note that by selecting coeicients in the curvelet domain, the reconstructed images would have natural tapering on the edges. In summary, the irst step consists of decomposing the two waveields in the curvelet domain, then of performing the local correlation. As before, depending on if the maximum value for d is obtained for a positive or negative , the selected coeicients are stored in two diferent arrays. After having investigated all spatial positions and all directions, the inverse curvelet transform is applied to the two intermediate arrays, providing two reconstructed images. For the curvelet-based decomposition, we can observe some artifacts of residual directional element as the curvelet is pseudo localized in the space domain. The artifacts get reduced when the number of directions considered in the curvelet decomposition increase.

We investigate a more complex waveield propagated in the Marmousi model. The decomposed results for the snapshot in Figure 3-30 using the Fourier transform method, the Poynting vector method and the curvelet method are shown in For this waveield, the drawbacks of the Fourier transform method and the Poynting vector method are obvious. We note that the Fourier transform method ilters out low frequency data and waveields propagating horizontally near the surface, and the Poynting vector method sufers from strong artefacts at intersections of diferent directions. The curvelet method yields the best result of the three.

From the point of view of computational cost, the Poynting vector method is the most eicient one, as the partial derivative waveields could be obtained during the wave propagation and the computational complexity is ( * ), where and are the model dimensions. The Fourier transform method is more expensive than the Poynting vector method as we have to compute the Fourier transform of the waveield in thedomain for each horizontal coordinate . If we use the FFT algorithm, the computational complexity is (***(*)), where and are the model dimensions, and is the number of time samples. The curvelet method is also expensive as we have to compute the curvelet transform for each snapshot of the waveield. Recall that the curvelet transform consists of one Fourier transform, one iltering for each direction and one inverse Fourier transform for each direction. The computational complexity is then (***(*)),

where and are the model dimensions, and is the number of time samples. Note that the Poynting vector method is the most eicient one. The diference of computational complexity between the Fourier transform method and the curvelet method is not obvious and should be evaluated case by case. In general, the decomposition based on the curvelet transform is more expensive than the decomposition based on the Fourier transform as the curvelet transform is about 10 times more expensive than the Fourier transform.

From the point of view of memory, for the Fourier transform method, we have to store the whole waveield on the hard disk. For the Poynting method, as the partial derivative waveields can be obtained during the wave propagation, we do not have to keep in memory the whole waveield. For the curvelet transform method, at each time, we compare two neighbouring snapshots, thus the decomposition can be done during the wave propagation and we do not need to store the whole waveield. Furthermore, since the decomposition is carried out for each snapshot separately, the curvelet decomposition method is more easily parallelizable. Further investigations on the curvelet implementation are needed before being able to apply it to large datasets.

For the tests shown in Chapters 4 and 5, we use the 2D Fourier transform decomposition method for computational cost and accuracy considerations. At the time these tests were performed, the correlation was done in the space domain but not in the curvelet domain, and the method was even more expensive. For this reason, the 2D Fourier domain was preferred. For 3D applications, the Poynting vector method may be the best choice, as the need of the Fourier transform method for storing the whole 3D waveield is diicult to meet. For the curvelet method, the signiicant computational cost of the 3D curvelet transform makes this method diicult to be widely used for large dataset applications.

Conclusion

In this chapter, we propose a relection-based waveform inversion method. The principle of the proposed method is to decompose the gradient of FWI into a long-wavelength part and a short-wavelength part. The inversion is performed in an alternating fashion between these two parts. This method allows to automatically separate the long and short wavelengths of the model without having a-priori knowledge of the model. The gradient decomposition is achieved by the waveield decomposition into their one-way components. Diferent waveield decomposition methods, especially those based on twoway wave equation, are presented. The 2D FT-based waveield decomposition method will be used for the 2D synthetic model in Chapter 4 and for the real data set in Chapter 5.

Résumé du chapitre 4

Dans ce chapitre, j'applique l'inversion des formes d'ondes (Full Waveform Inversion, FWI) et l'inversion basée sur la décomposition des champs d'ondes (Decomposition-based waveform inversion, DWI) sur des modèles synthétiques 1D et 2D. Rappelons que DWI consiste à décomposer le gradient en une partie de grandes longueurs d'ondes et une partie de courtes longueurs d'ondes après la décomposition des champs d'ondes. L'inversion est efectuée d'une manière alternée entre ces deux parties. Pour la modélisation numérique, nous utilisons un schéma aux diférences inies d'ordre 8 en espace et d'ordre 4 en temps. Des conditions aux limites absorbantes sont imposées aux bords du modèle pour simuler un milieu inini. Nous utilisons la méthode d'optimisation de quasi-Newton, plus précisément la L-BFGS avec 5 itérations conservées en mémoire. D'abord, nous testons FWI et DWI sur un modèle de couche 1D. Pour ce test, la source contient les basses fréquences, mais comme il y un grand écart entre le modèle exact et le modèle initial, la FWI classique n'arrive pas à retrouver les basses fréquences du modèle.

Ensuite, nous étudions un modèle synthèse 2D et comparons les résultats de l'inversion avec FWI, FWI multi-échelle et DWI. Nous faisons quatre tests avec deux choix de la source et deux choix du modèle de vitesse initial. Pour les deux modèles de vitesse initiale, l'un est proche, et l'autre est loin du modèle de vitesse exacte pour évaluer les résultats de l'inversion. Celui proche du modèle exact est un modèle de gradient constant, bien qu'il ne soit pas si proche du modèle exact, car il n'y a pas d'information a priori dans ce modèle. Celui loin du modèle exact est un modèle homogène. Pour la source, nous utilisons une source sans les basses fréquences pour être plus réaliste. Dans ce cas, la fonction objectif de FWI est moins linéaire, et l'inversion a plus de mal à converger. Cependant, nous avons également testé FWI en utilisant une source avec les basses fréquences comme une référence. La source qui contient les basses fréquences est une ondelette de Ricker avec une fréquence centrale de 8 Hz. Celle sans les basses fréquences est la même ondelette Ricker, mais les fréquences au dessous de 5 Hz sont enlevées. Notons que la densité reste constante pendant l'inversion. A travers des quatre tests, nous observons que la FWI multi-échelle est plus robuste que la FWI classique, et la DWI est plus robuste que la FWI classique et que la FWI multi-échelle face à l'absence de basses fréquences dans la source et aux erreurs importantes dans le modèle initial.

Nous avons également analysé l'inluence de plusieurs paramètres clés pendant l'inversion, comme l'application de l'inversion itérative non-linéaire dans l'étape de mise à jour des courtes longueurs d'ondes, et en utilisant les données de hautes fréquences et zéro ofset dans cette étape. Grâce à cette analyse, nous pouvons conclure que lorsque ces paramètres clés sont tous proprement appliqués à l'inversion, la mise à jour des grandes longueurs d'ondes est plus proche de la perturbation exacte du modèle que dans les cas où les paramètres sont partiellement appliqués. Du point de vue du coût de calcul, DWI est très cher, car l'inversion itérative est efectuée à chaque étape de la mise à jour des courtes longueurs d'ondes pour ajuster les données au lieu de juste imager les rélecteurs. Toutefois, cette étape itérative est importante, car elle facilite la mise à jour des grandes longueurs d'ondes du modèle. Cette conclusion est compatible avec l'importance d'introduire la migration itérative dans l'analyse de vitesse par migration (Migration Velocity Analysis, MVA) [Symes, 2015, Lameloise et al., 2015].

Introduction

In this chapter, I apply the classical FWI and the Decomposition-based Waveform Inversion (DWI) on 1D and 2D synthetic models. Recall that the decomposition-based waveform inversion consists of decomposing the gradient of FWI into a long-wavelength part and a short-wavelength part. The inversion is performed in an alternating fashion between the short-wavelength update and the long-wavelength update. For the numerical modeling, we use a inite-diference scheme with 8 th order in space and 4 th order in time. Absorbing boundary conditions are imposed at the borders of the model to simulate an ininite medium. We use the quasi-Newton optimization method, more precisely the L-BFGS with 5 iterations kept in memory.

We irst test FWI and DWI on a 1D layer model. Then we investigate a 2D synthetic model and compare the inversion results of FWI, multiscale FWI and DWI. We also analyse the inluence of key inversion parameters on the inversion result.

1D layer model with low frequency

We use the simple 1D layer model that has been analyzed in Chapter 3. Figure 4-1 shows the true velocity model. The model dimensions are 1.7 km large by 0.85 km deep.

The model consists of a strong velocity contrast, the velocities for the top layer, middle layer and bottom layer are 2 km/s, 3 km/s, 2 km/s, respectively. The initial model is a homogeneous model at 2 km/s, thus there is no a priori information in the initial model except the velocity at the surface. A Ricker wavelet with a central frequency of 12 Hz and temporal sampling of 1.5 ms is used as the source function. The sources and receivers are evenly distributed at the surface, with the source spacing being 17 m and receivers spacing being 8.5 m. Figure 4-2 on the top displays the result after running 80 iterations of conventional FWI with L-BFGS scheme. All ofsets available at the surface (up to 1.7 km) are used. It fails into a local minimum when starting from this homogeneous initial model. We note that it does not recover the long wavelengths of the model when all frequencies are inverted at the same time. Figure 4-2 in the second panel displays the inal result of DWI for 20 global iterations, and for each global iteration associated to the long-wavelength update, we run 20 iterations for the short-wavelength update using zero-ofset data. For the long-wavelength update, we use a mean ilter of size 9x9 grid points to remove the potential residual short-wavelength components. In this step, we use all ofsets available at the surface. Figure 4-2 in the third panel displays the result after 20 iterations of FWI with the result in the second panel as the initial model. We observe that DWI can better invert the long-wavelength components of the middle layer than the conventional FWI. Figure 4-2 on the bottom is the comparison of the vertical proile taken at x = 0.85 km. Compared to the conventional FWI, the decomposition-based waveform inversion better recovers the long-wavelength part of the middle layer, and as a result the bottom of the middle layer is better positioned. After the test on the simple 1D layer model, we investigate a more complex 2D synthetic model. The exact velocity model is shown in Figure 4-3. This model is inspired from [Perrone, 2013]. The dimensions of the model are 3.12 km by 1.212 km. The model consists of several layers and a velocity anomaly in the central part. Sources and receivers are evenly distributed at the surface, with source spacing being 48 m and receiver spacing being 12 m. The maximum ofset we use is 1.2 km in order to mainly simulate relected events. This maximum ofset is about the same dimension as the model depth. The density is considered constant with = 1000 kg/m 3 .

For the initial velocity model, we consider two models, one that is close, and the other that is far from the exact velocity model to evaluate convergence results in FWI and DWI. The one close to the exact model is a gradient-constant model (Figure 4-4 on the top), however it is not so close to the exact model, as there is no a prior information in the initial model. The one far from the exact model is a homogeneous model at 1.5 km/s (Figure 4-4 on the bottom). For the source wavelet, we use a source without low frequencies to be more realistic. In this case, the objective function of FWI is more nonlinear for data without low frequencies, and the inversion is more diicult to converge. However, we also test FWI using source wavelet with low frequencies as a reference. The source wavelet containing low frequencies is a Ricker wavelet with central frequency of 8 Hz (black line in Figure 45). The one without low frequencies is the same Ricker wavelet, but low frequencies lesser than 5 Hz are iltered (red line in Figure 45). The source wavelet is iltered from the initial Ricker wavelet with a high-pass Butterworth ilter of order 9 and with a cut frequency of 8 Hz . With these two initial models and two source wavelets we test the four combinations (Table 4.1). We start from the optimal case: good initial model and source wavelet with low frequencies, and gradually degrade the inversion conditions. In this case, the source wavelet and the initial model are favorable to the success of FWI.

The result of FWI after 100 iterations is shown in Figure 4-6 and 4-7. The comparison of the exact data, the initial data and inal data is displayed in Figure 4-8. We observe that the conventional FWI succeeds in this case. The result is considered as a reference for the inversion without low frequencies.

FWI

Firstly, we run FWI on this model. The FWI result after 100 iterations is shown in Figure 4-12. Compared to the irst test, as low frequencies are missing from the source wavelet, the problem becomes more nonlinear, and the FWI yields a poorer result, especially for the long-wavelength part in the central anomaly zone. In order to retrieve the long-wavelength components of the model, multi-scale FWI is widely used [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF]. We perform the multi-scale inversion to see if it could improve the inversion result. Multi-scale FWI

The multi-scale FWI is performed here in the time domain [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF], starting from the lowest available frequency band (5 Hz) and then gradually adding higher frequency bands. Here, we perform 4 FWI steps. The source wavelets are iltered from the initial Ricker wavelet with a Butterworth ilter of order 9 and in the frequency band of [8][9] Hz, [8][9][10][11][12] Hz, [8][9][10][11][12][13][14][15] Hz respectively (Figure 4-13). We perform 100 iterations of FWI at each step. The FWI result at each step is considered as the initial model for the next step. The results at each step are shown in Figure 4-14. Compared to the FWI result in Figure 4-12, the multi-scale FWI provides a better result and allows to better retrieve the main features of the model. However, the central velocity anomaly between 0.72 km and 0.9 km in depth is not well retrieved. As a result, two relectors at 0.72 km and 1.06 km are not positioned correctly due to the background velocity error.

DWI

We test the same model using DWI. As presented in Chapter 3, we use a nested algorithm. Within each global iteration, we perform 20 iterations of FWI using zero-ofset data and quasi-Newton (L-BFGS) algorithm for the model update [Nocedal, 1980]. The longwavelength update step is then performed using the gradient described in equation (3.4).

A single iteration of inversion is performed using the gradient descent algorithm [Sirgue,

FWI

Firstly, we apply the FWI on the initial model. The result after 100 iterations is shown in Figure 4-19. As low frequencies are absent in the source wavelet, and the initial model is too far from the exact model, the problem becomes quite nonlinear in this case, and the FWI falls into a local minimum . We perform the multi-scale inversion to see if it could improve the inversion result.

Multi scale FWI

The parameters for multi-scale FWI remain the same as in the third test, and the inversion results at each step are shown in Figure 4-20. Compared to the conventional FWI, multiscale FWI improves the long-wavelength part of the model (Figure 4-21). But with this initial model and this source wavelet, the multi-scale FWI is still unable to retrieve the long wavelengths of the velocity model and falls into a local minimum. The sign is opposite. As by using largeofset data, move-out information is introduced and long-wavelength components are also updated in the irst step, which is contrary to the idea of decoupling long wavelengths and short wavelengths, the inversion falls into a local minimum after the irst global iteration. Figure 4-28 shows the model perturbation after a single iteration of the long-wavelength update step of inversion when running 1 iteration (6 iterations for L-BFGS) instead of 15 iterations (20 iterations for L-BFGS) of zero-ofset FWI in the short-wavelength update step of inversion. Similarly, as shown in Figure 4-29, the result is far from the true model perturbation as iterations help to correct the relectivity. If just one iteration is performed, the amplitude of the relectivity is not correctly recovered and the amplitude error in the data residual will be interpreted as background velocity error in the second step. Figure 4-30 and 4-31 shows the model perturbation when full-spectrum source is used in the irst step. The image is close to the true model perturbation except on the edge. This is caused by the fact that low frequencies in the source will bring long wavelength update in the irst step.

Figure 4-32 shows the model perturbation when using a linearized inversion based on the Born approximation in the short-wavelength update step as [START_REF] Symes | Velocity inversion by diferential semblance optimization[END_REF] and [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] proposed. Note that the nonlinear inversion is performed in the second step as discussed previously. The poorer result (Figure 4-33) highlights the importance of the compatibility of waveform between these two steps even though the linearized inversion requires less iterations to converge and provides mainly short wavelength update. Figure 4-34 shows the comparison of a vertical slice taken on the relectivity models obtained by FWI and linearized inversion respectively, both using zeroofset data. We observe that the results are quite similar. But as the velocity contrast at interfaces is quite large, about 15% in this case, the Born approximation is no longer strictly valid. As a result, small diferences in relectivity could result in large diferences in the long-wavelength model update. If we use the Born approximation in the irst step, more iterations may be required to get a similar long-wavelength model.

Through these four comparisons we conclude that the three key ingredients: high frequency source, zero-ofset data and iterative inversion in the irst step, help to build a more correct short-wavelength relectivity which is favorable to the update of the longwavelength model in the second step. Besides, we note from the comparison of vertical proiles , that these key issues also help to better balance the amplitude for diferent parts of the model. Figure 4-35 shows the model perturbation when performing the irst step (shortwavelength update) with the high-frequency source and then using full-spectrum source in the second step (long-wavelength update). The similarity of the results shown in Figure 4-36, when compared to the one displayed in Figure 4-25 on the bottom, conirms that the proposed method is capable of retrieving long-wavelength components without low frequency in the source. This property is similar to that of DSO. As shown by [van Leeuwen and Mulder, 2010], the DSO functional does not depend on the frequency content of the data.

There are two possible interpretations to explain the success of the decomposed gradi- ent to reduce the nonlinearity of the inversion problem. The irst one is that the decomposed gradient is a iltered version of the original gradient of FWI. The ilter serves as a projector of the original gradient direction into the direction of the true model perturbation. The second interpretation is that since we use only one part of the FWI gradient, the objective function in equation ( 2.1) is no longer valid. The new objective function corresponding to the new gradient is less nonlinear and has a larger basin of attraction.

Conclusion and discussion

In this chapter, we apply the FWI and DWI on two synthetic models and show that compared to the conventional FWI and the multi-scale FWI, DWI is more robust with respect to the absence of low frequencies in the source wavelet and to the large errors in the initial model. We have also analyzed the inluence of three key parameters during the inversion and show that when these three key parameters are all applied in the inversion, the update of the long-wavelength part is closer to the true velocity perturbation.

From the computation cost point of view, the new method is computational expensive, as iterative inversion is performed at the short-wavelength update step to it the data. However this iterative step is important, since it facilitates the long-wavelength update of the model. For migration velocity analysis, [Symes, 2015, Lameloise et al., 2015] also underline the importance of iterative migration to it the data in the short-wavelength update, rather than merely image the relectivity. A possible remedy to reduce the computation time is to perform the iterative inversion in the pseudo-time domain to avoid the re-computation of the relectivity at each step [Plessix, 2012, Brossier et al., 2015].

We have also tested the idea that in the second step of inversion, we do not decompose waveields, but directly smooth the FWI gradient to update the long-wavelength part of the model. However, this method fails to give a satisfactory result. Our analysis is that the limit between the short wavelengths and the long wavelengths is not ixed at each iteration, their spectrum may overlap or vary during the iterations, thus it is not trivial to determine the limit at each iteration.

Finally, all the examples shown in this chapter use models with constant and known density. We have also tested the more realistic case with the exact model having variable density, and the model used for the inversion having constant density. In this case, similarly as for conventional FWI, DWI fails to retrieve the correct model. For future work, we would like to extend the method to be suitable for elastic models with variable density. This aspect is described in the perspectives (Chapter 6).

Résumé du chapitre 5

Dans ce chapitre, nous appliquons la FWI classique et la DWI au jeu de données réelles 2D Brunei fourni par Total. L'acquisition des données Brunei a été réalisée en 2011, et la géométrie d'acquisition est montrée dans la Figure 5-1. Nous présentons d'abord les paramètres d'acquisition et les traitements de données que nous efectuons, y compris la vériication de la cohérence entre traces, le pointé des premières rélexions, la suppression des ondes directes, le iltrage passe-bas des données, le ré-échantillonnage et le débruitage des données. En particulier, nous testons et comparons trois méthodes de pointé et observons que celle proposée par [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF] est la plus eicace.

Nous estimons la source avec la méthode proposée par [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF]. La bande de fréquence estimée pour la source correspond à la bande de fréquence des données. En théorie, l'ondelette de la source doit être estimé à chaque itération de l'inversion, en même temps que la mise à jour du modèle de vitesse. Toutefois, pour simpliier, l'ondelette de la source que nous utilisons reste inchangée au cours de l'inversion. Néanmoins, cette hypothèse pourrait pénaliser le résultat d'inversion.

Nous soulignons aussi les déis que nous avons rencontrés. Tout d'abord, nous étudions le cas de la densité constante. Nous observons que avec l'hypothèse de densité constante la variation de l'amplitude avec ofset difère considérablement entre les données observées et calculées. En outre, pour les grands ofsets, nous pouvons aussi constater un déphasage entre ces deux données. Comme nous utilisons l'approximation acoustique de la propagation des ondes avec densité constante, l'inluence de la densité et d'autres paramètres physiques du sous-sol sur l'amplitude des données de rélexion ne pouvait pas être incorporé correctement. Notons que pour des ondes plongeantes, l'inluence de la variation de densité n'a pas beaucoup d'efet sur l'amplitude des données, mais que pour les rélexions, cet efet est important. Nous essayons d'utiliser un modèle de densité variable pour compenser la diférence d'amplitude. Nous testons la formule de Gardner et observons que l'amplitude est mieux expliquée que dans le cas de densité constante. Mais le problème de déphasage aux grands ofsets existe toujours. Il semble donc diicile de trouver un modèle de densité qui pourrait mieux expliquer la variation d'amplitude avec ofset. Nous pensons que l'atténuation est un autre facteur important, surtout dans les zones supericielles en présence de gaz. De ce fait, nous avons décidé de prendre en compte les ofsets limités dans les données.

Le deuxième test préparatoire nous efectuons est d'utiliser le modèle de vitesse fourni par Total pour faire la migration. Nous observons que l'image migrée dans ce modèle de vitesse n'est pas focalisée. Si l'on migre avec les très courts ofsets, les images sont bien focalisées mais difèrent les unes des autres et la somme des contributions n'est plus cohérente, surtout dans les zones de gaz.

Avec les diicultés que nous avons rencontrées, nous montrons quand même les résultats d'inversion en utilisons FWI classique et DWI. D'après le résultat de FWI, nous observons que la partie des grandes longueurs d'ondes du modèle n'est pas mis à jour. En outre, en raison de l'hypothèse de densité constante, la mise à jour de la densité se retrouve dans la mise à jour de vitesse, visible à la fond de l'eau. Pour DWI, nous observons que la poche de gaz est mieux récupérée que dans FWI. Mais le problème de l'amplitude existe toujours, qui est aussi visible sur le fond de l'eau. D'autres études sur ce sujet vont certainement aborder l'utilisation du modèle visco-élastique. En outre, plus de travail serait nécessaire

Introduction

In this chapter, we present an application on the 2D Brunei real dataset as provided by Total. The acquisition of the Brunei data set was carried out in 2011, and the acquisition geometry is displayed in Figure 5-1. We irst present the seismic acquisition parameters and data processings we performed. We point out the challenges we met when dealing with this dataset. Then we show the irst inversion results using the conventional Full Waveform Inversion (FWI) and Decomposition-based Waveform Inversion (DWI). More work is still needed to correctly image the subsurface.

Seismic dataset

The seismic dataset consists of a sub-line (that mimics a 2D acquisition) with 608 shot gathers with sources spaced 50 m apart and receivers spaced 12.5 m apart. Within each shot gather, there are 516 channels. The irst channel is situated at 150 m from the source and the last channel is situated at 6587.5 m. The temporal sampling is 3 ms and the total recording time is 9 s. Several preprocessings were already done by the provider of the data:

• Multiple removal using 3D SRME and Radon-based methods.

• Butterworth high-pass ilter with low cut frequency of 3 Hz.

• Swell noise and linear noise attenuation.

Figure 5-2 displays one shot gather after the preprocessings described above, and the spectrum of the data corresponding to this shot gather. To meet our needs, we have also applied other preprocessings, including:

• Seismic event coherency check. In this step, we check the quality of each shot gather and noisy traces are dropped.

• Picking and direct wave removal. The direct waves are not completely removed from the initial data, as shown in Figure 5-2. As we do not use direct waves for the inversion, we apply a picking on the data and remove the direct waves. For the irst-relection picking we have irst tested two methods that are classically used for picking of the irst-breaks: 1) the method proposed by [Coppens, 1985], that consists in calculating the ratio of energy of seismogram for two windows, one before the current time point, and one after the current time point. However, with this dataset, this approach fails to detect the irst-relection arrival because of the presence of the residual direct wave (Figure 5-3); and 2) the STA/LTA (Short Term Average / Long Term Average) ratio method [Allen, 1978], that consists in calculating the ratio of two averages of energy between a short-term window and a long-term window. The picking results using diferent sizes of windows are shown in Figure 5-4. It fails to pick the irst relection and indicates the second relection. We also tested an algorithm that is widely used in image processing, but less standard in geophysics. The method is based on the contour detection [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF].

It is similar to Coppens's algorithm. This approach calculates the absolute value of the diference, instead of the ratio of energy of seismogram for two windows. This alternative improves the picking result, as the diference operator is more stable than the ratio operator. The picking result is shown in Figure 5-5. It successfully picks the irst relection. The result after the muting of the direct wave is shown in Figure 5-6 on the left.

• Low pass iltering of the data. As shown in Figure 5-2 on the right, the spectrum of the original data is up to 100 Hz. This frequency is very high for FWI. Moreover, with this frequency, the spatial and the temporal samplings for the wave propagation modeling would be very small and the size of model would exceed our computational capacity. Thus, we apply a low pass ilter on the data, by using a 10 th order Butterworth ilter with high cut frequency of 10 Hz. The iltered shot gather is shown in Figure 5-6 on the right.

• Data resampling. The data is resampled from 3 ms to 2 ms in the time domain, and from 25 m to 12.5 m in the space domain, to satisfy the criteria for the spatial and the temporal samplings for the waveield propagation.

• Data de-noising. In the lower part of the shot gather (Figure 5-6), we can observe some artefacts that intersects with each other and have slopes diferent from the slopes of the signal. We ilter these artefacts in the f-k domain. The denoised shot gather and its spectrum are displayed in Figure 5-7.

Note that there is no 3D to 2D amplitude correction in the initial data. This correction consists in multiplying the data by √ , where is the recorded time. Since there is a lot of noise in the bottom of the data, rather than amplifying later arrivals in the observed data, we instead correct the calculated data by dividing it by √ . This behaves as a preconditioning on the data.

Seismic inversion

A FWI model is provided by Total using long-ofset data (Figure 5678). The dimensions of the model are 16.7625 km large by 3.0125 km deep. The velocity model building has to contend with challenging geology, including steeply dipping structure and the presence of large shallow gas pockets. These gas pockets exhibited extremely low seismic velocities, and had a negative impact on deeper relector amplitudes. The objective is to localize these gas traps. When we began to deal with this data set, one of the irst tests we wanted to perform was to compare the data modeled on this FWI model with the processed observed data. For this, we need irst to estimate the source wavelet.

Estimation of the source wavelet

We use the source-estimation method proposed by [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF], and described in Chapter 2. This method was initially applied on direct waves, but as the direct waves are partially removed in our initial data, the residual direct waves are not suicient for the source estimation (Figure 5-2). Instead we use the irst relection at the water bottom to estimate the source wavelet. We have tested this approach at three source locations: at 2.5 km, 7.5 km and 10 km respectively. The estimated source wavelets are not exactly the same for these three locations. However, we decided to use an average source wavelet.

The averaged source wavelet and its spectrum are displayed in Figure 5-9. The bandwidth of the estimated source corresponds to the bandwidth of the data (Figure 567). In theory, the source wavelet should be estimated at each iteration of the inversion, along with the velocity update. However, for simplicity, the source wavelet we use remains unchanged during the inversion. This may penalize the inversion result.

Preparatory tests

With the estimated source wavelet, we compute the data using the model in Figure 5-8. First, we investigate the constant density case. Let consider the shot gather for the source situated at 7.5 km. The comparison of the observed data and the calculated data is displayed in Figures 5-10 and 5-11. We note that the variation of the amplitude with ofset difers greatly between these two data. Besides, for large ofset, the phase is also shifted. As we use the acoustic wave propagation approximation with constant density, the inluence of density and other physical parameters of the subsurface on the relection data amplitude could not be incorporated correctly. Note that for diving waves, the inluence of the variation of density does not have much efect on data amplitude, but for relections, this efect is important. We try to use a variable density model to compensate the amplitude diference. We rely on the Gardner's formula [START_REF] Gardner | Formation velocity and density-the diagnostic basics for stratigraphic traps[END_REF]]:

= β .
(5.1) The Gardner's formula is an empirically derived equation that relates seismic P-wave velocity to the bulk density of the lithology. We tested diferent choices of and and ind none of them could match the data amplitudes satisfactorily. We cite one example here. To compensate the amplitude diference, a large density contrast at the water-bottom layer is needed. The velocity of the water here is 1495 m/s. We choose = 7.74 -12 and = 3.5. The obtained density model is displayed in Figure 5-12. The source wavelet is re-estimated. We then recompute the shot gather for the source situated at 7.5 km, using the velocity model in Figure 5-8 and the density model in Figure 5-12. The comparison of the calculated and the observed data is shown in Figure 5-13 and 5-14. We observe that the amplitude is better matched than in the constant density case. But the phase shift problem still exists. Besides, due to the density change, the amplitude of the events below the irst relection is attenuated; refer for example to the event at =5 s and =6 km. It seems therefore diicult to ind a density model that could better explain the real data amplitudes with ofset variations.

Offset (km) Time (s)
As the decomposition-based waveform inversion does not need large ofset data, we decide to only consider the data in the ofset range of m to avoid the large variations in amplitude and phase.

The second preparatory test we carry out is to use the velocity model in Figure 5-8 to get a migrated image. We use the constant density model. The migrated image is displayed in Figure 5-15. We note that this image is not focused, especially in the red circled zone. To analyse this problem, we use data with smaller ofset ranges. We test 5 ofset ranges: respectively [150-275] m, [275-400] m, [400-525] m, [525-650] m, and [650-775] m. For each ofset range, the migrated image is focused . But if we compare carefully, we ind that these images are shifted horizontally (indicated by red vertical lines in Figure 5-16). Due to this shift, the sum of the ive migrated images, providing the migrated image for data with ofset range of [150-775] m, is still unfocused at = 6 km (see Figure 5-17). Three vertical slices at = 3 km, 6 km, 10 km respectively in these ive gradients are shown in Figure 5-18. We can observe the moveout at diferent ofsets. For the following inversion tests, we only consider the right part of the model, starting from = 8 km, to avoid the diicult zone near = 6 km.

Seismic inversion

In this section, we apply FWI and DWI on the data. We consider the constant density case. The initial model (Figure 5-19) we use is the smoothed version of the FWI model in Figure 5-8.

FWI

We apply the conventional FWI on the data in the time domain, with the estimated source wavelet (Figure 56789). All frequencies are involved at the same time. The sources are spaced 100 m apart and the receivers are spaced 12.5 m apart. The FWI result after 100 iterations is displayed in Figure 5-20. We note that the long-wavelength part of the model is not updated. Besides, due to the use of constant density, the density update is leaked into the velocity update (see the bottom of the water layer).

DWI

We apply the DWI on the same data, with the same source wavelet. We use a nested algorithm. For the short-wavelength update step, we perform 20 iterations of conventional FWI with short-ofset data. The ofset range is [150-400] m. For the long-wavelength update, we decompose the waveield and only consider the long-wavelength part in the gradient. We perform 3 global iterations. The long-wavelength gradient at the irst iteration is displayed in Figure 5-21. We compare it with the true velocity perturbation in Figure 5-22. The "true model" here means the FWI model in Figure 5-8. The true velocity perturbation thus is the diference between the FWI model and the initial model . We can observe that the gas pocket at = 2 km is better retrieved than in the FWI case. But the problem of density still exists as visible below the water-bottom layer. The long-wavelength velocity model after 3 global iterations is displayed in Figure 5-23 on the top, and the FWI result after 100 iterations, starting from this model, is displayed in Figure 5-23 on the bottom. The result is not yet satisfactory. Compared to conventional FWI, one improvement may be found in the gas pocket at = 2 km, however, the bottom of the water layer pose more problems in this case.

Discussion and conclusion

In this chapter, we apply the FWI and DWI on the Brunei real data set and show that compared to the conventional FWI, DWI better recovers the gas pocket. However, due to the acoustic approximation with constant density, the data is not correctly simulated, especially for dealing with the data amplitude at the water-bottom layer (Figure 5-24).

Attenuation with low quality factor is certainly visible due to the presence of gas pockets [START_REF] Operto | A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice[END_REF]. Further investigations are needed to incorporate elastic, density and attenuation efect. Besides, we had wanted to consider more frequencies in the data to get a more focused relectivity in the short-update step of the inversion, which is favorable for the long-wavelength update of the model. However, adding higher frequencies in the data also means a iner spatial and temporal sampling, yielding larger model which is beyond our actual computational capacity.

Chapter 6

Conclusions and perspectives

Résumé du chapitre 6

Ce chapitre est consacré aux conclusions et aux perspectives de la thèse. L'inversion des formes d'ondes (Full Waveform Inversion, FWI) est largement utilisé pour l'imagerie sismique. Cependant, la fonction objective de FWI est fortement non linéaire par rapport aux paramètres du modèle et présente beaucoup de minima locaux [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF]. Ces minima locaux empêchent les techniques basées sur le gradient de trouver le minimum global si le modèle initial est loin de la solution globale. Si un ensemble de données contient principalement des rélexions, les méthodes basées sur le gradient souvent n'arrivent pas à reconstruire les grandes longueurs d'ondes du modèle [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF]. Dans cette thèse, nous avons présenté et testé une nouvelle méthode d'inversion basée sur la séparation du gradient après la décomposition du champ d'onde (Decomposition-based Waveform Inversion, DWI). L'objectif est de reconstruire les grandes longueurs d'ondes du modèle de vitesse. La recherche est principalement axée sur trois aspects, comme décrit dans la suite.

Dans le chapitre 2, nous avons revu les formulations classiques de FWI, y compris la fonction objective, les méthodes de modélisation numérique pour la propagation des ondes, le calcul du gradient et du Hessien de la fonction objective, et les méthodes de mise à jour du modèle de vitesse. L'analyse de la résolution du gradient de FWI indique que la résolution dépend des diférents types d'ondes et de la coniguration d'acquisition. Les ondes plongeantes, les transmissions liée aux rélexions et aux réfractons, et les multiples particulièrement aident à reconstruire les grandes longueurs d'ondes du modèle. Les inversions avec diférents paramètres d'acquisition indiquent que les données de grands ofsets et de basses fréquences sont favorables pour la mise à jour des grandes longueurs d'ondes du modèle. Les analyses efectuées dans ce chapitre et la littérature nous ont amenés à proposer la nouvelle méthode d'inversion. Pour mieux reconstruire les grandes longueurs d'ondes du modèle, dans le chapitre 3, nous avons proposé une méthode d'inversion basée sur la décomposition. Cette approche consiste à séparer le gradient de FWI en une partie de grandes longueurs d'ondes et une partie de courtes longueurs d'ondes, puis de les inverser séparément et alternativement. Cette séparation du gradient est obtenue par la décomposition des champs d'ondes en une partie descendante et une partie montante. Au lieu de mettre des diférents poids dans le gradient [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF], où le poids est diicile à déterminer, nous proposons d'inverser ces deux parties séparément. L'inversion est efectuée avec un algorithme imbriqué. La méthode est illustrée au travers d'un modèle de couche 1D.

Nous avons également comparé DWI avec d'autres méthodes d'inversion basée sur la rélexion. Nous comparons en particulier avec la méthode proposée par [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF]. La comparaison à travers le modèle de couche 1D indique que les mises à jour du gradient pour ces deux méthodes sont très similaires. La diférence réside dans la façon de calculer les champs d'ondes partiels. En outre, comme la méthode de [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] est basée sur la modélisation linéarisé, elle suppose que les perturbations de vitesse sont petites et ne considère que les rélexions primaires. En revanche, comme la FWI classique, DWI résout l'équation des ondes non linéarisée, donc les autres types d'ondes pourraient être pris en considération. Comparée à d'autres méthodes qui iltrent le gradient pour obtenir les grandes longueurs d'ondes, DWI récupère les grandes longueurs d'ondes automatiquement. Ceci est un avantage quand on connaît peu sur le modèle. Une autre méthode similaire est proposé par [START_REF] Zhou | Combining diving and relected waves for velocity model building by waveform inversion[END_REF]. Cette approche consiste à combiner la partie tomographique dans le gradient de FWI obtenu par les rélexions avec le gradient obtenu par les ondes plongeantes pour extraire les grandes longueurs d'ondes du modèle. Cette méthode présente l'avantage d'être plus rapide que DWI, mais elle est plus sensible à la rélectivité. Elle fonctionne de manière optimale pour la rélectivité avec une fréquence ininie, mais en pratique, en raison des données de bande limitée, la rélectivité est également limité en fréquence.

La décomposition du gradient est obtenue par la décomposition du champ d'ondes. Nous avons présenté et comparé trois méthodes pour décomposer le champ d'ondes. Ces méthodes sont illustrées à travers un champ d'ondes qui se propage dans le modèle de couche 1D. Les avantages et les inconvénients de chaque méthode sont également détaillée. La méthode par la transformation de Fourier est la plus populaire, et a été largement utilisée pour supprimer des artefacts de haute fréquence dans les images migrées, mais avec cette méthode les basses fréquences sont iltrées dans le résultat inal. La méthode de Poynting est la plus rapide car elle peut être calculée au cours de la propagation des ondes, mais elle soufre du problème d'artefacts numériques forts. La méthode de décomposition par les curvelet donne le meilleur résultat car il peut prendre en compte les diférentes directions en même temps, avec le coût de calcul le plus élevé des trois. Cependant, contrairement à la méthode de Fourier qui nécessite de garder en mémoire l'ensemble des champs d'ondes, la méthode de curvelet est facilement parallélisable parce que la décomposition est efectuée séparément pour chaque pas de temps.

Dans le chapitre 4, nous appliquons FWI, FWI multi-échelles et DWI sur deux modèles synthétiques et montrons que par rapport à la FWI conventionnel et la FWI multi-échelle, la DWI est plus robuste face à l'absence de basses fréquences dans la source et à des erreurs importantes dans le modèle initial. Nous avons également analysé l'inluence de plusieurs paramètres clés lors de l'inversion, comme l'application de l'inversion itérative et non linéaire en utilisant des données de hautes fréquences et zéro ofset dans l'étape de mise à jour des courtes longueurs d'ondes. Grâce à cette analyse, nous pouvons conclure que lorsque ces paramètres clés sont tous appliqués dans l'inversion, la mise à jour des grandes longueurs d'ondes est plus proche de la vraie perturbation du modèle de vitesse que dans les cas où les paramètres clés sont partiellement appliqués.

Dans le chapitre 5, nous avons appliqué FWI et DWI sur un jeu de données réelles de Brunei. Nous montrons que, avec l'hypothèse acoustique, densité constante et milieu non atténuant, il est diicile de simuler les formes d'onde dans ce jeu de données réelles, en particulier en ce qui concerne le comportement de l'amplitude avec ofset. Avec ces hypothèses limitées, ni FWI ni DWI peuvent fournir un résultat satisfaisant. Cependant, DWI donne un résultat légèrement meilleur pour retrouver les poches de gaz.

A part le problème de l'amplitude, pour DWI, la variation de densité inluence également la rélectivité, ce qui est un facteur important pour DWI. Vue dans le chapitre 4, une petite diférence dans la rélectivité peut provoquer une grande diférence dans le résultat de l'inversion. D'autres investigations sur ce sujet vont certainement aborder l'utilisation du modèle visco-élastique. En outre, des recherches supplémentaires seraient nécessaires pour analyser l'inversion multi-paramétrique avec la méthode DWI. Modèles synthétiques doivent être validés avant de traiter les données réelles.

Pour DWI, l'inversion itérative est efectuée pour la mise à jour des courtes longueurs d'onde pour ajuster les données. Cette étape est très coûteuse en temps de calcul. Toutefois, cette étape itérative est importante car elle facilite la mise à jour des grandes longueurs d'onde du modèle. Cette conclusion est compatible avec l'importance d'introduire la migration itérative dans l'analyse de vitesse par migration (Migration-based Velocity Analysis, MVA) [Symes, 2015, Lameloise et al., 2015]. Pour les futures études, nous aimerions mettre en oeuvre cette étape dans le domaine pseudo-temps pour éviter la calcul à chaque itération globale [Plessix, 2012, Brossier et al., 2015]. Notons que pour DWI, la longueur du pas pour la mise à jour des grandes longueurs d'ondes est déinie au début. Ainsi, la décroissance de la fonction objective n'est pas assurée au cours des itérations. Si nous utilisons l'interpolation quadratique pour déterminer la longueur du pas à chaque itération, nous avons besoin d'efectuer deux problèmes directs supplémentaires, et ici le problème direct correspond à l'inversion itérative, donc rend le calcul très cher. Si cette étape d'inversion itérative pourrait être réalisée dans le domaine pseudo-temps, le coût de calcul pourrait être largement réduit.

Par ailleurs, dans cette thèse, DWI principalement utilise les rélexions. Cependant, en théorie, DWI n'est pas limitée aux rélexions, ce qui est diférent des autres méthodes d'inversion basées sur la rélexion, et pourrait intégrer d'autres types d'ondes en même temps. Pour les futures études, nous tenons à combiner ces diférents types d'ondes et voir si la combinaison pourrait faciliter l'inversion. [START_REF] Biondi | Tomographic full waveform inversion (tfwi) by combining full waveform inversion with wave-equation migration velocity anaylisis[END_REF] introduisent l'ofset en sous-sol pour atténuer la non linéarité du problème. Nous aimerions étudier de plus près sur ce point, en introduisant l'ofset en sous-sol pour DWI.

Ain d'atténuer la grande dépendance vis-à-vis de l'amplitude des données, nous aimerions tester d'autres fonctions objectives, qui sont moins sensibles à l'amplitude, telles que la fonction objective basée sur la corrélation [van Leeuwen and Mulder, 2008] et fonction objective normalisée.

Enin, avec le développement des ordinateurs de haute performance, l'application de la méthode de décomposition basée sur le curvelet aux données réelles pourrait devenir bientôt possible. L'élément important est la corrélation dans le domaine des curvelets et non pas en temps pour chacun des directions.

Conclusions

Full Waveform Inversion (FWI) is widely used for seismic imaging. However, the objective function of FWI is highly nonlinear with respect to the model parameters and has many local minima [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF]. These local minima prevent gradient-based techniques from inding the global minimum if the initial model is far from the global solution. If a dataset mainly contains relections, the gradient-based methods does not recover the long wavelengths of the model [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by waveield inversion[END_REF]. In this thesis, I have presented and tested a new inversion method based on the gradient separation after waveield decomposition.

The objective is to retrieve the long wavelengths of the velocity model. The research mainly focused on three aspects, as described in the following.

6.1.1 Basic FWI theory and resolution analysis

In Chapter 2, we have reviewed the classical formulations of FWI, including the objective function, the methods for the numerical modeling of waveield propagation, the computation of the gradient and the Hessian of the objective function, and the velocity model update methods. The resolution analysis of the FWI gradient indicates that the resolution depends on the types of waves and on the acquisition setup. The diving waves, relection-related transmission, refraction-related transmission and multiples particularly help to recover the long-wavelength part of the model. The inversions with diferent acquisition parameters indicate that long-ofset and low-frequency data are favorable for long-wavelength update of the model. The analyses performed in this chapter and the revised literature brought us to propose the DWI inversion method.

6.1.2

The DWI method

To better recover the long-wavelength part of the model, in Chapter 3, we have proposed a Decomposition-based Waveform Inversion (DWI) method. The method consists in separating the gradient of FWI into a long-wavelength part and a short-wavelength part, and then inverting them separately. This separation of gradient is based on the decomposition of the waveield into their down-going and up-going part. By introducing this waveield decomposition into the gradient of FWI, a recombination of them gives the long-and short-wavelength parts of the gradient. Instead of setting diferent weights ( [START_REF] Tang | Tomographically enhanced full waveield inversion[END_REF]), for which the weight is diicult to determine, we propose to separately invert these two parts. The inversion is performed using a nested algorithm. The method is illustrated through a 1D layer model. The decomposed gradients using the 1D layer model show that one part of the gradient of FWI indicates the positions of relectors and thus provides the short wavelengths of the model, and the other part indicates the wavepath that connecting the surface and the relector and provides the long wavelengths of the model.

We have also compared DWI with other relection-based waveform inversion methods. We compare particularly with the method proposed by [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF]. The comparison through the 1D layer model indicates that the gradient updates for these two methods are very similar. The diference resides in the way to compute the partial waveields. Besides, as the method of [START_REF] Xu | Full waveform inversion for relected seismic data[END_REF] is based on the linearized modeling, it assumes that the velocity perturbations are small and the data contain only primaries. While similarly to classical FWI, DWI resolves the nonlinear wave equation; therefore refractions, diving waves and multiples could be included in the dataset. Compared to others methods that ilter the gradient to get the long wavelengths, DWI retrieves the long-wavelength part automatically. This is an advantage when we know little about the model.

Another similar method is proposed by [START_REF] Zhou | Combining diving and relected waves for velocity model building by waveform inversion[END_REF]. This method consists in combining the tomographic part in the FWI gradient using relection data with the gradient using diving waves to retrieve the long wavelengths of the model. The tomographic part in the FWI gradient with relection is the same as we use. Instead of decomposing waveields, the authors irst compute the gradient of FWI in the model with relectivity 0 + , and then they compute the gradient of FWI in the smooth model 0 , and the diference of these two cancels the migration part in the FWI gradient and provides the tomographic part. This method has the advantage of being faster than the waveield decomposition-based method, however it is more sensitive to the relectivity. It works optimally for relectivity with ininite frequency, but in practice, due to the band-limited data, the relectivity is also band-limited. This band-limited relectivity eventually change the kinematics of the model. Thus, the migration part in the FWI gradient calculated in the full model and in the smooth model may be diferent and the diference of these two gradients may contain residual short wavelengths.

The decomposition of the gradient is based on the waveield decomposition. We have presented and compared three methods to decompose the waveield into their one way component. These methods are illustrated through a waveield that propagates in the 1D layer model. The advantages and the disadvantages of each method are also detailed. The Fourier transform method is the most popular one, and was widely used for suppressing high frequency artefacts in migrated images, but low frequencies are iltered in the inal results. The Poynting method is the fastest one as it could be computed during the wave propagation, but it sufers from strong numerical artefacts. The curvelet decomposition method yields the best result as it can take into account diferent directions at the same time, with the expense of the highest computational cost of the three. However, diferently from the 2D Fourier transform method which requires to keep in memory the whole propagated waveield, the curvelet decomposition method is more easily parallelizable because the decomposition is carried out for each snapshot separately.

6.1.3

Applications to synthetic models and real data

In chapter 4, we apply FWI, multiscale FWI and DWI on two synthetic models and show that compared to the conventional FWI and the multiscale FWI, DWI is more robust with respect to the absence of low frequencies in the source wavelet and to large errors in the initial model. We have also analyzed the inluence of several key parameters during the inversion, such as performing the nonlinear iterative inversion in the short-wavelength update step, and using high frequency and zero-ofset data in the short-wavelength update step. Thanks to this analysis, we can conclude that when these key parameters are all applied in the inversion, the update of the long-wavelength part is closer to the true velocity perturbation than in cases where the key parameters are partly applied. From the point of view of computational cost, DWI is computational expensive, as iterative inversion is performed at the short-wavelength update step to it the data. However this iterative step is important, since it facilitates the long-wavelength update of the model. This conclusion is consistent with introducing iterative migration in MVA [Symes, 2015, Lameloise et al., 2015].

In chapter 5, we have applied FWI and DWI on the Brunei real data set. We show that with the acoustic and the constant density and non attenuative assumption, it is diicult to simulate the waveforms present in this real dataset, especially regarding the behavior of the amplitude with ofset. With this limiting assumption, neither FWI nor DWI could provide a satisfactory result. However, DWI yields a slightly better result to retrieve the gas pockets. For the wave propagation, we assume in this thesis the acoustic approximation, with constant density and no attenuation. This assumption is very limiting when dealing with real datasets, as they are better interpreted by visco-elastic modeling. The inadequacy of the acoustic assumption with constant density and without attenuation to interpret real data and especially to interpret the data amplitude, is observed in Chapter 5. We think that attenuation is a major factor. Further investigations on this topic will certainly address the use of visco-elastic model. Moreover, further research would be needed to analyze multiparameter inversion with the DWI method. Synthetic models need to be validated before dealing with real data.

Role of density

We have investigated the inluence of density on the data amplitude in Chapter 5. Apart from the data amplitude issue, for DWI, the variation of density also inluence the relectivity, which is an important factor for DWI. We analyze this efect through a 1D simple model (Figure 6-1). The true velocity for the irst and the third layer is 2 km/s and for the second layer is 2.5 km/s. The true density for the irst and the third layer is 1000 kg/m 3 and for the second layer is 1500 kg/m 3 . The initial velocity is homogeneous at 2 km/s. We have tested three density models. The irst one is the true density model. The second one is the homogeneous model at 1000 kg/m 3 , and for the third model, the density for the second layer is 2000 kg/m 3 and for the rest of the model is 1000 kg/m 3 . The result after running 20 FWI iterations with ixed density parameters are shown in Figure 6-2. We observe that the density error is leaked into the velocity update. For the initial density that is smaller than the true density, the velocity is overestimated to compensate the density error, and for the initial density that is larger than the true density, the velocity is underestimated to correct the overestimation of density. Besides, the sign of the relectivity is changed. As we have analyzed in chapter 4, DWI is quite sensitive to the relectivity. A small diference in the relectivity may cause a large diference in the inversion result. The inluence of the density on the inversion, especially on DWI needs further investigations.

decomposed waveields in the curvelet domain before the inverse curvelet transform to compute the gradient, rather than performing the inverse curvelet transform on each snapshot and then correlating the whole waveield in the space domain.

  . . . . . . . . . . . . . . . . . . . 2-1 An example of seismic data. We can observe the direct wave and the relected waves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-2 The original Ricker wavelet (black) and the estimated source (red) in the frequency domain (left) and in the time domain (right). . . . . . . . . . . . 2-3 Wavenumber illumination. One source-receiver pair, one scattering point and one frequency in the data provides one wavenumber in the image space [Huang and Schuster, 2014]. . . . . . . . . . . . . . . . . . . . . . . . 2-4 Approximate resolution limits for (a) diving waves, (b) and (c) relectionrelated transmission, (d) difraction-related transmission and (e) relection migration [Huang and Schuster, 2014]. . . . . . . . . . . . . . . . . . . . . 2-5 First Fresnel zone for the specular relection

  ]. 2-7 a) Ray diagram for interbed multiples generated by a difractor in a thin layer and b) the associated mirror sources diagram [Huang and Schuster, 2014]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-8 1D exact velocity model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-9 Four Ricker wavelets. From top to bottom and from left to right, the central frequency is respectively 4 Hz, 8 Hz, 12 Hz, and 16 Hz. . . . . . . . . . . . 2-10 (Left) Gradient of the irst iteration of FWI with source wavelets displayed in Figure 2-9. The central frequency is respectively 4 Hz, 8 Hz, 12 Hz, and 16 Hz. (Right) The bandwidth of the vertical slice taken at the middle of each model on the left. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-11 (Left) Gradient of the irst iteration of FWI for data with ofset range [0-0.17] km, [0.34-0.51] km, [0.68-0.85] km respectively. (Right) The bandwidth of the vertical slice taken at the middle of each model on the left. . . 2-12 Comparison of the convergence rate of local methods applied on a 2D linear inverse problem. The problem is solved using a) Gauss-Newton, b) gradient descent, c) preconditioned gradient and d) conjugate-gradient methods. The gradient in c) was preconditioned with the diagonal term of the inverse Hessian [Sirgue, 2003]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-9 Direct and backpropagated waveield in the model in Figure 2-8. . . . . . . 3-10 Decomposed waveields for the initial waveields in Figure 3-9. + (top left), + (top right), -(bottom left) and -(bottom right). . . . . . . . 3-11 Four waveields in equation (3.11). 0 (top left), (top right), (bottom left) and 0 (bottom right). . . . . . . . . . . . . . . . . . . . . . . . . . .

  -8 at = 0.525 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-14 (Left) The section of the waveield in thedomain for = 0.85 km in the middle of the model. (Right) FT of the igure on the left in thez domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-15 The down-going and up-going ilter in thez domain. . . . . . . . . . . 3-16 Filtered parts of the initial image in Figure 3-14. . . . . . . . . . . . . . . . 3-17 Decomposed waveields in thedomain. Notice the diferent color-scale between the two images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-18 Decomposed waveields in the space domain using the 2D FT-based decomposition. + (left) and -(right). . . . . . . . . . . . . . . . . . . . . . . 3-19 (Left) Partial derivative of the waveield with respect to . (Right) Partial derivative of the waveield with the respect to . . . . . . . . . . . . . . . . 3-20 The product of the two images in Figure 3-19. . . . . . . . . . . . . . . . . 3-21 The decomposed waveields using Poynting-vector method. + (left) and -(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-22 Original mask (left) and median-iltered mask (right) for the decomposition by the Poynting-vector method. . . . . . . . . . . . . . . . . . . . . . . . . 3-23 Decomposed waveields using the mask in Figure 3-22 on the right. + (left) and -(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-24 A scale example of curvelet

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-27 Two snapshots at diferent recording time for = 0.525 s (left) and = 0.54 s (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-28 One direction after the curvelet transform for the two snapshots in Figure 3-27. 3-29 Decomposed waveields using the curvelet method. + (left) and -(right). 3-30 One snapshot of the waveield propagated in the Marmousi model. . . . . . 3-31 Decomposed waveields using the 2D FT method (irst panel), Poynting method (second panel) and curvelet method (third panel). + (left) and -(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-1 1D layer exact velocity model. . . . . . . . . . . . . . . . . . . . . . . . . . 4-2 Final inversion results using a homogeneous initial model for FWI (irst panel), long-wavelength model obtained with DWI (second panel), and FWI result starting from model in the second panel (third panel). Velocity proiles at 0.85 km (fourth panel) for the exact model (black line), initial model (red line), FWI (blue line) and DWI + FWI (green line). . . . . . . . . . .4-20 Test 4. Multiscale FWI results at each steps. From top to bottom, the source bandwidths used are respectively from (a) to (d) those shown in Figure 4-13. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-21 Test 4. 1D vertical proiles of exact (black line), initial (red line) and FWI model (blue line) (left), or multiscale FWI model (blue line) (right). From top to bottom: = 0.4, 1, 1.6 km. . . . . . . . . . . . . . . . . . . . . . . 4-22 Test 4. Long-wavelength model obtained for DWI (top), and FWI result starting from model on the top (bottom). . . . . . . . . . . . . . . . . . . . 4-23 Test 4. 1D vertical proiles of exact (black line), initial (red line) and multiscale FWI model (blue line) (left), or decomposition-based waveform inversion + FWI model (blue line) (right). From top to bottom: = 0.4, 1, 1.6 km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-24 Test 4. Top panels: Shot gathers for observed data (left), initial data (right). Middle panels: Final data obtained with FWI (left), multi-scale FWI (middle) and DWI (right). Bottom panels: Observed (black), initial (red), FWI modeled (blue), multi-scale FWI modeled (green) and DWI modeled (magenta) seismic traces taken at three ofset positions: 0.5 km (top), 1.5 km (middle) and 2.5 km (bottom). . . . . . . . . . . . . . . . . . The diference between the exact velocity model (Figure4-3) and the gradient-constant initial velocity model (Figure 4-4 on the top). (Middle) The smoothed version of the image on the top. (Bottom) The model perturbation of the irst global iteration for the proposed decomposition-based inversion using zero-ofset data and 20 iterations of FWI in the irst step of inversion. Source wavelet without low frequencies less than 5 Hz is used for the two steps of inversion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-26 Model perturbation of the irst iteration of long-wavelength update for the proposed decomposition-based waveform inversion when using large ofset data (1.2 km) in the short-wavelength update step of the inversion. . . . . 4-27 1D vertical proiles of exact model perturbation (black line), model perturbation for DWI applying all the key parameters (red line) and model perturbation for DWI when large ofset data (1.2 km) is used in the shortwavelength update step of the inversion (blue line). From top to bottom: = 0.4, 1, 1.6 km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-28 Model perturbation of the irst iteration of the long-wavelength update for DWI when using 1 iteration instead of 15 iterations in the short-wavelength update step of the inversion. . . . . . . . . . . . . . . . . . . . . . . . . . . 4-29 1D vertical proiles of exact model perturbation (black line), model perturbation for DWI (red line) and model perturbation for DWI when 1 iteration instead of 15 iterations is used in the short-wavelength update step (blue line). From top to bottom: = 0.4, 1, 1.6 km. . . . . . . . . . . . . . . . . 4-30 Model perturbation of the irst global iteration for DWI when using source with low frequency in the short-wavelength update step of the inversion. . 4-31 1D vertical proiles of exact model perturbation (black line), model perturbation for DWI (red line) and model perturbation for DWI when using source with low frequency in the short-wavelength update step of the inversion (blue line). From top to bottom: = 0.4, 1, 1.6 km. . . . . . . . . . . 4-32 Model perturbation of the irst global iteration for DWI when using Born approximation in the short-wavelength update step of the inversion. . . . . 4-33 1D vertical proiles of exact model perturbation (black line), model perturbation for DWI (red line) and model perturbation for DWI when using Born approximation in the short-wavelength update step of the inversion (blue line). From top to bottom: = 0.4, 1, 1.6 km. . . . . . . . . . . . . . 4-34 (Top) Velocity proiles at 1.5 km for the model obtained after the irst step of inversion by using zero-ofset FWI (red line) and zero-ofset iterative migration based on the Born approximation (blue line). (Bottom) Zoom of the zone in green box in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . 4-35 Model perturbation of the irst global iteration for DWI when using source with low frequency in the short-wavelength update step of the inversion. . 4-36 1D vertical proiles of exact model perturbation (black line), model perturbation for DWI (red line) and model perturbation for DWI when using source with low frequency in the short-wavelength update step of the inversion (blue line). From top to bottom: = 0.4, 1, 1.6 km. . . . . . . . . . . 5-1 Acquisition geometry of the Brunei dataset. . . . . . . . . . . . . . . . . . 5-2 One original shot gather (left) and its spectrum (right). . . . . . . . . . . . 5-3 Picking results using the Coppen's algorithm (blue) and smoothed versions (red). From top to bottom, the size of windows are respectively 0.09 s, 0.15 s and 0.21 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-4 Picking results using the STA/LTA algorithm (blue) and smoothed versions (red). From top to bottom, the size of short-term windows are respectively 0.09 s, 0.15 s and 0.21 s, and the size of the long-term windows are respectively 0.9 s, 1.5 s and 2.1 s . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-5 Picking results using the diference-based algorithm (blue) and smoothed versions (red). The size of the window is 0.09 s. . . . . . . . . . . . . . . . 5-6 The shot gather after removal of direct waves (left) and the shot gather after low-pass iltering with cut frequency of 10 Hz (right). . . . . . . . . . 5-7 The shot gather after denoising in the f-k domain (left) and its spectrum (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-8 FWI model provided by Total. . . . . . . . . . . . . . . . . . . . . . . . . . 5-9 Estimated source wavelet in the time domain (left) and in the frequency domain (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-10 Constant density case. The observed data (left) and the calculated data in the model of Figure 5-8 (right). . . . . . . . . . . . . . . . . . . . . . . . . 5-11 Constant density case. Seismic traces of observed (black) and calculated (red) data taken at four ofset positions, from top to bottom: 0.15 km, 2 km, 4 km and 6 km respectively. . . . . . . . . . . . . . . . . . . . . . . . . 5-12 Density model obtained with the Gardner's formula in equation 5.1. . . . . 5-13 Variable density case. The observed data (left) and the calculated data (right) using the density model in Figure 5-12. . . . . . . . . . . . . . . . . 5-14 Variable density case. Seismic traces of observed (black) and calculated (red) data taken at four ofset positions, from top to bottom: 0.15 km, 2 km, 4 km and 6 km respectively. . . . . . . . . . . . . . . . . . . . . . . . . 5-15 The migrated image using data within ofset range of [150-2500] m. Note that the image is not focused in the red circled zone. . . . . . . . . . . . . 5-16 From top to bottom: the migrated image using data with ofset range of [150-275] m, [275-400] m, [400-525] m, [525-650] m and [650-775] m respectively. Note that from top to bottom, the migrated images shift horizontally. The red vertical lines are ixed at = 5.2 km in these images. . . . . . . . 5-17 The migrated image using data with ofset range of [150-775] m. . . . . . . 5-18 Vertical slices taken in Figure 5-16. From left to right, = 3 km, 6 km and 11 km respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-19 The initial velocity model for inversion. . . . . . . . . . . . . . . . . . . . . 5-20 FWI result after 100 iterations. . . . . . . . . . . . . . . . . . . . . . . . . 5-21 The model perturbation of the irst global iteration with DWI. . . . . . . .

Figure 2 - 1 :

 21 Figure 2-1: An example of seismic data. We can observe the direct wave and the relected waves.

Figure 2 - 2 :

 22 Figure 2-2: The original Ricker wavelet (black) and the estimated source (red) in the frequency domain (left) and in the time domain (right).

Figure 2 - 3 :

 23 Figure 2-3: Wavenumber illumination. One source-receiver pair, one scattering point and one frequency in the data provides one wavenumber in the image space [Huang and Schuster, 2014].

Figure 2 - 4 :

 24 Figure 2-4: Approximate resolution limits for (a) diving waves, (b) and (c) relectionrelated transmission, (d) difraction-related transmission and (e) relection migration [Huang and Schuster, 2014].

Figure 2 - 5 :

 25 Figure 2-5: First Fresnel zone for the specular relection [Huang and Schuster, 2014].

Figure 2 -

 2 Figure 2-7: a) Ray diagram for interbed multiples generated by a difractor in a thin layer and b) the associated mirror sources diagram [Huang and Schuster, 2014].

Figure 2

 2 Figure 2-8: 1D exact velocity model.

Figure 2

 2 Figure 2-9: Four Ricker wavelets. From top to bottom and from left to right, the central frequency is respectively 4 Hz, 8 Hz, 12 Hz, and 16 Hz.

Figure 2 -Figure 2 -

 22 Figure 2-10: (Left) Gradient of the irst iteration of FWI with source wavelets displayed in Figure 2-9. The central frequency is respectively 4 Hz, 8 Hz, 12 Hz, and 16 Hz. (Right) The bandwidth of the vertical slice taken at the middle of each model on the left.

Figure 2 -

 2 Figure 2-12: Comparison of the convergence rate of local methods applied on a 2D linear inverse problem. The problem is solved using a) Gauss-Newton, b) gradient descent, c) preconditioned gradient and d) conjugate-gradient methods. The gradient in c) was

Figure 2 -

 2 14 to Figure 2-19 show an example that compares the result of conventional FWI and multiscale FWI through a 2D synthetic model. The true model (Figure 2-14) is a constant gradient model with a circular anomaly. The initial model (Figure 2-15) is just the constant gradient model. The circular velocity anomaly is missing. The source wavelet used for FWI is a Ricker wavelet with central frequency of 8 Hz. The result of FWI after 100 iterations is displayed in Figure 2-16, and the model perturbation with respect to the initial model is displayed in Figure 2-17. Compared to the exact model, only the edges of the circular anomaly are retrieved. FWI fail to recover the long-wavelength part inside the anomaly. For the multiscale FWI, we use 7 singular frequencies: 1.8 Hz, 2.5 Hz, 3.5 Hz, 5 Hz, 7 Hz, 10 Hz, 14 Hz.

Figure 2 -

 2 Figure 2-13: Principle of cycle-skipping problem[Virieux and Operto, 2009]. If the phase mismatch is less than one half of the period, the local optimization methods could adjust the correct phase. In the contrary case, the local optimization methods adjust these two signals with one period shift and the inversion falls into a local minimum.[Virieux and Operto, 2009]

Figure 2 -

 2 Figure 2-14: Exact velocity model.

Figure 2 Figure 2 -

 22 Figure 2-16: Conventional FWI result.

Figure 2 -FWIFigure 2 -

 22 Figure 2-18: Multiscale FWI results. From top to bottom, the frequency used is 1.8 Hz, 2.5 Hz, 3.5 Hz, 5 Hz, 7 Hz, 10 Hz, 14 Hz respectively. For each frequency, 25 iterations of FWI is performed.

Contents 3 .

 3 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3 Inversion strategy . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.1 Choice of ofset range . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3.2 Iterative inversion . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3.3 Gradient smoothing . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.4 1D layer model case . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Comparison with similar methods . . . . . . . . . . . . . . . . 57 3.5 Waveield decomposition methods . . . . . . . . . . . . . . . . 62 3.5.1 Methods based on one-way wave equation . . . . . . . . . . . . . 63 3.5.2 Methods based on two-way wave equation . . . . . . . . . . . . . 63 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Résumé du chapitre 3

  d'ondes et une partie de grandes longueurs d'ondes après la décomposition des champs d'ondes. Nous appelons cette nouvelle méthode en anglais Decomposition-based Waveform Inversion (DWI). Notons que la formule du gradient de FWI somme le terme des courtes longueurs d'ondes et le terme des grandes longueurs d'ondes avec le même poids. Pour les données de rélexion , le terme des grandes longueurs d'ondes est généralement beaucoup plus faible que le terme des courtes longueurs d'ondes. C'est la raison que la FWI classique souvent n'arrive pas à retrouver les grandes longueurs d'ondes du modèle. Au lieu de mettre des poids diférents, nous proposons d'inverser séparément les deux parties. Supposons que l'on commence par un modèle lisse, d'abord nous inversons les courtes longueurs d'onde en utilisant les données zéro ofset. Dans cette étape, nous efectuons l'inversion itérative pour retrouver une rélectivité bien focalisée, qui est essentielle pour l'étape suivante. Notons que cette étape se comporte comme une minimisation au sens des moindres carrés (en anglais, Least-Squares Migration, LSM). L'utilisation des données courts ofsets est importante pour ne pas introduire les grandes longueurs d'ondes dans cette étape. Après l'inversion itérative, le champ d'ondes direct et le champ d'onde rétropropagé maintenant contient les champs descendants et montants. Dans la deuxième étape, nous efectuons la mise à jour des grandes longueurs d'ondes. Dans cette étape, les données avec plus d'ofsets sont utilisées pour obtenir les informations moveout. Un lissage sur le gradient est appliqué ain d'éliminer les résidus de courtes longueurs d'ondes, causée par les limitations d'acquisition. Une fois que nous avons mis à jour les grandes longueurs d'ondes du modèle, nous recommençons la mise à jour des courtes longueurs d'ondes du modèle puisque les positions des rélecteurs ne sont pas compatibles avec le modèle de vitesse actuel et devraient être mis à jour en conséquence. Puis nous continuons à alterner entre ces deux étapes. Le point clé de cette méthode d'inversion est ce schéma alternatif. La limitation d'ofset et l'application de l'inversion itérative dans la première étape et le lissage du gradient dans la deuxième étape sont les points importants pour assurer le découplage entre les grandes longueurs d'onde et les courtes longueurs d'onde.

Figure 3 - 1 :

 31 Figure 3-1: Raypaths and illumination in the model space for diferent components of the gradient. Schematic wave propagation in 1D layer model. + : direct forward waveield (red solid line). -: relected forward waveield (red dashed line). -: direct backward waveield (blue solid line) and + : relected backward waveield (blue dashed line).

Figure 3 - 2 :

 32 Figure 3-2: Comparison of conventional FWI (left) and the decomposition-based waveform inversion (right).

Figure 3 - 3 :

 33 Figure 3-3: The initial velocity model. It is the sum of the homogeneous model at 2 km/s and the true relectivity model. The true relectivity is obtained by subtracting the smoothed true velocity model from the exact velocity model of Figure 2-8.

Figure 3 - 4 :Figure 3 - 5 :Figure 3 Figure 3

 343533 Figure 3-4: Correlation of decomposed forward and backward waveields for one sourcereceiver pair with the model in Figure 3-3. The source is located at 0.85 km. From top to bottom, the receivers are at 0.085 km, 0.85 km and 1.615 km respectively. Correlation of + and -(left); -and + (right).

Figure 3 - 8 :Figure 3 - 9 :

 3839 Figure3-8: Correlation of decomposed forward waveield and backward waveield using all the available sources at the surface, starting from the model in. From top to bottom, correlation of + and -; -and + ; + and + ; -and -respectively. (right) Spectra of the vertical velocity proiles taken at 0.85 km for igures on the left respectively.

Figure 3 -

 3 Figure 3-13: One snapshot of the waveield propagated in the model in Figure 2-8 at = 0.525 s.

Figure 3 -Figure 3 -Figure 3 -

 333 Figure 3-16: Filtered parts of the initial image in Figure 3-14.

Figure 3 -

 3 Figure 3-19: (Left) Partial derivative of the waveield with respect to . (Right) Partial derivative of the waveield with the respect to .

Figure 3 -

 3 Figure 3-20: The product of the two images in Figure 3-19.

Figure 3 -Figure 3 -

 33 Figure 3-21: The decomposed waveields using Poynting-vector method. + (left) and - (right).

Figure 3 -

 3 Figure 3-23: Decomposed waveields using the mask in Figure 3-22 on the right. + (left) and -(right).

Figure 3 -

 3 Figure 3-24: A scale example of curvelet [Chauris et al., 2006].

Figure 3 -Figure 3 -

 33  The direction (left) and scale (right) ilter in the curvelet domain[START_REF] Chauris | Seismic demigration/migration in the curvelet domain[END_REF].Surface (km)

Figure 3 - 28 :

 328 Figure 3-28: One direction after the curvelet transform for the two snapshots in Figure 3-27.

Figure 3 - 29 :

 329 Figure 3-29: Decomposed waveields using the curvelet method. + (left) and -(right).

Figure 3 - 30 :

 330 Figure 3-30: One snapshot of the waveield propagated in the Marmousi model.

Figure 3 - 31 :

 331 Figure 3-31: Decomposed waveields using the 2D FT method (irst panel), Poynting method (second panel) and curvelet method (third panel). + (left) and -(right).

  Figure 4-1: 1D layer exact velocity model.

Figure 4 - 3 :Figure 4 - 4 :Figure 4 - 5 :

 434445 Figure 4-3: 2D exact velocity model [Perrone, 2013].

  Figure 4-6: Test 1. FWI result after 100 iterations.

  with low frequency and poor initial modelIn the second test, we apply FWI starting from the homogeneous model shown in Figure4-4 on the bottom. The source wavelet remains the same as in the irst test. In this test, as the initial model is very far from the true model, the conventional FWI fails to converge to the exact model. The inal FWI result is shown in Figure4-9 and Figure4-10, and is considered as a reference for the fourth test with a poor initial model and without low frequencies. The comparison of the exact data, the initial data and inal data is displayed in Figure4-11. We note that for this initial model FWI falls into a local minimum. 4.3.3 Test 3: inversion without low frequency and good initial model In the third test, the initial model is the same as in the irst test (Figure 4-4 on the top). But low frequencies are iltered in the original Ricker wavelet to be more realistic (Figure 4-5).

Figure 4 - 7 :Figure 4 Figure 4 Figure 4 -

 47444 Figure 4-7: Test 1. 1D vertical proiles of exact (black), initial (red) and FWI model (blue). From top to bottom: = 0.4, 1, 1.6 km.

Figure 4 -

 4 Figure 4-11: Test 2. Top panels: shot gathers for observed data (left), initial data (middle) and inal data obtained with FWI (right). Bottom panels: observed (black), initial (red) and modeled (blue) seismic traces taken at three ofset positions: 0.5 km (top), 1.5 km (middle) and 2.5 km (bottom).

Figure 4 -

 4 Figure 4-12: Test 3. FWI result after 100 iterations.

Figure 4 -Figure 4 -Figure 4 -Figure 4 -

 4444 Figure 4-13: Test 3. Source bandwidths used for the four steps in multiscale FWI. The black lines represent the initial Ricker bandwidth with central frequency of 8 Hz and the red lines represent the source bandwidth after the Butterworth ilter. The cut frequencies for the band-pass ilter are respectively (a) 8-9 Hz, (b) 8-12 Hz and (c) 8-15 Hz. A high-pass ilter with cut frequency of 8 Hz is used in (d).

  Figure 4-19: Test 4. FWI result

Figure 4 - 4 Figure 4 -Figure 4 -

 4444 Figure 4-21: Test 4. 1D vertical proiles of exact (black line), initial (red line) and FWI model (blue line) (left), or multiscale FWI model (blue line) (right). From top to bottom: = 0.4, 1, 1.6 km.

Figure 4 -Figure 4 - 28 :

 4428 Figure 4-27: 1D vertical proiles of exact model perturbation (black line), model perturbation for DWI applying all the key parameters (red line) and model perturbation for DWI when large ofset data (1.2 km) is used in the short-wavelength update step of the inversion (blue line). From top to bottom: = 0.4, 1, 1.6 km.

Figure 4 - 29 :Figure 4 - 30 :Figure 4 - 31 :Figure 4 -Figure 4 - 33 :Figure 4 -

 42943043144334 Figure 4-29: 1D vertical proiles of exact model perturbation (black line), model perturbation for DWI (red line) and model perturbation for DWI when 1 iteration instead of 15 iterations is used in the short-wavelength update step (blue line). From top to bottom: = 0.4, 1, 1.6 km.

Figure 4 - 35 :Figure 4 -

 4354 Figure 4-35: Model perturbation of the irst global iteration for DWI when using source with low frequency in the short-wavelength update step of the inversion.

Figure 5 - 1 :

 51 Figure 5-1: Acquisition geometry of the Brunei dataset.

Figure 5 - 2 :

 52 Figure 5-2: One original shot gather (left) and its spectrum (right).

Figure 5 - 3 :Figure 5 - 4 :Figure 5 - 5 :Figure 5

 5354555 Figure 5-3: Picking results using the Coppen's algorithm (blue) and smoothed versions (red). From top to bottom, the size of windows are respectively 0.09 s, 0.15 s and 0.21 s.

Figure 5 - 7 :

 57 Figure 5-7: The shot gather after denoising in the f-k domain (left) and its spectrum (right).

Figure 5 Figure 5

 55 Figure 5-8: FWI model provided by Total.

Figure 5 -

 5 Figure 5-10: Constant density case. The observed data (left) and the calculated data in the model of Figure 5-8 (right).

Figure 5 -

 5 Figure 5-11: Constant density case. Seismic traces of observed (black) and calculated (red) data taken at four ofset positions, from top to bottom: 0.15 km, 2 km, 4 km and 6 km respectively.

Figure 5 -

 5 Figure 5-14: Variable density case. Seismic traces of observed (black) and calculated (red) data taken at four ofset positions, from top to bottom: 0.15 km, 2 km, 4 km and 6 km respectively.

Figure 5 -

 5 Figure 5-15: The migrated image using data within ofset range of [150-2500] m. Note that the image is not focused in the red circled zone.

Figure 5 -

 5 Figure 5-17: The migrated image using data with ofset range of [150-775] m.

Figure 5 -Figure 5 -

 55 Figure 5-18: Vertical slices taken in Figure 5-16. From left to right, = 3 km, 6 km and 11 km respectively.

  de données réelles, et présentons les pré-traitements efectués sur les données et les diicultés que nous rencontrons avec l'hypothèse acoustique de la propagation des ondes.

	1.1	Seismic imaging
	1.1.1	General principles

4, nous appliquons la FWI, la FWI multi-échelle et la nouvelle méthode d'inversion proposée sur deux modèles synthétiques et comparons les résultats. Nous analysons aussi l'inluence des paramètres clés pendant l'inversion sur le résultat d'inversion. Dans le chapitre 5, nous appliquons la FWI et la nouvelle méthode d'inversion proposée sur un jeu

Table 4 .

 4

	N ∘	Frequency content	Initial model
	Test1	with low frequency	gradient-constant model
	Test2	with low frequency	homogeneous model
	Test3 without low frequency gradient-constant model
	Test4 without low frequency	homogeneous model

1: List of tests for the 2D model.

  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The DWI method . . . . . . . . . . . . . . . . . . . . . . . . . . . Physics of the earth . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Role of density . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.3 Inversion techniques . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents
	6.1
		6.1.1	Basic FWI theory and resolution analysis	. . . . . . . . . . . . .
		6.1.2
		6.1.3	Applications to synthetic models and real data	. . . . . . . . . .
	6.2	Perspectives . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1

Mots Clés: Inversion des formes d'ondes, décomposition des champs d'ondes, imagerie sismique.

• Multi-parameter inversion. The principle of FWI is to use all types of waves (direct waves, relections, refractions, and multiples) to resolve diferent model parameters

Some authors[Lee and Kim, 2003, Zhou andGreenhalgh, 2003] propose to normalize the waveield to make the inversion independent of the source function, such that the source estimation step is not needed for these methods.

Acknowledgements

List of tables

We apply this method on a 2D waveield (, , ). The waveield is computed in the model of Figure 2-8 with the source at = 0.85 km at the surface. A snapshot at = 0.525 s is shown in Figure 3-13. The algorithm for the vertical decomposition is described as follows.

For the central trace ( = 0.85 km) of the image:

1. The image indomain (, ) is shown in Figure 3-14 (left).

2. The 2D FT of (, ) gives  (, z ) in thez domain, as shown in Figure 3-14 (right).

3. We select appropriate quadrants using the ilters in Figure 3-15.  (, z ) is separated into a down-going part  z+ (, z ) and an up-going part  z-(, z )

(Figure 3-16).

4. Then we apply the 2D inverse FT on  z+ (, z ) and  z-(, z ) to get the decomposed images in thedomain (Figure 3-17).

5. We repeat steps (1)-( 4) for each horizontal coordinate , and we sum all the contributions to get + (, , ) and -(, , ), shown in Figure 3-18.

The drawbacks of the FT-based decomposition are: 2003], choosing in equation ( 2.11) such that the product of and the maximum value of the gradient remains a constant ixed at the irst iteration. We smooth the gradient in this step to remove the residual short-wavelength components that are caused by acquisition limitations and decomposition artefacts. The smoothing is performed in the frequency domain as shown in Chapter 3. Once the background velocity is updated, a new global iteration is started.

In the long-wavelength update step, as we perform only one iteration of inversion, to compensate the geometrical spreading, we need to weight the gradient by the inverse of the Hessian to speed up convergence. The inverse of the Hessian helps to deconvolve the gradient, and to compensate for the geometrical spreading and for illumination. But the computation of the inverse of the Hessian is extremely expensive [START_REF] Clément | Migration-based traveltime waveform inversion of 2-D simple structures: A synthetic example[END_REF]. As the test model is relative simple, we use a linear preconditioning to compensate the geometrical spreading assuming that the velocity increases linearly with depth. We also tested the preconditioning presented in [START_REF] Plessix | Frequency-domain inite-diference amplitude -preserving migration[END_REF] which calculates an approximation of the inverse of the Hessian by the inverse of the norm of the waveield. However, the linear preconditioning is suicient in this case and provides a similar result.

The long-wavelength velocity model we obtain after 30 global iterations is shown in Figure 4-16 on the top. Starting from this model, we run 100 iterations of the conventional FWI (inverting for all frequencies at once) and the inal result is shown in Figure 4-16 on the bottom. The whole procedure could be regarded as a two-step FWI worklow, consisting of the DWI followed by the conventional FWI. Compared to the multi-scale FWI (Figure 4-17 and 4-18), we can observe that the relectors are better positioned and the central velocity anomaly, which is absent in the initial model, is successfully recovered. As presented in Chapter 3, the construction of the short-wavelength relectivity model is a critical step of the proposed algorithm, since a more correct relectivity allows for a better long-wavelength update. The proposed inversion scheme combines several elements in the short-wavelength update step to build a true-amplitude, focused and sharp relectivity model: (1) zero-ofset data; (2) nonlinear iterative inversion (20 iterations in our example) and ( 3) a source with high-frequency content. In this section, we show the importance of these elements by analyzing and comparing the model perturbations that we obtain after the irst iteration of the long-wavelength update. The analysis of the model perturbation is similar to the analysis of the sensitivity kernel in the gradient.

In the following tests, the initial velocity model is the gradient-constant velocity model (Figure 4-4 on the top). The true model perturbation (the diference between the true velocity model and the initial velocity model) is displayed in Figure 4-25 on the top. The model perturbation retrieved after the irst global iteration (short-wavelength update followed by long-wavelength update) by using the three above-mentioned elements is shown in Figure 4-25 on the bottom. This is the reference for comparison. The result is close to the smoothed version of the true model perturbation shown in Figure 4-25 in the middle as if the new method reduces the nonlinearity of the problem and speeds up the convergence.

In the following, we analyze the efect of partly applying all of the above-mentioned elements and compare the velocity perturbation to the one displayed in Figure 4-25 on the bottom.

Figure 4-26 shows the model perturbation after the irst iteration of long-wavelength update when full-ofset data (1.2 km ofset) instead of zero-ofset data is used in the short- Offset (km)

Amplitude

Figure 5-24: Amplitude variation with ofset for constant density case. Blue: observed data; red: calculated data; green: the ratio of calculated data and observed data; magenta:

1/ √ as a reference.

Inversion techniques

For DWI, iterative inversion is performed at the short-wavelength update step to it the data. This step is quite computational expensive. However this iterative step is important, since it facilitates the long-wavelength update of the model. For further tests, we would like to implement this step in the pseudo-time domain to avoid the re-computation of the short-wavelength model at each step [Plessix, 2012, Brossier et al., 2015]. It should be noted that for DWI, in the long-wavelength update step, the step length of the gradient descent method is deined at the beginning. Thus the decrease of the objective function is not assured during all the steps of the inversion. If we use the quadratic interpolation to determine the step length at each iteration, we need to carry out two extra direct problems, however here the direct problem corresponds to the iterative inversion, thus this is quite computational expensive. If this iterative inversion step could be performed in the pseudo-time domain, the computational cost could be largely reduced. Besides, in this thesis, DWI is mainly used for relections. However, in theory, DWI is not limited to work with relections only, which is diferent from other relection-based waveform inversion methods, and could incorporate other types of waves at the same time, such as diving waves and difractions. For future research, we would like to combine these diferent types of waves and to see if the combination could facilitate the inversion.

[ [START_REF] Biondi | Tomographic full waveform inversion (tfwi) by combining full waveform inversion with wave-equation migration velocity anaylisis[END_REF] introduce the subsurface ofset to the relection-based waveform inversion to mitigate the nonlinearity of the problem. We would like to investigate more this point, by introducing the subsurface ofset into DWI.

In order to mitigate the problem of data amplitude, we would like to test other objective functions, that are less sensitive to amplitude, such as the correlation-based objective function [van Leeuwen and Mulder, 2008] and normalized objective function. We need to investigate on the feasibility of combining the new objective functions with waveield decomposition.

Finally, with the development of high performance computers, the application of the curvelet-based decomposition method on real datasets could become soon feasible. To further reduce the computational cost, we would like to test the idea of correlating the In this thesis, I propose a variant of FWI based on the wavefield separation, typically between up-and down-going waves, to mitigate the nonlinearity of the problem. The new method consists of decomposing the gradient of FWI into a short-wavelength part and a longwavelength part after wavefield decomposition. The inversion is performed in an alternating fashion between these two parts. We apply this method to several case studies and show that the new method is more robust especially for constructing the long-wavelength model.