$Q$-compensated LSRTM using the GMRES method

To compensate for attenuation in seismic images, Zhu et al. (2014) 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(\mathbf{x},t)$ with $Q$-compensation by solving:
   

   $${\frac{1}{c^2}} {\frac{\partial^2 S_i(\mathbf{x},t)}{\part... ...i(\mathbf{x},t) = \delta(\mathbf{x}-\mathbf{x}_i) f(t) \; .$$ (20)

   
   

2. Backward propagate the receiver wavefield $R_i(\mathbf{x},t)$ with $Q$-compensation by solving:
   

   $${\frac{1}{c^2}} {\frac{\partial^2 R_i(\mathbf{x},t)}{\part... ...{\partial}{\partial t} \mathbf{H} R_i(\mathbf{x},t) = 0 \;,$$ (21)

   
   
   with the boundary condition: $R_i(\mathbf{x}_r,t) = d_i(\mathbf{x}_r,t)$.
3. Apply the imaging condition (equation 15).

Notice the sign reversal in front of $\tau$ 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; Zhu et al., 2014), 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 $\mathbf{A}^*$ with $Q$-RTM $\mathbf{A}^*_c$ as the backward operator. The true model defined in equation 19 can be equivalently expressed as:

$$\mathbf{m} = (\mathbf{A}^*_c \mathbf{A})^{-1}\mathbf{A}^*_c \mathbf{d} \; .$$ (22)

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

$$\mathbf{m} = (\mathbf{L}\mathbf{A}^*_c\mathbf{A})^{-1} \mathbf{L}\mathbf{A}^*_c \mathbf{d}\;,$$ (23)

where $\mathbf{L}$ denotes the Laplacian operator. Since the operator $\mathbf{L}\mathbf{A}^*_c\mathbf{A}$ is closer to an identity operator than $\mathbf{A}^*\mathbf{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

$$J_p(\mathbf{m}) = \frac{1}{2}\Vert \mathbf{P} (\mathbf{A} \mathbf{m} - \mathbf{d}) \Vert^2_2 \; ,$$ (24)

which leads to the solution

$$\mathbf{m} = (\mathbf{A}^* \mathbf{P}^* \mathbf{P} \mathbf{A})^{-1} \mathbf{A}^* \mathbf{P}^* \mathbf{P} \mathbf{d} \;.$$ (25)

Instead of looking for the preconditioner $\mathbf{P}$, we simply replace $\mathbf{A}^* \mathbf{P}^* \mathbf{P}$ with $\mathbf{L} \mathbf{A}^*_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 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, produces the desired image within only a few iterations.