Yves-Marie Schilke

Tiago Batany

Jihane Barros

Julien Belhadj

Keurfon Cotton

Alexandre Luu

Tianyou Kazantsev

Michelle Zhou

Ahmed Almakari

Maxime Jabrane

Zufri Go- Dano

Ibrahim Malik

Mohamadhasan Joshua Lartey

Mohomadian Sarvandani

Rita Abou

Milad Farshad

Pranav Au- Dhkhasi

Asma Berrhouma

Fengjuan Chen

Kaiwen Chang

Naixin Chen

Qinglin Deng

Cecilia Dip

Sebatien Drouyer

Minyu Fan

Mario Fernandez

Yuxin Fu

Fang Gao

Xin Gao

Weiguang He

Yanhui He

Hanqin Huang

Shun Huang

Song Jin

Bo Kang

Eleni Kolokyta

Dingan Liao

Zimin Li

Fang Lu

Dac-Thuong Ngo

Dandan Niu

Oleg Ovcharenko

Meng Pang

Zilin Peng

Bin Qiao

Yanfang Qin

Sara Rachdi

Yi Rong

Yang Shen

Xiangdong Song

Yang Song

Adam Sukma Putra

Yuwei Sun

Yuxin Sun

Chaojing Wang

Dandan Wang

Haisheng Wang

Jialan Wang

Shuaitao Wang

Tianyuan Wang

Wei Wang

Zhuangzhuang Wang

Xiaojing Wu

Mingguan Yang

Shuran Yang

Tiantian Yin

Xi Yu

Shihao Yuan

Fan Zhang

Yizhuo Zhang

Zezhong Zhang

I express my thanks to Yubing for showing me around when I first arrived at Fontainebleau and many research discussions on the trains between Paris and Fontainebleau. I would like to thank Chen Tang for inspiring me to work on seismic imaging when I first began to learn geophysics. I would like to thank my best friend Dezhi Wang for keeping encouraging me to pass through the hard time. I can not write down all the helps from my friends in every details, but I will keep them always in my heart.

Finally, I would like to express my deepest gratitude to my family. Without their supports and encouragements, I could never

I would like to thank all the people who have contributed to this thesis work. The process of research during the PhD is a journey of exciting the curiosity, pursuing the knowledges and answering the unknowns. I always remember all the people's helps during the wonderful journey.

I would like to thank my supervisor Prof. Hervé Chauris with my sincere gratitude. I began my research with him from my internship until the end of my PhD study. He always supported and encouraged me to conduct the research with great patiences. He always had an insight on the problems and gave me good suggestions. It is the most impressive that he always can describe the geophysical problems in a very precise mathematical language. During the time of writing my thesis, he spent much time on correcting my manuscript lines by lines. He guided me to present my works in an easy understanding way with good timing and rhythm for my presentation at EAGE and my final thesis defense. His passion and rigorous attitude on research inspires me a lot, his optimistic philosophy, healthy life style and humor influence me a lot. I appreciate him for all those helps and positive impacts from him.

I would like to thank the juries of my thesis committee. I thank Alessandra Ribodetti and Romain Brossier for reviewing my thesis in details with professional comments and recommendations. I thank Jean-Pierre Vilotte for chairing the thesis defense. I thank Mark Noble and Fatiha Gamar for examining my work.

I would like to thank my school, MINES ParisTech for funding me on the research, and HPC of our Geoscience Center for implementing the computations.

I would like to thank Baotao Qin, Nobuaki Fuji, Patrick Goblet and Laurent De Windt for examining my work during the annual qualification meeting with fruitful advise. I also express my thanks to Feng Yang for his encouragements and helps for preparation of the annual qualification meeting.

Many thanks to the members of the geophysics team at MINES ParisTech. I especially thank Véronique Lachasse for her helps to facilitate the administrative affairs. I thank the permanent researchers Daniela Donno, Alexandrine Gesret and Mark Noble who are also my master teachers, and Pierre Dublanchet, and the scientific visitors Nidhal Belayouni and Pierre-Francois Roux, for their kindness. I also would like to thank all my colleagues Fang Wang, Charles-Antoine Lameloise, Emmanuel Cocher, Sven 1.9 Workflow of the sequential (left) simultaneous (middle) and hybrid (right) strategies. τ approximates to the inverse of Q. v k = f (v p , τ, ω c ) is the kinematic relationship, ω c is the central frequency of the source (more details in Section 2.6.3 and Chapter 4). . . . . . . . . . . . . . . . . . . 16 1.10 The sketch of searching direction of strategies: the kinematic curves (dashed green line) superimposed to the contour lines of the misfit distribution (gray solid) with the true values being v p = 2000 m/s and τ = 0.023 (Q = 50). The blue dots and red hollow circles are the searching direction of the simultaneous and hybrid strategies, respectively. For the hybrid strategy: v k updates (1 to 3) and the v p updates for fixed v k (3 to 5), more details in Chapter 4. . . . . . . . . . . . . . . . . . . . . . . 17

2.1 (Left) Generalized Maxwell Body: MB denotes Maxwell body, formed by a spring with elastic modulus M l and a dashpot with viscosity coefficient η l . n pairs of Maxwell bodies and one spring with elastic modulus M H in parallel make up the GMB, here denoted as GMB-EK to refer to the model from paper of [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]. (Right) Generalized Zener Body: ZB denotes the Zener body. n pairs of Zener bodies form the GZB. There are two kinds of Zener bodies. M Rl , δM l , M 1l and M 2l are elastic modulus, and η l is the viscosity coefficient (adapted from [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2 The stress relaxation function of a SLS. The strain is unchanged, the imposed stress gradually decreases from unrelaxed M U to relaxed M R (adapted from [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF]. . . . . . . . . . . . . . . . . . . . . . . 26 2.3 The optimization frequency range is from 1.5 Hz to 100 Hz for desired Q = 15, with empirically given Q pre with 15 (top), 13 (middle) and 12 (bottom). The τ method is the analytical solution with different predefined Q pre , the modified method is the global searching method. . . . 31 2. 4 The optimization frequency range is from 1.5 Hz to 100 Hz for Q 0 = 50 with empirically given Q pre with 48. . . . . . . . . . . . . . . . . . . . 31 2.5 Same as Figure 2.4, but for Q = 100 with empirically given Q with 98. . 32 2.6 Discretization of the media on a staggered grid, p is the pressure wavefields, v x and v z are particle velocity components. . . . . . . . . . . . . 37 2.7 One snapshot of wavefields: stable wave propagation without numerical dispersion (a), unstable wave propagation (b), stable condition with numerical dispersion (c) (adapted from [START_REF] Thorbecke | 2D Finite-Difference wavefield modelling[END_REF]. . . . . . . . . 39 2.8 Snapshot without CPML at 1.4 s (a) and 1.7 s (b). . . . . . . . . . . . . 40 2.9 Same as Figure 2.8, but snapshot with CPML at 1.4 s (a) and 1.7 s (b). . Dans ce chapitre, j'introduis d'abord l'atténuation sismique et discute des effets de dissipation et de dispersion intrinsèque et liée au "scattering". Je donne aussi des valeurs pour des roches communes. J'introduis plusieurs modèles de Q. Dans la deuxième partie de ce chapitre, je donne un état de l'art sur l'imagerie des structures avec de l'atténuation. Enfin, je conclus avec les motivations pour cette recherche et avec l'organisation de ce travail. Les détails techniques sont introduits dans le chapitre 2.

Seismic attenuation is a useful physical parameter to describe the properties of specific geological bodies, e.g., saturated rocks and gas clouds. The introduction of attenuation can improve the seismic migration image and remove the blurring of images under the gas clouds. Seismic attenuation is also more sensitive than velocity changes for fluid detection. Thus, attenuation combined with velocity analysis provide extra information on the subsurface.

In this chapter, I first introduce seismic attenuation including dissipation and dispersion phenomenon, intrinsic and scattering attenuation. I also give attenuation values for common rocks. Based on rock experiments on attenuation, I introduce different Q models. The second part of this chapter is an overview of the methods on attenuation imaging. Finally, I conclude with the motivation for my research and with the organization of the thesis. In this chapter, I only introduce a few notations and equations. More technical details are presented in Chapter 2.

Seismic attenuation 1.Attenuation phenomenon and impacts

Seismic attenuation is the term to characterize the seismic energy loss through the subsurface media, including the intrinsic and scattering attenuation (Section 1.1.2). It is a useful property to characterize geological bodies as lithology discriminator. Due to attenuation, absorption impacts waves, resulting in change in the shape of waveform after passing through attenuation media. The energy absorption can be expressed by the quality factor Q. The reciprocal of Q is defined as the proportion of energy loss per wave propagation cycle [START_REF] Aki | Quantitative seismology[END_REF][START_REF] Sheriff | Exploration seismology[END_REF][START_REF] Wang | Seismic Inverse Q Filtering[END_REF].

1 Q = ∆E 2πE . (1.1)
On the one hand, seismic attenuation information can supply extra information on the subsurface media solely or combined with velocity information. Attenuation is more sensitive to pore pressure, saturation, and fracture compared to velocity changes [START_REF] Gardner | Effects of pressure and fluid saturation on the attenuation of elastic waves in sands[END_REF][START_REF] Winkler | Seismic attenuation: Effects of pore fluids and frictional?sliding[END_REF][START_REF] Klimentos | Attenuation of P-and S-waves as a method of distinguishing gas and condensate from oil and water[END_REF][START_REF] Muller | Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks -A review[END_REF].

For reservoir detection and characterization, the information of attenuation is crucial. For example, the large attenuation or small Q values on seismic image indicate gas clouds in Figure 1.1. Further more, different ratios of Q p and Q s associated to P and S waves, can determine the reservoir contents (Winkler andNur, 1979, 1982;[START_REF] Cadoret | Fluid distribution effect on sonic attenuation in partially saturated limestones[END_REF]. Lab experiment of sandstone samples are resonated at frequency mainly from 500 Hz to 1700 Hz, to study the attenuation related to fractional sliding and fluid flow mechanism [START_REF] Winkler | Seismic attenuation: Effects of pore fluids and frictional?sliding[END_REF]. Compared to dry rock, the P and S wave attenuation of partial water saturation rock increases significantly, which leads to a larger P wave attenuation than S wave attenuation (i.e., Q -1 p > Q -1 s ). As the water saturation increases, the S wave attenuation reaches the peak value until the full saturation, meanwhile the P wave attenuation increases to maximum value as water saturation is between 60-90 percent. Those phenomenons of lab experiment further extend to the frequency range of exploration seismology [START_REF] Murphy | Effects of partial water saturation on attenuation in Massilon Sandstone and Vycor Porous Glass[END_REF][START_REF] Walls | Seismic attenuation for reservoir characterization[END_REF] and well logging [START_REF] Klimentos | Attenuation of P-and S-waves as a method of distinguishing gas and condensate from oil and water[END_REF][START_REF] Qi | Sonic Qp/Qs ratio as diagnostic tool for shale gas saturation[END_REF]. Partly under the exploration seismology frequency range from 25 Hz to 400 Hz, [START_REF] Murphy | Effects of partial water saturation on attenuation in Massilon Sandstone and Vycor Porous Glass[END_REF] demonstrates the same conclusion that P and S wave attenuation are sensitive to partial water saturation, leading to

Q -1 p > Q -1
s . The sonic well logging (i.e., 10-20 kHz) results indicate that Q -1 p is larger than Q -1 s when the gas or condensate is inside the sandstone [START_REF] Klimentos | Attenuation of P-and S-waves as a method of distinguishing gas and condensate from oil and water[END_REF]. Similar values of P wave attenuation and S wave attenuation indicate the full saturation of water or water + oil. In global seismology, the attenuation model retrieved from Rayleigh wave provides extra information on the temperature, composition, melt content, and volatile content of the rocks of Earth [START_REF] Dalton | Global models of surface wave attenuation[END_REF][START_REF] Bao | Imaging Rayleigh wave attenuation with USArray[END_REF]. Attenuation is more sensitive than velocity to temperature increasing for fixed subsurface pressure [START_REF] Anderson | The anelasticity of the mantle[END_REF][START_REF] Sato | The use of laboratory velocity data for estimating temperature and partial melt fraction in the low velocity zone: comparison with heat flow and electrical conductivity studies[END_REF][START_REF] Karato | Importance of anelasticity in the interpretation of seismic tomography[END_REF][START_REF] Goes | Shallow mantle temperature under Europe from P and S wave tomography[END_REF][START_REF] Jackson | Grain-size sensitive seismic wave attenuation in polycrystalline olivine[END_REF]. Besides, increased attenuation and lower velocity can be caused by partial melting of mantle materials (Hammond and Humphreys, 2000b,a;[START_REF] Faul | Shear wave attenuation and dispersion in melt-bearing olivine polycrystals: 2. Microstructral interpretation and seismological implications[END_REF]. Combined with velocity information, attenuation helps to reduce the ambiguity for the interpretation of media state of Earth interiors.

On the other hand, seismic attenuation information can be used to improve the resolution of seismic images. For the scope of exploration seismology, it supplies the extra information to improve the resolution of seismic migration image. For linearized inversion, Q model can be applied for Q-compensated migration to remove the image blurring under the gas clouds (Figure 1.2) [START_REF] Zhang | Compensating for visco-acoustic effects in reverse-time migration[END_REF][START_REF] Zhou | Compensating attenuation due to shallow gas through Q-tomography and Q-PSDM, a case study from Brazil[END_REF][START_REF] Fletcher | Wavepath-consistent effective Q estimation for Q-compensated reverse time migration[END_REF]Zhu et al., 2014;[START_REF] Gamar | Success of highresolution volumetric Q-tomography in the automatic detection of gas anomalies on offshore Brunei data[END_REF][START_REF] Gamar | Image Quality Enhancement Using Volumetric Q-tomography and Q-PSDM -Martin Linge Case Study[END_REF]. Q model also can be used as a passive parameter for visco-elastic least-squares migration to improve the resolution of migration images [START_REF] Dutta | Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation[END_REF]. For non-linear visco-elastic full waveform inversion, the Q-compensated inversion can accelerate the convergence speed [START_REF] Xue | Accelerating full-waveform inversion with attenuation compensation[END_REF].

Intrinsic and scattering attenuation

Seismic wave attenuation includes intrinsic attenuation [START_REF] Voigt | Ueber innere Reibung fester Korper, insbesondere der Metalle[END_REF][START_REF] Ricker | The form and laws of propagation of seismic wavelets[END_REF][START_REF] Kolsky | The propagation of stress pulses in viscoelastic solid[END_REF][START_REF] Jackson | Physical mechanisms of seismic?wave attenuation[END_REF][START_REF] Aki | Quantitative seismology[END_REF] and scattering attenuation (Aki, 1980a,b;[START_REF] Shapiro | Seismic Attenuation By Scattering: Theory and Numerical Results[END_REF][START_REF] Shapiro | Elastic waves in random media[END_REF][START_REF] Hong | Scattering attenuation ratios of P and S waves in elastic media[END_REF][START_REF] Bongajum | Investigating seismic wave scattering in heterogeneous environments and implications for seismic imaging[END_REF][START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous earth[END_REF], which cause wave dissipation and dispersion. The total quality factor [START_REF] Aki | Scattering and attenuation[END_REF][START_REF] Mukhopadhyay | Variation of intrinsic and scattering attenuation with depth in NW Himalayas[END_REF][START_REF] Fielitz | Intrinsic and scattering attenuation as derived from fluid induced microseismicity at the German Continental Deep Drilling site[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF] is given as

1 Q = 1 Q int + 1 Q sct , (1.2)
where Q is the total quality factor, Q int is the intrinsic quality factor and Q sct is the scattering quality factor. Scattering is a common wave phenomenon and can happen in any kind of media with scatters or highly heterogeneous media. Wave energy is transferred to later arrivals or to waves from other directions caused by scattering attenuation (Figure 1.3). Meanwhile, intrinsic attenuation is the property of anelastic material due to the internal friction. The internal friction transfers wave energy to heat, which results in energy loss and lack of high frequency waves (Figure 1 The typical Q values of rocks are shown in Table 1.1. Q p is the quality factor of P waves and Q s is the quality factor of S waves. Seismic data lack the high frequency parts and have lower wave energy after traveling through the Earth. To explain the phenomenon, people consider the Earth as a frequency filter or a visco-elastic media [START_REF] Ricker | The form and laws of propagation of seismic wavelets[END_REF]. Visco-elastic media is related to intrinsic attenuation. Different from elastic and acoustic media, the ratio between stress and strain is variable with time in visco-elastic media. To simulate the effects of intrinsic attenuation, there are mainly three Q model, including: Voigt-Ricker model [START_REF] Voigt | Ueber innere Reibung fester Korper, insbesondere der Metalle[END_REF][START_REF] Ricker | The form and laws of propagation of seismic wavelets[END_REF], the nearly constant Q model [START_REF] Kolsky | The propagation of stress pulses in viscoelastic solid[END_REF][START_REF] Futterman | Dispersive body wave[END_REF][START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF][START_REF] Day | Numerical simulation of wave fields using a Pade approximation method[END_REF] and the constant Q theory [START_REF] Strick | The determination of Q, dynamic viscosity and creep curves from wave propagation measurements[END_REF][START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF].

For Voigt-Ricker model, the attenuation effect is regarded as an energy absorption term, which is described by adding a single damping term to the acoustic wave equation in the model [START_REF] Ricker | The form and laws of propagation of seismic wavelets[END_REF]. Because of this simplicity, the theory of wave propagation with absorption has been widely adopted to compute the synthetic seismograms in early researches [START_REF] Boore | Comparison of two independent methods for the solution of wave-scattering problems: Response of a sedimentary basin to vertically incident SH wave[END_REF][START_REF] Munasinghe | Finite difference analysis of Rayleigh wave scattering at vertical discontinuities[END_REF][START_REF] Deng | True-amplitude prestack depth migration[END_REF]. However, the Q varies with frequency. The observations contradict the frequency-dependence of Q: other models have thus been developed.

In exploration seismology, within the specific frequency bandwidth, the nearly constant Q or constant Q theory means that the effect of attenuation of seismic waves for various frequencies are different (Knopoff, 1964;[START_REF] Toksoz | Attenuation of seismic waves in dry and saturated rocks: I. Laboratory measurements[END_REF][START_REF] Toksoz | Seismic Wave Attenuation[END_REF][START_REF] Barton | Rock quality, seismic velocity, attenuation and anisotropy[END_REF]. The attenuation of seismic waves can be defined by the attenuation factor or attenuation coefficient or absorption factor [START_REF] Aki | Quantitative seismology[END_REF][START_REF] Sheriff | Exploration seismology[END_REF][START_REF] Wang | Seismic Inverse Q Filtering[END_REF], given as

α = ω 2cQ , (1.3)
where c is the speed of waves. With constant Q assumption, the high frequency part attenuates faster than lower frequency one, which is consistent with many laboratory experiments [START_REF] Peselnick | Internal friction of fine grained limestones at ultrasonic frequencies[END_REF][START_REF] Peselnick | Internal friction in shear and shear modulus of Solenholfen limestone over a frequency range of 10 7 cycles per second[END_REF][START_REF] Toksoz | Attenuation of seismic waves in dry and saturated rocks: I. Laboratory measurements[END_REF] and field measurements [START_REF] Mcdonal | Attenuation of shear and compressional waves in Pierre shale[END_REF][START_REF] Tullos | Seismic attenuation of Gulf coast sediments[END_REF][START_REF] Hamilton | Compressional-wave attenuation in marine sediments[END_REF]. One example of field measurement of Q was conducted by [START_REF] Mcdonal | Attenuation of shear and compressional waves in Pierre shale[END_REF] in Pierre shale (Figure 1. 5). The attenuation factor α is linear within the 100-500 Hz frequency range, the equivalent Q is around 32.5. A constant or nearly constant Q model is widely used in exploration seismology. For nearly constant Q model, attenuation coefficients are determined such that Q is constant within the desired frequency bandwidth. [START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF] firstly connected the nearly constant Q model with a physically rheological model. The nearly constant Q is obtained by superimposing L pairs of Zener model. A classical Zener body is called standard linear solid (SLS). [START_REF] Day | Numerical simulation of wave fields using a Pade approximation method[END_REF] first incorporated realistic attenuation into the time-domain computation using the Padé approximation. Many authors developed the nearly constant Q models and visco-elastic or visco-acoustic wave equations based on different rheological models, mainly including generalized Maxwell body (GMB) [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]Carcione et al., 1988a,b;[START_REF] Emmerich | PSV-wave propagation in a medium with local heterogeneities: A hybrid formulation and its application[END_REF]Fah, 1992;[START_REF] Moczo | Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures[END_REF][START_REF] Kay | Applying finite element analysis to the memory variable formulation of wave propagation in anelastic media[END_REF] and generalized Zener body (GZB) [START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF][START_REF] Xu | Composite memory variables for viscoelastic synthetic seismograms[END_REF][START_REF] Robertsson | A numerical free-surface condition for elastic/ viscoelastic finite-difference modelling in the presence of topography[END_REF][START_REF] Hestholm | Composite memory variable velocity-stress viscoelastic modelling[END_REF]. [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF] prove both GMB and GZB are equivalent to implement realistic attenuation in the time-domain. The rheological models will be further discussed in Section 2.3.

The constant Q model corresponds to the description of wave propagation and attenuation with Q exactly independent of frequency [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF]. However, the purely constant Q model leads to a infinite phase velocity, as the frequency increases to infinity [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF]. This contradicts the physical reality. Therefore, the application of this model should be used within a limited frequency bandwidth, e.g, under the frequency range of exploration seismology. Some authors incorporated the constant Q theory in the visco-elastic or visco-acoustic wave equations [START_REF] Carcione | Time-domain modelling of constant-Q seismic waves using fractional derivatives[END_REF][START_REF] Carcione | A generalization of the Fourier pseudospectral method[END_REF][START_REF] Zhang | Compensating for visco-acoustic effects in reverse-time migration[END_REF]Zhu and Harris, 2014). [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] constant Q model involves fractional derivative wave equation in the time-domain. The fractional derivatives can be implemented in the wavenumber domain using extra space Fourier transform (Zhu and Harris, 2014).

Compared with constant Q, nearly constant Q needs several pairs of rheological models to obtain the Q being constant over a specific frequency range. The nearly constant Q model is widely applied in visco-acoustic / visco-elastic modeling due to a simple implementation in the time domain.

Estimation of attenuation

In exploration seismology, there are different kinds of methods to estimate Q. Mainly, they can be divided into ray-based methods under the high frequency approximation and wave equation methods. Ray-based methods can either use spectral amplitude loss or spectrum distortion to estimate the attenuation along the wave travel path. Wave equation methods rely on the full waveform inversion to update the Q. Independently of this classification, the estimation of attenuation can be applied in the time and frequency domains. Ray-based methods can both measure the wavelet broadening time in the time domain or estimate the spectrum changes in the frequency domain. Wave equation methods can also be implemented in the time and frequency domains. Furthermore, Q estimation also can be derived in data domain or image domain. The data domain methods update the Q model by fitting the calculated data and observed data. The image domain methods retrieve the Q model by focusing the image panels, with the construction in subsurface image domain. We review those three approaches in this section.

Ray-based methods

Seismic dissipation caused by attenuation, leads to spectrum amplitude loss and spectrum distortion of seismic data. Apart from the scattering, geometric spreading, source and receiver coupling, radiation patterns, and transmission/reflection effects, the recorded amplitude loss and spectrum distortion are proportional to the integration of attenuation over the distance between source and receiver positions. With a known velocity, this information links the attenuation and traveltime. Thus, the Q model can be estimated by the ray-based method. The ray theory is based on the high frequency asymptotic approximation. Traveltime and amplitude are solutions of the Eikonal and transport equations, respectively. In this context, the Q estimation by the usage of traveltime is mainly into two main categories, the amplitude method and spectrum distortion method.

Spectrum amplitude loss method [START_REF] Brzostowski | 3D tomographic imaging of near-surface seismic velocity and attenuation[END_REF][START_REF] Leggett | Study of traveltime and amplitude time-lapse tomography using physical model data[END_REF][START_REF] Cavalca | Estimation of the near-surface Q distribution from pre-stack surface seismic data[END_REF][START_REF] Cavalca | Raybased Tomography for Q Estimation and Q Compensation in Complex Media[END_REF] is first needed to estimate a velocity model, commonly by the standard traveltime tomography. Then the spectrum amplitude loss is estimated along the rays to obtain Q values. The relation between source and observed received amplitude along ray path [START_REF] Futterman | Dispersive body wave[END_REF][START_REF] Aki | Quantitative seismology[END_REF][START_REF] Bregman | Traveltime and amplitude inversion in crosshole seismology[END_REF][START_REF] Tonn | The determination of the seismic quality factor Q from VSP data: A comparison of different computational methods[END_REF] is

A k (f ) = A s (f ) exp[-πf ray k ds Q(x, y, z)v(x, y, z) ], (1.4)
where A s (f ) is the source amplitude, and A k (f ) is the observed amplitude along the ray k. The velocity model is obtained by a prior velocity tomography, then the Q model is determined by equation 1.4. However, the spectrum amplitude of seismic data is easily contaminated by many factors such as scattering, geometric spreading, source and receiver coupling, radiation patterns, and transmission/reflection effects. The assumption of this method is that geometrical spreading and scattering effects are first well removed.

The other way is to use the spectrum distortion method. Because of attenuation, the high frequency waves attenuate faster than low frequency waves, which leads to the downshift of centroid frequency. [START_REF] Quan | Seismic attenuation tomography using the frequency shift method[END_REF] first estimated Q based on centroid shift method. The advantage of the method is that only the attenuation is responsible for centroid frequency shift, compared with amplitude loss method contaminated by many other factors. However, this method is based on the assumption that the spectrum of incident waves is known. Based on centroid frequency method, Plessix (2006a) applied the method to estimate Q using limited frequency band data in crosswell case. [START_REF] Zhang | Estimation of quality factors from CMP records[END_REF] applied the method to obtain quality factors from common midpoint records using a layer-stripping approach.

Additionally, another way for estimating perturbed Q is to combine the ray theory and Born approximation (Ribodetti and Virieux, 2000;Ribodetti et al., 2000). The Geen's function has an asymptotic expression. The perturbed wavefields caused by the perturbed Q are obtained according to the Born approximation. The perturbed Q is reconstructed by minimizing the data misfit between perturbed wavefields and recorded data. This linearized inversion is fast compared with full waveform inversion, and is devoted to the situation of the short wavelength components of the attenuation model.

Wave equation methods

Full waveform inversion (FWI) exploits all the waveform information by minimizing the data difference [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations: Conference on Inverse Scattering, Theory and Application[END_REF]Tarantola, 1984a). It is a powerful method to retrieve physical models of the subsurface, including pressure and shear velocity as well as models. [START_REF] Tarantola | Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation[END_REF] extended the FWI algorithm to retrieve the elastic parameters and attenuation parameter. As the computation ability increases, the theory and applications of FWI have been widely developed, mainly after the 2000's with a better understanding of the role of low frequency in the inversion [START_REF] Pratt | Seismic waveform inversion in the frequency domain, Part I: Theory and verification in a physical scale model[END_REF]. At the same time, visco-elastic or visco-acoustic FWI also has been developed for more realistic physical models, with applications on both synthetic and real data cases. Visco-elastic or visco-acoustic FWI can be implemented in the time domain (Liao and McMechan, 1995;[START_REF] Charara | The state of affairs in inversion of seismic data: an OVSP example[END_REF][START_REF] Charara | Full waveform inversion of seismic data for a viscoelastic medium[END_REF][START_REF] Barnes | Borehole seismic data inversion for attenuating media[END_REF][START_REF] Charara | Crosswell seismic data tomography for heterogeneous viscoelastic media[END_REF]Barnes andCharara, 2009, 2010;[START_REF] Bai | Waveform inversion with attenuation[END_REF][START_REF] Kurzmann | Acoustic full waveform tomography in the presence of attenuation: A sensitivity analysis[END_REF][START_REF] Bai | Viscoacoustic waveform inversion of velocity structures in the time domain[END_REF][START_REF] Barnes | [END_REF][START_REF] Kurzmann | Viscoacoustic full waveform inversion for spatially correlated and uncorrelated problems in refection seismics[END_REF][START_REF] Belahi | Elastic versus viscoelastic full waveform inversion of near-offset and wide-angle data in the presence of attenuation[END_REF][START_REF] Cheng | A new approach of viscoacoustic waveform inversion in the time domain[END_REF][START_REF] Li | 2D time-domain visco-acoustic full waveform inversion method[END_REF][START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF][START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF][START_REF] Jiang | A modified 2D visco-acoustic waveform inversion scheme[END_REF][START_REF] Wang | Visco-acoustic fullwaveform inversion in the presence of complex gas clouds[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF] and frequency domain [START_REF] Song | Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part II-Inversion method, synthetic experiments and real-data results[END_REF][START_REF] Causse | Finite-difference inversion in viscoacoustic media[END_REF][START_REF] Liao | Multifrequency viscoacoustic modeling and inversion[END_REF][START_REF] Causse | Preconditioning of full-waveform inversion in viscoacoustic media[END_REF][START_REF] Hicks | Reflection waveform inversion using local descent methods: estimating attenuation and velocity over a gas-sand deposit[END_REF][START_REF] Pratt | Waveform tomography images of velocity and inelastic attenuation from the Mallik 2002 crosshole seismic surveys[END_REF][START_REF] Kamei | Waveform tomography strategies for imaging attenuation structure for cross-hole data[END_REF][START_REF] Malinowski | High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full waveform inversion[END_REF][START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF][START_REF] Operto | Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation[END_REF][START_REF] Gao | Simultaneous inversion for velocity and attenuation by waveform tomography[END_REF][START_REF] Operto | On the role of density and attenuation in threedimensional multiparameter viscoacoustic VTI frequency-domain FWI: an OBC case study from the North Sea[END_REF]) (more details in Chapter 2).

The advantage of FWI is that it applies for all the kinds of wave and exploits all the waveform information, with higher resolution compared with ray-based methods. However, it is a highly nonlinear problem, and it is hard to retrieve the velocity and Q due to different parameter sensitivities. Because of the low sensitivity to Q, some researches focus on the data representation to study when it is acceptable to neglect the Q model [START_REF] Kurzmann | Acoustic full waveform tomography in the presence of attenuation: A sensitivity analysis[END_REF][START_REF] Barnes | [END_REF][START_REF] Belahi | Elastic versus viscoelastic full waveform inversion of near-offset and wide-angle data in the presence of attenuation[END_REF]. Their results show it is worth to implement the visco-acoutic or visco-elastic FWI. Visco-elastic FWI leads to better results of elastic parameters than purely elastic FWI [START_REF] Barnes | [END_REF]. Furthermore, researches investigate the sensitivities of parameters and inversion strategies [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF][START_REF] Operto | Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation[END_REF][START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF][START_REF] Jiang | A modified 2D visco-acoustic waveform inversion scheme[END_REF] to better retrieve the velocity and Q models. The different inversion strategies recommended by different authors mainly include the sequential inversion [START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF], the simultaneous inversion [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF] and the hierarchical inversion [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF][START_REF] Operto | Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation[END_REF]. Sequential inversion first updates the velocity model, then continues to update the Q model. It can handle the different parameter sensitivities; however it may over-estimate the parameters [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF]. Simultaneous inversion updates the velocity and Q at same time, but suffers from the differences between parameter sensitivities without the estimation of the Hessian matrix [START_REF] Métivier | Full waveform inversion and the truncated Newton method[END_REF][START_REF] Métivier | Full waveform inversion and truncated Newton method: quantitative imaging of complex sub-surface structures[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. Hierarchical inversion first retrieves a good velocity model, then follows a subsequent simultaneous inversion. [START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF] applied the sequential inversion to estimate the Q factor, which is used to monitor the CO 2 injection in Cambrian-Ordovician sedimentary basin of the St. Lawrence Platform in southern Quebec, Canada. [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF] used the simultaneous inversion to obtain the high resolution velocity and Q models. They precisely detected the gas cloud in offshore Malaysia. [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF]; [START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF][START_REF] Operto | Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation[END_REF] conducted synthetic study for imaging density, velocity and Q and recommended the hierarchical inversion strategy. They applied this strategy to a real data from Valhall field in North sea. The seismic images and Q model reveal the precise location of gas clouds. All those successful real data examples show FWI is a powerful tool to obtain high resolution images of velocity and Q. Therefore, this thesis continues the research on velocity and Q inversion under the framework of FWI, to investigate a new inversion strategy to guide the inversion and mitigate cross-talks be-tween parameters without relying on the Hessian matrix (see motivations in Section 1.3 and methodology and results in Chapter 4 and 5).

Data domain and image domain methods

As mentioned in the introduction, the methods of estimating Q also can be divided into data [START_REF] Brzostowski | 3D tomographic imaging of near-surface seismic velocity and attenuation[END_REF][START_REF] Quan | Seismic attenuation tomography using the frequency shift method[END_REF][START_REF] Tarantola | Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation[END_REF] and image domains [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF][START_REF] Shen | Image-based Q tomography using reverse time Q migration[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF]. Data domain methods are based on the data fitting to update the Q model. Image domain methods retrieves the Q model by focusing the image panel [START_REF] Chauris | Analyse de vitesse par migration pour l'imagerie des structures complexes en sismique réflexion[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF][START_REF] Cocher | Seismic iterative migration velocity analysis: two strategies to update the velocity model[END_REF]Chauris and Cocher, 2017;Chauris et al., 2017;[START_REF] Li | Subsurface seismic imaging based on inversion velocity analysis in both image and data domains[END_REF]. Image domain can be implemented in common-shot and subsurface gathers [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF][START_REF] Symes | Extended Waveform Inversion[END_REF]. This is not a common approach for estimating Q [START_REF] Shen | Image-based Q tomography using reverse time Q migration[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF].

In data domain, ray-based methods (section 1.2.1) and wave-equation methods (section 1.2.2) measure the misfit of observed and calculate data to update Q model. For the ray-based methods, the observed data is exactly the spectrum amplitude or central frequency. For wave-equation methods, the observed data are the wavefields recorded at specific positions. The data can be all kinds or specific kinds of waves with whole or limited time window, e.g., direct waves, transmission waves, reflections.

The image domain methods back-project the image residuals to update the macro models [START_REF] Sava | Wave-equation migration velocity analysis[END_REF][START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF][START_REF] Sava | Extended imaging conditions for wave-equation migration[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF]. Inaccurate velocity or Q leads to a defocused energy of image panel (Figure 1.6). Based on this criteria, the macro Q model can be updated with a known fixed velocity [START_REF] Shen | Image-based Q tomography using reverse time Q migration[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF]. Data domain methods are straightforward. However, it is sensitive to the data quality, e.g., signal-to-noise ratio. Data domain methods, e.g., FWI can produce high resolution image. Image domain methods are based on Born approximation, dividing the model into a macro and perturbed model. Image domain method focusing on reconstructing the macro model. Methods in both domains suffer from the leakage between velocity and Q. In image domain, the image panel is focused when the background kinematic model is correct. For example, attenuation leads to dispersion phase velocity. Different velocity and Q can have same kinematic information. In this situation, the inversion of correct Q relies on prior known correct velocity model.

Motivation, strategy and thesis organization

In the previous sections, I introduce the attenuation effects. On the one side, attenuation information can solely be used as the lithology discriminator, e.g., detection of the gas clouds and saturated layers. On the other side, attenuation model can be applied to improve the velocity images. Attenuation can be represented as quality factor Q. Different Q models can be used to explain the attenuation effects. In the exploration seismology, the constant or nearly constant Q model is widely applied and consistent with field observations. Based on this Q model, ray-based and wave-equation methods are applied to retrieve velocity and Q models. Wave equation methods have the advantages of high resolution. The thesis focuses on the reconstruction of velocity and Q models under the frame of visco-acoustic FWI, without the estimation of the Hessian matrix. We encounter the well-known difficulties related to cross-talks between velocity and Q and different parameter sensitivities. Our objectives of the thesis are: (1) investigating the inversion strategies to handle different sensitivities of parameters, and (2) introducing extra information to guide the inversion and mitigate cross-talks between parameters.

Motivation: Inversion strategies to reconstruct velocity and attenuation

Visco-acoustic FWI like many other multi-parameter inversion suffers from cross-talks between the parameters [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF]. Additionally, due to the limited observed information on the surface, similar radiation patterns for velocity and Q [START_REF] Malinowski | High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full waveform inversion[END_REF], and different sensitivities of parameters, those difficulties lead to the insufficient resolution of velocity and Q of visco-acoustic FWI, referring to cross-talks. The cross-talks can be observed in the pattern of the gradient of objective function with respect to parameters (Figure 1.7). The lack of resolution of data disparity can also lead to cross-talks (Figure 1.8).

Velocity and Q have the different sensitivities due to physical properties. The changes of v p and Q can both generate a kinematic phase shift and amplitude changes. However, v p changes greatly impact the kinematic phase, leading to a large data difference. A scaling factor is needed to balance the magnitude. The incorporation of inverse of the Hessian matrix [START_REF] Métivier | Full waveform inversion and the truncated Newton method[END_REF][START_REF] Métivier | Full waveform inversion and truncated Newton method: quantitative imaging of complex sub-surface structures[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF] in the minimization process is one solution, but computation costs are very heavy, even with the truncated Newton approaches.

To mitigate cross-talks and balance the different gradient sensitivities, conventional inversion strategies include sequential inversion [START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF], simultaneous inversion [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF] and hierarchical inversion [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF][START_REF] Operto | Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation[END_REF]. Sequential inversion (Figure 1.9, left) usually first retrieves the velocity and then Q, but it may over-estimate the parameter [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF]. Simultaneous inversion (Figure 1.9, middle) reconstructs velocity and Q at the same time. It has to overcome the different sensitivities between parameters without the estimation of the Hessian matrix. Hierarchical inversion first retrieves velocity model then continues with a simultaneous inversion. There is not yet a clear statement about the optimal recommendation. Therefore, we investigate different inversion strategies and propose a hybrid inversion scheme without the estimation of Hessian matrix.

Strategy: Kinematic relationship to guide the inversion

As indicated above, one difficulty of visco-acoustic FWI is cross-talks between two parameters. In this sense, the inversion process of fitting the same observed data can lead to different inversion results. To improve the inversion, more information is needed to guide the inversion and reduce the searching model parameter space, without necessarily relying on the information of the Hessian matrix.

Attenuation effects lead to phase velocity dispersion [START_REF] Futterman | Dispersive body wave[END_REF][START_REF] Aki | Quantitative seismology[END_REF]. Phase velocity increases as frequency and attenuation increase. It means that different combinations of velocity and attenuation result in the same kinematic velocity v k for band-limited seismic data. We call this relation the kinematic relationship. The v k model controls the kinematic of wave propagation. It means the traveltime of recorded data is same, but the amplitudes are different. It is different from the P-velocity model v p because of attenuation.

We propose a two-step process (hybrid strategy) by incorporating the kinematic relationship (Figure 1.9 right): (1) For a fixed Q, we iteratively update v k . (more details in Chapter 4). ( 2) In a second-step, we update v p , but here for a fixed v k model. The gradient of the objective function is modified accordingly. During the second phase, Q is the passive parameter, i.e. is automatically updated when v p is modified. Unlike conventional inversion strategies with more flexibilities to search the parameters in whole model space, the model searching space of the hybrid inversion strategy in the secondstep is only limited to be along kinematic curve (Figure 1.10). Therefore, the key factor is to obtain the kinematic velocity for the inversion. We will discuss if it is easier to obtain the kinematic velocity for transmitted waves compared to reflection waves. We first apply this inversion strategy on models with mainly transmitted waves to obtain kinematic velocity. With this information, we then apply our hybrid inversion strategy to retrieve the velocity and Q. Those two steps are alternated with each other to further update the kinematic velocity, then decouple the two parameter based on updated kinematic velocity. We will discuss how this extra information can guide the inversion and mitigate cross-talks between parameters, especially on the model of Gaussian anomaly case in Chapter 4. Then we further extend this hybrid strategy to reflections dominated model, i.g., Marmousi model (Chapter 5).

Figure 1.9: Workflow of the sequential (left) simultaneous (middle) and hybrid (right) strategies. τ approximates to the inverse of Q. v k = f (v p , τ, ω c ) is the kinematic relationship, ω c is the central frequency of the source (more details in Section 2.6.3 and Chapter 4).

Thesis organization

The thesis is organized as followed: [START_REF] Komatitsch | An unspoilt convolutional perfectly matched layer improved at grazing incident for the seismic wave equation[END_REF][START_REF] Martin | An unspoilt convolutional perfectly matched layer technique improved at grazing incident for the viscoelastic wave equation[END_REF]. Then I compare the numerical seismogram with analytical results to evaluate the accuracy of the numerical solution. Finally, I analyze the visco-elastic effects. The visco-elastic waves have dissipation and dispersion. The dispersion effect lead to an equivalent kinematic velocity. The kinematic velocity will be used as an extra information to guide the inversion.

• In Chapter 3, the objective of this chapter is to better relate the different forms of visco-elastic wave equations mainly based on nearly constant, constant Q models and power law. Those models can be rewritten as a general relation between stress and strain with either ordinary or fractional derivatives. The general relationship between stress and strain can turn into different rheological models. Those rheological models include Generalized Maxwell body, Generalized Zener body and Scott-Blair model. Starting from the unification form, we can obtain the commonly used visco-acoustic wave equations in exploration geophysics. This unification form also can be used as a way to derive other forms of visco-elastic wave equation. Those equations are not limited to constant Q and can have complex relationship between Q and frequency. We demonstrate the feasibility of the modified inversion strategy. We also apply the frequency continuation inversion when higher frequencies are progressively included in the inversion scheme.

• In Chapter 6, I conclude the thesis and propose some directions for further study: Q-compensated reverse-time migration (RTM), Q least-squares RTM (LS-RTM), and migration velocity analysis (MVA) with Q. The three cases incorporate Q as a passive parameter to obtain macro velocity model or velocity perturbation based on linearized inversion.

Contributions

I reviewed visco-acoustic theories, forward modeling and inversion. I implemented the codes of 2D visco-acoustic wave simulation and inversion using MPI, a modified τ method to optimize nearly constant Q, forward modeling based on staggered grid finite difference (Virieux, 1986) with CPML [START_REF] Komatitsch | An unspoilt convolutional perfectly matched layer improved at grazing incident for the seismic wave equation[END_REF][START_REF] Martin | An unspoilt convolutional perfectly matched layer technique improved at grazing incident for the viscoelastic wave equation[END_REF], validation with analytical solution (Chapter 2), investigating a selfjoint back-propagation formula to obtain the gradient (Plessix, 2006b), stabilization for back-propagation, validation the gradient with comparison with FD method, conjugate gradient optimization and MPI implementation (part of Chapter 4). I investigated a unification formalism of visco-elastic wave equations, including fractional-order, stan-dard linear solids based and power-law based visco-elastic wave equations (Chapter 3).

Based on those modeling and inversion code packages, I propose a hybrid visco-acoustic FWI scheme, including kinematic velocity analysis, a total gradient by introducing the relation between Q and velocity, compare with different inversion strategies, and test on simple models, the synthetic Middle-East model (part of Chapter 4) and Marmousi model (Chapter 5). I draw some perspectives on Q migration velocity analysis and Qcompensated reverse-time migration (Chapter 6). Main contributions are:

• A modified τ method to optimize nearly constant Q in Chapter 2;

• A unification formalism to link different forms of the visco-elastic wave equations in Chapter 3;

• A hybrid visco-acoustic FWI strategy to guide the inversion by incorporating the kinematic relationship between velocity and Q in Chapter 4;

• A new term ∂τ (v p )/∂v p acting as a scaling factor to balance the gradient sensitivities in Chapter 4;

• Comparisons of different inversion strategies to better understand the visco-acoustic FWI for transmitted and reflected waves in Chapter 4 and Chapter 5.

Part of the results has been presented / accepted in

Introduction

This chapter is dedicated to introducing the numerical simulation of attenuation and visco-elastic wave equations based on constant or nearly constant Q model. Constant Q models are consistent with the field observations. In other word, the absorption coefficient is proportional to frequency, which means the Q is independent of frequency. High frequency waves attenuate faster than low frequency ones, leading to central frequency downshifts. This chapter first introduces the physical models to represent viscoelasticity, then extends to different constant Q models. Based on those Q models, different forms of visco-elastic or visco-acoustic wave equations are obtained. In the thesis, we use the visco-acoustic wave equation based on Generalized Zener Solids, implemented in time domain. Then this chapter introduces the numerical solution to simulate visco-acoustic waves, including finite-difference solver, numerical dispersion and stability criteria, and absorbing boundary condition. Finally, this chapter analyzes the visco-elastic effects in time and frequency domain, the dispersion curve, kinematic velocity and kinematic relationship.

Visco-elasticity

The linear visco-elasticity represents the phenomenon that the stress at any moment depends on the whole strain history. The stress-strain relationship or constitutive relationship can be built by a time-dependent function [START_REF] Boltzmann | Theorie der elastiche Nachwirkung[END_REF][START_REF] Christensen | Theory of Viscoelasticity[END_REF]. In contrast to visco-elastic media, the ratio of stress and strain in elastic media is invariable with time [START_REF] Fung | Foundations of solid mechanics[END_REF]. This ratio is called the elastic modulus. The simplest elastic media can be represented with a spring (Hooke element), with the relation σ(t) = M (t). tis the time, σ(t) means stress, (t) is strain, and M is the elastic modulus.

A dashpot is an damper, combined with a spring to exhibit visco-elastic behavior. The stress-strain relationship of a dashpot is σ(t) = η d (t) dt . η is the coefficient of viscosity. It is a linear relationship between the stress and the time derivative of the strain. One simple visco-elastic media is a Maxwell body; it consists of a spring and a dashpot in parallel. Different combinations of springs and dashpots compose Generalized Maxwell Body (GMB) and Generalized Zener Body (GZB). They are commonly used for viscoelastic media in exploration geophysics [START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF][START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]Carcione et al., 1988a,b;[START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF][START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF]. n pairs of Maxwell bodies and one spring in parallel make up the GMB. A spring in parallel with a Maxwell body form the one kind of Zener bodies, also called Stand Linear Solid (SLS). n pairs of Zener/SLS bodies in parallel are called Generalized Zener Body (GZB) or Generalized SLS (Figure 2.1).

To obtain the stress-stain relation for visco-elastic media, we start from the elastic media for comparison. The elastic stress-strain relationship is given as

σ(t) = M (t). (2.1) σ(t)
is the stress and (t) is the strain. M is the elastic modulus. It is constant and invariant with time.

For a linear isotropic visco-elastic media, the stress depends on the previous history of strain. The stress-strain relationship is obtained in a time convolutional form [START_REF] Boltzmann | Theorie der elastiche Nachwirkung[END_REF] given as

σ(t) = t -∞ ψ(t -τ ) d (τ ) dτ , = ψ(t) * ˙ (t), = ψ(t) * (t), (2.2)
where σ(t) and (t) are same as before, representing the stress and strain. ˙ (t) is the first order time derivative of strain. A time integration can be written into a convolution form, where * is the time convolutional operator. ψ(t) is the stress relaxation function, which is the stress response to a Heaviside step function of strain. It can be measured by laboratory experiments [START_REF] Zener | Elasticity and anelasticity of metals[END_REF][START_REF] Christensen | Theory of Viscoelasticity[END_REF][START_REF] Sasso | Application of fractional derivative models in linear viscoelastic problems[END_REF].

To compare with the elastic stress-strain relation, we write the third formula of equation 2.2 in frequency domain. The convolution in time turns into a product in frequency,

σ(ω) = M (ω) (ω), (2.3)
where M (ω) is denoted as the complex modulus of visco-elastic media. It is equal to

F[ ψ(t)],
where F is the Fourier transform operator. The transform pair of modulus and stress relaxation function is

M (ω) = F ψ(t) , (2.4
)

ψ(t) = F -1 M (ω) iω , (2.5)
where F -1 is the inverse Fourier transform operator. Frequency domain M (ω) is easy to be obtained by the analyzing the stress-strain relation based on rheological models, e.g. SLS. ψ(t) is therefore deduced by the inverse Fourier transform of M (ω)/iω. Taking the SLS as example in Figure 2.1, SLS is composed of a spring and a dashpot in parallel, then with another spring in series. Complex modulus in the frequency domain and stress relaxation function in time domain respectively, are given by

M (ω) = M R 1 + iωτ 1 + iωτ σ , (2.6)
ψ(t) = M R 1 -1 - τ τ σ exp (-t/τ σ ) H(t).
(2.7)

H(t) is the Heaviside step function. τ σ = η/(M 1 + M 2 ) is called
τ = η/(M 2 ) is called strain relaxation time, which is related to strain creep rate. M R is called relaxed modulus obtained when ω = 0 or t = ∞, defined as M 1 M 2 /(M 1 +M 2 )
. It means the visco-elastic media reaches the relaxed status after long time continuous stress imposing. When ω = ∞ or t = 0, the visco-elastic media is imposed by an instantaneous stress M U , which is called unrelaxed modulus equal to M R τ /τ σ . During the stress imposing history, the strain is unchanged (Figure 2.2). The complex modulus of visco-elasticity media varies with time. The complex modulus can be used for further discussion on Q models.

Constant Q approximation

The reciprocal of Q is defined as the proportion of energy loss per wave propagation cycle. It is obtained by ratio of imaginary part and real part of the visco-elastic modulus [START_REF] Boit | Theory of stress-strain relations in anisotropic visco-elasticity and relaxation phenomena[END_REF][START_REF] Knopoff | The seismic pulse in materials possessing solid friction[END_REF][START_REF] Fung | Foundations of solid mechanics[END_REF][START_REF] O'connell | Measures of dissipation in viscoelastic media[END_REF][START_REF] Ben-Menahem | Seismic waves and sources[END_REF], yielded by

1 Q(ω) = [M (ω)]/ [M (ω)]. (2.8) [M (ω)
] and [M (ω)] are related to the time-averaged strain-energy and dissipated energy, respectively. The constant Q models are consistent with the field observations. It is widely applied in exploration seismology. There are two kinds of constant Q models. The first one is based on the model proposed by [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF]. The visco-elastic wave equation is deduced with fractional derivative. The other one is called nearly constant Q model. It mainly is based on rheological models [START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF][START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]Carcione et al., 1988a,b;[START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF], e.g, GZB and GMB. The nearly constant Q model has to be approximated over a frequency range. The fractional derivative is hard to solve by a finite-difference method in the time domain. Thus, the latter nearly constant Q model is widely applied. We first introduce the constant Q model, and continue to discuss nearly constant Q model.

Kjartansson constant Q

Kjartansson (1979) Q model is independent with frequency. Based on laboratory experiments, it is reasonable to have relaxation function in a power law, similar to the shape displayed in Figure 2.2. Constant Q model is obtained by giving a specific stress relaxation function

ψ(t) = M 0 Γ(1 -2γ) ( t t 0 ) -2γ , 0 < 2γ < 1, (2.9)
where M 0 is the modulus, and Γ is the gamma function. For t < 0, relaxation function is zero. According to equation 2.5, the complex modulus M (ω) is the Fourier transform of the time derivative of ψ(t).

M (ω) = M 0 ( iω ω 0 ) 2γ = M 0 | ω ω 0 | 2γ e iπγ sgn(ω) , (2.10)
where ω 0 is the reference frequency. At reference frequency, all waves have same velocity. γ is a constant, with 0 < 2γ < 1. When ω > 0, sgn(ω) = 1. For ω < 0, sgn(ω) = -1.

Substituting the complex modulus into equation 2.8, the constant quality factor is

1 Q(ω)
= tan(πγ).

(2.11)

Q is independent of frequency and represented by the parameter γ. When γ trends 0, the visco-elasticity turns into elasticity. Oppositely, the viscoelastic media is highly attenuated when γ trends 0.5.

Based on [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] Q model, the scalar visco-acoustic wave equation is obtained as

∂ 2-2γ p ∂ 2-2γ t - M 0 ρ ω -2γ 0 ∇ 2 p = 0, 0 < 2γ < 1, (2.12)
where p is the pressure wavefields. The wave equation involves a fractional derivative.

It is hard to solve by finite-difference in the time domain. γ → 0 turns the time derivative of pressure wavefields into ∂ 2 p/∂ 2 t, the equation becomes an elastic propagation function. Meanwhile, γ → 0.5, the time derivative of pressure wavefields is a diffusion term ∂p/∂t. It indicates the fractional visco-acoustic wave-equation is between the propagation and diffusion equations. Constant Q is easily achieved by using a single γ over all the frequency band. However, it not a realistic due to the infinitely large modulus when frequency increases to infinity (Equation 2.10). The application should be limited to a finite frequency range. Nearly constant Q models is optimized to find best value to have Q constant within a frequency bandwidth.

Nearly Constant Q model

Due to the difficulty to solve fractional derivatives, alternative nearly Q models are widely used. The nearly constant Q can be estimated based on rheological models [START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF][START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]Carcione et al., 1988a,b;[START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF], e.g., GMB and GZM. Analyzing the stress-strain relation of GMB or GZB, the complex modulus M (ω) is obtained (Equation 2.6). As before, Q is equal to the ratio of real and imaginary term of M (ω). The nearly constant Q is approximated within a frequency bandwidth. We first introduce the Q models based on GZB and GMB, then show the equivalence between the two approaches.

The modulus of a Zener body / SLS is given in equation 2.6. The total modulus of a GZM is the sum of modulus of L pairs of Zener body / SLS in parallel [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF][START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF], given by

M (ω) = L l=1 M Rl 1 + iωτ l 1 + iωτ σl , (2.13) = M R L l=1 1 L 1 + iωτ l 1 + iωτ σl , (2.14)
where M Rl is the elastic modulus of spring inside each Zener body. To be modified in accordance to GMB, M Rl is defined as 1 L M R , where M R is the total relaxed modulus of a GZB [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF]. τ σl and τ l are the stress and strain relaxation time inside each Zener body [START_REF] Pipkin | Lectures on viscoelasticity theory[END_REF]. According to the relation between M (ω) and relaxation time function ψ(t) in equation 2.5, ψ(t) is obtained as

ψ(t) = M R (1 - 1 L L l=1 (1 - τ l τ σl )e -t/τ σl )H(t), (2.15)
The relaxation function is used for constructing the wave equations in next section 2.4.

The ratio between the real and imaginary part of complex modulus is defined as Q (equation 2.8). The Q based on a GZB model is given as

1 Q(ω) = [M (ω)]/ [M (ω)], = L l=1 ω(τ l -τ σl ) 1+ω 2 τ 2 σl L l=1 1+ω 2 τ l τ σl 1+ω 2 τ 2 σl (2.16)
The nearly constant Q model can also be based on GMB [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]. A Maxwell body is composed of a spring and a dashpot in parallel. The GMB is obtained by L pairs of Maxwell bodies and a spring in parallel. The complex modulus is given by

M (ω) = M R + L l=1 M l iωτ σ l 1 + iωτ σ l , (2.17)
M R is the relaxed modulus, is related to single parallel spring inside the GMB. τ σ l = η l /M l = 1/ω l is called stress relaxation time of Maxwell body. M l is the elastic modulus of each spring and η l is the viscosity inside the same Maxwell body. Therefore, the qualify factor is achieved as

1 Q(ω) = L l=1 M l M R ωτ σl 1+ω 2 τ 2 σl 1 + L l=1 M l M R ω 2 τ 2 σl 1+ω 2 τ 2 σl (2.18)
Q formulas based on GZB (Equation 2.16) and GMB (Equation 2.18) may appear to be different. By defining

α l = 1 L τ l τ σl -1 (2.19)
for Equation 2.16, and

α l = M l /M R (2.20)
for Equation 2.18, an equivalent Q formula is derived as

1 Q(ω) = L l=1 α l ωτ σl 1+ω 2 τ 2 σl 1 + L l=1 α l ω 2 τ 2 σl 1+ω 2 τ 2 σl . (2.21)
α l is a positive parameter to determine the constant Q value within a frequency bandwidth [ω min , ω max ]. τ σl is determined in advance, empirically given by 1/ω l , where ω l is chosen evenly over the frequency bandwidth in a logarithm scale to optimize the nearly constant Q.

The same Q format based on GZB and GMB turns into a rheology-independent Q model [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF]. However, rheology models are still useful to understand and analyze the stress-strain relationship and visco-elastic behaviors. The term α l in total number L is optimized to have nearly constant Q values. Generally, when L is larger, the constant Q trends to be more flat within the frequency bandwidth. However, L number of α l leads to computationally expensive numerical schemes. Taking L equal to 1 is not enough to approximate constant Q. We will introduce the τ method [START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF] with a single value (Section 2.3.3).

Modified optimization of nearly constant Q

The nearly constant Q is approximated within a frequency bandwidth. The τ method [START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF] with a single term is convenient for inversion, first assuming all α l in Equation 2.21 are replaced by a single τ , turning into

1 Q(ω) = L l=1 τ ωτ σl 1+ω 2 τ 2 σl 1 + L l=1 τ ω 2 τ 2 σl 1+ω 2 τ 2 σl , (2.22)
where τ is the parameter for representing Q. To achieve the nearly constant Q, the least-squares objective function of determining optimal τ , given by

J[τ ] = ωmax ω min (Q -1 (ω, τ ) -Q -1 0 ) 2 dω, (2.23)
where ω min and ω max are upper and bottom bounds of the frequency range. Q 0 is the desired value. To minimize the misfit function, the τ method finds the optimal τ by setting the derivative of J[τ ] with respect to τ as zero. To obtain a fast analytical solution, the second assumption of the τ method is τ 1, leading to 1 + τ ≈ 1. Under this approximation, the Equation 2.21 turns into

1 Q(ω) ≈ L l=1 τ ωτ σl 1 + ω 2 τ 2 σl .
(2.24)

Finally, the τ is determined by a analytical solution,

τ = 1 Q 0 ωmax ω min F (ω)dω ωmax ω min F (ω) 2 dω , (2.25)
where

F (ω) is defined as L l=1 ωτ σl 1+ω 2 τ 2 σl .
However, the τ method over-estimates the Q value. Because the approximation of 1 + τ ≈ 1, for same τ value, the Q defined in Equation 2.24 is smaller than the Q defined in Equation 2.22. Equation 2.22 is the actual Q formula and connected to visco-elastic wave equation. To solve the over-estimation, the way is to give a lower pre-determined value Q pre , other than the desired Q 0 to approximate the constant Q. For example, the desired quality factor is 15, the Q pre may be around 13 or 12 (Figure 2.3).

Therefore, we propose a modified τ method based on grid search for Equation 2.23, to find more precise τ to simulate constant Q within the desired frequency bandwidth. In practice, we search the τ values from 0 to 1 with a fine interval 0.00005. Each constant Q value corresponds to a best τ value. Here, we present some results for comparison between the τ method and the modified method (Figure 2.4 and 2.5). The modified method is more accurate and without estimating a empirical initial Q pre . 5 SLS can simulate a constant Q over 1.5 Hz to 100 Hz, and we use L = 5 for our later viscoacoustic wave equation. The τ method [START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF] is widely applied due to its simplicity with analytical solution. However, this method over-estimates the desired Q value. Our modified τ method is more precise without estimating a lower Q pre to approach the desired Q 0 . Based on this modified method, we obtain nearly Q model and the associated viscoacoustic wave equation. This is a fast procedure compared with the computation of the wave equations.

Visco-acoustic wave equations

Visco-elastic / visco-acoustic wave equations have various formats. We split the wave equations into three categories corresponding to Q models in Section 1.1.3, including Voigt-Ricker model, constant Q model and nearly constant Q model.

First, based on the second-order acoustic wave equation, a first-order time derivative of pressure wavefields describes the energy dissipation [START_REF] Ricker | The form and laws of propagation of seismic wavelets[END_REF]. Because of this simplicity, this wave equation theory with absorption has been widely adopted to compute the synthetic seismograms in early researches [START_REF] Boore | Comparison of two independent methods for the solution of wave-scattering problems: Response of a sedimentary basin to vertically incident SH wave[END_REF][START_REF] Munasinghe | Finite difference analysis of Rayleigh wave scattering at vertical discontinuities[END_REF][START_REF] Deng | True-amplitude prestack depth migration[END_REF]. However, the Q varies with frequency, which contradicts observations. Under the seismic frequency bandwidth, a constant Q is realistic. Based on constant Q [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF], the wave equation has to solve a fractional derivative in time domain, leading to few applications (Section 2.3.1).

Based on the nearly constant Q model (Section 2.3.2), visco-elastic / visco-acoustic wave equations are obtained by incorporating the complex modulus or complex velocity in frequency domain [START_REF] Aki | Quantitative seismology[END_REF][START_REF] Song | Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part II-Inversion method, synthetic experiments and real-data results[END_REF][START_REF] Liao | Multifrequency viscoacoustic modeling and inversion[END_REF][START_REF] Stekl | Accurate viscoelastic modeling by frequency-domain finite difference using rotated operators[END_REF][START_REF] Hicks | Reflection waveform inversion using local descent methods: estimating attenuation and velocity over a gas-sand deposit[END_REF][START_REF] Pratt | Waveform tomography images of velocity and inelastic attenuation from the Mallik 2002 crosshole seismic surveys[END_REF]Operto et al., 2007b;[START_REF] Brossier | Frequency-domain numerical modelling of visco-acoustic waves based on finite-difference and finite-element discontinuous Galerkin methods[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Gosselin-Cliche | 3D frequency-domain finite-difference viscoelastic-wave modeling using weighted average 27-point operators with optimal coefficients[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF][START_REF] Operto | Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation[END_REF] or time convolution in time domain [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]Carcione et al., 1988a,b;[START_REF] Emmerich | PSV-wave propagation in a medium with local heterogeneities: A hybrid formulation and its application[END_REF]Fah, 1992;[START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF][START_REF] Xu | Composite memory variables for viscoelastic synthetic seismograms[END_REF][START_REF] Robertsson | A numerical free-surface condition for elastic/ viscoelastic finite-difference modelling in the presence of topography[END_REF][START_REF] Moczo | Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures[END_REF][START_REF] Kay | Applying finite element analysis to the memory variable formulation of wave propagation in anelastic media[END_REF][START_REF] Hestholm | Composite memory variable velocity-stress viscoelastic modelling[END_REF]. Frequency domain methods are easy to implement by incorporating complex modulus or complex velocity. However, modeling of one frequency wave with a finite-difference method requires solving a large, sparse system of linear equations (Operto et al., 2007b;[START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF], especially for 3D applications. The time domain methods incorporate the time convolution into equations. The stress is presented by visco-elastic or visco-acoustic modulus convoluting with strain. Based on nearly constant Q model and the methods of optimizing constant Q, the viscoelastic modulus is written as nth-order rational function, turning time convolution into nth-order differential equation.

In this thesis, we apply the visco-elastic / visco-acoustic wave equation based on nearly constant Q model. The nearly Q is approximated by modified τ method, representing the attenuation parameter (Section 2.3.3).

First, we start from the stress relaxation function based on GZB (Equation 2.15), replaced with parameter τ , the stress relaxation functions related to P and S waves in a isotropic visco-elastic media (Carcione et al., 1988a,b;[START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF] are given by

ψ P (t) = (λ + 2µ)(1 + L l=1 τ P e -t/τ σl )H(t), (2.26)
ψ S (t) = µ(1 + L l=1
(1 + τ S e -t/τ σl )H(t), (2.27) where ψ P (t) and ψ S (t) are the stress relaxation functions related to P and S waves, respectively. As before, H(t) is the Heaviside function. τ P and τ S are attenuation parameters for P and S waves. λ and µ are the Lamé parameters. λ + 2µ is the relaxed modulus of P waves and µ is the relaxed modulus of S waves. l means the lth SLS, and total number is L to simulate constant Q.

For an isotropic elastic media, the constitutive relation is σ ij = λδ ij kk + 2µ ij . The index i, j, k denote x, y in 2D media and x, y, z in 3D media. According the convolution relation between stress and strain in a linear visco-elastic media in Equation 2.2, the stress-strain in isotropic visco-elastic media is

σ ij = ψλ * δ ij kk + 2 ψµ * ij , (2.28)
where ψ denotes the time derivative of ψ, * denotes the time convolution, σ is stress and is strain. The term ψ λ and ψ µ are defined as

ψ µ = ψ S , (2.29) ψ λ = ψ P -2ψ S .
(2.30)

The strain-displacement relationship is

ij = 1 2 (∂ i u j + ∂ j u i ), (2.31)
where u is displacement. Replacing the displacement u with velocity of particles, the time derivative of strain is

˙ ij = 1 2 (∂ i v j + ∂ j v i ), (2.32)
v is velocity of particle, equal to u. Substituting Equation 2.32 into time derivative of Equation 2.28, the relation between stress and velocity of particles is

σij = ( ψP -2 ψS ) * ∂ k v k + 2 ψS * ∂ i v j , i = j, (2.33) σij = ψS * (∂ i v j + ∂ j v i ), i = j, (2.34)
Introducing the relaxation function (Equation 2.27) into the convolution of Equation 2.33 and 2.34, the convolution equation turns into a ordinary differential form

σij = (λ + 2µ) (1 + Lτ P ) -2µ 1 + Lτ S ∂ k v k + 2µ 1 + Lτ S ∂ i v j (2.35) + L l=1 r l ij , i = j, σij = µ(1 + Lτ S )(∂ i v j + ∂ j v i ) + L l=1 r l ij , i = j, (2.36)
with memory variables r l ij , denoted as

ṙ l ij = - 1 τ σl (r l ij + (λ + 2µ)τ P -2µτ S ∂ k v k + 2µτ S ∂ i v j , i = j, (2.37) ṙ l ij = - 1 τ σl (r l ij + µτ S (∂ i v j + ∂ j v i )), i = j, (2.38)
According to the Newton's second law, the relation between velocity of particles and

stress is vi = 1 ρ ∂ j σ ij . (2.39)
Finally, the visco-elastic wave-equation in stress-velocity format is composed by Equations 2.36, 2.36, 2.37, 2.38 and 2.39. If τ P = τ S = 0, the visco-elastic wave equation becomes the elastic one. Setting µ and τ S as zero, the visco-acoustic wave-equation is

∂p ∂t + κ(Lτ + 1)(∇ • v) + L l=1 r l p = f (x s , t), (2.40) ∂v ∂t + 1 ρ ∇p = 0, (2.41) ∂r l p ∂t + 1 τ l σ (r l p + τ κ(∇ • v)) = 0. l =∈ [1, L] (2.42)
p is the pressure wavefield, v is the particle velocity. κ = ρv 2 p denotes for relaxed bulk modulus. This modulus is obtained as the limit of the visco-acoustic complex modulus when t → ∞ or ω → 0 (Section 2.2 and Figure 2.2). τ is the term to represent Q model. The visco-acoustic wave equation has two parameters v p for P wave velocity and τ for attenuation. Based on the visco-acoustic wave equation, we introduce the numerical solution.

Numerical soltuion

Wave equations can be solved numerically in many schemes, including the finite difference method (FDM) [START_REF] Alterman | Propagation of elastic waves in layered media by finite difference methods[END_REF][START_REF] Kelly | Synthetic seismograms: A finite-difference approach[END_REF]Virieux, 1986) and the finite element method (FEM) [START_REF] Smith | The application of finite element analysis to body wave propagation problems[END_REF][START_REF] Marfurt | Accuracy of finite-difference and finite-element modelling of the scalar and elastic wave equations[END_REF]. FEM is more flexible to deal with complex geological body interfaces, as well as with the topography. Among the FEM family, spectral element method (SEM) [START_REF] Komatitsch | Méthodes spectrales et elements spectraux pour l'équation de l'élastodynamics 2D et 3D en milieu hétérogène[END_REF][START_REF] Komatitsch | The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF] combines the flexibility of FEM with the accuracy of the analytical / semi-analytical spectral method, which is gaining the popularities in global seismology and strong motion seismology [START_REF] Komatitsch | Méthodes spectrales et elements spectraux pour l'équation de l'élastodynamics 2D et 3D en milieu hétérogène[END_REF]Tromp, 1999, 2002;[START_REF] Chaljub | Spectral element analysis in seismology[END_REF]. However, FEM is expensive in terms of computational costs. FDM provides a relatively more efficient and cheaper implementation, as a compromise between the accuracy and computational costs. It is widely adopted in exploration seismology, can be implemented in the time (Virieux, 1986;[START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF][START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[END_REF] or in the frequency [START_REF] Pratt | Frequency domain elastic wave modeling by finite differences: A tool for crosshole seismic imaging[END_REF][START_REF] Jo | An optimal 9 point, finite difference, frequency space, 2-D scalar wave extrapolator[END_REF][START_REF] Pratt | Seismic waveform inversion in the frequency domain, Part I: Theory and verification in a physical scale model[END_REF] domains. The efficiency and complexity of the different numerical modeling methods are reviewed by [START_REF] Virieux | A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging[END_REF]. In this thesis, I use FDM to simulate the visco-acoustic in the time domain.

Finite difference solver

FDM is a natural way to solve partial differential equations. The method directly estimates each differential term using the Taylor expansion on a regular grid. The numerical solutions by FDM for visco-elastic wave equtions is more straightforward to solve in frequency domain by incorporating a complex and frequency-dependent viscoelastic modulus in the constitutive relation [START_REF] Stekl | Accurate viscoelastic modeling by frequency-domain finite difference using rotated operators[END_REF]Operto et al., 2007a,b;[START_REF] Brossier | Frequency-domain numerical modelling of visco-acoustic waves based on finite-difference and finite-element discontinuous Galerkin methods[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF]. However, modeling of one frequency wave with a finite-difference method requires solving a large, sparse system of linear equations (Operto et al., 2007b;[START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. A sparse direct solver needs large memories to store the matrix during LU decomposition, especially for 3D appliations. Therefore, visco-elastic wave equation can be solved in time domain by efficiently discretizing the convolution integral constitutive equations [START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF][START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF][START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF][START_REF] Carcione | Time-domain modelling of constant-Q seismic waves using fractional derivatives[END_REF]Saenger andBohlen, 2002, 2004).

The 2D spatial samplings are ∆x and ∆z and the time step is ∆t, the first-order derivatives in temporal coordinates and in the spatial coordinates (lateral position x and depth position z) and are both approximated by second-order approximation (Figure 2.6), given as

∂v x ∂t | k i,j = v k+ 1 2 x i,j -v k-1 2 x i,j ∆t , ∂v z ∂t | k i+ 1 2 ,j+ 1 2 = v k+ 1 2 z i+ 1 2 ,j+ 1 2 -v k-1 2 z i+ 1 2 ,j+ 1 2 ∆t , ∂p ∂t | k+ 1 2 i+ 1 2 ,j = p k+1 i+ 1 2 ,j -p k i+ 1 2 ,j ∆t , (2.43)
and

∂v x ∂x | k+ 1 2 i+ 1 2 ,j = v k+ 1 2 x i+1,j -v k+ 1 2 x i,j ∆x , ∂v z ∂z | k+ 1 2 i+ 1 2 ,j = v k+ 1 2 z i+ 1 2 ,j+ 1 2 -v k+ 1 2 z i+ 1 2 ,j-1 2 ∆z , ∂p ∂x | k i,j = p k i+ 1 2 ,j -p k i-1 2 ,j ∆x , ∂p ∂z | k i+ 1 2 ,j+ 1 2 = p k i+ 1 2 ,j+1 -p k i+ 1 2 ,j ∆z , (2.44)
where i, j are the 2D spatial coordinates and k is the temporal coordinate. ± 1 2 means that the position of the grid is shifted by half a grid point, corresponding to staggered grids (Virieux, 1986). The first-order velocity-pressure wave equation is solved under the staggered grid scheme. v x ,v z and p are at the grids, shifted by half a grid point with each other (Figure 2.6). The memory variables r l p are related to p, and have the same positions as p.

Substituting the Taylor expansion in the visco-acoustic wave equation, the components v x , v z , p and r pl are updated as

v k+ 1 2 x i,j = v k-1 2 x i,j - 1 ρ ∆t ∆x (p k i+ 1 2 ,j -p k i-1 2 ,j ), (i, j) . (i + 1, j) . (i, j + 1) z x v x , ρ v z , ρ p, r l p , k(Lτ + 1)
Figure 2.6: Discretization of the media on a staggered grid, p is the pressure wavefields, v x and v z are particle velocity components.

v k+ 1 2 z i+ 1 2 ,j+ 1 2 = v k-1 2 z i+ 1 2 ,j+ 1 2 - 1 ρ ∆t ∆z (p k i+ 1 2 ,j+1 -p k i+ 1 2 ,j ), p k+1 i+ 1 2 ,j = p k i+ 1 2 ,j -κ(Lτ + 1) ∆t ∆x (v k+ 1 2 x i+1,j -v k+ 1 2 x i,j ) + ∆t ∆z (v k+ 1 2 z i+ 1 2 ,j+ 1 2 -v k+ 1 2 z i+ 1 2 ,j-1 2 ) -∆t L l=1 r k pl i+ 1 2 ,j , r k+1 pl i+ 1 2 ,j = r k pl i+ 1 2 ,j - ∆t τ σl r k pl i+ 1 2 ,j - κτ τ σl ∆t ∆x (v k+ 1 2 x i+1,j -v k+ 1 2 x i,j ) - κτ τ σl ∆t ∆z (v k+ 1 2 z i+ 1 2 ,j+ 1 2 -v k+ 1 2 z i+ 1 2 ,j-1 2 ), l =∈ [1, L] (2.45)
where the physical parameters κ, τ and ρ have the same space positions as the updating component.

The last equation for updating r pl has an explicit form. The implicit form is obtained by replacing r k pl i+ 1 2 ,j , leading to

r k+1 pl i+ 1 2 ,j = (1 - 1 2 ∆t τ σl )r k pl i+ 1 2 ,j - κτ τ σl ∆t ∆x (v k+ 1 2 x i+1,j -v k+ 1 2 x i,j ) - κτ τ σl ∆t ∆z (v k+ 1 2 z i+ 1 2 ,j+ 1 2 -v k+ 1 2 z i+ 1 2 ,j-1 2 ) /(1 + 1 2 ∆t τ σl ). (2.46)
With the staggered grid discretization, all the components of visco-acoustic wave equa-tion are updated inside the time loop.The spatial discretization of FDM is usually implemented with a regular rectangle grid. Irregular boundaries and shapes limit the application and resolution of modeling by FDM, especially dealing with topography. Finite element methods [START_REF] Marfurt | Accuracy of finite-difference and finite-element modelling of the scalar and elastic wave equations[END_REF] is more flexible to solve irregular boundaries. Due to the computation efficiency reason and for regular boundary, I apply FDM for forward modeling.

Stability and dispersion

FDM is discretized by finite steps. The spatial steps and temporal step should satisfy the Courant-Friedrichs-Lewy (CFL) condition to meet the numerical stability [START_REF] Courtant | Uber die partiellen differenzengleichungen der mathematischen physik[END_REF]. The theoretical analysis assuming a homogeneous media (Virieux, 1986;[START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF][START_REF] Saenger | Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[END_REF][START_REF] Martin | An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave propagation in poroelastic media[END_REF] gives the following condition for a second-order approximation both in spatial and temporal space

max(v p )∆t 1 ∆x 2 + 1 ∆y 2 + 1 ∆z 2 ≤ 1. (2.47)
Similar to the analog elastic schemes to satisfy the stability criteria, maximum v p has to be adjusted to the highest phase velocity along the frequency bandwidth in viscoacoustic media [START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Igel | Anisotropic wave propagation through finitedifference grids[END_REF][START_REF] Saenger | Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[END_REF].

In practice, the maximum values of frequency of source and attenuation leads to highest phase velocity. For 2D simulation with ∆x = ∆z = ∆h, stability condition becomes

max(v p ) ∆t ∆h ≤ 1 √ 2 .
(2.48)

If the criteria is not satisfied, numerical simulation is unstable generating results as shown in Figure 2.7b. For purely elastic media, body waves are not dispersive waves. However, finite difference discretization generates numerical dispersion during the simulation of wave propagation. The spatial sampling ∆h should be defined such that

∆h ≤ λ min n λ = v min f max 1 n λ , (2.49)
where λ min is shortest wavelength depending on the maximum frequency f max and the minimum velocity v min . n λ is grid point number per wavelength, it is related to spatial and temporal orders in the Taylor expansion. n λ = 10 or more is adequate for second-order approximation in space and time, versus n λ = 5 or more for fourth-order in space and second-order in time [START_REF] Alford | Accuracy of finite-difference modeling of the acoustic wave equation[END_REF], within a propagation up to 100 wavelengths through the media [START_REF] Sei | Dispersion analysis of numerical wave propagation and its computational consequences[END_REF]. It the criterion is not satisfied, the simulation generates numerical dispersion (Figure 2.7c). One important point mentioned is that numerical dispersion is different with the physical dispersion of viscoacoustic or visco-elastic waves.

Boundary condition

Different from the real seismic propagation, the numerical simulation encounters with a finite space due to the limited computation capacities. When waves reach the boundaries of the model, the waves will reflect back from the rigid boundaries, which generates the non-physical waves. To solve the problem, absorbing boundaries should be adopted to eliminate the artificial reflected waves. The intuitive way is adding extra layers around the computational space as sponges to absorb incident waves. It is effective with normal incident waves. Non-physical spurious reflections appear from other incident angles [START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF].

A new kind of absorbing boundary called perfectly matched layer (PML) was firstly applied to electromagnetic wave-equation to absorb the boundary reflections [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. PML is widely used in the seismic wave equation [START_REF] Chew | Perfectly Matched Layers for elastodynamics: a new absorbing boundary condition[END_REF][START_REF] Collino | Applications of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Komatitsch | A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[END_REF]. Furthermore, the behavior of PML at grazing incidence can be improved using an unsplit convolutional approach called Convolutional PML (CPML) proposed by [START_REF] Komatitsch | An unspoilt convolutional perfectly matched layer improved at grazing incident for the seismic wave equation[END_REF]. [START_REF] Martin | An unspoilt convolutional perfectly matched layer technique improved at grazing incident for the viscoelastic wave equation[END_REF] introduced CPML to visco-elastic wave equation. The principal is to damp the wave energy until it is weak enough, by introducing the complex spatial derivatives inside the CMPL. We implemented CMPL for numerical simulation with limited model extension. For example, a homogeneous model has 941 grid points in length and 941 grid points in depth, with the grid size ∆h = 5 m. The simulation time is 2 seconds with ∆t = 0.001 s. The CPML has 30 grid points with ∆h = 5 m around four each side. The central frequency for the Ricker source is 10 Hz. Spurious reflections appear when the wave arrives at the boundary without CPML (Figure 2.8). With CMPL boundary for visco-acoustic wave equation modeling, boundary reflections are absorbed (Figure 2.9). The CMPL works well and absorbs the boundary reflections. In this thesis, we put CMPL around each side of the model for simulation. We now compare the numerical solution with analytical solution to validate accuracy our FDM program.

Analytical and numerical comparison

To validate the accuracy of the finite difference solution, we compare the analytical solution with the staggered grid finite difference method in the pure acoustic media, by setting visco-acoustic parameter τ = 0. The 2D homogeneous acoustic analytical solution [START_REF] Aki | Quantitative seismology[END_REF] is

G 2D (x, t) = 1 2π H(t -r c ) t 2 -r 2 c 2 , (2.50)
where x is spatial coordinate, r is the distance between source and recorded location, c is the velocity. In frequency domain, it becomes

G 2D (x, ω) = i 4 H 1 0 (ω r c ), (2.51)
where H is the Heaviside function and H 1 0 is the Hankel's function of order 1. According to the representation theorem [START_REF] Aki | Quantitative seismology[END_REF], the seismic trace is obtained as the convolution of source term with the Green's function. We use a Ricker wavelet as source term with central frequency 20 Hz. To compare the analytical solution with the numerical solution, the source term for first order velocity-pressure wave equation is the time integration of Ricker wavelet. The acoustic velocity is 2000 m/s. The spatial grid size is 5 m in both horizontal and vertical directions. Two solutions are compared at distance 200 m and 800 m (Figure 2.10), corresponding to 5 and 20 wavelengths for the maximum frequency waves, respectively. The maximum relative difference is 2%. It validates the accuracy of the finite difference method for solving the wave equation in the acoustic case. As for the visco-elastic case, the analytical solution is obtained by introducing the complex modulus.

Visco-acoustic effects

The seismic attenuation effects include the dissipation and dispersion, leading to the wave energy loss, central frequency downshifting and kinematic changes. We first illustrate the spectrum distortion of 2D propagation of acoustic waves. After clarifying the spectrum distortion caused by 2D propagation, then we analyze the visco-elastic effects. Finally, we rewrite the visco-acoustic wave equation in the frequency domain, with a complex velocity. The real term of the complex velocity is related to the kinematic information, defined as kinematic velocity.

Distortion of 2D propagation

Seismic attenuation causes the spectrum distortion of the visco-elastic waves. The recorded data has lower central frequency than the source wavelet. However, the spectrum distortion also occurs for the acoustic waves in 2D propagation. It is necessary to check the acoustic propagation distortion in 2D media by both numerical and analytical solutions. We use a central 20 Hz Ricker wavelet as the source time function

S(t) = (1 -2π 2 f 2 c (t -t 0 ) 2 ) exp (-π 2 f 2 c (t -t 0 ) 2 ), (2.52)
where t is the time, t 0 is the time shift, and f c is the central frequency, equal to 20 Hz. The analytical solution is obtained by the source term convoluting with 2D Green's function (Equation 2.50). In 2D media, because of Hankel's transform, the central frequency is not preserved. The central frequency of a recorded trace is 17.5 Hz computed by analytical, and also by numerical solutions, lower than the central frequency of the source wavelet (Figure 2.12). The central frequency is invariant for 3D acoustic propagation. Therefore, one should pay attention to the central frequency downshifting of visco-acoustic waves. The phenomenon is caused by the effects of both attenuation and 2D geometric spreading. This kind of frequency downshifting should be compensated in 2D media for central frequency shift method for estimating attenuation. Figure 2.12: Spectrum of the recorded trance by analytical and numerical solutions, the central frequency of trace is 17.5 Hz computed by analytical solution (solid red) and finite difference method (dashed blue), lower than the source central frequency with 20 Hz.

Visco-acoustic effects in the time and frequency domains

To analyze the visco-elastic effects, we simulate one trace of acoustic and visco-acoustic waves at a distance of 400 m from the source (Figure 2.13). The visco-acoustic model has Q = 15 and same velocity as the last subsection. The central frequency of the source is 20 Hz. Because of 2D propagation distortion, the central of the acoustic wave at receiver position is 17.5 Hz. The central frequency of the visco-acoustic is 15 Hz due to attenuation. For the snapshots (Figure 2.14), the left half snapshots are acoustic waves, the right half ones are visco-acoutic waves. We notice the visco-acoustic waves have lower amplitude and faster velocity because the phase dispersion. We now introduce the kinematic changes and the kinematic velocity.

Kinematic velocity

Seismic attenuation leads to dispersion, leading to phase velocity changes. The kinematic information of visco-acoustic waves changes. We call the actual traveling speed of visco-acoustic waves the kinematic velocity. Here, we use the visco-acoustic wave equation to analyze the kinematic changes and introduce the kinematic velocity. First, the visco-acoustic wave equation (Equations 2.42) in frequency domain reads

iωp + κ(Lτ + 1)(∇ • v) + L l=1 r l p = f (x s , ω), (2.53) iω • v + 1 ρ ∂p ∂x = 0, (2.54)
(iω + 1 τ l σ )r l p + 1 τ l σ τ κ(∇ • v) = 0, l =∈ [1, L], (2.55)
where κ = ρv 2 p is called relaxed modulus, and v p is the P wave velocity. We substitute the term r l p from the third equation into the first equation, then taking divergence on the both sides of the second equation to replace the term ∇ • v. The second-order wave equation with constant density gives as

(iω) 2 p - 1 ρ κ(Lτ + 1) - L l=1 κτ 1 + iωτ l σ ∇ 2 p = F (x s , ω), l =∈ [1, L]. (2.56)
The new equation is equal to a second-order wave equation with a complex velocity

(iω) 2 p -v 2 c ∇ 2 p = F (x s , ω), ( 2
v 2 c (v p , τ, ω) = v 2 p (Lτ + 1) - L l=1 τ 1 + iωτ l σ , = v 2 p (Lτ + 1) - L l=1 τ 1 + ω 2 τ 2 σl + i • v 2 p L l=1 τ ωτ σl 1 + ω 2 τ 2 σl .(2.58)
The complex velocity can be spitted into a real term plus a imaginary part

v c (v p , τ, ω) = v k (v p , τ, ω) + i • v l (v p , τ, ω), (2.59) with v k = [v 2 c ] 2 + [v 2 c ] 2 + [v 2 c ] 2 , (2.60) v l = [v 2 c ] 2 + [v 2 c ] 2 -[v 2 c ] 2 . (2.61)
where v k is corresponding to kinematic effects, controlling the travel time. We call it the equivalent kinematic velocity. v l is related to the energy loss (Appendix 2.8). The v k is the function of ω, v p and τ , whose explicit form is

v k = v p Lτ + 1 -L l=1 τ 1+ω 2 τ 2 σl 2 + L l=1 τ ωτ σl 1+ω 2 τ 2 σl 2 + Lτ + 1 -L l=1 τ 1+ω 2 τ 2 σl 2 .
(2. For the band-limited seismic waves, we choose one single reference frequency to determine the kinematic velocity for the propagation of seismic waves. Generally, v k is determined for the central frequency of the source. For example, we simulate three kinds of waves with a 28 Hz central frequency Ricker wavelet (Figure 2.16). We set the kinematic velocity at 28 Hz, which is 2250 m/s. The three waves at the near and far offset have different amplitudes in the time domain and spectrum in the frequency domain in Figure 2.16. The correct traveltime is obtained at the maximum amplitude. The different combinations of velocity and attenuation lead to the same kinematic velocity. This relation is called kinematic relationship, it will be applied to guide the inversion to retrieve the velocity and Q models (Chapter 4).

Summary

This chapter first reviews visco-elasticity represented by rheological models and introduces the concept of complex modulus and stress relaxation function of visco-elastic media. To be consistent with field observations, different kinds of Q models including constant Q and nearly constant Q models are introduced. Constant Q is deduced by assuming a relaxation function in the shape of power function. Nearly constant Q models represented by rheological models are intuitive, e.g., GMB and GZB. A unification formalism of Q models represented by rheological description will be introduced in the next chapter. Nearly constant Q is widely applied in exploration seismology due to the simplicity of the implementation. The visco-elastic / visco-acoustic wave equation in the time domain is obtained based on the method of nearly constant Q optimization. The τ method reduces the number of attenuation parameters to a single term τ . We propose a modified τ method to avoid an empirical initial guess for nearly constant Q optimization. Based on 2D visco-acoustic wave equation, we have implemented the forward simulation based FDM and analyzed the accuracy by comparing with an analytical solution. Finally, we have analyzed the visco-elastic effects including dispersion and dissipation. The former leads to phase shift, and the latter causes the amplitude loss. In frequency domain, the visco-acoustic wave equation can be represented by a complex velocity. The real and imaginary term of complex velocity corresponding to kinematic effect and amplitude loss respectively. For a fixed frequency, visco-acoutic waves with various combinations of v p and Q (or τ ) have the same kinematic velocity except that the amplitudes are different. The information is important for building a relation between v p and Q (or τ ) to guide the inversion. This aspect will be introduced in Chapter 4.

Appendix: Analysis of the wave propagation

We analyze the plane wave propagation in a homogeneous visco-acoustic media. A plan wave at distance x and time t is represented by

u(x, t) = u 0 exp iω t - x v , (2.63)
where u(x, t) denotes the plane wave propagation. u 0 is the initial amplitude. ω is the angular frequency. v is the velocity of the waves. This is a constant real value for a homogenous elastic media. For a homogenous attenuated media, the amplitude of visco-acoustic waves decays along with the propagation distance and velocity disperses with frequency, caused by attenuation. The wave propagation is written as

u(x, t) = u 0 exp [-α(ω)x] exp iω t - x v(ω) , (2.64)
where v(ω) is real and changes with frequency. The phase of waves is determined by iω(t -x/v(ω)). α(ω) is the attenuation coefficient, equal to ω/2v(ω)Q. The wavefield with attenuation is

u(x, t) = u 0 exp iω t - x v(ω) 1 - i 2Q , (2.65)
with a complex velocity [START_REF] Aki | Quantitative seismology[END_REF])

1 v c (ω) = 1 v(ω) 1 - i 2Q . (2.66)
With the first order approximation [START_REF] Aki | Quantitative seismology[END_REF][START_REF] Toksoz | Seismic Wave Attenuation[END_REF][START_REF] Tarantola | Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation[END_REF][START_REF] Malinowski | High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF] with Q 1, it is equivalent to

v c (ω) = v(ω) 1 + i 2Q .
(2.67)

The complex velocity has a real and imaginary term, related to kinematic effects and energy loss. Recall the complex velocity (Equation 2.59)

v c (v p , τ, ω) = v k (v p , τ, ω) + i • v l (v p , τ, ω),
and the definition of Q (Equation 2.8)

1 Q(ω) = [M (ω)] [M (ω)] ,
Q is rewritten as

Q(ω) = [v 2 c ] [v 2 c ] , (2.68) = v 2 k -v 2 l 2v k v l , (2.69) = 1 2 v k v l - v l v k . (2.70)
Reorganize the third term, the ratio of v k and v l is

v k v l = 2Q + v l v k > 2Q.
(2.71)

We can neglect the small value term v l /v k due to

v l v k < 1 2Q 1.
(2.72) leading to the Q defined as [START_REF] Song | Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part II-Inversion method, synthetic experiments and real-data results[END_REF][START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF])

Q = v k 2v l . (2.73) Finally, we have v l = v k 2Q , (2.74)
which means that Equation 2.59 is consistent with Equation 2.67. v k is related to the kinematic effects, controls the phase of waves, and v l is related to the energy loss.

Chapter 3

Unification of visco-elastic wave equations

Abstract

Visco-elasticity is the essential ingredient for quantitative seismic imaging and geological interpretation in a number of contexts, such as in the presence of gas clouds. Decades of developments of numerical simulation of visco-elastic wave equations in seismology are mainly based on constant Q model, leading to numerous different forms of time-domain visco-elastic wave equations. Based on rheological models, [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] adopted the Generalized Maxwell body (GMB) to implement viscoelastic wave equations in time domain. Carcione et al. (1988a) incorporated the Generalized Zener body (GZB) into the time-domain visco-elastic wave equation. [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF] proved that visco-elastic complex modulus based on GMB and GZB are equivalent. However, from the rheological point of view, this formalism can not incorporate the fractional visco-elastic wave equations based on the constant Q model [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF]. [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity[END_REF] first mentioned that the constant Q model is based on a fractional Scott-Blair model. The stress-strain relationship of the Scott-Blair model is between a spring and a dashpot. Therefore, we review the various visco-elastic wave equations in the text of seismology. Based on the stress-strain constitutive law of rheological models, we propose a unification way to describe the existed visco-elastic wave equations. The unification formalism indicates that each kind of visco-elastic wave equation is composed by the combination of basic rheological elements, e.g., GMB, GZB.

In this chapter, we gather knowledge usually available in separated papers in the fields of fractional calculus, rheology, mechanics, and seismology. By unification formalism, we can establish more clearly links between different approaches used in seismology.

Introduction

Linear visco-elasticity means that the stress is represented by a time relaxation function convoluted with strain. There are mainly 3 ways to describe the visco-elasticity. First, from the rheological point of view, the stress-strain relationship for a linear viscoelastic system can be represented by a spring-dashpot mechanical model (Figure 3.0a, b). Different combinations of springs and dashpots compose the family of visco-elasticity models. The general form of constitutive relation is expressed by an ordinary differential equation [START_REF] Fung | Foundations of solid mechanics[END_REF]) (a): spring:

N n=0 a n d n σ(t) dt n = M m=0 b m d m (t) dt m , ( 3
σ(t) = M 0 (t) (b): dashpot: σ(t) = η ∂ (t) ∂t (c): spring-pot: σ(t) = G ∂ α (t) ∂ α t (d): Maxwell: σ(t) + η M 0 ∂σ(t) ∂t = η ∂ (t)
∂t (e): Kelvin:

σ(t) = M 0 (t)+η ∂ (t) ∂t (f): Burgers: σ(t) + ( η 1 M 1 + η 1 M 2 + η 2 M 2 ) ∂σ(t) ∂t + η 1 η 2 M 1 M 2 ∂ 2 σ(t) ∂ 2 t = η 1 ∂ (t) ∂t + η 1 η 2 M 2 ∂ 2 (t) ∂ 2 t
(g): Zener: However, there is another kind of equation describing the relation between stress and strain. The idea behind is that the material properties are determined by various states between an elastic solid and a viscous fluid characterized as spring-pot (Figure : 3.0c), rather than a combination of an elastic and a viscous element represented by previous ordinary differential formula [START_REF] Gemant | On fractional differentials[END_REF][START_REF] Blair | Limitation of the Newtonian time scale in relation to non-equilibrium rheological states and a theory of quasiproperties[END_REF]. It can be represented by a fractional-order differential form [START_REF] Park | Analytical modelling of viscoelastic dampers for structural and vibration control[END_REF])

σ(t) + η M 1 +M 2 ∂σ(t) ∂t = M 1 M 2 M 1 +M 2 (t) + M 1 η M 1 +M 2 ∂ (t) ∂t (h): GMB: σ(t) + η l M l ∂σ(t) ∂t = M 0 (t) + ( M 0 η l M l + η l ) ∂ (t) ∂t (i): GZB: σ(ω) = L l=1 M 1l M 2l +iωη l M 1l +M 2l +iωη l
N n=0 a n d pn σ(t) dt pn = M m=0 b m d qm (t) dt qm , 0 ≤ p n , q m ≤ 1 (3.2)
where p n and q m are fractional derivative. In analogy with ordinary differential models, [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity[END_REF] and [START_REF] Mainardi | Creep, relaxation and viscosity properties for basic fractional models in rheology[END_REF] displayed the fractional rheological models by replacing integer differential coefficients with fractional ones, e.g., fractional Maxwell model and fractional Zener body. Fractional differential models are widely applied in fields of materials and microstructure. Finally, another widely used theory to explain the visco-elasticity behaviors is called power law based on phenomenological models. The relaxation or creep time function has the shape of power function, which can well fit the observations of relaxation or creep behaviors of visco-elastic media [START_REF] Nutting | A new general law of deformation[END_REF]. The classic stress relaxation function gives as

ψ(t) = M 0 t -α , 0 ≤ α ≤ 1, (3.3)
where, ψ(t) is the stress relaxation function, M 0 is the elastic modulus, and α is the power of the function. The stress is the time convolution between the relaxation function with the time derivative of strain in the time domain (Equations 2.2, 2.5), leading to

σ(t) = M 0 Γ(1 -α) d α (t) dt α , (3.4)
where Γ(1 -α) is the gamma function, α is the time differential coefficient. The first introduction in seismology of constant Q model by [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] is based on the phenomenological power law. The power law with fractional coefficient actually belongs to fractional derivative models [START_REF] Bagley | Power law and fractional calculus model of viscoelasticity[END_REF][START_REF] Park | Analytical modelling of viscoelastic dampers for structural and vibration control[END_REF], which is also applied for the visco-elastic simulation and inversion in seismology [START_REF] Hanyga | Power-law attenuation in acoustic and isotropic anelastic media[END_REF][START_REF] Ribodetti | Some effects of the memory kernel singularity on wave propagation and inversion in viscoelastic media -II. Inversion[END_REF].

After analyzing the three formulas, we can write Equations 3.1, 3.2 and 3.4 as a general form

N n=0 a l n d pn σ(t) dt pn = M m=0 b l m d qm (t) dt qm , 1 ≤ l ≤ L (3.5)
where l is associated to a rheological model among total number L, p n = n and q m = m correspond to ordinary order representing, e.g., GZB, 0 ≤ p n , q m ≤ 1 corresponds to fractional order representing fractional derivative models, and models with N = 0, M = 0 can be explained by the classic power law.

To clarify the relations of different theories and the connections between different visco-elastic wave equations, we review different rheological models and applications in seismology. Most commonly used vicso-elastic wave equations obtained by the combinations of springs and dashpots, e.g. GZB and GMB models, are based on nearly constant Q model. [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] constant Q model and associated visco-elastic wave equation are based on power law. Other applications in seismology are based on power law [START_REF] Hanyga | Power-law attenuation in acoustic and isotropic anelastic media[END_REF][START_REF] Ribodetti | Some effects of the memory kernel singularity on wave propagation and inversion in viscoelastic media -II. Inversion[END_REF][START_REF] Wang | Generalized viscoelastic wave equation[END_REF]. The power law belongs to the family of fractional derivative models.

Two prominent papers should be mentioned. The first one is written by [START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF], who gave different mechanical visco-elastic models and relevant wave equations. Those mechanical models can be obtained from ordinary differential equation (Equation 3.1). The second is written by [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity[END_REF], who mentioned that [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] constant Q is a specific example of fractional differential models called Scott-Blair model. Those different visco-elastic wave equations and their rheological models need a unification formalism to be represented in exploration seismology. Thus we combine the three equations (Equation 3.1, 3.2 3.4) into a unification formula (Equation 3.5). This unification formula is a easy way to identify each specific form of viscoelastic wave equations and the related rheological models (Table 3.1). Furthermore, it is a way to obtain other visco-elastic wave equation to describe more complicated attenuation behaviors, which is not limited to seismology. We begin from the unified stress-

Stress-strain relationship

Unification formalism (Here) ⇓

Specific example

One spring with L pairs of M = N = 1, p 0 = 0, q 0 = 0 and p 1 = 1, q 1 = 1 M = N = 0, p 0 = 0, q 0 = α and L = 1

Others, e.g., fractional Zener M = N = 1, p 0 = 0, q 0 = 0 and p

1 = α, q 1 = β L = 1 Mainardi (2010) Q model Nearly constant Q Constant Q Others Rheological elements

Complex modulus and visco-elastic wave equations

For L pairs of rheological models in parallel, the stress and strain relation in frequency domain is given by

N n=0 a l n (iω) pn σ(ω) = M m=0 b l m (iω) qm (ω), 1 ≤ l ≤ L. (3.6)
Then the complex modulus [START_REF] Bland | The Theory of Linear Viscoelasticity[END_REF][START_REF] White | Seismic waves: Radiation, transmission and attenuation[END_REF] reads

M (ω) = L l=1 M m=0 b l m (iω) qm N n=0 a l n (iω) pn . (3.7)
The imaginary part corresponds to energy loss. The Q factor [START_REF] Boit | Theory of stress-strain relations in anisotropic visco-elasticity and relaxation phenomena[END_REF][START_REF] Knopoff | The seismic pulse in materials possessing solid friction[END_REF][START_REF] Fung | Foundations of solid mechanics[END_REF][START_REF] O'connell | Measures of dissipation in viscoelastic media[END_REF][START_REF] Ben-Menahem | Seismic waves and sources[END_REF] is defined as

1 Q = [M (ω)] [M (ω)] . (3.8)
From the constant Q or nearly constant Q model based on the observations of seismic waves and rock experiments [START_REF] Mcdonal | Attenuation of shear and compressional waves in Pierre shale[END_REF], one first needs to have the proper number for M , N , p n and q m to dertermine the basic rheological element. Then the objective is to determine the a n and b m corresponding to elastic modulus and viscosity coefficients and L to obtain constant Q within the frequency range [START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF].

With the constant Q optimization, then the general form of visco-ealstic wave equation is obtained with the complex modulus in the frequency domain. The general form of visco-acoustic wave equation with constant density is

(iω) 2 p - 1 ρ M (ω)∇ 2 p = 0, (3.9)
where p is the pressure and ρ is the density. The complex modulus M (ω) leads to visco-elastic effects including dispersion and dissipation. Different forms of M (ω) are obtained from different ways, e.g., phenomenological observation of relaxation or creep function, leading to the existed various forms of visco-elastic wave equations. From a rheological point of view, M (ω) can be derived from various rheological models. Thus we will derive visco-acoustic wave equations from the unification formalism including nearly constant Q based on GMB and [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] constant Q. Compared with the visco-acoustic wave equation, visco-elastic wave equation only requires to introduce Q p for P wave and Q s for S wave separately. Both Q p and Q s are determined by the same way. Thus we only illustrate the visco-acoustic wave equations from the unification formalism based on different rheological models, instead of visco-elastic wave equations.

Example I: Nearly constant Q and associated wave equations

Nearly constant Q model is approximated by several rheological models like GZB and GMB. [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] first introduced the GMB rheological model to obtain nearly constant Q. Carcione et al. (1988a,b) incorporated the GZB into the timedomain visco-elastic wave equation. [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF] proved the equivalence between GZB and GMB and put forward a material-independent formula. However, this material-independent formula is still based on the ordinary differential constitutive relation (Equation 3.1). We choose the integer values to obtain nearly constant Q from equation (Equation 3.5) and visco-acoustic wave equations.

The GMB is composed by a spring and L pairs of Maxwell body in parallel. A spring is represented as M = N = 0, p 0 = 0, q 0 = 0: the constitutive relation is

σ(ω) = M 0 (ω).
(3.10)

A Maxwell body is represented by M = N = 1, p 0 = 0, q 0 = 0 and p 1 = 1, q 1 = 1: the constitutive relation leads to

σ(t) + η l M l σ (1) (t) = M l η l M l (1) (t). (3.11)
In frequency domain is expressed as

σ(ω) + iω ω l σ(ω) = M l iω ω l (ω).
(3.12) ω l is equal to 1/t l , in which t l is referred to the relaxation time from the stress relaxation function of Maxwell body ψ(t) = M l exp(-t/t l )H(t). t l is defined as η l /M l . Complex modulus for l Maxwell body is

M l (ω) = M l iω ω l + iω . (3.13)
Therefore the total complex modulus of GMB is

M (ω) = M 0 + L l=1 M l iω ω l + iω , = M 0 + δM L l=1 α l iω ω l + iω , (3.14)
with α l δM = M l . M l is the elastic modulus of each Maxwell body and M 0 is the elastic modulus of the spring. Based on rheological models, one needs to optimize a n and b m corresponding to elastic modulus and viscosity coefficients and L to obtain constant Q within the frequency range (Equation 3.7, 3.8). Here a n and b m have been parameterized as ω l and α l . ω l is usually determined logarithmically over the frequency range [START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF]. Different α l optimizes different Q values.

The wave-equation reads

(iω) 2 p - 1 ρ (M 0 + δM L l=1 α l iω ω l + iω )∇ 2 p = 0 (3.15)
Therefore the form in time domain is

∂p ∂t + 1 ρ (M 0 + δM L l=1 α l )(∇ • v) + L l=1 r l p = f (x s , t), ∂v ∂t + 1 ρ ∇p = 0, ∂r l p ∂t + ω l (r l p + δM α l ∇ • v) = 0, l =∈ [1, L]. (3.16)
v is the particle velocity vector. r l p is the memory variables. From the general formula of stress-strain relationship, we get the visco-acoustic wave-equation based on GMB rheological model. The visco-equation can be obtained by GZB rheological model. Next section, we introduce the [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] constant Q visco-acoustic wave equation the under frame of the unification formalism.

3.4

Example II: Kjartansson constant Q and associated wave equations [START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF] gives a relaxation function and complex modulus defined as

M (ω) = M 0 ( iω ω 0 ) 2γ , 0 < 2γ < 1 (3.17)
where ω 0 is the reference frequency. The complex modulus is derived from a power creep function. [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity[END_REF] pointed out this form can be obtained from the rheological model called Scott-Blair model. In fact, the Scott-Blair model called by [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity[END_REF] is a spring-pot model (Figure 3.0). The stress of a spring-pot is linear, with the fractional derivative of strain. Taking M = N = 0, p 0 = 0, q 0 = α and L = 1 in Equation 3.5, the constitutive relation turns into

σ(t) = M 0 ( η 0 M 0 ) α (α) (t) (3.18)
by choosing the b 0 /a 0 as M 0 (η 0 /M 0 ) α . M 0 stands for the elastic modulus, η is the viscous coefficient of a dashpot. Defining t 0 = 1/ω 0 as η 0 /M 0 , t 0 is referred as the relaxation time. The stress-strain relation in the frequency domain is

σ(ω) = M 0 ( 1 ω 0 ) α (iω) α (ω). (3.19)
Thus the complex modulus is defined as

M (ω) = M 0 ( iω ω 0 ) α = M 0 ( iω ω 0 ) 2γ , (3.20)
with α = 2γ, we have the complex modulus equivalent to Equation 3.17 which is obtained through the power law. The parameter α depends the Q value. The wave equation is derived by introducing the complex modulus as

(iω) 2 p - 1 ρ M 0 ( iω ω 0 ) 2γ ∇ 2 p = 0. (3.21)
Combing the (iω) 2γ with (iω) 2 and transferring into time domain, the visco-acoustic wave equation reads

∂ 2-2γ p ∂t 2-2γ - M 0 ρ ω -2γ 0 ∇ 2 p = 0. (3.22)
This function is consistent with the result of [START_REF] Carcione | Time-domain modelling of constant-Q seismic waves using fractional derivatives[END_REF]. This form of viscoacoustic wave equation is obtained with a spring-pot rheological model. The constant Q visco-acoustic wave equation in time domain involves with time-domain fractional derivatives. It is hard to be solved by finite difference method. Thus the commonly applied visco-acoustic wave equation is based on GMB / GZB rheological models by introducing memory variables instead of solving fractional time derivative [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]Carcione et al., 1988a).

Conclusions

The general constitutive relation of stress and strain combines ordinary differential system, fractional order differential system and power law together. From the constitutive relation, the unification formalism is a way to describe and summarize visco-elastic wave equations in seismology. Depending on the different strategies of optimizing constant Q, the rheological models of visco-elastic wave equations contain different combinations of basic rheological elements with springs, dashpots and spring-pots. The widely applied visco-elastic wave equations in seismology can be characterized as different specific cases. Furthermore, it is a way to obtain other kinds of visco-elastic wave equations, e.g. based on fractional Zener model. Those visco-elastic wave equations can be summarized under the unification formalism.

Résumé du chapitre 4

Abstract

Seismic attenuation is a very useful property to detect fluid and gas clouds combined with the conventional velocity-based seismic imaging in exploration geophysics. To obtain attenuation and velocity models, visco-acoustic full waveform inversion (FWI) is an advanced and accurate method compared to traditional ray-based methods. Viscoacoustic FWI has more realistic rheology characterizing the Earth media and exploits all the wave information to retrieve attenuation and velocity. However, attenuation and velocity inversion suffers cross-talks between parameters and different parameter sensitivities. Many researches focus on different inversion strategies to mitigate cross-talks.

Those strategies mainly include sequential inversion, simultaneous inversion and hierarchical inversion. Because of dispersion, different combinations of attenuation and velocity may have the same kinematic information for band-limited data, called kinematic relationship. The real travel speed of visco-acoustic waves is called kinematic velocity. We propose a hybrid strategy to explicitly use the kinematic relationship during full waveform inversion to retrieve the parameters. In a first step, a classical visco-acoustic inversion is applied for a fixed quality factor to obtain kinematic velocity. Then a subsequent FWI with a modified gradient of the objective function is put forward by introducing the relation between attenuation and velocity based on the kinematic velocity. In case of inaccurate kinematic velocity derived in the first step, we perform alternate steps between the determination of kinematic model and visco-acoustic models. The hybrid FWI leads to a faster convergence compared to the conventional inversion strategies. This is first validated on simple 2D synthetic data sets mainly dominated by transmitted waves. Finally, the hybrid strategy is applied to a synthetic Middle-East model with surface acquisition geometry; the inversion has a fast convergence and lower final misfit.

The retrieved attenuation has closer exact average values compared with conventional inversion strategies. Both simple and complex models validate that the hybrid inversion strategy can mitigate cross-talks between two parameters and retrieve attenuation and velocity more effectively.

Introduction

Conventional seismic inversion relies on an acoustic Earth model only depending on velocity. As the demands of high-resolution inversion, more accurate rheology models have been taken into consideration, such as elasticity, anisotropy, and attenuation.

For fluid and gas clouds detection, attenuation is more sensitive than velocity of pure acoustic assumption [START_REF] Toksoz | Seismic Wave Attenuation[END_REF][START_REF] Best | The relationships between the velocities, attenuations and petrophysical properties of reservoir sedimentary rocks[END_REF][START_REF] Wang | Seismic Inverse Q Filtering[END_REF]. The property of attenuation is very useful to compensate for energy loss and subsequently to obtain accurate seismic images, such as removing seismic image blurring beneath the gas clouds [START_REF] Zhou | Compensating attenuation due to shallow gas through Q-tomography and Q-PSDM, a case study from Brazil[END_REF]. Seismic attenuation is divided into two categories, one is the scattering attenuation [START_REF] Marion | From ray to effective medium theories in stratified media: an experimental study[END_REF]) and the other one is the intrinsic attenuation [START_REF] Aki | Quantitative seismology[END_REF][START_REF] Klimentos | Relationships among compressional wave attenuation, porosity, clay content, and permeability in sandstones[END_REF][START_REF] Thomsen | Converted wave imaging of Valhall reservoir[END_REF]. The scattering attenuation can happen in any kind of media with multi scatters.

The intrinsic attenuation is related to the property of media. Waves transfers the energy to heat due to the internal friction. The seismic attenuation usually refers to intrinsic attenuation, which can be represented by quality factor Q [START_REF] Knopoff | The seismic pulse in materials possessing solid friction[END_REF]. The reciprocal of Q stands for the wave energy loss per traveling cycle.

Traditional methods for estimating Q include ray-based methods, based on the spectral amplitude loss or spectrum distortion due to attenuation. Spectral amplitude loss method [START_REF] Brzostowski | 3D tomographic imaging of near-surface seismic velocity and attenuation[END_REF][START_REF] Leggett | Study of traveltime and amplitude time-lapse tomography using physical model data[END_REF][START_REF] Cavalca | Estimation of the near-surface Q distribution from pre-stack surface seismic data[END_REF][START_REF] Cavalca | Raybased Tomography for Q Estimation and Q Compensation in Complex Media[END_REF] retrieves the velocity and Q model using traveltime tomography. Q is estimated by integrating the amplitude loss along the ray path [START_REF] Futterman | Dispersive body wave[END_REF][START_REF] Aki | Quantitative seismology[END_REF][START_REF] Bregman | Traveltime and amplitude inversion in crosshole seismology[END_REF][START_REF] Tonn | The determination of the seismic quality factor Q from VSP data: A comparison of different computational methods[END_REF]. However, amplitude of seismic waves is easily contaminated by many factors such as scattering, geometric spreading, source and receiver coupling, radiation patterns, and transmission/reflection effects. It is hard to eliminate those effects. To overcome this problem, spectrum distortion methods include pulse broadening method [START_REF] Zucca | Seismic velocity and attenuation structure of the Geysers geothermal field[END_REF] and centroid frequency shift method [START_REF] Quan | Seismic attenuation tomography using the frequency shift method[END_REF]Plessix, 2006a). Those methods compare the spectrum difference between incident and recorded data to estimate the Q model. This kind of methods is based on the assumption that the spectrum of incident waves is known, and the observed spectrum also has the similar shape.

Beyond the traditional methods using ray information, full waveform inversion (FWI) [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations: Conference on Inverse Scattering, Theory and Application[END_REF]Tarantola, 1984a) is a powerful method to retrieve physical models of subsurface, including velocity, density, attenuation and anisotropy. FWI is effectively applied to non-linear multi-parameter inversion by exploiting full wave information [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]. [START_REF] Tarantola | Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation[END_REF] first derived the full waveform inversion algorithm to retrieve the elastic and attenuation parameters. Visco-elastic or viscoacoustic FWI can be implemented in frequency domain [START_REF] Song | Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part II-Inversion method, synthetic experiments and real-data results[END_REF][START_REF] Causse | Finite-difference inversion in viscoacoustic media[END_REF][START_REF] Liao | Multifrequency viscoacoustic modeling and inversion[END_REF][START_REF] Causse | Preconditioning of full-waveform inversion in viscoacoustic media[END_REF][START_REF] Hicks | Reflection waveform inversion using local descent methods: estimating attenuation and velocity over a gas-sand deposit[END_REF][START_REF] Pratt | Waveform tomography images of velocity and inelastic attenuation from the Mallik 2002 crosshole seismic surveys[END_REF][START_REF] Kamei | Waveform tomography strategies for imaging attenuation structure for cross-hole data[END_REF][START_REF] Malinowski | High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full waveform inversion[END_REF][START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF][START_REF] Operto | Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation[END_REF][START_REF] Gao | Simultaneous inversion for velocity and attenuation by waveform tomography[END_REF][START_REF] Operto | On the role of density and attenuation in threedimensional multiparameter viscoacoustic VTI frequency-domain FWI: an OBC case study from the North Sea[END_REF] and time domain (Liao and McMechan, 1995;[START_REF] Charara | The state of affairs in inversion of seismic data: an OVSP example[END_REF][START_REF] Charara | Full waveform inversion of seismic data for a viscoelastic medium[END_REF][START_REF] Barnes | Borehole seismic data inversion for attenuating media[END_REF][START_REF] Charara | Crosswell seismic data tomography for heterogeneous viscoelastic media[END_REF]Barnes andCharara, 2009, 2010;[START_REF] Bai | Waveform inversion with attenuation[END_REF][START_REF] Kurzmann | Acoustic full waveform tomography in the presence of attenuation: A sensitivity analysis[END_REF][START_REF] Bai | Viscoacoustic waveform inversion of velocity structures in the time domain[END_REF][START_REF] Barnes | [END_REF][START_REF] Kurzmann | Viscoacoustic full waveform inversion for spatially correlated and uncorrelated problems in refection seismics[END_REF][START_REF] Belahi | Elastic versus viscoelastic full waveform inversion of near-offset and wide-angle data in the presence of attenuation[END_REF][START_REF] Cheng | A new approach of viscoacoustic waveform inversion in the time domain[END_REF][START_REF] Li | 2D time-domain visco-acoustic full waveform inversion method[END_REF][START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF][START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF][START_REF] Jiang | A modified 2D visco-acoustic waveform inversion scheme[END_REF][START_REF] Wang | Visco-acoustic fullwaveform inversion in the presence of complex gas clouds[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. Frequency-domain FWI is implemented by incorporating a complex modulus or complex velocity into elastic wave equation. Time-domain FWI needs to introduce extra memory variables instead of solving the time convolution. For large-scale 3D inversion, the preferred choice is applied in the time-domain [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. For computational reasons, the applications rely on the visco-acoustic assumption and not on the visco-elasticity [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF][START_REF] Wang | Visco-acoustic fullwaveform inversion in the presence of complex gas clouds[END_REF].

Visco-acoustic FWI with two parameters still encounters many difficulties. First, due to attenuation, energy of waves in visco-elastic media is weaker than pure elastic waves. With a simple but fairly correct Q model, mono-parameter velocity FWI regards Q as a passive parameter. Dissipation happens twice during source propagation and back-propagation of adjoint source. The weak energy leads to the difficulty of retrieving velocity. [START_REF] Kurzmann | Acoustic full waveform tomography in the presence of attenuation: A sensitivity analysis[END_REF] applied a smooth Q as a passive parameter to retrieve v p using visco-acoustic FWI. The shallow part is well retrieved, however resolution of the deep part decreases as the attenuation increases with traveling distance of wave.

Second, the limited observation leads to the problem of non-uniqueness [START_REF] Barnes | [END_REF][START_REF] Kurzmann | Viscoacoustic full waveform inversion for spatially correlated and uncorrelated problems in refection seismics[END_REF]. Different solutions can represent the data well. As the two parameters v p and Q are characterized as a complex velocity in the frequency domain, v p and Q have the same radiation patterns [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. The gradients of the objective function with respect to v p and Q are hard to distinguish, especially with limited acquisition [START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. It leads to the low resolution of parameters, which can be explained as cross-talks [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF].

Third, parameter gradients have different sensitivities [START_REF] Malinowski | High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full waveform inversion[END_REF]. Different inversion strategies have to be considered, which include sequential inversion, simultaneous inversion and hierarchical inversion. Sequential inversion consists of first determining the velocity parameter for a fixed Q factor, and in a second phase of inverting the Q parameter as the velocity model is better reconstructed. [START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF] suggested this inversion strategy for a cross-hole model to decrease the footprints of simultaneous inversion. However, this sequential approach may lead to data over-fitting [START_REF] Kurzmann | Acoustic full waveform tomography in the presence of attenuation: A sensitivity analysis[END_REF][START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF]. To avoid this, simultaneous inversion updates v p and Q at the same time with gradient optimizations. The gradients of v p and Q have different magnitudes without the introduction of the Hessian matrix (second-order derivative of the objective function). However, the inverse of Hessian to normalize the gradients of different parameters is hard due to the ill-posed Hessian [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF]. Therefore, proper magnitude of damping factors have to be added to the Hessian matrix to scale the parameter gradients and reconstruct the models. The hierarchical inversion [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF] combines a mono-parameter velocity inversion with simultaneous inversion to focus on the dominant information v p . A good estimation of v p contributes to the prior information for latter simultaneous inversion [START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF].

Besides those, prior information, smoothing operator and regularization have to be considered [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. [START_REF] Pratt | Waveform tomography images of velocity and inelastic attenuation from the Mallik 2002 crosshole seismic surveys[END_REF] and [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF] suggested a stronger Q smoothing factor than velocity smoothing factor to retrieve a better results in a highly heterogeneous attenuated media. A strong soothing can effectively overcome the oscillatory artefacts on the attenuation model. Regularization is applied to handle the ill-posed problem, but the choice of damping factor is still challenging [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF].

Despite recent progresses, visco-acoustic FWI still needs to find a way to decouple the two parameters. The previous inversion strategies have more flexible degrees of freedom to update each parameter. On the one side, the inversion could reach to a global minimum with the correct updating direction. On the other side, the inversion may be stuck into a local minimum due to cross-talks between parameters. We propose to introduce a constrained relation to guide the inversion.

The dispersion of visco-acoustic waves leads to phase velocity modifications. Because of the band-limited aspect of seismic waves, different combinations of velocity and Q parameters can have the same equivalent kinematic velocity. v p is the dominant information compared with Q, FWI preferentially retrieves the velocity information due to the impacts on kinematics [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF]. The mono-parameter velocity inversion for visco-acoustic data tends to update the kinematic velocity. In this work, we propose a hybrid FWI algorithm to take into account the kinematic information to retrieve the two parameters. The objective is to handle cross-talks between velocity and attenuation, without relying on the information of the Hessian matrix. The paper is organized as follows: we first indicate how we model attenuation, then we present the modified FWI approach. It mainly leads to a modified gradient of the objective function. We compare the behaviors of the classical and new hybrid strategies on 2D inversions on simple synthetic data sets, and draw some first conclusions on the recommended strategy for inverting velocity and attenuation. Finally, we apply the modified algorithm to a more complex 2D synthetic Middle-East model.

Visco-acoustic FWI

Attenuation effects lead to phase velocity dispersion [START_REF] Futterman | Dispersive body wave[END_REF][START_REF] Aki | Quantitative seismology[END_REF]. It means that different combinations of velocity and attenuation result in the same kinematic velocity v k for band-limited seismic data. The v k model controls the kinematic of wave propagation. The v k model, precisely defined in the next section, is different from the P-velocity model v p because of attenuation. Based on v k , we put forward a hybrid inversion strategy for visco-acoustic FWI.

Visco-acoustic wave equation and kinematic velocity

The time-domain visco-acoustic wave-equation based on nearly constant Q theory [START_REF] Liu | Velocity dispersion due to anelasticity; Implications for seismology and mantle composition[END_REF][START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]Carcione et al., 1988a,b;[START_REF] Robertsson | Viscoelastic finite-difference modeling[END_REF][START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF] involves a time convolution in the stress-strain relation. The constant Q theory can be obtained by L pairs of Standard Linear Solid (SLS) within desired frequency range, which is given as

1 Q(ω) = [M (ω)] [M (ω)] = L l=1 ω(τ l -τ σl ) 1+ω 2 τ 2 σl L l=1 1+ω 2 τ l τ σl 1+ω 2 τ 2 σl . (4.1)
The reciprocal of quality factor Q is defined as the ratio of imaginary and real part of complex modulus M (ω). τ σl and τ l are the stress and strain relaxation times for the lth SLS within total number L. To obtain constant Q within a frequency range, τ σl are given in advance and distributed logarithmically over the frequency range of interest, which yields good approximations to a constant Q with τ σl = 1 /ω l . τ l in total number L are optimized to have constant Q. To simplify the implementation, the τ method [START_REF] Blanch | Modelling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique[END_REF] is applied to reduce the L parameters into one parameter τ = τ l /τ σl -1.

With the single τ parameter, the visco-acoustic wave-equation reads

∂p ∂t + κ(Lτ + 1)(∇ • v) + L l=1 r l p = f (x s , t), (4.2)
∂v ∂t + 1 ρ ∇p = 0, (4.3)
∂r l p ∂t + 1 τ l σ (r l p + τ κ(∇ • v)) = 0, l =∈ [1, L] (4.4)
where v is the velocity vector, p is acoustic pressure field, and κ = ρv 2 p is the relaxed modulus. The wavefields depend on the source position s, the space coordinate x and the time t, written as p(s, x, t) and v(s, x, t). Model parameters are related to model space x, wring as v p (x) and τ (x). To analyse the visco-acoustic kinematic velocity in details, we rewrite the constant density wave equation in the frequency domain as a second-order differential equation

(iω) 2 p -v 2 c ∆p = F (x s , ω), (4.5)
with a complex velocity

v 2 c = v 2 p (Lτ + 1) - L l=1 τ 1 + iωτ l σ .
(4.6)

The complex velocity can be split into a real term and an imaginary part

v c = v k (v p , τ, ω) + iv l (v p , τ, ω), (4.7)
where v k and v l are the real and imaginary parts of the complex velocity, associated to the kinematic part and the energy loss, respectively. We call v k the equivalent kinematic velocity. It stands for real travel speed of visco-acoustic waves. v k is the function of v p , τ and frequency ω. We call it the kinematic relationship, written as

v k = f (v p , τ, ω) = v p Lτ + 1 -L l=1 τ 1+ω 2 τ 2 σl 2 + L l=1 τ ωτ σl 1+ω 2 τ 2 σl 2 + Lτ + 1 -L l=1 τ 1+ω 2 τ 2

FWI Strategies

FWI exploits all the wave information to retrieve the model parameters. The basic principle is matching the waveform to reduce the data misfit between the calculated and observed data (Tarantola, 1984a;[START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]. The objective function is given as

J (m) = 1 2 d cal (s, r, t) -d obs (s, r, t) 2 .
(4.9) d cal (s, r, t) is the selected from the whole wavefields at the source s and receiver r. m is the model parameter. m is composed of v p (x) and τ (x) for visco-acoustic FWI.

Theoretically, m can be retrieved from a global optimization method through the whole model space, e.g. Monte-Carlo method [START_REF] Sambridge | MONTE CARLO METHODS IN GEO-PHYSICAL INVERSE PROBLEMS[END_REF]. However, it is almost impossible due to the large model space and computation limitations. The practical solution is to apply FWI based on gradient updating method, e.g. conjugate gradient method (details of FWI is in Appendix 4.6).

Conventional visco-asouctic full waveform inversion suffers from cross-talks between velocity and τ (or Q) parameters. To overcome cross-talks, the inversions strategies [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF] include: sequential inversion, simultaneous inversion and hierarchical inversion. The key point of the different inversion strategies based on the conjugate gradient method is to compute the gradients of the objective function in equation 4.25. The derivatives of J with respect to v p and τ can be obtained with the adjoint state method (details in Appendix 4.7). We review different inversion strategies (Strategies A and B) and introduce the new approach (Strategy C and the hybrid strategy).

Strategy A

Visco-acoustic FWI preferentially retrieves the velocity information due to the kinematic effects [START_REF] Kamei | Inversion strategies for visco-acoustic waveform inversion[END_REF][START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Fabien-Quellet | Time domain viscoelastic full waveform inversion[END_REF]. If τ is a passive parameter, the mono-parameter velocity inversion tends to retrieve the kinematics. The kinematics are the combinations of the retrieved v p and fixed τ (Equation 4.8). The gradient of the objective function is

∂J (v p , τ ) ∂v p τ , (4.10)
where | τ means the derivative is obtained for fixed τ . This strategy is called Strategy A to retrieve v p . If τ is zero, Strategy A is a purely acoustic FWI. In this case, the acoustic FWI is also used to examine the Q impacts on the data. For a fixed τ , the iteratively updating of v p leads to a modified v k according to Equation 4.8.

Strategy B

Simultaneous inversion updates the velocity and attenuation at the same time. It is applied to avoid parameter over-estimating [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF] where | τ and | vp mean the derivatives are obtained for fixed τ and v p , respectively. v p and τ are considered as independent parameters.

Strategy C -new approach

Different combinations of v p and τ , visco-acoustic waves have same the kinematic velocity at a given frequency in equation 4.8. When the kinematic velocity v k is fixed, v p and τ are linked by the kinematic relationship. A modified inversion strategy is applied to decouple the v k as the combinations of updated v p and τ to match the amplitudes. For fixed v k , τ is the function of v p . The gradient of v p turns into a full differential formula

dJ (v p , τ (v p )) dv p v k = ∂J (v p , τ (v p )) ∂v p v k + ∂τ (v p ) ∂v p v k ∂J (v p , τ (v p )) ∂τ v k , (4.13)
where

∂τ (v p ) ∂v p v k = - ∂f (v p , τ ) ∂v p ∂f (v p , τ ) ∂τ v k . (4.14) | v k means the derivative is obtained for fixed v k . f (v p , τ
) is the function of kinematic velocity (Equation 4.8). The derivative of τ (v p ) with respect to v p for a fixed v k (Equation 4.14) can be obtained analytically (details in Appendix 4.8). The first step is to update v p based on Equation 4.13. Then τ is updated form the kinematic relationship (Equation 4.8). This strategy is called Strategy C, which guides the non-linear visco-acoustic FWI.

Strategy C needs a kinematic velocity model. The first step consists of updating the kinematic velocity v k and then further decouple the v p and τ parameters.

Hybrid strategy

When the kinematic velocity is not correct enough, v k can be updated for fixed τ . The objection function is function of v p . Then the updated v k can be used for the kinematic velocity for Strategy C. The gradient is

dJ (v p (v k )) dv k τ = ∂v p (v k ) ∂v k τ ∂J (v p (v k )) ∂v p τ , (4.15)
where

∂v p (v k ) ∂v k τ = 1 ∂f (v p , τ ) ∂v p τ . ( 4

.16)

| τ means the derivative is obtained for fixed τ . The derivative of v p (v k ) with respect to v k is explained in Appendix 4.9.

This hybrid strategy includes a first step for retrieving v k , and a subsequent step of Strategy C to update v p and τ for a fixed v k . Actually, the first step is equivalent to the Strategy A (details in Appendix 4.9). In practice, we propose the two-step hybrid strategy: (1) For a fixed τ , we iteratively update v p (Strategy A). This leads to a modified v k model according to kinematic relationship (Equation 4.8). ( 2) In a second-step, we still update v p , but here for a fixed v k model. The gradient of the objective function is modified accordingly. During the second phase, τ is the passive parameter, i.e. is automatically updated when v p is modified (Strategy C).

In this paper we have 2 schemes of the hybrid inversion strategies include A → C, and A → C A (→: followed by, : alternate scheme). For the simple models, the A → C scheme is enough. For the complex models, the scheme of A → C

A is preferred. The first several iterations of A is to obtain an initial v p model and v k model for Strategy C, then alternate C

A is applied to retrieve v p and τ and continuously improve the kinematic velocity. Additionally, without updating the kinematic velocity v k by Strategy A, the hybrid strategy turns into pure Strategy C. We also test the Strategy C with known v k for simple models. Among the schemes of the hybrid inversion strategies, Strategy C is the key algorithm to guide the visco-acoustic FWI.

In summary (Table 4.1):

• Strategy A is visco-acoustic inversion, where only v p is updated, while τ (or Q) remains fixed (Equation 4.10);

• Strategy B is the classical visco-acoustic FWI where v p and τ (or Q) are simultaneously inverted with conjugate gradient (Equation 4.11 and 4.12);

• Strategy C is the modified viso-acoustic FWI based on fixed v k , v p is updated by Equation 4.13, then τ (or Q) is derived from Equation 4.8.

• Hybrid Strategy is composed of Strategies A and C alternately; Strategy A is used for updating v k , and Strategy C is applied to guide the inversion to retrieve v p and τ .

Numerical examples

We apply Strategy C and the hybrid strategy to investigate the feasibilities on different models, and compare with Strategy A and Strategy B (Table 4

τ n = τ 0 fixed v k , v k = f (v n p , τ n ) Hybrid strategy A → C, or A → C A A given initial v 0 p A given initial τ 0 v k is updated by Strategy A, v k = f (v A p , τ A )
retrieve the parameters. The velocity model is homogeneous, but the attenuation model has a Gaussian anomaly. The last category (example 4.3.4) is applied to 2D synthetic Middle-East model. The objective is to show the feasibility of Strategy C and the hybrid strategy on complex models. The key factors are: (1) how to obtain an accurate kinematic velocity for complex models, and (2) how to guide the inversion based on fixed v k . Finally, we discuss if the method is a way of decoupling the parameters by explicitly introducing the kinematic relationship. To analyze the features of different strategies, we project the mean values of results of each iterations to the kinematic curve, to obtain the the searching direction for the three strategies (Figure 4.5).

Homogeneous model -transmission acquisition

First, for the results of mono-parameter velocity inversion (Strategy A) and simultaneous inversion (Strategy B) after first several iterations, results are located on the kinematic curve. v p is close to the correct v k 2060 m/s, τ is close to zero. Thus the first several inversions focus on correcting the kinematic information. It indicates that the visco-acoustic FWI preferentially retrieves the v p to match the phase of the observed and calculated data.

Second, after several iterations of Strategy B, the inversion results show that v p and τ are on the kinematic curve (e.g., 2rd iteration). The next searching direction (e.g., 3rd iteration) do not follow the kinematic curve to decouple the v p and τ . Once the result is located on kinematic curve, it means that the phase is matched well (Figure 4.6a). The updating space for velocity is very limited, due to good match of phase, leading to small values of v p gradient under the area of acquisition coverage (Figure 4.8a). The τ gradient (Figure 4.8b) is obviously positive. In this iteration, the kinematic information is accurate enough, the velocity updating depends on τ updating. Then the updated velocity will distort the phase match (Figure 4.6b), that is the reason why the searching direction is off the kinematic curve. The next followed iteration (e.g., 4th iteration) updates v p to correct the kinematic information (Figure 4.6c). The v p gradient in Figure 4.8b is large, leading to large decreases to approach the true value. v p is the dominant information to update.

Third, the searching direction of Strategy B indicates velocity updating is not always enough. The searching direction of Strategy C strictly follows the kinematic curve to decrease the amplitude difference and keep the phase matching (Figure 4.7). Thus, the modified inversion scheme by building the relation between velocity and τ compensates and corrects the velocity updating, with obvious positive full derivative dJ /dv p (Figure 4.9). The positive dJ /dv p denotes the large decreases for velocity updating.

Fourth, without the estimation of the inverse of Hessian matrix, the scaling between the ∂J /∂v p and ∂J /∂τ has to be estimated and given in advance, under conjugate gradient scheme. During the inversion, the fixed scaling factor may constrain the parameter updating. For example, the retrieved results which are on the kinematic curve of Strategy B, the increase scaling factor on ∂J /∂v p will improve the velocity updating. Thus, the velocity derivative with respect to τ acts as a scaling factor between velocity and τ gradients.

In summary, we apply Strategy C for a homogeneous model, and compare with the mono-parameter inversion (Strategy A) and simultaneous inversion (Strategy B). With a given kinematic velocity, Strategy C is a guided inversion and keeps phase matching. In the end, it converges faster than simultaneous inversion. In next section, we still apply the three strategies, but for a model with blocky anomalies to further investigate the inversion strategies.

Blocky model -transmission acquisition

In this section, we further test the Strategy C for a slightly more complex model with blocky v p and τ anomalies, and compare with the other two strategies. Following the same inversion sequence as previous test 4.3.1, we first apply the Strategy A for mono-parameter inversion. A few iterations (3 iterations in Figure 4.11a,b) leads to close value between mean velocity of v p and true v k . This is consistent with the previous discussion that the first a few iterations of visco-acoustic FWI preferentially update the kinematic information. More iterations (30 iterations in Figure 4.11c, d) generate more artifacts, leading to the over-estimation of v p without introducing τ . It is necessary to update τ , as well as v p .

Strategy B updates the two parameters simultaneously. 3 iterations over v p with Strategy B (Figure 4.11e) has the similar result of Strategy A, meanwhile, the retrieved τ after 3 iterations (Figure 4.11f) is not obvious. This indicates that the first several iterations focus on the kinematic information. After more iterations with Strategy B (Figure 4.11g, h), the anomaly is gradually retrieved, but not very clear.

The Strategy C is also applied to reconstruct the v p and τ models. The initial model of Strategy C is from the v p result of Strategy A after 3 iterations (Figure 4.11a). The initial τ is zero, also same as Strategy A. The kinematic velocity (Figure 4.12) is obtained by smoothing v p result of Strategy A after 3 iterations. Results of Strategy C after 3 iterations (Figure 4.11i, j) are comparable with the results of Strategy B after 30 iterations. After more iterations with Strategy C, the v p and τ anomalies are more obvious (30 iterations in Figure 4.11k,l).

From the misfit curve ( We analyze the inversion results under the area of effective acquisition coverage, noticing the strong edge effects (Figure 4.11). The inverted values of v p and τ change gradually from the sources to the receivers, relying on the information derived from transmitted waves. According to the test and analysis of the previous section, Strategy B converges slowly because the inversion is always off the kinematic curve, but Strategy C is guided and follows the kinematic curve. It leads to that Strategy C has better and more obvious inversion results on v p and τ anomalies compared with Strategy B.

In this section, as well as the previous section, the kinematic velocity v k of the two tests is homogeneous and obtained by the results of Strategy A after 3 iterations. In the next section, we extend the different inversion strategies to a more complex models, with the updating of v k model.

Gaussian model -transmission acquisition

In the previous sections (Section 4.3.1, 4. Strategy B updates both v p and τ models. The inverted τ is only slightly updated after 3 iterations, but with a strong v p artifacts (Figure 4.15e, f). With more iterations, the inverted v p has stronger artifacts. The inverted τ anomaly (Figure 4.15h) increases, but is much lower than the true value. The misfit curve of Strategy B can not decrease any more (Figure 4.16 green).

To mitigate the v p artifacts, we apply the hybrid strategy for this model, defining as A3 → C A (A3 denotes the 3 iterations of Strategy A). After 3 iterations of Strategy A, we start with Strategy C for 2 iterations, then continue another 2 iterations with Strategy A. To improve the kinematic velocity v k , we perform each 2 iterations of Strategy A and Strategy C in an alternate way until the end. Following this sequence with 3 iterations, the τ anomaly (Figure 4.15i) is slightly updated. The v p model (Figure 4.15j) still has artifacts. With more iterations, the final inverted τ anomaly (30 iterations in Figure 4.15k) is more obvious than the one obtained with Strategy B, and the v p artifacts are much lower. The final misfit of the hybrid strategy is also lower than Strategy B (Figure 4.16 blue).

Due to cross-talks between parameters, the final results of simultaneous inversion (Strategy B) leave strong v p artifacts. The hybrid strategy guides the inversion and mitigates cross-talks. Although, the final results still have v p artifacts, but those artifacts are much lower than the one obtained with Strategy B. The τ anomaly is better retrieved. Without the introduction of the Hessian matrix, the parameter gradients of the objective function need a scaling factor to balance the sensitivities. The extreme situation is that the scaling factor of v p gradient is fixed at zero, leading to no contributions on v p up-dating. The inversion only focuses on retrieving τ anomaly. The term ∂τ (v p )/∂v p | v k of the hybrid strategy acts as the scaling factor to balance the sensitivities. The τ (v p ) derivative with respect to v p is based on the fixed v k . Accurate v k is important for Strategy C to mitigate cross-talks. We apply the Strategy C with known aproximated v k (Figure 4.17 In the next section, we apply the hybrid strategy for the 2D synthetic Middle-East model, to test the feasibility on more complex models.

Middle-East model -surface acquisition

The previous tests are based on simple models to investigate the features and feasibility of the hybrid strategy and Strategy C. The first several iterations with a mono-parameter inversion could basically retrieve the kinematic velocity. Based on a good approximate kinematic velocity, this modified strategy decouples the parameters and mitigates crosstalks, with faster convergence than the classical simultaneous inversion. However, the question is how to retrieve the kinematic model in complicated models. In this section, we apply the hybrid strategy to iteratively update v k and decouple v p and τ in a 2D synthetic Middle-East model. As for the previous sections, we also compare the hybrid strategy with Strategy A, Strategy B and an additional sequential inversion (Figure 1.9). The model is 800 m laterally and 212 m in depth. 20 m water layer is on the top of the model. 81 sources are located on the surface from 80 m to 720 m, with 361 receivers evenly distributed on the surface from 40 m to 760 m. The source is a Ricker wavelet with central frequency being 40 Hz. The true velocity is from 1500 m/s to 4500 m/s (Figure 4.20a). τ is set to 0 in the water layer (Q is infinite). The the true τ except in water layer is from 0.015 to 0.033 (Q is from 75 to 35) (Figure 4.20b). We start the inversion, with an initial velocity model obtained by smoothing the true model ( The hybrid strategy starts with the results of Strategy A after 3 A-iterations, and continues with each 4 C-iterations and 4 A-iterations in an alternate way until the end. We count iteration times as the sum of times of A-iterations and C-iterations in sequence. v p and τ models (Figure 4.23a, b) are slightly updated after 3 iterations. The inverted τ has slightly obvious result than the result of Strategy B. With more iterations, the inverted v p has more detailed structures, and the inverted τ is significantly improved. The hybrid strategy has higher misfit for the first 18 iterations compared with the misfit of Strategy B (Figure 4.24), because the kinematic velocity v k needs to be iteratively updated. With good v k , the final misfit of the hybrid strategy is lower than the misfit of Strategy B.

We perform more iterations for the three strategies Additionally, we also apply the sequential inversion. This strategy includes a first step for updating v p with fixed τ , and a subsequent step for updating τ with fixed v p . The hybrid strategy is composed of a first step for updating v p with fixed τ (v k updating) and the alternate step of Strategy C. We perform the sequential inversion same as the hybrid strategy, except that we replace the alternate each 4 C-iteration with 4 mono-parameter τ inversion with fixed v p . The misfit curves (Figure 4.31) shows that the hybrid strategy has faster convergence and a lower final misfit, compared with the sequential inversion. It indicates that Strategy C guiding the inversion with fixed v k , leads to faster convergence.

Discussion

We incorporate the kinematic relationship to the visco-acoustic inversion. The results show that the hybrid strategy can guide the inversion and mitigate cross-talks. We need to discuss the following aspects: the model searching space, scaling factor between gradients, the Hessian matrix and the regularization terms. First, the explicit kinematic relationship imposes strong constrains on the inversion, which reduces the searching space of model to guide the inversion. The simultaneous inversion has more flexibilities to explore the model space, but still suffers from crosstalks between parameters (e.g., Figure 4.15g, h). The accurate kinematic velocity is important for the hybrid inversion strategy to search in the less model space. The simple models test show that the retrieved parameters have higher resolution (e.g., Figure 4.11), fast convergence and lower misfit values (Figure 4.4,4.13) based on the correct kinematic velocity models. For the 2D synthetic Middle-East model, the main difficulty is to update the kinematic velocity. The result shows that the modified inversion has lower misfit values (Figure 4.26). The retrieved v p has comparable resolution as the one inverted by the simultaneous inversion (Figure 4.25a, c), the retrieved τ has higher average values (Figure 4.25b,d). However, the kinematic velocity for the hybrid strategy is determined for a single frequency, which is the same as the central frequency for a Ricker wavelet source. From the frequency spectrum of the source, the energy of central frequency is dominant. The inversion results of our examples with Ricker wavelet sources show that the central frequency is a proper choice to determine the kinematic velocity.

Second, the scaling factor of the gradients for simultaneous inversion based on conjugate gradient method are fixed before the inversion. The scaling factor α m (Equation 4.25) is used to determine the unitary searching step of each parameter [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF] with a conjugate gradient optimization scheme. For synthetic studies, α m can be obtained by the difference of lower and upper bounds of each model parameter divided by the iteration number [START_REF] He | Parametrization study of the land multiparameter VTI elastic waveform inversion[END_REF]. Therefore, the different choices of α vp and α τ adjust the different sensitivities of parameters. However, the scaling factor is fixed during the inversion. As we discussed before for simple model cases, when the kinematic information is well retrieved, the v p sensitivity is low. Increasing the α vp , the inversion accelerates the v p updating, leading to a faster convergence speed. The factor of ∂τ (v p )/∂v p | v k is a dynamic scaling factor to adjust the sensitivities between parameters.

Third, our comparison between different strategies is based on the conjugate gradient method. Newton method is more accurate in the case for quadratic behaviors, but is very expensive. Visco-acoustic FWI updates two parameters v p and τ , the inverse of Hessian matrix can act as the scaling factor of the gradients of v p and τ and mitigate cross-talks [START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. Due to the ill-posed problem, the Hessian matrix has singular values and it is hard to obtain the inverse Hessian matrix. The Hessian matrix can be solved by second-order adjoint state method, but the difficulty lies in the negative eigenvalues [START_REF] Métivier | Full waveform inversion and truncated Newton method: quantitative imaging of complex sub-surface structures[END_REF]. Another way is to add a regularization term to objective function [START_REF] Asnaashari | Regularized seismic full waveform inversion with prior model information[END_REF]. Then the inverse of Hessian matrix including a damping factor to scale the gradients. [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF] show the different damping factors significantly impact the inversion results. The explicit kinematic relation provides a way to connect the two parameters during the inversion by the term ∂τ (v p )/∂v p | v k . It is a possible alternative to replace the empirical scaling factor.

To illustrate the impacts of the scaling factors, we select different scaling factors for Strategy B. We also apply the term ∂τ (v p )/∂v p for Strategy B as a dynamic scaling factor during each iteration. The true v p model is a homogeneous 2000 m/s with a blocky anomaly 2500 m/s (Figure 4.32a), and true τ model is homogeneous 0.023 (Q = 50) with a blocky anomaly 0.0923 (Q = 15) (Figure 4.32b). The initial v p is 2000 m/s and initial τ is 0.023 (Figure 4.32c,d). Shots are on the top and receivers are on bottom.

The inversion only focuses on reconstructing the blocky v p and τ anomalies. The lower and upper bounds of v p are 2000 m/s and 2500 m/s, and the lower and upper bounds of τ are 0.023 and 0.0923. The α vp and α τ are given by the differences of lower and upper bounds dividing the a chosen a iteration number. For example, we select α vp = 10 m/s and α τ = 0.005, the ratio α τ /α vp is 5 × 10 -4 . The other different ratios can be obtained by choosing different α vp and α τ . The different scaling factors have different convergence speeds (Figure 4.33). The scaling 1 represents the term ∂τ [v p (x)]/∂v p (x) as a dynamic scaling factor, and the scaling 2 represents the spatial mean value of ∂τ [v p (x)]/∂v p (x) (x is the spatial coordinates). Scaling 1 changes with spatial coordinates, but scaling 2 is a scalar value (Figure 4.34).

As we analyzed before, the first iteration mainly focuses on updating the velocity to fit the kinematic information. For the same step of updating v p , a large scaling factor (e.g., 50 × 10 -4 and 15 × 10 -4 ) increases the updating of τ , leading to a smaller misfit. However, a over-estimated large scaling factor constrains v p updating. The scaling factor 50 × 10 -4 (Figure 4.33 dashed dark blue) has the largest misfit after 30 iterations. On the opposite, a over-estimated small scaling factor 0.5 × 10 -4 also constrains τ updating (Figure 4.33 dashed blue). A proper scaling factor (e.g., 2.5 × 10 -4 and 5 × 10 -4 ) is more consistent with the sensitivities between parameters, the inversion has lower misfit values. Among those proper scaling factors, a smaller one, e.g., 2.5 × 10 -4 contributes more on velocity updating when the kinematic information is well matched, which has the lowest misfit value (Figure 4.33 dashed blown). A dynamic scaling factor (Figure 4.33 solid dark and dashed green) has comparable low misfit as fixed scaling factor, e.g., 5 × 10 for each iteration, where x is the spatial coordinates.

Conclusion

We have proposed a hybrid inversion strategy based on modification of the parameterization leading to a different gradient, by explicitly introducing the kinematic relationship. The total gradient is a combination between the partial derivatives of the objective function with v p and τ . ∂τ (v p )/∂v p is the linking term between v p gradient and τ gradient.

The term acts as a dynamic scaling factor. It is derived for a fixed v k and modifies the v p gradient to retrieve the parameters.

To apply the hybrid strategy, we first update the kinematic velocity, then decouple the parameters based on a fixed v k . This strategy is a guided inversion with faster convergence speed. The four tested models are mainly dominated by transmitted waves. The kinematic velocity is relatively easy to be obtained from the transmitted waves, compared with the reflections. With a good approximated v k , the hybrid strategy mitigates cross-talks and has a lower final misfit than conventional strategies, including the sequential and simultaneous inversions.

Appendix I: Introduction to FWI

The theory of full waveform inversion is based on Tarantola (1984a); Plessix (2006b); [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]. Based on Taylor expression, m practically can be decomposed into an initial model m 0 with a perturbation of model ∆m using a local optimization method based on gradient searching.

m = m 0 + ∆m (4.17)
The objective function can be written as a Taylor expansion in vicinity of m 0

J (m 0 +∆m) = J (m 0 )+∇ m J (m 0 )∆m+ 1 2 ∇ 2 m J (m 0 )∆m 2 +O(∆m 3 ). (4.18)
∇ m is the derivative with respect to m. Taking the derivatives with respect to m on the both side of equation 4.18, with the minimum of the misfit function in vicinity of m 0 reached, the left part vanishes and the right part neglect small perturbations, leading to

0 = ∇ m J (m 0 ) + ∇ 2 m J (m 0 )∆m. (4.19) From 4.
19, the perturbed model can be obtain from

∆m = -∇ 2 m J (m 0 ) -1 ∇ m J (m 0 ), (4.20)
where ∇ m J (m 0 ) is the first-order derivative of objective function (gradient), and ∇ 2 m J (m 0 ) is second-order derivative of the objective function called Hessian matrix. In an explicit format, ∇ m J (m 0 ) is written as

∇ m J (m 0 ) = ∂ ∂m 1 2 d cal -d obs 2 , (4.21)
= ∂d cal ∂m † (d cal -d obs ), = J 0 † δd.
δd is the data difference between the calculated and observed data. † is the transpose operator. J 0 is the data derivative with respect to model m, which is also called sensitivity or Fréchet derivative matrix. With the explicit format of ∇ m J (m 0 ), the Hessian matrix is written as

∇ 2 m J (m 0 ) = J 0 † J 0 + ∂J 0 ∂m δd. (4.22)
The Hessian matrix involves a second-order derivative with respect to m, it is computationally expensive. The perturbed model with the explicit formats of ∇ m J (m 0 ) and

∇ 2 m J (m 0 ) is express as ∆m = -J 0 † J 0 + ∂J 0 ∂m δd -1 J 0 † δd. (4.23)
The best solution of m can be obtained by the iteratively updating the ∆m until the objective function converges or satisfies with convergence criterion.

m k+1 = m k -J k † J k + ∂J k ∂m δd k -1 J k † δd k (4.24)
is the recursive optimization. However, the Hessian matrix is hard to obtain due to computation limitation. The inverse of Hessian matrix contains many singularities because of the ill-posed problem. The different ways of incorporating of inverse of Hessian matrix divide the optimization methods into conjugate-gradient, Quasi-Newton, and Newton methods. Conjugate gradient approximates the inverse of Hessian as a positive scalar λ with a conjugate searching direction.

Based on conjugate gradient method, the model updating is obtain by

m k+1 = m k -λ k α m max (|∇ m J (m 0 )|) p k m=m k . (4.25)
α m are scalars called unitary step length, which determine the basic unitary searching length for each parameters. α m are estimated by empirically. For different parameters, α m depend on the gradient scaling. max (|∇ m J (m 0 )|) is maximum of absolute value of first gradient of each parameter to normalize the gradients during each iteration. The λ k is the searching step length which is determined by a parabolic curve fitting method [START_REF] Nocedal | Numerical Optimization[END_REF]. p k is the searching direction determined by conjugate gradient given as

p k = ∇ m J (m k ) + β k p k-1 , (4.26)
β k is used to determine the conjugate direction [START_REF] Polak | Note sur la convergence de méthodes de directions conjuguées[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF] given as

β k = ∇ m J (m k ) -∇ m J (m k-1 )|∇ m J (m k ) ∇ m J (m k ) 2 . (4.27)

Appendix II: Gradient of the objective function

We use L = 5 standard linear solid (SLS) to simulate constant Q with one τ and five τ l σ . As we presented before, 5 SLS produces a good approximation for constant Q from 1.5 Hz to 100 Hz. We rewrite the visco-acoustic wave equation in frequency domain.

         iωp + κ(Lτ + 1) ∂vx ∂x + ∂vz ∂z + L l=1 r l p = f (x s , ω), iωv x + 1 ρ ∂p ∂x = 0, iωv z + 1 ρ ∂p ∂z = 0, (iω + 1 τ l σ )r l p + 1 τ l σ τ κ ∂vx ∂x + ∂vz ∂z = 0, l =∈ [1, L], (4.28)
where κ is the relaxed modulus defined as the ρv 2 p . v p is the P wave velocity and ρ is the density. p, v x v z and associated memory variables r l p are wave fields. The forward modeling is written as a operator format

S(m)u(x, ω) = F, (4.29)
where S(m) is the wave-equation operator in terms of the model parameters m, u(x, ω) is the whole wave fields in whole spatial space x, and F is the source term. The observed data d cal (s, r, ω) is projected to source and receiver positions from the whole u(x, ω).

Mu = d cal (s, r, ω), (4.30)
where M is projection operator defined as dxδ(x -s)δ(xr).

The objective function can be extended as a Lagrangian formulation depending on (m, u, λ). The model parameters include v p and τ , the wave fields include p, v x , v z and r l p . λ are the adjoint state variables associated with p, v x and v z , denoted as λ p , λ vx , λ vz , λ l rp , which are all independent. The Lagrangian formulation is

L(m, u, λ) = 1 2 Mu(x, ω) -d obs (s, r, ω) 2 (4.31) -< λ p , iωp + κ(Lτ + 1)( ∂v x ∂x + ∂v z ∂z ) + L l=1 r l p -f (x s , ω) > -< λ vx , iωv x + 1 ρ ∂p ∂x > -< λ vz , iωv z + 1 ρ ∂p ∂z > - L l=1 < λ l rp , (iω + 1 τ l σ )r l p + 1 τ l σ τ κ( ∂v x ∂x + ∂v z ∂z ) > .
= 0 = M † (Mu(x, ω) -d obs )(s, r, ω) + iωλ p + ∂ ∂x ( λv x ρ ) + ∂ ∂z ( λv z ρ ), ∂L(m,u,λ) ∂vx = 0 = ∂ ∂x [κ(Lτ + 1)λ p ] + iωλ vx + ∂ ∂x L l=1 1 τ l σ τ κλ l rp , ∂L(m,u,λ) ∂vz = 0 = ∂ ∂z [κ(Lτ + 1)λ p ] + iωλ vz + ∂ ∂z L l=1 1 τ l σ τ κλ l rp , ∂L(m,u,λ) ∂r l p = 0 = (iω -1 τ l σ )λ l rp -λ p , l =∈ [1, L].
.

(4.
34) The adjoint equation can be solved in time domain as

             ∂λp ∂t + 1 ρ ( ∂λv x ∂x + ∂λv z ∂z ) = -M † (Mu(x, ω) -d obs ), ∂λv x ∂t + ∂ ∂x [κ(Lτ + 1)λ p ] + ∂ ∂x L l=1 1 τ l σ τ κλ l rp = 0, ∂λv z ∂t + ∂ ∂z [κ(Lτ + 1)λ p ] + ∂ ∂z L l=1 1 τ l σ τ κλ l rp = 0, ∂λ l rp ∂t - λ l rp τ l σ -λ p = 0, l =∈ [1, L]. (4.35)
The adjoint state equation is solved for λ t=T = 0. -M † (Mu(x, ω) -d obs ) is the adjoint source. To satisfy the initial condition, the adjoint state equation is solved with a back-propagation way, from the last time t = T to the initial time t = 0. The last equation in equation 4.35 increases exponentially from time t = 0 to t = T . The back-propagation turns the equation as a exponentially decreasing function, which satisfies the physical meaning due to attenuation.

Finally, the gradient of the objective function becomes Compared with adjoint state approach, the finite difference approach is very expensive. The initial model is a homogeneous model with v p = 2000 m/s. The true model is a two-layer model, 2000 m/s and 2500 m/s, the interface is at the depth of 195 m. The Q model is homogeneous with Q = 50 (τ = 0.023). For the finite difference approach, the perturbation of velocity is 1 m/s. The comparison between the velocity gradients obtained from the adjoint state method and from the finite difference method (Figure 4.35) validates the accuracy of gradients computed from adjoint-state method. 4.9 Appendix IV: Equivalence between two ways for updating v k

As we discussed in Section 4.2.2.4, there are two solutions for updating v k . The first one is that v k is updated by a modified gradient for fixed τ (Equations 4.15,4.16). The other way first starts with mono-parameter v p inversion for fixed τ (Strategy A), and a subsequent v k is obtained based on Equation 4.8. First, we derive the relation between v k and v p for fixed τ . We use the same approach as Appendix 4.8. v p is the function of v k for a fixed τ based on Equation 4.8. Taking perturbations of v k and v p on the two sides, the equation leads to

δv k = ∂f ∂v p δv p + 0 (4.42)
The ratio between perturbations of v k and v p is given by .43) The relation between v k and v p for fixed τ can be written by an explicit way 4.44) where x is spatial locations. The linear relation between v k and v p for fixed τ leads to

δv k δv p = ∂f ∂v p . ( 4
v k (x) = f [v p (x), τ (x)] = v p (x) f [τ (x)], (
v k v p = ∂v k ∂v p τ = f (τ ) = δv k δv p = ∂f ∂v p τ , (4.45)
where | τ means that the τ is fixed. The two ways for updating v k (x) need a initial v 0 p (x), and initial

v 0 k (x) = v 0 p (x) f [τ (x)
]. The solution for updating v k for fixed τ is

v k = v 0 k -H -1 v k ∂J ∂v k τ , = v 0 p f [τ (x)] -f [τ (x)]H -1 vp ∂J ∂v p τ , = f [τ (x)] v 0 p -H -1 vp ∂J ∂v p τ , (4.46)
where we have

∂J ∂v k τ = ∂v p ∂v k τ ∂J ∂v p τ , = 1 f [τ (x)] ∂J ∂v p τ , (4.47)
and

H v k τ = ∂ ∂v k ∂J ∂v k τ , = 1 f [τ (x)] ∂ ∂v p 1 f [τ (x)] ∂J ∂v p τ , = 1 f 2 [τ (x)] H vp τ . (4.48)
The third formula of Equation 4.46 is composed of updating v p for fixed τ

v p = v 0 p -H -1 vp ∂J ∂v p τ , (4.49)
and the following step for updating v k based on Equation 4.44

v k = f [τ (x)]v p . (4.50)
The two solutions are equivalent. In this thesis, we apply the conjugate gradient method without estimation of the Hessian matrix. The model updates are simplified as the searching step multiplying the gradient of objective function.

Introduction

FWI deals with all kinds of waves for inversion, including transmissions and reflections. Information contained in different waves combines the tomography mode and migration mode to build macro and fine structures [START_REF] Mora | Inversion = migration + tomography[END_REF]. Transmissions / diving waves and post-critical reflections are applied to retrieve the long and intermediate wavelengths of the structures, regarded as an efficient scheme of FWI [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF][START_REF] Zhou | Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation[END_REF]. This information is utilized to retrieve velocity and attenuation models using visco-acoustic / visco-elastic FWI [START_REF] Prieux | Multiparameter full waveform inversion of multicomponent OBC data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation[END_REF][START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF]. The 2D synthetic Middle-East model in Chapter 4 mainly focuses on the diving waves to build macro models. The traveltime of transmitted waves are connected with traveling distance and kinematic velocity straightforwardly. It is relatively easier to retrieve the kinematic velocity from transmitted waves, then further update v p and τ . Compared with transmitted waves, the phases of reflections are more strongly affected by the kinematic velocity, depth of interfaces and interference overlaps of reflected waves on thin layers. The phase information explains the kinematic velocity as the combined effects of interface depth, and velocity and τ values inside each layer.

In this chapter, we focus on the reflections to retrieve the velocity and Q in the Marmousi model. The Marmousi model is composed by faults with strong horizontal and vertical velocity variations. Except for direct waves, reflections are the dominated wavefields on this model due to limited offset acquisition. We test the three strategies (Section 5.3) on the Marmousi model to validate the feasibility of the hybrid inversion strategy proposed in Chapter 4.

Marmousi model and acquisition system

The Marmousi model in this section is extracted from the entire Marmousi (Figure 5.1a). The model size is 6 km laterally and 2.2 km in depth with a 200 m water layer at the top. The space grid size is 20 m in both directions. 76 sources are located all over the surface.

The receivers are deployed at the interval of 20 m on the both sides of the source with the maximum offset of 3 km. The true τ model has two anomalies (Figure 5.1b). In terms of Q, the anomaly in light red has a Q value around 22 (τ = 0.052). The other anomaly in dark red has a Q value around 17 (τ = 0.0689). The initial velocity model is obtained by smoothing the true model (Figure 5.1c). The initial τ model is homogeneous (Q is 1000) (Figure 5.1d). The source is a Ricker wavelet with central frequency 4 Hz. The forward simulation time lasts 4 seconds.

Results

Three inversion strategies in this part include:

• Strategy A: mono-velocity inversion for fixed τ value.

• Strategy B: simultaneous inversion for v p and τ .

• Hybrid Strategy (Strategy H): combination of Strategy A and Strategy C. Strategy

A is applied to update the kinematic velocity v k . Strategy C is the modified inversion to decouple v p and τ based on the fixed v k . The Strategy A and C are applied in an alternate sequence. Every 4 A-inversions are alternated with 4-C inversions (the workflow in Chapter 4).

After 5 iterations, the shallow part of v p is the main updated areas for three inversion strategies (Figures 5.2a,5.3a,5.4a). τ is not updated for Strategy A, because it is fixed (Figure 5.2c). The shallow part τ is the main updated area for Strategy B (Figure 5.3b) and the hybrid strategy (Figure 5.4b). More iterations are needed to further analyze the inversions.

After 80 iterations with the three strategies, the shallow part of v p is well-retrieved (Figures 5.2b,5.3c,5.4c). But v p value of the bottom part is lower than the true value, due to weak energy of reflections from the deep layers. A proper pre-conditioner for enhancing gradients of the deep layers could improve the inversion.

We first analyze the retrieved v p of Strategy A. Taking one vertical section from the results with Strategy A at the position x = 2500 m (Figure 5.6), the velocity profile matches well the true model only up to 1 km depth, but some obvious interface shifts occur at depth of 1.4 km and 1.8 km in the v p profile. From the data comparison (Figure 5.10), the phase matches well, but the amplitude of the calculated data shows difference with the observed data without considering τ . This kind of shifts means that the kinematic velocity is translated into the combination of shallow depth and lower v p (for the detailed analysis of retrieved v p with Strategy B and the hybrid strategy, we discuss the results after 400 iterations in later paragraphs.). τ result of Strategy B after 80 iterations shows that the inversion successfully detects the anomalies, except for large artifacts on the bottom (Figure 5.3d). The hybrid strategy after 80 iterations locates the τ anomalies with lower values compared to the true values (Figure 5.4d). Obvious artifacts are observed on the top of the retrieved τ model. More iterations are needed to improve the kinematic velocity v k , and further update v p and τ .

After 400 iterations, v p and τ results (Figure 5.5a, c) are similar to the results after 80 iterations (Figure 5.3c,d) with Strategy B. Inverted v p has a good resolution on the shallow part, but under-estimated values on the bottom. Inverted τ has strong artifacts on the bottom. The v p result of the hybrid strategy after 400 iterations (Figure 5.5b) is similar to the previous result after 80 iterations. The τ result (Figure 5.5d) shows that the two anomalies are retrieved, but some artifacts are distributed on both the top and bottom parts. The τ anomalies on the top part is lightened.

We extract 3 sections from the results of Strategy B and the hybrid strategy at the position x = 2500 m, 3320 m and 4000 m, respectively (Figure 5.7). The three positions correspond to the areas, without τ anomaly, with the left τ anomaly and the right τ anomaly, respectively. We first analyze the v p results. The v p profiles of Strategy B (Figure 5.7a, b, c) and the hybrid strategy (Figure 5.7d, e, f) show the retrieved v p matches well the true model on the shallow part. On the bottom part, the retrieved v p is lower than the true value. The thin layers are not reconstructed due to lacking high frequencies. The τ result of Strategy B (Figure 5.7g,h,i) shows severe variations along the depth, although the structure of true τ is smooth except for two anomalies. The inverted two τ anomalies have larger values and a shift on the bottom interface of the anomalies. We observe that the lower under-estimated v p value is consistent with the larger τ value on the bottom, e.g., the τ artifacts around 1.3-1.5 km at position x = 2500 m (Figure 5.7g), 1.6-1.8 km at position x = 3320 m (Figure 5.7h) and 1.7-1.9 km at position x = 4000 m (Figure 5.7i). The structure of inverted τ has relevance to v p structure, which shows cross-talks between two parameters. The τ result of the hybrid strategy (Figure 5.7j,k,l) shows that the two anomalies are located around the correct position but with a vertical deeper shift. The anomaly value is close to the true model, especially for the anomaly at position x = 3320 m. The vertical shift is caused by the under-estimated low kinematic velocity v k . The τ structure is related to the ratio of v k and v p based on kinematic relationship in Equation 4.8 (Figure 5.7m,n,o). For example, the anomaly at position x = 3320 m, the τ anomaly is located around 1.3-1.5 km (Figure 5.7k). The interfaces of v p and v k at depth of 1.3 km match well (Figure 5.7n), leading to the matches between observed and calculated waves reflected on this interface. To match the waves reflected on the bottom interface, the v k has larger value than v p , and a downshifted interface, leading to data fitting (the arrows in Figure 5.7n). We simulate one shot gather at position x = 3320 m by using v p = v k and by setting τ as 0. The phase of calculated data matches well the observed data, except for the amplitude difference (Figure 5.13).

From the misfit curves for the 3 strategies after 400 iterations (Figure 5.8), the misfit of Strategy A stops decreasing after 80 iterations due to the lack of τ updating. Strategy B has the smallest final misfit. The misfit of the hybrid strategy decreases slowly until a close misfit as Strategy B. The slow decreasing is due to the v k updating. One shot data comparison at horizontal position x = 3320 m shows the large misfit of Strategy A (Figures 5.9 and 5.10). The final data misfit for Strategy B (Figure 5.11) and the hybrid strategy (Figure 5.12) is lower than 1%, and the data difference is mainly distributed between the later arrivals. The good data matches of both Strategy B and the hybrid strategy indicate that there are non-unique solutions, due to cross-talks between parameters, limited acquisition system, resolution of limited frequencies and so on. We compare the two simulated shot gathers (Figures 5.15 and 5.16). The data difference is very small. We need to investigate the τ resolution with small data difference.

Results

Velocity

To estimate the τ resolution [START_REF] Malinowski | High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full waveform inversion[END_REF], we add checkerboard perturbations on the retrieved τ model by Strategy B after 400 iterations (Figure 5.5b). The retrieved τ by simultaneous inversion contains strong artifacts on the bottom of the model, thus we add τ perturbations on that part, with the value being 0.01 (Q is around 100) (Figure 5.17), is lower than the τ anomalies (τ is around 0.06). We apply the mono-parameter τ inversion for fixed v p . The retrieved τ model below the dashed line lacks resolution due to the limited acquisition coverage. Those less-resolution checkerboard artifacts can partly explain the artifacts of the reconstructed τ below the dashed line. Other reasons for strong artifacts in retrieved τ model (Figure 5.5b) are due to cross-talks between parameters, as we explained before.

In summary, the data difference is very small between the shot gathers simulated with / without the artifacts removed (Figures 5.15 and 5.16). But, the checkerboard test validates the high resolution of τ under the acquisition coverages by mono-parameter τ inversion (Figure 5.17) on small data difference. We conclude that the small data difference is difficult to reconstruct both v p and τ correctly, due to cross-talks with the

Over-estimation and cross-talks

From the results of the three strategies in section 5.3, v p and τ have cross-talks. The under-estimated lower v p transfers the phase difference to τ , leading to a over-estimated τ (τ artifacts) in Strategy B and the hybrid strategy. Meanwhile, Strategy A overestimates v p without updating τ . We need to understand the inversion behaviors on over-estimation and cross-talks between parameters. We use different v p models as input for mono-parameter τ inversion, and apply a simultaneous inversion over v p and τ , including:

• Mono-τ 1: The initial v p is the inverted result with Strategy A after 100 iterations.

• Mono-τ 2: The initial v p is the true v p model.

• Mono-τ 3: The initial v p is the inverted result with Strategy B after 100 iterations.

• Simultaneous: The initial v p is the inverted result with Strategy A after 100 iterations.

Mono-τ denotes the mono-parameter τ inversion for fixed v p . Three v p models are inverted v p with Strategy A after 100 iterations, true v p model and inverted v p with Strategy B after 100 iterations. τ is almost not updated by mono-τ 1 (Figure 5.24), due to the incorrect over-estimated v p model with a coupling of incorrect depth and velocity values (see the v p profile in Figure 5.6). The kinematic information is mainly dominated by the velocity model v p . There is not updating space for τ , because the inversion is stuck into the local minimum. With the same v p model as input, the simultaneous inversion can not update v p and τ (Figure 5.27) towards the true models. We observe that the retrieved v p is larger than the true model at the locations of the two τ anomalies. The τ anomalies are not retrieved. To match the data, the under-estimated τ values transfer the kinematic information to over-estimated larger v p . The over-estimated incorrect v p makes the inversion be stuck into the local minimum, leading to no updating space for τ . The mono-τ 2 uses the true v p model as input, the retrieved τ locates the anomalies and has identical structure as the true τ model (Figure 5.25). Meanwhile, the mono-τ 3 has the similar τ result (Figure 5.26) as Strategy B with many artifacts due to the under-estimated lower v p .

From the three mono-τ inversions, τ updating is strongly affected by the v p models. The incorrect v p leads to the limited updating space for τ , even with a simultaneous inversion. The under-estimated v p maps the phase difference into the over-estimations or artifacts on the τ result. The cross-talks between v p and τ make the inversion get stuck into a local minimum.

Frequency continuation inversion

The kinematic velocity v k is the function of frequency, v p and τ (Equation 4.8). The kinematic velocity v k is determined at the central frequency of the source. The central frequency is defined as the peak frequency in all the previous examples. Concerning the broadband Ricker wavelet, we filter the Ricker wavelet into two narrow bandwidth sources. The Ricker wavelet in Section 5.3 is a 4 Hz central frequency wavelet. The maximum frequency is 10 Hz (Figure 5.28a). We apply the two band-pass filters to have a low and high frequency source term. The spectrum of the low frequency source is from 1.5-5 Hz (Figure 5.28b). The central frequency is 3.8 Hz. The spectrum of the high frequency source is from 4.5-6 Hz (Figure 5.28c). The central frequency is 4.8 Hz. We choose 3.8 Hz and 4.8 Hz for the determination of v k . The low frequency and high frequency data both contain bandwidth of several frequencies. The high frequency data has 0.5 Hz bandwidth overlapping with low frequency data. This overlapping frequency continuation scheme is more stable than single frequency continuation scheme [START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF]. We first apply the low frequency data with both Strategy B and the hybrid strategy for 40 iterations. Then the high frequency data is exploited to obtain the final results. Strategy B after first 40 iterations retrieves the main v p structures (Figure 5.29a). With higher frequency data after another 40 iterations, more fine structures are retrieved, especially on the shallow part (Figure 5.29c). For the result of τ at low frequency, two anomalies are retrieved (Figure 5.29b). There are strong artifacts on the bottom for high frequency data (Figure 5.29d). For the hybrid inversion, v p retrieved by high frequency data (Figure 5.30c) has more more detailed structure than v p retrieved by low frequency data (Figure 5.30a). The two τ anomalies are retrieved by low frequency data (Figure 5.30b). There are some artifacts on the shallow part. After another 40 iterations with high frequency data, two anomalies are more clearly located (Figure 5.30d). τ artifacts on the shallow part is relieved with more iterations. Some of the artifacts are explained by the τ checkerboard result, showing the lack of resolution below the dashed line (Figure 5.17), due to the insufficient acquisition coverage. With the same iterations, τ anomalies are more obvious than τ result without data filtering (Figure 5.4). The misfits of Strategy B and the hybrid strategy both decease to be below 1% (Figure 5.31).

To obtain the results with 1% data misfit left (Figure 5.8), we need 170 iterations with the hybrid strategy without data filtering. The frequency continuation scheme needs 50 iterations to have the same data misfit. The frequency continuation scheme is around 3 time faster than previous one without data filtering to retrieve the τ anomalies with the same data misfit until 1%. The frequency continuation scheme has more accurate v k , which is determined at the central frequency of each narrow frequency bandwidth, leading to a faster convergence.

Conclusions

We have applied the three strategies to the Marmousi model. The data is mainly composed of reflections, except for direct waves in the water layer. The v p is over-estimated with Strategy A without considering τ , ending with the combination of incorrect depth and velocity values. The incorrect over-estimated v p makes the inversion get stuck into a local minimum. The mono-parameter v p cannot perfectly match the data without updating τ .

Strategy B updates v p and τ simultaneously. The two τ anomalies are retrieved, but the structure of inverted τ is affected by the structure of inverted v p due to cross-talks. The shallow part of v p is well-retrieved, meanwhile the bottom part of v p is lower than true model. The under-estimated low v p has cross-talks with τ , leading to strong artifacts on τ . This inversion behavior is validated by three mono-τ inversions. The deeper τ artifacts appear on the locations where v p is insufficiently retrieved. Meanwhile, the mono τ inversion with the true v p model as input has few artifacts. On the one side, the impacts of τ artifacts on the data are small, but this small data difference cannot correctly reconstruct v p and τ at the same time. The simultaneous inversion ends with the final result of combination between strong τ artifacts and under-estimated low v p values on the bottom part. On the other side, the checkerboard test indicates that the monoparameter τ inversion has high resolutions on τ perturbations under the illustration zone controlled by acquisition coverage. The two aspects indicate that there are cross-talks between parameters. The inversion has to be guided to explain the small data difference.

Compared with results of Strategy B, τ artifacts with the hybrid strategy are lower, especially on the bottom part. The two anomalies are retrieved by the hybrid strategy, but with a vertical depth shift. The kinematic of final calculated is controlled by the v k model, it matches well the kinematic of observed data. To match the waves reflected on the bottom interfaces of the τ anomalies, v k has large values and deeper interfaces. τ is further updated based on the ratio of v k and v p to match the phase and amplitude, leading to a vertical depth shift. This kinematic relationship assumes that v k is a function of a single frequency, which is determined at the central frequency of source during the inversion. From the inversion results, this assumption proves to be proper and effective to retrieve v p and τ . Frequency continuation scheme is also recommended to obtain more accurate v k for the hybrid inversion strategy.

Finally, the hybrid inversion strategy incorporates the kinematic relationship to guide the inversion. Compared with simultaneous inversion, the hybrid strategy retrieves a good velocity model (mainly focus on the top) and locates the attenuation anomalies with less artifacts. This is an alternative way to retrieve v p and τ models on complex models, e.g., Marmousi model that dominated by reflections, and insensitive to Gaussian noise with L 2 norm objective fucntion.

Chapter 6

Conclusions and perspectives

Conclusions

In this thesis, I have investigated different inversion strategies for visco-acoustic full waveform inversion to reconstruct velocity and attenuation models. Visco-acoustic FWI belonging to the multi-parameter inversion suffers from cross-talks between parameters. Conventional inversion strategies include sequential and simultaneous inversion strategies. Sequential inversion includes a first step of inverting the velocity model, and a subsequent step for updating the attenuation model. This approach may lead to parameter over-estimation. In this thesis, we apply the mono-parameter velocity inversion to investigate the inversion behaviors. The simultaneous inversion inverts velocity and attenuation at the same time. But it suffers from cross-talks between parameters without the introduction of the Hessian matrix.

In this thesis, I have proposed a hybrid inversion strategy to guide the inversion without relying on the information of the Hessian matrix and to mitigate cross-talks. This inversion strategy includes two key steps: (1) the determination of the kinematic velocity and (2) the decoupling between velocity and attenuation based on the kinematic relationship.

Hybrid strategy -A guided inversion

The hybrid inversion strategy incorporates the kinematic relationship to guide the inversion and mitigate cross-talks. To apply the hybrid strategy, we first update the kinematic velocity, then decouple the parameters based on a fixed v k . For the fixed kinematic velocity, the total derivative of the objective function with respect to v p is a combination between the partial derivatives of the objective function with respect to v p and τ . ∂τ (v p )/∂v p is the linking term between v p and τ gradients. The term acts as a dynamic scaling factor to balance the gradient sensitivities without the information of the Hessian matrix. The hybrid strategy follows the kinematic curves to decouple the parameters.

I applied the hybrid strategy to transmitted-wave-dominated simple models and to the 2D Middle-East model (Chapter 4) and to reflected-wave-dominated Marmousi model (Chapter 5). I have compared the approach with conventional inversion strategies. The key factor of the hybrid strategy is to update the kinematic velocity. The kinematic velocity is relatively easy to invert in transmitted-wave-dominated models, leading to a faster converge speed and lower final data misfit, compared with conventional simultaneous inversion. The reflection-dominated Marmousi model needs more iterations to update the kinematic velocity. The hybrid inversion strategy successfully reconstructs the attenuation anomalies with lower artifacts, and has a comparable low final data misfit, compared with the conventional simultaneous inversion. For the 2D Middle-East model, 18 iterations with the hybrid strategy leads to a lower misfit than the misfit of the simultaneous inversion. But for the Marmousi model, more than 100 iterations with the hybrid strategy are needed to obtain a comparable low misfit as the one provided by the simultaneous inversion. The frequency continuation scheme is more efficient to decouple the parameters, with the kinematic velocity being determined within narrow bandwidth.

In conclusion, the hybrid inversion strategy is an alternate inversion strategy to reconstruct velocity and attenuation models without the information of the Hessian matrix.

Perspectives

In this thesis, we focus on inverting velocity and attenuation models. The hybrid inversion strategy could extend to 3D and apply to real data. Furthermore, with the attenuation models, we also can incorporate Q model a passive parameter to build macro velocity through image-domain method and invert for velocity perturbation by reversetime migration.

Extension to 3D

We implemented the time domain visco-acoustic FWI in 2D media. Considering the 3D media, the computation costs and memories consuming are significant large. First, the visco-acoustic wave equation involves with extra equations in a total number L for computing memory variable r l p , compared to pure acoustic wave equation. L is commonly chosen as 3 to have a constant Q [START_REF] Plessix | Visco-acoustic full waveform inverision[END_REF][START_REF] Yang | A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion[END_REF]. Second, the 3D media has one extra dimension, leading to significant longer computation time. Third, the gradient is obtained by time integration of the cross-correlation between forward-propagation and back-propagation wavefields at each time step. The intuitive way is to store the forward-propagation wavefields, which leads to very large memories. Careful memory optimization is needed, e.g., checkpointing strategy [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation[END_REF][START_REF] Griewank | Algorithm 799: Revolve: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation[END_REF][START_REF] Symes | Reverse time migration with optimal checkpointing[END_REF][START_REF] Anderson | Time-reversal checkpointing methods for RTM and FWI[END_REF][START_REF] Komatitsch | Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion[END_REF]. The checkpointing strategy allows the forward-propagation wavefields to be reconstructed within a time window with stored checkpointing snapshots. [START_REF] Yang | Wavefield reconstruction in attenuating media: A checkpointing-assisted reverse-forward simulation method[END_REF] combined the checkpointing strategy with the reconstruction of forwardpropagation wavefields by reverse propagation in time. The forward-propagation wavefields are recomputed from final snapshot with the stored boundary conditions. This method can efficiently stabilize the reverse propagation in the attenuating media to correctly reconstruct forward-propagation wavefields.

Application to real data

For the synthetic study, we assume that we know the source in advance. But for real data, we need to estimate the source. Due to the attenuation, high frequency parts are absorbed, leading to the difficulty for estimating accurate source. Besides, the more accurate physical models are needed, e.g., the role of density and the anisotropy parameters. Limited acquisition system and noise should also be considered for the real data case.

Image domain macro velocity building with Q model

Image domain methods retrieves the macro velocity model by focusing the image panel [START_REF] Chauris | Analyse de vitesse par migration pour l'imagerie des structures complexes en sismique réflexion[END_REF][START_REF] Sava | Wave-equation migration velocity analysis[END_REF][START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF][START_REF] Cocher | Seismic iterative migration velocity analysis: two strategies to update the velocity model[END_REF]Chauris and Cocher, 2017;Chauris et al., 2017;[START_REF] Li | Subsurface seismic imaging based on inversion velocity analysis in both image and data domains[END_REF]. Image domain can be implemented in common-shot or subsurface gathers [START_REF] Symes | Extended Waveform Inversion[END_REF]. The image domain methods back-project the image residues to update the macro models [START_REF] Sava | Wave-equation migration velocity analysis[END_REF][START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF][START_REF] Sava | Extended imaging conditions for wave-equation migration[END_REF][START_REF] Shen | Wave-equation migration Q analysis[END_REF]. Inaccurate velocity leads to a defocused energy of image panel (Figure 1.6). Based on this criteria, the macro velocity model can be updated for a known fixed Q model.

Migration velocity analysis (MVA) [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF] and inversion velocity analysis (IVA) (Chauris and Cocher, 2017) optimize the common image gathers (CIG, image panels) to update the macro velocity model. Compared with MVA, IVA uses a direct pseudo inverse of the Born operator in shot-domain [START_REF] Li | Coupling direct inversion to common-shot image-domain velocity analysis[END_REF] or subsurface domain (Chauris and Cocher, 2017), instead of the adjoint Born operator. Only when the background velocity model is accurate, energy of CIGs will be focused. When the data is affected by attenuation, the attenuation and inaccurate macro velocity lead to dispersed CIGs. However, different combinations of velocity and Q values lead to the same kinematic velocity at the reference frequency. MVA or IVA only takes account of kinematic information. CIGs can be focused with the equivalent kinematic velocity when data is affected by attenuation, which means velocity and Q can not be updated simultaneously by MVA or IVA. To validate this, we apply IVA for macro velocity building with the visco-acosutic data to invert the kinematic velocity model. Then we incorporate the Q model as a passive parameter to invert the macro velocity model.

IVA for visco-acoustic data

We apply IVA for three datasets in subsurface domain (Chauris and Cocher, 2017). The criterion for evaluating focusing or dispersion of the CIG panel is

J (v 0 ) τ = 1 2 hI(x, h) 2 , (6.1)
where, v 0 is the background velocity model, h is the horizontal shift in subsurface offset, and x is the model space with horizontal distance x and depth z for 2D case. I(x, h) is the CIG. The accurate background macro velocity leads to focused CIG. Three observed datasets are acoustic data, visco-acoustic data with Q = 50 and viscoacoustic data with Q = 15. The central frequency of source is 20 Hz. The true velocity model is a two-layer model, the first layer velocity is 2000 m/s and the second layer velocity is 2100 m/s. The kinematic velocity of the first layer for the three datasets are 2000 m/s, 2060 m/s and 2234 m/s, respectively. We start IVA with the three datasets in a initial background velocity v p = 2500 m/s.

We perform 6 iterations in pure acoustic media to update the background macro velocity. After 6 iterations, the energy of CIGs is focused (Figure 6.1). The background macro velocity models are updated with the 3 datasets (Figure 6.2). Without the information from the second layer, IVA only can retrieve the kinematic velocity of first layer. The inverted average kinematic velocity of the first layer is 2042 m/s, 2098 m/s and 2345 m/s, respectively. The results coincide with the predicted kinematic values for three models. It shows that MVA or IVA can not update the macro velocity and Q models independently. Thus, we need to incorporate the Q model as a passive parameter to update the macro velocity model.

MVA with Q model

We incorporate Q model as a passive parameter to reconstruct the macro velocity model by MVA (Equation 6.1). But for the MVA in visco-acoustic media, we need to derive the equations for fixed Q.

First, the common image gather I(x, h) in the subsurface offset is obtained from adjoint Born operator with a horizontal shift h, given by I(x, h) = -1 2ρv 2 0 xs xr ω (Lτ + 1)∇ • v(x s , x -h, ω)q(x r , x + h, ω)

+ L l=1 τ τ l σ ∇ • v(x s , x -h, ω)s l (x r , x + h, ω) , (6.2)
where, p, v x , v z and r l p are forward-propagation wave fields, and q, u x , u z and s l are back-propagation wavefields. Forward-propagation wave fields are obtained by viscoacosutic wave equation, in the frequency, given as (6.4)

The gradient of ∂J /∂v 0 | τ is obtained by the adjoint state method (Plessix, 2006b). The objective function is extended as a Lagrangian formulation depending on (m, U , λ).

The model parameter m is v 0 , U is the wavefield, which include forward-propagation wave fields p, v x , v z and r l p , back-propagation wave fields q, u x , u z and s l and CIG I. λ are the adjoint state variables associated with p, v x , v z , r l p , q, u x , u z , s l and I, which are all independent. The Lagrangian formulation is L(m, U , λ) = )s l -q = 0, l =∈ [1, L] (6.19) where λ 1 , λ 2 , λ 3 and λ l 4 are written as q, u x , u z and s l respectively. Therefore, λ are obtained by solving the adjoint state equations (Equations 6.17,6.18 and 6.19) (Lτ + 1)∇ • v(x -h)q(x + h)λ 9 (x, h)

- 2 v 4 0 h L l=1 τ τ l σ ∇ • v(x -h)s l (x + h)λ 9 (x, h) .
For example, we have a two-layer model with velocity being 2100 m/s and 2200 m/s. The interface is at the depth of 51. Q model is homogeneous with the value being 1000. Initial velocity is 2500 m/s and input Q is 1000. We obtain the gradient with respect to the macro velocity (Figure 6.3). The blue area indicates the negative gradient to decrease the velocity values.

Inverting velocity perturbation with Q model

The velocity model can be divided into a background macro model and a velocity perturbation under the Born approximation. The velocity perturbation is also called seismic image in reflection seismology. Seismic image characterizes the structures of the geological bodies, which is useful for interpretation. Based on Born approximation, the velocity perturbation leads to perturbed data. The velocity perturbation is inverted with the adjoint Born operator applied to the perturbed data [START_REF] Claerbout | Earth soundings analysis: Processing versus inversion[END_REF]. It can be implemented by wave-equation reverse-time migration (RTM) [START_REF] Mcmechan | Migration by extrapolation of time-dependent boundary values[END_REF][START_REF] Baysal | Reverse time migration[END_REF]. To compensate for insufficient illumination and geometrical spreading effects, true-amplitude migration is applied to retrieved the velocity perturbation [START_REF] Zhang | True amplitude wave equation migration arising from true amplitude one-way wave equations[END_REF][START_REF] Zhang | Theory of true-amplitude one-way wave equations and true-amplitude common-shot migration[END_REF][START_REF] Zhang | Practical issues in reverse time migration: true amplitude gathers, noise removal and harmonic source encoding[END_REF]Zhang et al., 2014b). Seismic attenuation absorbs wave energy, leading to imaging blurring under the attenuation area. There are two ways to improve the seismic images. The first way is Q compensated RTM, which is implemented on the pre-stacked data by enhancing the wave energy [START_REF] Zhang | Compensating for visco-acoustic effects in reverse-time migration[END_REF]. The other way is to apply the least-squares RTM (LS-RTM) to iteratively update the velocity perturbations with Q being a passive parameter (Zhang et al., 2014a;[START_REF] Dutta | Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation[END_REF][START_REF] Zeng | A guide to least-squares reverse time migration for subsalt imaging: Challenges and solutions[END_REF].

Q least-squares RTM

The least-squares RTM is an iterative scheme to update the seismic images (Zhang et al., 2014a;[START_REF] Dutta | Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation[END_REF][START_REF] Zeng | A guide to least-squares reverse time migration for subsalt imaging: Challenges and solutions[END_REF]. Q is incorporated as a passive parameter. LS-RTM is similar to the mono-parameter velocity inversion for visco-acoustic FWI with fixed true Q model, except for a linearized inversion [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations: Conference on Inverse Scattering, Theory and Application[END_REF][START_REF] Lines | Tutorial: A review of least-squares inversion and its application to geophysical problems[END_REF]Tarantola, 1984b). This aspect is not discussed in detail.

6.2.4.2 Q compensated RTM Conventional methods of compensating attenuation are implemented in the data domain by an inverse Q-filter [START_REF] Wang | Seismic Inverse Q Filtering[END_REF]. This Q-filter is based on a 1D Q model, which is hard to handle complicated geological models. As the visco-elastic / viscoacoustic simulation becomes achievable, the Q compensated RTM was first proposed by [START_REF] Zhang | Compensating for visco-acoustic effects in reverse-time migration[END_REF]. The basic idea is to compensate for the energy absorptions during the wave extrapolation (Figure 6.4). The seismic images are obtained by the zero-lag cross-correlation between source side forward-propagation with receiver side back-propagation wave fields [START_REF] Claerbout | Toward a unified theory of reflector mapping[END_REF]. Zhu and Harris (2014) and Zhu et al. (2014) extended the fractional derivative visco-acoustic wave equation into a dispersionabsorption-split wave equation. Through changing the absorption term, the Q compensated extrapolation is achieved. However, the simulation needs to solve the fractional derivative with space Fourier transform, which is computational expensive. Based on our previous analysis of kinematic velocity (Section 2.6.3), we modifies the viscoacoustic wave equation (Equation 2.42) to fix the kinematic velocity and to compensate the absorption. We rewrite the visco-acoustic wave equation into the second-order form with constant density, given as (6.22) where φ is equal to 1 for the visco-acoustic wave equation. As we analyzed before, the real term of the complex velocity is related to the kinematic effects. The imaginary term of the complex velocity corresponds to the energy absorption with a positive φ. Therefore, a negative φ corresponds to energy compensation. To obtain the compensated image condition, φ turns into -1 to amplify the energy. The first-order velocity-pressure equations for forward modeling is where q, u x , u z and s l are back-propagation wave fields. The r l p in the third line of Equation 6.23 is exponentially increasing with time, and s l in the third line of Equation 6.24 is exponentially increasing with reverse-time. The increases of r l p and s l stand for the Q compensation. The amplifying leads to unstable computations, especially for higher frequency noise during extrapolation. A low-pass filter is needed to stabilize the computations [START_REF] Zhang | Compensating for visco-acoustic effects in reverse-time migration[END_REF]Zhu et al., 2014). Finally, we implement the Q compensated RTM with the visco-acoustic wave equation based on nearly constant Q model (the results are not presented here).

  attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Attenuation phenomenon and impacts . . . . . . . . . . . . . . 1.1.2 Intrinsic and scattering attenuation . . . . . . . . . . . . . . . . 1.1.3 Different Q models . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Estimation of attenuation . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Ray-based methods . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Wave equation methods . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Data domain and image domain methods . . . . . . . . . . . . 1.3 Motivation, strategy and thesis organization . . . . . . . . . . . . . . . 1.3.1 Motivation: Inversion strategies to reconstruct velocity and attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Strategy: Kinematic relationship to guide the inversion . . . . . 1.3.3 Thesis organization . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . Modeling 2 Visco-elasticity and visco-acoustic wave equation 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Visco-elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Constant Q approximation . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Kjartansson constant Q . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Nearly Constant Q model . . . . . . . . . . . . . . . . . . . . . 2.3.3 Modified optimization of nearly constant Q . . . . . . . . . . . 2.4 Visco-acoustic wave equations . . . . . . . . . . . . . . . . . . . . . . vii List of Figures1.1 Q image of subsurface superimposed to velocity model in an offshore Malaysia reservoir; the purple area indicates the gas clouds (adapted from Plessix et al. (2016)) . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Migrated image without (left) and with compensation for attenuation (right) (adapted from Zhou et al. (2011)) . . . . . . . . . . . . . . . . . 1.3 Scattering attenuation effects in highly heterogeneous media: the amplitude decays and traveltime fluctuates (Adapted from Sato et al. (2012)); the coda waves are made up of incoherent scattered waves. . . . . . . . 1.4 Intrinsic attenuation effects for one seismic trace in time domain (left) and normalized spectrum in frequency domain (right), the solid red line is the acoustic wave, and the solid blue line is the visco-acoustic wave with Q = 15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Q measurements in Pierre shale. The open circles are field experiments measured by McDonal et al. (1958), the solid line is the fitting curve between α and frequency, the dotted and dashed lines are computed by numerical modeling method from Carcione et al. (2002) (Adapted from Carcione (2011). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Image panel (i.e, common image gather) in subsurface domain, a lower background velocity (a), a correct background velocity (b) and a higher background velocity (c) (Adapted from Lameloise et al. (2015). . . . . . 1.7 Effects of cross-talks between P-wave velocity v p and τ (approximates to the inverse of Q), J is the FWI objective function. The true v p model is homogenous v p = 2000 m/s with an low value anomaly 1800 m/s (a), and the true τ model is homogeneous τ = 0.001 with a high τ anomaly 0.07 (b). Both gradients with respect to v p (c) and τ (d) contain information on τ and v p , respectively (circles). . . . . . . . . . . . . . . . . . . . . 1.8 Two normalized recorded seismic traces with velocity and attenuation being (v p = 2234 m/s, Q = infinity) (red) and (v p = 2168 m/s, Q = 50) (blue). The associated data are almost identical. . . . . . . . . . . . . .

  40 2.10 Numerical solution compares with analytical solution and at distance 200 m (left) and 800 m (right) with the spatial grid space being 5 m. Top: trace comparisons, bottom: trace difference. . . . . . . . . . . . . 42 List of Figures xiii 2.11 20 Hz central Ricker wavelet in the time (left) and frequency (right) domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Spectrum of the recorded trance by analytical and numerical solutions, the central frequency of trace is 17.5 Hz computed by analytical solution (solid red) and finite difference method (dashed blue), lower than the source central frequency with 20 Hz. . . . . . . . . . . . . . . . . . . . 2.13 Trace comparison between acoustic (red) and visco-acoustic (blue) waves in time domain (left panel) and frequency domain (right panel). . . . . . 2.14 Snapshot comparison for acoustic (left parts of the images) and viscoacoustic (right parts) waves at 200 ms (a), 500 ms (b), 800 ms (c) and 1000 ms (d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Left: different kinematic velocity v k at 28 Hz for three kinds of waves: (infinite Q, v p = 2000 m/s) in red, (Q = 50, v p = 2000 m/s) in black and (Q = 15, v p = 2000 m/s) in green. Right: same kinematic velocity v k = 2250 m/s at 28 Hz for three kinds of waves: (infinite Q, v p = 2250 m/s) in red, (Q = 50, v p = 2180 m/s) in black and (Q = 15, v p = 2000 m/s) in green. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Traces for three different Q and velocity combinations (the same configurations as Figure 2.15 right): (infinite Q, v p = 2250 m/s) in green, (Q = 50, v p = 2180 m/s) in black and (Q = 15, v p = 2000 m/s) in red at the near offset (left) and far offset (right), in the time (top) and frequency domains (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.0 Elastic / rheological models and constitutive relations in the time domain, except GZB in the frequency domain due to the complexity in the time domain, M l is the elastic modulus of spring, η l is the viscosity coefficient and G is the coefficient of spring-pot. . . . . . . . . . . . . . 4.1 Different combinations of v p and τ , leading to the same v k = 2250 m/s at reference frequency 4, 20 and 40 Hz (left) and comparison of three shot gathers for different combinations of v p and τ at 20 Hz (right), the traveltime is same and amplitude becomes weaker as attenuation increases. 4.2 True velocity (a) and τ (b) models, initial velocity (c) and τ (d) models. 4.3 Updated velocity (a) and τ (b) models after 3 iterations (Strategy A), updated velocity (c) and τ (d) models after 30 iterations (Strategy A), updated velocity (e) and τ (f) models after 3 iterations (Strategy B), updated velocity (g) and τ (h) models after 3 iterations (Strategy B), updated velocity (i) and τ (j) models after 3 iterations (Strategy C). . . . 4.4 Misfit comparison for three different strategies: Strategy A (red), Strategy B (green), and Strategy C (blue). . . . . . . . . . . . . . . . . . . . xiv List of Figures 4.5 Search direction for Strategies A, B and C, the blue line is the kinematic curve for all combinations of v p and τ equal to v k 2060 m/s, red solid line is the searching direction of Strategy B. The dots are the mean values of each iteration results. . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.6 Data comparison of Strategy B after 2 (a), 3 (b), 4 (c) iterations, in each subfigure, the upper left is the real data, upper middle is the calculated data, upper right is the data difference, the bottom left is one trace comparison at 201 m, the bottom middle is one trace comparison at 321 m, the bottom right the one trace difference at 201 m and 321 m respectively. The phase of data after 2 (a) and 4 iterations (c) are matched, meanwhile the phase of data after 3 iterations (b) has a shift. . . . . . . 80 4.7 Same as Figure 4.6, except for data comparison of Strategy C after 1 (a), 2 (b), 3(c), the phase of data is kept during the inversion. . . . . . . . . 82 4.8 Gradients of v p and τ for Strategy B, the left panels (a) and (c) are gradients after 2 iterations, and the bottom two panels (b) and (d) are gradient after 3 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.9 Same as Figure 4.8, but for gradients of v p and τ of Strategy C. The left panels (a) (c) and (e) are gradients after 1 iteration, and the right panels (b) (d) and (f) are gradient after 2 iterations. The full derivative dJ /dv p (e, f) are positive for the area under the acquisition coverage. . . . . . . 84 4.10 True v p (a) and τ (b) models, initial v p (c) and τ (d) models. . . . . . . . 86 4.11 Inversion results: updated v p (a) and τ (b) models after 3 iterations (Strategy A), updated v p (c) and τ (d) models after 30 iterations (Strategy A), updated v p (e) and τ (f) models after 3 iterations (Strategy B), updated v p (g) and τ (h) models after 3 iterations (Strategy B), updated v p (i) and τ (j) models after 3 iterations (Strategy C), updated v p (k) and τ (l) models after 30 iterations (Strategy C). The dashed lines in (i), (j), (k) and (l) denote the position of the v p and τ anomalies. . . . . . . . . 88 4.12 Kinematic velocity v k model for Strategy C, it is obtained by smoothing v p (Figure 4.11a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.13 Misfit curve of 3 strategies (left), and enlarge in details (right). . . . . . 88 4.14 True v p (a) and τ (b) models, initial v p (c) and τ (d) models. . . . . . . . 90 4.15 Inversion results: updated v p (a) and τ (b) models after 3 iterations (Strategy A), updated v p (c) and τ (d) models after 30 iterations (Strategy A), updated v p (e) and τ (f) models after 3 iterations (Strategy B), updated v p (g) and τ (h) models after 3 iterations (Strategy B), updated v p (i) and τ (j) models after 3 iterations (the hybrid strategy), updated v p (k) and τ (l) models after 30 iterations (the hybrid strategy). . . . . . . . 92 4.16 Misfit curves for Strategy A (solid red), Strategy B (solid green), the hybrid strategy (dashed blue). . . . . . . . . . . . . . . . . . . . . . . . 92 4.17 The kinematic velocity v k model, obtained by mono-parameter v p inversion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Inversion results of Strategy C, updated v p (a) and τ (b) models after 30 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.19 Misfit curves for Strategy A (solid red), Strategy B (solid green), Strategy C (dashed blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 True and initial models: true v p (a), true τ (b), initial v p (c) and initial τ (d). In (a), the stars denote the source positions and the reversed triangles represent the receiver positions; the dashed white lines denote the interfaces where the true v p rapidly changes. . . . . . . . . . . . . . . . 4.21 Inversion results of Strategy A. . . . . . . . . . . . . . . . . . . . . . . 4.22 Inversion results of Strategy B. . . . . . . . . . . . . . . . . . . . . . . 4.23 Inversion results of the hybrid strategy. . . . . . . . . . . . . . . . . . . 4.24 Misfit curves for Strategy A (solid red line), Strategy B (solid green line), the hybrid strategy (dashed black line). . . . . . . . . . . . . . . . 4.25 v p and τ results after 80 iterations of Strategy B (a, b) and the hybrid strategy (c, d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.26 Misfit curves for Strategy A (solid red), Strategy B (solid green), and the hybrid strategy (dashed black) with 120 iterations. . . . . . . . . . . . . 4.27 One observed shot gather (direct waves in water are removed) (a) and one shot gather of initial data (b) at position x = 400 m. . . . . . . . . . 4.28 Data comparison for Strategy A at position x = 400 m, with final calculated data after 80 iterations (a) and final data difference (b), two traces of observed data (dashed blue and dark) and calculated data (solid green and red) are compared at near offset (c) and far offset (d). . . . . . . . . 4.29 Same as Figure 4.28, but for data comparison for Strategy B. . . . . . . 4.30 Same as Figure 4.28, but for data comparison for the hybrid strategy. . . 4.31 Misfit curves for Strategy A (solid red), Strategy B (solid green), the strategy of sequential inversion (dashed blue) and the hybrid strategy (dashed black) with 120 iterations. . . . . . . . . . . . . . . . . . . . . 4.32 True velocity (a) and τ (b) models, initial velocity (c) and τ (d) models. 4.33 The misfit curves of different scaling factors, the scaling 1 uses the term ∂τ (v p (x))/∂v p (x) as a dynamic scaling factor, and the scaling 2 uses the spatial mean value of ∂τ (v p (x))/∂v p (x), where x is the spatial coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.34 The scaling factor obtained by the spatial mean value of ∂τ (v p (x))/∂v p (x) for each iteration, where x is the spatial coordinates. . . . . . . . . . . 4.35 Velocity gradient comparison between velocity gradient ∂J (v p , τ )/∂v p obtained by finite difference (red) and the adjoint state (blue) methods, the left panel are the gradients without attenuation, the right panel are the gradients with Q = 50. . . . . . . . . . . . . . . . . . . . . . . . . . 4.36 ∂τ /∂v p for fixed v k of 2000 m/s (red), 2250 m/s (green) and 2500 m/s (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.1 v p and τ models: the true velocity values (a) go from 1500 m/s (water layer) to 4500 m/s (dark red on the bottom). The true τ model (b) is from 0.0011 (Q = 1000) to 0.0689 (Q = 17). τ in the water layer is set as 0.0011 (Q = 1000). The τ anomaly in light red is around 0.052 (Q = 22) and the right τ anomaly is around 0.0689 (Q = 17). The initial velocity model (c) is a smooth model. The initial τ model (d) is homogeneous with τ = 0.0011 (Q = 1000). . . . . . . . . . . . . . . . . . . . . . . . 123 5.2 Inversion results of Strategy A after 5 iterations (a) and 80 iterations (b) over v p , with τ always being fixed (c). . . . . . . . . . . . . . . . . . . 124 5.3 Inversion results of Strategy B after 5 iterations (a) and 80 iterations (c) over v p , and after 5 iterations (b) and 80 iterations (d) over τ . . . . . . . 125 5.4 Same as Figure 5.3, except for the hybrid strategy. . . . . . . . . . . . . 126 5.5 Final inversion results for (v p ,τ ) with Strategy B (a, b) and with the hybrid strategy (c, d). . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.6 One vertical section extracted at position x = 2500 m from the v p and τ models obtained with Strategy A (dashed red), compared with the initial model (dashed black) and the true model (solid blue). . . . . . . . . . . 128 5.7 Three vertical sections extracted at position x = 2500 m (left column: a, d, g, j), 3320 m (middle column: b, e, h, k) and 4000 m (right column: c, f, i, l) from inverted v p and τ . B denotes Strategy B, and H denotes the hybrid strategy. Three profiles extracted from the inverted v p and τ (dashed red) for Strategy B and for the hybrid strategy, compared with initial model (dashed black) and true model (solid blue). Inverted v p (dashed red) and final v k (dashed green) are compared at the three positions (m, n, o). The arrows in (n) denote the positions of the vertical shift between v p and v k profiles. . . . . . . . . . . . . . . . . . . . . . 130 5.8 Misfit curves for Strategy A (solid red line), Strategy B (solid green line), and the hybrid strategy (dashed blue line) up to 550 iterations (left), and the enlarged misfit curves (right). . . . . . . . . . . . . . . . . . . . . . 130 5.9 One observed shot gather (direct waves in water are removed) (a) and one shot gather of initial data (b) at position x = 3320 m. . . . . . . . . 131 5.10 Data comparison for Strategy A at position x = 3320 m, with final calculated data after 400 iterations (a) and final data difference (b), two traces of observed data (dashed blue and dark) and calculated data (solid green and red) are compared at near offset (c) and far offset (d). . . . . . . . . 132 5.11 Same as Figure 5.10, but for data comparison for Strategy B. . . . . . . 133 5.12 Same as Figure 5.10, but for data comparison for the hybrid strategy. . . 134 5.13 Same as Figure 5.10, except that the calculated data are simulated with v p = v k and τ = 0. Comparison between one observed shot gather (direct waves in water are removed) at position at x = 3320 m (a), calculated data over v k (v p = v k ) (b) and data difference (c). Comparison between observed data (dashed blue and dark) and calculated data over v k (solid green and red) at near offset (d) and far offset (e) (enlarged panel f). . . 5.14 The τ model with the artifacts removed from 2.7-3.4 km laterally and 1.56-1.82 km in depth (area inside the white circle). . . . . . . . . . . . 5.15 Comparison between the calculated shot gather for Strategy B (a), and the calculated shot gather in the τ model with one artifact removed (b), the data difference (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Following Figure 5.16, the trace comparison at near offset (d), and far offset (e), and enlarged comparison at far offset 800 m. . . . . . . . . . 5.17 The checkerboard consists of rectangular blocky zones, with fixed τ perturbations as 0.01 (Q around 100) from the depth of 1 km to the bottom part (a). τ is reconstructed for fixed v p (b). The dashed line indicates the limit of the constrained inversion zone. . . . . . . . . . . . . . . . . . . 5.18 Inversion results of Strategy B after 80 iterations over v p (a) and τ (b). . 5.19 Same as Figure 5.18, but for inversion results of the hybrid strategy. . . 5.20 Same as Figure 5.9, except for the noisy observed data. . . . . . . . . . 5.21 Same as Figure 5.11, except for the noisy observed data. . . . . . . . . 5.22 Same as Figure 5.12, except for the noisy observed data. . . . . . . . . 5.23 Same as Figure 5.8: misfit curves for Strategy B (solid green line), and the hybrid strategy (dashed blue line) up to 100 iterations (left), and the enlarged misfit curves (right). . . . . . . . . . . . . . . . . . . . . . . . 5.24 Mono-τ 1: the final τ inversion result; input v p is from the Strategy A after 100 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.25 Same as Figure 5.24, but for mono-τ 2 inversion; input v p is true v p model.148 5.26 Same as Figure 5.24, but for mono-τ 3 inversion; input v p is from Strategy B after 100 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . 5.27 Simultaneous inversion: the final v p (a) and τ (b); input v p is same as in Figure 5.24. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.28 Different frequency sources in the time (left) and frequency domains (right), 4 Hz central frequency Ricker wavelet (a), the source after 1.5-5 Hz band-pass filtering (b), and the source after 4.5-6 Hz band-pass filtering (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.29 Inversion results with Strategy B: 40 iterations (a) and 80 iterations (c) over v p , and 40 iterations (b) and 80 iterations (d) over τ . The low frequency data are used in first 40 iterations, and the high frequency data are used in the subsequent 40 iterations. . . . . . . . . . . . . . . . 5.30 Same as Figure 5.29, but for the inversion results with the hybrid strategy. xviii List of Figures 5.31 Misfit curves of Strategy B (solid green line), and the hybrid strategy (dashed blue line) up to 100 iterations (left), and the enlarged misfit curves (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 6.1 CIGs of three datasets obtained from updated velocity with 6 iterations: acoustic data (top), visco-acoustic data with Q = 50 (middle), viscoacoustic data with Q = 15 (right). . . . . . . . . . . . . . . . . . . . . . 163 6.2 IVA results of three datasets. . . . . . . . . . . . . . . . . . . . . . . . 164 6.3 The gradient of the objective function with respect to the macro velocity. 169 6.4 The image condition for non-attenuating (a, c) and attenuating media (b, d) by cross-correlating the source side forward-propagation with receiver side back-propagation wave fields (c, d) at zero-lag. For attenuating media, the forward data are attenuated by e -αL D e -αL U (b).

Figure 1

 1 Figure 1.1: Q image of subsurface superimposed to velocity model in an offshore Malaysia reservoir; the purple area indicates the gas clouds (adapted from Plessix et al. (2016))

Figure 1 . 2 :

 12 Figure 1.2: Migrated image without (left) and with compensation for attenuation (right) (adapted from Zhou et al. (2011))

  Figure 1.3: Scattering attenuation effects in highly heterogeneous media: the amplitude decays and traveltime fluctuates (Adapted from Sato et al. (2012)); the coda waves are made up of incoherent scattered waves.

Figure 1 . 4 :

 14 Figure 1.4: Intrinsic attenuation effects for one seismic trace in time domain (left) and normalized spectrum in frequency domain (right), the solid red line is the acoustic wave, and the solid blue line is the visco-acoustic wave with Q = 15.

Figure 1 . 5 :

 15 Figure 1.5: Q measurements in Pierre shale. The open circles are field experiments measured by McDonal et al. (1958), the solid line is the fitting curve between α and frequency, the dotted and dashed lines are computed by numerical modeling method from Carcione et al. (2002) (Adapted from Carcione (2011).

Figure 1 . 6 :

 16 Figure 1.6: Image panel (i.e, common image gather) in subsurface domain, a lower background velocity (a), a correct background velocity (b) and a higher background velocity (c) (Adapted from Lameloise et al. (2015).

Figure 1 . 7 :

 17 Figure 1.7: Effects of cross-talks between P-wave velocity v p and τ (approximates to the inverse of Q), J is the FWI objective function. The true v p model is homogenous v p = 2000 m/s with an low value anomaly 1800 m/s (a), and the true τ model is homogeneous τ = 0.001 with a high τ anomaly 0.07 (b). Both gradients with respect to v p (c) and τ (d) contain information on τ and v p , respectively (circles).

Figure 1 . 8 :

 18 Figure 1.8: Two normalized recorded seismic traces with velocity and attenuation being (v p = 2234 m/s, Q = infinity) (red) and (v p = 2168 m/s, Q = 50) (blue). The associated data are almost identical.

Figure 1

 1 Figure1.10: The sketch of searching direction of strategies: the kinematic curves (dashed green line) superimposed to the contour lines of the misfit distribution (gray solid) with the true values being v p = 2000 m/s and τ = 0.023 (Q = 50). The blue dots and red hollow circles are the searching direction of the simultaneous and hybrid strategies, respectively. For the hybrid strategy: v k updates (1 to 3) and the v p updates for fixed v k (3 to 5), more details in Chapter 4.

Figure 2 .

 2 Figure 2.1: (Left) Generalized Maxwell Body: MB denotes Maxwell body, formed by a spring with elastic modulus M l and a dashpot with viscosity coefficient η l . n pairs of Maxwell bodies and one spring with elastic modulus M H in parallel make up the GMB, here denoted as GMB-EK to refer to the model from paper of Emmerich and Korn (1987). (Right) Generalized Zener Body: ZB denotes the Zener body. n pairs of Zener bodies form the GZB. There are two kinds of Zener bodies. M Rl , δM l , M 1l and M 2l are elastic modulus, and η l is the viscosity coefficient (adapted from Moczo and Kristek, 2005).

  stress relaxation time, where M 1 and M 2 are the elastic modulus of the springs inside the SLS in series (Figure 2.1). It depends on the stress relaxation rate.

Figure 2 . 2 :

 22 Figure 2.2: The stress relaxation function of a SLS. The strain is unchanged, the imposed stress gradually decreases from unrelaxed M U to relaxed M R (adapted from Carcione (2011).

FrequencyFigure 2 . 3 :

 23 Figure 2.3: The optimization frequency range is from 1.5 Hz to 100 Hz for desired Q = 15, with empirically given Q pre with 15 (top), 13 (middle) and 12 (bottom). The τ method is the analytical solution with different pre-defined Q pre , the modified method is the global searching method.

FrequencyFigure 2 . 4 :FrequencyFigure 2 . 5 :

 2425 Figure 2.4: The optimization frequency range is from 1.5 Hz to 100 Hz for Q 0 = 50 with empirically given Q pre with 48.

Figure 2

 2 Figure 2.7: One snapshot of wavefields: stable wave propagation without numerical dispersion (a), unstable wave propagation (b), stable condition with numerical dispersion (c) (adapted from Thorbecke, 2015).

Figure 2

 2 Figure 2.8: Snapshot without CPML at 1.4 s (a) and 1.7 s (b).

Figure 2 .

 2 Figure 2.9: Same as Figure 2.8, but snapshot with CPML at 1.4 s (a) and 1.7 s (b).

Figure 2 .

 2 Figure 2.10: Numerical solution compares with analytical solution and at distance 200 m (left) and 800 m (right) with the spatial grid space being 5 m. Top: trace comparisons, bottom: trace difference.

AmplitudeFigure 2 .

 2 Figure 2.11: 20 Hz central Ricker wavelet in the time (left) and frequency (right) domain.

Figure 2 .

 2 Figure 2.13: Trace comparison between acoustic (red) and visco-acoustic (blue) waves in time domain (left panel) and frequency domain (right panel).

Figure 2 .

 2 Figure 2.14: Snapshot comparison for acoustic (left parts of the images) and viscoacoustic (right parts) waves at 200 ms (a), 500 ms (b), 800 ms (c) and 1000 ms (d).

  62) For example, three kinds of waves have the same v p = 2000 m/s with Q being infinite (acoustic), 50 and 15. The corresponding τ values are 0, 0.023 and 0.0923, respectively. Due to the dispersion, the kinematic velocity of three kinds of waves is different. Accordingly, v k of three waves at 28 Hz is 2000 m/s, 2065 m/s and 2250 m/s. Attenuation is stronger, the kinematic velocity is higher (Figure 2.15 left). On the other side, waves with different combinations of velocity and attenuation can have same kinematic velocity. For example, with the same kinematic velocity v k = 2250 m/s at 28 Hz, we have three kinds of waves with v p = 2250 m/s and infinite Q (acoustic), v p = 2180 m/s and Q = 50 and v p = 2000 m/s and Q = 15 (Figure 2.15 right).

Figure 2 .

 2 Figure 2.15: Left: different kinematic velocity v k at 28 Hz for three kinds of waves: (infinite Q, v p = 2000 m/s) in red, (Q = 50, v p = 2000 m/s) in black and (Q = 15, v p = 2000 m/s) in green. Right: same kinematic velocity v k = 2250 m/s at 28 Hz for three kinds of waves: (infinite Q, v p = 2250 m/s) in red, (Q = 50, v p = 2180 m/s) in black and (Q = 15, v p = 2000 m/s) in green.

Figure 2 .

 2 Figure 2.16: Traces for three different Q and velocity combinations (the same configurations as Figure 2.15 right): (infinite Q, v p = 2250 m/s) in green, (Q = 50, v p = 2180 m/s) in black and (Q = 15, v p = 2000 m/s) in red at the near offset (left) and far offset (right), in the time (top) and frequency domains (bottom).

  .1) where n and m are integer coefficients. Taking different orders of derivative and coefficients of a n and b m can obtain the constitutive relation between stress and strain, which represent different rheological models, such as: Maxwell body (MB), Kelvin body (KB), Zener body (ZB), Burgers body, Generalized Maxwell body (GMB) or Generalized Zener Body (GZB) (Figure 3.0d, e, f, g, h, i).

Figure 3

 3 Figure 3.0: Elastic / rheological models and constitutive relations in the time domain, except GZB in the frequency domain due to the complexity in the time domain, M l is the elastic modulus of spring, η l is the viscosity coefficient and G is the coefficient of spring-pot.

  strain relationship to obtain the general visco-elastic wave equation, and give specific examples with their associated rheological models and visco-ealstic wave equations.

Figure 4 . 1 :

 41 Figure 4.1: Different combinations of v p and τ , leading to the same v k = 2250 m/s at reference frequency 4, 20 and 40 Hz (left) and comparison of three shot gathers for different combinations of v p and τ at 20 Hz (right), the traveltime is same and amplitude becomes weaker as attenuation increases.

  We consider sources on the surface and receivers on the bottom of the model at 200 m depth. 61 sources are located between 200 and 440 m laterally, every 4 m, while receivers are deployed between 0 and 640 m. The exact model is homogeneous with v p = 2000 m/s and τ = 0.023 (or Q = 50)(Figures 4.2a, b). The equivalent kinematic velocity for a central frequency at 20 Hz is 2060 m/s. The initial model is also homogeneous, with v p = 2150 m/s and τ = 0 (Q is infinite)(Figures 4.2c, d).

Figure 4 . 2 :

 42 Figure 4.2: True velocity (a) and τ (b) models, initial velocity (c) and τ (d) models.

Figure 4 . 4 :

 44 Figure 4.3: Updated velocity (a) and τ (b) models after 3 iterations (Strategy A), updated velocity (c) and τ (d) models after 30 iterations (Strategy A), updated velocity (e) and τ (f) models after 3 iterations (Strategy B), updated velocity (g) and τ (h) models after 3 iterations (Strategy B), updated velocity (i) and τ (j) models after 3 iterations (Strategy C).

Figure 4 . 5 :

 45 Figure 4.5: Search direction for Strategies A, B and C, the blue line is the kinematic curve for all combinations of v p and τ equal to v k 2060 m/s, red solid line is the searching direction of Strategy B. The dots are the mean values of each iteration results.

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: Data comparison of Strategy B after 2 (a), 3 (b), 4 (c) iterations, in each subfigure, the upper left is the real data, upper middle is the calculated data, upper right is the data difference, the bottom left is one trace comparison at 201 m, the bottom middle is one trace comparison at 321 m, the bottom right the one trace difference at 201 m and 321 m respectively. The phase of data after 2 (a) and 4 iterations (c) are matched, meanwhile the phase of data after 3 iterations (b) has a shift.

Figure 4 . 9 :

 49 Figure 4.9: Same as Figure 4.8, but for gradients of v p and τ of Strategy C. The left panels (a) (c) and (e) are gradients after 1 iteration, and the right panels (b) (d) and (f) are gradient after 2 iterations. The full derivative dJ /dv p (e, f) are positive for the area under the acquisition coverage.

  The model size is 200 m in depth and 640 m laterally. 61 sources are located on the top from 200 m to 420 m. Each source has 201 receivers on the bottom. The true velocity model is 2000 m/s embedded with a blocky v p anomaly being 2060 m/s (Figure 4.10a). The true τ model is 0.023 (Q = 50) embedded by a non-attenuated block (Q is infinite) (Figure 4.10b). The initial models are homogenous, with initial v p = 2150 m/s (Figure 4.10c) and initial τ = 0 (Figure 4.10d). The source is a Ricker wavelet with central frequency being 20 Hz. The kinematic velocity is 2060 m/s at 20 Hz for the combination of v p = 2000 m/s and τ = 0.023, which means that the kinematic velocity is homogeneous 2060 m/s.

Figure 4 .

 4 13), the misfits of Strategy A and Strategy B have similar values for the first several iterations, focusing on the match of the kinematic information. The misfit of Strategy A could not decrease further without including the attenuation. Strategy B has a lower misfit with more iterations. Compared with Strategy B, Strategy C decreases faster during the first several iterations and has lower misfit up to 30 iterations.

Figure 4 .Figure 4 .Figure 4 .

 444 Figure 4.10: True v p (a) and τ (b) models, initial v p (c) and τ (d) models.

Figure 4 .

 4 Figure 4.13: Misfit curve of 3 strategies (left), and enlarge in details (right).

  3.2), we tested two simple models. In this section, we test Strategy C and the hybrid strategy on the model containing a τ Gaussian anomaly. The v p and τ models have uncorrelated structures. The model size is 220 m in depth and 600 m laterally. The geometry is same as previous test 4.3.2. 76 sources are located on the top from the left to the right edge. 301 receivers are distributed on the bottom. The source is a Ricker wavelet with central frequency being 20 Hz. The true v p is homogeneous 2000 m/s (Figure 4.14a). The true τ is 0.001 (Q is around 1000) with a Gaussian anomaly, the maximum value of τ is 0.0702 (Q is around 18) (Figure 4.14b). The initial v p is homogeneous 2000 m/s (Figure 4.14c), same as the true model. The initial τ is homogeneous 0.001 (Figure 4.14d). Same as in the previous sections, we perform 3 iterations with Strategy A (Figure 4.15a, b). The inverted v p has artifacts at the same position as the τ anomaly inside the true τ model. The τ anomaly turns into the v p artifacts due to cross-talks. More iterations with Strategy A (30 iterations in Figure 4.15c, d) lead to the more localized v p artifacts with higher values. It appears that the v p inversion tends to mach the kinematic information without updating τ . The misfit curve of Strategy A stops to decrease with 10% final data difference left (Figure 4.16 red).

  ). The final inversion results have much less v p artifacts and more obvious τ anomaly (Figure 4.18a, b) with lower misfit, compared with the results associated to Strategy B (Figure 4.19).

Figure 4 .

 4 Figure 4.14: True v p (a) and τ (b) models, initial v p (c) and τ (d) models.

Figure 4 .Figure 4 .Figure 4 .Figure 4 .

 4444 Figure 4.16: Misfit curves for Strategy A (solid red), Strategy B (solid green), the hybrid strategy (dashed blue).

Figure 4 .

 4 20c) and an initial τ = 0 (Figure 4.20d). As for the previous sections, we have inversion results of Strategy A, Strategy B and the hybrid strategy. 3 iterations with Strategy A briefly reconstruct the shallow part of v p model (Figure 4.21a, b). With 30 iterations, inverted v p contains more detailed structures, especially on the shallow part (Figure 4.21c, d). At the same time, 3 (Figure 4.22a) and 30 (Figure 4.22c) iterations over v p with Strategy B leads similar inversion results as Strategy A, but in some parts, the v p values in the layer in the depth of 50 m between 500 m-700 m laterally are lower than v p values obtained with Strategy A. The inverted τ has very lower values after 3 iterations (Figure 4.22b). With more iterations, the inverted τ has increased values (30 iterations in Figure 4.22 d), but the retrieved τ values are still much lower than the true ones. The misfit of Strategy B keeps lower than the misfit of Strategy A during the inversion (Figure 4.24).

  . The misfit curve of Strategy A stops to decrease after 50 iterations (Figure 4.26 red). We compare the results of Strategy B and the hybrid strategy after 80 iterations. The v p result of Strategy B (Figure 4.25a) is close to the inverted v p of the hybrid strategy C (Figure 4.25c). The inverted τ of Strategy B (Figure 4.25b) is still lower than the true model, but the retrieved area is almost on the shallow part and the deeper part is less reconstructed. Meanwhile, the inverted τ of the hybrid strategy (Figure 4.25b) has more close average values to the true τ . With more iterations, the hybrid strategy (Figure 4.26 black) has smaller misfit than the misfit of Strategy B (Figure 4.26 green). The two strategies both have small final misfit values. The final data residuals of the hybrid strategy (Figures 4.27 and 4.29) are smaller than Strategy B final data residuals (Figure 4.30).

Figure 4 .

 4 Figure 4.20: True and initial models: true v p (a), true τ (b), initial v p (c) and initial τ (d). In (a), the stars denote the source positions and the reversed triangles represent the receiver positions; the dashed white lines denote the interfaces where the true v p rapidly changes.

Figure 4 .

 4 Figure 4.23: Inversion results of the hybrid strategy.

Figure 4 .

 4 Figure 4.24: Misfit curves for Strategy A (solid red line), Strategy B (solid green line), the hybrid strategy (dashed black line).

Figure 4 .

 4 Figure 4.25: v p and τ results after 80 iterations of Strategy B (a, b) and the hybrid strategy (c, d).

Figure 4 .

 4 Figure 4.26: Misfit curves for Strategy A (solid red), Strategy B (solid green), and the hybrid strategy (dashed black) with 120 iterations.

Figure 4 .Figure 4 .

 44 Figure 4.27: One observed shot gather (direct waves in water are removed) (a) and one shot gather of initial data (b) at position x = 400 m.

Figure 4 .

 4 Figure 4.29: Same as Figure 4.28, but for data comparison for Strategy B.

Figure 4 .Figure 4 .

 44 Figure 4.30: Same as Figure 4.28, but for data comparison for the hybrid strategy.

Figure 4 .Figure 4 .

 44 Figure 4.32: True velocity (a) and τ (b) models, initial velocity (c) and τ (d) models.

Figure 4 .

 4 Figure 4.34: The scaling factor obtained by the spatial mean value of ∂τ (v p (x))/∂v p (x) for each iteration, where x is the spatial coordinates.

  m k+1 is updated by the previous model m k plus the model perturbation. J k denotes the d k cal derivative with respect to model and δd k represents the data difference between d k cal and d obs , where d k cal is generated by forward modeling based on model m k .

Figure 4 .Figure 4 .

 44 Figure 4.35: Velocity gradient comparison between velocity gradient ∂J (v p , τ )/∂v p obtained by finite difference (red) and the adjoint state (blue) methods, the left panel are the gradients without attenuation, the right panel are the gradients with Q = 50.

Figure 5

 5 Figure 5.1: v p and τ models: the true velocity values (a) go from 1500 m/s (water layer) to 4500 m/s (dark red on the bottom). The true τ model (b) is from 0.0011 (Q = 1000) to 0.0689 (Q = 17). τ in the water layer is set as 0.0011 (Q = 1000). The τ anomaly in light red is around 0.052 (Q = 22) and the right τ anomaly is around 0.0689 (Q = 17). The initial velocity model (c) is a smooth model. The initial τ model (d) is homogeneous with τ = 0.0011 (Q = 1000).

Figure 5 . 2 :Figure 5 . 3 :Figure 5 . 4 :Figure 5 . 5 :Figure 5 . 6 :

 5253545556 Figure 5.2: Inversion results of Strategy A after 5 iterations (a) and 80 iterations (b) over v p , with τ always being fixed (c).

Figure 5 . 8 :

 58 Figure 5.7: Three vertical sections extracted at position x = 2500 m (left column: a, d, g, j), 3320 m (middle column: b, e, h, k) and 4000 m (right column: c, f, i, l) from inverted v p and τ . B denotes Strategy B, and H denotes the hybrid strategy. Three profiles extracted from the inverted v p and τ (dashed red) for Strategy B and for the hybrid strategy, compared with initial model (dashed black) and true model (solid blue). Inverted v p (dashed red) and final v k (dashed green) are compared at the three positions (m, n, o). The arrows in (n) denote the positions of the vertical shift between v p and v k profiles.

Figure 5 . 9 :Figure 5 .

 595 Figure 5.9: One observed shot gather (direct waves in water are removed) (a) and one shot gather of initial data (b) at position x = 3320 m.

Figure 5 .Figure 5 .Figure 5 .

 555 Figure 5.11: Same as Figure 5.10, but for data comparison for Strategy B.

Figure 5 .

 5 Figure 5.14: The τ model with the artifacts removed from 2.7-3.4 km laterally and 1.56-1.82 km in depth (area inside the white circle).

Figure 5 .Figure 5 .Figure 5 .Figure 5 .Figure 5 .Figure 5 . 21 :Figure 5 .Figure 5 .

 5555552155 Figure 5.15: Comparison between the calculated shot gather for Strategy B (a), and the calculated shot gather in the τ model with one artifact removed (b), the data difference (c).

Figure 5 .Figure 5 .Figure 5 .Figure 5 .

 5555 Figure 5.24: Mono-τ 1: the final τ inversion result; input v p is from the Strategy A after 100 iterations.

Figure 5 .Figure 5 .Figure 5 .

 555 Figure 5.28: Different frequency sources in the time (left) and frequency domains (right), 4 Hz central frequency Ricker wavelet (a), the source after 1.5-5 Hz band-pass filtering (b), and the source after 4.5-6 Hz band-pass filtering (c).

  Contents 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Hybrid strategy -A guided inversion . . . . . . . . . . . . . 6.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Extension to 3D . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Application to real data . . . . . . . . . . . . . . . . . . . . 6.2.3 Image domain macro velocity building with Q model . . . . 6.2.4 Inverting velocity perturbation with Q model . . . . . . . . 157

  Figure 6.1: CIGs of three datasets obtained from updated velocity with 6 iterations: acoustic data (top), visco-acoustic data with Q = 50 (middle), visco-acoustic data with Q = 15 (right).

  Figure 6.2: IVA results of three datasets.

  propagation wave fields are the adjoint state variables of the forward-propagation wavefields, obtained by the adjoint state equation ∂z ) = -δd(x g , ω; x s ) iωu x + κ(Lτ + 1) ∂q ∂x + l -q = 0, l =∈ [1, L].

-

  (x g , ω; x s ) > -< λ 6 , iωu x + κ(Lτ + 1) < λ 7 , iωu z + κ(Lτ + 1) ∇ • v(x s , x -h, ω)q(x r , x + h, ω) (x s , x -h, ω)s l (x r , x + h,ωcalled the state equation, which has to be satisfied by Equations 6.2, 6.3, 6.4. Equation6.7 is called adjoint state equation, which gives as∂L(m, U , 1)∇ • v(x -2h)λ 9 (x -h) ∂L(m, U , λ) (x -2h)λ 9 (x -h) ∂L(m, U , λ) ∂I = h 2 I -λ 9 (x, h) (x -2h)λ 9 (x -h), l =∈ [1, L],(6.18) where, λ 5 , λ 6 , λ 7 and λ l 8 are written as p, v x ,v z , r l p respectively.

  Finally, the gradient of objective function becomes ∂L(m, U , λ) U (m), m) is forward functions, defined as the group of Equations 6.2, 6.3, 6.4. As m is defined as v 0 , the macro velocity is updated by∂J (v 0 )∂v 0 τ = -2ρv 0 (Lτ + 1)∇ • vλ 1 +

Figure 6

 6 Figure 6.3: The gradient of the objective function with respect to the macro velocity.

2 σl∇ 2 p

 22 = F (x s , ω),

+Figure 6 . 4 :

 64 Figure 6.4: The image condition for non-attenuating (a, c) and attenuating media (b, d) by cross-correlating the source side forward-propagation with receiver side backpropagation wave fields (c, d) at zero-lag. For attenuating media, the forward data are attenuated by e -αL D e -αL U (b). To obtain the same quality of seismic images, the source side and receiver side wave fields should be compensate the energy by e αL D and e αL U , respectively (d) (Adapted from Zhu et al. (2014)).

  To obtain the same quality of seismic images, the source side and receiver side wave fields should be compensate the energy by e αL D and e αL 'atténuation sismique est un paramètre physique utile pour détecter la présence de fluides et de gaz en profondeur. L'introduction de l'atténuation peut améliorer la qualité des images migrées et résoudre des problèmes d'amplitudes faibles sous le gaz. Avec l'information de la vitesse, c'est un facteur important pour caractériser la sub-surface.

U , respectively (d) (Adapted from Zhu et al. (2014)). . . . . . . . . . . . . 171 Chapter 1 Introduction Contents 1.1 Seismic attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Attenuation phenomenon and impacts . . . . . . . . . . . . 1.1.2 Intrinsic and scattering attenuation . . . . . . . . . . . . . . 1.1.3 Different Q models . . . . . . . . . . . . . . . . . . . . . . 1.2 Estimation of attenuation . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Ray-based methods . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Wave equation methods . . . . . . . . . . . . . . . . . . . . 1.2.3 Data domain and image domain methods . . . . . . . . . . . 1.3 Motivation, strategy and thesis organization . . . . . . . . . . . . . 1.3.1 Motivation: Inversion strategies to reconstruct velocity and attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Strategy: Kinematic relationship to guide the inversion . . . 1.3.3 Thesis organization . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . 1 Résumé du chapitre 1 L

Table 1 .

 1 1: Typical Q values of common rocks from[START_REF] Thorne | Modern Global Seismology[END_REF].

	Rock types Shale Sandstone Granite Peridotite Midmantle Lowermantle
	Q p	30	58	250	650	360	1200
	Q s	10	31	70-250 280	200	520
	1.1.3 Different Q models

  , through updating v k for fixed v p . We update the kinematic velocity and continue to retrieve v p and Q, alternately. Finally, we compare the results with conventional inversion strategies.

	trieve velocity v p and Q. This chapter first starts with the review of current re-
	search status on visco-acoustic or visco-elastic FWI. Main issues include cross-
	talks between parameters and different gradient sensitivities. To deal with these
	problems, people investigate different inversion strategies, the role of Hessian ma-
	trix and pre-conditioning operators. In this chapter, I incorporate the kinematic
	relationship into the inversion to obtain a hybrid inversion strategy (indicated in
	Section 1.3.2) based on conjugate gradient method, without the estimation of in-
	verse of Hessian matrix and pre-conditioning operators. I test this strategy on
	4 different models and compare with mono-velocity inversion and simultaneous
	inversion strategies. The key factor of this hybrid strategy is to obtain the kine-
	matic velocity. First, I apply this strategy in a homogeneous model and in another
	blocky anomaly model with known kinematic velocity. Then I test this strategy on
	a Gaussian anomaly model. Finally, I test this strategy on 2D synthetic Middle-
	East model, mainly dominated by transmitted waves. I develop the way to retrieve
	kinematic velocity v k • In Chapter 5, I apply the hybrid inversion strategy to Marmousi model. The Mar-
	mousi model is different from the models in Chapter 4. The dominated waves
	in Marmousi model are reflections.
	• In Chapter 4, I develop a hybrid inversion strategy for visco-acoustic FWI to re-

  La visco-élasticité est un élément essentiel pour l'imagerie quantitative. Les équations reposent principalement sur un modèle de Q constant. Avec les modèles rhéologiques,[START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] a proposé les "Generalized Maxwell Bodies" (GMB) avec une implémentation dans le domaine temporel.Carcione et al. (1988a) a utilisé les GZB également dans le domaine temps.[START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF] a montré que le module complexe visco-élastique était équivalent dans les deux approches. Cependant, du point de vue rhéologique, ce formalisme ne peut pas prendre en compte les équations d'onde visco-élastiques fractionnelles avec un Q constant[START_REF] Kjartansson | Constant Q-wave propagation and attenuation[END_REF].[START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity[END_REF] le premier indiqua que le modèle de Q constant est basé sur un modèle de Scott-Blair. La relation entre la contrainte et la déformation est à mi-chemin entre un ressort et un patin. Ainsi, nous proposons ici un regard différent et plus large sur le lien entre équations d'onde et modèles rhéologiques.Ce chapitre rassemble des approches existantes dans différents papiers et différents champs d'applications que sont le calcul fractionnel, la rhéologie, la mécanique et la sismologie. Avec cette approche d'unification, des liens plus clairs sont mis en avant.

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Complex modulus and visco-elastic wave equations . . . . . . . . 3.3 Example I: Nearly constant Q and associated wave equations . . . 3.4 Example II: Kjartansson constant Q and associated wave equations 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Résumé du chapitre 3

Table 3 .

 3 1: A review of visco-elasticity, ⇓ &⇐ lead to the contents in the arrow direction.

	In seismology		Viscoelasticiy
		ZB, MB, GZB, GMB ...
	Existing equations	(Nearly constant Q) Emmerich and Korn (1987) Carcione et al. (1988a)	Constant Q Kjartansson (1979) ⇓	⇐ Power law Bagley (1989) ⇓
		⇓
	Differential equation	Ordinary order Moczo and Kristek (2005) Carcione (2011) ⇓	Fractional order Mainardi (2010) ⇓

  L'inversion des formes d'ondes (FWI) est une méthode avancée et potentiellement plus puissante que les approches qui utilisent des rais. FWI exploite le champ d'onde complet pour estimer la vitesse et l'atténuation. Cependant, il est difficile de retrouver indépendamment les deux paramètres. Beaucoup de pistes ont été proposées, en particulier l'inversion séquentielle, l'inversion conjointe et l'inversion hiérarchique. À cause de la dispersion, plusieurs combinaisons de paramètres conduisent à la même vitesse cinématique. Nous proposons une inversion hybride. Dans un premier temps, l'inversion classique est appliquée pour un facteur de qualité fixe. Ensuite, une version modifiée de la FWI est menée. Dans ce cas, la vitesse cinématique est maintenue fixe. Ceci conduit à un gradient différent. Si la vitesse cinématique doit être remise à jour, plusieurs cycles sont répétés. La FWI hybride est plus rapide que les stratégies conventionnelles. Les validations sont réalisées sur exemples synthétiques (exemple gaussien et du Moyen-Orient) avec chacun ses caractéristiques. L'inversion hybride mitige les interférences entre vitesse et atténuation.

  . We call this simultaneous inversion as Strategy B. The gradients of the objective function are

	∂J (v p , τ ) ∂v p	τ	(4.11)
	and
	∂J (v p , τ )
	∂τ

vp , (4.12)

Table 4 .

 4 2: Inversion strategies for 4 different models.

	Tests	Models			Strategies
	Test1: 4.3.1	Homogeneous	A	B	C
	Test2: 4.3.2	Blocky	A	B	C
	Test3: 4.3.3	Gaussian	A	B	C	A → C	A
	Test4: 4.3.4	Middle-East	A	B		A → C	A

  The gradient of objection function(Equations 4.11, 4.12) is obtained by cross-correlating the forward propagation wave fields with the back-propagation wave fields. obs -d cal (v p + ∆v p ) 2 -1 2 d obs -d cal (v p -∆v p ) 2 2∆v p

				∆vp→0 = lim	1 2 d
	∂L(m,u,λ) ∂vp ∂L(m,u,λ) ∂τ	= ∂J (vp,τ ) ∂vp = ∂J (vp,τ ) ∂τ	+ 2ρv p (Lτ + 1)( ∂vx ∂x + ∂vz ∂z )λ p + 2ρv p + (Lκ + 1)( ∂vx ∂x + ∂vz ∂z )λ p + L l=1 1 τ l σ κ( ∂vx L l=1 ∂x + ∂vz τ τ l σ ( ∂vx ∂x + ∂vz ∂z )λ l rp , ∂z )λ l rp .	.
					(4.36)
	∂J (vp,τ ) ∂vp ∂J (vp,τ ) ∂τ	= -= -	T 0 2ρv p (Lτ + 1)( ∂vx ∂x + ∂vz ∂z )λ p + 2ρv p T 0 (Lκ + 1)( ∂vx ∂x + ∂vz ∂z )λ p + L l=1 1 τ l σ κ( ∂vx L l=1 ∂x + ∂vz τ τ l σ ( ∂vx ∂x + ∂vz ∂z )λ l rp , ∂z )λ l rp .	.
					(4.37)
	To verify the accuracy of the gradient, we compare the gradient with a finite differ-
	ence method, given by
	∂J (v p , τ ) ∂v p	= lim ∆vp→0	J (v p + ∆v p ) -J (v p -∆v p ) 2∆v p	(4.38)

321(b) Data comparison after 3 iterations

Acknowledgments

Abstract

In this chapter, we apply and compare three strategies, including mono-velocity inversion, simultaneous inversion and the hybrid inversion strategy described in Chapter 4, on the Marmousi model. During the discussion, Gaussian noise is added to the data to test the robustness of the hybrid strategy. We evaluate the mono-parameter τ inversion to discuss cross-talks between parameters and the reason of artifacts on τ . Finally, frequency continuation inversion is applied when higher frequencies are progressively included in the inversion scheme.