Layout 6


1

ANNALS OF GEOPHYSICS, 63, 3, DM329, 2020; doi:10.4401/ag-8197

A staggered-grid lowrank finite-difference 
method for elastic wave extrapolation 
Qizhen Du1,2,3, Jing Ba*,4, Dong Han5, Pengyuan Sun6, Jianlei Zhang6 
 
(1) Key laboratory of deep oil and gas, China University of Petroleum (East China), Changjiang West Road 66th, Qingdao, 

China, 266580 
(2) Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao, 

China, 266071 
(3) Key Laboratory of Geophysical Prospecting, CNPC, China University of Petroleum (East China) Changjiang West Road 

66th, Qingdao, China, 266580 
(4) School of Earth Sciences and Engineering, Hohai University, Nanjing, Jiangsu, China, 211100 
(5) Sinopec Geophysical Research Institute, Nanjing, Jiangsu, China, 211103  
(6) China National Petroleum Corp Bureau of Geophysical Prospecting Inc. Zhuozhou, Hebei, China, 072750 
 
Article history: received May 8, 2019; accepted December 18, 2019 
 
 
 

Abstract  
 
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. 
 
 
Keywords: Lowrank finite-difference, Lowrank approximation, Elastic wavefield extrapolation. 
 
 
 

 



1. 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 P-
wave 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., 2012, 
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; Han et al., 2016]. 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-difference-
based 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; Liu and 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 [Fang et al., 2014]. 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 first-
order 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 [Fomel et al., 2013; Cheng et al., 2016] and the lowrank FD method [Song et al., 2013; Fang et al., 2014; Han 
et al., 2016]. 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 mixed-
domain 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.

2



[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 S-
wavefield 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. 

 
 

2. Methodology 
 
2.1 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 second-
order 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 S-
wave 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 

(3) 
 

 
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 

𝜕𝑡²𝜕²𝐮 ‒𝑐   𝛻(𝛻 � 𝐮) + 𝑐   𝛻 × (𝛻 × 𝐮) = 0² � ² �

𝜕𝑡²𝜕²𝐔� +𝑐  𝐤 �𝐤 � 𝐔�� ‒ 𝑐   𝐤 × �𝐤 × 𝐔�� = 0² � ² �

𝐔� =            ∫𝐮 � 𝑒‒𝑖𝐤�𝐱d𝐱 .1 (2𝜋)²

3

Lowrank elastic wave extrapolation



Qizhen Du et al.

4

           𝐔�p = 𝐤� (𝐤� � 𝐔�) (4) 
       𝐔�𝑠 = ‒𝐤� × (𝐤� × 𝐔�) (5) 
 

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: 
             𝐔� = 𝐔�p + 𝐔�𝑠  (6) 

 
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) characterize the P- and S-wave propagations, respectively. The derivation details are given in 
Appendix A. 

According to the plane wave theory, analytical solutions to equations (7) and (8) are expressed as 
    𝐔�p = 𝐤� (𝐤� � 𝐔�₀)𝑒±𝑖𝑐�𝐤𝑡 (9) 
and                 𝐔�𝑠 = ‒ 𝐤� × (𝐤� × 𝐔�₀)𝑒±𝑖𝑐�𝐤𝑡 . (10) 
 
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:  
 
             𝐔�𝑝 (𝑡 +𝛥𝑡)+ 𝐤�𝑝 (𝑡‒𝛥𝑡) ‒ 2𝐔�𝑝 (𝑡)=𝐤 ��𝐤� · 𝐔�(𝑡)� · �2(cos(|𝑘|𝑐𝑝𝛥𝑡)‒1)�    (11)          𝐔�𝑠 (𝑡 +𝛥𝑡)+ 𝐤�𝑠 (𝑡‒𝛥𝑡) ‒ 2𝐔�𝑠 (𝑡)=‒ 𝐤 �×�𝐤� × 𝐔��(𝑡)� · �2(cos(|𝑘|𝑐𝑠𝛥𝑡)‒1)� (12) 
 
 
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:  

 
              . (13) 

 
 
Indices 𝑗,𝑙 = 1,3 denote x- and z-axes of the Cartesian coordinate in a 2D medium and the mixed-domain terms 

𝐤 |𝐤|

𝜕²𝐔�p 𝜕𝑡² +𝑐   𝐤 �𝐤 � 𝐔�� = 0² �
𝜕²𝐔�𝑠 𝜕𝑡² ‒𝑐   𝐤 ×�𝐤 × 𝐔�𝑠� = 0 .² �

� 𝑘�𝑗𝑘�𝑙 �2(cos(|𝑘|𝑐𝑝𝛥𝑡)‒1)�𝑘�𝑗𝑘�𝑙 �2(cos(|𝑘|𝑐𝑠𝛥𝑡)‒1)�



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. 

 
 
2.2 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 1986]. 
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 first-
order system.  

