A staggered-grid lowrank finite-difference method for elastic wave extrapolation

Elastic wave extrapolation in the time domain is significant for an elastic wave equation-based processing. To improve the simulation reliability and accuracy of decoupled elastic P- and S- waves, we propose the staggered-grid lowrank finite-difference method based on the elastic wave decomposition. For elastic wave propagation, a lowrank finite-difference method based on the staggered grid is derived to improve the accuracy. Regarding the application of the decoupled elastic wave equation, we derive the finite-difference scheme coefficients which are dependent on velocity. Based on the elastic wave decomposition and plane wave theories, we formulate the elastic wave-extrapolation operators, which contain trigonometric adjustment factors. Accordingly, by applying the lowrank method to approximating the operators, the finite- difference scheme is designed to discretize the decoupled wave equation. The derivation processing implies the combination of elastic wave-mode decomposition and extrapolation. The proposed method enables elastic P- and S-waves to extrapolate in the time-space domain separately and produces accurate P-and S-wave components simultaneously. Dispersion analysis suggests that our proposed method is reliable and accurate in a wide range of wavenumber. Numerical simulation tests on a simple model and the Marmousi2 model validate the accuracy and effectiveness of the method, showing its ability in handling complex structures. Although the operators are accurate only when the medium is homogeneous, they are of high accuracy when the velocity gradient is quite small and are applicable when the velocity gradient is large. The subsequent results of reverse time migration for the Marmousi2 model also suggest that the proposed method is enough to serve as an extrapolator in elastic reverse time migration.


