Elastic Reverse Time Migration for Weakly Illuminated Structure

Abstract

One of the most effective techniques to obtain PP and PS images is elastic reverse time migration which employs multi-component seismic data. The two types of complementary images play an important role in reducing blind spots in seismic exploration. However, the migration image of deep structures is always blurred due to the shielding effect of overburden rock on seismic waves. To overcome this issue, we develop an elastic reverse time migration approach for insufficient illumination. This approach contains two crucial elements. The first is that we derive an elastic wave equation to extract the wavefields associated with the exploration target using the staining algorithm. Secondly, we develop an inner product imaging condition with a filter to mute migrated artifacts. The filter, consisting of two vectors, determines which part of the wavefield is contributed to imaging. Synthetic examples exhibit that the proposed elastic reverse time migration method can improve the signal-to-noise ratio of PP and PS images of weakly illuminated structures.

1. Introduction

Elastic reverse time migration (ERTM) possesses the ability to reposition the multi-component seismic data into underground structure information, thus providing a scientific basis for further seismic interpretation and hydrocarbon development [1,2]. Particularly, subsalt has become a significant hydrocarbon exploration target due to the appearance of salt beds as regional overburden rock [3]. However, owing to the shielding effect generated by the complex overburden rock on seismic signals, the imaging results of the subsurface geological structure are often unsatisfactorily [4,5]. Therefore, it is essential to apply the elastic migration technique for insufficient illumination areas. The wide-azimuth measurement, illumination compensation and acquisition aperture correction approaches have been developed to enhance the imaging quality of shadow regions [6,7,8]. In addition, the staining algorithm can highlight the wavefields associated with weakly illuminated structures among the full wavefields [9,10,11], which is a powerful means to realize target-oriented imaging within the ERTM framework.

ERTM images produced by multi-component seismic data always involve crosstalks due to the reason that P- and S-waves are concurrently excited when the waves propagate to a reflecting interface [12]. To suppress the crosstalk artifacts, the elastic wavefield should be adequately decoupled before imaging. A common approach is to apply the curl and divergence operators to the coupled wavefield [13]. However, the decoupled S-wave presenting the polarity reversal phenomenon ultimately destroys the continuity of the reflection lineups of multi-shot stacked profiles. Du et al. [14] argued that the polarity of the decoupled S-wave is reversed on both sides of the vertical incident direction of the P-wave. With the assistance of the Poynting vector, a sign factor can be estimated to correct the polarity of PS imaging [15]. However, the divergence and curl operators will change the amplitude and phase of wavefields. As a consequence, it is difficult for ERTM to obtain a definite physical meaning from the converted wave imaging profile. To address this problem, the pure P- and S-wave can be effortlessly decoupled by solving the vector-decoupled wave equation [16].

When applying the correlation imaging condition, the noises inevitably accompanying valuable information will debase the quality of migration profiles, especially on the interface where the velocity varies dramatically. Although the high-pass filter can suppress migration noises, it alters the waveform of the original image [17]. Another image denoising approach whereby the high angle information of the common imaging gathers in the angle domain is eliminated has advantages in amplitude-preserving [18,19]. Liu et al. [20] first suggested that the migration noises are related to the wave propagation path through the scattering point, and finally proposed an imaging condition with wavefield decomposition to remove the noises. Unfortunately, the initial wavefield decomposition method is realized in the wavenumber-frequency domain, which may result in a huge burden in computation, storage and input/output. Fei et al. [21] imported a cost-effective scheme to decompose the wavefield based on the Hilbert transform. Further, Wang et al. [22] decomposed the wavefield into several components according to the propagating directions. Particularly, the Poynting vector can indicate the direction of energy transfer, thereby flexibly realizing wavefield decomposition [23].

The paper is arranged as follows: after the introduction section, the vector-decoupled wave equation is reviewed. On this basis, a novel equation is derived to extract the wavefield associated with insufficient illumination. Next, an inner product imaging condition with filters is described. Finally, some synthetic examples are employed to demonstrate the effectiveness of our approach.

2. Methodology

2.1. Vector-Decoupled Wave Equation

The velocity-stress formula in 2D isotropic medium is expressed as,

$$\{\begin{array}{c}\rho \frac{\partial \mathit{V}}{\partial t}=ℒ\mathbf{\sigma }\\ \frac{\partial \mathbf{\sigma }}{\partial t}=\rho \mathit{A}{ℒ}^T\mathit{V}\end{array}$$

(1)

where $\rho$ denotes the density of media, $\mathit{V}={\left[v_x,v_z\right]}^T$ is the particle-velocity vector, $t$ signifies the time variable, $\mathbf{\sigma }={\left[{\tau }_{xx},{\tau }_{zz},{\tau }_{xz}\right]}^T$ denotes the stress vector and $ℒ$ is a spatial partial derivative operator expressed as:

$$ℒ=\left[\begin{array}{ccc}{\partial }_x& 0& {\partial }_z\\ 0& {\partial }_z& {\partial }_x\end{array}\right]$$

(2)

$x$ and $z$, respectively represent the horizontal and vertical directions, and $\mathit{A}$ is a parameter matrix as shown below:

$$\mathit{A}=\left[\begin{array}{ccc}{c}_P^2& c_P^2-2c_S^2& 0\\ c_P^2-2c_S^2& c_P^2& 0\\ 0& 0& c_S^2\end{array}\right]$$

(3)

where $c_P$ and $c_S$ are P-wave and S-wave velocity in the propagation medium. Equation (1) can be named as vector-coupled wave equation (VCWE).

When the wavefield is simulated according to Equation (1), the P- and S-waves are coupled in the wavefield of particle velocity. If the coupled wavefield is directly applied to the ERTM, the imaging profile will suffer from crosstalk. To tackle this problem, Wang et al. [24] introduced a scalar P-wave stress ${\tau }^P$ to decouple the wavefield, and the formulation is written as,

$$\{\begin{array}{c}\frac{\partial {\tau }^P}{\partial t}=\rho c_P^2\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_z}{\partial z}\right)\\ \rho \frac{\partial v_x^P}{\partial t}=\frac{\partial {\tau }^P}{\partial x}\\ \rho \frac{\partial v_z^P}{\partial t}=\frac{\partial {\tau }^P}{\partial z}\end{array}$$

(4)

where $v_x^P$ and $v_z^P$ are the decoupled P-waves in $x$- and $z$-directions. The residual S-waves are defined as $v_x^S=v_x-v_x^P$ and $v_z^S=v_z-v_z^P$. The decoupled wavefield can preserve the same physical properties as the coupled wavefield. The combination of Equations (1) and (4) is the vector-decoupled wave equation (VDWE).

2.2. Stained VDWE

To extract the wavefield associated with the weakly illuminated structure, we establish a governing equation. According to the principle of staining algorithm [9], the variables $\mathit{V}$, $\mathbf{\sigma }$ and $\mathit{A}$ in VDWE are extended to complex domain, and possess the following forms:

$$\mathit{V}=\overline{\mathit{V}}+i\tilde{\mathit{V}}={\left[\overline{v_x},\overline{v_z}\right]}^T+i{\left[\widetilde{v_x},\widetilde{v_z}\right]}^T,$$

(5)

$$\mathbf{\sigma }=\overline{\mathbf{\sigma }}+i\tilde{\mathbf{\sigma }}={\left[\overline{{\tau }_{xx}},\overline{{\tau }_{zz}},\overline{{\tau }_{xz}}\right]}^T+i{\left[\widetilde{{\tau }_{xx}},\widetilde{{\tau }_{zz}},\widetilde{{\tau }_{xz}}\right]}^T,$$

(6)

$$\mathit{A}=\overline{\mathit{A}}+i\tilde{\mathit{A}}=\left[\begin{array}{ccc}{\left(\overline{c_P}+i\widetilde{c_P}\right)}^2& {\left(\overline{c_P}+i\widetilde{c_P}\right)}^2-2{\left(\overline{c_S}+i\widetilde{c_S}\right)}^2& 0\\ {\left(\overline{c_P}+i\widetilde{c_P}\right)}^2-2{\left(\overline{c_S}+i\widetilde{c_S}\right)}^2& {\left(\overline{c_P}+i\widetilde{c_P}\right)}^2& 0\\ 0& 0& {\left(\overline{c_S}+i\widetilde{c_S}\right)}^2\end{array}\right]$$

(7)

where $i=\sqrt{-1}$, variables with $-$ and $~$ symbols, respectively, represent the real and imaginary parts. Note that setting $c_P$ and $c_S$ is the core step in extracting the wavefield. $\overline{c_P}$ and $\overline{c_S}$ are consistent with the P- and S-wave velocity. $\widetilde{c_P}$ and $\widetilde{c_S}$ are obtained by a dimensionless multiplied $\overline{c_P}$ and $\overline{c_S}$ near the exploration target, while $\widetilde{c_P}$ and $\widetilde{c_S}$ are equal to 0 in the unconcerned area. The location where $\widetilde{c_P}$ and $\widetilde{c_S}$ are non-zero is referred to as the stained target.

By substituting Equations (5)–(7) into the VDWE, we obtain,

$$\{\begin{array}{c}\rho \frac{\partial \left(\overline{\mathit{V}}+i\tilde{\mathit{V}}\right)}{\partial t}=ℒ\left(\overline{\mathbf{\sigma }}+i\tilde{\mathbf{\sigma }}\right)\\ \frac{\partial \left(\overline{\mathbf{\sigma }}+i\tilde{\mathbf{\sigma }}\right)}{\partial t}=\rho \left(\overline{\mathit{A}}+i\tilde{\mathit{A}}\right){ℒ}^T\left(\overline{\mathit{V}}+i\tilde{\mathit{V}}\right)\end{array}$$

(8)

$$\{\begin{array}{c}\frac{\partial \left(\overline{{\tau }^P}+i\widetilde{{\tau }^P}\right)}{\partial t}=\rho {\left(\overline{c_P}+i\widetilde{c_p}\right)}^2\left(\frac{\partial \left(\overline{v_x}+i\widetilde{v_x}\right)}{\partial x}+\frac{\partial \left(\overline{v_z}+i\widetilde{v_z}\right)}{\partial z}\right)\\ \rho \frac{\partial \left(\overline{v_x^P}+i\widetilde{v_x^P}\right)}{\partial t}=\frac{\partial \left(\overline{{\tau }^P}+i\widetilde{{\tau }^P}\right)}{\partial x}\\ \rho \frac{\partial \left(\overline{v_z^P}+i\widetilde{v_z^P}\right)}{\partial t}=\frac{\partial \left(\overline{{\tau }^P}+i\widetilde{{\tau }^P}\right)}{\partial z}\end{array}$$

(9)

Taking the real parts of Equations (8) and (9), we acquire an equation expressed as,

$$\{\begin{array}{c}\rho \frac{\partial \overline{\mathit{V}}}{\partial t}=ℒ\overline{\mathbf{\sigma }}\\ \frac{\partial \overline{\mathbf{\sigma }}}{\partial t}=\rho \left(\overline{\mathit{A}}{ℒ}^T\overline{\mathit{V}}-\tilde{\mathit{A}}{ℒ}^T\tilde{\mathit{V}}\right)\\ \frac{\partial \overline{{\tau }^P}}{\partial t}=\rho \alpha \left(\frac{\partial \overline{v_x}}{\partial x}+\frac{\partial \overline{v_z}}{\partial z}\right)-2\rho \beta \left(\frac{\partial \widetilde{v_x}}{\partial x}+\frac{\partial \widetilde{v_z}}{\partial z}\right)\\ \rho \frac{\partial \overline{v_x^P}}{\partial t}=\frac{\partial \overline{{\tau }^P}}{\partial x}\\ \rho \frac{\partial \overline{v_z^P}}{\partial t}=\frac{\partial \overline{{\tau }^P}}{\partial z}\end{array}$$

(10)

where $\alpha ={\overline{c_P}}^2-{\widetilde{c_P}}^2$ and $\beta =\overline{c_P}\widetilde{c_P}$. The real parts of S-waves are defined as $\overline{v_x^S}=\overline{v_x}-\overline{v_x^P}$ and $\overline{v_z^S}=\overline{v_z}-\overline{v_z^P}$. The imaginary parts of Equations (8) and (9) can be written as,

$$\{\begin{array}{c}\rho \frac{\partial \tilde{\mathit{V}}}{\partial t}=ℒ\tilde{\mathbf{\sigma }}\\ \frac{\partial \tilde{\mathbf{\sigma }}}{\partial t}=\rho \left(\overline{\mathit{A}}{ℒ}^T\tilde{\mathit{V}}-\tilde{\mathit{A}}{ℒ}^T\overline{\mathit{V}}\right)\\ \frac{\partial \widetilde{{\tau }^P}}{\partial t}=\rho \alpha \left(\frac{\partial \widetilde{v_x}}{\partial x}+\frac{\partial \widetilde{v_z}}{\partial z}\right)+2\rho \beta \left(\frac{\partial \overline{v_x}}{\partial x}+\frac{\partial \overline{v_z}}{\partial z}\right)\\ \rho \frac{\partial \widetilde{v_x^P}}{\partial t}=\frac{\partial \widetilde{{\tau }^P}}{\partial x}\\ \rho \frac{\partial \widetilde{v_z^P}}{\partial t}=\frac{\partial \widetilde{{\tau }^P}}{\partial z}\end{array}$$

(11)

The imaginary parts of S-waves are defined as $\widetilde{v_x^S}=\widetilde{v_x}-\widetilde{v_x^P}$ and $\widetilde{v_z^S}=\widetilde{v_z}-\widetilde{v_z^P}$. Equations (10) and (11) jointly constitute the stained VDWE (SVDWE). According to Equation (11), we find that $\tilde{\mathit{V}}$ is only excited when $\overline{\mathit{V}}$ arrives at the stained target. The imaginary part $\tilde{\mathit{V}}$ is named as stained wavefield. The term $\tilde{\mathit{A}}{ℒ}^T\overline{\mathit{V}}$ can be considered as the source signature of $\tilde{\mathit{V}}$. $\tilde{\mathit{A}}$ and $\beta$ are dimensionless due to the principle of $\widetilde{c_P}$ and $\widetilde{c_S}$ setting. It is evident that Equation (10) is approximately equivalent to the VDWE.

2.3. Inner Product Imaging Condition with Filters

Imaging condition can convert seismic data into structure images through the focus of the source and receiver wavefield. For the ERTM, the inner product imaging condition has an advantage in terms of definite physical meanings. However, applying such imaging condition inevitably results in migration artifacts because of backscatter. To suppress the unwanted seismic events in the migration profile, we add filters to the imaging conditions. The final images can be expressed as:

$$\widetilde{I^{PP}}={{\displaystyle \sum }}_{t=1}^T\left(F_{So}^P\widetilde{{\mathit{V}}_{So}^P}\right)·\left(F_{Re}^P\widetilde{{\mathit{V}}_{Re}^P}\right)$$

(12)

$$\widetilde{I^{PS}}={{\displaystyle \sum }}_{t=1}^T\left(F_{So}^P\widetilde{{\mathit{V}}_{So}^P}\right)·\left(F_{Re}^S\widetilde{{\mathit{V}}_{Re}^S}\right)$$

(13)

where $\widetilde{I^{PP}}$ and $\widetilde{I^{PS}}$ are PP and PS images, the symbol · denotes the inner product operation. $\tilde{\mathit{V}}$ with the subscripts $So$ and $Re$ denotes the stained wavefields of seismic sources and receivers. The superscripts $P$ and $S$ represent P- and S-wave, respectively. $T$ is the maximum time of seismic records. $F_{So}^P$ is a filter with the mathematical form shown below:

$$F_{So}^P=\{\begin{array}{c}1, \widetilde{{\mathit{D}}_{So}^P}·{\mathit{n}}_{So}^P>0\\ 0, \widetilde{{\mathit{D}}_{So}^P}·{\mathit{n}}_{So}^P<0\end{array}$$

(14)

where $\widetilde{{\mathit{D}}_{So}^P}=\left(\widetilde{d_x^P,}\widetilde{d_z^P}\right)$ is the Poynting vector, ${\mathit{n}}_{So}^P$ is the direction vector. If the result of the inner product between $\widetilde{{\mathit{D}}_{So}^P}$ and ${\mathit{n}}_{So}^P$ is positive (negative), then the wavefield is retained (muted). Other filters have a similar expression. The Poynting vector can be calculated by:

$$\{\begin{array}{c}\widetilde{d_a^P}=-\widetilde{{\tau }^P}\widetilde{v_a^P}\\ \widetilde{d_a^S}=-\left[\widetilde{{\tau }_{ab}}-\widetilde{{\tau }^P}\delta \left(a-b\right)\right]\widetilde{v_b^S}\end{array}$$

(15)

where $a$ and $b$ are the $x$- or $z$-component of the Poynting vectors, and $\delta$ is the Dirac function. From formula (15), the Poynting vector is 0 if either particle-velocity and stress is 0. There is a position between the crest and the trough where the wavefield is 0, but the Poynting vector at that position should not be 0. To improve the stability of equation (15), the local least squares strategy [25] is applied.

3. Numerical Examples

In the first example, we employ a layered velocity model to probe the performance of the above equations. The P-wave migration velocity model is sketched in Figure 1a. The S-wave migration velocity is obtained using an empirical formula $c_S=c_P/\sqrt{3}$. The density is set as 2000 kg/m³. The model is dispersed by 650 × 400 grid numbers, with a grid size of 10 × 10 m². A Ricker wavelet with the dominant frequency of 25 Hz is set as the source signature at the location of (1 km, 0). The wavefields are simulated by the finite difference method with the second-order in time and tenth-order in space. In addition, we add a perfectly matched layer around the simulated area to absorb spurious reflections from the truncated boundary.

Based on the VCWE, the $x$-component of the particle-velocity is shown in Figure 2a, in which the P- and S-waves are coupled together. When the coupled wavefield is extrapolated by solving Equation (4), it can be completely decoupled into pure P- and S-waves (Figure 2b,c). Notably, the seismic energy of the deep structure information is relatively weak in the decoupled wavefield.

To reduce the shielding effect of the overburden rock on seismic waves, we simulate the wavefield based on SVDWE. Here, the second interface is set as the stained target. $\widetilde{c_P}$ is obtained by $\overline{c_P}$ multiplied 10⁻⁶ in the target area while is zero in other locations. The setting of $\widetilde{c_S}$ is similar to that of $\widetilde{c_P}$. Considering the similarity between $\widetilde{c_P}$ and $\widetilde{c_S}$ models, only one of them is shown in Figure 1b. By solving SVDWE, the stained wavefield $\tilde{\mathit{V}}$ is generated after the $\overline{\mathit{V}}$ hits the stained target. The stained wavefield $\tilde{\mathit{V}}$ only contains the reflection and transmission about the second interface as exhibited in Figure 3. Consequently, we infer that imaging with the stained wave can reduce interference from non-target structures.

Furthermore, we utilize a contrastive approach to illustrate the adaptability of the ERTM for different equations. We evenly deploy 10 shots and 650 receivers on the surface of the layered velocity model. The total length of seismic records is 3 s, with a step time of 1ms. Figure 4 shows the multi-shot stacked images based on VDWE. Moreover, we implement the ERTM using the SVDWE, the PP and PS images are sketched in Figure 5. The settings of the stained target are the same as those of previous numerical experiments. Thanks to the stained wavefields, the primary information in the images is the second interface. However, some unwanted migration artifacts may interfere with subsequent interpretations.

Next, we demonstrate that the proposed imaging condition plays an important role in eliminating migration artifacts. The filters in the imaging condition need to be prepared in advance, and their properties are determined by the direction vector and Poynting vector. For the snapshot shown in Figure 3a, the Poynting vector (Figure 6) can be obtained using Equation (15). Here, we exhibit four filtered wavefields (Figure 7) that, respectively, correspond to different direction vectors: ${\mathit{n}}_{So}^P=\left(1, 0\right)$, ${\mathit{n}}_{So}^P=\left(-1, 0\right)$, ${\mathit{n}}_{So}^P=\left(0,-1\right)$, and ${\mathit{n}}_{So}^P=\left(0,1\right)$.

For flat strata, images with fewer artifacts can be obtained in the inner product of the up-going stained source wavefield and the down-going stained receiver wavefield (Figure 8). Figure 9 quantitatively compares PP and PS images, and the two sets of single-channel data are extracted from Figure 8 (from 2 km to 3 km in depth), and then mapped to the wavenumber domain through the Fourier transform. It is prominent that PS images are characterized by abundant high wavenumber information, robustly confirming that PS images have a higher spatial resolution than PP images. The basic reason for such a phenomenon is that S-wave travels slower than P-wave in the subsurface. Namely, the converted S-waves possess a shorter wavelength.

In this example, a salt model is applied to evaluate the efficiency of the proposed approach in complex geologic regions. The true P-wave velocity model is exhibited in Figure 10a, and S-wave velocity is also calculated using the previous empirical formula. The model size is 400 × 400 grids, with a grid size of 10 × 10 m² in both horizontal and vertical directions. A total of 40 sources are evenly excited on the surface of the model. Each source has 400 receivers with a spacing of 10 m. The maximum recording time is 6 s. Figure 11 shows the PP and PS images obtained by the conventional ERTM, and it is obvious that the shallow structure and the salt boundary can be effectively identified. However, the high-velocity overburden rock leads to insufficient illumination for the deep subsurface, which ultimately blurs the subsalt image. We set the stained target in Figure 10b. Figure 12 displays the PP and PS images with the proposed imaging condition. We can see that the two complementary images are favorable for the identification of subsalt faults.

In the last experiment, a field data is utilized to confirm the effectiveness of our method. Figure 13 displays the migration velocity model that has 940 × 300 grid points with a 20 × 20 m² element. There are 240 sources unevenly distributed at the free surface. Figure 14 show some shot gathers whose maximum recording time is 5 s. Since we only have the vertical component of the seismic record where P-wave information dominates, the valuable PS images are not obtained. The PP images based on VDWE are sketched in Figure 15, from which we can see the main structures. To enhance the quality of imaging in weakly illuminated areas, we set the box in Figure 13 as the stained targets. By using the proposed method, the migration events (Figure 16) relevant to the stained targets are highlighted, and the reflections of nontarget structures are muted.

4. Discussion

Our approach shows the satisfying performance for imaging weakly illuminated structures. When an overall migrated profile is acquired, it is proper to achieve locally enhanced imaging by artificially setting exploration targets in advance. Since the SVDWE contains two systems of equations, the algorithm of such an approach is twice as computationally intensive as the conventional one. The staining strategy is intrinsically related to the imaging results and is worthy of further investigation. The filters in the inner product imaging condition screen different wavefields according to the direction of wave propagation. Theoretically, it is possible to achieve wavefield decomposition in an arbitrary direction. However, in this paper, we decompose the wavefield into up- and down-going categories under the assumption that the target structure is relatively flat. To enhance the imaging adaptability for complex geological situations, the wavefield screening method according to stratigraphic dip angle needs to be persistently explored.

5. Conclusions

We have presented a new ERTM technique to produce the weakly illuminated structure. To extract the elastic wavefield related to the target, a stained vector-decoupled wave equation is adopted. The filter in the inner product imaging condition can suppress the migration artifact, improving the resolution of ERTM results. Numerical results show that our method can generate better amplitude for the weakly illuminated structure, providing possible technical support for reducing the risk of hydrocarbon development.