If we define 𝜈� and 𝜏�� as the particle velocity and stress components, respectively, the velocity-stress formulas 
[Leaney, 1988] follow that: 
 
 
  (14) 
 

 
 
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 ‒ 
2sin2(𝑥/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. 

⎧ ⎪  ⎨ ⎪ ⎩

= 1 𝜌(𝐱) �𝜏�� �𝑥��𝜈� �𝑡
�𝜏�� �𝑥� = 𝜆(𝐱)𝛿��          +𝜇(𝐱)�        +        �

,�𝑣� �𝑥� �𝑣� �𝑥��𝑣� �𝑥� �𝑣� �𝑥�

⎧ ⎪  ⎨ ⎪ ⎩

�𝑣�𝑙 �𝑡 1 𝜌= 
�𝜏��� �𝑥� = 𝜆𝛿�� 𝑖𝑘𝑚𝜈�𝑚 + 𝜇 (𝑖𝑘𝑗𝜈�𝑙 + 𝑖𝑘𝑙𝜈�𝑗�

,
𝑖𝑘𝑗𝜏���

⎧  ⎨ ⎩ ² �𝑖²𝑘𝑗𝑘𝑙 𝑐   𝛥𝑡 sinc²�|𝐤|𝑐�𝛥𝑡� 
.   ²

² �𝑖²𝑘𝑗𝑘𝑙 𝑐   𝛥𝑡 sinc²�|𝐤|𝑐�𝛥𝑡�   ²

�|𝐤|𝑐�,� 𝛥𝑡� ₂

𝑖𝑘𝑙sinc²�|𝐤|𝑐�,�𝛥𝑡� .     ²

5

Lowrank elastic wave extrapolation



These variables follow the linear principle, 𝜈𝑗 = 𝜈�𝑗 + 𝜈�𝑗. The second-order decoupled wave equations are transformed 
to the first-order ones (Li et al., 2007):  
 
 
 

(18) 
 
 

 
 
 
 
 
 

(19) 
 
 

 
 

 
 

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. 

�
�

𝜕𝑧𝜕𝜏�𝜌1=𝜕𝑡𝜕𝜈�
𝜕𝑥𝜕𝜏�𝜌1=𝜕𝑡𝜕𝜈�

�𝜕𝑧𝜕𝜈�+𝜕𝑥𝜕𝜈��= (𝜆 +2𝜇)𝜕𝑡𝜕𝜏�⎧ ⎪⎪  ⎨ ⎪⎪ ⎩
𝜕𝜈�𝜕𝑡 = 1𝜌 � 𝜕𝜏��𝜕𝑧 + 𝜕𝜏��𝜕𝑥 �

𝜕𝜏��𝜕𝑡 = −2𝜇 𝜕𝜈�𝜕𝑧
=

𝑠 𝑠 𝑠
𝜕𝜈�𝜕𝑡 = 1𝜌 � 𝜕𝜏��𝜕𝑧 + 𝜕𝜏��𝜕𝑥 �

𝑠 𝑠 𝑠
𝜕𝜏��𝜕𝑡 = 𝜇 � 𝜕𝜈�𝜕𝑧 + 𝜕𝜈�𝜕𝑥 � ,

𝑠
𝑠

𝜕𝜏��𝜕𝑡 −2𝜇 𝜕𝜈�𝜕𝑥
𝑠

⎧ ⎪⎪⎪⎪⎪  ⎨ ⎪⎪⎪⎪⎪ ⎩

~𝜈�(𝑡)‒𝜈� (𝑡‒𝛥𝑡)~𝛥𝑡 1𝜌 𝑖sinc(c�|k|𝛥𝑡/2)𝑘�𝜏�~~ ~
=

𝛥𝑡

� � =
1𝜌 𝑖sinc(c�|k|𝛥𝑡/2)𝑘�𝜏�~𝜈�(𝑡)‒𝜈� (𝑡‒𝛥𝑡)� �𝛥𝑡~ ~𝜏�(𝑡+𝛥𝑡/2)‒𝜏�(𝑡‒𝛥𝑡/2) = (𝜆+2𝜇)𝑖sinc(c�|k|𝛥𝑡/2)(𝑘�𝑣�+𝑘�𝑣�)~ ~

⎧ ⎪⎪⎪  ⎨ ⎪⎪⎪ ⎩
~��=�� ~𝑖sinc(c�|k|𝛥𝑡/2)�𝑘�𝜏��+𝑘�𝜏���𝜌1𝛥𝑡

~𝜈�(𝑡)‒𝜈� (𝑡‒𝛥𝑡)~
~��=�� ~𝑖sinc(c�|k|𝛥𝑡/2)�𝑘�𝜏��+𝑘�𝜏���𝜌1𝛥𝑡

~𝜈�(𝑡)‒𝜈� (𝑡‒𝛥𝑡)~
~= ~𝜇𝑖sinc(c�|k|𝛥𝑡/2)�𝑘�𝑣�+𝑘�𝑣�� ,��𝜏��(𝑡+𝛥𝑡/2)‒𝜏��(𝑡‒𝛥𝑡/2)

~~
𝛥𝑡

= ~‒2𝜇𝑖sinc(c�|k|𝛥𝑡/2)𝑘�𝑣���𝜏��(𝑡+𝛥𝑡/2)‒𝜏��(𝑡‒𝛥𝑡/2)
~~

𝛥𝑡
= ~‒2𝜇𝑖sinc(c�|k|𝛥𝑡/2)𝑘�𝑣���𝜏�(𝑡+𝛥𝑡/2)‒𝜏�(𝑡‒𝛥𝑡/2)

~~
𝛥𝑡

⎧ ⎪⎪⎪⎪⎪⎪  ⎨ ⎪⎪⎪⎪⎪⎪ ⎩

Qizhen Du et al.

6



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 𝑐�(𝐱) [Fang, et al., 2014]. Applying reverse Fourier transforms to equations 
(20) and (21) yields:  
 
 
 
 
 
 
 
 
 
 

(22) 
 
 
 
 
 
 
 
 
 
 

 
 
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 (Fang et al., 2014] 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 
[Fomel et al., 2013; Fang et al., 2014]. In the next section, we will demonstrate the details of the lowrank 
approximation scheme and the lowrank finite-difference method. 

��
��

𝜈𝑧 = 𝜈𝑧 + 𝜈𝑧
𝜈𝑥 = 𝜈𝑥 + 𝜈𝑥�

⎧ ⎪⎪  ⎨ ⎪⎪ ⎩ � ,𝜕𝑧‒
𝜕�𝑣�+𝜕𝑥‒𝜕�𝑣��= �𝜆(𝐱)+2𝜇(𝐱)�𝜏�(𝑡+𝛥𝑡/2)‒𝜏�(𝑡‒𝛥𝑡/2)𝛥𝑡

𝜕𝑧+𝜕𝜏𝜌= 𝜌(𝐱)1�� 𝛥𝑡𝜈�(𝑡+𝛥𝑡/2)‒𝑣� (𝑡‒𝛥𝑡/2)
𝜕𝑥+𝜕𝜏𝜌=�� 𝜌(𝐱)1𝛥𝑡𝜈�(𝑡+𝛥𝑡/2)‒𝜈� (𝑡‒𝛥𝑡/2)

⎧ ⎪⎪⎪⎪⎪  ⎨ ⎪⎪⎪⎪⎪ ⎩

�� �𝜕𝑧+𝜕�𝜏��+𝜕𝑥+𝜕�𝜏���= 𝜌(𝐱)1�� 𝛥𝑡𝜈�(𝑡+𝛥𝑡/2)‒𝜈� (𝑡‒𝛥𝑡/2)� � 𝜕𝑧‒𝜕�𝑣�= ‒2𝜇(𝐱)𝜏��(𝑡+𝛥𝑡/2)‒𝜏��(𝑡‒𝛥𝑡/2)𝛥𝑡� � 𝜕𝑥‒𝜕�𝑣�= ‒2𝜇(𝐱)𝜏��(𝑡+𝛥𝑡/2)‒𝜏��(𝑡‒𝛥𝑡/2)𝛥𝑡
𝜕𝑥+𝜕𝑧+ �𝜕𝜈�+𝜕�𝑣��� � = 𝜇(𝐱)𝜏��(𝑡+𝛥𝑡/2)‒𝜏��(𝑡‒𝛥𝑡/2)𝛥𝑡

�� �𝜕𝑥‒𝜕𝜏��+𝜕𝑧‒𝜕�𝜏���=�� 𝜌(𝐱)1𝛥𝑡𝜈�(𝑡+𝛥𝑡/2)‒𝜈� (𝑡‒𝛥𝑡/2)

⎧ ⎪⎪⎪⎪  ⎨ ⎪⎪⎪⎪ ⎩

˖𝜕𝑥�𝜕� = 𝐹‒1�𝑖𝑘�𝑒+𝑖𝑘�𝛥�/2 sinc�c�(𝐱)|k|𝛥𝑡/2�𝐹{·}�
𝜕𝑥�̠𝜕� = 𝐹‒1�𝑖𝑘�𝑒−𝑖𝑘�𝛥�/2 sinc�c�(𝐱)|k|𝛥𝑡/2�𝐹{·}�

˖𝜕𝑥�𝜕� = 𝐹‒1�𝑖𝑘�𝑒+𝑖𝑘�𝛥�/2 sinc�c�(𝐱)|k|𝛥𝑡/2�𝐹{·}�
˗𝜕𝑥�𝜕� = 𝐹‒1�𝑖𝑘�𝑒−𝑖𝑘�𝛥�/2 sinc�c�(𝐱)|k|𝛥𝑡/2�𝐹{·}�

7

Lowrank elastic wave extrapolation



2.3 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 [Fomel et al., 2013; 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 𝑊𝑗� (𝐱,𝐤) = 𝑘𝑗 sinc(𝑐𝑝(𝐱) ǀ𝐤ǀ 𝛥𝑡 / 2) which is an 𝑁� by 𝑁� matrix, 𝐱 = {𝑥1, ..., 𝑥𝑁�} represent samples 
in the space domain and 𝐤 = {𝑘1, ..., 𝑘𝑁�} represent samples in the wavenmber domain. The exact operator, 𝑊𝑗� (𝐱,𝐤), 
could be factorized into the product of three submatrices, which has the form as 
 
 

(24) 
 

 
where 𝑊𝑗� (𝐱,𝐤�) is a submatrix of the original matrix 𝑊𝑗� (𝐱,𝐤) consisting of M columns selected from 𝑊𝑗� (𝐱,𝐤) 
related to 𝐤�, while 𝑊𝑗� (𝐱�,𝐤) is another submatrix of the original matrix 𝑊𝑗� (𝐱,𝐤) consisting of N rows selected from  𝑊𝑗� (𝐱,𝐤) related to 𝐱�. 𝑎�� represents the middle coefficients which have M-by-N entries, and M and N stand for the 
rank of the decomposition. The way to construct the decomposition representations follows the method by Engquist 
and Ying [2007, 2009]. The computational complexity of the algorithm is 𝑂(𝑁𝑁� log𝑁�), and the computational cost 
of the Fourier transform is 𝑁� log𝑁�. M and N are usually far less than the spatial number, 𝑁�. The ranks are up to 
the complexity of the model. For homogeneous media, the mixed-domain operators are equivalent to the pseudo 
spectral operators. There is an actual trade-off in the selection of ranks: larger values lead to a more accurate wave 
representation but require a longer computational time. In Fomel et al. [2013], the authors select the ranks based 
on an estimate of the approximation accuracy of 10−4. For the isotropic heterogeneous media, ranks between 2 and 
4 can lead to satisfying accuracy [Fang et al. 2014]. In this paper, we recommend 𝑀, 𝑁 = 4. If we utilize equation (24), 
the approximate form of   could be given by 
 
 

(25) 
 

 
2.4 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 

𝜕�𝑣� 𝜕𝑥 ̟�

𝑊𝑗� (𝐱, 𝐤) ≈�   � 𝑊𝑗� (𝐱, 𝐤�)𝑎��𝑊𝑗� (𝐱�, 𝐤),� �=1
� �=1

𝜕�𝑣� 𝜕𝑥 ̟�
𝜕�𝑣�𝜕𝑥�̟ ≈�   � 𝑊𝑗� (𝑥, 𝑘�)𝑎��𝐹‐1�𝑖𝑒𝑖𝑘�𝛥�/2𝑊𝑗� (𝐱�, 𝐤)𝐹[𝜈�(𝐱,𝑡)]�.� �=1

� �=1

Qizhen Du et al.

8

Figure 1. A schematic of a 2D staggered grid scheme. 



submatrix and obtain the staggered-grid lowrank finite-difference. We further decompose the third term on the 
right-hand side of equation (25) into: 
 
 

(26) 
 

 
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 
 

(29) 

 
 

(30) 
 

 

(31) 
 

 
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). 

Figure 2b, 2c and 2d illustrate the exact matrix with the k-space method, the approximate matrix with the 
ESGLFD method and the absolute error between these two matrices, respectively. There is almost no difference 
between Figures 2b and 2c at the macroscopic level. By inspection, the absolute error, shown as Figure 2d, is almost 
constrained on the interval (-0.04, 0.04). The error is tolerant in the processing of elastic wave extrapolation. 

Compared to the existing works, the elastic lowrank finite-difference operators are optimal finite-difference 
operators and constructed in accordance with the spectral response in certain range of frequency. 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. It extends the application of lowrank 
finite-difference schemes from acoustic case to elastic situation, which provides a more accurate approximation in 
both temporal and spatial domains. 

𝜕�𝑣�𝜕𝑥�̟ ≈      � 𝐺𝑗� (𝐱, 𝑙)𝐹‐1[𝑖�𝑒𝑖𝑘�𝛥��/2𝐁(𝑙, 𝐤)�𝐹(𝜈�)]1 2
≈� 𝐺𝑗� (𝐱, 𝑙)𝐹‐1[𝑖(𝑒𝑖𝑙𝑘�𝛥�� ‒ 𝑒‒𝑖(𝑙‒1)𝑘�𝛥��)𝐹(𝜈�)].
≈      � 𝐺𝑗�(𝐱, 𝑙)𝐹‐1[𝜈�(𝑥� + 𝑙)‒𝜈�(𝑥� ‒ (𝑙‒1)]1 2

� �=1� �=1 � �=1

𝑊𝑗� (𝐱�, 𝐤) ≈ � 𝐶𝑗� (𝐱�, 𝑙)𝐵 (𝑙, 𝐤),� �=1
2𝑙‒1 ₂ 

𝐺𝑗� (𝐱, 𝑙) =�   � 𝑊𝑗� (𝐱, 𝐤�)𝑎��𝐶𝑗� (𝐱�, 𝑙).� �=1
� �=1

𝜕�𝑣�𝜕𝑥𝑗 ≈      � 𝐺𝑗�(𝐱, 𝑙)𝐹‐1[𝜈�(𝑥� + 𝑙 ‒ 1)‒𝜈�(𝑥� ‒ 𝑙)]1 2 � �=1
𝜕�𝑣�𝜕𝑥𝑗 ≈      � 𝐺𝑗�(𝐱, 𝑙)𝐹‐1[𝜈�(𝑥� + 𝑙)‒𝜈�(𝑥� ‒ 𝑙 + 1)]1 2 � �=1
𝜕�𝑣�𝜕𝑥𝑗 ≈      � 𝐺𝑗�(𝐱, 𝑙)𝐹‐1[𝜈�(𝑥� + 𝑙 ‒ 1)‒𝜈�(𝑥� ‒ 𝑙)] .1 2 � �=1

9

Lowrank elastic wave extrapolation



2.5 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: 
 

(32) 
 

and 
 
                                                                        

,                                                               (33) 

where,  
 𝜂�� = 𝐺𝑥� (𝐱, 𝑙)sin((2𝑙‒1)𝑘�𝛥ℎ / 2),  

 
 𝜂�� = 𝐺𝑧� (𝐱, 𝑙)sin((2𝑙‒1)𝑘�𝛥ℎ / 2),  

 
 

 � 
� �=1

𝜉𝑝 = 2arcsin (�(𝜂𝑝𝑥 + 𝜂𝑝𝑧)𝜈�𝛥𝑡²/4)𝑘𝛥𝑡𝜈�² ² ²

𝜉𝑠 = 2arcsin (�(𝜂𝑠𝑥 + 𝜂𝑠𝑧)𝜈�𝛥𝑡²/4)𝑘𝛥𝑡𝜈�² ² ²

 � 
� �=1

Qizhen Du et al.

10

Figure 2. (a) A linearly increasing P-wave velocity model, (b) An exact matrix for 𝜕� / 𝜕�𝑗 with the k-space method, (c) An 
approximation matrix for 𝜕� / 𝜕�𝑗 with the ESGLFD method, (d) The absolute error between the exact matrix and 
the approximate matrix.

a) b)

c) d)



