**1. Introduction**

Elastic reverse time migration (ERTM) possesses the ability to reposition the multicomponent 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 wideazimuth measurement, illumination compensation and acquisition aperture correction approaches have been developed to enhance the imaging quality of shadow regions [6–8]. In addition, the staining algorithm can highlight the wavefields associated with weakly illuminated structures among the full wavefields [9–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

**Citation:** Song, L.; Shi, Y.; Liu, W.; Zhao, Q. Elastic Reverse Time Migration for Weakly Illuminated Structure. *Appl. Sci.* **2022**, *12*, 5264. https://doi.org/10.3390/ app12105264

Academic Editors: Guofeng Liu, Xiaohong Meng and Zhifu Zhang

Received: 26 March 2022 Accepted: 20 May 2022 Published: 23 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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{cases} \begin{array}{c} \rho \frac{\partial V}{\partial t} = \mathcal{L}\sigma\\ \frac{\partial \sigma}{\partial t} = \rho \mathcal{A} \mathcal{L}^T V \end{array} \end{cases} \tag{1}$$

where *ρ* denotes the density of media, *V* = [*vx*, *vz*] *<sup>T</sup>* is the particle-velocity vector, *t* signifies the time variable, σ = [*τxx*, *τzz*, *τxz*] *<sup>T</sup>* denotes the stress vector and L is a spatial partial derivative operator expressed as:

$$\mathcal{L} = \begin{bmatrix} \partial\_x & 0 & \partial\_z \\ 0 & \partial\_z & \partial\_x \end{bmatrix} \tag{2}$$

*x* and *z*, respectively represent the horizontal and vertical directions, and *A* is a parameter matrix as shown below:

$$A = \begin{bmatrix} 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{bmatrix} \tag{3}$$

where *cP* and *cS* 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 *τ<sup>P</sup>* to decouple the wavefield, and the formulation is written as,

$$\begin{cases} \begin{array}{l} \frac{\partial \mathbf{r}^P}{\partial t} = \rho c\_P^2 \left( \frac{\partial \upsilon\_x}{\partial x} + \frac{\partial \upsilon\_z}{\partial z} \right) \\ \rho \frac{\partial \upsilon\_x^P}{\partial t\_p} = \frac{\partial \tau^P}{\partial x} \\ \rho \frac{\partial \upsilon\_z^P}{\partial t} = \frac{\partial \tau^P}{\partial z} \end{array} \tag{4}$$

where *v<sup>P</sup> <sup>x</sup>* and *v<sup>P</sup> <sup>z</sup>* are the decoupled P-waves in *x*- and *z*-directions. The residual S-waves are defined as *v<sup>S</sup> <sup>x</sup>* = *vx* − *v<sup>P</sup> <sup>x</sup>* and *v<sup>S</sup> <sup>z</sup>* = *vz* − *v<sup>P</sup> <sup>z</sup>* . 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).
