**1. Introduction**

At present, reverse-time migration imaging simulation is an important means to explore the morphology of underground media. As the reverse-time migration imaging method is based on the two-way wave equation, which can accurately describe the dynamic and kinematic characteristics of a seismic wavefield propagating underground, reversetime migration has no inclination angle limit and can adapt to the imaging of complex structural areas, especially for structures with clear lateral velocity changes [1–4]. The finite difference scheme is widely used in the numerical simulation of the elastic-wave equation because of its simplicity and flexibility, high calculation efficiency, and small memory requirement [5–13]. On the one hand, with the development of multicomponent seismic exploration in recent years, particularly shear-wave seismic exploration, in order to minimize computational costs, there is a need to increase the time and spatial steps used in finite difference modeling while maintaining sufficient accuracy during numerical simulation [14]. On the other hand, if the traditional finite difference scheme is used for numerical simulation, small time and spatial steps are required to achieve sufficient accuracy. Compared with the traditional difference scheme, the compact difference scheme has the same accuracy and high stability. The development of compact difference schemes can be traced back to 1989 when Dennis and Hudson (1989) first proposed spatial fourth-order compact schemes for Navier–Stokes equations [15]. Lele (1992) studied the Pade scheme and proposed a symmetric compact difference scheme for solving the first and second derivatives,

**Citation:** Zhou, C.; Wu, W.; Sun, P.; Yin, W.; Li, X. The Combined Compact Difference Scheme Applied to Shear-Wave Reverse-Time Migration. *Appl. Sci.* **2022**, *12*, 7047. https://doi.org/10.3390/ app12147047

Academic Editor: Filippos Vallianatos

Received: 18 May 2022 Accepted: 6 July 2022 Published: 12 July 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/).

and the accuracy of the scheme can reach that of the spectral method [16]. Chu et al. (2000) used the combined compact difference scheme (CCD) for the convection–diffusion equation [17].

Since the mid-1990s, reverse-time migration has been applied to multicomponent wave seismic data excited from P-wave sources [18–20], having overcome calculation problems and interference artifacts in P- and S-wave simulations. During reverse-time migration, the P- and S-wave vectors of the source wavefield and the P- and S-wave vectors of the receiver wavefield are obtained and then imaged. In general, P-and S-wave wavefields are obtained either via the Helmholtz decomposition or by using the pled elastic-wave equation. These approaches are all designed for P-wave and elastic-wave sources, focusing mainly on how to retain the phase, amplitude, wavefield and vector characteristics of the wavefield efficiently during wavefield separation [21–28]. In recent years, because of a breakthrough in S-wave vibrator technology, pure S-wave seismic data can be obtained via artificial excitation [29,30]. As the wavelet length of an S-wave is shorter than that of a P-wave for the same frequency bandwidth, the resolution of its wavefield is higher, and its advantage for imaging beneath gas clouds area is unmatched by a P-wave [31]. Using the S-wave wavefield generated by an S-wave source for reverse-time migration imaging (RTM) is a key step for processing shear-wave data. Due to the characteristics of shear-wave wavelength, higher accuracy is required in numerical simulation. To ensure the accuracy of shear-wave RTM, a smaller time and spatial step than P-wave is needed, which reduces the calculation efficiency.

To improve the accuracy and efficiency of shear-wave RTM, based on the characteristics of the shear-wave velocity model, we used the combined compact difference (CCD) and combined supercompact difference scheme (CSCD) for shear-wave (SH) RTM, with larger spatial grid conditions. Wang Shuqiang (2002) applied the compact difference scheme to the numerical simulation of seismic wavefields [32]. Most recent studies of the finite difference scheme focus on the compressional wavefield from an explosion source, but few have studied the finite difference scheme based on the shear wavefield of the shear-wave source (SH-wave).

This paper aims at the problem of low shear-wave velocity by introducing the supercompact difference scheme to suppress the numerical dispersion caused by the large spatial step. We introduce the basic concept of shear-wave reverse-time migration, followed by the implementation methods of the combined compact difference scheme (CCD) and combined supercompact difference scheme (CSCD); the numerical simulation accuracy of CCD and CSCD is also discussed, and the accuracy of CCD and CSCD is compared with the traditional finite difference scheme. Finally, the method is applied to synthetic data to verify the accuracy and efficiency of the algorithm. Furthermore, extending the present research to take into the viscoelastic behavior of media will have a great potential for imaging oil and gas reservoirs [33].

### **2. Principle of RTM and Combined Compact Difference Scheme**

### *2.1. Principle of RTM*

The technique of reverse-time migration imaging (RTM) is composed of the following three steps [34]:


$$I(\mathbf{x}, z) = \int\_0^T \mathcal{S}(\mathbf{x}, z, t) \mathcal{R}(\mathbf{x}, z, t) dt \tag{1}$$

where *S*(*x*, *z*, *t*) is the source wavefield obtained via forward modeling, *R*(*x*, *z*, *t*) is the receiver wavefield obtained at the same time via reverse continuation simulation under the same velocity model, and *t* is the total propagation time.

It can be seen from Equation (1) that the final result of reverse-time migration is affected by the accuracy of the source wavefield *S*(*x*, *z*, *t*) and the receiver wavefield *R*(*x*, *z*, *t*). Notably, the method to improve the accuracy of the source wavefield *S*(*x*, *z*, *t*) can also be applied to the receiver wavefield *R*(*x*, *z*, *t*). Thus, a high-precision finite difference scheme is applied to generate the source and receiver wavefields *S*(*x*, *z*, *t*) and *R*(*x*, *z*, *t*), which will give an ideal result of reverse-time migration. In this paper, we adopted cross-correlation imaging conditions; in addition, the method in this paper can also be applied to some imaging conditions developed by researchers in recent years [35–37].
