Next Article in Journal
A Review of Models for Photovoltaic Crack and Hotspot Prediction
Previous Article in Journal
Properties of Lightweight Controlled Low-Strength Materials Using Construction Waste and EPS for Oil and Gas Pipelines
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Finite Difference Pre-Stack Depth Migration of One-Way Wave Equation in Isotropic Visco-Acoustic Media of a Ray-Centered Coordinate System

School of Ocean and Earth Science, Tongji University, Siping Road, Shanghai 200092, China
*
Author to whom correspondence should be addressed.
Energies 2022, 15(12), 4302; https://doi.org/10.3390/en15124302
Submission received: 6 May 2022 / Revised: 2 June 2022 / Accepted: 10 June 2022 / Published: 12 June 2022
(This article belongs to the Section H: Geo-Energy)

Abstract

:
In seismic exploration, obtaining accurate bandwidth of the reflected signal is essential for determining seismic resolution. A portion of the stored energy is attenuated during signal propagation, narrowing the received seismic signal bandwidth. Therefore, compensating the energy attenuation is important for improving the seismic resolution. The current method of compensating absorption and attenuation based on a single channel can only compensate the post-stack data (self-exciting and self-receiving), whereas in practice, seismic waves do not propagate along the self-exciting and self-receiving seismic wave path; the propagation path is complex. The absorption and attenuation depend on the propagation path. The primary methods used for Q-compensation along the propagation paths are one-way wave extrapolation in the Cartesian coordinate system and Gaussian beam Q-compensation migration in the ray-centered coordinate system. However, the large angle limits the one-way wave method, and the Gaussian beam method refers to the high-frequency approximation solution of the two-way wave equation. Therefore, a 15-degree equation in the ray-centered coordinate system is proposed. Seismic waves extrapolate along the ray, which compensates the absorption and attenuation along the real propagation path. The 15-degree equation in the ray-centered coordinate system does not perform high-frequency approximation in the ray beam and has no large angle limit, facilitating the accurate description of local wavefields in the ray beam.

1. Introduction

