An explicit stable Q-compensated reverse time migration scheme for complex heterogeneous attenuation media

Prestack reverse-time migration (RTM) is a popular imaging technique for complex geological conditions, since the amplitude attenuation and velocity dispersion are common in seismic recordings. To image attenuated seismic recordings accurately, a robust migration algorithm with a stable attenuation compensation approach should be considered. In the context of the Q-compensated RTM approach based on the decoupled fractional Laplacians (DFLs) viscoacoustic wave equation, amplitude compensation can be implemented by flipping the sign of the dissipation term. However, the non-physical magnification of image amplitude could lead to a well-known numerical instability problem. The explicit stabilization operator can rectify the amplitude attenuation and suppress the numerical instability. However, limited by the inconvenient mixed-domain operator, the average Q value rather than the real Q value is often used in the compensation operator, lowering the compensated accuracy of the migration image. To overcome this problem, we propose a novel explicit Q-compensation scheme. The main advantage of the proposed compensation operator is that its order is space-invariant, making it more suitable for handling complex heterogeneous attenuation media. Several two-dimensional (2D) and three-dimensional (3D) synthetic models are used to verify the superiority of the proposed approach in terms of amplitude fidelity and image resolution. Field data further demonstrates that our approach has potential applications and can greatly enhance the resolution of seismic images.