Introduction
Over the past few decades, seismic exploration technique has been extended to the elastic wavefield. This remarkable stride benefits from the development of multi-component seismic acquisition technique [Granger et al., 2005;Robertsson et al., 2008], which provides the accessibility to not only the primary P-wave but also the S-wave.
A great deal of information from diverse wave modes exist in multicomponent data, which portray the prominent subsurface structure and discriminate the lithological variation as compared to conventional single-component Pwave data [Stewart et al., 2003;Yan and Xie, 2012]. In such case, S-wave becomes an important factor to calibrate the bright-spot and to estimate the fractures [Du et al., 2012].
The processing based on the elastic wave equations act as a bridge between multicomponent data and elegant profiles. The subsurface wavefield consists of the P-and S-waves. For a physical meaning of P-and S-waves, the processing is commonly observed to contain either the wavefield separation [Yan and Sava, 2008;Du et al., 2012Du et al., , 2014 or the wavefield decomposition [Zhang and McMechan, 2010;Wang and McMechan, 2015]. As for the separation in isotropic media, classic curl and divergence operators are efficient and practical [Aki and Richards, 2002;Yan and Sava, 2009;. However, Sun et al. [2001] showed that they induce amplitude distortion and phase shift. Such problems can be fully avoided by using the decomposition operators [Zhang and McMechan, 2010] or the decoupled wave equations [Ma and Zhu, 2003;Li et al., 2007;Du et al., 2017]. The decomposition method can fully retain the amplitude and phase of the P-and S-waves. The decomposition operators are applied to the subsurface wavefield and are not related to the wavefield extrapolation. Therefore, the decoupled wave equations need improved extrapolation methods.
The elastic wavefield extrapolation algorithm is the core of processing such as elastic wave simulation, reverse time migration (RTM), and full waveform inversion. Those techniques require an efficient and accurate algorithm, in which the large time step is a key to the efficiency. Conventional finite-difference methods [Virieux, 1984] and pseudo-spectral methods [Reshef et al., 1988] are the two popular methods of wave extrapolation. However, conventional finite different methods are constrained to the Courant-Friedrch-Lewy (CFL) condition and require small time steps to guarantee stability [Finkelstein and Kastner, 2007]. The pseudo-spectral methods, equivalent to infinite order finite-difference, require even smaller time step. A solution to this issue lies in compensating error of approximating the temporal derivatives approximation in spatial derivative computation. For finite-differencebased methods, the optimized coefficients of finite-difference in joint time-space domain could become a good choice to improve the accuracy [Tan and Huang, 2014;Sen,2009, 2011]. For the spectral-based method, incorporation of the temporal compensation factors into the extrapolation operators in the mixed domain is a desirable way to reduce the artifacts caused by the temporal dispersion [Soubaras and Zhang,2008;Zhang and Zhang, 2009;Stoffa and Pestana, 2009;Etgen and Brandsherg-Dahl, 2009;Du et al., 2013].
The k-space method is another type of the spectral domain wave extrapolation method. It is essentially an extension of the pesudo-spectral method and reaches high accuracy by using a temporal iteration procedure in combination with the accurate spectral evolution of spatial derivatives [Tabei, et al., 2002]. The dispersion errors of the second-order time integration operators are compensated by this procedure . Early attempts of the k-space method include Fornberg and Whitham [1978]'s application to nonlinear wave equations and Bojarski [1982]'s research on the scattering problem. It is then extended by many researches including the application to the elastic scalar wave equation [Compani-Tabrizi, 1986] and the application to the acoustical wave equation [Pourjavid and Tretiak, 1992]. The first-order wave equations are capable of handling velocity and density variations and are easier to apply the perfectly matched layer (PML) boundary condition. Hence, k-space method on the basis of firstorder wave equations attracts many attentions [Mast et al, 2001;Tabei et al., 2002;Fang et al., 2014].
However, due to its velocity-dependency [Etgen and Brandsberg-Dahl, 2009], solving the k-space operators in the spectral domain involves high computational cost. To reduce computational costs, the flexible-control lowrank methods established a great interest toward further investigation which contains the lowrank decomposition method Cheng et al., 2016] and the lowrank FD method Fang et al., 2014;. The lowrank decomposition method generates a sparse separable decomposition of the operator matrix and balance the accuracy and efficiency through the rank. As a branch of the lowrank decomposition method, Sun et al. [2015] developed a lowrank one-step method. Cheng et al. [2016] proposed a series of approximate mixeddomain integral operators in the anisotropic media for elastic wave extrapolation. The lowrank FD method (LFD) is now widely used to design the FD scheme in relation to the spectral response in the mixed-domain. Song et al.

Qizhen Du et al.
[2013] proposed the LFD method for the propagation of the scalar wave. It is extended to a coupled first-order system to handle the variable density and velocity by Fang et al. [2014]. However, the current literature on the LFD method exerts their efforts to the scalar wave extrapolation based on the acoustic assumption and ignores the elastic nature of the earth. It is essential to develop a new LFD-based method to depict the structures of the subsurface from elasticity and linearity standpoints.
In this paper, our prime goal is to design a lowrank FD method for elastic wavefield extrapolation based on decoupled wave equations. The LFD-based elastic wave extrapolation is available to provide accurate P-and Swavefield vectors. Firstly, we begin with the common elastic displacement equations to derive second-order extrapolation operators in the mixed-domain and deduce first-order operators from these second-order operators.
Then, for heterogeneous media, we formulate a discretization form of the first-order decoupled system. Next, the lowrank approximation is applied in the first-order extrapolation operators and the corresponding staggered-grid LFD scheme is designed. Finally, to validate the proposed method, we implement the elastic wave simulation on the simple models and on the complex ones, respectively. We also apply our proposed method to elastic RTM on the Marmousi2 model.

Second-order Elastic Wave Extrapolation Operators in Mixed-domain
The LFD-based elastic wave extrapolation matters to the construction of the decoupled P-and S-wavefield. In consideration of the decoupled wave equations, we derive second-order extrapolation operators to obtain secondorder elastic wave extrapolation in the mixed-domain. According to the elastic wavefield vector decomposition method [Zhang and McMechan, 2010], we could obtain the decoupled elastic wave equations from the common second-order elastic wave equations [Aki and Richard, 2002]. Given that the velocity terms are constant, we employ the analytical solutions to these decoupled elastic wave equations and formulate the second-order extrapolation operators.
In a 2D infinite, homogeneous, and isotropic medium, the elastic wave equation is [Aki and Richards, 2002]: (1) where = ( , ) denotes a displacement vector, represents the time, is the P-wave particle velocity, is the Swave particle velocity, and = ( , ) is the spatial location vector. Note that in this paper, the discussions and numerical tests are in 2D cases. Given that the velocity terms are constant, we could derive the following equations in the wavenumber domain: (2) where = ( , ) represents the vector of wavenumber. The vector of represents the displacement vector in the domain of wavenumber, and the forward Fourier transform is expressed as Equations (1) and (2) are the coupled elastic wave formulas in which the displacement vector contains not only the P-wave but also the S-wave. For a physical interpretation of seismic processing results (e.g., seismic imaging profiles), it is necessary to implement wavefield separation or decomposition. The derivation of decoupled wave equations is based on the condition that P-and S-waves propagate independently in homogeneous media. The decomposed P-and S-wavefields [Zhang and McMechan, 2010] could be expressed as where = is the unit vector of wavenumber, p and represent the P-and S-wave vectors in the wavenumber domain, respectively. Under the linear assumption, the vector wavefield is expressed as the superposition of P-wave (irrotational) part and S-wave (non-divergence) part: Substituting the vector wavefield expressed by equation (6) into equation (2), and applying the identities in connection with the divergence and the curl, the decoupled wave equations can be given as (7) and (8) equations (7) and (8)  According to the plane wave theory, analytical solutions to equations (7) and (8) are expressed as Herein, ₀ is the amplitude vectors. The sign ± denotes that the solution stands for the outgoing or the ingoing wavefield.
We merge the analytical solutions into the temporal discretization equation of the leapfrog form and obtain the mixed-domain integral solutions for the elastic case as follows: Note that the equation related to P-wave only describes the P-wave propagation, which means that if the terms related to S-wave are added to the P-wave equation, the products of the interaction between these terms and the propagation operators are zeros. Therefore, under the condition that no effects of the S-wave exist in the P-wave formulas, we can replace the P-wavefield symbol with the wavefield symbol at the right-hand side of equation (11). Similarly, we can determine the S-wave propagation equation as shown by equation (12).
equations (11) and (12) are the exact discretization representations of the second-order time derivatives. When describing equations (11) and (12) at the component level, we could obtain the mixed-domain terms: .
Indices , = 1,3 denote xand z-axes of the Cartesian coordinate in a 2D medium and the mixed-domain terms are the kernels of the mixed-domain recursive time integration operators [Pestana and Stoffa, 2010;Du et al., 2013] for the elastic case. We use the Taylor expansions to approximate these exact mixed-domain operators and find that the pseudo-spectral operators, in the elastic case, are only the second-order Taylor polynomial approximations of these operators. In this sense, it is the high-order terms of the mixed-domain operators' Taylor expansion polynomials that compensate the temporal accuracy loss caused by the second-order polynomials.

First-order Elastic Wave Extrapolation Operators in Mixed-Domain
In practice, the wavefield extrapolation methods based on the first-order velocity-stress formulas catch higher ability to handle a variant density situation than the second-order displacement formulas [Virieux, 1984 and.
It is easier to combine the PML boundary conditions [Bérenger, 1994;Yuan et al., 2014] with the first-order system of elastic wave equations [Tabei et al., 2002]. Therefore, we extend the method demonstrated above to the firstorder system.
If we define and as the particle velocity and stress components, respectively, the velocity-stress formulas [Leaney, 1988] follow that: where ( ) and ( ) are the lamé coefficients, and ( ) is density. The subscripts ,, represent the components of the Cartesian coordinate system, and is the Kronecker delta.
For homogeneous media, the formulas in the time-wavenumber domain are determined as (15) where i is the imaginary unit, and are Lame's constants, and , are the 2D spatial forward Fourier transforms of particle velocity and stress components, respectively.
To derive the first-order operator in the mixed-domain, we employ the trigonometric identities cos( ) = 1 -2sin 2 ( /2) and sin( ) = sinc( )/ to reshape operators (13), then we obtain (16) Note that, in the equations above, replacing the pseudo spectral operator ² with ² ² , ² sinc² improves the accuracy of the numerical simulation. In a first-order system, the pseudo-spectral operator is changed into the following form at the large time-step without lowering the accuracy (17) However, the velocity and stress components in equation (14) are not only related to P-wave but also to S-wave.
We introduce the variables = / and = / as P-and S-wave particle velocity components, respectively.
These variables follow the linear principle, = + . The second-order decoupled wave equations are transformed to the first-order ones (Li et al., 2007): where is P-wave stress, and are S-wave normal stresses and is S-wave shear stress.
We use the forward 2D spatial Fourier transform to derive the representations in the wavenumber domain (20)   (21) where , , and are the Fourier transform of , , , and . After that, we will design the decoupled first-order system of the elastic wave on staggered grids. Figure 1 describes a schematic of the 2D staggered grid scheme. We use and to denote the space grid intervals along the x-and z-axes, respectively, and i and j to denote the grid cell indices, respectively. , and are discretized at space locations ( , ), is discretized at space locations, (( +0.5) , ( +0.5) ) and are discretized at space locations (( +0.5) , ), and and are discretized at space locations ( ,( +0.5) ). In the discussion, the staggered-grid schemes are used for theoretical analyses and numerical tests.

Qizhen Du et al.
In the wavenumber domain, the staggered concepts conveyed by Figure 1 are achieved with the shift properties of Fourier transform. When the gradient of velocity is small, we can replace constant velocities, and , with the space-dependent velocities, ( ) and ( ) . Applying reverse Fourier transforms to equations (20) and (21) yields: where (23) These derivative operators are named as the first-order k-space operators [Tabei, 2002]. Note that equations (23) are actually the approximations of the exact ones which use constant velocities in the sinc operators. Its accuracy is still relatively high  in media with small velocity gradient. And it also allows us to testify this method in the complex model with large velocity gradient.
For numerical solutions to the spatial derivatives in equation (22), we could introduce the localized Fourier transform methods [Wards, 2008;Etgen and Brandsberg-Dhl, 2009;Zhang and Zhang, 2009]. However, these methods have access to all nodes in the mixed-domain and require several multidimensional inverse and forward Fourier transforms. Let stand for the total number of the spatial grids, the computational complexity of the localized Fourier transform methods amounts to ( ² ). From the mathematical perspective, the mixed-domain operators usually exhibit an oscillatory behavior and approximate lowrank [Candes et al., 2007;Engquist and Ying, 2009]. Therefore, one choice for reducing the computational complexity is the application of the lowrank methods Fang et al., 2014]. In the next section, we will demonstrate the details of the lowrank approximation scheme and the lowrank finite-difference method.

Lowrank Approximation to First-order operators in the Heterogeneous Case
It is not desirable to extrapolate the elastic wave using the recursive time integration operators because of the high computational cost. In this section, we will derive the decomposition representations of the mixed-domain operators with the lowrank decomposition method Cheng et al., 2016]. As equation (23) expressed, the mixed-domain operators can be categorized into two groups. One is associated with P-wave and the other is associated with S-wave. For a brief description, we will take as an example to introduce the approximation. Let 4 can lead to satisfying accuracy . In this paper, we recommend , = 4. If we utilize equation (24), the approximate form of could be given by

