Implementation of Elastic Prestack Reverse-Time Migration Using an Efficient Finite-Difference Scheme

Elastic reverse-time migration (RTM) can reflect the underground elastic information more comprehensively than single-component Pwave migration. One of the most important requirements of elastic RTM is to solve wave equations. The imaging accuracy and efficiency of RTM depends heavily on the algorithms used for solving wave equations. In this paper, we propose an efficient staggered-grid finite-difference (SFD) scheme based on a sampling approximation method with adaptive variable difference operator lengths to implement elastic prestack RTM. Numerical dispersion analysis and wavefield extrapolation results show that the sampling approximation SFD scheme has greater accuracy than the conventional Taylor-series expansion SFD scheme. We also test the elastic RTM algorithm on theoretical models and a field data set, respectively. Experiments presented demonstrate that elastic RTM using the proposed SFD scheme can generate better images than that using the Taylor-series expansion SFD scheme, particularly for PS images. FurH. thermore, the application of adaptive variable difference operator lengths can effectively improve the computational efficiency of elastic RTM.

ment RTM with angle-domain imaging formulated for multicomponent elastic data. Du et al. (2012Du et al. ( , 2014 followed the procedure presented by Yan and Sava (2008) to perform elastic RTM in the common-shot domain and discussed the polarity reversal correction by using the energy flux density vector. Chung et al. (2012) proposed a frequency domain elastic RTM algorithm in the common shot domain based on the wavefield separation technology. Yan and Xie (2012) presented an angle-domain imaging condition using the local slant stack method for multicomponent elastic RTM.
Elastic RTM requires forward extrapolation of the source wavefield and backward extrapolation of the recorded receiver wavefield in time (Yan and Sava 2008), and the essence of wavefield extrapolation is the numerical solution of wave equations. The imaging accuracy and efficiency of elastic RTM depends heavily on the algorithms used for solving wave equations. Hence, how to solve wave equations is very important to the implementation of elastic RTM. In seismic exploration, the most popular implementation for solving equations in wavefield extrapolation is finite-difference (FD) methods, especially staggered-grid FD (SFD) methods , Liu and Sen 2011b, Yang et al. 2014. However, numerical approximations of the FD methods are required for spatial and time derivatives, and the conventional FD methods often suffer from serious numerical errors (numerical dispersion) (Li et al. 2013). Traditionally, accurate numerical approximations are achieved by using either relatively very fine computation grids or very long FD operators (Dablain 1986). Otherwise, numerical dispersion will be present in the image data, which contaminates the signals. However, either approach will increase the computational cost dramatically ). Another alternative is to optimize the FD coefficients to suppress the numerical dispersion in seismic modeling and imaging (e.g., Li et al. 2015.
In this paper, we propose an efficient SFD scheme based on a sampling approximation method with adaptive variable difference operator lengths to implement elastic prestack RTM. First, we simply review the conventional SFD scheme based on Taylor-series expansion method, and introduce a high-accuracy SFD scheme based on a sampling approximation method and adaptive variable difference scheme. Then, we analyze the accuracy of the numerical solution of elastic wave equations using the SFD schemes, including the Taylor-series expansion SFD and sampling approximation SFD schemes. Next, we introduce the basic workflow of the elastic RTM. Finally, we adopt the SFD scheme based on the sampling approximation method with adaptive variable difference operator lengths to implement forward and backward wavefield extrapolations for elastic prestack RTM, testing the algorithm with two theoretical models and a field data set, respectively.

METHOD 2.1 SFD schemes for the numerical solutions of elastic wave equations
In 2D heterogeneous media, the elastic velocity-stress equations are (Virieux 1986) where t is the time, x and z are the space coordinates, (v x ,v z ) is the particle velocity vector, (τ xx ,τ zz ,τ xz ) is a vector containing three components of stress, ρ(x,z), is the density, and λ(x,z), and μ(x,z) are Lame's constants.
To solve the 2D elastic velocity-stress equations, the high-order SFD schemes are generally used to calculate the spatial first-order derivatives. The (2M)th-order SFD scheme for the first-order derivative of the function p(x) can be expressed as (Kindelan et al. 1990) where h is the grid size and a m are difference coefficients on the first-order derivative.
Using plane wave theory, we obtain: where p 0 is a constant, 1 i = − , and k is the wavenumber. Substituting Eq. 3 into Eq. 2, the dispersion relation can be obtained as follows (e.g., Yang et al. 2014: where: β = kh/2 and 0 ≤ β ≤ π/2. When the Taylor-series expansion method (e.g., Dong et al. 2000, Pei 2004, Liu and Sen 2009) is employed to compute the SFD coefficients, Eq. 4 can be expanded as Comparing the coefficients of β in Eq. 5, we obtain the SFD coefficients based on Taylor-series expansion method as follows (Liu and Sen 2009): When the conventional classic coefficients determined by the Taylorseries expansion method are used to solve wave equations for a larger frequency content, serious numerical dispersion will occur, which will affect the seismic modeling accuracy and imaging quality .
We introduce a high-accuracy SFD scheme with weak-dispersion to solve elastic wave equations. First, we take M sampling points for β, which can be expressed as β 1 ,β 2 ,…β M , and are distributed evenly over the range from 0 to a given u, where u is a constant, and u ∈ (0,π/2]. Then we approximate the dispersion relation (Eq. 4) at the M sampling points. When β 1 = 0, the equality of Eq. 4 can be unconditional. The basic idea is similar to the optimization method of Yang et al. (2015), but their optimization method is for rotated SFD scheme. The near zero-wavenumber constraint condition is usually adopted to improve the accuracy of FD (e.g., Zhou and Zhang 2011, Liu 2014. We introduce the constraint condition by using the following formula: Then, we obtain: Subsequently, we can use the rest of the (M − 1) sampling points (β 1 , β 2 , β M ) to construct a system of linear equations from Eq. 9, and then obtain the SFD coefficients for first-order derivatives by solving these linear equations and Eq. 8. The difference coefficients can be simply obtained by this method. We call the method sampling approximation SFD scheme.
According to Eq. 4, the parameter δ(β) is defined to describe the numerical dispersion in SFD modeling: If δ(β) is close to zero, the numerical dispersion is weak. When δ(β) is much greater or less than zero, large numerical dispersion will exist. To compare the sampling approximation SFD with Taylor-series expansion SFD schemes, we analyze the numerical dispersion relations by showing the dispersion curves. Figure 1 shows variations of δ(β) with β for different M by the Taylor-series expansion SFD and the sampling approximation SFD for numerical dispersion, respectively. From Fig.1, we can see that the sampling approximation SFD scheme has greater accuracy than the Taylor-series expansion SFD scheme over a wider range of wavenumbers. Namely, the sampling approximation SFD scheme suffers much less numerical dispersion than the Taylor-series expansion SFD scheme for the same spatial difference operator length. In addition, the difference operator length of sampling approximation SFD scheme is shorter than that of Taylor-series expansion SFD scheme with the same numerical accuracy.
Generally, the fixed difference operator length is used to compute the spatial derivatives, which leads to more computing time. Based on the idea a) b) of the adaptive variable-length FD operators presented by Liu and Sen (2011a), we propose the sampling approximation SFD scheme with adaptive variable difference operator lengths. It uses long operators in regions of low velocity and short operators in regions of high velocity to improve the computational efficiency for the solution of elastic wave equations in elastic prestack RTM.
According to Eq. 10, the numerical error ε is defined as follows: where f is the frequency and v is velocity. Here, we adopt the S-wave velocity to calculate the numerical error because the S-wave dispersion is dominant in elastic wave propagation.
From Eq. 11, we find that ε is a function of v, M and f. Therefore, for the given maximum frequency f max and maximum numerical error η, the following inequality is satisfied (e.g., Liu and Sen 2011a) where f ≤ f max . The adaptive variable SFD operator lengths can be determined by the inequality (Eq. 12).

The basic workflow of elastic prestack RTM
The implementation of elastic prestack RTM consists of three main parts: forward extrapolation of the source wavefield in time, backward extrapolation of the recorded receiver wavefield in time, and the application of an imaging condition to the pure waves (Yan and Sava 2008). The workflow includes the same three typical parts found in acoustic RTM. However, more operations are performed. After reconstructing the source wavefield and the receiver wavefield, we should separate the elastic wavefields into P-and Swaves, and then perform the phase correction, the amplitude balancing and the polarity reversal correction prior to applying the imaging condition. Figure 2 is the workflow chart for elastic prestack RTM. The extrapolated wavefield should be separated into P-and S-wave components after extrapolation and crosscorrelation of the vector and scalar potentials are used for imaging (Dellinger and Etgen 1990). In a 2D isotropic medium, the P-and S-wave components can be obtained by applying Helmholtz decomposition to the elastic extrapolated wavefield (e.g., Yan and Sava 2008 The pure S-waves are represented with the curl of the wavefield, We can use the high-order SFD method to compute spatial derivatives for Eqs. 13 and 14.
The phase and amplitude of P-waves and S-waves separated by Helmholtz decomposition have changed, so we need to correct the phase and amplitude for the wavefield. Here, we adopt the method proposed by Sun et al. (2001Sun et al. ( , 2011 to solve this issue: Implementing the phase correction by Hilbert transform for separated P-and S-waves; and performing the amplitude balancing by taking the ratio of the P-wave velocity to the S-wave velocity times the S-wave amplitude. Additionally, an S-image polarity reversal occurs because the S-wave changes its polarity when crossing normal incidence (Balch and Erdemir 1994). So the polarities must be corrected before applying the imaging conditions. Sun and McMechan (2001) specifically discussed the problem of polarity reversal for converted S-wave, and proposed that the polarity reversal could be corrected simply by multiplying the S-wave amplitudes by -1 on one side of the source position while keeping the other side unchanged after P-S wave separation. In this paper, we use this simple method (Sun and McMechan 2001) to perform polarity reversal correction in elastic RTM. Yan and Sava (2008) proposed crosscorrelation imaging conditions to obtain elastic images with clear physical meanings. To reduce migration artifacts and improve image quality, Du et al. (2012Du et al. ( , 2014 recommended source-normalized crosscorrelation imaging conditions. In this paper, the source-normalized crosscorrelation imaging conditions are adopted to perform the PP and PS imaging as follows (Du et al. 2012): , , , , where S(x,z,t) is the forward extrapolating source wavefield; R(x,z,t) is the backward extrapolating receiver wavefield; the subscripts P and S denote P and S components, respectively; and I pp (x,z) and I ps (x,z) are PP-and PSmode migrated images, respectively. The cross-correlation imaging conditions usually produce low frequency migration noises. In this paper, we adopt a Laplace filtering presented by Zhang and Sun (2009) to suppress the low frequency noises.

EXAMPLES 3.1 Elastic prestack RTM for a groove model
We use a groove model as shown in Fig. 3 to demonstrate the effects of elastic modeling and prestack RTM. The model is discretized into a 501 by 301 grid points, with a grid interval of 15 m. The P-and S-wave velocities and the density for the upper layer are 2500 m/s, 1500 m/s and 2000 kg 3 m -3 , respectively, and for the lower layer are 2800 m/s, 1800 m/s and 2100 kg 3 m -3 , respectively. We use the sampling approximation SFD scheme with 2ndorder accuracy in time and 16th-order accuracy in space to generate a synthetic multicomponent data on the groove model. We choose a Ricker wavelet with a peak frequency of 25 Hz as the compressional source function.
There are 42 sources in total, and the sources are placed at a horizontal interval of 150 m, the first source being located at (225 m, 30 m). The receivers are distributed on the surface, with an interval of 15 m. The sample interval is 0.001 s with a total propagation time of 5.0 s. We perform forward and backward wavefield extrapolations for the elastic prestack RTM using the SFD methods with 2nd-order accuracy in time and 8th-order accuracy in space. Here, the perfectly matched layer (PML) absorbing boundary condition (Bérenger 1994) is used to reduce unwanted reflections from artificial boundaries. Figures 4 and 5 show forward elastic wavefield snapshots at t = 1.35 s computed with the Taylor-series expansion SFD and the sampling approximation SFD schemes, respectively, when the source is located at (3525 m, 30 m). Because we use an explosive source, the snapshots contain a mix of P-and S-wave modes, as can be seen form the x and z components. Furthermore, as highlighted by the white arrows, the numerical dispersion is very serious in the snapshots computed by the Taylor-series expansion SFD scheme. However, there is very small numerical dispersion in the snapshots computed by the sampling approximation SFD scheme, which further demonstrates that the sampling approximation SFD scheme has better accuracy in solving elastic wave equations. We implement elastic prestack RTM using the Taylor-series expansion SFD scheme and the sampling approximation SFD scheme, respectively. Figures 6 and 7 show the final image results including PP images and PS images. Comparing Figs. 6 and 7, the differences among images generated by the Taylor-series expansion SFD scheme and the sampling approximation SFD scheme are obvious. And it can be observed that the events in Fig. 7 are clearer and more focused (highlighted by the white arrows). However, some artifacts along the reflector interfaces in the images migrated by the Taylorseries expansion SFD scheme are very serious (highlighted by the white arrows in Fig. 6), particularly for the PS image in Fig. 6b (highlighted by the white arrows), which results from numerical dispersion effects in the wave extrapolation. The serious dispersion in the extrapolating wavefield contaminates the images. Therefore, the numerical tests indicate the groove model is better imaged by the sampling approximation SFD scheme than the Taylor-series expansion SFD scheme.

Elastic prestack RTM for a salt model
To further demonstrate the imaging quality and examine the efficiency of difference scheme with adaptive variable difference operator lengths, we test the elastic prestack RTM algorithms on the modified 2D Society of Exploration Geophysicists/European Association of Geoscientists and Engineers (SEG/EAGE) salt model. The salt model has 600 (in x) by 200 (in depth) grid points, with grid intervals 20 m (in x) and 20 m (in depth). The P-wave velocity information, as shown in Fig. 8, varies 2000 m/s from 4981 m/s. The S-wave velocities (v s ) are created from P-wave velocities (v p ) following v p /v s = 1.7. The density is constant. We simulate a synthetic multicomponent data set which has 75 shots with a 80 m shot spacing from the leftmost location (3000 m, 40 m) using the sampling approximation SFD scheme with 2nd-order accuracy in time and 16th-order accuracy in space. The compressional source function is represented by a Ricker wavelet with a peak frequency of 30 Hz. Each shot is recorded by 300 receivers spaced 20 m apart. The time sampling is interval 0.001 s and the recording time is 5 s. We implement the elastic prestack RTM using the Taylor-series expansion SFD scheme with the fixed difference operator lengths, the sampling approximation SFD scheme with the fixed difference operator lengths, and the sampling approximation SFD scheme with the adaptive variable difference operator lengths, respectively. Figure 9 shows the fixed and the variable sampling approximation SFD operator lengths M used in elastic RTM. Figure 10 shows the final PP and PS images migrated by the Taylor-series expansion SFD scheme with fixed difference operator lengths. Figure 11 shows the final PP and PS images migrated by the sampling approximation SFD scheme with the fixed difference operator lengths. Figure 12 shows the final PP and PS images migrated by the sampling approximation SFD scheme with the adaptive variable difference operator lengths.  50191 s, respectively. The application of adaptive variable difference operator lengths can improve computational efficiency compared with that of fixed difference operator lengths under the same imaging accuracy.

Field data application
In the final example, we use a field data set from the northeast China to test elastic RTM algorithm based on the sampling approximation SFD scheme with adaptive variable difference operator lengths, and compare it with that based on Taylor-series expansion SFD scheme. Figure 13 displays the migrated PP and PS images using the Taylor-series expansion SFD scheme with fixed difference operator lengths. Figure 14 displays the migrated PP and PS images using the sampling approximation SFD scheme with adaptive variable difference operator lengths. As shown in Figs. 13 and 14, two elastic RTM algorithms produce similar images, except in the areas marked by the black arrows. Some fuzzy events in the PP and PS images generated by the conventional elastic RTM algorithm based on Taylor-series expansion SFD scheme are visible (highlighted by the black arrows in Fig. 13a and b), which reduces the seismic imaging accuracy. In contrast, the events in the PP and PS images (highlighted by the black arrows in Fig. 14a and b) generated by the proposed elastic RTM algorithm are clearer and more focused. The field data test shows that the proposed RTM algorithm can produce better images than the conventional algorithm. Furthermore, this field data test also shows computational efficiency of the RTM with the adaptive variable b) a) a) b) difference operator lengths improves by about 28%, compared with that with the fixed difference operator lengths.

CONCLUSIONS
We have put forward and used an efficient SFD scheme based on a sampling approximation method with adaptive variable difference operator lengths to implement elastic prestack RTM. The sampling approximation SFD scheme is used to solve elastic wave equations for forward and backward wavefield extrapolations. Numerical dispersion analysis and the wavefield extrapolation results show that the sampling approximation SFD scheme has greater accuracy than the Taylor-series expansion SFD scheme. We tested the elastic prestack RTM algorithm based on the sampling approximation SFD scheme with adaptive variable difference operator lengths. Experiments on two theoretical models and a field data set demonstrate that the elastic prestack RTM using the sampling approximation SFD scheme can generate better images than that using the Taylor-series expansion SFD scheme, particularly for PS images. Compared with the conventional fixed difference operator lengths, the application of the adaptive variable difference operator lengths can effectively improve computational efficiency in elastic prestack RTM.