Introduction
Seismic attenuation is commonly observed in the rock matrix and pores such as sandstones and gas clouds (Aki and Richards, 1980;Dutta and Schuster, 2014), which is commonly characterized by the quality factor Q (Tselentis, 1998). The inherent attenuation properties of the Earth media significantly affect the characteristics of seismic waves in terms of amplitude, phase, and frequency band (Carcione, 1990;Zhu et al., 2013;Yang et al., 2015). As accurate attenuation compensation is critical for understanding the mechanisms of seismic data and mapping the Earth's interior, ignoring the anelasticity would result in blurred migrated images, dimmed amplitudes, and ultimately reduce the reliability of seismic interpretation (Wang and Guo, 2004;Zhang and Gao, 2021). Currently, the Q-compensated schemes are divided into two categories. The first method is usually implemented on the poststack seismic profile, such as the inverse-Q filtering method (Bickel and Natarajan, 1985;Wang, 2002;Wang, 2003), time-variant spectral whitening (Yilmaz, 2001), and time-varying deconvolution (Margrave et al., 2011). Although these methods are computationally efficient, they are usually suitable for simple structures because of the lack of consideration for the internal physical connotation in the process of wave propagation (Hargreaves and Calvert, 1991). The second approach is based on wave equations and compensates for attenuation along the propagation path (Dai and West, 1994;Mittet, 2007;Li et al., 2016b). Since the attenuation occurs during wave propagation, it is more reasonable to implement Q-compensated as part of the migration process (Zhang et al., 2010;Zhu and Sun, 2017;Zhao Y. et al., 2018).
Earth materials usually exhibit nearly frequency-independent Q behavior over the seismic frequency band (Korneev et al., 2004;Dvorkin and Mavko, 2006;Zhu, 2017). Mechanical models are the most commonly used approaches for describing frequencyindependent Q behavior (Mainardi, 2010;Rossikhin, and Shitikova, 2010). Among these, the standard linear solid (SLS) and the generalized standard linear solid (GSLS) models have been extensively used in modeling and imaging (Xu and McMechan, 1995;Causse and Ursin, 2000;Deng and McMechan, 2008;. Alternatively, some mathematicalmodel-based schemes, such as the Kolsky-Margravechen model (Kolsky, 1956;Futterman, 1962), the power-law model, and Kjartansson's constant-Q model (Kjartansson, 1979;Zhu and Bai, 2018;Wang et al., 2020), are also gradually applied to the characterization of Earth attenuation. Recently, the decoupled fractional Laplacians (DFL) viscoacoustic wave equation has attracted attention (Chen et al., 2016;Yao et al., 2016;Wang N. et al., 2018;Xing and Zhu, 2019;Liu and Luo, 2021) for the following reasons. First, it incorporates the spatial fractional Laplacians to avoid the memory issue often encountered by conventional anelastic modeling (Dutta and Schuster, 2014) and fractional time derivative approaches (Carcione, 2008). Second, compared with SLS, this strategy is more attractive for Q-RTM because it realizes amplitude compensation by flipping the sign of the amplitudeloss term without changing phase information (Guo and McMechan, 2015;Guo et al., 2016).
Attenuation compensation is critical for improving imaging quality in complex attenuation structures. However, Q-RTM usually suffers from numerical instability due to the inversephysical energy amplification of high-frequency noise. Several strategies have been proposed to improve numerical stability. For example, the schemes include regularization approaches (Zhang et al., 2010;Wang et al., 2012), filter-based approaches Wang Y. et al., 2018;Chen et al., 2020a), improved imaging conditions (Xie et al., 2015;Zhao X. H. et al., 2018;Sun and Zhu, 2018;Yang et al., 2021) and least-squares Q reverse-time migration (QLSRTM) (Chen et al., 2020b;. Among these, the implicit adaptive stabilization compensation scheme (Wang et al., 2017) that adjusts the truncation frequency according to the propagation time and Q value provides a better trade-off between numerical stability and imaging resolution. However, the implicit compensation operator is calculated in the wavenumber domain, implying that it requires additional Fourier transforms. Wang et al. (2019) further proposed an explicit compensation strategy, which adjusts compensation parameters more conveniently and simplifies the workflow of Q-RTM. Nevertheless, this strategy is not straightforward for dealing with the complex heterogeneous Q media (Sun and Zhu, 2015) because it involves the spatial variable-order Laplacians (mixeddomain operators). To solve the mixed-domain problem, the average Q value is usually used to replace the real Q value (Wang et al., 2019), which would introduce significant errors for the sharp Q-contrast models (Chen et al., 2016;Yang and Zhu, 2018a;Xing and Zhu, 2019).
In this study, we aim to derive a new explicit Q-compensated scheme so as to accurately image complex heterogeneous attenuation structures. Our derivation seeks an equation that avoids calculating spatial variableorder Laplacian operators starting from the explicit compensation wave equation (Wang et al., 2019). To accomplish this, we used Taylor expansion to approximate the spatial variable-order Laplacian operators to the spatial constant-order form and then integrate the constant-order compensated scheme into the Q-RTM framework. The following are the advantages of the proposed method. First, the proposed scheme enables us to compensate for amplitude loss in Q-RTM without changing the phase because the dispersion and dissipation effects are naturally separated. Second, the explicit Q-compensated scheme is free from frequent Fourier transforms, so it is expected to simplify the workflow of Q-RTM. Third, compared with the original compensation operator with average Q, the proposed algorithm visibly improves the imaging quality in the heterogeneous Q media.
The organization of this study is as follows. First, we review the explicit Q-compensation wave equation with spatial variable-order Laplacian operators (Wang et al., 2019;Wang et al., 2022). Second, we demonstrate the derivation of the proposed spatial constant-order Q-compensated wave equation, followed by detailing its numerical implementation and validating its accuracy. Then, much synthetic and field data are used to demonstrate its advantages. Finally, we conduct a discussion and draw some conclusions.

Materials and methods
Constant-order compensated viscoacoustic wave equation. According to a previous study (Wang et al., 2019;Wang et al., 2021), the wave equation based on the explicit compensated operators can be expressed as follows: and where c 0 is the phase velocity at the reference frequency ω 0 . p and ∇ 2 represent the pressure wavefield and Laplacian operator, respectively. η a and τ a represent phase dispersion and amplitude compensation terms, respectively. η c1 and η c2 is introduced to keep the phase unchanged during Q compensation. τ c is a stabilization term to keep the high-wavenumber component from uncontrolled amplification, in which β and σ are the stable compensation parameters. The stable compensation parameters usually depend on the media's physical properties. Generally, A larger β can reserve more higher frequency, and the more high-wavenumber signal is eliminated with σ increases (Wang et al., 2021). Equation 1 can compensate for amplitude without introducing high-frequency noise because it includes the stabilization term. However, subject to the thorny mixed-domain (spatialwavenumber domain) operators, the average Q value is often used when extrapolating the wavefields Mu et al., 2021). Even though the average scheme works well for the homogeneous or smooth Q model, it cannot accurately describe the characteristics of wave propagation in the heterogeneous media (Yang and Zhu, 2018b). To overcome this shortcoming, we propose a new explicit viscoacoustic wave equation with the constant fractional-order Laplacians using a truncated Taylor expansion algorithm (Chen et al., 2016).
If seismic waves propagate in a homogeneous medium, the plane wave equation is substituted into Eq. 1 to obtain the frequencywavenumber domain viscoacoustic wave equation wherep represents the wavefield in the wavenumber domain. Here, the fractional Laplacian c 3 is taken as an example of an approximation where f d and k d represent the dominant frequency and dominant wave number, respectively, 1 + 2γ ln k k d can be approximated by Taylor expansion again where ζ is a empirical constant coefficient (Chen et al., 2016), which is introduced to guarantee ( k kd ) ζ − 1 close to zero. Therefore, the lefthand side of Eq. 5 can be approximated as follows: The other spatial variable-order Laplacians in Eq. 4 can be approximated in the same way. Therefore, Eq. 3 can be expressed as follows: By transforming Eq. 9 back to the space domain, we obtain the following equation In Eq. 11, the fractional Laplacians is the constant-order form, and it naturally adapts to sharp Q media. Note that for Q → ∞ case, Eq. 11 reduces to the classic acoustic-wave equation. We introduce the general principle of Q-RTM in the framework of the proposed wave equation (Zhu, 2016;Wang et al., 2021).

Stable constant-order Q-compensated RTM framework
The complete procedure of the proposed Q-compensated RTM is summarized as follows.
Step 1. Forward extrapolating the source wavefield.
Through a given source wavelet S(x s , t), we solve Eq. 11 to extrapolate the source wavefield, then checkpoint wavefields (Chen et al., 2020c) are stored. For source excitation at source positions where x s indicate the source positions.
Step 2. Perform backward propagation at each receiver position. We solve Eq. 11 to extrapolate the receiver wavefields where the backward source is recorded data R(x r , t) in time as the boundary condition. The receivers's excitation at the record positions is given as follows: where x r and T max represent the receiver positions and record duration, respectively.
Step 3. Compute the final imaging results using an imaging condition.
Finally, the zero-lag crosscorrelation imaging condition is used to obtain the image of the subsurface structure (Claerbout, 1971) as follows:  where I p (x) represents the imaging result at position x. p s (x, t) and p r (x, t) express the compensated source and receiver wavefields, respectively. To ensure the time consistency between the source and receiver wavefields, the source wavefields should be read from the disk or reconstructed (Ren et al., 2022).

Numerical implementation
The pseudo-spectral method is widely applied to solve the fractional Laplacian operator (Xue et al., 2017;Wang et al., 2018;Wang et al., 2020;Wang et al., 2022). We used the second-order central finite-difference approach and fast Fourier transform (FFT) to calculate the temporal derivative and fractional Laplacian operators, respectively (Bai et al., 2019;Xing and Zhu, 2020;. The detailed numerical implementation is summarized in the following three steps. 1) Calculate the fractional Laplacians. Take (−∇ 2 ) 1+0.5ξ p as an example, and it can be expressed as follows: where F and F −1 denotes the forward and inverse Fourier transform, respectively, and |k| is the norm of the complex wavenumber vector. ζ is an invariable coefficient defined ζ 1/16 in this study (Chen et al., 2016).
2) Calculate the time partial derivative and next moment wavefields.
where p it+1 , p it , p it−1 represents the wavefields at the next, current, and previous timestep, respectively. Lap denotes all Laplacians.
3) Update the wavefields and enter the next time cycle until the maximum simulation time.
Therefore, the Q-compensated wave equation based on the constant fractional-order Laplacians can be implemented.

