Lowrank seismic wave extrapolation on a staggered grid

Second-order and first-order mixed-domain methods

We consider first-order acoustic wave equations for a medium of variable velocity and density. For a lossless 2-D medium,

$$\begin{array}{l} \displaystyle \rho(\mathbf{x})\frac{\partial... ...artial t} = -\nabla \cdot \mathbf{u}(\mathbf{x},t), \end{array}$$ (1)

where $\mathbf{u}(\mathbf{x},t)$ is acoustic particle velocity with components $u_x(\mathbf{x},t)$ and $u_z(\mathbf{x},t)$ ; $p(\mathbf{x},t)$ is the acoustic pressure; $\rho(\mathbf{x})$ is density of the medium; $v(\mathbf{x})$ is seismic wave velocity of the medium; and $\mathbf{x}=(x,z)$ denotes the space location in vector coordinate.

The second-order wave equation corresponding to equation 1 is

$$\nabla\cdot \frac{1}{\rho(\mathbf{x})}\nabla p(\mathbf{x},t) - \f... ...\mathbf{x})v^2(\mathbf{x})}\frac{\partial^2 p(\mathbf{x},t)}{\partial t^2} = 0.$$ (2)

In the case of homogeneous velocity and density, equation 2 can be written in the spatial-frequency domain as

$$\frac{\partial^2 \hat{p}(\mathbf{k},t)}{\partial t^2 } = -v^2_0\mathbf{k}^2\hat{p}(\mathbf{k},t),$$ (3)

where $\hat{p}(\mathbf{k},t)$ is the 2-D spatial Fourier transform of $p(\mathbf{x},t)$ . Equation 3 has an analytical solution,

$$\hat{p}(\mathbf{k}, t+\Delta t) = e^{\pm iv_0\left\vert\mathbf{k}\right\vert\Delta t}\hat{p}(\mathbf{k},t).$$ (4)

Applying the second-order time-marching scheme leads to the k-space scheme (Tabei et al., 2002),

$$\begin{array}{c} \displaystyle \frac{\hat{p}(\mathbf{k},t+\De... ...thbf{k}\right\vert\Delta t/2)\hat{p}(\mathbf{k},t), \end{array}$$ (5)

where sinc$(x)=\sin(x)/x$ . In general, velocity and density vary in space. When both the gradient of velocity and the time step are small, replacing $v_0$ with $v(\mathbf{x})$ in equation 5 provides a new approximation. Applying inverse Fourier transform to this approximation leads to the scheme

$$\begin{array}{c} \displaystyle \frac{p(\mathbf{x},t+\Delta t)... ... t/2)\mathbf{F}\left[p(\mathbf{x},t)\right]\right], \end{array}$$ (6)

where $\mathbf{F}$ denotes a spatial Fourier transform. The operator on the right-hand side of equation 6 depends on both $\mathbf{x}$ and $\mathbf{k}$ . Following Tabei et al. (2002) we call it second-order $kx$ -space operator,

$$\left[\nabla^{v(\mathbf{x})\Delta t}\right]^2 p(\mathbf{x},t)\equ... ...t\mathbf{k}\right\vert\Delta t/2)\mathbf{F}\left[p(\mathbf{x},t)\right]\right],$$ (7)

where the operator $\left[\nabla^{v(\mathbf{x})\Delta t}\right]^2$ is analogous to the standard gradient operator, but it is a function of parameter $v(\mathbf{x})\Delta t$ . Similar to the definition of the standard gradient operator, we can define

$$\left[\nabla^{v(\mathbf{x})\Delta t}\right]^2p(\mathbf{x},t) \equ... ...\frac{\partial}{\partial^+z}\frac{\partial}{\partial^-z}\right)p(\mathbf{x},t).$$ (8)

Tabei et al. (2002) suggested a factorization, which can factor their second-order $k$ -space operator into parts associated with each spatial direction. This factorization can also be applied to $kx$ -space operator by replacing constant velocity with variable velocity. The factored operators are called first-order $kx$ -space operators:

$$\begin{array}{l} \displaystyle \frac{\partial p(\mathbf{x},t)... ...rt\Delta t/2)\mathbf{F}[p(\mathbf{x},t)]\right]. \\ \end{array}$$ (9)

The spatial frequency components $k_x$ and $k_z$ are defined so that $k^2=k_x^2+k_z^2$ . Application of the exponential coefficient in equation 9 requires the corresponding wavefield to be evaluated on grid points staggered by distance of $\Delta x/2$ along the positive or negative $x$ direction and $\Delta z/2$ along the positive or negative $z$ direction. The spatial staggering in equation 9 is implicitly incorporated into the spatial derivative by the shift property of the Fourier transform. Using operators in equation 9 within equation 1 enables a new construction of spectral method. The first-order coupled equations for acoustic-wave extrapolation in variable velocity and density media with staggered spatial and temporal grids are therefore:

$$\begin{array}{l} \displaystyle \frac{u_x(\mathbf{x}_1,t^+) - ... ...rac{\partial u_z(\mathbf{x}_2,t^+)}{\partial^- z}), \end{array}$$ (10)

where $\mathbf{x}_1\equiv(x+\Delta x/2, z)$ , $\mathbf{x}_2\equiv(x, z+\Delta z/2)$ , $t^+ \equiv t+\Delta t/2$ , $t^- \equiv t-\Delta t/2$ .

The partial derivative operators in equation 10 are defined by equation 9. Note that the ordering of $\partial /\partial ^+x$ and $\partial/\partial^-x$ is arbitrary depending on the configuration of staggered grid. However, these operators should be used in pairs, such that the spatial shifting cancel out over any temporal interval length $\Delta t$ . The spatial and temporal staggered grids used in equation 10 are analogous to staggered scheme employed in previous finite-difference methods(Virieux, 1986; Madariaga, 1976; Virieux, 1984), which are know to increase accuracy and stability by halving spatial and time interval without increasing the number of computational points. Equation 9 can be solved with a localized Fourier transform (Wards et al., 2008). However, this kind of solution has a high computational cost. Song et al. (2012) proposed to apply FFD method to calculate the first-order $kx$ -space operators in equation 9, which can handle the variable velocity and density accurately and efficiently. Another possible ways to speed up the computation is to represent the extrapolation operator with a lowrank matrix.

Lowrank seismic wave extrapolation on a staggered grid