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

Lowrank Approximation

The following acoustic-wave equation is widely used in seismic modeling and reverse-time migration (Etgen et al., 2009):

$$\frac{\partial^2p}{\partial t^2} = v(\mathbf{x})^2\,\nabla^2p\;,$$ (1)

where $\mathbf{x}\,=\,(x_1,x_2,x_3)$, $p(\mathbf{x},t)$ is the seismic pressure wavefield and $v(\mathbf{x})$ is the propagation velocity. Assuming a constant velocity, $v$, after Fourier transform in space, we could obtain the following explicit expression,

$$\frac{d^2\hat{p}}{dt^2} = -v^2\vert\mathbf{k}\vert^2\hat{p}\;,$$ (2)

where

$$\hat{p}(\mathbf{k},t)=\int^{+\infty}_{-\infty}{p(\mathbf{x},t)e^{i\mathbf{k}\cdot\mathbf{x}}d\mathbf{x}}\;,$$ (3)

and $\mathbf{k}\,=\,(k_1,k_2,k_3)$.

Equation 2 has an explicit solution:

$$\hat{p}(\mathbf{k},t+\Delta t) = e^{\pm i\vert\mathbf{k}\vert v\Delta t}\hat{p}(\mathbf{k},t)\;.$$ (4)

A second-order time-marching scheme and the inverse Fourier transform lead to the well-known expression (Etgen, 1989; Soubaras and Zhang, 2008) :

$$p(\mathbf{x},t+\Delta t)+p(\mathbf{x},t-\Delta t) = 2\int^{+... ...\vert v\Delta t)e^{-i\mathbf{k}\cdot\mathbf{x}}d\mathbf{k}}\;.$$ (5)

Equation 5 provides an efficient solution in the case of a constant-velocity medium with the aid of the fast Fourier transform (FFT). When velocity varies in space, equation 5 can provide an approximation by replacing $v$ with $v(\mathbf{x})$. In such a case, a mixed-domain term, $\cos(\vert\mathbf{k}\vert v(x)\Delta t)$, appears in the expression. As a result, the computational cost of a straightforward application of equation 5 is $O(N_x^2)$, where $N_x$ is the total size of the three-dimensional space grid.

Fomel et al. (2010,2012) showed that the mixed-domain matrix,

$$W(\mathbf{x},\mathbf{k}) = \cos(\vert\mathbf{k}\vert v\Delta t),$$ (6)

can be efficiently decomposed into a separate representation of the following form:

$$W(\mathbf{x},\mathbf{k}) \approx \sum\limits_{m=1}^M \sum... ... W(\mathbf{x},\mathbf{k}_m) a_{mn} W(\mathbf{x}_n,\mathbf{k}),$$ (7)

where $W(\mathbf{x},\mathbf{k}_m)$ is a submatrix of $W(\mathbf{x},\mathbf{k})$ that consists of selected columns associated with ${\mathbf{k}_m}$, $W(\mathbf{x}_n,\mathbf{k})$ is another submatrix that contains selected rows associated with ${\mathbf{x}_n}$, and $a_{mn}$ stands for the middle matrix coefficients. The construction of the separated form 7 follows the method of Engquist and Ying (2009). The main observation is that the columns of $W(\mathbf{x},\mathbf{k}_m)$ need to span the column space of the original matrix and that the rows of $W(\mathbf{x}_n,\mathbf{k})$ need to span the row space as well as possible.

Representation (7) speeds up the computation of $p(\mathbf{x},t+\Delta t)$ because

$$p(\mathbf{x},t+\Delta t) + p(\mathbf{x},t-\Delta t) = 2 \int e^{-... ...hbf{x} \mathbf{k}} W(\mathbf{x},\mathbf{k}) \hat{p}(\mathbf{k},t) d \mathbf{k}$$

$$\approx 2 \sum\limits_{m=1}^M W(\mathbf{x},\mathbf{k}_m) \left( ... ... W(\mathbf{x}_n,\mathbf{k}) \hat{p}(\mathbf{k},t) d\mathbf{k} \right) \right).$$ (8)

Evaluation of equation 8 requires N inverse FFTs. Correspondingly, the lowrank approximation reduces the cost to $O(NN_x\log N_x)$, where $N$ is a small integer, which is related to the rank of the above decomposition and can be automatically calculated at some given error level with a pre-determined $\Delta t$. Increasing the time step size $\Delta t$ may increase the rank of the approximation ($M$ and $N$) and correspondingly the number of the required Fourier transforms.

As a spectral method, the lowrank approxmation is highly accurate. However, its cost is several FFTs per time step. Our goal is to reduce the cost further by deriving an FD scheme that matches the spectral response of the output from the lowrank decomposition.

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