Exploring three-dimensional implicit wavefield extrapolation with the helix transform

Implicit versus Explicit extrapolation

The difference between implicit and explicit extrapolation is best understood through an example. Following Claerbout (1985), let us consider, for instance, the diffusion (heat conduction) equation of the form

$${\frac{\partial T}{\partial t}} = {a (x)\,{\frac{\partial^2 T}{\partial x^2}}}\;.$$ (1)

Here $t$ denotes time, $x$ is the space coordinate, $T (x,t)$ is the temperature, and $a$ is the heat conductivity coefficient. Equation (1) forms a well-posed boundary-value problem if supplied with the initial condition

$$\left.T\right\vert _{t=0} = T_0 (x)$$ (2)

and the appropriate boundary conditions. Our task is to build a digital filter, which transforms a gridded temperature $T$ from one time level to another.

It helps to note that when the conductivity coefficient $a$ is constant and the space domain of the problem is infinite (or periodic) in $x$, the problem can be solved in the wavenumber domain. Indeed, after the Fourier transform over the variable $x$, equation (1) transforms to the ordinary differential equation

$${\frac{d \hat{T}}{d t}} = {- a k^2\, \hat{T}}\;,$$ (3)

which has the explicit analytical solution

$$\hat{T} (k,t) = \hat{T}_0 (k) e^{- a k^2 t}\;,$$ (4)

where $\hat{T}$ denotes the Fourier transform of $T$, and $k$ stands for the wavenumber. Therefore, the desired filter in the wavenumber domain has the form

$$H (k) = e^{- a k^2}\;,$$ (5)

where for simplicity the coefficient $a$ is normalized for the time step $\triangle t$ equal to $1$.

Returning now to the time-and-space domain, we can approach the filter construction problem by approximating the space-domain response of filter (5) in terms of the differential operators $\frac{\partial^2}{\partial x^2} = - k^2$, which can be approximated by finite differences. An explicit approach would amount to constructing a series expansion of the form

$$H_{\mbox{ex}} (k) \approx a_0 + a_1 k^2 + a_2 k^4 + \ldots\;,$$ (6)

and selecting the coefficients $a_j$ to approximate equation (5). For example, the three-term Taylor series expansion around the zero wavenumber yields

$$H_{\mbox{ex}} (k) = 1 - a\,{k^2} + {\frac{{{a }^2}\,{k^4}}{2}} \;.$$ (7)

The error of approximation (7) as a function of $k$ for two different values of $a$ is shown in the left plot of Figure 1.

exerror,imerror
Figure 1. Errors of second-order explicit (a) and implicit (b) approximations for the heat extrapolation.

An implicit approach also approximates the ideal filter (5), but with a rational approximation of the form

$$H_{\mbox{im}} (k) \approx \frac{b_0 + b_1 k^2 + b_2 k^4 + \ldots} {1 + c_1 k^2 + c_2 k^4 + \ldots} \;.$$ (8)

One way of selecting the coefficients $b_i$ and $c_i$ is to apply an appropriate Padé approximation (Baker and Graves-Morris, 1981). For example the $[2/2]$ Padé approximation is

$$H_{\mbox{im}} (k) = \frac{1 - \frac{a}{2}\,k^2}{1 + \frac{a}{2}\,k^2} \;.$$ (9)

This approximation corresponds to the famous Crank-Nicolson implicit method (Crank and Nicolson, 1947). The error of approximation (9) as a function of $k$ for different values of $a$ is shown in the right plot of Figure 1. Not only is it significantly smaller than the error of the same-order explicit approximation, but it also has a negative sign. It means that the high-frequency numerical noise gets suppressed rather than amplified. In practice, this property translates into a stable numerical extrapolation.

The second derivative operator $-k^2$ can be approximated in practice by a digital filter. The most commonly used filter has the $Z$-transform $D_2 (Z) = -Z^{-1} + 2 - Z$, and the Fourier transform

$$D_2 (k) = e^{-ik} - 2 + e^{-ik} = 2 (\cos{k} - 1) = -4 \sin^2{\frac{k}{2}}\;.$$ (10)

Formula (10) approximates $-k^2$ well only for small wavenumbers $k$. As shown in Appendix A, the implicit scheme allows the accuracy of the second-derivative filter to be significantly improved by a variation of the ``1/6-th trick'' (Claerbout, 1985). The final form of the implicit extrapolation filter is

$$H_{\mbox{im}} (k) = \frac{1 + \frac{a+\beta}{2}\,D_2 (k)}{1 - \frac{a-\beta}{2}\,D_2 (k)} \;,$$ (11)

where $\beta$ is a numerical constant, found in Appendix A.

heat
Figure 2. Heat extrapolation with explicit and implicit finite-different schemes. Explicit extrapolation appears stable for $a=2/3$ (left plot) and unstable for $a=4/3$ (middle plot). Implicit interpolation is stable even for larger values of $a$ (right plot).

A numerical 1-D example is shown in Figure 2. The initial temperature distribution is given by a step function. The discontinuity at the step gets smoothed with time by the heat diffusion. The left plot shows the result of an explicit extrapolation with $a=2/3$, which appears stable. The middle plot is an explicit extrapolation with $a=4/3$, which shows a terribly unstable behavior: the high-frequency numerical noise is amplified and dominates the solution. The right plot shows a stable (though not perfectly accurate) extrapolation with the implicit scheme for the larger value of $a=2$.

The difference in stability between explicit and implicit schemes is even more pronounced in the case of wave extrapolation. For example, let us consider the ideal depth extrapolation filter in the form of the phase-shift operator (Claerbout, 1985; Gazdag, 1978)

$$W (k) = e^{i \sqrt{a^2 - k^2}}\;,$$ (12)

where $a = \omega / v$, $\omega$ is the time frequency, and $v$ is the seismic velocity (which may vary spatially); we assume for simplicity that both the depth step $\triangle z$ and the space sampling $\triangle x$ are normalized to $1$. A simple implicit approximation to filter (12) is

$$W_{\mbox{im}} (k) = e^{i a}\, \frac{1 -4 a^2 + i a\,k^2}{1 - 4 a^2 - i a\,k^2} = e^{i \phi}\;,$$ (13)

where $\phi = a - 2 \arctan{\frac{a\,k^2}{4 a^2-1}}$. We can see that approximation (13) is again a pure phase shift operator, only with a slightly different phase. For that reason, the operator is unconditionally stable for all values of $a$: the total wave energy from one depth level to another is preserved. Operator (12) corresponds to the Crank-Nicolson scheme for the 45-degree one-way wave equation (Claerbout, 1985). Its phase error as a function of the dip angle $\theta = \arcsin{\frac{k}{a}}$ for different values of $a$ is shown in Figure 3.

phase Figure 3. The phase error of the implicit depth extrapolation with the Crank-Nicolson method.

The unconditional stability property is not achievable with the explicit approach, though it is possible to increase the stability of explicit operators by using relatively long filters (Hale, 1991b; Holberg, 1988).

Exploring three-dimensional implicit wavefield extrapolation with the helix transform