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

Lowrank Finite Differences

In a matrix notation, the lowrank decomposition problem takes the following form:

$$\mathbf{W} \approx \mathbf{W}_1 \cdot \mathbf{A} \cdot \mathbf{W}_2,$$ (9)

where $\mathbf{W}$ is the $N_x\times N_x$ matrix with entries $W(\mathbf{x},\mathbf{k})$, $\mathbf{W}_1$ is the submatrix of $\mathbf{W}$ that consists of the columns associated with ${\mathbf{k}_m}$, $\mathbf{W}_2$ is the submatrix that consists of the rows associated with ${\mathbf{x}_n}$, and $\mathbf{A} = \{a_{mn}\}$.

Note that $\mathbf{W}_2$ is a matrix related only to wavenumber $\mathbf{k}$. We propose to further decompose it as follows:

$$\mathbf{W_2} \approx \mathbf{C} \cdot \mathbf{B},$$ (10)

where we determine $\mathbf{B}$ to be an $L\times N_x$ matrix, and the entry, $\mathbf{B}(\mathbf{\xi},\mathbf{k})$, has the form of $\cos(\sum\limits_{j=1}^3{\xi^j k_j\Delta x_j})$, in which $\mathbf{\xi}$ is a 3-D integer vector, $\mathbf{\xi}=(\xi^1,\xi^2,\xi^3)$, $k_j$ is the jth component of wavenumber $\mathbf{k}$, $\Delta x_j$ is the space grid size in the jth direction, $j=1,2,3$, and $\mathbf{C}$ is the matrix product of $\mathbf{W}_2$ and the pseudo-inverse of $\mathbf{B}$. In practice, we apply a weighted-inversion to achieve the pseudo-inverse: putting a larger weight on the low-wavenumber part and a smaller weight on the high-wavenumber part to enhance the stability. Now we have a new decompostion for the mixed-domain matrix:

$$\mathbf{W} \approx \mathbf{G} \cdot \mathbf{B},$$ (11)

where $\mathbf{G}$ is an $N_x\times L$ matrix,

$$\mathbf{G} = \mathbf{W}_1 \cdot \mathbf{A} \cdot \mathbf{C},$$ (12)

and

$$p(\mathbf{x},t+\Delta t) + p(\mathbf{x},t-\Delta t) =$$

$$2 \int e^{-i \mathbf{x} \mathbf{k}} W(\mathbf{x},\mathbf{k}) \hat... ...hbf{k}} B(\mathbf{\xi}_m,\mathbf{k}) \hat{p}(\mathbf{k},t) d\mathbf{k} \right)$$

$$\approx \sum\limits_{m=1}^L G(\mathbf{x},m) \left(\int e^{-i \ma... ...mits_{j=1}^3{\xi_m^j k_j\Delta x_j}) \hat{p}(\mathbf{k},t) d\mathbf{k} \right)$$

$$\approx \sum\limits_{m=1}^L G(\mathbf{x},m) \left[\int e^{-i \mat... ...mits_{j=1}^3{\xi_m^j k_j\Delta x_j}})\hat{p}(\mathbf{k},t) d\mathbf{k} \right].$$ (13)

According to the shift property of FFTs, we finally obtain an expression in the space-domain

$$p(\mathbf{x},t+\Delta t) + p(\mathbf{x},t-\Delta t) = \su... ...m=1}^L G(\mathbf{x},m) [p(\mathbf{x}_L,t)+p(\mathbf{x}_R,t)],$$ (14)

where $\mathbf{x}_L\,=\,(x_1-\mathbf{\xi}_m^1\Delta x_1,x_2-\mathbf{\xi}_m^2\Delta x_2,x_3-\mathbf{\xi}_m^3\Delta x_3)$, and $\mathbf{x}_R\,=\,(x_1+\mathbf{\xi}_m^1\Delta x_1,x_2+\mathbf{\xi}_m^2\Delta x_2,x_3+\mathbf{\xi}_m^3\Delta x_3)$.

