## ABSTRACT

Attenuation of seismic waves needs to be taken into account to improve the accuracy of seismic imaging. In viscoacoustic media, reverse time migration (RTM) can be performed with $Q$-compensation, which is also known as $Q$-RTM. Least-squares RTM (LSRTM) has also been shown to be able to compensate for attenuation through linearized inversion. However, seismic attenuation may significantly slow down the convergence rate of the least-squares iterative inversion process without proper preconditioning. We have found that incorporating attenuation compensation into LSRTM can improve the speed of convergence in attenuating media, obtaining high-quality images within the first few iterations. Based on the low-rank one-step seismic modeling operator in viscoacoustic media, we have derived its adjoint operator using nonstationary filtering theory. The proposed forward and adjoint operators can be efficiently applied to propagate viscoacoustic waves and to implement attenuation compensation. Recognizing that, in viscoacoustic media, the wave-equation Hessian may become ill-conditioned, we propose to precondition LSRTM with $Q$-compensated RTM. Numerical examples showed that the preconditioned $Q$-LSRTM method has a significantly faster convergence rate than LSRTM and thus is preferable for practical applications.

## INTRODUCTION

Seismic attenuation is caused by the effective anelastic properties of the earth (Aki and Richards, 2002; Carcione, 2007) and may lead to poor illumination and misplacement of reflectors in a migration image. To directly compensate for seismic attenuation during reverse time migration (RTM) (Baysal et al., 1983; McMechan, 1983; Whitmore, 1983), Zhang et al. (2010) propose a viscoacoustic wave equation involving a pseudodifferential operator based on the constant-$Q$ model (Kjartansson, 1979) with decoupled effects of amplitude loss and velocity dispersion. Suh et al. (2012) extend the operator to vertically transversely isotropic (VTI) media. Bai et al. (2013) adopt a similar approach for attenuation compensation in RTM, but use a viscoacoustic wave equation without memory variables. Using fractional Laplacians, Zhu and Harris (2014) propose a constant-$Q$ viscoacoustic wave equation with separate terms accounting for amplitude loss and velocity dispersion, which is further applied for $Q$-compensated RTM using synthetic and field data (Zhu et al., 2014; Zhu and Harris, 2015). Fletcher et al. (2012) and Sun and Zhu (2015) investigate stable approaches for *Q*-compensation in RTM.

The imaging problem can also be cast as an inverse problem, with the objective of minimizing the $L2$ norm of the difference between recorded data and predicted data (Ronen and Liner, 2000). Such approaches are known as least-squares migration (Nemeth et al., 1999; Tang, 2009; Dai et al., 2011) and more specifically least-squares RTM (LSRTM) in the context of RTM (Wong et al., 2011; Dai et al., 2012; Dai and Schuster, 2013; Liu et al., 2013; Zhang et al., 2013; Sun et al., 2014a; Xue et al., 2014; Hou and Symes, 2015). LSRTM is capable of mitigating imaging artifacts and enhancing subsurface illumination and may have a correlative objective function to relax the amplitude matching requirement (Zhang et al., 2014). Pioneering works of linearized inversion in viscoacoustic and viscoelastic media are done by Ribodetti et al. (1995, 2000) using an asymptotic theory and by Blanch and Symes (1994, 1995) using the wave equation. Recently, Dutta and Schuster (2014) and Sun et al. (2014b) show that LSRTM can be applied for attenuation compensation in viscoacoustic media. Dutta and Schuster (2014) use the standard linear solid model (Robertsson et al., 1994; Blanch et al., 1995), with a simplified stress-strain relation and incorporated a single relaxation mechanism (Blanch and Symes, 1995). Sun et al. (2014b) use the low-rank one-step method to solve the constant-$Q$ wave equation, which allows for an efficient formulation involving fractional Laplacians (Carcione, 2010; Zhu and Harris, 2014). The computational cost of LSRTM depends on the number of iterations, which hinges on the conditioning of the wave-equation Hessian that it tries to invert. In acoustic media, RTM is an efficient approximation to the inverse of reverse time demigration (RTDM), the forward modeling operator, and it provides accurate kinematic information of subsurface structures (Symes, 2008). However, in viscoacoustic media, RTM is a poor approximation to the inverse of RTDM because the wave amplitude suffers from attenuation during forward and backward propagation (Zhu and Harris, 2014; Sun et al., 2015). As a result, the wave-equation Hessian becomes ill-conditioned and iterative LSRTM suffers from a slow convergence rate. To improve the convergence rate of LSRTM in viscoacoustic media, we propose to seek for a proper preconditioner that mitigates the effect of attenuation in the inversion.