Lowrank Finite-difference Schemes on Staggered Grids for the Elastic Case
In this section, we will further decompose the recursive time integration operators and generate optimal coefficients for elastic wave extrapolation. The submatrix ( , ) is only related to the wavenumber and exhibits the behavior of oscillation [Engquist and Ying, 2009]. Thus, we could adopt Fourier sine series to represent this ˖ submatrix and obtain the staggered-grid lowrank finite-difference. We further decompose the third term on the right-hand side of equation (25) into: where ( , ) is an by submatrix and determined as sin( ), L is the half size of the staggered grid stencil, ( , ) is an by submatrix calculated as the product of ( , ) and the pseudo-inverse of ( , ). We introduce the auxiliary matrix (27) After that, we could overwrite the spatial derivative of as (28) equation (28) is analogous to the one expressed by the standard staggered grid scheme. It is a striking feature of the method that the deduction of the coefficients containing lowrank decomposition and the coefficients depend upon the media. This method could be entitled as an elastic staggered-grid lowrank FD(ESGLFD) method. The coefficients for the different media are supposed to be stored before the implementation of elastic wave simulation because of the velocity-dependency. Through the same processing, the other spatial derivatives of , containing / ˗ , / ˖ , / ˗ , could be deduced as Herein, we use a 2D simple homogeneous media to analyze the accuracy of the ESGLFD method. Also, we choose the spatial derivative operator, / ˗ as an example. The spatial derivative operator, / ˗ shown in equation (22), is calculated by the k-space method [Tabei et al., 2002] and the ESGLFD method, respectively. The velocity is given by ( , ) = 20 + 25 + 1200. The velocity value increases to 3900 m/s (Figure 2a).

