Elastic Wavefield Decomposition for Reverse-Time Migration in 3D Transverse Isotropic Media

Elastic reverse-time migration (ERTM), which utilizes the advantages of both P- and S-wave modes, is a widely used application for imaging in 3D anisotropic media. However, crosstalk due to intrinsically coupled P- and S-wavefields may degrade the image quality. To solve this problem, this study presents an effective vector P- and S-wavefield decomposition scheme in ERTM that can improve the images of 3D transversely isotropic (TI) media. The proposed method consists of four steps: (1) rotating the observation coordinate system to align its vertical axis with the symmetry axis of 3D TI media; (2) deriving the formulations of the 3D TI decomposition operator by applying the VTI P/S wave-mode decomposition strategy based on eigenform analysis in the new coordinate system; (3) implementing vector P- and S-wavefield decomposition by constructing the 3D TI Poisson equation, and introducing a novel and efficient method based on the first-order Taylor expansion to accelerate the computational efficiency of the decomposition; and (4) applying a vector-based dot-product imaging condition to generate PP and PS images. Compared with previous studies, the algorithm of our proposed method in 3D TI media is both numerically stable and computationally efficient. The 3D TI decomposition operator generates vector P- and S-wavefields showing the correct amplitude/phase with the input ones. Several numerical examples illustrate the satisfactory performance of the proposed 3D TI decomposition operator and the effective image improvement.


Introduction
Acoustic reverse-time migration (ARTM) loses shear waves by assuming zero shear velocity in the symmetry directions (Duveneck et al., 2008;Fletcher et al., 2009;Liu et al., 2019;Zhan et al., 2012;Zhang & Zhang, 2008).However, the assumption does not prevent the shear waves from propagating in other directions.It generates diamond-shaped pseudo-shear wavefronts that make the ARTM unstable, which jeopardizes the image's quality (Grechka et al., 2004;Yoon et al., 2010;Zhang et al., 2009).With the development of multi-component data acquisition technology, elastic reverse-time migration (ERTM) has been gradually applied to oil and gas exploration and production.It considers the characteristics of converted P-to-S waves which contain more information related to anisotropy (Jin et al., 2010;Weibull & Arntsen, 2014).The application of both P-and S-waves can accurately image complex reservoir structures such as gas chimneys compared with ARTM (Caldwell, 1999;Zhao & Li, 2018;Zhao et al., 2018).
However, coupled P-and S-wavefields produce crosstalk noise that adversely influences image quality.In anisotropic media, the popular elastic wavefield separation is based on the polarization directions from the solution of the Christoffel equation (Dellinger & Etgen, 1990).Non-stationary filters are used as the separator for heterogeneous vertical transverse isotropy (VTI)/tilted transverse isotropy (TTI) media to estimate the normalized polarization directions (Yan & Sava, 2009a, 2009b).Nevertheless, the separation of scalar P-waves and vector S-waves produces multiple PS images impeding its 3D application.
Another effective method based on vector elastic wavefield decomposition can obtain vector P-and S-wavefields which show consistent amplitude, phase, and physical units with original inputs (Zhang & McMechan, 2010).However, this method requires a local homogeneous model assumption for stable Fourier transforms.Low-rank approximation simplifies the vector elastic wavefield decomposition in the form of matrix multiplication in the spacewavenumber domain (Cheng & Fomel, 2014).However, the 3D computational cost of the vector elastic wavefield decomposition is prohibitively high because of multiple fast Fourier transforms (FFTs) corresponding to the selected ranks of low-rank approximation.
Recently, based on the eigenform analysis of the Christoffel equation, a pseudo-Helmholtz decomposition method is developed (Yang et al., 2019).The method mitigates the crosstalk artifacts and generates high-quality PP and PS images in 2D heterogeneous VTI/TTI media.The algorithm improves the computational efficiency by using LU factorization to solve Poisson's equation.But the pseudo-Helmholtz decomposition operators lead to incorrect amplitudes of vector P-and S-wavefields because they possess a mispositioned coefficient related to the local elastic parameters.Later, a corrected pseudo-Helmholtz decomposition operator is implemented in 3D VTI and generates the vector decomposition wavefields with the same amplitudes/phases as the input elastic wavefields (Zuo et al., 2022).However, the corrected pseudo-Helmholtz decomposition operator is not suitable for 3D TI media.Because the corrected operator overlooks the information of the tilted symmetry axis, it cannot project the P-and S-wavefields in their polarization directions.
In this paper, we propose an efficient vector elastic wavefield decomposition method for 3D TI media based on the strategy of the corrected 3D VTI pseudo-Helmholtz operator and apply the new decomposition to 3D TI ERTM.First, we review the corrected pseudo-Helmholtz decomposition method for 3D VTI media by solving the eigenvalues and eigenvectors of the Christoffel equation.Second, we derive the 3D TI decomposition operator in a rotational coordinate where its vertical axis is aligned with the tilted symmetry axis.In the new coordination, the corrected pseudo-Helmholtz decomposition for VTI can be directly applied in TI media.Then, we implement the vector P-and S-wavefield decomposition based on Poisson's equation using the derived 3D TI operator.A fast algorithm is applied to reduce the computation involved in solving Poisson's equation.Finally, the 3D TI ERTM procedure based on the derived decomposition method generates PP and PS images using a vector-based dot-product imaging condition.
The main differences between our derived decomposition method for 3D TI media and the previous work extend to 3D are: 1.The formula derivation of the 3D TI decomposition operator based on the Christoffel equation using a rotational coordinate is superior to using the Bond transform (Winterstein, 1990).2. The derived decomposition operator for 3D TI media shows the correct position of the local elastic parameters, compared with the 2D pseudo-Helmholtz decomposition operator (Yang et al., 2019).3. The proposed decomposition method for 3D TI media outperforms the previous methods (Cheng & Fomel, 2014) in terms of computational efficiency.

Review of Estimating Polarization Directions in 3D VTI
The eigenvalues and eigenvectors of the elastic wave equation represent the different wave modes (Tsvankin, 2001).The eigenvectors refer to the polarizations of P-, SV-, and SH-waves.The elastic decomposition operator based on their polarizations can generate pure P-and S-wavefields.Based on the elastic wavefield decomposition method in 3D VTI (Zuo et al., 2022), the second-order elastic wave equation in the frequency-wavenumber domain is where k x , k y , and k z are the wavenumbers, U x , U y , and U z are the three components of the particle displacement in the frequency-wavenumbers domain, q is the density, and x is the angular frequency.
where C 12 ¼ C 11 À 2C 66 .The stiffness parameters are in terms of Thomsen parameters (Thomsen, 1986): where v p and v s are the velocities of P-and S-waves along the symmetry axis.e and c refer to the differences of P and S velocities along the vertical and horizontal directions, respectively.d refers to the rate of change of P velocity along the vertical direction.Dividing both sides by the norm of wavenumber 1) is rewritten as where and x 2 k 2 refers to the square of the phase velocity for each wave mode.qx 2 and U are the eigenvalue and eigenvector of the matrix A. Applying the eigenform analysis to Eq. ( 4), the eigenvectors are where a 1 , a 2 and a 3 are the polarizations of P-, SH-, and SV-waves.In 3D VTI, the polarization directions are mutually perpendicular.The corrected pseudo-Helmholtz decomposition operator for 3D VTI can be constructed using the polarization direction of the P wave.In the space domain, the 3D VTI decomposition operator is where o x , o y , and o z are the spatial derivatives.The Detailed formulation derivation of the corrected pseudo-Helmholtz decomposition can be found in Zuo et al. (2022).
Elastic Wavefield Decomposition for Reverse-Time Migration

The Decomposition Operator for 3D TI Media
In this section, we derive the vector decomposition operator in 3D TI media.The geometrical elements shown in the observation coordinate are defined in Fig. 1b.h and u are the tilt and azimuthal angles of the anisotropic symmetry axis.Because of the symmetry of TI media, the rotated VTI media can be seen as TI.In Fig. 1c, we rotate the observation coordinate to align its vertical axis along the tilted anisotropic symmetry axis.
Based on the geometrical relation, the rotational coordinate is In the rotational coordinate, the 3D VTI elastic wave decomposition method above can be directly applied in TI.The new spatial partial derivatives are The 3D TI decomposition operator is When h ¼ u ¼ 0, the operator is simplified to VTI.For 3D isotropic media, the operator in Eq. ( 10) becomes a gradient operator, i.e., r ¼ Any elastic wavefield u can be decomposed into a curl-free u p and divergence-free u s field (Aki & Richards, 2002).For 3D TI media, with a proper decomposition operator based on the polarization of each wave mode, the decomposed vector P-and Swaves are given by where w is a virtual wavefield that satisfies Poisson's equation, and is related to vector displacements by Note that the decomposition formulations are simplified to 2017) for isotropic media.Substituting the 3D TI operator r 3d ti into Eq.( 12), Poisson's equation is given by In the frequency-wavenumber domain, Poisson's equation is where U and W are the Fourier transforms of the original u and virtual w wavefields.LU factorization is an effective approach to solving Poisson's equation (Zuo et al., 2022).But for complex anisotropic structures, the decomposition approach still suffers from expensive computational costs.Lately, a fast algorithm is proposed to solve the virtual wavefield w and improve the computational efficiency of the vector P-and S-wavefield decomposition for 2D VTI (Zhang et al., 2022).They use the first-order Taylor expansion to approximate the partial derivative term of Poisson's equation and obtain an efficient decomposition approach.In this paper, we extend the 2D fast algorithm to 3D and derive the formulation for TI media.Equation ( 14) is simplified as It is a function of the elastic parameters, the tilt angle and azimuth of the symmetry axis, and the propagation direction.Based on the weak anisotropic assumption, we estimate Eq. ( 16) using first-order Taylor expansion around e ¼ 0 and d ¼ 0. It is approximated as Substitute Eq. ( 17) into Eq.( 15), and Poisson's equation is rewritten as Transforming Eq. ( 19) into the space domain, Poisson's equation is Elastic Wavefield Decomposition for Reverse-Time Migration À u 1 À f P sin 2 h cos 2 uu 2 À f P sin 2 h sin 2 uu 3 where u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , and u 7 are the forms of U 1 , U 2 , U 3 , U 4 , U 5 , U 6 , and U 7 in the time-space domain.
We construct Poisson's equation for 3D TI media using Eq. ( 21).The approximated formulation consists of three parts: the function f P associated with the elastic parameters; the trigonometric function associated with the tilt angle and azimuth of the symmetry axis; the functions u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , and u 7 associated with the propagation direction of the elastic wavefields.The first two functions related to the elastic parameters are calculated in the space domain.The functions with propagation direction should be estimated in the frequency-wavenumber domain.The workflow of our method is summed up in Table 1.

3D TI ERTM Using the Vector Elastic Wavefield Decomposition
We introduce the vector elastic wavefield decomposition into the 3D ERTM frame.Based on the decomposed vector P-and S-wavefields, vector-based dot-product imaging conditions are used to construct PP and PS images (Zhu, 2017) Figure 2 presents a flowchart showing the modified ERTM frame for 3D TI media.The process of our proposed vector elastic wavefield decomposition is marked in the green box.

