Anisotropic Attenuation Compensated Reverse Time Migration of Pure qP-Wave in Transversely Isotropic Attenuating Media

Abstract

The absorption (anelastic attenuation) and anisotropy properties of subsurface media jointly affect the seismic wave propagation and the quality of migration imaging. Anisotropic viscoelastic model can effectively describe seismic velocity and attenuation anisotropy effects. To reduce the computational cost and complexity of elastic wave modes decoupling for seismic imaging in anisotropic attenuating media, we have developed a pure-viscoacoustic transversely isotropic (TI) wave equation starting from the complex-valued velocity dispersion relation of quasi-compressional (qP) wave. The wave equation involving fractional Laplacians has advantages of being able to describe the constant-Q (frequency-independent quality factor) attenuation, arbitrary TI velocity and attenuation, decoupled amplitude loss and velocity dispersion effects. Numerical analyses showed that the simplified equation can accurately hold the velocity and attenuation anisotropy of qP-wave in viscoelastic anisotropic media in the range of moderate anisotropy. Compared to previous pseudo-viscoacoustic equations, the pure-viscoacoustic equation can be completely free from undesirable S-wave artifacts and behaves good numerical stability in tilted transversely isotropic (TTI) attenuating media. There are obvious wavefield differences between isotropic attenuation and anisotropic attenuation cases especially in the direction perpendicular to the axis of symmetry. Furthermore, to mitigate the influences of velocity and attenuation anisotropy on migrated seismic images, we have developed an anisotropic attenuation (Q) compensated reverse time migration (AQ-RTM) approach based on the new propagator. The compensation can be implemented by reversing the sign of the dissipation terms and keeping the dispersion terms unchanged during wavefields extrapolation. Synthetic example from a Graben model illustrated that the anisotropic Q-compensated RTM scheme can produce images with more balanced amplitude and accurate position of reflecters compared with conventional RTM methods under assumptions of acoustic anisotropic (uncompensated) and isotropic attenuating media. Results from a Marmousi-II model demonstrated that the new methodology is applicable for complicated geological model to significantly improve imaging resolution of the target area and deep layers.

Appendices

Appendix 1: 3D Pure-Viscoacoustic TI Wave Equation

Under acoustic approximation, the dispersion relation of qP-wave in 3D VTI attenuating media can be simplified as

$$\omega^{2} \approx \rho^{ - 1} \left[ {M_{11} \left( {k_{x}^{2} + k_{y}^{2} } \right) + M_{33} k_{z}^{2} + \frac{{\left( {M_{13} \sqrt {1 + 2\delta } - M_{11} } \right)\left( {k_{x}^{2} + k_{y}^{2} } \right)k_{z}^{2} }}{{\left( {1 + 2\varepsilon } \right)\left( {k_{x}^{2} + k_{y}^{2} } \right) + k_{z}^{2} }}} \right],$$

(27)

The corresponding time-wavenumber domain wave equation for dispersion relation is,

$$\frac{{\partial^{2} p}}{{\partial t^{2} }} = - \left[ {L_{11} \left( {k_{x}^{2} + k_{y}^{2} } \right) + L_{33} k_{z}^{2} + \left( {L_{13} \sqrt {1 + 2\delta } - L_{11} } \right)\frac{{\left( {k_{x}^{2} + k_{y}^{2} } \right)k_{z}^{2} }}{{k^{2} }}} \right]p,$$

(28)

where L_ij denote fractional Laplace operators in Eq. (10). where $k^{2} = k_{x}^{2} + k_{y}^{2} + k_{z}^{2}$, here we assume $\varepsilon = 0$ for the denominator of the fractional term. The wave equation can be conveniently solved by using the pseudospectral method as follow