Accuracy analysis
We conduct a numerical test to analyze this approximation accuracy because the spatial constant-order fractional Laplacian is approximated using Taylor expansion. Here, we take Eq. 8 as an example. This test is performed using a 2D homogeneous model with a grid size of 200 × 200 nodes and spatial intervals of 10 m. The reference velocity is 4,000 m/s. For simplicity, we set λ 1, ζ 1/16, β 1. The degree of match between our approximation solution (the right-hand side of Eq. 8) and the original solution (left-hand side of Eq. 8) is evaluated using the relative root mean square error (RRMS) defined as follows: where x ai and x bi represent the left and right sides of Eq. 8, respectively. Figure 1 shows the RRMS of different simulation parameters, in which Q ranges from 10 to 100 and f d (dominant frequency) ranges from 5 to 50 Hz. The value of RRMS decreases with an increase in Q (Figure 1), and it is almost unchanged as the frequency varies. Furthermore, our approximation has a small RRMS even for the lower Q, validating the accuracy of our approximation. Furthermore, we compare several numerical simulations of different wave equations to investigate the accuracy of Eq. 11 in wavefield extrapolation. A Ricker wavelet is used as the source and is located at the center of a homogeneous model. The single traces of the different wave equation wavefields at x = 1 km are shown in

Numerical examples
Two-layer model Different from Wang et al. (2021) method, the proposed method avoids calculating the spatial variable-order fractional Laplacians. Therefore, it is more straightforward for dealing with the sharp Q model, which can be verified by the two-layer model (Figure 3). The grid is 400 × 400 cells with a unified 5 m interval. Ricker wavelet with a peak frequency of 30 Hz is the location in the center of the model FIGURE 3 Model parameters for the two-layer models.