𝜂�� = 𝐺𝑥� (𝐱, 𝑙)sin((2𝑙‒1)𝑘�𝛥ℎ / 2),    
and                                                                𝜂�� = 𝐺𝑧� (𝐱, 𝑙)sin((2𝑙‒1)𝑘�𝛥ℎ / 2).  

 
 

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 finite-
difference (SGFD) scheme and the ESGLFD scheme in terms of difference orders, velocity, time intervals and 
propagation angles. 

Figure 3 depicts the relation between 𝜉 and 𝑘ℎ for different orders of SGFD and ESGLFD methods. Clearly, 
when  𝐿 increases, the dispersion error decreases for both SGFD and ESGLFD methods. However, the proposed 
ESGLFD method covers a wider wavenumber range than the SGFD method. The dispersion curves are given in 
Figure 4. Note that the errors change significantly with the velocities when employing the SGFD method. By 
comparison, the errors equal or approximately equal the limit value on the interval [0,2.2] when the ESGLFD 
method is used. The ESGLFD method keeps highly accurate on the wavenumber range mentioned above. 

From Figure 5, we can estimate the impact of time intervals on the dispersion. The dispersion error by the 
SGFD method becomes larger when the time interval increases. However, as shown in Figure 5b, all the 
dispersion errors are confined in a small neighborhood of 1 with 𝑘ℎ∊[0,2.0]. It proves that the ESGLFD method 
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. 

 

 � 