A Homogeneous Model
We consider a homogeneous 3D TI model to illustrate the performance of our proposed decomposition method.The cubic model is 2 Â 2 Â 2 km with a spatial spacing of 10 m.The elastic parameters are v p ¼ 3:25km=s, v s ¼ 1:90km=s, q ¼ 2:00g=cm 3 , e ¼ 0:22, d ¼ 0:22, c ¼ 0:00, h ¼ 30 and u ¼ 45 .The extrapolated wavefields are calculated using the second-order accuracy in time and a staggered-grid finite-difference scheme with eighth-order accuracy in space.
We set a directional source along the three coordinate axes with a frequency of 25 Hz located at the center of the model to excite elastic wavefields.Figure 3 shows the input elastic wavefields from the x, y, and z components of displacement.Decomposed vector P-and SV-waves are marked by the blue and white arrows.Figure 4 shows the decomposition using the 3D isotropic operator.Strong energies are observed for SV-and P-wave modes (black arrows) in the vector P-and SV-wavefields, respectively.Based on the local elastic parameters, the 3D VTI operator is still unable to remove the energies completely (Fig. 5). Figure 6 shows the decomposed results using the 3D TI operator based on the tilt angle and azimuth of the symmetry axis.The numerical examples illustrate that only the operator with correct elastic parameters can produce pure vector P-and S-wave modes.
To compare the amplitude of wavefields before (the black lines) and after (the red lines) decomposition using different decomposition operators, we extract a trace located at x ¼ 1000m and y ¼ 1000m (Fig. 7).The decomposed P and SV amplitudes using the isotropic operators are not correct.These traces still display visible residuals shown in the green boxes.The amplitudes of P and SV arrivals using anisotropic decomposition operators are consistent with the original ones.Only the 3D TI decomposition

Numerical Examples
We apply the decomposition method to 3D TI ERTM on a two-layer model and a complex model modified from the SEAM Arid model.The PP and PS images both perform filtering.

3D Two-Layer TI Model
The 3D TI two-layer model is discretized by a grid of 3 Â 3 Â 2 km with a spatial interval of 10m.To test the performance of the 3D TI operator in the strong shear anisotropy, we increase the value of c.The elastic parameters are shown in Table 2.The layer is located at 1500 m.The displacement source is a Ricker wavelet with a peak frequency of 20 Hz.There are 85 shot gathers and 9801 receivers are evenly distributed at a depth of 100 m.The spatial interval of shots is approximately 40 m.In the horizontal square, there are 99 receivers on each side, with a spatial sampling interval of 30 m. Input multi-component records are simulated by the exact two-layer model.The elastic parameters in the migration model are the same as that in the first layer.We consider different ERTM procedures to illustrate the performance of our proposed 3D TI operator.
Figure 5 The elastic wavefield decomposition using the 3D VTI operator.The panels are in the same setting as Fig. 4 Elastic Wavefield Decomposition for Reverse-Time Migration Figure 8 shows the multi-component simulated wavefields from one shot.The propagation of SH mode (red arrow) is elliptical because of the enhanced shear anisotropy.Figure 9 depicts the vector P-and S-wavefields decomposed by the 3D isotropic operator showing strong crosstalk artifacts (black arrows).Without considering the tilt angle and azimuth of the symmetry axis, the 3D VTI operator fails to generate clear P-and S-wave modes (Fig. 10).As expected, Fig. 11 shows the well-decomposed vector P-and S-wavefields using the proposed 3D TI operator.SH and SV wave modes are coupled in the S-wavefields.
Figure 12a and b show the PP and PS images from 16 shots using 3D isotropic ERTM.There are four shots on each side with an interval of 64 m.
Because of different phase velocities due to anisotropy, the isotropic elastic wave equation causes wrong travel times during wavefield extrapolation.The wrong reflector image is located at about 1700 m.The isotropic operator produces crosstalk as shown in Fig. 9 that causes the obvious artifacts marked by the black arrow.Horizontal reflectors have poor spatial resolution, especially for the PS image.Figure 12c and d show the results from VTI ERTM.The imaging artifacts (black arrows) are caused by the inaccurate anisotropic decomposition operator.In contrast, the 3D TI ERTM removes most of the migration artifacts and produces clear PP and PS images (Fig. 12e and f).With sufficient stacking of 85 shots (Fig. 13), the 3D VTI operators can suppress most crosstalk artifacts, but the images using the 3D Figure 6 The elastic wavefield decomposition using the 3D TI operator.The panels are in the same setting as Fig. 4 J. Zuo et al.
TI operator are clearer.The reflections from the TI ERTM are more focused at the correct locations with better illumination (yellow arrows).

