1 Introduction

Geophysical and geological exploration of ocean-bottom environments has attracted considerable research efforts in petroleum and geological exploration. Ocean-Bottom Node (OBN) acquisition using acoustic sources near sea- surface and four-component receivers at the sea floor offer obvious advantages over conventional tower-streamer acquisition, such as improved imaging in the presence of obstructions such as drilling platforms, full-azimuth acquisition, and wider seismic bandwidth. Additionally, seafloor acquisition also allows the imaging of seismic energy converted from P-wave to S-wave at subsurface reflection points, which helps in overcoming limitations in high-resolution seismic images (Ker et al., 2010). However, density variations are often ignored in many prestack migration algorithms (Li et al., 2016). S-wave reflection amplitude is a function of both velocity and density (Jin et al., 2000), which contains information about the velocity and density distributions (Mora, 1987). Furthermore, converted-wave amplitude variation patterns reveal the changes in S-wave velocity and density, thus providing important lithologic information that is useful for hydrocarbon discrimination. Conventional methods of PS imaging have been operated in time-space domains based on second-order scalar acoustic wave equations that are not staggered grids of a finite difference. As shown in Fig.1, the P-wave can be imaged using a P-velocity model for data extrapolation, and the S-wave can be imaged using an S-velocity model for reverse-time extrapolation (Sun and McMechan, 1986). Because of different ray paths, travel time is also not the same. S-waves have been obtained by elastic wavefield separation after numerical simulation, which may cause distortion on amplitude and phase.

In this paper, we put forward new elastic equations transformed from traditional elastic wave equations for converted S-wave simulation combined with a staggered grid finite difference. This represents great progress in computational efficiency, accuracy and dispersion suppression (Dablain, 1986). Then, after phase correction, we illustrated the method of the PS image on a 2-D synthetic data set.

2 Theory 2.1 Converted S-wave Simulation

The traditional 2-D isotropy elastic wave equation can be written as (Ma and Zhu, 2003):

$ \frac{{{\partial ^2}S}}{{\partial {t^2}}} = {v_p}^2\nabla (\nabla \cdot S) - {v_S}^2\nabla \times (\nabla \times S). $

(1)

Here, S is the elastic wavefield consisting of P-waves and S-waves, t is time, v_p and v_S represents velocities of P-wave and S-wave, $\bigtriangledown$ is the gradient operator and $\bigtriangledown$× is a curl operator. The gradient operator can be described as the following:

$\nabla {\rm{ = }}\frac{\partial }{{\partial {\rm{x}}}}w + \frac{\partial }{{\partial z}}v.$

(2)

Here, w and v are the unit vectors along the x- and z-axes, respectively. It is worth pointing out that conventional elastic waves only provide a mixed wavefield of P- and S-waves.

It is important to perform wavefield separation if we want to image P and S reflection energy independently, where it is difficult to avoid distortion of amplitude and phase. Therefore, we introduced two variables – S_p and S_s. Eq. (1) can then be expressed as:

$\begin{align} & \frac{{{\partial }^{2}}{{S}_{p}}}{\partial {{t}^{2}}}={{v}_{p}}^{2}\nabla (\nabla \cdot S) \\ & \frac{{{\partial }^{2}}{{S}_{s}}}{\partial {{t}^{2}}}=-{{v}_{S}}^{2}\nabla \times (\nabla \times S). \\ & S={{S}_{p}}+{{S}_{s}} \\ \end{align}$

(3)

Thus, new Eq. (3) has been obtained from the conventional elastic wave equation. Tang and He (2016) gave proof that S_p was the P-wave displacement vector and S_s the S-wave displacement vector. The formula can thus be written as (Tang and He, 2016):

$\begin{gathered} \frac{{{\partial ^2}}}{{\partial {t^2}}}(\nabla \times {S_p}) = 0 \\ \frac{{{\partial ^2}}}{{\partial {t^2}}}(\nabla \cdot {S_s}) = 0 \\ \end{gathered} .$

(4)

Here, $\bigtriangledown$ is the gradient operator and $\bigtriangledown$× is the curl operator.

2.2 Reverse Time Migration

There are many migration approaches that can be used to produce an image of the subsurface by extrapolating both the source and the receiver wavefields into the earth (Mora, 1987; Yoon and Marfurt, 2006; Lu et al., 2009). RTM, as a full wave-equation algorithm, is a powerful tool to effectively handle multi-arrival seismic waves and accurately image overhanging reflectors and steep dips. For PS imaging, Sun and McMechan (2001) and Sun et al., 2011) put forward a preliminary method of scalar elastic wave pre-stack reverse-time depth migration. In this method, converted S-waves are treated as acoustic waves when performing reverse-time wavefield extrapolation. However, the algorithm needs another wavefield extrapolation and compared with SSG, is lacking in computational efficiency and accuracy.

So, in order to take density into account, in 2-D isotropy media, the first-order velocity-stress equation is (Zhang and Shi, 2017):

$\begin{align} & \frac{\partial P}{\partial t}=-\rho {{v}^{2}}(\frac{\partial {{v}_{x}}}{\partial x}+\frac{\partial {{v}_{z}}}{\partial z}) \\ & \frac{\partial {{v}_{x}}}{\partial t}=\frac{1}{\rho }\frac{\partial P}{\partial x} \\ & \frac{\partial {{v}_{z}}}{\partial t}=\frac{1}{\rho }\frac{\partial P}{\partial z} \\ \end{align}.$