Seismic data obtained from exploration areas with strong absorption and attenuation properties feature low frequencies and weak amplitudes. In addition, the imaging results of traditional imaging algorithms, without considering the absorption and attenuation compensation in such exploration areas (such as gas-bearing reservoirs), decrease regarding resolution (Traynin et al., 2008 [1]). To improve the imaging accuracy of this type of exploration area, the absorption and attenuation effect during the imaging process must be compensated.
The most commonly used seismic wave absorption and attenuation compensation methods include the Q deconvolution and visco-acoustic medium pre-stack depth migration (PSDM) methods. The former is used for the zero-offset seismic channel (or post-stack data), whereas the latter is used for the change in Q value and the absorption and attenuation compensation along the wave propagation path. The inverse Q migration processing of pre-stack seismic data is conducted in the following ways: Inverse Q-compensation is conducted when the upward received wavefield is extended backward, and the inverse Q-compensation is also conducted when the source wavefield is extended downward; the imaging is then conducted for the primary wave under certain imaging conditions.
Deconvolution is a standard method for improving the resolution of seismic imaging in the traditional seismic data processing. The optimal Wiener filter is a traditional deconvolution method used to identify the optimal solution regarding the least square. However, it assumes the time invariance theory. When the seismic wave propagates through an underground medium, the wavelet’s shape is time-varying. In areas with severe absorption and attenuation, the achievement of ideal results using the general Wiener deconvolution method is difficult. Futterman (1962) [2] proposed a constant Q medium model based on the linear absorption theory. Kjartansson (1979) [3] introduced the frequency-dependent complex modulus to improve the linear seismic wave absorption model. In this model, absorption is proportional to a certain power of a frequency function, which better reflects the principle of causality in seismic wave absorption and attenuation. Wang (2002, 2008) [4,5] studied seismic wave attenuation and inverse Q filter compensation through the numerical wavelet test, and the Q value obtained from the Vertical Seismic Profiling (VSP) downward wavefield was applied to the surface data to improve resolution. Chopra et al. (2003) [6] proposed a high-frequency recovery method, using the Wiener shaping filter to obtain an inverse filter operator to describe the formation absorption effect from VSP records, and applied it to the post-stack data to improve the profile resolution. The time-varying Wiener filter proposed by van der Baan (2012) [7] can also achieve the same absorption and attenuation compensation effect as the inverse Q filter. Braga and Moraes (2013) [8] applied and realized the inverse Q filtering method in the wavelet transform domain and determined the amplitude damping coefficient according to the signal-to-noise ratio of seismic data for amplitude attenuation compensation and velocity dispersion correction. Ren et al. (2007) [9] proposed a method of Q-compensation along the ray path defined by the root-mean-square velocity in the Common Middle Point (CMP) channel gather, which improved the deficiency of the Q-compensation for the zero-offset seismic channel. Li et al. (2015) [10] proposed a compensation method that can compensate the horizontal and vertical depth absorption attenuation at the same time. Starting from the absorption attenuation operator, this method separates the attenuation operator into the attenuation operator related to the horizontal and vertical absorption attenuation, respectively, and then applies the separated operator to the horizontal and vertical absorption attenuation compensation, respectively, so as to realize the absorption attenuation compensation algorithm.
Single channel inverse Q filtering bears certain advantages, including simplicity and improved wavelet property after processing. However, the absorption and attenuation of seismic waves occur along the propagation path, as a result, a better compensation method is based on the PSDM method in visco-acoustic media (or even visco-elastic media). Mittet et al. (1995) [11] proposed the use of the MeClellan transformation to conduct PSDM in visco-acoustic media in the frequency-space domain, which increases the calculation amount. Wang (2004) [12] achieved PSDM imaging of visco-acoustic media for a common shot point gather.
In addition, reverse time migration (RTM) is a widely used imaging method. When RTM performs Q-compensation migration, it is necessary to modify the two-way wave equation (Suh et al., 2012 [13]; Dutta et al., 2014 [14]; Zhu et al., 2014 [15]). Chen et al. (2020) [16] developed a robust stabilization strategy for Q-RTM by incorporating a time-variant filter into the amplitude-boosted wavefield extrapolation step. Yang et al. (2021) [17] developed a robust space-wavenumber compensation operator to overcome numerical instabilities for seismic data with strong high-frequency noise and applied it to visco-acoustic RTM. Fathalian et al. (2021) [18] set up a procedure for solving the time-domain visco-acoustic wave equation for tilted transversely isotropic (TTI) media, based on a standard linear solid model and, from this, develop a visco-acoustic reverse time migration (Q-RTM) algorithm. However, numerical solution of two-way wave equation has some problems: for example, for high-frequency data, numerical dispersion is especially difficult to overcome and so on. Sun et al. (2017) [19] developed two attenuation compensation operators that correct for both amplitude loss and velocity dispersion in a stable manner and when applied to RTM, it leads to two kinds of Q-compensation imaging conditions.
This paper uses the 15-degree equation in the ray-centered coordinate system. In the ray-centered coordinate system, the wavefield propagates along the ray. In isotropic media, the ray direction represents the propagation direction of seismic wave. As a result, we can carry out Q-compensation along the propagation path of seismic wave. In addition, since the 15-degree equation in the ray-centered coordinate system is solved in the frequency-space domain, Q attenuation and Q-compensation can be realized as long as the real velocity is changed to complex velocity, while we need to modify the two-way wave equation accordingly for RTM. Finally, using the 15-degree equation can accurately describe the local wavefield in the beam, while the wavefield obtained by Gaussian beam is the high-frequency approximation solution of the two-way wave equation, which can not accurately describe the local wavefield (Zhang et al., 2021 [20]).
Wang (2004) [12] emphatically discussed the control challenge associated with amplitude compensation stability in the wavefield continuation process. This method uses the quality factor of the overlying strata of the seismic trace, where the imaging point is located to determine the total amplitude compensation of each frequency component and the amplitude stabilization compensation within the current extrapolation step. Considering that only the absorption effect of the single channel overlying strata is used, the method does not strictly follow the propagation path to compensate the absorption effect for the post-stack wavefields. In addition, although the amplitude stabilization compensation in the current extrapolation step is determined by the total amplitude compensation, the wavefield is recursively extended; therefore, the influence of stability control (lowpass gain effect) is cumulative. When the extrapolation depth is large, the amplitude compensation for the high-frequency components may be overly suppressed and thus may affect the imaging resolution and have a negative effect.
Theoretically, the inverse Q migration compensates and corrects the attenuation effect of inelastic absorption in the media in strict accordance with the wavefield propagation path. However, the compensation operator of the attenuation amplitude in the inverse Q migration increases exponentially as frequency increases, and the calculation stability issue occurs when the quality factor is small with a long extrapolation time. Two methods are commonly used to control high-frequency amplitude compensation instability—the gain cut-off frequency and stability factor methods. Both methods control the stability of the inverse Q-compensation by limiting the effective gain cut-off frequency.
In this study, based on the 15-degree equation beam propagation in the ray-centered coordinate system, absorption and attenuation compensation along the wave propagation path is conducted, offering three advantages:
(1)
This method provides Q-compensation along the real propagation path of seismic waves in isotropic media;
(2)
In terms of mathematical derivation, this method is simpler than the Q-compensation imaging method of RTM;
(3)
Compared with the Gaussian beam method, this method does not introduce high-frequency asymptotic (Zhang et al., 2021 [20]).
(4)
Numerical experiments demonstrate the effectiveness of the proposed method.

2. Seismic Wave Absorption and Attenuation Mechanism

