Seismic Modeling, Imaging and Inversion in Viscoacoustic Media

During wave propagation, seismic energy is dissipated by the geometrical spreading, heterogeneity scattering and lattice internal friction. The energy decay related to internal friction is known as intrinsic attenuation, which reflects the viscosity (anelastic) property of subsurface minerals and rocks. Particularly, the saturations of gas give rise to strong seismic intrinsic attenuation. Incorporating attenuation into seismic modeling, imaging and inversion enables accurate detection of hydrocarbon reservoir and characterization of fluid properties. To date, a number of wave equations have been developed to describe the intrinsic attenuation effects. For example, in the frequency-domain, the viscous behavior can be described using a complex-valued velocity. It explicitly incorporates the quality factor (Q) into the wave equation and therefore is easy to utilize to compensate attenuation effects in reverse-time migration (RTM) and full-waveform inversion (FWI). But its requirements for solving the Helmholtz equation include large computer memory cost, which is still challenging for largescale 3D models. On the other hand, the attenuation can be incorporated in the time-domain wave equation based upon the standard linear solid (SLS) theory. Since the dispersion and dissipation are coupled in the SLS, and the quality factor Q has to be transformed to stress and strain relaxation times, it is difficult to use in seismic imaging and inversion. Another popular time-domain wave equation is formulated based on constant-Q theory. Although this approach has the advantage that the dispersion and dissipation terms are decoupled, it is necessary to calculate a mixed-domain operator using complicated numerical solvers, such as the low-rank approximation. In this study, starting from the frequency-domain viscoacoustic wave equation, I first use a second-order polynomial to approximate the dispersion term, followed by a pseudo-differential operator to approximate the dissipation term. These two approximations make it possible to transform the frequency-domain equation into the time domain, and derive a new complexvalued viscoacoustic wave equation. The advantages of the new wave equation include: (1) the dispersion and dissipation effects are naturally separated, which can be used to compensate amplitude loss in seismic migration by reversing the sign of the dissipation term; (2) Q is explicitly incorporated into the wave equation, which makes it easy to directly derive the misfit gradient with respect to Q and estimate subsurface attenuation models using Q-FWI; and (3) this new viscoacoustic wave equation can be numerically solved using finite-difference time marching and a Fourier transform, which does not require mixed-domain solvers as required in the constant-Q method, and has lower memory cost than the frequency-domain approach. Based on the new complex-valued wave equation, I develop a viscoacoustic RTM workflow to correct the attenuation-associated dispersion and dissipation effects. A time-reversed wave equation is derived to extrapolate receiver-side wavefields, in which the sign of the dissipation term is reversed while the dispersion term remains unchanged. In wavefield extrapolation, both source and receiver wavefields are complex-valued and their real and imaginary parts satisfy the Hilbert transform. This analytic property helps to explicitly decompose up- and down-going waves. Then, a causal imaging condition, which crosscorrelates the down-going source-side wavefield and the up-going receiver-side wavefield, is utilized to suppress lowwavenumber artifacts in migrated images. Furthermore, with limited recording apertures, finite-frequency source functions, irregular subsurface illuminations, viscoacoustic RTM is still insufficient to produce satisfactory reflectivity images with high resolution and amplitude fidelity. By incorporating the complexvalued wave equation into a linear waveform inversion scheme, I develop a viscoacoustic least-squares reverse-time migration (LSRTM) scheme. Based on the Born approximation, I first linearize the wave equation and derive a viscoacoustic demigration operator. Then, using the Lagrange multiplier method, I derive the adjoint viscoacoustic wave equation and the corresponding sensitivity kernels. With the forward and adjoint operators, a linear inverse problem is formulated to estimate the subsurface reflectivity model. A total-variation regularization is introduced to enhance the robustness of the proposed viscoacoustic LSRTM, and a diagonal Hessian is used as a preconditioner to accelerate convergence. Traditional waveform inversion for attenuation is commonly based on the SLS wave equation, in which case the quality factor Q has to be converted to the stress and strain relaxation times. When using multiple attenuation mechanisms in the SLS method, it is difficult to directly estimate these relaxation time parameters. Based on the new time-domain complex-valued viscoacoustic wave equation, I present an FWI framework for simultaneously estimating subsurface P-wave velocity and attenuation distributions. Since Q is explicitly incorporated into the wave equation, I directly derive sensitivity kernels for P-wave velocity and attenuation using the adjoint-state method, and simultaneously estimate their distributions. By analyzing the Gauss-Newton Hessian, I observe strong inter-parameter crosstalk artifacts, especially the leakage from velocity to Q. I approximate the Hessian inverse using a preconditioned L-BFGS method in FWI, which significantly reduces inter-parameter crosstalk, and produces accurate velocity and attenuation models.