Dispersion Analysis
We analyze the accuracy of the ESGLFD method. According to the plane wave theory, we could acquire a set of ratios containing the ratio of the P-wave numerical velocity to the real P-wave velocity, , and that of the S-wave, .
These ratios indicate the errors: When = 1, no dispersion exists. If keeps a long way off 1, serious dispersion will occur. Also, it is clear that the error depends on the propagation angle. Suppose that the propagation angle is in a 2D case, let = cos( ) and = sin( ). Due to the symmetry of equations (32) and (33), we only investigate the properties of the error when ∊[0, /4]. Thereafter, we will compare the difference between the staggered grid finitedifference (SGFD) scheme and the ESGLFD scheme in terms of difference orders, velocity, time intervals and propagation angles. is effective for large time-step simulation. As mentioned above, we only investigate the dispersion relation with ∊[0, /4] and select the representative angles including /16, /8, /6 and /4 to illustrate the variation of the dispersion errors. From Figures 6a and 6b we can see that the trend of the dispersion curves by the two methods for the designed propagation angle is comparable. It is the wider coverage of wavenumber that gives the advantage of the ESGLFD method over the SGFD method. After the analysis of these Figures, we can find that the ESGLFD method exhibits higher accuracy. Especially, the ESGLFD method is able to yield the promise results when the time intervals become larger.

