We consider the problem of constructing a wave extrapolation operator in a variable and possibly anisotropic medium. Our construction involves Fourier transforms in space combined with the help of a lowrank approximation of the spacewavenumber wave-propagator matrix. A lowrank approximation implies selecting a small set of representative spatial locations and a small set of representative wavenumbers. We present a mathematical derivation of this method, a description of the lowrank approximation algorithm, and numerical examples which confirm the validity of the proposed approach. Wave extrapolation using lowrank approximation can be applied to seismic imaging by reverse-time migration.
Wave extrapolation in time plays an important role in seismic imaging (reverse-time migration), modeling, and full waveform inversion. Conventionally, extrapolation in time is performed by finite-difference methods (Etgen, 1986). Spectral methods (Tal-Ezer et al., 1987; Reshef et al., 1988) have started to gain attention recently and to become feasible thanks to the increase in computing power. The attraction of spectral methods is in their superb accuracy and, in particular, in their ability to suppress dispersion artifacts (Chu and Stoffa, 2008; Etgen and Brandsberg-Dahl, 2009). Theoretically, the problem of wave extrapolation in time can be reduced to analyzing numerical approximations to the mixeddomain space-wavenumber operator (Wards et al., 2008). In this paper, we propose a systematic approach to designing wave extrapolation operators by approximating the space-wavenumber matrix symbol with a lowrank decomposition. A lowrank approximation implies selecting a small set of representative spatial locations and a small set of representative wavenumbers. The optimized separable approximation or OSA (Song, 2001) was previously employed for wave extrapolation (Zhang and Zhang, 2009; Du et al., 2010) and can be considered as another form of lowrank decomposition. However, the decomposition algorithm in OSA is significantly more expensive, especially for anisotropic wave propagation, because it involves eigenfunctions rather than rows and columns of the original extrapolation matrix. Our algorithm can also be regarded as an extension of the interpolation algorithm of Etgen and Brandsberg-Dahl (2009), with optimally selected reference velocities and weights. Another related method is the Fourier finite-difference (FFD) method proposed by Song and Fomel (2010). FFD may have an advantage in efficiency, because it uses only one pair of forward and inverse Fast Fourier Transforms per time step. However, it does not offer flexible controls on the approximation accuracy.
The algorithm does not require, at any step, access to the full matrix W, only to its selected rows and columns. Once the decomposition is complete, it can be used at every time step during the wave extrapolation process.
Our next example (Figure 3) corresponds to wave extrapolation in a 2-D smoothly variable isotropic velocity field. As shown by Song and Fomel (2010), the classic finite-difference method tends to exhibit dispersion artifacts with the chosen model size and extrapolation step, while spectral methods exhibit high accuracy. The wavefield snapshot (Figure 5) confirms the ability of our method to handle complex models and sharp velocity variations.