3D TI Complex Model
The second model is a small 3D portion extracted from the SEAM Arid model (Oristaglio, 2012).The portion size is 2 Â 2 Â 2 km with the same spatial spacing.The spatial interval of shots is approximately The interval is 20 m.To validate our proposed method for 3D TI media, we transfer the elastic model from HTI to TI.We increase the values of anisotropic parameters to enhance the anisotropic intensity (Fig. 14).The tilt angles and azimuth of the symmetry axis are estimated from the v p and v s models, respectively (Fig. 15).This modified model is used to simulate the multi-component records and compute the 3D TI operator.The elastic parameters of the migration models are generated by smoothing the true model using a box window with a length of 200m.The elastic wavefields are simulated using the same space-time staggered-grid finite-difference solution to the elastic wave equation.
Figure 16 shows the simulated multi-component records using a Ricker wavelet with a peak frequency of 20 Hz located at the depth of 100 m.For the comparison of different modeling engines, we use both the finite-difference and finite-element methods to generate the observation records.Because of the regular size of the model, their records are similar.So we only show the records calculated by the finitedifference method.The coupled P-, SV-, and SHmode are marked by the blue, white, and red arrows, respectively.The 3D isotropic and VTI operators fail and display visible crosstalk artifacts (black arrows in Figs. 17 and 18).In contrast, the 3D TI operator is capable of removing most of the crosstalk energies and thus generates clearer P and S records (Fig. 19).There are still some residual energies, especially the SV-and SH-wavefields in the P records.The large difference between e and d is the main reason for the residual crosstalk energies.With the first-order Taylor expansion there is some mismatch between the estimated and exact eigenvalues.
Figure 20 shows the PP and PS images from 85 shot gathers via 3D ERTM with the different decomposition methods.To improve the resolution, source normalization is used for the ERTM imaging condition.Compared with the anisotropic ERTM, the images using the isotropic operator are unfocused (black arrow shown in Fig. 20a).As expected, the anisotropic ERTM using the 3D TI operator (Fig. 20e  and f) gives the best image at the locations where medium parameters change rapidly.The PP and PS stacking results are consistent with the true elastic model (yellow arrows).These images provide more structural information than VTI ERTM (Fig. 20c and  d).To compare clearly the performance of TI and VTI ERTM, we extract 2D profiles from their 3D images.In the x À z plane with y ¼ 1000m (Fig. 21), for example, at the locations of the fault, the highangle structures are clearer from both PP and PS images; at a depth of 1700 m, the horizontal layers in the PP image are more coherent.In the horizontal x À y plane at a depth of 960 m (Fig. 22), the horizontal images display more detailed structures, especially for PP images.For both 3D VTI and TI operators, the quality of the deep PS image profile is lower than that of the PP image because of the different illumination angles between the two propagation modes for the given geometry.

