Lowrank seismic wave extrapolation on a staggered grid

Lowrank FD for first-order extrapolation operators

Approximation 12 can also be used to design accurate finite-difference schemes. Here we extend the lowrank finite-difference method (Song et al., 2013) to first-order $kx$ -space operators. Note that $W_x(\mathbf{x}_n, \mathbf{k})$ in equation 13 is a matrix related only to wavenumber $\mathbf{k}$ . It can be further decomposed as follows:

$$\displaystyle W_x(\mathbf{x}_n,\mathbf{k}) \approx \sum\limits_{l=1}^LC(\mathbf{x}_n,\mathbf{\xi}_l)B(\mathbf{\xi}_l,\mathbf{k}),$$ (14)

where $B$ is an $L \times N_x$ matrix. Specifically, we can define $B(\mathbf{\xi}_l,\mathbf{k})$ to take the form of $\sin(\sum_{j=1}^3\xi_l^jk_j\Delta_j)$ , in which $\xi_l^j$ is the $j$ -th component of a $3$ -D vector, $\mathbf{\xi}_l=(\xi_l^1,\xi_l^2,\xi_l^3)$ , $k_j$ is the $j$ -th component of wavenumber $\mathbf{k}$ , $\Delta_j$ is the space grid size in the $j$ -th direction, $j=1,2,3$ corresponds to $x,y,z$ direction in space. $C$ is the matrix product of $W_x$ and the pseudo-inverse of $B$ . If we define

$$\displaystyle G(\mathbf{x},l) = \sum\limits_{m=1}^M\sum\limits_{n=1}^NW_x(\mathbf{x},\mathbf{k}_m)a_{mn}C(\mathbf{x}_n,\mathbf{\xi}_l),$$ (15)

then equation 13 can be described as

$$\begin{array}{l} \displaystyle \frac{\partial p(\mathbf{x},t)... ...\Delta_j})\mathbf{F}\left[p(\mathbf{x},t)\right]] . \end{array}$$ (16)

For a staggered grid, in which pressure is defined on main grid points and partial velocity on half-grid points, we can choose $\mathbf{\xi}_l$ as $\mathbf{\xi}_l^1 = (2l^1-1)/2, \: \mathbf{\xi}_l^2 = l^2, \: \mathbf{\xi}_l^3 = l^3$ to calculate partial derivative in $x$ -direction, where $l^1$ , $l^2$ , $l^3 = 1,2,\cdots,L$ , and $L$ is the length of the stencil. According to the shift property of FFTs, we can finally obtain the following expression in the space domain:

$$\displaystyle \frac{\partial p(\mathbf{x},t)}{\partial^+x} \appro... ...c{1}{2}\sum\limits_{l=1}^LG(\mathbf{x},l)[P(\mathbf{x}_R,t)-P(\mathbf{x}_L,t)],$$ (17)

where $\mathbf{x}_R=(x_1+l^1\Delta_1,x_2+l^2\Delta_2,x_3+l^3\Delta_3)$ and $\mathbf{x}_L=(x_1-(l^1-1)\Delta_1,x_2-l^2\Delta_2,x_3-l^3\Delta_3)$ .

Equation 17 corresponds to the procedure of finite-difference scheme for calculating the $kx$ -space operator. The vector $\mathbf{\xi}_l=(\xi_l^1,\xi_l^2,\xi_l^3)$ provides stencil information, and $G(\mathbf{x},l)$ stores the corresponding coefficients. A similar derivation can be applied to the remaining partial derivative operators in equation 9. We call this staggered grid lowrank finite-differences (SGLFD) method.

While the SGL method (equation 13) is proposed by applying lowrank approximation to the $kx-$ space method on a staggered grid (equation 9), SGLFD (equation 17) is a further approximation of SGL. Theoretically, the SGLFD method using a longer stencil reaches higher accuracy. It is hard to derive stability condition for SGL. However, applying Von Neumann stability analysis, we can easily obtain a sufficient condition of stability for SGLFD as

$$\left\vert\Delta tv_{max}\sum\limits_{l=1}^{L}G(x,l)sin(\sum\limits_{j=1}^3\xi_l^jk_j\Delta_j) \right\vert \leq 1,$$ (18)

