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 post-stack 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 frequency-independent 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; Zhang and Gao, 2022). Alternatively, some mathematical-model-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 amplitude-loss 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 inverse-physical 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 (Zhu and Harris, 2014; 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; Zhang et al., 2022; Zhang and Gao, 2022). 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 (mixed-domain 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 variable-order 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.

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: 1 c 2 ∂ 2 p ∂ t 2 = η a + τ a + η c 1 + η c 2 + τ c (1) and η a = − c 0 2 γ ω 0 − 2 γ ⁡ cos π γ − ∇ 2 γ + 1 p , τ a = c 0 2 γ − 1 ω 0 − 2 γ ⁡ sin π γ ∂ ∂ t − ∇ 2 γ + 0.5 p , η c 1 = c 0 2 γ − 1 ω 0 − 2 γ ⁡ sin π γ σ − ∇ 2 2 γ + 1 + β / 2 p , η c 2 = − σ 2 c 2 − ∇ 2 β p , τ c = − 2 σ c 2 ∂ ∂ t − ∇ 2 β / 2 p , γ = 1 π Q , c = c 0 ⁡ cos π γ / 2 , (2) 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. η c 1 and η c 2 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 (spatial-wavenumber domain) operators, the average Q value is often used when extrapolating the wavefields (Zhu et al., 2014; Zhu and Harris, 2014; 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 frequency-wavenumber domain viscoacoustic wave equation 1 c 2 ⁡ cos π γ ∂ 2 p ∼ ∂ t 2 = η ∼ a + τ ∼ a + η ∼ c 1 + η ∼ c 2 + τ ∼ c , (3) and η ∼ a = − c 0 2 γ ω 0 − 2 γ k 2 γ + 2 p ∼ , τ ∼ a = c 0 2 γ − 1 ω 0 − 2 γ 1 Q k 2 γ + 1 ∂ p ∼ ∂ t , η ∼ c 1 = c 0 2 γ − 1 ω 0 − 2 γ σ Q k 2 γ + 1 + β p ∼ , η ∼ c 2 = − σ 2 k 2 β c 2 ⁡ cos π γ p ∼ , τ ∼ c = − 2 σ k β c 2 ⁡ cos π γ ∂ p ∼ ∂ t (4) where p ∼ represents the wavefield in the wavenumber domain. Here, the fractional Laplacian c 0 2 γ − 1 ω 0 − 2 γ σ Q k 2 γ + 1 + β in Eq. 3 is taken as an example of an approximation c 0 2 γ − 1 ω 0 − 2 γ σ Q k 2 γ + 1 + β = λ σ Q c 0 k 1 + β k k d 2 γ (5) where f d and k d represent the dominant frequency and dominant wave number, respectively, λ = ω d / ω 0 2 γ , ω d = 2 π f d . Because 2 γ ln k k d ≪ 1 , we use Taylor expansion to approximate Eq. 3 as λ σ Q c 0 k 1 + β k k d 2 γ ≈ λ σ Q c 0 k 1 + β 1 + 2 γ ⁡ ln k k d , (6) 1 + 2 γ ⁡ ln k k d can be approximated by Taylor expansion again 1 + 2 γ ⁡ ln k k d = 1 + 2 γ ζ ln 1 + k k d ζ − 1 ≈ 1 + 2 γ ζ k k d ζ − 1 , (7) where ζ is a empirical constant coefficient (Chen et al., 2016), which is introduced to guarantee k k d ζ − 1 close to zero. Therefore, the left-hand side of Eq. 5 can be approximated as follows: c 0 2 γ − 1 ω 0 − 2 γ σ Q k 2 γ + 1 + β ≈ λ σ Q c 0 k 1 + β 1 + 2 γ ζ k k d ζ − 1 (8)

The other spatial variable-order Laplacians in Eq. 4 can be approximated in the same way. Therefore, Eq. 3 can be expressed as follows: 1 λ c 2 ⁡ cos π γ ∂ 2 p ∼ ∂ t 2 = η ∼ a n e w + τ ∼ a n e w + η ∼ c n e w + 1 λ ⁡ cos π γ η ∼ c 2 + 1 λ ⁡ cos π γ τ ∼ c (9) where η ∼ a n e w = − 1 − 2 ξ π Q k 2 p ∼ − 2 ξ π Q c 0 ω d ξ k 1 + 0.5 ξ p ∼ , τ ∼ a n e w = 1 Q c 0 ∂ ∂ t 1 − 2 ξ π Q k 0.5 + 2 ξ π Q c 0 ω d ξ k 0.5 + 0.5 ξ p ∼ , η ∼ c n e w = σ Q c 0 1 − 2 ξ π Q k 1 + β / 2 + 2 σ ξ π Q c 0 ω c ζ k ξ + 1 + β / 2 p ∼ , (10)

By transforming Eq. 9 back to the space domain, we obtain the following equation 1 λ c 2 ⁡ cos π γ ∂ 2 p ∂ t 2 = η a n e w + τ a n e w + η c n e w + 1 λ ⁡ cos π γ η c 2 + 1 λ ⁡ cos π γ τ c (11) where η a n e w = 1 − 2 ξ π Q ∇ 2 p − 2 ξ π Q c 0 ω d ξ − ∇ 2 1 + 0.5 ξ p , τ a n e w = 1 Q c 0 ∂ ∂ t 1 − 2 ξ π Q − ∇ 2 0.5 + 2 ξ π Q c 0 ω d ξ − ∇ 2 0.5 + 0.5 ξ p , η c n e w = σ Q c 0 1 − 2 ξ π Q − ∇ 2 1 + β / 2 + 2 σ ξ π Q c 0 ω d ξ − ∇ 2 ξ + 1 + β / 2 p . (12)

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).

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.

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.