Equation 14 indicates a procedure of finite differences for wave extrapolation: the integer vector, $\mathbf{\mathbf{\xi}_m}\,=\,(\xi_m^1,\xi_m^2,\xi_m^3)$ provides the stencil information, and $G(\mathbf{x},m)$ stores the corresponding coefficients. We call this method lowrank finite differences (LFD) because the finite-difference coefficients are derived from a lowrank approximation of the mixed-domain propagator matrix. We expect the derived LFD scheme to accurately propagate seismic-wave components within a wide range of wavenumbers, which has advantages over conventional finite differences that focus mainly on small wavenumbers. In comparison with the Fourier-domain approach, the cost is reduced to $O(L\,N_x)$, where $L$, as the row size of matrix $\mathbf{B}$, is related to the order of the scheme. $L$ can be used to characterize the number of FD coefficients in the LFD scheme, shown in equation 14. Take the 1-D 10th order LFD as an example, there are 1 center point, 5 left points ($\mathbf{x}_L$) and 5 right ones ($\mathbf{x}_R$). So $\xi_m^1={0,1,2,3,4,5}$, and $\mathbf{\xi}_m=(\xi_m^1,0,0)$. Thanks to the symmetry of the scheme, coefficients of $\mathbf{x}_L$ and $\mathbf{x}_R$ are the same, as indicated by equation 14. As a result, one only needs 6 coefficients: $L=6$.

Mexact,Mlrerr,Mapperr,Mfd10err
Figure 1. (a) Wavefield extrapolation matrix for 1-D linearly increasing velocity model. Error of wavefield extrapolation matrix by:(b) lowrank approximation, (c) the 10th-order lowrank FD (d) the 10th-order FD.

We use a one-dimentional example shown in Figure 1 to demonstrate the accuracy of the proposed LFD method. The velocity linearly increases from 1000 to 2275 m/s. The rank is 3 ($M=N=3$) for lowrank decomposition for this model with 1 ms time step. The propagator matrix is shown in Figure 1a. Figure 1b-Figure 1d display the errors corresponding to different approximations. The error by the 10th-order lowrank finite differences (Figure 1c) appears significantly smaller than that of the 10th-order finite difference (Figure 1d). Figure 2 displays the middle column of the error matrix. Note that the error of the LFD is significantly closer to zero than that of the FD method.

slicel Figure 2. Middle column of the error matrix. Solid line: the 10th-order LFD. Dash line: the 10th-order FD.

To analyze the accuracy, we let

$$p(\mathbf{x},t) = e^{i(\mathbf{k}\cdot\mathbf{x}-\omega t)},$$ (15)

by using the plane wave theory. Inserting 15 into equation 14 and also adopting the dispersion relation $\omega\,=\,\vert k\vert\,v$, defines the phase velocity of LFD ($v_{LFD}$) as follows:

$$v_{LFD}=\frac{1}{\vert k\vert\Delta t} \arccos(\sum\limit... ...\cos(\xi_m^2\,k_2\Delta x_2)+\cos(\xi_m^3\,k_3\Delta x_3))) ,$$ (16)

For 1-D 10th order LFD, $L=6$, $\mathbf{\xi}_m=(\xi_m^1,0,0)$ and $\xi_m^1={0,1,2,3,4,5}$. With equation 16, we can calculate phase-velocities ($v_{LFD}$) by 1-D 10th order LFD with different velocities ( $v=2500,3000,3500,4000$), and we use the ratio $\delta=v_{LFD}/v$ to describe the dispersion of FD methods. Figure 3a displays 1D dispersion curves by 1-D 10th order LFD, and Figure 3b shows those by conventional FD method. Note that compared with the conventional FD method, LFD is accurate in a wider range of wavenumbers (up to 70% of the Nyquist frequency).

app,fd10
Figure 3. Plot of 1-D dispersion curves for different velocities, $v=2500$ (red), $3000$ (pink), $3500$ (green), $4000$ (blue) $m/s$, $\Delta t=1$ ms, $\Delta x=10$ m by: (a) the 10th-order LFD (b) the 10th-order conventional FD.

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