$$\begin{aligned} \frac{{\partial^{2} p}}{{\partial t^{2} }} & = L_{11} \mathcal{F}^{ - 1} \left[ { - \left( {k_{x}^{2} + k_{y}^{2} } \right)\mathcal{F}\left[ p \right]} \right] + L_{33} \mathcal{F}^{ - 1} \left[ { - k_{z}^{2} \mathcal{F}\left[ p \right]} \right] \\ & \,\,\,\, + \left( {L_{13} \sqrt {1 + 2\delta } - L_{11} } \right)\mathcal{F}^{ - 1} \left[ { - \frac{{\left( {k_{x}^{2} + k_{y}^{2} } \right)k_{z}^{2} }}{{k^{2} }}\mathcal{F}\left[ p \right]} \right], \\ \end{aligned}$$

(29)

The similar dispersion relation of qP-wave for TTI attenuating media in 3D case can be deduced from Eq. (27) through coordinate transformation as

$$\begin{aligned} \rho \omega^{2} & \approx M_{11} \left( {\tilde{k}_{x}^{2} + \tilde{k}_{y}^{2} } \right) + \left( {M_{33} + M_{13} \sqrt {1 + 2\delta } - M_{11} } \right)\tilde{k}_{z}^{2} \\ + \left( {M_{11} - M_{13} \sqrt {1 + 2\delta } } \right)\frac{{\tilde{k}_{z}^{4} }}{{k^{2} }}, \\ \end{aligned}$$

(30)

where $\tilde{k}_{x} ,\tilde{k}_{y} ,\tilde{k}_{z}$ are spatial wavenumbers in the rotated coordinate system.

Substituting the following transformation formula from Eq. (13)

$$\begin{gathered} \tilde{k}_{x} = k_{x} \cos \tilde{\theta }\cos \phi + k_{y} \cos \tilde{\theta }\sin \phi + k_{z} \sin \tilde{\theta }, \hfill \\ \tilde{k}_{y} = - k_{x} \sin \phi + k_{y} \cos \phi , \hfill \\ \tilde{k}_{z} = - k_{x} \sin \tilde{\theta }\cos \phi - k_{y} \sin \tilde{\theta }\sin \phi + k_{z} \cos \tilde{\theta }, \hfill \\ \end{gathered}$$

(31)

into Eq. (30) and after some algebraic manipulations, the dispersion relation can be futher formulate as.

$\rho \omega^{2} \approx M_{11} w_{1} + \left( {M_{33} + M_{13} \sqrt {1 + 2\delta } - M_{11} } \right)w_{2} + \left( {M_{11} - M_{13} \sqrt {1 + 2\delta } } \right)w_{3} ,$, (32).where

$$\begin{aligned} w_{1} & = \left( {\cos^{2} \tilde{\theta }\cos^{2} \phi + \sin^{2} \phi } \right)k_{x}^{2} + \left( {\cos^{2} \tilde{\theta }\sin^{2} \phi + \cos^{2} \phi } \right)k_{y}^{2} + \sin^{2} \tilde{\theta }k_{z}^{2} \\ \, & \,\,\,\,\, - \sin^{2} \tilde{\theta }\sin 2\phi k_{x} k_{y} + \sin 2\tilde{\theta }\cos \phi k_{x} k_{z} + \sin 2\tilde{\theta }\sin \phi k_{y} k_{z} , \\ \\ w_{2} & = \sin^{2} \tilde{\theta }\cos^{2} \phi k_{x}^{2} + \sin^{2} \tilde{\theta }\sin^{2} \phi k_{y}^{2} + \cos^{2} \tilde{\theta }k_{z}^{2} \\ \, & \,\,\,\,\,{ + }\sin^{2} \tilde{\theta }\sin 2\phi k_{x} k_{y} - \sin 2\tilde{\theta }\cos \phi k_{x} k_{z} - \sin 2\tilde{\theta }\sin \phi k_{y} k_{z} , \\ \\ w_{3} & = \sin^{4} \tilde{\theta }\cos^{4} \phi \frac{{k_{x}^{4} }}{{k^{2} }} + \sin^{4} \tilde{\theta }\sin^{4} \phi \frac{{k_{y}^{4} }}{{k^{2} }} + \cos^{4} \tilde{\theta }\frac{{k_{z}^{4} }}{{k^{2} }} \\ \, & \,\,\,\,\, + 2\sin^{4} \tilde{\theta }\sin 2\phi \cos^{2} \phi \frac{{k_{x}^{3} k_{y} }}{{k^{2} }} - 2\sin 2\tilde{\theta }\sin^{2} \tilde{\theta }\cos^{3} \phi \frac{{k_{x}^{3} k_{z} }}{{k^{2} }} + 2\sin^{4} \tilde{\theta }\sin 2\phi \sin^{2} \phi \frac{{k_{x} k_{y}^{3} }}{{k^{2} }} \\ \, & \,\,\,\,\, - 2\sin 2\tilde{\theta }\sin^{2} \tilde{\theta }\sin^{3} \phi \frac{{k_{y}^{3} k_{z} }}{{k^{2} }} - 2\sin 2\tilde{\theta }\cos^{2} \tilde{\theta }\cos \phi \frac{{k_{x} k_{z}^{3} }}{{k^{2} }} - 2\sin 2\tilde{\theta }\cos^{2} \tilde{\theta }\sin \phi \frac{{k_{y} k_{z}^{3} }}{{k^{2} }} \\ \, & \,\,\,\,\, + \frac{3}{2}\sin^{4} \tilde{\theta }\sin^{2} 2\phi \frac{{k_{x}^{2} k_{y}^{2} }}{{k^{2} }} + \frac{3}{2}\sin^{2} 2\tilde{\theta }\cos^{2} \phi \frac{{k_{x}^{2} k_{z}^{2} }}{{k^{2} }} + \frac{3}{2}\sin^{2} 2\tilde{\theta }\sin^{2} \phi \frac{{k_{y}^{2} k_{z}^{2} }}{{k^{2} }} \\ \, & \,\,\,\,\, - 3\sin 2\tilde{\theta }\sin^{2} \tilde{\theta }\sin 2\phi \cos \phi \frac{{k_{x}^{2} k_{y} k_{z} }}{{k^{2} }} - 3\sin 2\tilde{\theta }\sin^{2} \tilde{\theta }\sin 2\phi \sin \phi \frac{{k_{x} k_{y}^{2} k_{z} }}{{k^{2} }} + \frac{3}{2}\sin 2\tilde{\theta }\sin 2\phi \frac{{k_{x} k_{y} k_{z}^{2} }}{{k^{2} }}. \\ \end{aligned}$$

(33)

The corresponding 3D pure-viscoacoustic TTI wave equation in the time-wavenumber domain is

$$\frac{{\partial^{2} p}}{{\partial t^{2} }} = \left[ {L_{11} w_{1} + \left( {L_{33} + L_{13} \sqrt {1 + 2\delta } - L_{11} } \right)w_{2} + \left( {L_{11} - L_{13} \sqrt {1 + 2\delta } } \right)w_{3} } \right]p.$$

(34)

When the dip angle $\tilde{\mu } = 0$, the Eq. (34) reduces to the 3D viscoacoustic VTI case. Similarly, the wavenumber terms and the fractional Laplacians in the equation can be implemented numerically using the pseudospectral method.

Appendix 2: Parameterization for TI Attenuating Media

The anisotropic attenuation can be described by quality factors matrix. As follows from Eq. (3), the Q matrix inherits the structure of the stiffness matrix. For the VTI media with VTI attenuation, the Q matrix has the form

$${\text{Q}} = \left[ {\begin{array}{*{20}c} {Q_{11} } & {Q_{12} } & {Q_{13} } & 0 & 0 & 0 \\ {Q_{12} } & {Q_{11} } & {Q_{13} } & 0 & 0 & 0 \\ {Q_{13} } & {Q_{13} } & {Q_{33} } & 0 & 0 & 0 \\ 0 & 0 & 0 & {Q_{44} } & 0 & 0 \\ 0 & 0 & 0 & 0 & {Q_{44} } & 0 \\ 0 & 0 & 0 & 0 & 0 & {Q_{66} } \\ \end{array} } \right],$$

