Accurately Stable Q-Compensated Reverse-Time Migration Scheme for Heterogeneous Viscoelastic Media

Abstract

The development of multi-component seismic acquisition technology creates new possibilities for the high-precision imaging of complex media. Compared to the scalar acoustic wave equation, the elastic wave equation takes the information of P-waves, S-waves, and converted waves into account simultaneously, enabling accurate description of actual seismic propagation. However, inherent attenuation is one of the important factors that restricts multi-component high-precision migration imaging. Its influence is mainly reflected in the following three ways: first, the attenuation of the amplitude energy makes the deep structure display unclear; second, phase distortion introduces errors to the positioning of underground structures; and third, the loss of high frequency components reduces imaging resolution. Therefore, it is crucial to fully consider the absorption and attenuation characteristics of the real Earth during seismic modeling and imaging. This paper aims to develop an accurate attenuation compensation reverse-time migration scheme for complex heterogeneous viscoelastic media. We first utilize a novel viscoelastic wave equation with decoupled fractional Laplacians to depict the Earth’s attenuation behavior. Then, an adaptive stable attenuation compensation operator is developed to realize high-precision attenuation compensation imaging. Several synthetic and field data analyses verify the effectiveness of the proposed method.

1. Introduction

Multi-component seismic acquisition and processing technology has come to the foreground in recent years. Compared with traditional acoustic scalar seismic data, multi-component data contain richer wave field information, which enables detailed description of oil and gas reservoirs [1,2,3,4,5,6,7]. To improve the imaging accuracy of multi-component seismic data, considerable attention has been paid to data preprocessing [8], wave field separation [9,10], velocity model estimation [11,12], imaging conditions [13,14] and high-performance computing [15]. In addition to the above factors, the inherent attenuation of the Earth is also a major challenge faced by multi-component imaging [16,17,18]. Its impact is mainly reflected in the following three effects: first, the attenuation of amplitude energy makes the deep structure display unclear; second, phase distortion introduces errors to the positioning of underground structures; and third, the loss of high-frequency components reduces the imaging resolution. Therefore, it is important to fully consider absorption and attenuation when undertaking seismic modeling and migration imaging.

Both laboratory experiments and field rock measurements have revealed that the real Earth shows almost frequency-independent attenuation behavior within the seismic frequency band [19]. The most common approach to considering attenuation in the time domain is based on mechanical (rheological) models, in which the frequency-independent attenuation is approximated by stacking mechanical units, e.g., the standard linear solid or Kelvin–Voigt elements [20]. However, such methods may not be suitable for large-scale Q-compensation reverse-time migration (Q-RTM) due to their large number of characterization parameters and high computational cost [21,22,23]. Alternatively, a viscous wave equation with decoupled fractional Laplacians (DFL) has unique advantages in many respects. The amplitude term and phase term in the DFL viscous equation are expressed separately, and accurate Q-RTM results can be obtained by reversing the sign of the amplitude term [24,25,26,27]. Currently, DFL viscous wave equations are widely used in Q-compensation seismic imaging [28,29,30,31] and inversion [32,33,34,35].

In the context of Q-ERTM, one of the fundamental essentials is to reverse-amplify the absorbed energy to ameliorate imaging accuracy and the precision of deep structures. Zhu and Sun [16] were the first to implement QERTM based on the DFL viscoelastic wave equation. Guo and McMechan [18] proposed compensating for Q effects by adjoint-based least-squares reverse-time migration. However, the high-frequency noise contained in the seismic records is also exponentially amplified in this process, which can easily cause numerical instability. Although the low-pass filtering method ensures numerical stability [16,36,37], its high-frequency cutoff threshold is fixed, so cannot flexibly adapt to medium changes and it is prone to excessive damage to high-frequency components. In recent years, some optimized stability compensation strategies have been proposed, e.g., regularization-operator-based schemes [38], stabilization strategies by improving imaging conditions [39,40], and approaches based on iterative algorithms [18,41,42]. Among them, Wang et al. [43,44,45] proposed an adaptive stability compensation operator to maximize the high-frequency components while ensuring numerical stability by which the cutoff threshold is automatically determined according to the propagation time and the Q values of the medium. However, this adaptive stability term contains hybrid-domain operators (a function of space and wavenumber). For heterogeneous media, the average Q method (the fractional power terms are computed using an average Q instead of the real Q) is usually used to solve the hybrid-domain operator problem but this method actually has the effect of blurring Q, resulting in insufficient compensation accuracy.

In our previous study [46], we proposed a novel viscoelastic wave equation with constant-order DFLs. This equation decouples the phase distortion and the amplitude attenuation effects, which facilitate the Q-RTM. More importantly, the hybrid-domain issue that plagued the original DFL equations is naturally avoided in our DFL viscoelastic wave equation, in which the order of the fractional Laplacian is fixed. Here, we use it as a forward engine, and further derive a mode-dependent adaptive stabilization operator. Although a large number of studies have contributed to Q-ERTM applying different concepts, some methods are still limited by compensation accuracy [16], computational efficiency [18] or model complexity [43]. In contrast, the proposed method achieves a more effective tradeoff between numerical stability and imaging accuracy. In contrast to the original adaptive stability approach [43], the power term in the proposed stability operator is spatial-independent, thus it meets the imaging requirements of complex heterogeneous media in an efficient manner.

The paper is arranged as follows: We start by reviewing the constant-order DFL viscoelastic wave equation [46]. Next, we present a vector-separation-based P- and S-wave decouple method. We then describe the derivation of the proposed adaptive stabilization operator in detail and embed it into the framework of Q-compensation viscoelastic reverse time migration (Q-ERTM). In a section on numerical examples, two synthetic examples and a field data application are employed to demonstrate the accuracy and the advantages of the proposed Q-ERTM scheme. Finally, we provide a summary and conclusions.

2. Theory and Methods

2.1. Viscoelastic Wave Equation with Constant-Order DFLs

Starting from an approximated dispersion relation of constant Q model [26], a novel viscoelastic wave equation with constant-order fractional Laplacians [46] is established, which is written as:

(1)$\left\{\begin{array}{l}\rho \frac{\partial v_x}{\partial t}=\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xz}}{\partial z}\\ \rho \frac{\partial v_z}{\partial t}=\frac{\partial {\sigma }_{xz}}{\partial x}+\frac{\partial {\sigma }_{zz}}{\partial z}\\ \frac{\partial {\sigma }_{xx}}{\partial t}=H_P\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z}\right)-2H_S\frac{\partial v_z}{\partial z}\\ \frac{\partial {\sigma }_{zz}}{\partial t}=H_P\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z}\right)-2H_S\frac{\partial v_x}{\partial x}\\ \frac{\partial {\sigma }_{xz}}{\partial t}=H_S\left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)\end{array}\right.$

where ${\sigma }_{xx}$, ${\sigma }_{zz}$ and ${\sigma }_{xz}$ denote the stress components, $v_x$ and $v_z$ represent the particle velocity, and

(2)$\begin{array}{c}{H}_{\alpha }=\rho {\gamma }_{\alpha }{c}_{\alpha }{\omega }_0{({\nabla }^2)}^{-\frac{1}{2}}+\rho c_{\alpha }^2-\rho {\gamma }_{\alpha }\frac{c_{\alpha }^3}{{\omega }_0}{({\nabla }^2)}^{\frac{1}{2}}\\ -\rho \pi {\gamma }_{\alpha }{c}_{\alpha }{({\nabla }^2)}^{-\frac{1}{2}}\frac{\partial }{\partial t}+\rho \pi {\gamma }_{\alpha }^2{c}_{\alpha }^2\frac{1}{{\omega }_0}\frac{\partial }{\partial t}\end{array}$

where $\alpha =P$ or $s$ corresponds to P- or S-wave, and attenuation strength ${\gamma }_{\alpha }=\arctan Q_{\alpha }{}^{-1}/\pi$. $c_{\alpha }$ refers to the propagation velocity that is defined as $c_{\alpha }=c_{0\alpha }\cos \left(\pi {\gamma }_{\alpha }/2\right)$, where $c_{0\alpha }$ represents the phase velocity at frequency ${\omega }_0$. It is worth mentioning that, in Equation (2), two terms including $\partial /\partial t$ control the amplitude decay, while the rest are phase-dominated terms. Unlike the original, for the DFL viscoelastic wave equation that contains hybrid-domain Laplacian operators, as shown in Equation (2), the wavenumber k has constant power (i.e., does not vary across space). Thus, taking it as the forward engine in Q-ERTM, this can guarantee simulation accuracy when handling heterogeneous Q media. The derivation and advantages of this equation can be found in Wang et al. [46].

2.2. Separation of P- and S- Wavefields

Multi-component seismic data simultaneously contain P-wave, S-wave and converted-wave components. If the coupled wavefields are directly introduced into the migration process for imaging, serious crosstalk noise between different wave components will be generated. Thus, it is essential to separate the P- and S-wavefields during seismic extrapolation. The early separation scheme is based on Helmholtz theory, in which the divergence and curl operators are applied to separate the coupled wavefields [47,48]. However, this scheme will introduce amplitude distortion and phase-reversal to the separation results, causing the final migration results to lose physical meaning. To overcome this phenomenon, a vector-separation-based P- and S-wave decouple method is proposed to ensure that the actual vector characteristics remain unchanged [49,50]. Similar to Wang et al. [50], we carry out wave field separation by introducing an intermediate variable ${\sigma }^p$, which satisfies:

(3)$\frac{\partial {\sigma }^p}{\partial t}=H_p\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z}\right)$

It has the same physical meaning as the scalar acoustic pressure in the viscoacoustic wave equation. From the above equation, we can obtain a series of equations to perform P- and S-wave vector-decomposition:

(4)$\rho \frac{\partial v_x^p}{\partial t}=\frac{\partial {\sigma }^p}{\partial x}$

(5)$\rho \frac{\partial v_z^p}{\partial t}=\frac{\partial {\sigma }^p}{\partial z},$

(6)$v_x^s=v_x-v_x^p,$

(7)$v_z^s=v_z-v_z^p.$

The superscripts p and s represent the P- and S-components of stress or particle velocity, respectively. Using Equations (3)–(7), we realize the independent propagation of P- and S- waves during wavefield extrapolation.

2.3. Adaptive Stable Q-Compensation Scheme

As stated in previous studies, reversing the amplitude operator while keeping the dispersion operator unchanged produces a correct compensated wavefield when Q-ERTM is implemented using the DFL viscous wave equation [28]. Using a similar operation, Equation (2) becomes:

(8)$\begin{array}{c}{H}_{\alpha }^C=-\rho {\gamma }_{\alpha }{c}_{\alpha }{\omega }_0{(-{\nabla }^2)}^{-\frac{1}{2}}+\rho c_{\alpha }^2+\rho {\gamma }_{\alpha }\frac{c_{\alpha }^3}{{\omega }_0}{(-{\nabla }^2)}^{\frac{1}{2}}\\ -\rho \pi {\gamma }_{\alpha }{c}_{\alpha }{(-{\nabla }^2)}^{-\frac{1}{2}}\frac{\partial }{\partial t}-\rho \pi {\gamma }_{\alpha }^2{c}_{\alpha }^2\frac{1}{{\omega }_0}\frac{\partial }{\partial t}.\end{array}$

Equations (2) and (8) restore the correct amplitude and phase. The last major problem faced by Q-ERTM is numerical stability. Since seismic records usually contain high-frequency noise, they will also be exponentially amplified in the process of energy recovery, and eventually lead to numerical instability. Compared to the conventional filter-based stability scheme [28], an adaptive stabilization operator [43,44,45] has the potential to gain a better balance between numerical stability and imaging resolution. According to Wang et al. [43], a mode-dependent adaptive stabilization has the following generalized form:

(9)$S_{\alpha }\left(k,l\Delta t\right)=\left\{\begin{array}{c}\frac{1}{1+{\sigma }^2{e}^{2{\varsigma }_{\alpha }\left(k\right)\Delta t}},l=1\\ \frac{1+{\sigma }^2{e}^{2{\varsigma }_{\alpha }\left(k\right)\left(l-1\right)\Delta t}}{1+{\sigma }^2{e}^{2{\varsigma }_{\alpha }\left(k\right)l\Delta t}},l=2,3,\cdots ,n\end{array}\right.$

where $e^{-{\varsigma }_{\alpha }\left(k\right)t}$ represents the amplitude-attenuated operator. Combined with Equation (2), the amplitude-attenuated operator can be given by:

(10)$e^{-{\varsigma }_{\alpha }\left(k\right)t}=\exp \left[-\frac{1}{2}\left(\pi {\gamma }_{\alpha }{c}_{\alpha }k+\pi {\gamma }_{\alpha }{}^2{c}_{\alpha }{}^2{\omega }_0{}^{-1}{k}^2\right)t\right]$

In Equation (9), ${\sigma }^2$ is the introduced stabilization factor [45], which can be empirically expressed as:

(11)${\sigma }^2=e^{-0.23G_{\lim }-1.39}$

where $G_{\lim }$ is a man-made upper limit of compensation. The proposed scheme does not need to blur Q to avoid the hybrid-domain operator because the power terms of k in Equation (10) are constant. Thus, it has higher potential to enable “point-to-point” compensation for the attenuation over the existing scheme.

2.4. Implementation of Adaptive Stable Q-ERTM

The following process is used to achieve an adaptive stable Q-ERTM:

• (a). Forward propagating the source-wavefield.

With a wavelet estimated from seismic data (or a given wavelet), solving Equations (2)–(7), we obtain the attenuated pure P- and S- forward wavefield.

• (b). Backward propagating the receiver-wavefield.

With the reversed attenuated multi-components recorded data as the boundary condition to solve Equations (3)–(8), we obtain the compensated receiver wavefield. Note that the adaptive stabilization operator (Equation (9)) is implemented at each time-step iteration to maintain numerical stability.

• (c). Applying the imaging condition.

The Q-compensated image of the subsurface structure can be produced using the following source-normalized cross-correlation imaging conditions:

(12)$\begin{array}{c}{I}_{PP}\left(x\right)={{\displaystyle \int }}_0^T\frac{v_p^S\left(x,t\right)\cdot v_p^R\left(x,t\right)}{v_p^S\left(x,t\right)\cdot v_p^S\left(x,t\right)}dt,\\ I_{ps}\left(x\right)={{\displaystyle \int }}_0^T\frac{v_p^S\left(x,t\right)\cdot v_s^R\left(x,t\right)}{v_p^S\left(x,t\right)\cdot v_p^S\left(x,t\right)}dt.\end{array}$

where the superscripts S and R correspond to the attenuated and compensated wavefield variables; and the subscripts P and S represent the P- and S-wave components, respectively.

3. Numerical Examples

3.1. Sample Sag Model

First, we conduct the proposed attenuation compensation algorithm for a simple sag model to verify its effectiveness. Figure 1 displays the P-wave velocity and quality factor $Q_p$ model, respectively. We use the conversion formula $v_s=v_p/1.73$, and $Q_s=Q_p/2$ to obtain the S-wave velocity and $Q_s$. The model contains 400 × 360 cells with a grid interval of 10 m. The simulation time is 1.5 s with a time step of 1.0 ms. There are 40 shots distributed evenly at a depth of 10 m. Each shot is recorded by 200 double-sided geophones. The source used is a Ricker wavelet with a dominant frequency of 15 Hz. Figure 2 shows the synthetic shot gathers corresponding to the elastic (Figure 2a) and viscoelastic (Figure 2b) media. As expected, the deep structural energy in the viscoelastic medium is obviously weakened due to the existence of attenuation. Figure 3a–c display the coupled wavefield, the decoupled P−component wavefield and the decoupled S−component wavefield of the 20th shot, respectively. We observe that the different components are separated accurately.

Figure 4 and Figure 5 show the migrated PP− and PS-images, where we consider three scenarios: conventional ERTM from lossless seismic data (Figure 4a and Figure 5a), attenuated seismic data without compensation (Figure 4b and Figure 5b), and attenuated seismic data with the proposed Q-ERTM scheme (Figure 4c and Figure 5c). Note that there is a wedge-shaped high attenuation area in Figure 1b, and the imaging quality below this area is significantly worse (see Figure 4b and Figure 5b), which confirms the necessity of attenuation compensation. The Q-ERTM images match the reference well and the clearer structure and better horizontal continuity verify its effectiveness.

3.2. Marmousi Model

We further employ a Marmousi model to show the compensation effect of the proposed method for complex heterogeneous media. Figure 6 displays the P−wave velocity and Q model. We also obtain the velocity and Q of S-wave through the conversion formula. This Marmousi model size is 234 × 640 with a spatial sampling interval of 10 m. There are 80 shots distributed evenly on the depth of 10 m. Each shot point has 160 bilateral records. The source we use is a Ricker wavelet with a dominant frequency of 30 Hz. The total calculation time is 2 s with a time-step of 1 ms. In this example, two different stabilization schemes are implemented: the low-pass-filter-based method and the proposed adaptive stabilization approach.

We show the migrated PP− and PS-images in Figure 7 and Figure 8 to demonstrate the compensation effect of the proposed method, where panels a–d correspond to the images using ERTM from lossless seismic data (reference solution), non-compensated images, compensated images with the low-pass filter, and compensated images with the proposed adaptively stabilized scheme, respectively. It is observed that the non-compensated migration (Figure 7b and Figure 8b) results exhibit obviously decreased resolution, weak energy of deep structure, and blurred structural details. For the low-pass-filter stabilization scheme, we only retain the components below 50 Hz to ensure numerical stability. Although the imaging quality is better than the non-compensated results, the resolution reduction caused by the loss of high-frequency information is obvious. The imaging results of the proposed scheme (Figure 7d and Figure 8d) look very close to the reference solution (Figure 7a and Figure 8a), which verifies the accuracy of compensation. It is worth mentioning that, due to the lower S-wave velocity, we observed that the PS images resolution is better than that of the PP images, which, to some extent, reflects the advantages of S-wave exploration. Figure 9 compares the traces extracted from the images generated by different migration algorithms at x = 4 km. It can be observed that the traditional low-pass filtering method (green line) is not sufficiently accurate to describe the deep structure, and both the amplitude and phase deviate from the reference trace (black line). In contrast, the proposed method (red line) can compensate for both P−wave and S-wave attenuation more accurately.

3.3. Apply to Field Data

The above synthetic data fully verify the compensation effectiveness of the proposed method. In the last example, we demonstrate its applicability to field data. Figure 10a,b show the velocity and Q model of P−wave, respectively. We also obtain the velocity and Q of S-wave through the conversion formula. The model size is 10 km × 6 km with a spatial sampling interval of 20 m. There are 250 shots distributed laterally at the depth of 20 m. Before implementing migration, we preprocess the initial data with de-noise and surface wave suppression. The source wavelet we use is estimated from the seismic records. Figure 11a displays the 200th seismogram. As represented in the averaged amplitude spectrum shown in Figure 11b, the dominant frequency is almost 25 Hz.

We show the PP− and PS-images of non-compensated migration and Q-ERTM obtained by the proposed scheme in Figure 12. For the PP−images, we find poor lateral continuity and ambiguous deep-structure in the non-compensated migration (Figure 12a); while the imaging quality in the Q-ERTM profile is improved (see Figure 12b). As indicated by the arrows, the lateral continuity of the structure is significantly enhanced, and the deep-structure becomes clear. Speaking frankly, the PS−images are definitely more challenging due to the absence of accurate S−wave velocity and Q models. Although the PS−images are inferior to the PP-images, we still observe the improvement in the imaging quality of the proposed method (Figure 12d) over the traditional ERTM (Figure 12c). Figure 13 shows a comparison of the wavenumber spectra (selected from Figure 12 at the depth of 3 km); clearly, the proposed adaptive stable Q-ERTM scheme is able to retain more high-frequency information over the non-compensated ERTM migration.

4. Conclusions

We developed an accurate attenuation compensation RTM scheme for complex heterogeneous viscoelastic media. We utilized a novel viscoelastic wave equation with constant-order fractional Laplacians as the forward engine to ensure simulation accuracy in complex attenuation media. Then, we derived a mode-dependent adaptive stabilization operator, in which the power term of the stability operator is spatial-independent, thus meeting the imaging requirements of complex heterogeneous media in an efficient manner. Two synthetic data analyses fully verify the compensation effectiveness of the proposed method, while a field data application further demonstrates its fidelity and stability. These numerical examples suggest that this adaptive, stable Q-ERTM scheme facilitates improvement in the description of complex structures in viscoelastic media.

Author Contributions

Conceptualization, N.W. and Y.S.; formal analysis, N.W. and H.Z.; methodology, N.W.; software, N.W.; supervision, Y.S. and H.Z.; writing—original draft, N.W.; writing—review and editing, H.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data used in this study is available by contacting the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Sag model. (a) P−wave velocity model and (b) P−wave quality factor [Forumla omitted. See PDF.] model.

Figure 2. Synthetic shot gathers corresponding to (a) the elastic media and (b) the loss media. Note that every 4th shot is shown.

Figure 3. The 20th snapshots. (a) the coupled wavefield, (b) the decoupled P−component wavefield and (c) the decoupled S−component wavefield.

Figure 4. PP−images of the Sag model obtained by (a) the conventional ERTM from lossless seismic data, (b) attenuated seismic data without compensation, (c) attenuated seismic data with the proposed Q-ERTM scheme.

Figure 5. PS-images of the Sag model obtained by (a) the conventional ERTM from lossless seismic data, (b) attenuated seismic data without compensation, (c) attenuated seismic data with the proposed Q-ERTM scheme.

Figure 6. The Marmousi (a) P−wave velocity and (b) quality factor [Forumla omitted. See PDF.] model.

Figure 7. PP−images generated by (a) ERTM from lossless seismic data (reference solution), (b) non-compensated images, (c) Q-compensated with low-pass filter, (d) Q-compensated with the proposed a stable scheme.

Figure 8. PS-images generated by (a) ERTM from lossless seismic data (reference solution), (b) non-compensated images, (c) Q-compensated with low-pass filter, (d) Q-compensated with the proposed a stable scheme.

Figure 9. Comparison of the traces at x = 4 km. (a) PP−component (b) PS-component.

Figure 10. The (a) P-wave velocity and (b) quality factor [Forumla omitted. See PDF.] models correspond to the field data.

Figure 11. (a) The observed seismogram selected from the 200th shot, (b) the corresponding amplitude spectrum.

Figure 12. PP-images of the field data obtained by (a) ERTM method (without Q-compensation), (b) the proposed Q-ERTM scheme; PS−images of the field data obtained by (c) ERTM method (without Q-compensation), (d) the proposed Q-ERTM scheme.

Figure 13. Comparison of wavenumber spectra that are selected from Figure 12 at the depth of 3 km, (a) PP-component, (b) PS−component.

References

1. Du, Q.; Zhu, Y.; Ba, J. Polarity reversal correction for elastic reverse time migration. Geophysics; 2012; 77, pp. S31-S41. [DOI: https://dx.doi.org/10.1190/geo2011-0348.1]

2. Duan, Y.; Sava, P. Scalar imaging condition for elastic reverse time migration. Geophysics; 2015; 80, pp. S127-S136. [DOI: https://dx.doi.org/10.1190/geo2014-0453.1]

3. Nguyen, B.D.; McMechan, G.A. Five ways to avoid storing source wavefield snapshots in 2D elastic prestack reverse time migration. Geophysics; 2015; 80, pp. S1-S18. [DOI: https://dx.doi.org/10.1190/geo2014-0014.1]

4. Luo, J.; Wang, B.; Wu, R.S.; Gao, J. Elastic full waveform inversion with angle decomposition and wavefield decoupling. IEEE Trans. Geosci. Remote Sens.; 2020; 59, pp. 871-883. [DOI: https://dx.doi.org/10.1109/TGRS.2020.2994959]

5. Zhang, Y.; Liu, Y.; Xu, S. Viscoelastic Wave Simulation with High Temporal Accuracy Using Frequency-Dependent Complex Velocity. Surv. Geophys.; 2021; 42, pp. 97-132. [DOI: https://dx.doi.org/10.1007/s10712-020-09607-3]

6. Zhang, W.; Gao, J.; Cheng, Y.; Su, C.; Liang, H.; Zhu, J. 3-D Image-Domain Least-Squares Reverse Time Migration with L1 Norm Constraint and Total Variation Regularization. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5918714. [DOI: https://dx.doi.org/10.1109/TGRS.2022.3196428]

7. Huang, X.; Greenhalgh, S.; Han, L.; Liu, X. Generalized Effective Biot Theory and Seismic Wave Propagation in Anisotropic, Poroviscoelastic Media. J. Geophys. Res. Solid Earth; 2022; 127, e2021JB023590. [DOI: https://dx.doi.org/10.1029/2021JB023590]

8. Wang, C.; Wang, Y.; Wang, X.; Xun, C. Multicomponent seismic noise attenuation with multivariate order statistic filters. J. Appl. Geophys.; 2016; 133, pp. 70-81. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2016.07.023]

9. Xiao, X.; Leaney, W.S. Local vertical seismic profiling (VSP) elastic reverse-time migration and migration resolution: Salt-flank imaging with transmitted P-to-S waves. Geophysics; 2010; 75, pp. S35-S49. [DOI: https://dx.doi.org/10.1190/1.3309460]

10. Li, Z.; Ma, X.; Fu, C.; Liang, G. Wavefield separation and polarity reversal correction in elastic reverse time migration. J. Appl. Geophys.; 2016; 127, pp. 56-67. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2016.02.012]

11. Fang, J.; Chen, H.; Zhou, H.; Zhang, Q.; Wang, L. Three-dimensional elastic full-waveform inversion using temporal fourth-order finite-difference approximation. IEEE Geosci. Remote. Sens. Lett.; 2021; 19, pp. 1-5. [DOI: https://dx.doi.org/10.1109/LGRS.2021.3125310]

12. Zhang, W.; Gao, J. Deep-learning full-waveform inversion using seismic migration images. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 5901818. [DOI: https://dx.doi.org/10.1109/TGRS.2022.3195858]

13. Du, Q.; Guo, C.; Zhao, Q.; Gong, X.; Wang, C.; Li, X.Y. Vector-based elastic reverse time migration based on scalar imaging condition. Geophysics; 2017; 82, pp. S111-S127. [DOI: https://dx.doi.org/10.1190/geo2016-0146.1]

14. Yan, R.; Xie, X.B. An angle-domain imaging condition for elastic reverse time migration and its application to angle gather extraction. Geophysics; 2012; 77, pp. S105-S115. [DOI: https://dx.doi.org/10.1190/geo2011-0455.1]

15. Zhang, Q.; Mao, W.; Chen, Y. Attenuating crosstalk noise of simultaneous-source least-squares reverse time migration with GPU-based excitation amplitude imaging condition. IEEE Trans. Geosci. Remote Sens.; 2018; 57, pp. 587-597. [DOI: https://dx.doi.org/10.1109/TGRS.2018.2858850]

16. Zhu, T.; Sun, J. Viscoelastic reverse time migration with attenuation compensation. Geophysics; 2017; 82, pp. S61-S73. [DOI: https://dx.doi.org/10.1190/geo2016-0239.1]

17. Yang, J.; Zhu, H. Viscoacoustic reverse time migration using a time-domain complex-valued wave equation viscoacoustic RTM. Geophysics; 2018; 83, pp. S505-S519. [DOI: https://dx.doi.org/10.1190/geo2018-0050.1]

18. Guo, P.; McMechan, G.A. Compensating Q effects in viscoelastic media by adjoint-based least-squares reverse time migration. Geophysics; 2018; 83, pp. S151-S172. [DOI: https://dx.doi.org/10.1190/geo2017-0235.1]

19. McDonal, F.; Angona, F.; Mills, R.; Sengbush, R.; Van Nostrand, R.; White, J. Attenuation of shear and compressional waves in Pierre shale. Geophysics; 1958; 23, pp. 421-439. [DOI: https://dx.doi.org/10.1190/1.1438489]

20. Liu, H.P.; Anderson, D.L.; Kanamori, H. Velocity dispersion due to anelasticity; implications for seismology and mantle composition. Geophys. J. Int.; 1976; 47, pp. 41-58. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1976.tb01261.x]

21. Robertsson, J.O.A.; Blanch, J.O.; Symes, W.W. Viscoelastic finite-difference modeling. Geophysics; 1994; 59, pp. 1444-1456. [DOI: https://dx.doi.org/10.1190/1.1443701]

22. Savage, B.; Komatitsch, D.; Tromp, J. Effects of 3D attenuation on seismic wave amplitude and phase measurements. Bull. Seismol. Soc. Am.; 2010; 100, pp. 1241-1251. [DOI: https://dx.doi.org/10.1785/0120090263]

23. Zhu, T.Y.; Carcione, J.M.; Harris, J.M. Approximating constant-Q seismic propagation in the time domain. Geophys. Prospect.; 2013; 61, pp. 931-940. [DOI: https://dx.doi.org/10.1111/1365-2478.12044]

24. Wang, N.; Zhou, H.; Chen, H.M.; Xia, M.M.; Wang, S.C.; Fang, J.W.; Sun, P.Y. A constant fractional-order viscoelastic wave equation and its numerical simulation scheme. Geophysics; 2018; 83, pp. T39-T48. [DOI: https://dx.doi.org/10.1190/geo2016-0609.1]

25. Wang, N.; Zhu, T.; Zhou, H.; Chen, H.; Zhao, X.; Tian, Y. Fractional Laplacians viscoacoustic wavefield modeling with k-space-based time-stepping error compensating scheme. Geophysics; 2020; 85, pp. T1-T13. [DOI: https://dx.doi.org/10.1190/geo2019-0151.1]

26. Xing, G.; Zhu, T. Modeling frequency-independent Q viscoacoustic wave propagation in heterogeneous media. J. Geophys. Res. Solid Earth; 2019; 124, pp. 11568-11584. [DOI: https://dx.doi.org/10.1029/2019JB017985]

27. Zhu, T.Y.; Bai, T. Efficient modeling of wave propagation in a vertical transversely isotropic attenuative medium based on fractional Laplacian. Geophysics; 2019; 84, pp. T121-T131. [DOI: https://dx.doi.org/10.1190/geo2018-0538.1]

28. Zhu, T.Y.; Harris, J.M.; Biondi, B. Q-compensated reverse time migration. Geophysics; 2014; 79, pp. S77-S87. [DOI: https://dx.doi.org/10.1190/geo2013-0344.1]

29. Li, Q.Q.; Zhou, H.; Zhang, Q.C.; Chen, H.M.; Sheng, S.B. Efficient reverse time migration based on fractional Laplacian viscoacoustic wave equation. Geophys. J. Int.; 2016; 204, pp. 488-504. [DOI: https://dx.doi.org/10.1093/gji/ggv456]

30. Wang, Y.; Harris, J.M.; Bai, M.; Saad, O.M.; Yang, L.; Chen, Y. An explicit stabilization scheme for Q-compensated reverse time migration. Geophysics; 2022; 87, pp. F25-F40. [DOI: https://dx.doi.org/10.1190/geo2021-0134.1]

31. Mu, X.; Huang, J.; Li, Z.; Liu, Y.; Su, L.; Liu, J. Attenuation Compensation and Anisotropy Correction in Reverse Time Migration for Attenuating Tilted Transversely Isotropic Media. Surv. Geophys.; 2022; 43, pp. 737-773. [DOI: https://dx.doi.org/10.1007/s10712-022-09707-2]

32. Xue, Z.; Sun, J.; Fomel, S.; Zhu, T. Accelerating full-waveform inversion with attenuation compensation. Geophysics; 2018; 83, pp. A13-A20. [DOI: https://dx.doi.org/10.1190/geo2017-0469.1]

33. Chen, H.; Zhou, H.; Rao, Y. Source Wavefield Reconstruction in Fractional Laplacian Viscoacoustic Wave Equation-Based Full Waveform Inversion. IEEE Trans. Geosci. Remote Sens.; 2021; 59, pp. 6496-6509. [DOI: https://dx.doi.org/10.1109/TGRS.2020.3029630]

34. Yang, J.; Zhu, H.; Li, X.; Ren, L.; Zhang, S. Estimating p wave velocity and attenuation structures using full waveform inversion based on a time domain complex-valued viscoacoustic wave equation: The method. J. Geophys. Res. Solid Earth; 2020; 125, e2019JB019129. [DOI: https://dx.doi.org/10.1029/2019JB019129]

35. Xing, G.; Zhu, T. Decoupled Fréchet kernels based on a fractional viscoacoustic wave equation. Geophysics; 2022; 87, pp. T61-T70. [DOI: https://dx.doi.org/10.1190/geo2021-0248.1]

36. Ammari, H.; Bretin, E.; Garnier, J.; Wahab, A. Time-reversal algorithms in viscoelastic media. Eur. J. Appl. Math.; 2013; 24, pp. 565-600. [DOI: https://dx.doi.org/10.1017/S0956792513000107]

37. Wang, Y. A stable and efficient approach of inverse Q filtering. Geophysics; 2002; 67, pp. 657-663. [DOI: https://dx.doi.org/10.1190/1.1468627]

38. Qu, Y.; Huang, J.; Li, Z.; Guan, Z.; Li, J. Attenuation compensation in anisotropic least-squares reverse time migration. Geophysics; 2017; 82, pp. S411-S423. [DOI: https://dx.doi.org/10.1190/geo2016-0677.1]

39. Zhao, X.; Zhou, H.; Wang, Y.; Chen, H.; Zhou, Z.; Sun, P.; Zhang, J. A stable approach for Q-compensated viscoelastic reverse time migration using excitation amplitude imaging condition. Geophysics; 2018; 83, pp. S459-S476. [DOI: https://dx.doi.org/10.1190/geo2018-0222.1]

40. Sun, J.Z.; Zhu, T.Y.; Fomel, S. Viscoacoustic modeling and imaging using low-rank approximation. Geophysics; 2015; 80, pp. A103-A108. [DOI: https://dx.doi.org/10.1190/geo2015-0083.1]

41. Zhang, W.; Gao, J. Attenuation compensation for wavefield-separation-based least-squares reverse time migration in viscoelastic media. Geophys. Prospect.; 2022; 70, pp. 280-317. [DOI: https://dx.doi.org/10.1111/1365-2478.13161]

42. Mu, X.; Huang, J.; Yang, J.; Guo, X.; Guo, Y. Least-squares reverse time migration in TTI media using a pure qP-wave equation. Geophysics; 2020; 85, pp. S199-S216. [DOI: https://dx.doi.org/10.1190/geo2019-0320.1]

43. Wang, Y.; Zhou, H.; Zhao, X.; Zhang, Q.; Chen, Y. Q-compensated viscoelastic reverse time migration using mode-dependent adaptive stabilization scheme. Geophysics; 2019; 84, pp. S301-S315. [DOI: https://dx.doi.org/10.1190/geo2018-0423.1]

44. Wang, Y.; Zhou, H.; Zhao, X.; Zhang, Q.; Zhao, P.; Yu, X.; Chen, Y. CuQ-RTM: A CUDA-based code package for stable and efficient Q-compensated reverse time migration CUDA-based Q-RTM. Geophysics; 2019; 84, pp. F1-F15. [DOI: https://dx.doi.org/10.1190/geo2017-0624.1]

45. Wang, Y.F.; Zhou, H.; Chen, H.M.; Chen, Y.K. Adaptive stabilization for Q-compensated reverse time migration. Geophysics; 2017; 83, pp. S15-S32. [DOI: https://dx.doi.org/10.1190/geo2017-0244.1]

46. Wang, N.; Xing, G.; Zhu, T.; Zhou, H.; Shi, Y. Propagating seismic waves in VTI attenuating media using fractional viscoelastic wave equation. J. Geophys. Res. Solid Earth; 2022; 27, e2021JB023280. [DOI: https://dx.doi.org/10.1029/2021JB023280]

47. Dellinger, J.; Etgen, J. Wave-field separation in two-dimensional anisotropic media. Geophysics; 1990; 55, pp. 914-919. [DOI: https://dx.doi.org/10.1190/1.1442906]

48. Yan, J.; Sava, P. Isotropic angle-domain elastic reverse-time migration. Geophysics; 2008; 73, pp. S229-S239. [DOI: https://dx.doi.org/10.1190/1.2981241]

49. Zhang, Q.; McMechan, G.A. 2D and 3D elastic wavefield vector decomposition in the wavenumber domain for VTI media. Geophysics; 2010; 75, pp. D13-D26. [DOI: https://dx.doi.org/10.1190/1.3431045]

50. Wang, W.; McMechan, G.A.; Zhang, Q. Comparison of two algorithms for isotropic elastic P and S vector decomposition. Geophysics; 2015; 80, pp. T147-T160. [DOI: https://dx.doi.org/10.1190/geo2014-0563.1]