The acoustic and elastic wave equations are used to describe wave propagation, and seismic waves are assumed to propagate in adiabatic media—the wave continues to propagate indefinitely upon generation by a particular source. Therefore, the wave may attenuate in space as it spreads from its source, but the total energy of the particle motion remains constant.
The propagation of seismic waves in real media differs from that in an ideal medium. The seismic wave amplitude and geometric diffusion factors are attenuated, and accompanied by chromatic dispersion, resulting in a change in spectrum components and the seismic waveform. The absorption and dispersion characteristics of seismic waves in real media are closely related to physical conditions, such as material composition, saturation and porosity of rock strata; the overall effect of such absorption and attenuation can be described by a dimensionless quantity (i.e., the quality factor, Q, of medium).
To evaluate the effect of attenuation, the effect of Q in frequency domain on wave propagation can replace the real angular frequency ω with a complex angular frequency; the influence of Q in the wavenumber domain on wave propagation can be regarded as a substitution for the complex wavenumber of real wavenumber k . In the first case, the complex angular frequency, ω , must have an imaginary part, i ω 2 Q t ; in the second case, the complex wavenumber, k , must have an imaginary part i k 2 Q s .
The complexity of real media complicates the description of absorption and attenuation of seismic waves for the viscous property of the media with a mathematical physical model. The absorption, attenuation and dispersion curves of the K-F and Kjartansson (K-J) models are the most commonly used in the compensation of absorption and attenuation in exploration seismology.
The mathematical expression of the K-J model is defined in Equations (1)–(4):
γ = 1 π arctan 1 Q
α x , z , ω = tan π γ 2 ω v 0 ω ω 0 γ
K c ω = K c R ω i K c I ω
K c R ω = ω v 0 ω ω 0 γ cos π γ 2 K c I ω = ω v 0 ω ω 0 γ sin π γ 2
where, v 0 represents the reference velocity and is used as the acoustic velocity; ω 0 refers to the reference frequency, which is considered to be the main frequency of the wavelet; Q represents the quality factor; α is the amplitude attenuation factor; γ is the loss angle; K c R ω is the real part of the complex wavenumber K c ω ; and K c I ω is the imaginary part of the complex wavenumber K c ω .
We used the K-J model in this study. The complex wavenumber expressed in Equation (3) can be derived, and the absorption and attenuation of seismic wave propagation can be described.
Wave propagation in the absorption and attenuation media is irreversible. The absorption and attenuation compensation by the K-J model is commonly used in seismic exploration as its effects are evident. However, no strict mathematical physics theorem guarantees the absolute rationality of absorption and attenuation compensation. The effect of inverse Q filtering and the effect of absorption and attenuation compensation of PSDM are used to judge the appropriateness of the absorption and attenuation compensation.
In addition, the scattering phenomena alter the amplitude attenuation and phase change in the observed seismic waves. O’Doherty and Anstey (1971) [21] noted that stratigraphic filtering would occur when seismic waves pass through fine stratigraphic strata, resulting in the attenuation of seismic wave amplitudes and phase changes, which vary with different offsets. However, the amplitude attenuation and phase change in the seismic waves caused by these two phenomena are not easily distinguishable. Therefore, the Q value estimated from seismic wave amplitude attenuation should be the equivalent Q value of these two effects.

3. Stable Inverse Q PSDM in the Ray-Centered Coordinate System

The core concept of this study is the compensation of the inelastic absorption and attenuation along the wave propagation path; therefore, the wavefield extrapolation equation in the visco-acoustic medium is generated first.

3.1. Wavefield Continuation with Attenuation Compensation in the Ray-Centered Coordinate System

The ray-centered coordinate system is centered on the ray, which is an orthogonal coordinate system at each point on the ray (Zhang et al., 2021 [20]).
The dispersion equation (Equation (5)) of the two-way wave in the ray-centered coordinate system of the isotropic acoustic medium is expressed as follows (Cheng et al., 2012 [22]):
ω 2 V 2 = K ω 2 = c s s k s 2 + c q q k q 2 i c s k s i c q k q
According to the above expression, the real wavenumber K ω is replaced by the complex wavenumber K c ω , and the dispersion equation of the two-way wave in the ray-centered coordinate system of isotropic visco-acoustic media can be obtained using Equation (6):
K c ω 2 = c s s k s 2 + c q q k q 2 i c s k s i c q k q
where s and q are the ray-centered coordinates; K c ω is the complex wavenumber, and it is assumed to be independent of the spatial position. In addition,
c s s = 1 h 2 c q q = 1 c s = 1 h s 1 h c q = 1 h h q h = 1 + 1 V V q q = 0 q v = V s , 0
Equation (6) follows the expression form of Equation (5) in the acoustic medium; however, the wavenumber changes from real to complex. There are many expressions of the complex wavenumber K c ω , and different mathematical models of absorption and attenuation can define different types of complex wavenumbers.
To obtain the wavefield compensation continuation operator in the ray-centered coordinate system, wavenumber, k s , must be determined (Equation (7))
k s = i c s 2 c s s ± K c ω 2 c s s c s 2 c s s 2 c q q c s s k q 2 i c q c s s k q
The 15-degree equation in the frequency-wavenumber domain can be expressed as Equation (8):
u ˜ s = i k s u ˜
Substituting Equation (7) into Equation (8) yields Equation (9):
u ˜ s = i i c s 2 c s s ± K c ω 2 c s s c s 2 c s s 2 c q q c s s k q 2 i c q c s s k q u ˜
Equation (9) is the one-way wave equation in the ray-centered coordinate system, where the positive sign represents the extrapolation along the negative direction of s ; the minus sign represents the extrapolation along the positive direction of s .
Equation (9) can be broken down in Equations (10) and (11):
u ˜ 1 s = c s 2 c s s u ˜ 1
u ˜ 2 s = ± i K c ω 2 c s s c s 2 c s s 2 c q q c s s k q 2 i c q c s s k q u ˜ 2
The analytic solution of the wavefield in Equation (10) can be obtained directly using Equation (12):
u ˜ 1 k q , s , ω = exp c s 2 c s s s
Equation (12) can then be substituted into Equation (11), which can be solved to control the propagation of the one-way wave in the ray-centered coordinate system. Therefore, we only numerically solve and focus on Equation (11).
Subsequently, in order to eliminate the influence of subscript “2”, we mark u 2 as p .
The analytic solution of Equation (11) is Equation (13):
p ˜ k q , s , ω = exp ± i K c ω 2 c s s c s 2 c s s 2 c q q c s s k q 2 i c q c s s k q s
Therefore, the decay continuation operator is expressed by Equation (14):
p ˜ k q , s + Δ s , ω = exp ± i K c ω 2 c s s c s 2 c s s 2 c q q c s s k q 2 i c q c s s k q Δ s p ˜ k q , s , ω
Equation (14) is a wavefield attenuation continuation operator considering the attenuation effect of the absorption media.
By back propagating the attenuation continuation operator (i.e. replacing Δ s in the right end term of equation 10 with Δ s ), the compensation continuation operator can be obtained (Zhang et al., 2010 [23])
p ˜ k q , s + Δ s , ω = exp ± i K c ω 2 c s s c s 2 c s s 2 c q q c s s k q 2 i c q c s s k q Δ s p ˜ k q , s , ω