Two-Layers Model
To investigate the validity of the ESGLFD method, we implement an experiment on a two-layer model (Figure 7).
Size of the model is defined as , with horizontal and vertical sampling intervals equal 15m. The time sampling interval is 2.0ms. The contrast of the model appears at the depth of 4500m. The density, P-and S-wave velocities of the upper layer are set as 2.0g/cm 3 , 1700m/s and 1000m/s, while the parameters of the lower layer are 2.2g/cm 3 , 2550m/ s and 1500m/ s, respectively. A Ricker wavelet is utilized as an explosive source at the location (4500m, 3975m), with 35 Hz the dominant frequency.
The x-component snapshots for P-and S-waves at 0.8 s are shown in Figure 8. Visible and apparent dispersion occurs in Figures 8a and 8b. When time and space intervals become larger, it is obvious that the waveform simulated by the SGFD method is distorted. The ESGLFD method performs well and yields nearly dispersion-free results.

The Marmousi2 Model
To validate our proposed method at the complex model, we select the modified elastic Marmousi2 model [Martin et al., 2006], and set the grid size of the model as 1700 × 350. The grid interval is redefined as 10 m. We divide the resampled P-wave velocity by 3 to calculate the S-wave velocity. For the elastic wave simulation and the subsequent RTM, we smooth the velocity model. And Figure 10a and 10b describe the original and smoothed P-wave velocities, respectively.
In Figure 10b  The snapshots corresponding to P-and S-wave vector components at 0.8 s (as shown in Figure 11) are generated by the sixth-order ESGLFD method. The distinction between P-and S-waves can be clearly delineated with these snapshots. We could also see that the information of the S-wave is a very valuable part to depict the subtle structure of the complex velocity-variation model. Figure 12 shows two representative elastic common-source gathers corresponding to horizontal and vertical components, respectively. The direct waves have been removed and we can see that the reflection waves are without any visible numerical dispersion. The experiment verifies that the ESGLFD method can handle the complex media and provide the accurate wavefield for the subsequent RTM.

Elastic Reverse Time Migration
The 2D Marmousi2 model, as is shown in Figure 10, between the source-side wavefields and the receiver-side wavefields. This method eliminates the artifacts resulting from the crosstalk between different types of wave modes and avoids the polarity reversal for PS images, which may lead to destructive interferences in the processing of the adjacent source stack.
The final migration images, including the PP image component and PS component, are described in Figure 13.
It can be observed that not only the stratoid and the flat structure but also the anticline and the faults appear clearly in all of these image profiles. Generally, these comparative image results shown in Figure 13 verify that the proposed method is sufficient to the function as an extrapolator to reconstruct the forward and the backward wavefields for elastic RTM.

Conclusions
In this work, a method of staggered-grid lowrank finite-difference is proposed to obtain a better representation of the elastic wave extrapolation. Numerical experiments on a two-layers model and the complex Marmousi2 model validate the precision and efficiency of the method. The following conclusions are achieved.
1. The basic assumption of the method is elasticity and linearity, which is composed of the dynamic and kinetic information of the primary and shear waves.
2. The method is derived on the basis of decomposition and plane wave theory, oriented to decompose the elastic wavefield to the pure-mode wavefield. The wavefield decomposition is the prerequisite of this method. 3. The elastic recursive time integration operator of the proposed method depends on the parameters of the elastic media and the simulation which can efficiently compensate the temporal accuracy loss.
4. The elastic lowrank finite-difference operators are optimal finite-difference operators and constructed in accordance with the spectral response in certain range of frequency. And the numerical tests indicate that the operators are capable of handling media with large velocity gradient.
5. The proposed method still follows the finite-difference scheme. In implementation, its efficiency will be a little lower than that of the regular staggered-grid finite-difference. That is because more memory space is needed to store the space-variant finite-difference coefficient and also more memory access will be introduced.