Discussion
The developed 3D TI operator is effective to suppress the crosstalk artifacts and improve the imaging accuracy of ERTM in 3D TI media.Our proposed decomposition method is based on the firstorder Taylor expansion of f e; d ð Þ around e ¼ 0 and d ¼ 0. For the medium with a high degree of anisotropy, there is a large difference between the exact and approximated value of f e; d ð Þ.But the misfit mainly comes from the high wavenumber areas, which contain very low reflection energy (Zhang et al., 2022).Compared with previous elastic Figure 9 The elastic wavefields decomposed using the 3D isotropic operator.The panels are in the same setting as Fig. 4 Elastic Wavefield Decomposition for Reverse-Time Migration wavefield decomposition methods, the innovations of our method include: 1. Based on the eigenform analysis of the Christoffel equation constructed by only one matrix rotation of the coordinate system, our formula derivation of the 3D TI decomposition operator is concise.Another method of constructing the Christoffel equation is using the Bond transform (Winterstein, 1990) in the observed coordinate system.The stiffness matrix of 3D TI media is calculated by the matrix multiplications between the Bond rotation matrix and the stress/strain tensors.The process of twice matrix rotations increases the difficulty of formula derivation.2. We derive the 3D TI decomposition operator based on the theory of the corrected pseudo-Helmholtz decomposition strategy in VTI (Zuo et al., 2022) and the rotational coordinate according to the tilted symmetry axis.The derived 3D TI decomposition operator shows the correct position of the coefficient related to local elastic parameters, compared with the 2D pseudo-Helmholtz decomposition operator (Yang et al., 2019).The amplitudes/phases of P-and S-wavefields using our derived 3D TI decomposition operator are consistent with the original ones.3. Our derived decomposition algorithm has the advantage of lower computational efficiency compared with other methods.The algorithm only uses one FFT and seven inverse FFTs with the cost 8N log 2 N. N is the grid point of a 3D TI model.However, the corrected pseudo-Helmholtz decomposition method for 3D VTI based on LU factorization has a cost N 3 (Zuo et al., 2022).The traditional vector elastic wavefield decomposition method based on low-rank approximation has a cost of R Â 2N log 2 N where R is the selected rank number (Cheng & Fomel, 2014).4. We also compare the computation time of the proposed 3D TI operator with other two methods.
Under the same computation environment, our proposed decomposition method requires the least time.The pseudo-Helmholtz decomposition needs the most time because of LU factorization.Take the Arid model as an example, the length of simulated multi-component records is 1.8 s with a time interval of 0.6 ms.Their computation times are shown in Table 3.

Conclusion
We present a vector elastic wavefield decomposition method for 3D TI media by solving the Christoffel equation and applying it to ERTM.By rotating the observation coordinate, we derive the formulation of the 3D TI operator related to the tilt angle and azimuth of the symmetry axis.It obtains the amplitude-preserved P-and S-wavefields.To improve the computational efficiency, we introduce a fast method to implement the vector elastic wavefield decomposition.Our algorithm has a lower cost 8N log 2 N compared with other elastic wavefield decomposition methods in 3D TI.

Figure 1
Figure 1 Schematics show the 3D VTI/TI media with an anisotropic symmetry axis and isotropic plane.a 3D VTI media in the observation coordination.b 3D TI media with a tilted angle h and azimuthal angle u. c 3D TI media in the rotational coordinate with its vertical axis along the tilted symmetry axis

Figure 2
Figure 2The flowchart of the ERTM using our proposed decomposition for 3D TI media

Figure 7
Figure 7 The comparison of displacement traces extracted from the multi-component wavefields.a, c, and e P arrivals.b, d, and f S arrivals.The input elastic wavefields are marked by black lines.The decomposed P-and S-wavefields are marked by red lines

Figure 8
Figure 8 The original 3D multiple wavefields shown at 0.5 s.The blue arrow indicates the reflected P-wave.The white and red arrows indicate the direct SV-wave and SH-wave, respectively.The panels are in the same setting as Fig. 2

Figure 10
Figure 10The elastic wavefields decomposed using the 3D VTI operator.The panels are in the same setting as Fig.4

Figure 11
Figure 11The elastic wavefields decomposed using the 3D TI operator.The panels are in the same setting as Fig.4

Figure 13
Figure 13The 3D migrations from 85 shots using different decomposition methods.a and b PP and PS images using the VTI ERTM.c and d PP and PS images using the TI ERTM procedure

Figure 15 a
Figure 15 a and b Tilted angles h and azimuthal angle u of the symmetry axis, respectively

Figure 18
Figure 18The decomposed multi-component elastic records using the 3D VTI operator.The panels are in the same setting as in Fig.4

Figure 20
Figure 20 The comparison of migrations using different decomposition methods.a and b PP and PS images using the isotropic ERTM.c and d PP and PS images using the VTI ERTM.e and f PP and PS images using the TI ERTM procedure Figure 22The 2D horizontal sections extracted from Fig.20.The panels are in the same setting as in Fig.20 . For vector-based ERTM, the equations of imaging conditions are

Table 2
The elastic anisotropic parameters for a simple two-layer TI model