(35)

where $Q_{11}$ and $Q_{33}$ are P-wave quality factors in horizontal and vertical directions, $Q_{55}$ is SV-wave quality factor and aslo responsible for the SH-wave attenuation in the symmetry (vertical) direction, $Q_{66}$ is responsible for the SH-wave attenuation in the horizontal direction.

For the convenience of description, three dimensionless attenuation anisotropy parameters following the idea of the Thomsen notation for elastic anisotropy are defined as (Zhu and Tsvankin 2006)

$$\varepsilon_{Q} = \frac{{Q_{33} - Q_{11} }}{{Q_{11} }},$$

(36)

$$\delta_{Q} = \frac{{\frac{{Q_{33} - Q_{55} }}{{Q_{55} }}C_{44} \frac{{\left( {C_{13} + C_{33} } \right)^{2} }}{{\left( {C_{33} - C_{55} } \right)^{2} }} + 2\frac{{Q_{33} - Q_{55} }}{{Q_{13} }}C_{13} \left( {C_{13} + C_{55} } \right)}}{{C_{33} \left( {C_{33} - C_{55} } \right)}},$$

(37)

$$\gamma_{Q} = \frac{{Q_{55} - Q_{66} }}{{Q_{66} }}.$$

(38)

where $\varepsilon_{Q}$ describes the difference between the P-wave attenuation in the horizontal and vertical directions, $\gamma_{Q}$ is responsible for the attenuation anisotropy of SH-waves, $\delta_{Q}$ controls the curvature of the P-wave attenuation coefficient near the vertical direction. Under the acoustic VTI approximation (i.e., $C_{44} = 0$), Eq. (37) can be simplified to

$$\delta_{Q} = \frac{{2\left( {Q_{33} - Q_{13} } \right)C_{13}^{2} }}{{Q_{13} C_{33}^{2} }}$$

(39)

when $\varepsilon_{Q} = \delta_{Q} = 0$, the P-wave attenuation is isotropic (independent of directions) for arbitrary velocity anisotropy.

Appendix 3: Analytical Phase Velocity and Quality Factor

The frequency-dependent complex velocity is the key to describe seismic wave velocity and attenuation effects in attenuating media. According to the Eq. (7), the directionally dependent complex velocity of P-wave in pure-viscoacoustic VTI approximation case can be can be written as

$$V_{{\text{P}}} (\theta ) \approx \rho^{ - 1/2} \left[ {\left( {M_{11} \sin^{2} \theta + M_{33} \cos^{2} \theta } \right) + \frac{{\left( {M_{13} \sqrt {1 + 2\delta } - M_{11} } \right)\sin^{2} \theta \cos^{2} \theta }}{{\left( {1 + 2\varepsilon } \right)\sin^{2} \theta + \cos^{2} \theta }}} \right]^{1/2} ,$$

(40)

and the analytical expressions of anisotropic phase velocity and quality factor of P-wave are given by

$$v_{{\text{P}}} \left( \theta \right) = \left[ {{\text{Re}} \left( {1/V_{{\text{P}}} } \right)} \right]^{ - 1} ,$$

(41)

$$Q_{{\text{P}}} (\theta ) = \frac{{{\text{Re}} \left( {V_{{\text{P}}}^{2} } \right)}}{{{\text{Im}} \left( {V_{{\text{P}}}^{2} } \right)}},$$

(42)

where ${\text{Re}}$ and ${\text{Im}}$ are the real and imaginary parts of a complex variable. Similarly, we can obtain the corresponding analytical expressions of P-wave phase velocity and quality factor in viscoelastic VTI and approximate pure-viscoacoustic TTI cases based on Eqs. (4) and (15), respectively.