Q-compensated least-squares reverse time migration using low-rank one-step wave extrapolation

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 asQ-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) proposed a viscoacoustic wave equation involving a pseudo-differential operator based on the constant-Q model (Kjartansson, 1979) with decoupled effects of amplitude loss and velocity dispersion. Suh et al. (2012) extended the operator to vertically transversely isotropic (VTI) media. Bai et al. (2013) adopted a similar approach for attenuation compensation in RTM, but used a viscoacoustic wave equation without memory variables. Using fractional Laplacians,  proposed a constant-Q viscoacoustic wave equation with Lowrank one-step Q-LSRTM separate terms accounting for amplitude loss and velocity dispersion, which was further applied for Q-compensated RTM using both synthetic and field data Zhu and Harris, 2015). Fletcher et al. (2012); Sun and Zhu (2015) investigated stable approaches for Q-compensation in RTM.
The imaging problem can also be cast as an inverse problem, with the objective of minimizing the L 2 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;Zhang et al., 2013;Liu 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 have been done by Ribodetti et al. (1995Ribodetti et al. ( , 2000 using an asymptotic theory and by Symes (1994, 1995) using the wave equation. Recently, Dutta and Schuster (2014) and Sun et al. (2014b) have shown that LSRTM can be applied for attenuation compensation in viscoacoustic media. Dutta and Schuster (2014) used the standard linear solid (SLS) 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) employed the lowrank one-step method to solve the constant-Q wave equation, which allows for an efficient formulation involving fractional Laplacians (Carcione, 2010;. 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 de-migration (RTDM), the forward modeling operator, and provides accurate kinematic information of subsurface structures (Symes, 2008). In viscoacoustic media, however, RTM is a poor approximation to the inverse of RTDM, because the wave amplitude suffers from attenuation during both forward and backward propagation 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.
In this paper, to construct LSRTM in viscoacoustic media, we use the lowrank one-step wave extrapolation (Sun et al., 2016) and derive its adjoint operator based on non-stationary linear filtering theory (Margrave, 1998). Sun et al. (2015) have successfully applied the lowrank 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., Lowrank one-step Q-LSRTM 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. Since the inverted matrix is numerically non-Hermitian, we adopt the Generalized Minimum Residual (GMRES) algorithm, a Krylov subspace method (Saad and Schultz, 1986), for iterative inversion. Using a synthetic model, we test the ability of the proposed Q-LSRTM to dramatically enhancing image quality at a reasonable cost.

Wave extrapolation in viscoacoustic media
We first briefly review the basic theory of lowrank one-step wave extrapolation in viscoacoustic media, and derive its adjoint operator for applications to RTM and LSRTM.
A constant-Q model (Kjartansson, 1979) describes an attenuating medium whose quality factor Q is constant in frequency (but may vary in space), indicating that the attenuation coefficient is linear in frequency.  derived the following approximate constant-Q wave equation with decoupled fractional Laplacians for modeling and imaging in viscoacoustic media: where Here γ is a dimensionless parameter that ranges between 0 to 1/2. P (x, t) is the pressure wavefield, c 0 (x) is the acoustic velocity model defined at a reference frequency ω 0 . The β 1 and β 2 parameters act like on/off switches that control velocity dispersion and amplitude loss effects, respectively ). For simplicity of notation, in the rest of the paper the fractional Laplacian operators are denoted as L = (−∇ 2 ) γ+1 and H = (−∇ 2 ) γ+1/2 .
Setting both β 1 and β 2 to one, equation 1 leads to the viscoacoustic dispersion relation with fractional powers of the wave number: Lowrank one-step Q-LSRTM Solving for ω in equation 6 yields: where: The phase function φ(x, k, ∆t) that determines the phase shift of the wavefield for propagation in time is then defined as The one-step wave extrapolation provides an approximate solution to equation 1 by incorporating the phase function defined in equation 10 into the Fourier integral operator (FIO): whereP is the spatial Fourier transform of P . The accuracy of the approximation increases with smaller ∆t . The adjoint of operator in equation 11 can be expressed asP whereφ denotes the complex conjugate of φ.
The FIOs introduced in equations 11 and 12 can be efficiently applied using the lowrank one-step wave extrapolation Sun et al. (2015), which we also refer to as the lowrank 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 lowrank PSPI operator, as well as the derivation of its adjoint operator, the lowrank NSPS operator, is shown in the appendix. RTM and 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 record from the i-th shot d i (x r , t), where x r denotes the receiver location, viscoacoustic RTM can be carried out in the following three steps: 2. Backward propagate the receiver wavefield R i (x, t) by injecting the observed seismic record as the boundary condition R i (x r , t) = d i (x r , t) and solving 3. Apply the cross-correlation imaging condition (Claerbout, 1985): whereR denotes the complex conjugate of R (Sun and Fomel, 2013).
Reverse-time demigration (RTDM) in viscoacoustic media can be formulated as the adjoint of the RTM process: 1. Calculate the source wavefield S i (x, t) in the background velocity model in the same manner as RTM by solving equation 13.
2. Generate the receiver wavefield R i (x, t) by using the stacked image I(x) as a secondary source and solving: 3. Extract the predicted seismic record (denoted by the hat) at receiver locations Note that, in order 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 lowrank PSPI is used to solve equation 14, then lowrank NSPS (derived in the appendix) 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
LSRTM aims to minimize the misfit between the observed data and predicted data measured by the quadratic function: Since A is a linear operator, a gradient-based local optimization method, such as the Conjugate Gradient method (CG), is usually applied to iteratively update the image (Dai and Schuster, 2013;Xue et al., 2014). J(m) is minimized when m satisfies (Tarantola, 2005) The square matrix A * A is known as the wave-equation Hessian, and its condition number affects the convergence rate of LSRTM implemented as an iterative inversion (Plessix and Mulder, 2004). In acoustic media, RTM usually provides a good approximation to the inverse of RTDM, and the Hessian matrix is well-conditioned (Symes, 2008). However, in viscoacoustic media, because both RTM and RTDM operators attenuate seismic waves, the image obtained by the aforementioned algorithm suffers from twice the amplitude loss accumulated along the reflection wavepath. Therefore, differently from the pure acoustic case, viscoacoustic RTM provides a poor approximation to the inverse of viscoacoustic RTDM, which makes the Hessian matrix A * A ill-conditioned. In the presence of strong attenuation and without proper preconditioning, this could slow down the convergence rate of an iterative method like CG and, in practice, may require a prohibitively large number of iterations to achieve a satisfactory result.

Q-compensated LSRTM using the GMRES method
To compensate for attenuation in seismic images,  and Sun et al. (2015) proposed the Q-compensated RTM (Q-RTM). Q-RTM in general can be formulated as follows: 1. Forward propagate the source wavefield S i (x, t) with Q-compensation by solving: 2. Backward propagate the receiver wavefield R i (x, t) with Q-compensation by solving: with the boundary condition: R i (x r , t) = d i (x r , t).

Apply the imaging condition (equation 15).
Notice the sign reversal in front of τ in equations 20-21 in comparison with equations 13-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 was developed by Sun and Zhu (2015) based on stable division between wavefields. Both the source and receiver wavefields need to be compensated in order to accumulate compensation along the entire reflection wavepath. Since Q-RTM is capable of restoring the attenuated energy in the seismic image Sun et al., 2015), it is reasonable to expect that Q-RTM is better than viscoacoustic RTM in approximating the inverse of viscoacoustic RTDM .
We propose to replace the original viscoacoustic RTM A * with Q-RTM A * c as the backward operator. The true model defined in equation 19 can be equivalently expressed as: Additionally, an RTM image may contain low-frequency noise, which can be efficiently removed by a Laplacian filter (Zhang and Sun, 2009). We propose to cascade the Q-RTM operator with a Laplacian filter to help with the least-squares inversion and speed up the convergence rate. Correspondingly, the inverted model is expressed as where L denotes the Laplacian operator. Since the operator LA * c A is closer to an identity operator than A * A, the iterative inversion process will converge faster. Equation 23 can be viewed as the solution of the preconditioned (weighted) least-squares system that seeks to minimize which leads to the solution m = (A * P * PA) −1 A * P * P d .
Instead of looking for the preconditioner P, we simply replace A * P * P with LA * c . Note that, theoretically, the inverted matrix in equation 25 is Hermitian. The new formulation (equation 23), however, 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 generalized minimum residual algorithm, GMRES(m), which solves a leastsquares 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, produces the desired image within only a few iterations. Lowrank one-step Q-LSRTM

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   (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 m along both 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 m 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).
To measure the convergence rate, we calculate the model residual as the L 2 norm of the misfit between the model calculated at each iteration m k and the true model m * , normalized by the L 2 norm of the true model: (26) Figure 8 shows the comparison of convergence rates. With the help of a Laplacian filter, LSRTM is able to achieve a slightly faster convergence rate at early iterations than the original LSRTM. The proposed Q-LSRTM, on the other hand, converges significantly faster than the other two methods. Convergence is achieved by Q-LSRTM within approximately 50 iterations, while the other two methods have not converged even after 100 iterations. The fast convergence is an important property, because for large-scale 3D seismic imaging problems only a few iterations can be afforded in practice.

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, as well as other types of preconditioners, 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 un-attenuated events, and approximately correspond to smaller eigenvalues of the forward operator . Inversion routines based on Krylov subspace methods, such as CG and GMRES, will favor large eigenvalues, which approximately correspond to shallower and un-attenuated reflectors. This leads to the observed behavior of LSRTM without Q-compensation, which first focused on improving shallow reflectors above the attenuation zone. Blanch and Symes (1994) suggested a simple way of assigning more weights to deeper reflectors, by post-conditioning 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 wave path.

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 implicitly 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.
Since 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 is preferable in application to accurate seismic imaging in attenuating media.

ACKNOWLEDGMENTS
We thank Joakim Blanch, Aaron Stanton, and one anonymous reviewer for their suggestions that helped to improve the paper. We thank Mrinal Sen and Lexing 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 UT Austin, and the third author is supported by the Jackson School Distinguished Postdoctoral Fellowship at UT Austin. We thank Texas Advanced Computing Center (TACC) for providing computational resources used in this study.

APPENDIX A: DERIVATION OF THE LOWRANK PSPI AND NSPS OPERATORS
Let p(x, t) be the seismic wavefield at location x and time t, with the spatial Fourier transform denoted by P (k, t). The wavefield at the next time step t + ∆t can be approximated by the Fourier integral operator (Wards et al., 2008;Fomel et al., 2013): p(x, t + ∆t) = P (k, t) e i φ(x,k,∆t)+ix·k dk , where φ(x, k, ∆t) is the phase function. The mixed-domain operator in equation 27 is also referred to as the one-step wave extrapolation operator Sun and Fomel, 2013). Because the wave extrapolation matrix is complex, it propagates a complex wavefield with the imaginary part being the Hilbert transform of the real part: where P and P r , respectively, denote the complex wavefield and real wavefield, and F denotes spatial Fourier transform ).
The adjoint form of operator 27 can be written as: whereφ denotes the complex conjugate of φ. The −iφ term in equation 30 indicates stepping backward in time. Expressing the dual domain operator 30 in the space domain and stepping forward in time, we arrive at a different operator: p(x, t + ∆t) = p(y, t) e i φ(y,k,∆t)+ik·(x−y) dk dy .
The phase function in equation 29 depends on the output space x, and thus represents a kind of non-stationary combination filter (Margrave, 1998). In comparison, the phase function appearing in equation 31 depends on the input space y, and leads to a kind of non-stationary convolution filter. Both operators 29 and 31 apply the same wave-propagation phase function, φ(x, k, ∆t); the one-step low-rank wave extrapolation operator 29 applies the phase function in the wavenumber domain after forward Fourier transform, whereas the new operator 31 applies the phase function in the space domain before the inverse Fourier transform. The essential difference between the two is that the non-stationary convolution has the physical interpretation of scaled, linear superposition of the non-stationary filter impulse responses, as suggested by Huygens' principle, whereas non-stationary combination filters do not have such implications (Margrave, 1998). The mixed-domain operator 29 is an equivalent to the most accurate limiting case of phase shift plus interpolation (PSPI) method (Gazdag and Squazzero, 1984;Kesinger, 1992), which has been a popular choice for one-way wavefield extrapolators. The proposed operator 31 is analogous to the non-stationary phase shift (NSPS) method (Margrave, 1998;Margrave and Ferguson, 1999) for one-way wave extrapolation.
The low-rank algorithm introduced by Fomel et al. (2013) is a separable approximation that selects a set of N representative spatial locations and M representative wavenumbers, which correspond to rows and columns from the original wave-propagation matrix. The low-rank one-step wave extrapolation uses low-rank decomposition to approximate the mixed-domain phase function in equation 29: a mn e ix·k W (x n , k)P (k, t)dk , whose computational cost effectively equals that of applying N inverse fast Fourier transforms per time step, where N is the approximation rank and is typically a number less than ten.
With the help of low-rank decomposition, the computational effort for the new NSPS method can be made identical to that of the low-rank PSPI wave extrapolation, by approximating the wave propagation operator appearing in equation 31 with P (k, t + ∆t) ≈ a mn e −ix·k W (x, k m )p(x, t)dx .
Note that, for simplicity of notation, equations 32 and 33 do not include an operation of the forward and inverse Fourier transforms in place between p(x, t) and P (k, t).
In this appendix we have presented the adjoint operator to lowrank one-step wave extrapolation. Because the derivation of the adjoint operator is discrete, i.e., using the lowrank 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 calculated the adjoint state variable. This strategy is usually refered to as discretize then optimize (Betts, 2010).