� �=1

 � 
� �=1

11

Lowrank elastic wave extrapolation

Figure 3. 2D dispersion curves of the FD methods (a) and the ESGLFD methods (b) for different orders, fourth-order (𝐿 = 2, blue), six-order (𝐿 = 3, red), eight-order (𝐿 = 4, purple), sixteen-order (𝐿 = 8, green), 𝜈 = 3000m/s, D𝑡 = 1ms 
and 𝛳 = 𝜋/8.

a) b)



Qizhen Du et al.

12

Figure 4. 2D dispersion curves of the FD methods (a) and the ESGLFD methods (b) for different velocities, 𝜈 = 2500 m/s 
(blue), 𝜈 = 3500 m/s (red), 𝜈 = 4500m/s (purple), 𝜈 = 5000m/s (green),  𝐿 = 6, D𝑡 = 1ms and 𝛳 = 𝜋/8.

a) b)

Figure 5. 2D dispersion curves of the FD methods (a) and the ESGLFD methods (b) for different time intervals, D𝑡 = 0.5ms 
(red), D𝑡 = 1.0ms (blue), D𝑡 = 1.5ms (purple), D𝑡 = 2.0ms (green), 𝜈 = 3000m/s, 𝐿 = 6 and 𝛳 = 𝜋/8.