Frontiers in Earth Science
frontiersin.org and has a time step of 0.5 ms. The compensated parameter β is 2 and σ 0.5c 2γ−1 0 ω −2γ 0 sin(πγ)|k ref | 2γ+1−β . We run simulations for three different cases. In the first case, the fractional Laplacians in Eq. 1 are calculated by the pointwise strategy (the reference solution). In the second case, the fractional Laplacians are calculated by using the average Q value scheme (the average Q-compensated scheme). The third case is implemented by solving Eq. 11 (the proposed constant-order compensated method). The reference solution, average Q-compensated, and constant-order compensated schemes are shown in Figures 4A-C. The difference between the reference solution and the average Q-compensated is shown in Figure 4D, and the difference between the reference solution and the constant-order compensated schemes is shown in Figure 4E. Note that the color scales are the same for all snapshots. Compared to Figure 4D, the smaller residual in Figure 4E confirms the accuracy of the proposed method for sharp Q media.

Marmousi model
Furthermore, we use the Marmousi model to validate the stability and reliability of the proposed scheme in a strongly heterogeneous Q model. The velocity and Q models that are discretized into 400 × 200 points with 10 m spacing are shown in Figures 5A, B. A Ricker wavelet with a peak frequency of 25 Hz is selected as the point source. The simulation time is 2.5 s with a time  Figure 6 where the first row corresponds to the elastic media and the second row correspond to the loss media. Columns 1-4 represent to the 10th, 30th, 50th, and 70th shots, respectively.
All common-gathers are displayed at the same color scales. As expected, the deep structural energy in loss medium is obviously weakened due to the existence of attenuation. Figure 7 shows the migration profiles with different methods. All profiles are displayed at the same color scales. We tested four different imaging methods: the acoustic RTM to acoustic data (reference solution, Figure 7A), the acoustic RTM to

FIGURE 6
Common-gathers. The first row corresponds to the elastic media and the second row correspond to the loss media. Columns 1-4 represent to the 10th, 30th, 50th, and 70th shots, respectively.

Frontiers in Earth Science
frontiersin.org  Frontiers in Earth Science frontiersin.org 08 viscoacoustic data (non-compensated RTM, Figure 7B), the average Q-compensated RTM ( Figure 7C), and the proposed constant-order compensated RTM (our method, Figure 7D). Obviously, significant amplitude reduction (particularly for deep structures) is observed ( Figure 7B) because of the attenuation property of the medium. Two compensation methods improve the imaging results without numerical instability. Although both compensation results show clear anticlinal structures because of the enhanced high-frequency components, there is still some difference between the average

Frontiers in Earth Science
frontiersin.org 09 Q-compensated RTM and the reference solution. In contrast, Figure 7D is very similar to the reference solution.
We also test the possible effects of σ on Q-RTM results in the Marmousi example and show results with the different σ in Figure 8. We set ref c  Figure 8A) and migration results with the very large σ damages too much high-frequency information ( Figure 8D). Figures 8B, C show better imaging results. Hence, it means that the more high-wavenumber signal is eliminated with σ increases.
For comparison, Figures 9A, B show vertical profiles at 1.5 km and 3.5 km along the horizontal direction, respectively. The black line, blue dashed line and red dashed line represent the reference solution, average Q-compensated, and constant-order compensated RTM, respectively. The blue dashed line suffers from amplitude mismatch and phase distortion phenomenon, while the red dashed line matches the reference trace better (particularly for deep reflections). Figure 10 represents the corresponding amplitude spectra. Compared with the non-compensated RTM (green line), the high-wavenumber components of the average Q-compensated RTM (blue dash line) and the constant-order compensated RTM (red line) are significantly improved. Furthermore, confirming the superiority of the proposed method in a strongly heterogeneous Q model, the consistency between the red and black lines is slightly better than that between the blue and black lines.

3D overthrust model
We verify the applicability of the proposed method in the 3D case. The velocity and Q models are shown in Figures 11A, B, respectively. The model contains 550 × 200 × 150 cells with unified spatial intervals of 10 m. A total of 270 shots and 22,500 receivers are evenly located. We simulate records using a Ricker wavelet with a

Frontiers in Earth Science
frontiersin.org peak frequency of 20 Hz. The records have a duration of 2 s with a time interval of 1 ms, and the compensated parameter β is 8. We set a maximum stacking aperture of 1.5 km for each shot to eliminate the diffraction artifacts from the long offset. Different migrated images with acoustic RTM to acoustic data (reference solution), acoustic RTM to viscoacoustic data (non-compensated), average Q-compensated RTM, and constant-order compensated RTM are shown in Figures 12A-D. Unsatisfactory imaging results, particularly below the high dip normal fault, indicate that attenuation has a certain impact on 3D migrated images, as shown in    Figure 13, where the black line, blue dashed line and red dashed line represent the reference solution, average Q-compensated, and constant-order compensated RTM, respectively. The blue line has a shifted phase due to inaccurate attenuation compensation, affecting the precise identification of the target layer and horizon interpretation. In contrast, the red line maintains polarity consistency with the reference solution, demonstrating the effectiveness of the proposed method in the 3D case.

Field data from a land survey
To confirm the potential of our approach in field applications, we apply the proposed constant-order compensated RTM to the field data from a land survey. The velocity (Sava and Vlad, 2008) and Q models (Tonn, 1991) are shown in Figures 14A, B, respectively. The model has 900 × 300 cells with a grid size of 20 m 2 × 20 m 2 . The

FIGURE 18
Amplitude spectra corresponding to the migrated seismic traces (A) at x = 5 km; (B) at x = 9 km. The record duration is 5 s, with a time step of 1 ms. A 2D acquisition line contains 240 excitation sources. We set the maximum stacking aperture of 3.0 km for each shot to eliminate the diffraction artifacts from the long offset, and the compensated parameter β is 6. Figures 15A-D show the common-shot gathers of 50, 100, 150, and 200. All the shots were preprocessed, including static correction, surface wave attenuation, multiple attenuations, and bandpass filtering.
The acoustic RTM (non-compensated) and the constant-order compensated RTM are shown in Figures 16A, B. Compared to the acoustic RTM results ( Figure 16A), the proposed compensated RTM ( Figure 16B) shows a higher resolution and has better illumination of the deep reflections. Specifically, compared with Figure 16B, the continuity of the structures in Figure 16A is destroyed, and the lateral variation in the reflectors is worse (approximately 2 km depth denoted by the red arrow) since acoustic RTM ignores the attenuation compensation. This can be considered a fake fault, thus degrading the reliability of the interpretation. A zoomed-in section of the red box in Figure 16 is shown in Figure 17, with Figures 17A, B showing the acoustic and the constant-order compensated RTM, respectively. The reflections of Figure 17B are stronger than Figure 17A. The wavenumber spectra of the images at x = 5 and 9 km are shown in Figures 18A, B, respectively. Indicating that the proposed method can effectively recover the amplitude, the wavenumber components are increased visibly after Q compensation (Figure 18).

Discussion
The accuracy of migration imaging is essential for seismic data processing and structural interpretation, and calculation efficiency is the premise of industrial practicality. Here, we discuss the cost of several compensation schemes and their precision of imaging results. The numerical examples are implemented using the Compute Unified Device Architecture programming on an Nvidia Geforce RTX 2080Ti. We run simulations for three different cases: the exact solution, the average Q-compensated method, and the constant-order compensated scheme. In Table 1, we list the computation time with the different methods for the 2D Marmousi and 3D overthrust models. For quantitative comparison, we calculate the mean relative absolute error and mark them on the right side of Table 1. The model parameters are consistent with the previous test. Due to the huge amount of computation, the exact solution (calculated by pointwise FFT) only runs a shot simulation. As shown in Table 1, even though the exact solution can accurately simulate the propagation of seismic waves, its numerical implementation is more expensive. Compared with the pointwise FFT, the average Q scheme significantly improves computing efficiency under the condition of sacrificing computational accuracy. Although the numerical implementation of the constant-order compensated schemes is slightly more expensive (about 1.3 times slower) than the average Q scheme in the 3D case, the calculation accuracy was significantly improved. As mentioned above, the computational cost and accuracy test demonstrated that the proposed method provides a better trade-off between imaging accuracy and calculation efficiency. In this paper, we do not test the 3D field data because a single GPU card cannot afford its memory requirements. Hence, the future research direction on this topic is the 3D field data Q-compensated RTM by multi-GPU computing based on model segmentation.

Conclusion
We presented a stabilized Q-compensated RTM scheme which has the explicit stabilization terms in the time-space domain. The proposed Q-compensated RTM avoids the compensating error introduced by averaging spatially varying fractional orders, and it enables us to precisely compensate for the seismic attenuation in complex heterogeneous Q media. The proposed algorithm enhances the resolution in Q-RTM without significantly increasing computational cost. The explicit stabilization term is free from frequent Fourier transforms, so it is expected to simplify the Q-RTM workflow. Numerical simulation examples for homogeneous models demonstrated that the numerical solutions of the proposed wave equation agree with those of the original viscoacoustic wave equation. Furthermore, the synthetic and land field datasets demonstrate the superiority and effectiveness of the proposed approach. We anticipate that the proposed Q-compensated RTM will directly benefit imaging applications as well as enhance the reliability of seismic interpretations.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Funding
This work is partly supported by the National Natural Science Foundation of China (41930431, 42274171, 42204129), the Joint Guiding Project of the Natural Science Foundation of Heilongjiang Province (LH 2021D009).

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.