A filtering approach to remove finite-difference errors from wave equation simulations

Abstract

Seismograms (i.e., recordings of seismic waves that propagate through the earth) can be used to uncover information about the earth's subsurface. Such investigations require accurate numerical wave simulations. One of the most common techniques to carry out these simulations is the finite-difference (FD) method. In the FD method, (1) derivatives are replaced with approximations of limited accuracy, and (2) continuous space and time are discretized into finite steps. The FD method is a fast numerical method, but it also introduces inaccuracies. In this thesis, we propose four procedures that reduce these inaccuracies. The overarching aim is to provide fast FD simulations (using large steps in space and time) while yielding accurate solutions. The first proposed method is the use of a filter pair: the forward and inverse time-dispersion transforms. These transforms must be applied before the simulation (to modify the source wavelet) and after the simulation (to modify the recorded seismic signals). They correct for the inaccuracy induced by the approximation of the temporal derivative in the wave equation. We show that the method applies to acoustic and elastic wave simulations. Furthermore, we show that the method applies to viscoelastic FD simulations if they use standard memory variables. The second proposed method is the use and design of `optimal' FD operators. Such FD operators are highly accurate for a prescribed wavenumber range. We obtain these FD operators using the Remez exchange algorithm, a well-known algorithm in the field of filter design. Our work generalizes the existing literature drastically: (1) we consider arbitrary derivative orders, (2) we consider three cost-functions [the absolute error, the relative error, the group velocity error], (3) we consider arbitrary input locations, (4) we can compute solutions that are optimal in a least-squares or maximum norm sense. Optimal results in FD modeling are obtained with the FD operator designed for the relative error. The third proposed method concerns the implementation of point-sources in FD simulations. These sources are typically modeled by exciting the source on a single FD node. We show that such an implementation leads to wavenumber-varying amplitude errors. In effect, two artifacts are generated: (1) ringing is introduced, and (2) erroneous wave modes may be excited. We show how to correct for this error using a filter in the wavenumber domain. The `FD-consistent' point-source that we propose minimizes these artifacts. The fourth proposed method concerns the use of interface representation schemes in wave simulations. For this, we compare five interface representations from geophysical literature. We find that, in acoustic simulations, optimal results are obtained with anti-aliasing of the fine velocity model. Conversely, in isotropic and anisotropic elastic simulations, optimal results are obtained with the Schoenberg \& Muir (1989) calculus. The proposed methods have two attractive features: (1) they allow FD simulations with large steps in space and time, (2) they must only be applied before and after the simulation to improve the accuracy, and have a negligible computational cost. Hence, they allow for fast FD simulations with a minimal computational cost, while yielding excellent accuracy.