The Helmholtz equation

Starting from the full wave equation in three-dimensions:

$$-\nabla^2 p + \frac{1}{v^2} \frac{\partial^2 p}{\partial t^2}=0$$ (5)

we can Fourier transform the time axis, and look for $(\omega, {\bf x}$)solutions of the form:

$$p(\omega,{\bf x}) = q ({\bf x}) \; e^{-i \omega t}$$ (6)

For a single frequency, the wave equation therefore reduces to the Helmholtz (time-independent diffusion) equation  

$$\left(-\nabla^2 - \alpha^2\right) q({\bf x})=0$$ (7)

where $\alpha=\omega/v$.

We aim to factor this equation on a helix, as with the Poisson equation above. However, before we can, we need to ensure that it is a `level-phase' function Claerbout (1998), that is to say the spectrum of the operator does not touch the negative real axis on the complex plane. The spectrum of the Helmholtz operator can be obtained by taking the Fourier transform of equation (7).

$$S({\bf k}) = \vert{\bf k}\vert^2 - \alpha^2$$ (8)

$S({\bf k})$ clearly becomes negative real for small values of $\vert {\bf k} \vert$; so as it stands, this equation is not factorable. Fortunately, however, replacing $\alpha$ by $\alpha'= \alpha-\epsilon / i$, where $\epsilon$ is a small positive number, successfully stabilizes the spectrum, by pushing the function off the negative real axis. The physical effect of $\epsilon$ is to provide damping as the wave propagates, differentiating between the forward and backward extrapolation directions.

Before factorization, equation (7) should therefore be rewritten to include the stabilization term  

$$\left[ -\nabla^2 + (-i \alpha + \epsilon)^2 \right] q({\bf x}) = \left( -\nabla^2 - \alpha'^2 \right) q({\bf x}) = 0$$ (9)

Following the helix solution to Poisson's equation above, a simple finite-difference approximation to the Laplacian, $-\nabla^2 \approx {\bf D}$, produces the matrix equation:

$$\left( {\bf D} - \alpha'^2 {\bf I} \right) {\bf q} = {\bf 0}$$ (10)

Alternatively, and more accurately, we can form a rational approximation to the Laplacian operator,  

$$-\nabla^2 \approx \frac{\bf D}{{\bf I} + \beta {\bf D}}$$ (11)

where $\beta \; (\approx 1/6)$ is Claerbout's 1985 adjustable `one-sixth' parameter, and ${\bf D}$ again represents convolution with a simple finite-difference filter, d.

Inserting equation (11) into equation (9) yields a matrix equation of similar form, but with increased accuracy at high spatial wavenumbers:

$$\begin{eqnarray} \left( \frac{\bf D}{{\bf I} + \beta {\bf D}} - \alpha'^2 {\bf I... ...lpha'^2 {\bf I} \right] {\bf q} = {\bf H} \; {\bf q} & = & {\bf 0}\end{eqnarray}$$ (12) (13)

The operator on the left-hand-side of equation (13) represents a three-dimensional convolution matrix, that can be mapped to an equivalent one-dimensional convolution by applying helical boundary conditions. Although the complex $\alpha'$ coefficients on the main diagonal cause the matrix not to be Hermitian, the spectrum of the matrix is of level-phase. Therefore, for constant $\alpha'$, it can be factored into causal and anti-causal (triangular) components with any spectral factorization algorithm that has been adapted for cross-spectra Claerbout (1998).  

$${\bf H} \; {\bf q} = {\bf U L} \; {\bf q} = {\bf 0}$$ (14)

The challenge of extrapolation is to find ${\bf q}$ that satisfies both the above equation and our initial conditions, ${\bf q}_{z=0}$. Starting from ${\bf q}_{z=0}$, we can invert ${\bf L}$ recursively to obtain a function that satisfies both the initial conditions, and

$${\bf L} \; {\bf q} = {\bf 0}.$$ (15)

Hence ${\bf q}$ will also satisfy equation (14).