where $v_{max}$ is the maximum value of velocity. Once we obtain the finite difference coefficient $G(\mathbf{x},l)$ for certain velocity $v(\mathbf{x})$ and the predefined parameters $\Delta_j$ and $\Delta t$ , we can use condition 18 to estimate the stability of the SGLFD scheme in equation 17.

Next, we use the plane wave theory to evaluate numerical dispersion for SGLFD method. Inserting the plane wave solution,

$$\begin{array}{l} \displaystyle p(\mathbf{x},t)=p_0e^{i\mathbf... ...\mathbf{u}_0e^{i\mathbf{k}\cdot\mathbf{x}-\omega t} \end{array}$$ (19)

into equation 17, and also adopting the dispersion relation $\omega=\vert\mathbf{k}\vert v$ , the relative error of phase velocity is defined as

$$\displaystyle \varepsilon=\frac{v_{LFD}}{v}-1=\frac{1}{\omega\Del... ...its_{l=1}^LG(\mathbf{x},l)\sum\limits_{j=1}^3\sin(\xi_l^jk_j\Delta_j)\right)-1.$$ (20)

The relative error $\varepsilon$ describes the numerical dispersion of SGLFD method. If $\varepsilon$ equals 0 , there is no dispersion. If $\varepsilon$ is far from 0, a large dispersion will occur. Here we define the order of SGLFD as that of conventional FD, which has the same stencil length ($L$ ). Next we compare the conventional SGFD method with the SGLFD method by the dispersion curves for different orders, time intervals and velocities.

Figure 1 shows the variation of $\varepsilon$ with frequency for different order. This figure demonstrates that dispersion decrease with the increase of the order for both SGFD and SGLFD method. Note that for SGFD method increase of order decreases the magnitude of the dispersion error without increasing the area where $\varepsilon$ nearly equals 0. Compared with the SGFD method, the SGLFD method is high accurate in a wider range of wavenumber. Figure 2 shows the variation of the $\varepsilon$ with frequency for difference time interval. From this figure, we can see that the dispersion becomes stronger when the SGFD method uses larger time interval. Moreover, if a large time interval is used, like $\Delta t=2.5ms$ in this example, the SGFD method will be unstable. However, for the SGLFD method, its dispersion mainly depends on the frequency. Compared with the SGFD method, the SGLFD method keeps high accuracy for different time intervals (up to $70\%$ of the Nyquist frequency). Figure 3 illustrates the effect of velocity on dispersion. Note that for the SGFD method, its dispersion curves change greatly with the variation of velocity. Compared with the SGFD method, the SGLFD method is more stable and accurate in a wider range of frequency(up to $70\%$ of the Nyquist frequency). In the previous examples, we used least squares to fit to nearly the $67\%$ of the Nyquist frequency.

Mfd,Mlr
Figure 1. Plot of $1$ -D dispersion curves of (a) the conventional SGFD method and (b) the SGLFD method for different orders, $4th-$ order(red, $2L=4$ ), $6th-$ order(pink, $2L=6$ ), $8th-$ order(green, $2L=8$ ), $16th-$ order(blue, $2L=16$ ), time interval $\Delta t=1ms$ , space interval $\Delta _x=10m$ , velocity $v=3000m/s$ .

Mfdt,Mlrt
Figure 2. Plot of $1$ -D dispersion curves of (a) the conventional SGFD method and (b) the SGLFD method for different time interval, $\Delta t=1ms$ (red), $\Delta t=1.5ms$ (pink), $\Delta t=2ms$ (green), $\Delta t=2.5ms$ (blue), $\Delta x=10m$ , $v=3000m/s$ , $2L=16$ .

Mfdv,Mlrv
Figure 3. Plot of $1$ -D dispersion curves of (a) the conventional SGFD method and (b) the SGLFD method for different velocity, $v=2500m/s$ (red), $v=3500m/s$ (pink), $v=4000m/s$ (green), $\Delta t=4500m/s$ (blue), $\Delta t=1ms$ , $\Delta x=10m$ , $2L=16$ .

Lowrank seismic wave extrapolation on a staggered grid