Fourier finite-difference wave propagation

Theory

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

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

where $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 can get 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)

Equation 2 has the following 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)-2p(\mathbf{x},t) = }$$

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

Equation 5 provides an elegant and efficient solution in the case of a constant-velocity medium with the aid of FFT. In the case of a variable-velocity medium, equation 5 can provide an approximation by replacing $v$ with $v(\mathbf{x})$ . However, FFT can no longer be applied directly for the inverse Fourier transform from the wavenumber domain back to the space domain. To overcome this problem, Etgen and Brandsberg-Dahl (2009) propose a velocity interpolation method. They present an implementation for isotropic, VTI (vertical transversely isotropic) and TTI (tilted transversely isotropic) media. In the isotropic case, two FFTs can be sufficient. For anisotropic media, more than one velocity parameter must be used. Therefore, it is necessary to perform velocity interpolation by combining different parameters and computing the corresponding forward and inverse FFTs for each of the velocity parameters, thus increasing the computational burden. Other FFT-based solutions include the optimized separable approximation or OSA (Du et al., 2010; Song, 2001; Liu et al., 2009; Zhang and Zhang, 2009) and the lowrank approximation (Fomel et al., 2010).

We propose an alternative approach. First, we adopt the following form of the right-hand side of equation 5 in the variable velocity case:

$${2\left[\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1\right] = }$$

$$2\left[\cos(v_0\vert\mathbf{k}\vert\Delta t)-1\right]\left[\frac{... ...t\mathbf{k}\vert\Delta t)-1}{\cos(v_0\vert\mathbf{k}\vert\Delta t)-1}\right]\;,$$ (6)

where $v_0$ is the reference velocity, such as the RMS (root-mean-square) velocity of the medium. After that, we apply the following approximation:

$$\frac{\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1}{\cos(v_0... ...mathbf{k}\vert\Delta t)-1} \approx a + 2\sum^3_{n=1}{b_n\cos(k_n\Delta x_n)}\;,$$ (7)

where coefficients $a$ and $b_n$ are defined using the Taylor expansion around $k=0$ ,

$${\begin{array}{*{20}c} \displaystyle a=\frac{v^2(\mathbf{x})}{v_0... ...thbf{x})(v^2(\mathbf{x})-v_0^2)}{12(\Delta x_n^2)v_0^2} \, \\ \end{array} } \;,$$ (8)

and $\Delta x_n$ is the sampling in the n-th direction. We only need to calculate these coefficients once. After completing the calculation, they can be used at each time step during the wave extrapolation process.

Equation 6 consists of two terms: the first term is independent of $\mathbf{x}$ and only depends on $\mathbf{k}$ . For this part, we use inverse FFT to return to the space domain from the wavenumber domain. For the remaining part, however, we can avoid phase shift in the wavenumber domain by implementing space shifts through finite differences (approximation 7) with coefficients provided by equation 8. This approach is analogous to the FFD method proposed by Ristow and Ruhl (1994) for one-way extrapolation in depth.

Figure 1 (a) shows approximations for $[\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1]$ by the 4th-order FD method (dash line) and pseudo-spectral method (dotted line). Figure 1 (b) shows approximations by the FFD method (2nd-order: dash line, 4th-order: dotted line). The solid lines stand for the exact values for function $[\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1]$ with true velocity: $v=4.0 km/s$ (bottom solid line) and $v_0=2.0 km/s$ (top solid line), which indicates a significant velocity contrast ($100\%$ difference). In this situation, all the approximations deviate from the exact solution as the wavenumber $\vert k\vert$ becomes large. However, the 4th-order FFD method approximates the exact solution with the most accuracy, as shown in the error plot (Figure 2). In order to enhance the stability, one can suppress the wavefield at high wavenumbers for both pseudo-spectral and the FFD method.

cos-side1
Figure 1. Different approximations for $\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1$ . Solid lines: exact solution ( $\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1$ ) for $v=4.0$ km/s and $v_0=2.0$ km/s. (a) Dash line: the 4th-order FD. Dotted line: pseudo-spectral method. (b) Dash line: the 2nd-order FFD method. Dotted line: the 4th-order FFD with $v_0$ as reference velocity. $\Delta t=0.001$ s. $\Delta x=0.005\,km$ .

diff1
Figure 2. Errors for different approximations for $\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1$ . Solid lines: exact solution ( $\cos(v(\mathbf{x})\vert\mathbf{k}\vert\Delta t)-1$ ) for $v=4.0$ km/s and $v_0=2.0$ km/s. (a) Dash line: the 4th-order FD. Dotted line: pseudo-spectral method. (b) Dash line: the 2nd-order FFD method. Dotted line: the 4th-order FFD with $v_0$ as reference velocity. $\Delta t=0.001$ s. $\Delta x=0.005\,km$ .

Assuming that $p(\mathbf{x},t-\Delta t)$ and $p(\mathbf{x},t)$ are already known, the time-marching algorithm can be specified as follows:

1. Transform $p(\mathbf{x},t)$ to $\hat{p}(\mathbf{k},t)$ by 3-D FFT;
2. Multiply $\hat{p}(\mathbf{k},t)$ by $2\left[\cos(v_0\vert\mathbf{k}\vert\Delta t)-1\right]$ to get $\hat{q}(\mathbf{k},t)$ ;
3. Transform $\hat{q}(\mathbf{k},t)$ to $q(\mathbf{x},t)$ by inverse FFT;
4. Apply finite differences to $q(\mathbf{x},t)$ with coefficients in equation 8 to get $q(\mathbf{x},t+\Delta t)$ . Namely,
   

   $${q^{i,j,k}(t+\Delta t) = a^{i,j,k}q^{i,j,k}(t)}$$

   $$+b_1^{i,j,k}(q^{i-1,j,k}(t)+q^{i+1,j,k}(t))$$

   $$+b_2^{i,j,k}(q^{i,j-1,k}(t)+q^{i,j+1,k}(t))$$

   $$+b_3^{i,j,k}(q^{i,j,k-1}(t)+q^{i,j,k+1}(t))\;,$$ (9)

   
   where $i$ is the grid index of $x_i$ direction;
5. $p(\mathbf{x},t+\Delta t) \leftarrow q(\mathbf{x},t+\Delta t) + 2p(\mathbf{x},t) - p(\mathbf{x}, t-\Delta t)$ .

Here $q(\mathbf{x},t)$ and $\hat{q}(\mathbf{k},t)$ are temporary functions.

The FFD approach is not limited to the isotropic case. In the case of transversally isotropic (TTI) media, the term $v(\mathbf{x})\,\vert\mathbf{k}\vert$ on the left-hand side of equation 7, can be replaced with the acoustic approximation (Fomel, 2004; Alkhalifah, 1998,2000),

$$f(\mathbf{v},\mathbf{\hat{k}},\eta)=\sqrt{\frac{1}{2}(v_1^2\,\hat... ...\,\hat{k}_2^2)^2-\frac{8\eta}{1+2\eta}v_1^2v_2^2\,\hat{k}_1^2\,\hat{k_2^2}}}\;,$$ (10)

where $v_1$ is the P-wave phase velocity in the symmetry plane, $v_2$ is the P-wave phase velocity in the direction normal to the symmetry plane, $\eta$ is the anisotropic elastic parameter (Alkhalifah and Tsvankin, 1995) related to Thomsen's elastic parameters $\epsilon$ and $\delta$ (Thomsen, 1986) by

$$\frac{1+2\delta}{1+2\epsilon}=\frac{1}{1+2\eta}\,;$$

and $\hat{k}_1$ and $\hat{k}_2$ stand for the wavenumbers evaluated in a rotated coordinate system aligned with the symmetry axis:

$$\begin{array}{*{20}c} \hat{k}_1=k_1\cos{\phi}+k_2\sin{\phi}\;... ...eta}-k_2\cos{\phi}\sin{\theta}+k_3\cos{\theta}\;,\\ \end{array}$$ (11)

where $\theta$ is the tilt angle measured with respect to vertical and $\phi$ is the angle between the projection of the symmetry axis in the horizontal plane and the original X-coordinate. The symmetry axis has the direction of $\left\{-\sin\theta\sin\phi,-\sin\theta\cos\phi,\cos\theta\right\}$ .

Using these definitions, we develop a finite-difference approximation analogous to equation 7 for FFD in TTI media. The details of the derivation are given in the appendix. For the 2D TTI case, the corresponding FFD algorithm is as following:

1. Transform $p(\mathbf{x},t)$ to $\hat{p}(\mathbf{k},t)$ by FFT;
2. Multiply $\hat{p}(\mathbf{k},t)$ by $$\frac{2[\cos(f(\mathbf{v_0},\mathbf{\hat{k}},\eta_0)\Delta t)-1]}{\vert\mathbf{\hat{k}}\vert^2}$$ to get $\hat{q}(\mathbf{k},t)$ ;
3. Transform $\hat{q}(\mathbf{k},t)$ to $q(\mathbf{x},t)$ by inverse FFT;
4. Apply finite differences to $q(\mathbf{x},t)$ with coefficients in Table 1 to get $q(\mathbf{x},t+\Delta t)$ . Namely,
   

   $${q^{i,j,t+\Delta t} = a^{i,j}q^{i,j,t}}$$

   $$+b_1^{i,j}(q^{i-1,j,t}+q^{i+1,j,t})$$

   $$+b_2^{i,j}(q^{i,j-1,t}+q^{i,j+1,t})$$

   $$+d_1^{i,j}(q^{i-2,j,t}+q^{i+2,j,t})$$

   $$+d_2^{i,j}(q^{i,j-2,t}+q^{i,j+2,t})$$

   $$+c^{i,j}(q^{i-1,j-1,t}+q^{i-1,j+1,t}+q^{i+1,j-1,t}+q^{i+1,j+1,t})\;.$$

   
   where $i$ is the grid index of $x_i$ direction;
5. $p(\mathbf{x},t+\Delta t) \leftarrow q(\mathbf{x},t+\Delta t) + 2p(\mathbf{x},t) - p(\mathbf{x}, t-\Delta t)$ .

Here, $q(\mathbf{x},t)$ and $\hat{q}(\mathbf{k},t)$ are temporary functions.

Fourier finite-difference wave propagation