Jing Wang
School of Statistics and Mathematics, Hubei University of Economics, Wuhan, China.
DOI: 10.4236/gep.2026.145019   PDF    HTML   XML   3 Downloads   29 Views   Citations

Abstract

Elastic reverse time migration is a vector wave theory-based depth domain migration algorithm. Wavefields separation is a key step to remove crosstalk artifacts and improve the imaging quality. This operation allows us to produce vector P- and S-wavefields with the same phases and amplitudes as the input coupled wavefields, while significantly reducing computational costs. The test results illustrate that the amplitude and phase of the separated wavefields are changed and the polarity is reversed after the Helmholtz decomposition, so that further correction is needed. Aiming to solve these problems, decoupled elastic wave equations are utilized to simulate the wavefields propagation, which can maintain the vector property of original wavefields without polarity reversal. With the separated vector wavefields, we implement a modified dot-product imaging condition for elastic reverse time migration, which produces migrated images with accurate amplitudes. Several numerical examples are used to demonstrate its feasibility and robustness for imaging complex subsurface structures. In this paper, we firstly review the Helmholtz decomposition and analyze its shortcomings. Next, we utilize the decoupled elastic wave equations to obtain vector P- and S-wavefields with the same phases and amplitudes as the input coupled wavefields. And then, a dot-product imaging condition and corresponding elastic reverse time migration workflow are presented to produce PP- and PS-images. Finally, several numerical examples are used to illustrate the performance of this wavefield decomposition strategy and the corresponding workflow. Besides, the numerical simulation methods for the forward and backward wavefields’ extrapolation, the effects of multiple waves are also discussed in this report to further improve the imaging efficiency and accuracy.

Keywords

Reverse Time Migration, Decoupled Elastic Wave Equations, Imaging Efficiency and Accuracy

Share and Cite:

1. Introduction

Elastic reverse time migration (RTM) is a powerful tool for estimating lithology information and imaging subsurface structures by taking advantages of P- and S-waves. The S-wave images always have a higher resolution than those of the P-wave due to its shorter wavelength. In addition, the interpretation of P- and S-wave images has great potential to detect sweet spots in unconventional tight-sandstone reservoirs (Tang et al., 2009). Elastic RTM reconstructs forward and backward vector wavefields in the subsurface with multicomponent records as boundary conditions, which can be further used to produce accurate PP- and PS-reflectivity when proper imaging conditions are applied. One important step in elastic RTM is to extract pure-mode P- and S-waves from the coupled elastic wavefields to remove crosstalk artifacts and improve imaging quality. The P- and S-waves separation in isotropic media can be implemented by numerically solving P- and S-wave separated elastic wave equations (Li et al., 2007; Zhang et al., 2007) or wavenumber-domain wavefield vector decomposition (Zhang & McMechan, 2010). Yan and Sava (2008) implement the wavefields separation to compute the scalar and vector potential wavefields using divergence and curl operators. However, divergence and curl operators resulting in the scalar and vector potential wavefields do not have the same phases, amplitudes, and physical units as the extrapolated wavefields, migration results are not accurate. Besides, this wavefields separation approach leads to a polarity reversal problem for PS-imaging profiles. Without proper angle corrections, summing PS-images over different sources leads to nonconstructive stacking results (Du et al., 2012; Duan & Sava, 2015; Wang et al., 2016). In order to avoid using the divergence and curl operators, a vector wavefield decomposition approach (Zhang & McMechan, 2010) by solving a linear system in the wavenumber domain is developed. This method can maintain vector P- and S-wavefields that have the same phases, amplitudes and units as the input elastic wavefields. And then, Zhu (2017) proposes an equivalent strategy in the time-space domain, which requires solving a vector Poisson’s equation. But the computational cost is still prohibitive for 3D problems. Another way to implement the vector wavefield decomposition is to introduce an auxiliary P-wave equation (Xiao & Leaney, 2010; Wang & McMechan, 2015). The S-waves can be obtained by subtracting the P-waves from the coupled input wavefields. Although this strategy can maintain the separation results as accurately as the vector wavefield decomposition approach methods, solving an auxiliary P-wave equation requires additional computational costs. The other way is to separate the wavefields while the propagation process. Yang et al. develop an efficient vector wavefields decomposition strategy and implement a modified dot-product imaging condition for elastic RTM. This imaging condition retains the signs of the dot product but recomputes the amplitudes using the multiplication of the absolute values of the separated source and receiver wavefields.