a) b)

Figure 6. 2D dispersion curves of the FD methods (a) and the ESGLFD methods (b) for different propagation angles, 𝛳 = 𝜋/16 (blue),𝛳 = 𝜋/8 (red), 𝛳 = 𝜋/6 (purple), 4 (green), D𝑡 = 1ms, 𝜈 = 3000m/s and 𝐿 = 6.
a) b)



3. Numerical examples 
 
3.1 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/cm3, 1700m/s and 1000m/s, while the parameters of the lower layer are 2.2g/cm3, 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.  

Figure 9 displays the vertical slices extracted from the x-component snapshots shown in Figure 8. Notice that there 
are three waveforms in Figures 9a and 9b and two waveforms in Figures 9c and 9d. It is in accordance with the fact that 
the reflected P-wave, the transmission P-wave, and the direct P-wave exist in the snapshots, as is shown in Figures 8a 
and 8b, while the reflected converted S-wave and the transmission converted S-wave exist in the snapshots of Figures 
8c and 8d. Meanwhile, the results also demonstrate that the wave modes, which occur in the elastic media, could be 
depicted by the proposed ESGLFD method effectively. The waveforms obtained by the SGFD method do not preserve 
the wavelet shapes well. Vibration curves occur in the neighborhood of the Ricker-shaped waveforms. Comparatively, 
the waveforms generated by the ESGLFD method do comply with the Ricker-wavelet forms. The results prove that the 
ESGLFD method is accurate and effective for elastic wave simulation when the time interval becomes larger. 

 
 
3.2 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, subset of the elastic Marmousi2 model denoted by the dashed rectangle is used for elastic wave 
simulation in this experiment. The four corners are all located at the border of the resampled Marmousi2 model, the 

13

Lowrank elastic wave extrapolation

Figure 7. The geometry and the elastic parameters of a two-layer model. The pentagram denotes the excitation position 
for elastic wave simulation.



Qizhen Du et al.

14

Figure 8. Elastic wavefield snapshots at t=0.8s of a two-layer model by using the 4th-order SGFD method and the 4th-order 
ESGLFD method: (a) the x-component snapshot of the P-wave by the SGFD method; (b) the x-component snapshot 
of the P-wave by the ESGLFD method; (c) the x-component snapshot of the S-wave by the SGFD method; (d) the 
x-component snapshot of the S-wave by the ESGLFD method. Head waves are highlighted by the red arrows.

a) b)

c) d)

Figure 9. Vertical slices extracted from the x-component snapshots (a) related to P-wave by the SGFD method, (b) related 
to P-wave by the ESGLFD method, (c) related to S-wave by the SGFD method, and (d) related to S-wave by the 
ESGLF D method at the distance of 3.0 km.

a) b) c) d)



exact coordinates of which are (6.0, 0) km, (13.0, 0) km, (6.0, 3.5) km, and (13.0, 3.5) km. Both the horizontal and the 
vertical space increments are 10 m. The partial model is determined with a grid size of 700 × 350, and the Ricker 
wavelet is used as the explosive source term in this case, with the dominant frequency of 25 Hz. The time sampling 
increment is 1 ms. The PML boundary is implemented on the four sides of the model and has a thickness of 300 m. 
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. 

 
 
3.3 Elastic Reverse Time Migration 
 
The 2D Marmousi2 model, as is shown in Figure 10, is used in this section for further RTM, in which the 

ESGLFD method serves as the forward and backward wavefield extrapolators. The geometry of the survey is from 

15

Lowrank elastic wave extrapolation

a)

b)

Figure 10. The resampled (a) and smoothed (b) Marmousi2 model corresponding to P-wave velocity for elastic wave 
simulation and reverse time migration, respectively. The subset denoted by the dashed rectangle is used for 
elastic wave simulation.



0 km to 17.0 km, in which a total number of the split-spread source is 340 with 700 receivers per shot. For each 
shot, one-half of the receivers are located on either side of the shot. All the sources and receivers are situated on 
the surface of the model, with the source spacing 50 m and receiver spacing 10 m. Similarly, to the elastic wave 
simulation, the Ricker wavelet with the dominant frequency of 25 Hz is utilized. Total record time is 6 s, with time 
sample increment of 1 ms. Horizontal and vertical sampling intervals are both set at 10 m. Additionally, PML 
with the thickness of 300 m on each side of the model is implemented. With regards to image condition, we utilize 
the scalar imaging conditions [Du et al., 2017], which obtain the scalar images using the normalized cross-correlation 

Qizhen Du et al.

16

Figure 11. Snapshots corresponding to the horizontal component of the P-wave (a), the vertical component of the P-wave 
(b), the horizontal component of the S-wave (c) and the vertical component of the S-wave (d) at 0.8 s with the 
6th-order ESGLFD method. The model used for elastic wave simulation is the subset, denoted by the dashed 
rectangle, of the resampled elastic Marmousi2, shown in Figure 9.

a) b)

c) d)

Figure 12. Elastic common-source gathers corresponding to the horizontal component (a) and the vertical component (b).

a) b)



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. 

 
 

4. 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. 

17

Lowrank elastic wave extrapolation

a)

b)

Figure 13. The stacked PP image component (a), the stacked PS image component (b).



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. 

 
 

Acknowledgements  We appreciate the support of the National Science Foundation of China (41930429, 41774139) and 
the National Key Research & Development Program of China (2017YFB0202903). And we appreciate the enlightening 
discussion of our related works with Chengfeng Guo, Sanyi Yuan, Qamar Yasin, Guanchao Wang and Xingyu Zhang. 
We thank the developers of Madagascar for providing the open-source software. 
 
 
References 
 
