The Helmholtz equation

Starting from the acoustic wave equation in three-dimensions,

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

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

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

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

$$\left(-\nabla^2 - \frac{\omega^2}{v^2}\right) \; q({\bf x})=0.$$ (19)

We can solve the Helmholtz equation on a regular grid by approximating the differential operator with a finite-difference stencil. A simple (all-zero) convolutional approximation to the Laplacian, $-\nabla^2 \approx {\bf D}$, produces the matrix equation:  

$$\left( {\bf D} - \frac{\omega^2}{v^2} {\bf I} \right) {\bf q} = {\bf 0}$$ (20)

Unfortunately the direct solution of equation () requires the inversion of a multi-dimensional convolutional matrix. As with Poisson's equation above, the application of helical boundary conditions allows me to factor the convolutional filter into a minimum-phase causal and anti-causal pair that can be inverted rapidly by polynomial division.

Factoring the Helmholtz operator is not as simple as factoring the Poisson operator, since its spectrum is not positive definite. The spectrum of the differential Helmholtz operator can be obtained by taking the spatial Fourier transform of equation (), to give

$$\left(\vert{\bf k}\vert^2 - \frac{\omega^2}{v^2}\right) \; \hat{q}({\bf k}) = S({\bf k}) \; \hat{q}({\bf k}) = 0,$$ (21)

where $\hat{q}({\bf k})$ represents the spatial Fourier transform of $q({\bf x})$, and $S({\bf k})=\vert{\bf k}\vert^2 - \frac{\omega^2}{v^2}$ is the Fourier representation of the Helmholtz operator. $S({\bf k})$ clearly becomes negative real for small values of $\vert {\bf k} \vert$; so as it stands, the Helmholtz operator does not represent an autocorrelation, and is not factorable. This problem also exists for discrete operators [e.g. equation ()].

Fortunately replacing $\omega$ by ${\hat \omega}= \omega + \epsilon i$, where $\epsilon$ is a small positive number, successfully stabilizes the spectrum, by pushing the function off the negative-real axis. Equation (), therefore becomes  

$$\left( {\bf D} - \frac{{\hat \omega}^2}{v^2} {\bf I} \right) {\bf q} = {\bf 0}.$$ (22)

The physical effect of $\epsilon$ is to provide damping as the wave propagates, differentiating between the forward and backward extrapolation directions. With the complex ${\hat \omega}$, the Helmholtz operator is symmetric but not Hermitian, and so represents a crosscorrelation rather than an autocorrelation.

Generalizing the concept of spectral factorization to cross-spectral factorization, Claerbout (1998c) showed that any function whose Fourier transform does not wrap around the origin in the complex plane can be factored into the crosscorrelation of two minimum-phase factors. Since for this class of function, the phase of the Fourier component at the positive Nyquist equals the phase at the negative Nyquist (with no bulk phase shift), he termed this class of function `level-phase'. Although the Helmholtz operator is not strictly an autocorrelation, it does satisfy the `level-phase' criterion, and so it can still be factored into a pair of minimum-phase factors.

Rather than considering a simple convolutional approximation to the Laplacian, we can form a rational approximation,  

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

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. It is safe to write equation () in this rational form since the numerator and denominator commute with each other.

Inserting equation () into equation () 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}} - \frac{{\hat \ome... ...bf D} - {\bf I} \right] {\bf q} = {\bf H} \; {\bf q} & = & {\bf 0}\end{eqnarray}$$ (24) (25)

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

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

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}.$$ (27)

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