This report is organized as follows. First, we review the Helmholtz decomposition and analyze its shortcomings. Next, we utilize the decoupled elastic wave equations to obtain vector P- and S-wavefields with the same phases and amplitudes as the input coupled wavefields. And then, a dot-product imaging condition and corresponding elastic RTM workflow are presented to produce PP- and PS-images. Finally, several numerical examples are used to illustrate the performance of the proposed wavefield decomposition-based elastic RTM workflow. Besides, the other strategies to further improve the elastic RTM imaging quality and efficiency are also discussed.

2. Method

2.1. Helmholtz Decomposition-Based Vector Wavefield Separation

Based on the Helmholtz decomposition theory in isotropic elastic media, the vector displacement or particle velocity wavefields $u$ can be decomposed as (Morse & Feshbach, 1946; Aki & Richards, 2002)

$u=u^{\text{P}}+u^{\text{S}}=\nabla \phi +\nabla \Psi$ , (1)

where $u^P$ and $u^S$ denote the curl-free P-wavefield and divergence-free S-wavefield, respectively. $\phi$ and $\Psi$ are their corresponding scalar- and vector-potential wavefields.

Note that the P- and S-waves are coupled together during wave propagation. By adopting vector Poisson’s formular ( ${\nabla }^{2}w=u$ ) and the property of the Laplacian operator ( ${\nabla }^{2}w=\nabla \left(\nabla \cdot w\right)-\nabla \times \nabla \times w$ ), isotropic decoupled P- and S-wave formulations are

$u^{\text{P}}=\nabla \left(\nabla \cdot w\right)$ , $u^{\text{S}}=-\nabla \times \nabla \times w$ , (2)

in which $\nabla$ , $\nabla \cdot$ and $\nabla \times$ are the gradient, divergence and curl operators, respectively. At every time step of each wavefield extrapolation, the forward and backward multidimensional Fourier transforms had to be calculated. Under current computational capability, their costs are still prohibitive, especially for large-scale 3D problems. Note that the P- and S-waves are coupled together during wave propagation.

2.2. Decoupled Elastic Wave Equations-Based Elastic RTM

The P- and S-wave separated elastic wave formulae are

$\rho \frac{\partial {\tau }_{xx}^{\text{p}}}{\partial t}=\left(\lambda +2\mu \right)\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z}\right)$ , $\rho \frac{\partial {\tau }_{xx}^{\text{s}}}{\partial t}=-2\mu \frac{\partial v_z}{\partial z}$ ,

$\rho \frac{\partial {\tau }_{zz}^{\text{p}}}{\partial t}=\left(\lambda +2\mu \right)\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z}\right)$ , $\rho \frac{\partial {\tau }_{zz}^{\text{s}}}{\partial t}=-2\mu \frac{\partial v_z}{\partial x}$ ,

$\rho \frac{\partial {\tau }_{xz}^{\text{s}}}{\partial t}=\mu \left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)$ , (3)

$\rho \frac{\partial v_x^{\text{p}}}{\partial t}=\frac{\partial {\tau }_{xx}^{\text{p}}}{\partial x}$ , $\rho \frac{\partial v_x^{\text{s}}}{\partial t}=\frac{\partial {\tau }_{xx}^{\text{s}}}{\partial x}+\frac{\partial {\tau }_{xz}^{\text{s}}}{\partial z}$ ,

$\rho \frac{\partial v_z^{\text{p}}}{\partial t}=\frac{\partial {\tau }_{zz}^{\text{p}}}{\partial z}$ , $\rho \frac{\partial v_z^{\text{s}}}{\partial t}=\frac{\partial {\tau }_{xz}^{\text{s}}}{\partial x}+\frac{\partial {\tau }_{zz}^{\text{s}}}{\partial z}$ ,