3.2. Finite Difference Solution of the One-Way Wave Equation in Frequency-Space Domain of the Ray-Centered Coordinate System

We considered wavefield extrapolation in the positive direction of s (equivalent to the downward wave equation in the Cartesian coordinate system) in Equation (11) as an example. To simplify the calculation, we introduce Equation (16):
k s = ± K c ω 2 c s s c s 2 c s s 2 c q q c s s k q 2 i c q c s s k q
To obtain the one-way wave equation in the ray-centered coordinate system, we apply the Taylor expansion of k s to generate Equation (17):
k s k q = k s k q 0 = 0 + k s k q 0 k q 0 = 0 k q + 1 2 2 k s k q 0 2 k q 0 = 0 k q 2 + O k q 3
Substituting Equation (16) into Equation (17) yields Equation (18):
k s k q k 0 i c q 2 c s s k 0 k q 1 2 k 0 c q 2 c s s k 0 2 c q q c s s k q 2
where k 0 2 = K c ω 2 c s s c s 2 c s s 2 .
Therefore, Equation (11) can be approximated as follows (Equation (19)):
p ˜ D s = i k 0 i c q 2 c s s k 0 k q 1 2 k 0 c q 2 c s s k 0 2 c q q c s s k q 2 p ˜ D
The simplified Equation (19) can be expressed as Equation (20):
p D s = i k 0 p D i c q 2 c s s k 0 p D q + i 2 k 0 c q 2 c s s k 0 2 c q q c s s 2 p D q 2
Equation (20) can be abbreviated as Equations (21) and (22):
p D s = A p D + B p D q + C 2 p D q 2
where
A = i k 0 B = i c q 2 c s s k 0 C = i 2 k 0 c s q 2 c s s 2 c q q c s s
Equation (23) is a derivation of the finite difference solution of Equation (21) (Ma, 1989 [24]):
p p i j + 1 + p i j 2 p s p i j + 1 p i j Δ s p q p i + 1 j + p i j 2 p i j + p i 1 j 2 Δ q = p i + 1 j p i 1 j 2 Δ q 2 p q 2 p i + 1 j 2 p i j + p i 1 j Δ q 2
where i = 1 , , n q ; j = 1 , , n s .
Through derivation, we have produced an implicit discretization of Equation (21).
For the wavefield continuation along the negative direction of s , Equation (22) can be revised to Equation (24):
A = i k 0 B = i c q 2 c s s k 0 C = i 2 k 0 c s q 2 c s s 2 c q q c s s

3.3. Stability Control