To construct LSRTM in viscoacoustic media, we use the low-rank one-step wave extrapolation (Sun et al., 2016) and derive its adjoint operator based on nonstationary linear filtering theory (Margrave, 1998). Sun et al. (2015) successfully apply the low-rank one-step wave extrapolation operator to solve the constant-$Q$ wave equation with fractional Laplacians. To solve the problem of slow convergence of LSRTM in viscoacoustic media, we propose to construct a preconditioned formulation by replacing the viscoacoustic RTM operator, i.e., RTM based on the solution of the viscoacoustic wave equation forward and backward in time, with a better approximate inverse of the RTDM operator, i.e., the $Q$-compensated RTM or $Q$-RTM (Zhang et al., 2010; Suh et al., 2012; Bai et al., 2013; Zhu et al., 2014). $Q$-RTM involves a modeling operator with separate control over amplitude and phase, and is designed to compensate for the amplitude loss along the attenuated wavepaths. As a result, the preconditioned wave-equation Hessian is well conditioned, helping the new framework to quickly converge to the true amplitude solution within only a few iterations. Because the inverted matrix is numerically non-Hermitian, we adopt the generalized minimum residual (GMRES) algorithm (Saad and Schultz, 1986), for iterative inversion. Using a synthetic model, we test the ability of the proposed $Q$-LSRTM to dramatically enhance image quality at a reasonable cost.

## THEORY

### Wave extrapolation in viscoacoustic media

We first briefly review the basic theory of low-rank one-step wave extrapolation in viscoacoustic media.

The FIOs introduced in equations 11 and 12 can be efficiently applied using the low-rank one-step wave extrapolation (Sun et al., 2015), which we also refer to as the low-rank phase shift plus interpolation (PSPI) operator because of its resemblance to the well-known PSPI method for solving the one-way wave equation (Gazdag and Squazzero, 1984; Kesinger, 1992; Margrave and Ferguson, 1999). The detailed formulation of low-rank PSPI operator, as well as the derivation of its adjoint operator, the low-rank NSPS operator, is shown in the Appendix A. LSRTM in viscoacoustic media can therefore be constructed using the forward and adjoint operators.

### Viscoacoustic RTM and RTDM

To obtain a seismic image with an attenuated data from the $i$th shot $di(xr,t)$, where $xr$ denotes the receiver location, viscoacoustic RTM can be carried out in the following three steps:

- Forward propagate the source wavefield $Si(x,t)$ by solving$1c2\u22022Si(x,t)\u2202t2\u2212\eta LSi(x,t)\u2212\tau \u2202\u2202tHSi(x,t)=\delta (x\u2212xi)f(t).$(13)
- Backward propagate the receiver wavefield $Ri(x,t)$ by injecting the observed seismic record as the boundary condition $Ri(xr,t)=di(xr,t)$ and solving$1c2\u22022Ri(x,t)\u2202t2\u2212\eta LRi(x,t)\u2212\tau \u2202\u2202tHRi(x,t)=0.$(14)
- Apply the crosscorrelation imaging condition (Claerbout, 1985):where $R\u203e$ denotes the complex conjugate of $R$ (Sun and Fomel, 2013).$I(x)=\u2211i\u2211tSi(x,t)R\xafi(x,t),$(15)

The RTDM in viscoacoustic media can be formulated as the adjoint of the RTM process:

Calculate the source wavefield $Si(x,t)$ in the background velocity model in the same manner as RTM by solving equation 13.

- Generate the receiver wavefield $Ri(x,t)$ by using the stacked image $I(x)$ as a secondary source and solving$1c2\u22022Ri(x,t)\u2202t2\u2212\eta LRi(x,t)\u2212\tau \u2202\u2202tHRi(x,t)=I(x)Si(x,t).$(16)
- Extract the predicted seismic record (denoted by the hat) at receiver locations $xr$:$d^i(xr,t)=Ri(xr,t).$(17)

Note that, to make the RTDM adjoint to RTM, the wave extrapolation operator used to solve equation 16 needs to be the adjoint of the operator used to solve equation 14. For example, if low-rank PSPI is used to solve equation 14, then low-rank NSPS (derived in Appendix A) needs to be used to solve equation 16. If we write the RTM process symbolically as $m^=A*d$, where $m^$ is the stacked image, $A*$ is the viscoacoustic RTM operator, and $*$ denotes the adjoint, then the RTDM process corresponds to $d^=Am$, where $d^$ represents the predicted data and $A$ is the viscoacoustic RTDM operator.

### Least-squares RTM in viscoacoustic media

### $Q$-compensated LSRTM using the GMRES method

To compensate for attenuation in seismic images, Zhu et al. (2014) and Sun et al. (2015) propose the $Q$-compensated RTM ($Q$-RTM). $Q$-RTM in general can be formulated as follows: Notice the sign reversal in front of $\tau $ in equations 20 and 21 in comparison with equations 13 and 14. This reversal aims to recover the attenuated wavefield by undoing phase distortion and amplifying the amplitude. Practically, it may also exponentially increase noise through each time step. A low-pass filter can be applied to stabilize the extrapolation process. Another robust compensation strategy is developed by Sun and Zhu (2015) based on stable division between wavefields. The source and receiver wavefields need to be compensated to accumulate compensation along the entire reflection wavepath. Because $Q$-RTM is capable of restoring the attenuated energy in the seismic image (Zhu et al., 2014; Sun et al., 2015), it is reasonable to expect that $Q$-RTM is better than viscoacoustic RTM in approximating the inverse of viscoacoustic RTDM.

- Forward propagate the source wavefield $Si(x,t)$ with $Q$-compensation by solving$1c2\u22022Si(x,t)\u2202t2\u2212\eta LSi(x,t)+\tau \u2202\u2202tHSi(x,t)=\delta (x\u2212xi)f(t).$(20)
- Backward propagate the receiver wavefield $Ri(x,t)$ with $Q$-compensation by solvingwith the boundary condition $Ri(xr,t)=di(xr,t)$.$1c2\u22022Ri(x,t)\u2202t2\u2212\eta LRi(x,t)+\tau \u2202\u2202tHRi(x,t)=0,$(21)
Apply the imaging condition (equation 15).

Instead of looking for the preconditioner $P$, we simply replace $A*P*P$ with $LAc*$. Note that, theoretically, the inverted matrix in equation 25 is Hermitian. However, the new formulation (equation 23) makes the inverted matrix numerically non-Hermitian. One complication with equation 23 is that because the square matrix being inverted is no longer Hermitian, iterative methods for Hermitian positive-definite matrices are not optimal (Saad, 2003). Therefore, we implement a complex-valued restarted GMRES(m) with m iterations between restarts, which solves a least-squares system by searching for the vector in the Krylov subspace with minimum residual (Saad and Schultz, 1986). We refer to the method of solving equation 23 by GMRES(m) as $Q$-compensated LSRTM or $Q$-LSRTM. As demonstrated in the numerical examples of the next section, $Q$-LSRTM is capable of achieving a significantly faster convergence rate than conventional LSRTM and, in practice, can produce the desired image within only a few iterations.

## NUMERICAL EXAMPLES

To test the convergence rate of $Q$-LSRTM, we use a portion of the BP 2004 velocity model (Billette and Brandsberg-Dahl, 2004) and the corresponding $Q$ model suggested by Zhu et al. (2014) (Figure 1). The model features a low-velocity, low-$Q$ area, which is assumed to be caused by the presence of a gas chimney. The model has a spatial sampling rate of $12.5\u2009\u2009m$ along vertical and horizontal directions. A total of 31 shots with a spacing of 162.5 m have been modeled with attenuation, and the source is a Ricker wavelet with 22.5 Hz peak frequency (Figure 2). Performing RTM without compensating for amplitude loss, i.e., using the dispersion-only operator, leads to an image corresponding to Figure 3b, which suffers from poor illumination below the gas chimney. In contrast, $Q$-RTM appears capable of recovering the amplitude at deeper reflectors (Figure 3c), but the image still exhibits some differences from the true reflectivity. Note that the dispersion-only RTM image and $Q$-RTM image have the same phase but differ in amplitude. Next, we perform LSRTM (equation 19) and $Q$-LSRTM (equation 23) through a number of iterations. To test the separate effect of applying a Laplacian filter without compensating for attenuation, we also perform LSRTM with a Laplacian operator that removes low-frequency artifacts. For fairness of comparison, all three methods are driven by the GMRES method and because the tested model is small enough, they were not run in a restarted fashion. Using the original LSRTM (Figure 4), the inversion process attempts to remove low-frequency noise and improve the illumination of deeper reflectors. However, at 30th iteration, the reflector amplitude and sharpness beneath the attenuating area still have not been recovered. LSRTM with a Laplacian filter (Figure 5) achieves a somewhat sharper image, but because a Laplacian filter boosts high-frequency components in the image, the reflectors beneath the attenuating zone remain poorly illuminated. In contrast, the proposed $Q$-LSRTM method (Figure 6) produces sharper reflectors with well-balanced illumination, especially in the area beneath the gas chimney using the same number of iterations. Note that the color scales used in all the three cases are kept the same as that of the true model (Figure 3a). Figure 7 compares the image traces extracted at $X=2500\u2009\u2009m$ from the 30th iteration results against the true model. Clearly, the result obtained by the proposed $Q$-LSRTM best represents the true reflectivity, especially at deeper parts beneath the gas chimney (below 800 m depth).

## DISCUSSION

In this work, the GMRES method is used to invert the non-Hermitian matrices. Similar to the CG method, the full GMRES method (with no restarts) converges in no more than $n$ steps, where $n$ is the total size of the model. However, the GMRES method requires additional memory to store the previous stepping directions. In the numerical examples presented above, because the model space is small enough, no restarts are needed. However, for large 3D models, restarts might be required for a large number of iterations, which could compromise global optimality. Fortunately, the proposed method is designed to achieve a satisfying result within only a small number of iterations. In practical applications where each iteration consumes large computing resources, only a small number of iterations is usually afforded.

The goal of preconditioning LSRTM in viscoacoustic media using $Q$-RTM is to alleviate the computational burden on iterative inversion by compensating for attenuation explicitly in wave propagation. Therefore, the iterations can be spent on removing migration artifacts and compensating irregularities in subsurface illumination caused by other reasons, such as acquisition footprint. Due to attenuation, the events of the reflectors beneath the attenuating zone yield a smaller amplitude compared with unattenuated events and approximately correspond to smaller eigenvalues of the forward operator (Blanch and Symes, 1994). Inversion routines based on Krylov subspace methods, such as CG and GMRES, will favor large eigenvalues, which approximately correspond to shallower and unattenuated reflectors. This leads to the observed behavior of LSRTM without $Q$-compensation, which first focuses on improving shallow reflectors above the attenuation zone. Blanch and Symes (1994) suggest a simple way of assigning more weights to deeper reflectors, by postconditioning the seismic record with an increasing function of time. The proposed method is similar in spirit but more accurate, in that $Q$-compensation removes the true effect of attenuation in the gradient by accurately compensating for attenuation along the entire wavepath.

## CONCLUSIONS

We have introduced a novel way of preconditioning least-squares RTM to achieve a faster convergence rate in viscoacoustic media. The data-space preconditioner is defined by the $Q$-compensated RTM operator, with the goal of recovering amplitude loss due to attenuation and removing low-frequency artifacts in the gradient. Because the square matrix to be inverted becomes numerically non-Hermitian, we adopt the GMRES algorithm to perform iterative inversion. Our synthetic examples show that the proposed $Q$-LSRTM is capable of producing an accurate $Q$-compensated image within significantly fewer iterations than LSRTM, and thus it is preferable in application to accurate seismic imaging in attenuating media.

## ACKNOWLEDGMENTS

We thank J. Blanch, A. Stanton, and one anonymous reviewer for their suggestions that helped to improve the paper. We thank M. Sen and L. Ying for constructive discussions. We thank the sponsors of the Texas Consortium for Computation Seismology (TCCS) for financial support. The first author is supported additionally by the Statoil Fellows Program at the University of Texas at Austin, and the third author is supported by the Jackson School Distinguished Postdoctoral Fellowship at the University of Texas at Austin. We thank Texas Advanced Computing Center (TACC) for providing computational resources used in this study.

### DERIVATION OF THE LOW-RANK PSPI AND NSPS OPERATORS

In this appendix, we have derived the adjoint operator to low-rank one-step wave extrapolation. Because the derivation of the adjoint operator is discrete, i.e., using the low-rank approximation matrices instead of state and adjoint state equations, the result presented in the appendix has wider applications not limited to the scope of this paper. For example, in full-waveform inversion applications, where the adjoint state equation is difficult to obtain, the discrete adjoint operator described in this appendix can be efficiently applied to calculate the adjoint state variable. This strategy is usually referred to as discretize then optimize (Betts, 2010).