$v_x=v_x^{\text{p}}+v_x^{\text{s}}$ , $v_z=v_z^{\text{p}}+v_z^{\text{s}}$ ,

where superscripts p or s stand for the corresponding P- or S-wavefield components. After the separated wavefields are given, we can simultaneously solve the P- and S-wavefields with high accuracy through substituting P-wave coefficients into P-wavefield formulae, and utilizing S-wave coefficients to calculate the S-wavefield formulae, and then the combined velocity terms v_x and v_z are obtained. Figure 1 illustrates an experiment of the vector wavefield decomposition using Equation (3) for a homogenous model. The calculation area is 4.0 km * 4.0 km with grid spacing 18 m. The time sampling interval is 1.5 ms. A ricker wavelet with 20 Hz peak frequency is chosen as the source and excited at the centre of the model. Figure 1 shows that errors between the input and reconstructed wavefields are almost zero both in the horizontal and vertical directions, and the amplitudes, phases, and units of the separated wavefields are the same as the input wavefields. So that we can conclude that we can obtain the accurate P- and S-wavefields when adopt the decoupled elastic wave equations to do the forward and backward wave extrapolation, which can improve the accuracy of the elastic RTM imaging profile.

2.3. A Dot-Product Imaging Condition (DP IMC) for Vector-Based Elastic RTM

The polarity reversal problem for PS-images, a well-known problem in traditional elastic RTM (Rosales et al., 2008; Du et al., 2012; Duan & Sava, 2015). However, this imaging condition involves multiplication by a cosine function, where $\Delta \theta$ denotes the included polarization angle between the incident and reflected wavefields. This leads to an inaccurate estimation of reflection coefficients. To eliminate the amplitude effects of the cosine function, Du et al. (2017) suggest normalizing the dot-product imaging results with the absolute value of $\cos \left(\Delta \theta \right)$ . When multiple waves intersect in the areas with complicated structures, it is difficult to accurately estimate the propagation and polarity directions (Tang & McMechan, 2016). To simplify this angle-dependent normalization, we retain the signs of the dot product and recompute the amplitudes using the multiplication of the absolute values of the separated source and receiver wavefields. This gives us the following MDP IMC:

$I^{PP}\left(x,y,z\right)=\frac{{\displaystyle \sum _{t}sig{n}^{\text{PP}}\cdot \left[\left|v_p^{src}\left(x,y,z,t\right)\right|\cdot \left|v_p^{rcv}\left(x,y,z,t\right)\right|\right]}}{{\displaystyle \sum _t\left[v_p^{src}\left(x,y,z,t\right)\cdot v_p^{src}\left(x,y,z,t\right)\right]}}$ , (4)

$I^{PS}\left(x,y,z\right)=\frac{{\displaystyle \sum _{t}sig{n}^{\text{PS}}\cdot \left[\left|v_p^{src}\left(x,y,z,t\right)\right|\cdot \left|v_s^{rcv}\left(x,y,z,t\right)\right|\right]}}{{\displaystyle \sum _t\left[v_p^{src}\left(x,y,z,t\right)\cdot v_p^{src}\left(x,y,z,t\right)\right]}}$ , (5)

$sig{n}^{\text{PP}}=\text{sign}\left[v_p^{src}\left(x,y,z,t\right)\cdot v_p^{rcv}\left(x,y,z,t\right)\right]$ , (6)

$sig{n}^{\text{PS}}=\text{sign}\left[v_p^{src}\left(x,y,z,t\right)\cdot v_s^{rcv}\left(x,y,z,t\right)\right]$ . (7)

The Helmholtz decomposition and the corresponding correction methods are tested for a three-layers homogeneous model. The calculation area is 4.0 km and 2.0 km in the horizontal and vertical directions, respectively, with grid spacing 10 m. The time sampling interval is 1 ms. A ricker wavelet with 23 Hz peak frequency is chosen as the source. There are 120 geophones excited at the surface of the model. Seismic records are received by 200 receivers evenly arranged on the surface.

In Figure 2, we observe that on both sides of the vertical incident position of the P-wave, the event polarities are opposite and discontinuous, also the energy gradually diminishes. The gradual weakened energy is due to the fact that when the P-wave is vertically incident, almost no converted wavefield is generated, and the polarities of the converted wave are opposite on both sides of the vertical incidence P-wave, resulting in the discontinuous of the event. Based on the dot-product imaging condition, we can observe from Figure 1(c) that there is no polarity reversal, the event continuous natural, and the imaging profile quality is greatly improved.

3. Numerical Examples

In this section, we adopt the Marmousi model (P-wave velocity is depicted in Figure 3) to illustrate the performance of the different elastic RTM methods. The calculation area is 5 km and 3.52 km in the vertical and horizontal directions,

respectively, with grid spacing 10 m. The time sampling interval is 1 ms. An explosive source of Ricker wavelet with the peak frequency of 27 Hz is located 9.5 km horizontally and 10 m vertically for the P- and S-wavefields decomposition test. Figure 4 depicts the horizontal and vertical snapshots components at 1.5 s calculated by temporal high-order and spatial implicit stagged-gird finite-difference schemes (Wang et al. 2022) with 18^th-order accuracy in space and 6^th-order accuracy in time, respectively. One can observe that there are obvious errors on the snapshot calculated by the conventional Helmholtz decomposition method. Although the P- and S-wavefields are separated more clearly in Figure 4(b), Helmholtz decomposition with polarity correction still has visible error. When

we utilize the decoupled elastic wave equations to perform the P-and S-wavefields propagation, the error is zero, which means we can obtain the accuracy P- and S-waves during the wavefields forward and backward process of the elastic RTM to gain the most satisfactory imaging performance.

Next, we still use the truncated Marmousi model to further test the PP- and PS-imaging performance using different elastic RTM workflows. The space and time sampling interval are 10 m and 1 ms, respectively. 98 explosive sources are arranged on the 10 m vertically line and begin at 7.5 km in horizontally direction with 50 m interval. 40 receivers are arranged in the well (the vertical black line shown in Figure 5(a)) with 30 m interval. Figure 5(a) and Figure 5(b) are the PP imaging results based on the Helmholtz decomposition method and the decoupled elastic wave equation, respectively. Both imaging results can clearly describe the complex structure of the Marmousi model. The wavefields separation based on Helmholtz decomposition applies a divergence operator to the wavefields of the source and the receiver, resulting in a π/2 phase shift of the separated P-wavefield. Figure 5(c) is the PS imaging results based on Helmholtz decomposition, Figure 5(e) is the PS imaging result based on the wavefields separation of the decoupled elastic wave equation. Comparing the two methods, we can observe that the former has discontinuous events and occurs polarity reversal phenomenon, resulting in poor imaging quality, while the latter has clear and continuous events, no polarity reversal, clear imaging of complex structures. Besides, the resolution of the PS imaging profile is higher than that of the PP imaging profile. Figure 5(d) shows the imaging results of the PS imaging condition with polarity correction. Compared with Figure 5(c), the imaging quality is greatly improved.

From the above analysis, we can conclude that the PS imaging results based on

Helmholtz decomposition resulting in a phase shift problem, because the wavefields of the curl operator also has a phase shift of π/2, which will cause a phase shift of π after cross-correlation. The elastic RTM imaging based on the decoupled elastic wave equation wavefields separation method can maintain the vector characteristics of original wavefields without the phenomenon of polarity reversal, so the imaging performance is more ideal.

4. Discussion

In this section, we will discuss some strategies to further improve the elastic RTM imaging accuracy and efficiency.

4.1. Variable Length Temporal and Spatial Operator Finite-Difference Method

In order to decrease the forward and backward wavefields extrapolation calculation cost, we develop a finite difference strategy with variable lengths of temporal and spatial operators (Zhou et al., 2018). For a given velocity value, spatial and temporal discretization interval, the minimum lengths of these operators can be automatically determined. Compared with conventional fixed length operator method or variable spatial operator length method, our approach has superior modeling efficiency under the same accuracy demand. Here we provide procedures for fast search strategy.

Step 1. For each velocity value in the complex velocity model, set spatial operator length parameter M = N (temporal accuracy length parameter) first, then search for minimum $M_{\min }$ so that $\epsilon \left(v,M,N,f\right)<\eta$ .

Step 2. Gradually decrease N to a minimum value $N_{\min }$ that satisfies the accuracy constrain $-\eta <\epsilon \left(v,M,N,f\right)$ and stability constrain.

Step 3. Using searched $M_{\min }$ and $N_{\min }$ to perform wave extrapolation for each local velocity.

In practice implementation, one needs only to calculate the temporal of spatial operator lengths with a velocity interval of 1 m/s and store corresponding FD coefficients before wave extrapolation. We use the 2D salt model (Figure 6(a)) to compare efficiency and accuracy of fixed length temporal and spatial operators (FLTSO), variable length spatial operator (VLSO), and simultaneously variable length temporal and spatial operator (VLTSO) strategies. For paper concise, FD coefficients are obtained only by TEM. Figure 7 further shows the variation of temporal and spatial operators’ lengths with velocity in complex velocity example. With the increase of velocity, both temporal and spatial operators’ lengths of VLTSO method are shorter than those of VLSO method. It is, therefore, fundamentally explained that VLTSO method is always more efficient (as shown in Table 1) than VLSO under the same accuracy demand.

4.2. Primaries and First-Order Multiples Combination-Based RTM

The conventional seismic migration method employs primaries to perform

migration imaging only and the multiples, regarded as a kind of noise, has to be removed before imaging. However, the migration of multiples has a higher illuminance than that of the primaries, which can improve the accuracy and illumination of seismic imaging profile. Recently, many imaging methods that use free surface waves have been developed, but most of them need to predict multiples first. Such methods are more time-consuming, also it becomes difficult to accurately predict multiples when the structure is complicated. In this report, we combine the primaries and first-order multiples together to implement RTM, which can image directly without predicting and separating the primaries and multiples.

In the process of seismic data acquisition, the received seismic wavefield data not only contains the primary, but also contains the free surface multiples, which are expressed as:

$R\left(z_0\right)=P\left(z_0\right)+M\left(z_0\right)+\mathrm{...}$ , (7)

in which, is the receiver wavefields, $P\left(z_0\right)$ is the primary reflected wave, $M\left(z_0\right)$ is the free surface multiple, and $z_0$ is the depth of the receiver point.

The basic idea of multiple reverse time migration is to regard the received wavefield data as a virtual source firstly, which can be utilized to obtain forward propagation data with the migration velocity model, and then combine the synthesized forward propagation data with the received free surface multiple to implement the migration according to the proper imaging condition. But this method needs to predict multiples in advance, which increases the computational cost. Through the analysis, it can be seen that the source wavelet and the virtual source in the forward modeling process can be regarded as a new source. Based on that, the forward propagation data can be expressed as:

$P_F\left(z_0\right)+M_F\left(z_0\right)=-XW\left(z_0\right)-XR\left(z_0\right)$ , (8)

in which, $W\left(z_0\right)$ is the source function, $X$ is the migration velocity model.

Figure 8 shows the imaging profiles calculated the primaries, and the primaries and first-order multiples combination-based RTM, respectively. According to results, the primaries and first-order multiples combination-based RTM can not only broaden the width of the shallow imaging, but also have an advantage in imaging the deep structure of complex models.

5. Conclusion

Accuracy and efficiency are two important issues to be considered when implementing the RTM. In order to improve the imaging quality and reduce the time consumption, we have made a lot of attempts focusing on the numerical simulation of wavefields extrapolation, P- and S-wavefield separation methods, imaging condition, and take the advantage of the multiples. And then, we adopt different synthetic models to demonstrate the effects and imaging performance of these methods. In further works, we will further explore the application of these aspects for the field walkway VSP data.

Acknowledgements

This research is supported by the Youth Fund of Hubei Province Natural Science Foundation (Grant No. 2024AFB348) and the Research Foundation of Hubei University of Economics (Grant No. XJ23BS45). We would like to thank Prof. Kristopher A. Innanen and CREWES for their continued support and valuable comments.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Aki, K., & Richards, P. G. (2002). Quantitative Seismology (2nd ed.). University Science Books.
[2]	Du, Q., Guo, C., Zhao, Q., Gong, X., Wang, C., & Li, X. (2017). Vector-Based Elastic Reverse Time Migration Based on Scalar Imaging Condition. Geophysics, 82, S111-S127. [Google Scholar] [CrossRef]
[3]	Du, Q., Zhu, Y., & Ba, J. (2012). Polarity Reversal Correction for Elastic Reverse Time Migration. Geophysics, 77, S31-S41. [Google Scholar] [CrossRef]
[4]	Duan, Y., & Sava, P. (2015). Scalar Imaging Condition for Elastic Reverse Time Migration. Geophysics, 80, S127-S136. [Google Scholar] [CrossRef]
[5]	Jianming, T., Yue, H., Xiangrong, X., Tinnin, J., & Hallin, J. (2009). Application of Converted-Wave 3D/3-C Data for Fracture Detection in a Deep Tight-Gas Reservoir. The Leading Edge, 28, 826-837. [Google Scholar] [CrossRef]
[6]	Li, Z., Zhang, H., Liu, Q., & Han, W. (2007). Numerical Simulation of Elastic Wavefield Separation by High-Order Staggered-Grid Finite-Difference Method. Oil Geophysical Prospecting, 42, 510-515. (In Chinese)
[7]	Morse, P. M., & Feshbach, H. (1946). Methods of Theoretical Physics. Technology Press.
[8]	Rosales, D. A., Fomel, S., Biondi, B. L., & Sava, P. C. (2008). Wave-Equation Angle-Domain Common-Image Gathers for Converted Waves. Geophysics, 73, S17-S26. [Google Scholar] [CrossRef]
[9]	Tang, C., & McMechan, G. A. (2016). Multidirectional Slowness Vector for Computing Angle Gathers from Reverse Time Migration. Geophysics, 81, S55-S68. [Google Scholar] [CrossRef]
[10]	Wang, C., Cheng, J., & Arntsen, B. (2016). Scalar and Vector Imaging Based on Wave Mode Decoupling for Elastic Reverse Time Migration in Isotropic and Transversely Isotropic Media. Geophysics, 81, S383-S398. [Google Scholar] [CrossRef]
[11]	Wang, J., Liu, Y., & Zhou, H. (2022). High Temporal Accuracy Elastic Wave Simulation with New Time-Space Domain Implicit Staggered-Grid Finite-difference Schemes. Geophysical Prospecting, 70, 1346-1366. [Google Scholar] [CrossRef]
[12]	Wang, W., & McMechan, G. A. (2015). Vector-Based Elastic Reverse Time Migration. Geophysics, 80, S245-S258. [Google Scholar] [CrossRef]
[13]	Xiao, X., & Leaney, W. S. (2010). Local Vertical Seismic Profiling (VSP) Elastic Reverse-Time Migration and Migration Resolution: Salt-Flank Imaging with Transmitted Pto-S Waves. Geophysics, 75, S35-S49. [Google Scholar] [CrossRef]
[14]	Yan, J., & Sava, P. (2008). Isotropic Angle-Domain Elastic Reverse-Time Migration. Geophysics, 73, S229-S239. [Google Scholar] [CrossRef]
[15]	Zhang, J., Tian, Z., & Wang, C. (2007). P- and S-Wave Separated Elastic Wave Equation Numerical Modeling Using 2D Staggered-Grid. In 77th Annual International Meeting, SEG, Expanded Abstracts (pp. 2014-2019). Society of Exploration Geophysics
[16]	Zhang, Q., & McMechan, G. A. (2010). 2D and 3D Elastic Wavefield Vector Decomposition in the Wavenumber Domain for VTI Media. Geophysics, 75, D13-D26. [Google Scholar] [CrossRef]
[17]	Zhou, H., Liu, Y., & Wang, J. (2018). Finite-Difference Modeling with Adaptive Variable-Length Temporal and Spatial Operators. In SEG Technical Program Expanded Abstracts 2018 (pp. 4015-4019). Society of Exploration Geophysicists. [Google Scholar] [CrossRef]
[18]	Zhu, H. (2017). Elastic Wave-field Separation Based on the Helmholtz Decomposition. Geophysics, 82, S173-S183. [Google Scholar] [CrossRef]