According to Equation (9), the one-dimensional (1D) one-way wave equation in the ray-centered coordinate system of frequency-space domain of isotropic visco-acoustic media can be expressed by Equation (25):
p ˜ s = ± i K c ω p ˜
which has the same form as the 1D one-way wave equation in the Cartesian coordinate system in the frequency-space domain of isotropic visco-acoustic media.
Next, we consider the positive extrapolation of s to be an example of the stability of isotropic visco-acoustic media in the ray-centered coordinate system.
The wave equation along the positive extrapolation of s (downward extrapolation wavefield) is expressed as Equation (26):
p ˜ s = i K c ω p ˜
The analytic solution of Equation (26) is expressed as Equation (27):
p ˜ s + Δ s , ω = p ˜ s , ω exp i K c ω Δ s
Substituting Equation (3) into Equation (27) yields Equation (28):
p ˜ s + Δ s , ω = p ˜ s , ω exp K c I ω Δ s exp i K c R ω Δ s
The specific expression of complex wavenumber K c ω (or complex velocity) is determined by the absorption and dispersion model. This study uses the K-J model, and the specific expression of this model is shown in Equations (3) and (4).
Substituting Δ s for Δ s in Equation (27), the expression of wavefield absorption and attenuation compensation can be expressed using Equation (29):
p ˜ s + Δ s , ω = p ˜ s , ω exp K c I ω Δ s exp i K c R ω Δ s = p ˜ s , ω Λ Δ s , ω Ψ Δ s , ω
where Λ Δ s , ω is the amplitude compensation term and Ψ Δ s , ω is the phase dispersion correction term, the respective expressions of which are defined in Equations (30) and (31):
Λ Δ s , ω = exp K c I ω Δ s
Ψ Δ s , ω = exp i K c R ω Δ s
Taking time (quasi-depth) Δ τ = Δ s / v as the extrapolation step and assigning the reference velocity v 0 = v of the constant Q model yields:
p ˜ τ + Δ τ , ω = p ˜ τ , ω Λ Δ τ , ω Ψ Δ τ , ω
where Λ Δ τ , ω is the amplitude compensation term and Ψ Δ τ , ω is the phase continuation dispersion correction term whose respective expressions are defined in Equations (33) and (34):
Λ Δ τ , ω = exp ω ω ω 0 γ sin π γ 2 Δ τ
Ψ Δ τ , ω = exp i ω ω ω 0 γ cos π γ 2 Δ τ
When the Q value is large, γ 1 π Q 0   sin π γ 2 1 2 Q . Then, the amplitude compensation operator (33) can be expressed as Equation (35):
Λ Δ τ , ω exp ω Δ τ 2 Q
According to Equation (34), the correction of the phase velocity dispersion effect can be completed by a phase shift, and the modulus of its correction term is a constant set to 1, making it unconditionally stable. However, Equation (35) shows that the compensation operator of the attenuation amplitude is an exponential function wherein base e increases as the frequency increases. The computational stability problem occurs when the value of quality factor Q is small and the extrapolation time is large.
Considering the control issue associated with the amplitude compensation stability, the amplitude compensation cumulant should be estimated. For the wavefield extrapolation along the positive direction of s (downgoing wavefield), let τ = 0 in Equation (33), p ˜ 0 , ω be the receiving record, and Δ τ becomes the total extrapolation time from the surface in the local beam to the current point.
The gain cut-off frequency method or stability factor method is used to control the amplitude compensation operator. The amplitude compensation operator of the gain cut-off frequency method can be expressed by Equation (36):
Λ S t b _ 1 Δ τ , ω = exp ω ω ω 0 γ sin π γ 2 Δ τ ω < ω q Λ Δ τ , ω q ω ω q
where ω q is the gain cut-off frequency. This study uses the method proposed by Bickel and Natarajan (1985) [25] to calculate the gain cut-off frequency, ω q (Equation (37)):
ω q = C 2 Q Δ τ
where the gain cut-off coefficient C is a constant with a value ≥1. Wang (2002) [4] also presents the same set of criteria for the gain cut-off frequency based on the overall characteristics of the amplitude compensation operator. The gain cut-off frequency, ω q , is applied to the wavefield extrapolation process of absorption and attenuation compensation for PSDM.

3.4. Inverse Q PSDM

This study briefly lists several core formulas for the PSDM of ray beams:
Source wavefield extrapolation can be expressed according to Equation (38):
P S p x , p y ; x , y , z + Δ z ; ω = P S p x , p y ; x , y , z ; ω W + p x , p y ; x , y , z ; ω ; V x , y , z
where, W + p x , p y ; x , y , z ; ω ; V x , y , z is the wavefield estimator projected back to the Cartesian coordinate system from the ray-centered coordinate system. P S p x , p y ; x , y , z ; ω is the local wavefield at the source point. p x , p y specifies beam direction, x , y specifies transverse spreading range, and z represents the wavefield at a certain depth.
The wavefield extrapolation at the receiver can be expressed by Equation (39):
P R p x , p y ; x , y , z + Δ z ; ω = P R p x , p y ; x , y , z ; ω W p x , p y ; x , y , z ; ω ; V x , y , z
The meaning of each quantity in Equation (39) is similar to that of Equation (38).
The imaging formula is expressed by Equation (40):
I m a g e p x , p y ; x , y , z + Δ z = i ω = 1 N ω P S p x , p y ; x , y , z + Δ z ; i ω P R p x , p y ; x , y , z + Δ z ; i ω
where I m a g e p x , p y ; x , y , z + Δ z is the pre-stack imaging result. The sum of the right end term of Equation (40) on the frequency component can only be conducted on the positive frequency.
As the wavefield extrapolation equation is linear, the synthesized local wavefield can be directly extrapolated by the wavefield extrapolator and the imaging results of Equation (40) can complete the finite difference PSDM of the 15-degree equation beam in the ray-centered coordinate system of the visco-acoustic media.
The vertical resolution of the seismic imaging is primarily determined by the wavelet bandwidth and the spread of the dominant frequency band. The bandwidth of the wavelet and the spread of the dominant frequency band are affected by the inelastic factors of the media in the propagation process, causing the elimination of the absorption, attenuation and dispersion of the seismic wave caused by inelastic factors in the full path from the source point to the reflection point and then to the receiver point. Compared with the general migration method, the inverse Q migration focuses on the compression of the seismic wavelet to obtain higher vertical resolution imaging results. Meanwhile, the deep seismic wave’s attenuation energy can be compensated to a certain extent.
For post-stack migration, the seismic profile is assumed to be self-exciting and self-receiving, and the upward and downward paths coincide. With inverse Q migration, the mere compensation of the upward wave is sufficient; however, the migration velocity is halved.
For pre-stack migration, the amplitude compensation of the upward wave extrapolation is simpler when the common shot gather migration is conducted based on the single square root operator in the ray-centered coordinate system. However, this only compensates for the attenuation and dispersion of energy absorption from the reflection point to the receiver point, whereas the attenuation and dispersion from the source point to the reflection point are not compensated. In this study, forward modeling is conducted from the source point to the reflection point without Q, yielding wavefield P 1 , and then conducted with Q, yielding wavefield P 2 . The ratio of P 1 / P 2 is the quantity of the amplitude attenuation and the dispersion of the corresponding downward local spherical wave, which is multiplied with the extrapolation result of the upward wavefield for the corresponding local spherical wave to obtain the amplitude absorption, attenuation and dispersion compensation of the two-way path (as shown in Figure 1). Stability control involves the use of Equation (37) to calculate the appropriate gain cut-off frequency in each extrapolation step to limit the frequency involved in extrapolation.

