Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation

Introduction

Wave extrapolation in time is crucial in seismic modeling, imaging (reverse-time migration), and full-waveform inversion. The most popular and straightforward way to implement wave extrapolation in time is the method of explicit finite differences (FDs), which is only conditionally stable and suffers from numerical dispersion (Wu et al., 1996; Finkelstein and Kastner, 2007). In practice, a second-order FD for temporal derivatives and a high-order FD for spatial derivatives are often employed to reduce dispersion and improve accuracy. FD coefficients are conventionally determined using a Taylor-series expansion around zero wavenumber (Dablain, 1986; Kindelan et al., 1990). Therefore, traditional FD methods are accurate primarily for long-wavelength components.

More advanced methods have been applied previously to FD schemes in the case of one-way wave extrapolation (downward continuation). Holberg (1988,1987) designed the derivative operator by matching the spectral response in the wavenumber domain. Soubaras (1996) adopted the Remez exchange algorithm to obtain the $L_{\infty}$-norm-optimized coefficients for second-derivative filters. Mousa et al. (2009) designed stable explicit depth extrapolators using projections onto convex sets (POCS). These approaches have advantages over conventional FD methods in their ability to propagate shorter-wavelength seismic waves correctly. To satisfy the general criterion for optimal accuracy (Geller and Takeuchi, 1995), Geller and Takeuchi (1998) derived an optimally accurate time-domain finite difference method for computing synthetic seismograms for 1-D problems extended later to 2-D and 3-D (Takeuchi and Geller, 2000). Liu and Sen (2009) proposed FD schemes for two-way scalar waves on the basis of time-space dispersion relations and plane-wave theory. Later on, they suggested adaptive variable-length spatial operators in order to decrease computing costs significantly without reducing accuracy (Liu and Sen, 2011). The Liu-Sen scheme satifies the exact dispersion relation and has greater accuracy and better stability than a conventional one. However, it still uses an expansion around the zero wavenumber.

In sedimentary rocks, anisotropic phenomena are often observed as a result of layering lithification, which is described as transversely isotropic (TI). Tectonic movement of the crust may rotate the rocks and tilt the natural vertical orientation of the symmetry axis (VTI), causing a tilted TI (TTI) anisotropy. Wavefields in anisotropic media are well described by the anisotropic elastic-wave equation. However, in practice, seismologists often have little information about shear waves and prefer to deal with scalar wavefields. Conventional P-wave modeling may contain shear-wave numerical artifacts in the simulated wavefield (Zhang et al., 2009; Grechka et al., 2004; Duveneck et al., 2008). Those artifacts as well as sharp changes in symmetry-axis tilting may introduce severe numerical dispersion and instability in modeling. Yoon et al. (2010) proposed to reduce the instability by introducing elliptical anisotropy in regions with rapid tilt changes. Fletcher et al. (2009) suggested to include a finite S-wave velocity in order to enhance stability when solving coupled equations. These methods can alleviate the instability problem; however, they may alter the wave propagation kinematics or leave residual S-wave components in the P-wave simulation. A number of spectral methods are proposed to provide solutions which can completely avoid the shear-wave artifacts (Cheng and Kang, 2012; Fowler and Lapilli, 2012; Song and Fomel, 2011; Etgen and Brandsberg-Dahl, 2009; Chu and Stoffa, 2011; Fomel et al., 2012; Zhan et al., 2012) at the cost of several Fourier transforms per time step. These methods differ from conventional pseudo-spectral methods (Gazdag, 1981; Fornberg, 2002), because they approximate the space-wavenumber mixed-domain propagation matrix instead of a Laplacian operator.

Our goal is to design an FD scheme that matches the spectral response in the mixed space-wavenumber domain for a wide range of spatial wavenumbers. The scheme is derived from the lowrank approximation of the mixed-domain operator (Fomel et al., 2010,2012) and its representation by FD with adapted coefficients. We derive this kind of FD schemes which we call lowrank FD or LFD for both isotropic and TTI media. Using this approach, we only need to compute the FD coefficients once and save them for the whole process of wave extrapolation or reverse-time migration. The method is flexible enough to control accuracy by the rank of approximation and by FD order selection.

The paper is organized as follows. We first give a brief review of the lowrank approximation method. As a spectral method, it provides an accurate wave extrapolation, but it is not optimally efficient. Next, we present the derivation of LFD. LFD as an FD method can reduce the cost and is also more adaptable for parallel computing on distributed computer systems. We also propose lowrank Fourier FD (LFFD), by replacing the original FD operator in the two-way Fourier FD (FFD) (Song and Fomel, 2011) with the corresponding LFD. LFFD improves the accuracy of FFD, in particular in tilted transversely isotropic (TTI) media. A number of synthetic examples of increasing complexity validate the proposed methods. In this paper, we solve the acoustic wave equation in constant-density media, aiming at incorporating wave extrapolation with LFD and LFFD into seismic imaging by reverse-time migration. It is possible to extend LFD and LFFD to variable-density media by factoring the second-order k-space operator into first-order parts (Tabei et al., 2002; Song et al., 2012).

Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation