Lowrank seismic wave extrapolation on a staggered grid

Introduction

Wave extrapolation in time is an essential part of seismic imaging, modeling, and full waveform inversion. Finite-difference methods (Etgen, 1986; Wu et al., 1996) and spectral methods (Tal-Ezer et al., 1987; Reshef et al., 1988) are the two most popular and straightforward ways of implementing wave extrapolation in time. The finite-difference (FD) methods are highly efficient and easy to implement. However, the traditional FD methods are only conditionally stable and suffer from numerical dispersion (Finkelstein and Kastner, 2007). Thanks to advances in supercomputer technology, spectral methods have become feasible for large-scale problems. Compared with FD methods, spectral methods have superior accuracy and are able to suppress dispersion artifacts (Etgen and Brandsberg-Dahl, 2009).

Several spectral methods have been developed for seismic wave extrapolation in variable-velocity media (Fowler et al., 2010; Du et al., 2014). Zhang et al. (2007) and Zhang and Zhang (2009) proposed a one-step extrapolation algorithm, which is derived from an optimized separable approximation (OSA). This algorithm formulates the two-way wave equation as a first-order partial differential equation in time without suffering from numerical instability or dispersion problems. Soubaras and Zhang (2008) proposed a two-step extrapolation method that is based on a high-order differential operator, which allows for large time steps in extrapolation. However, the decomposition algorithm in OSA can be expensive, particularly in anisotropic media. Fomel et al. (2013b) presented an approach to approximating the mixed-domain operator using a lowrank decomposition, which reduces computational cost by optimally selecting reference velocities and weights. Song and Fomel (2011) developed a related method, Fourier finite-differences (FFD), by cascading a Fourier transform operator and a finite-difference operator to form a chain operator.

In practice, first-order wave equations are often involved in handling wave extrapolation in media with both velocity and density variations. Mast et al. (2001) provided a derivation of the $k$ -space method for solving the ultrasonic wave equation. Tabei et al. (2002) extended this method to solving coupled first-order differential equations for wave propagation, efficiently accounting for velocity and density variations. In the $k$ -space method, dispersion errors from a second-order time integration operator are compensated by a modified spectral operator in the wavenumber domain. This correction is exact for a medium with constant velocity in particular. Song et al. (2012) modified the $k$ -space method with a mixed-domain operator and applied FFD to hand these operators. This method has highly accurate for variable velocity and density.

In the FD methods, the FD coefficients are conventionally determined through a Taylor-series expansion around the zero wavenumber (Dablain, 1986; Kindelan et al., 1990). Traditional FD methods are therefore particularly accurate for long-wavelength components. Several approaches have been proposed to improve the performance of FD method in practice. Implicit FD operators (Chu and Stoffa, 2011; Liu and Sen, 2009) can be used to achieve high numerical accuracy. Another way to control numerical errors is to use optimized FD operators (Takeuchi and Geller, 2000; Liu, 2013; Chu et al., 2009). Song et al. (2013) derived optimized coefficients of the FD operator from a lowrank approximation (Fomel et al., 2013b) of the space-wavenumber extrapolation matrix.To improve the accuracy and stability, the FD methods have been developed on a staggered grid (Levander, 1988; Virieux, 1986,1984; Madariaga, 1976). Moczo et al. (2002) investigated the stability and grid dispersion in the 3D fourth-order staggered grid FD scheme. In the past years, the viscous wave modeling using staggered grid FD methods (Operto et al., 2007; Robertsson et al., 1994; Bohlen, 2002)have also been studied and reported.

In this paper, we use modified staggered grid $k$ -space method (Tabei et al., 2002; Song et al., 2012) to handle the derivative operator in mixed-domain for variable velocity and density. We introduce lowrank decomposition (Fomel et al., 2013b) to approximate the modified k-space extrapolation operator and reduce the computational cost. Inspired by lowrank finite-differences (Song et al., 2013), we derive optimized finite-difference coefficients for coupled first-order wave-propagation equations using staggered spatial and temporal grids (Levander, 1988; Virieux, 1986,1984). We apply dispersion analysis and use the method of manufactured solutions to evaluate the accuracy of the proposed methods. Numerical tests demonstrate that the proposed SGL(staggered grid lowrank) and SGLFD(staggered grid lowrank finite-differences) methods are highly accurate and applicable for variable velocity and density modeling and reverse time migration (RTM) in complicated models. Our implementation of the new methods and the numerical examples are based on Madagascar software (Fomel et al., 2013a) and can be reproduced using the latest version of Madagascar.

Lowrank seismic wave extrapolation on a staggered grid