4. Numerical Examples

The following numerical experiments verify the validity of the theory and implementation of inverse Q migration in the ray-centered coordinate system of the visco-acoustic media. We use theoretical data and field data, respectively.

4.1. Theoretical Data

First, a trapezoidal depression model is numerically tested. The velocity model is shown in Figure 2. We use a homogeneous Q model of 10. In this numerical experiment, acoustic data (without absorption and attenuation) and visco-acoustic data (with absorption and attenuation) are used for different migration methods. Figure 3a shows the acoustic migration result of the acoustic data, and this result is the potentially ideal reference result. Figure 3b shows the result obtained from the visco-acoustic wave simulation data using the acoustic migration algorithm. The wavelet after imaging exhibits severe energy attenuation and decreased resolution. Figure 3c shows the result of the method in this study (the visco-acoustic migration algorithm is used to process the visco-acoustic data). Figure 3c also shows that the event’s energy absorption and resolution reduction are compensated. Figure 4 shows a comparison of the imaging results of the three migration results at different Common Depth Point (CDP) points. The amplitude of the visco-acoustic data using acoustic migration exhibits severe attenuation, which is expected to affect the seismic geological interpretation result. However, the method presented in this study can compensate absorption and attenuation, making the imaging result comparable to the imaging quality with no attenuation.
To simulate geological anomalies with strong absorption and attenuation characteristics, such as a gas-bearing region, we design a bottom splitting model and modify the Hess model. The bottom splitting velocity and Q models are shown in Figure 5. The gas-bearing region typically has a local spatial distribution, and the corresponding Q model has a small quality factor (severe absorption and attenuation) in the designed gas-bearing region, while other regions have a large quality factor (non-absorption and non-attenuation). Figure 6a shows the result obtained by the acoustic migration algorithm based on the visco-acoustic wave simulation data. Figure 6a also shows severe energy attenuation and reduced resolution in the lower boundary and area of the gas-bearing region. Therefore, the amplitude of the strata under the cover is distorted, and the continuity of the event has deteriorated. Figure 6b shows the result of the method in this study (visco-acoustic migration algorithm from visco-acoustic data). The energy absorption and resolution decline caused by the gas-bearing region are compensated, and the continuity and energy consistency of the event are restored.
To further verify the robustness of the method proposed in this paper, similar to the above bottom splitting model, we adopt a more complex Hess model—to simulate geological anomalies such as a gas-bearing region. Thus, we reduce the velocity of the salt from the original Hess model. Figure 7 shows the velocity and Q models. Figure 8 shows the imaging results. In the gas-bearing region, the amplitude of the visco-acoustic data is severely attenuated by the acoustic migration. The method proposed in this study can compensate the absorption and attenuation.

4.2. Field Data

In order to further test the practicability of our Q-compensation method, this paper selects two-dimensional field data from an exploration area in China. Figure 9a shows the imaging result without Q-compensation and Figure 9b shows the imaging result with Q-compensation.
Comparing the results before and after Q-compensation, it can be seen that the Q-compensation strategy proposed in this paper can effectively compensate the attenuation and deep amplitude energy in the imaging process.
Collectively, these experiments show that the method proposed in this study can solve the problem of inverse Q migration imaging.

5. Conclusions and Discussions

Because of the gas overflow or existence of shallow gas, the imaging of fuzzy areas is common in many exploration areas. However, the imaging of absorption and attenuation compensation using visco-acoustic media PSDM effectively improves the imaging of these fuzzy areas. Amplitude gain control is key for optimal imaging. The proposed method yielded remarkable results in the numerical experiments of inverse Q migration. This study discusses post-stack migration and how to achieve the full path absorption and attenuation compensation from the source point to the reflection point (i.e., image point) and then to the receiver point in pre-stack migration based on a one-way wave propagator in the ray-centered coordinate system.
Since wavefield propagation is coincident with the extrapolation direction in isotropic media, we can use inexpensive 15-degree finite-difference extrapolators to achieve high accuracy extrapolation and perform Q-compensation imaging along the wavefield propagation path.
In the future, referring to the idea similar to the method in this paper, we can easily extend our method to realize Q-compensation imaging in anisotropic viscous quasi-acoustic media along a wavefield propagation path.

