Theory

A constant-$Q$ model assumes that the attenuation coefficient is linear in frequency (Kjartansson, 1979). The time-domain viscoacoustic wave equation for constant-$Q$ model can be written as (Carcione et al., 2002)

$$\frac{\partial^{2-2\gamma} P}{\partial t^{2-2\gamma}} = c^2\omega_0^{-2\gamma}\nabla^2 P \;,$$ (1)

where $\mathbf{x}$ is the spatial coordinate vector, $P(\mathbf{x},t)$ is the pressure wavefield,

$$\gamma(\mathbf{x})=\arctan(1/Q(\mathbf{x}))/\pi$$ (2)

is a dimensionless parameter, $\nabla^2$ is the Laplacian operator,

$$c^2(\mathbf{x})=c_0^2(\mathbf{x})\cos^2(\pi\gamma(\mathbf{x})/2) \; ,$$ (3)

and $c_0(\mathbf{x})$ is the velocity model defined at the reference frequency $\omega_0$. When $Q$ is finite, $\gamma$ is greater than zero, and the wave equation involves a fractional time derivative.

In a homogeneous medium, considering the plane-wave solution $e^{(i\omega t - i\tilde{\mathbf{k}}\cdot\mathbf{x})}$, where $\omega$ is the angular frequency and $\tilde{\mathbf{k}}$ is the complex wavenumber vector and substituting it into equation 1, leads to the dispersion relation

$$\frac{\omega^2}{c^2} = (i)^{2\gamma} \omega_0^{-2\gamma}\omega^{2\gamma}\tilde{\mathbf{k}}^2 \;.$$ (4)

Starting from equation 4, Zhu and Harris (2014) derived the following approximate dispersion relation:

$$\frac{\omega^2}{c^2} = -\old{\beta_1 }\eta \vert\mathbf{k}\vert^{2\gamma +2} - i \old{\beta_2 }\omega\tau \vert\mathbf{k}\vert^{2\gamma +1} \;,$$ (5)

which corresponds to the following constant-$Q$ wave equation with fractional Laplacians:

$$\frac{1}{c^2}\frac{\partial^2 P}{\partial t^2} = \nabla^2 P + \ol... ...+ \old{\beta_2} \tau \frac{\partial}{\partial t}(-\nabla^2)^{\gamma+1/2} P \; ,$$ (6)

where $\eta = -c_0^{2\gamma}\omega_0^{-2\gamma}\cos(\pi \gamma)$ and $\tau = -c_0^{2\gamma-1}\omega_0^{-2\gamma}\sin(\pi \gamma)$. Note that both $c_0$ (the phase velocity) and $\gamma$ (the fractional power) can be heterogeneous (depending on $\mathbf{x}$). Solving for $\omega$ in equation 5 yields:

$$\omega = \frac{-ip_1 + p_2}{2} \; ,$$ (7)

where:

$$p_1 = \tau c^2\vert\mathbf{k}\vert^{2\gamma+1} \; ,$$ (8)

$$p_2 = \sqrt{-\tau^2c^4\vert\mathbf{k}\vert^{4\gamma+2}-4\eta c^2\vert\mathbf{k}\vert^{2\gamma+2}} \;.$$ (9)

The phase function $\phi (\mathbf{x},\mathbf{k},\Delta t)$ that determines the phase shift of wavefield for propagation in time is then defined as

$$\phi_1 (\mathbf{x},\mathbf{k},\Delta t) = \mathbf{k} \cdot \mathbf{x} + \frac{-ip_1 + p_2}{2} \Delta t$$ (10)

and can be supplied to the low-rank one-step wave extrapolation algorithm (Sun and Fomel, 2013).

To compensate for attenuation, reversing the sign of the second term on the right hand side of equation 5 can amplify the amplitude, while the other term must be kept unchanged to counteract the dispersion effects (Zhu, 2014). Thus, the attenuation-compensated constant-$Q$ wave equation corresponds to the dispersion relation:

$$\frac{\omega^2}{c^2} = -\eta \vert\mathbf{k}\vert^{2\gamma +2} + i\omega\tau \vert\mathbf{k}\vert^{2\gamma +1} \; ,$$ (11)

which defines the complex-conjugate phase function:

$$\phi_2 (\mathbf{x},\mathbf{k},\Delta t) = \bar{\phi_1} (\mathbf{x... ...},\Delta t) = \mathbf{k} \cdot \mathbf{x} + \frac{ip_1 + p_2}{2} \Delta t \; .$$ (12)

Both $\phi_1$ and $\phi_2$ involve the fractional power of wavenumber, and depend on both $\mathbf{x}$ and $\mathbf{k}$, the Fourier transform pair. The low-rank one-step wave extrapolation operator (Sun and Fomel, 2013) provides a convenient way to utilize the phase function to extrapolate a viscoacoustic wavefield, while allowing $\gamma(\mathbf{x})$, the fractional power coefficient, to vary in space. The one-step mixed-domain operator has the form of the following Fourier integral operator:

$$P(\mathbf{x},t+\Delta t) = \int \hat{P}(\mathbf{k},t) e^{i \phi(\mathbf{x},\mathbf{k},\Delta t)} d\mathbf{k}\;,$$ (13)

where $\hat{P}$ is the spatial Fourier transform of $P$. Its adjoint form can be expressed as

$$\hat{P}(\mathbf{k},t) = \int P(\mathbf{x},t+\Delta t) e^{-i \bar{\phi}(\mathbf{x},\mathbf{k},\Delta t)} d\mathbf{x}\; ,$$ (14)

where $\bar{\phi}$ denotes the complex conjugate of $\phi$. Substituting equation 10 into equation 13, we can write the forward extrapolation operator explicitly as

$$P(\mathbf{x},t+\Delta t) = \int \hat{P}(\mathbf{k},t) e^{i \phi... ... e^{i\mathbf{k} \cdot \mathbf{x} + (p_1 + ip_2) \Delta t/2} d\mathbf{k} \; .$$ (15)

The corresponding adjoint operator takes the form:

$$\hat{P}_{adj}(\mathbf{k},t) = \int P(\mathbf{x},t+\Delta t) e^{-... ...,e^{-i\mathbf{k} \cdot \mathbf{x} + (p_1 - ip_2) \Delta t/2} d\mathbf{x} \; .$$ (16)

To make the computation feasible, we apply low-rank decomposition proposed by Fomel et al. (2013) to approximate the wave extrapolation symbol in equations 15 and 16.

Note that the adjoint operator compensates for velocity dispersion. Because the sign of $p_1$ remains unchanged, the amplitude of waves will be attenuated during backward propagation using equation 16. This means that the adjoint operator may not be well suited for RTM because the backward wavefield would be attenuated twice. However, the adjoint operator can be supplied into the framework of least-squares RTM which can then iteratively recover the amplitude loss caused by viscoacoustic attenuation (Sun et al., 2014).

An alternative strategy for backward wave propagation is to try compensating for the full attenuation effect (both amplitude loss and velocity dispersion) using the operator

$$P_{comp}(\mathbf{x},t+\Delta t) = \int \hat{P}(\mathbf{k},t) e^{... ...,e^{i\mathbf{k} \cdot \mathbf{x} + (-p_1 + ip_2) \Delta t/2} d\mathbf{k} \; ,$$ (17)

which is the adjoint of

$$\hat{P}_{comp}(\mathbf{k},t) = \int P(\mathbf{x},t+\Delta t) e^{... ...e^{-i\mathbf{k} \cdot \mathbf{x} + (-p_1 - ip_2) \Delta t/2} d\mathbf{x} \; .$$ (18)

In application to RTM, operator in equation 17 has an exponentially growing term which can amplify the energy in the forward-propagated source wavefield. In the same manner, operator in equation 18 can compensate for the amplitude loss in the backward-propagated receiver wavefield. In order to avoid magnifying high-frequency components during wave propagation, we employ a low-pass Tukey filter for the attenuation and dispersion operators in the wavenumber domain (Zhu et al., 2014). A complex-valued imaging condition that compensates for both velocity dispersion and amplitude loss can be obtained by cross-correlating the source and receiver wavefields modeled by equations 17 and 18 (Zhu et al., 2014).

In this way, the migrated viscoacoustic data gets corrected for attenuation due to viscoacoustic material encountered along the wave path. As will be demonstrated in the numerical examples, $Q$-compensated RTM is capable of enhancing illumination in areas where attenuating material is present. For a more accurate compensation of attenuation, operators 17 and 18 can be used to design a preconditioner in iterative least-squares RTM.