Aki K., Richards P.G. (2002). Quantitative seismology, University Science Books (2nd Edition). 
Bérenger J.P. (1994). A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114(2), 

185–200. 
Bojarski N.N. (1985). The k-space formulation of the scattering problem in the time domain: An improved single 

propagator formulation, J. Acoust. Soc. Am. 71(2), 97 - 100. 
Candes E., L. Demanet and L. Ying (2007). Fast computation of fourier integral operators, SIAM J. Sci. Comput. 29(6), 

2464-2494.  
Cheng J.B., T. Alkhalifah, Z.D. Wu, P. Zou, C. Wang (2016). Simulating propagation of decoupled elastic waves using 

low-rank approximate mixed-domain integral operators for anisotropic media, Geophysics 81(2), T63-T77. 
Compani-Tabrizi B. (1986). K-space scattering formulation of the absorptive full fluid elastic scalar wave equation 

in the time domain, J. Acoust. Soc. Am. 79(4), 901-905. 
Du Q.Z, Y.T. Zhu, J. Ba (2012). Polarity reversal correlation for elastic reverse time migration, Geophysics 77( 2), 

S31-S41. 
Du Q.Z., X.F. Gong, M.Q. Zhang, Y.T., G. Fang (2014). 3D PS-wave imaging with elastic reverse-time migration, 

Geophysics 79 (5), S173-S184. 
Du Q.Z., D. Han, S.A. Hou (2016). Lowrank Finite-difference Method for Elastic Wave Propagation, 78th EAGE 

Conference & Exhibition, Vienna, A.T., Extended Abstracts. 
Du Q, Guo C F, Zhao Q, et al. (2017). Vector-based elastic reverse time migration based on scalar imaging condition, 

Geophysics 82(2), S111-S127. 
Du X., P.J. Fowler, R.P. Fletcher (2013). Recursive integral time-extrapolation methods for waves: A comparative 

review, Geophysics 79(1), T9-T26. 
Engquist B., L. Ying (2007). Fast directional multilevel algorithms for oscillatory kernels, SIAM J. Sci. Comput. 29(4), 

1710-1737. 
Engquist B., L. Ying (2009). A fast-directional algorithm for high frequency acoustic scattering in two dimensions, 

Comm. Math. Sci. 7(2), 327-345. 
Etgen J., S. Brandsberg-Dhl (2009). The pseudo-analytical method: Application of Pseudo-Laplacians to acoustic 

acoustic and acoustic anisotropic wave propagation, 79th SEG Meeting, Houston, USA, Expanded Abstracts, 
2552-2556. 

Fang G., S. Fomel, Q.Z. Du, J. Hu (2014). Lowrank seismic wave extrapolation on a staggered grid, Geophysics 79(3), 
T157-T168. 

Finkelstein B., R. Kastner (2007). Finite-difference time domain dispersion reduction schemes, J. Comput. Phys. 
221(1), 422-438. 

Fornberg B., and G.B. Whitham (1978). A numerical and theoretical study of certain nonlinear wave phenomena, 
Philos. Trans. R. Soc. London 289, 373-404 

Qizhen Du et al.

18



Fomel S., L. Ying, X. Song (2013). Seismic wave extrapolation using low rank symbol approximation, Geophys. 
Prospect. 61(3), 526-536. 

Granger P.Y., M. Martin, J.L. Boelle, E. Ceragioli, F. Lefeuvre, E. Crouzy (2005). Autonomous 4C nodes used in infill 
areas to complement streamer data, deepwater case study, 75th SEG Meeting, Houston, USA, Expanded 
Abstracts, 84-87. 

Han D., Q.Z. Du, M.Q. Zhang, C.F. Guo, X.F. Gong (2016). Effects of staggered-grid and rotated staggered-grid on the 
wavefield separation in VTI media, Acta Seismol. Sin. 38(1), 111-121. 

Han D., Q.Z. Du, G. Fang (2016). Elastic wave propagation on a staggered-grid using lowrank finite-difference 
method. 86th SEG Meeting, Dallas, USA, Expanded Abstracts, 3900-3904. 

Li Z.C., H. Zhang, Q.M. Liu, W.G. Han (2007). Numerical simulation of elastic wavefield separation by staggered-grid 
high-order finite-difference algorithm (in Chinese), Oil Geophys. Prospect. 42(5), 510-515. 

Liu Y., M.K. Sen (2009). A new time-space domain high-order finite-difference method for the acoustic wave 
equation, J. Comput. Phys. 228(23), 8779-8806.  

Liu Y., M.K. Sen (2011). Scalar wave equation modeling with time-space domain dispersion-relation-based 
staggered-grid finite-difference schemes, Bull. Seismol. Soc. Am. 101(1), 141-159. 

Ma D.T., G.M. Zhu (2003). P-and S-wave separated elastic wave equation numerical modeling (in Chinese), Oil 
Geophys. Prospect. 38, 482-486.  

Martin G.S., R. Wiley, K.J. Marfurt (2006). Marmousi2: An elastic upgrade for Marmousi, The Leading Edge, 25: 156-
166. 

Mast T.D., L.P. Souriau, D. Liu, M. Tabei, A.I. Nachman, and R.C. Waag (2001). A k-space method for large-scale 
models of wave propagation in tissue, IEEE Trans. Ultrason Ferroelectr. Freq. Control 48, 341-354. 

Nguyen B.D., G.A. McMechan (2013). Excitation amplitude imaging condition for prestack reverse-time migration, 
Geophysics 78(1), S37-S46. 

Pestana R.C., P. L. Stoffa (2010). Time evolution of the wave equation using rapid expansion method, Geophysics 
75(4), T121-T131. 

Pourjavid S., O.J. Tretiak (1992). Numerical solution of the direct scattering problem through the transformed 
acoustical wave equation, J. Acoust. Soc. Am. 91(2), 639-45. 

Reshef M., D. Kosloff, M. Edwards, C. Hsiung (1988). Three-dimensional acoustic modeling by the Fourier method, 
Geophysics 53(9), 1175-1183. 

Robertsson J.O., I. Moore, M. Vassallo, A.K. , D.J. Manen van, A.  (2008). On the use of multicomponent streamer 
recordings for reconstruction of pressure wavefields in the crossline direction, Geophysics 73(5), A45-A49. 

Song X., S. Fomel, L. Ying (2013). Lowrank finite-differences and low-rank Fourier finite-differences for seismic 
wave extrapolation, Geophys. J. Int. 193, 960-969. 

Soubaras R., Y. Zhang (2008). Two-step explicit marching method for reverses time migration. 70th EAGE Conference 
& Exhibition incorporating SPE EUROPEC, Rome, I.T., 2214-4609.  

Stewart R.R., J.E. Gaiser, R.J. Brown, D.C. Lawton (2003). Converted-wave seismic exploration: Applications, 
Geophysics 68(1), 40-57. 

Stoffa P.L., R.C. Pestana (2009). Numerical solution of the acoustic wave equation by the rapid expansion method 
(REM)-A one step time evolution algorithm, 79th SEG Meeting, Houston, USA, Expanded Abstracts, 2272-2276. 

Sun J., S. Fomel, L. Ying (2015). Low-rank one-step wave extrapolation for reverse time migration, Geophysics 81(1), 
S39-S54. 

Sun R., J.D. Chow, K.J. Chen (2001). Phase correction in separation P-and S-waves in elastic data, Geophysics 66(5), 
1515-1518. 

Tabei M., T.D. Mast, R.C. Waag (2002). A k-space method for coupled first-order acoustic propagation equations, J. 
Acoust. Soc. Am. 111(1), 53-63. 

Tan S.R., L.J. Huang (2014). An efficient finite-difference method with high-order accuracy in both time and space 
domains for modelling scalar-wave propagation, Geophys. J. Int. 197, 1250-1267.  

Virieux J. (1984). SH-wave propagation in heterogeneous media: velocity-stress finite-difference method, 
Geophysics 49(11), 1933-1942. 

Virieux J. (1986). P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method, 
Geophysics 51(4), 889-901. 

Voynov O.Y., V.I. Golubev, M.S. Zhdanov, I.B. Petrov (2018). Migration of elastic wavefield using adjoint operator and 

19

Lowrank elastic wave extrapolation



Born approximation. Innovations in Wave Modelling and Decision Making, SIST Series, Volume 90, Springer 
Switzerland, Chapter 8, P. 219-240. 

Wang W., and G.A. McMechan (2015). Vector-based elastic reverse-time migration, Geophysics 80, S245-S258. 
Wards B., G. Margrave, M. Lamoureux (2008). Phase-shift time-stepping for reverse-time migration. 78thSEG 

Meeting, Las Vegas, USA, Expanded Abstracts, 2262-2266. 
Yan J., P. Sava (2008). Isotropic angle-domain elastic reverse-time migration, Geophysics 73(6), S229-S239. 
Yan J., P. Sava (2009). Elastic wave-mode separation for VTI media, Geophysics 74(5), WB19-WB32. 
Yan R., X.B. Xie (2012). An angle-domain imaging condition for elastic reverse time migration and its application 

to angle gather extraction, Geophysics 77(5), S105-S115.  
Yuan S., S. Wang, W. Sun, L. Miao, Z. Li (2014). Perfectly matched layer on curvilinear grid for the second-order 

seismic acoustic wave equation, Explor. Geophys. 45, 94-104. 
Zhang Q., G.A. McMechan (2010). 2D and 3D elastic wavefield vector decomposition in the wavenumber domain for 

VTI median, Geophysics 75(3), D13-D26. 
Zhang Y., G. Zhang (2009). One-step extrapolation method for reverse time migration. Geophysics 74(4), A29-A33. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

*CORRESPONDING AUTHOR: Jing BA, 
School of Earth Sciences and Engineering,  

Hohai University, Nanjing,  

Jiangsu, China, 211100,  

e-mail: jba@hhu.edu.cn 

© 2020 the Istituto Nazionale di Geofisica e Vulcanologia. 

All rights reserved

Qizhen Du et al.

20