Author Contributions

Conceptualization, H.W.; methodology, H.W. and B.Z.; validation, B.Z.; funding acquisition, H.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key R&D Program of China (grant number: 2018YFA0702503, 2019YFC0312004), National Natural Science Foundation of China (grant number: 41774126, 42074143) and Southern Marine Science and Engineering Guangdong Laboratory (Zhanjiang) (grant number: ZJW-2019-04).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

The authors appreciate the sponsors of the WPI (Wave Phenomena and Intelligent Inversion Imaging) group for their financial support and help.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Traynin, P.; Liu, J.; Reilly, J.M. Amplitude and bandwidth recovery beneath gas zones using Kirchhoff prestack depth Q-migration. In SEG Technical Program Expanded Abstracts 2008; Society of Exploration Geophysicists: Houston, TX, USA, 2008; pp. 2412–2416. [Google Scholar] [CrossRef]
  2. Futterman, W.I. Dispersive body waves. J. Geophys. Res. 1962, 67, 5279–5291. [Google Scholar] [CrossRef]
  3. Kjartansson, E. Constant Q-wave propagation and attenuation. J. Geophys. Res. 1979, 84, 4737–4748. [Google Scholar] [CrossRef]
  4. Wang, Y. A stable and efficient approach of inverse Q filtering. Geophysics 2002, 67, 657–663. [Google Scholar] [CrossRef]
  5. Wang, Y. Inverse- Q filtered migration. Geophysics 2008, 73, S1–S6. [Google Scholar] [CrossRef]
  6. Chopra, S.; Alexeev, V.; Sudhakar, V. High frequency restoration of surface seismic data. Lead. Edge 2003, 22, 730–738. [Google Scholar] [CrossRef]
  7. van der Baan, M.A. Bandwidth enhancement—Inverse Q deconvolution or time-varying wiener filtering? In Proceedings of the 74th EAGE Conference and Exhibition Incorporating Europec 2012, Copenhagen, Denmark, 4–7 June 2012. [CrossRef]
  8. Braga, I.L.S.; Moraes, F.S. High-resolution gathers by inverse Q filtering in the wavelet domain. Geophysics 2013, 78, V53–V61. [Google Scholar] [CrossRef]
  9. Ren, H.; Wang, H.; Zhang, L. Compensation for absorption and attenuation using wave equation continuation along ray path. Geophys. Prospect. Pet. 2007, 46, 557–563. [Google Scholar] [CrossRef]
  10. Li, G.F.; Liu, Y.; Zheng, H.; Huang, W. Absorption decomposition and compensation via a two-step scheme. Geophysics 2015, 80, V145–V155. [Google Scholar] [CrossRef]
  11. Mittet, R.; Sollie, R.; Hokstad, K. Prestack depth migration with compensation for absorption and dispersion. Geophysics 1995, 60, 1485–1494. [Google Scholar] [CrossRef]
  12. Wang, H. Seismic wave imaging in visco-acoustic media. Sci. China Ser. A 2004, 47, 146–154. [Google Scholar] [CrossRef]
  13. Suh, S.; Yoon, K.; Cai, J.; Wang, B. Compensating visco-acoustic effects in anisotropic reverse-time migration. In SEG Technical Program Expanded Abstracts 2012; Society of Exploration Geophysicists: Houston, TX, USA, 2012; pp. 1–5. [Google Scholar] [CrossRef]
  14. Dutta, G.; Schuster, G.T. Attenuation compensation for least-squares reverse time migration using the viscoacoustic-wave equation. Geophysics 2014, 79, S251–S262. [Google Scholar] [CrossRef]
  15. Zhu, T.; Harris, J.M.; Biondi, B. Q-compensated reverse-time migration. Geophysics 2014, 79, S77–S87. [Google Scholar] [CrossRef]
  16. Chen, H.M.; Zhou, H.; Rao, Y. An implicit stabilization strategy for Q-compensated reverse time migration. Geophysics 2020, 85, S169–S183. [Google Scholar] [CrossRef]
  17. Yang, J.D.; Huang, J.P.; Zhu, H.J.; Li, Z.C.; Dai, N.X. Viscoacoustic reverse time migration with a robust space-wavenumber domain attenuation compensation operator. Geophysics 2021, 86, S339–S353. [Google Scholar] [CrossRef]
  18. Fathalian, A.; Trad, D.O.; Innanen, K.A. Q-compensated reverse time migration in tilted transversely isotropic media. Geophysics 2021, 86, S73–S89. [Google Scholar] [CrossRef]
  19. Sun, J.; Zhu, T. Strategies for stable attenuation compensation in reverse-time migration. Geophys. Prospect. 2017, 66, 498–511. [Google Scholar] [CrossRef]
  20. Zhang, B.; Wang, H. One-way wave propagation in the ray-centred coordinate system for vertical transversely isotropic media. Acta Geophys. 2021, 69, 1225–1239. [Google Scholar] [CrossRef]
  21. O’Doherty, R.F.; Anstey, N.A. Reflections on amplitudes. Geophys. Prospect. 1971, 19, 430–458. [Google Scholar] [CrossRef]
  22. Cheng, L.; Wang, H.; Liu, S. Simulation and imaging with paraxial one-way wave equation in ray-centered coordinates. Geophys. Prospect. Pet. 2012, 51, 338–342. [Google Scholar] [CrossRef]
  23. Zhang, L.; Wang, H. A stable inverse Q migration method. Geophys. Prospect. Pet. 2010, 49, 115–120. [Google Scholar] [CrossRef]
  24. Ma, Z.T. Seismic Imaging Techniques-Finite Difference Migration; Petroleum Industry Press: Beijing, China, 1989; p. 84. [Google Scholar]
  25. Bickel, S.H.; Natarajan, R.R. Plane-wave Q deconvolution. Geophysics 1985, 50, 1426–1439. [Google Scholar] [CrossRef]