(5)

Here, P is acoustic wavefield and v_x represents the x- component of velocity and v_z represents the z-component. ρ is density. Imaging condition is thus:

$I(x) = \frac{{\int {{P_s}(x, t) \cdot {S_r}} (x, t){\rm{d}}t}}{{\int {{P_s}(x, t) \cdot {P_s}} (x, t){\rm{d}}t}}.$

(6)

Here, I represent image value, P_s is the P-wave of the source-side, and S_r is the S-wave of the receiver side (Chen and He, 2014). Compared with the scalar acoustic wave equation, the first-order velocity-stress equation takes density into consideration and offers another velocity component when performing finite difference simulations. This helps in multi-component joint processing. Perfectly Matched Layer is used to reduce reflections from boundaries. Poynting vectors are also used to suppress low-frequency noise after RTM.

3 Examples

Fig.2 Compares improved elastic wave equations with conventional equation The model dimensions were 500× 500, with 3 m grid intervals, with the source at (250, 250), the improved equations offer completely separated wavefield of P-wave and converted S-wave, which help to avoid waveform distortion when applying wavefield separation before P- and S-image independently.

The effect of density variations on synthetic data was first illustrated with a two-layer model as shown in Fig.3. Amplitude difference between the constant density finite simulations (in black) and variable density (in blue) were considerable (offset = 300 m). Density was calculated by the Gardiner formula and the velocity of the P-wave in the first layer was 1000 m s^-1. In practice, we usually tended to focus on parameters that are the main influences on seismic imaging (velocity), while the others (density) were permitted to be neglected. However, seismic reflections amplitude was predominantly affected by contrasts in acoustic impedance, which was a product of velocity and density (Yang et al., 2016). It was our initial objective to focus on imaging of converted waves in the presence of density variations.

It should also be pointed out that in Fig.4a, phase had been changed between different sides of source; phase correction was required before pre-stack migration to avoid constructive interference between data from adjacent sources (Sun et al., 2001; Du et al., 2012). The formula used to correct phase was:

$T = \left({\begin{array}{*{20}{c}} {{T_x}} \\ {{T_z}} \end{array}} \right) = \left({\begin{array}{*{20}{c}} {T{v_x}} \\ {T{v_z}} \end{array}} \right), $

(7)

where T is the poynting vector of the mixed wavefield, and Tv_x and Tv_z represent poynting vectors of wavefields of the dilatation-component and the rotation-component. Based on Eq. (7), we could calculate T_s which represents the direction of the source-side wavefield propagation, and T_r which represents the direction of the receiver-side wavefield propagation.

The synthetic data in Fig.4a was calculated in a 2-D homogeneous model using Eq. (3). During elastic wavefield propagation, if T_s × T_r > 0, the converted S-wave phase was taken to be the same, while, if T_s × T_r < 0, the phase was taken as opposite. Fig.4b shows that after phase correction, constructive interference from adjacent sources was averted.

Here, we first demonstrate the advantages of elastic VD- RTM using a part of the Marmousi2 model. The model dimensions were 1688 × 368, with 3 m grid intervals. The observed data were generated by solving the proposed Eq. (3), including 50 OBN nodes at a depth of 108 m, and 150 shots at a depth of 3 m below the seafloor. Shots and receivers intervals were, respectively, 30 m and 90 m, with a 30 Hz peak frequency Ricker wavelet. Time sampling interval was 0.5 ms. Velocity and density models are shown in Fig.5. After numerical simulation, imaging the PP and PS reflections energy independently thus became possible. Following conventional methods, density was treated as a constant when performing migration.

However, if the P- wave velocity was the only variable parameter in migration when, in reality, density varied, the layers were not correctly imaged. For comparison purposes, when applying reverse-time migration to converted S-wave with constant density (CD-RTM), we used both a real-density model (Fig.5c) and hypothetical constant-density model (1000 kg m^-3). Figs. 6a and 6b compare CD-RTM and VD-RTM images using the data set generated with density variations. The black arrows highlight the noticeable improvements in preserving event continuities with VD-RTM compared with CD-RTM (Fig.6e). In this case, migration velocity was precise. We would speculate that the event discontinuities were caused by density errors.

On the other hand, RTM is extremely dependent on its velocity model; while, with both of these two methods, migration velocity was precise, improvements in signal-to- noise ratio and vertical resolution were not very obvious. At same time, with errors due to a lack of iterations and approximations of the wave equation (Feng and Gerard, 2016; Yang et al., 2016; Youn and Zhou, 2001), we may not be able to obtain amplitude-preserved RTM images. However, the relative amplitudes of the reflectors can still be estimated correctly, which is beneficial for reservoir characterization.

4 Conclusion

In this paper, we focused on RTM of S-waves in the presence of density variations and proposed improved elastic equations transformed from traditional elastic wave equations for numerical simulation. We then imaged converted S-waves by RTM-based first-order velocity-stress equations. Synthetic data based on the Marmousi2 model demonstrates the superior quality of migration with variable density.

Acknowledgements

This study was supported by the National Science and Technology Major Project (No. 2016ZX05027-002), the National Natural Science Foundation of China (No. 41230 318), and Qingdao National Laboratory for Marine Science and Technology Innovation Project of Ao-Shan (No. 2015ASKJ03). We are also grateful to GeoEast software.