Figure 1. Schematic of the two-way path compensation region (blue region).
Figure 1. Schematic of the two-way path compensation region (blue region).
Energies 15 04302 g001
Figure 2. Depression velocity model.
Figure 2. Depression velocity model.
Energies 15 04302 g002
Figure 3. Comparison of imaging results: (a) acoustic migration result of acoustic data; (b) acoustic migration result from visco-acoustic data; (c) visco-acoustic migration result of visco-acoustic data.
Figure 3. Comparison of imaging results: (a) acoustic migration result of acoustic data; (b) acoustic migration result from visco-acoustic data; (c) visco-acoustic migration result of visco-acoustic data.
Energies 15 04302 g003
Figure 4. Comparison of imaging results at different CDP locations: (a) 120 CDP (i.e., the blue line on the left of each figure in Figure 3); (b) 200 CDP (i.e., the red line on the right of each figure in Figure 3). Each trace NO. 1 of Figure 4a,b represents the acoustic migration result of acoustic data (ideal result); each trace NO. 2 represents the acoustic migration result of visco-acoustic data; and each trace NO. 3 represents the visco-acoustic migration result in visco-acoustic data.
Figure 4. Comparison of imaging results at different CDP locations: (a) 120 CDP (i.e., the blue line on the left of each figure in Figure 3); (b) 200 CDP (i.e., the red line on the right of each figure in Figure 3). Each trace NO. 1 of Figure 4a,b represents the acoustic migration result of acoustic data (ideal result); each trace NO. 2 represents the acoustic migration result of visco-acoustic data; and each trace NO. 3 represents the visco-acoustic migration result in visco-acoustic data.
Energies 15 04302 g004
Figure 5. Bottom splitting model: (a) velocity; (b) Q model.
Figure 5. Bottom splitting model: (a) velocity; (b) Q model.
Energies 15 04302 g005
Figure 6. Comparison of imaging results: (a) acoustic migration result from visco-acoustic data; (b) visco-acoustic migration result of visco-acoustic data.
Figure 6. Comparison of imaging results: (a) acoustic migration result from visco-acoustic data; (b) visco-acoustic migration result of visco-acoustic data.
Energies 15 04302 g006
Figure 7. Modified Hess model: (a) velocity; (b) Q model.
Figure 7. Modified Hess model: (a) velocity; (b) Q model.
Energies 15 04302 g007
Figure 8. Comparison of imaging results: (a) acoustic migration result from visco-acoustic data; (b) visco-acoustic migration result of visco-acoustic data.
Figure 8. Comparison of imaging results: (a) acoustic migration result from visco-acoustic data; (b) visco-acoustic migration result of visco-acoustic data.
Energies 15 04302 g008
Figure 9. Comparison of imaging results: (a) acoustic migration result; (b) visco-acoustic migration result.
Figure 9. Comparison of imaging results: (a) acoustic migration result; (b) visco-acoustic migration result.
Energies 15 04302 g009
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, B.; Wang, H. Finite Difference Pre-Stack Depth Migration of One-Way Wave Equation in Isotropic Visco-Acoustic Media of a Ray-Centered Coordinate System. Energies 2022, 15, 4302. https://doi.org/10.3390/en15124302

AMA Style

Zhang B, Wang H. Finite Difference Pre-Stack Depth Migration of One-Way Wave Equation in Isotropic Visco-Acoustic Media of a Ray-Centered Coordinate System. Energies. 2022; 15(12):4302. https://doi.org/10.3390/en15124302

Chicago/Turabian Style

Zhang, Bohan, and Huazhong Wang. 2022. "Finite Difference Pre-Stack Depth Migration of One-Way Wave Equation in Isotropic Visco-Acoustic Media of a Ray-Centered Coordinate System" Energies 15, no. 12: 4302. https://doi.org/10.3390/en15124302

APA Style

Zhang, B., & Wang, H. (2022). Finite Difference Pre-Stack Depth Migration of One-Way Wave Equation in Isotropic Visco-Acoustic Media of a Ray-Centered Coordinate System. Energies, 15(12), 4302. https://doi.org/10.3390/en15124302

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop