THE HELIX AND FINITE DIFFERENCES

Discretize the (x,y)-plane to an $N\times N$ array and pack the array into a vector of N² components. Likewise pack the Laplacian operator $\partial_{xx}+\partial_{yy}$into a matrix. For N=4 that matrix is shown in (5).

 

$$\begin{array} {\vert cccc\vert cccc\vert cccc\vert cccc\vert... ...ot & \cdot & 1 & -4 \\ &&&& &&&& &&&& &&& \\ \hline\end{array}$$ (5)

The two-dimensional matrix of coefficients for the Laplacian operator is shown in (5), where, on a cartesian space, h=0, and in the helix geometry, h=1. (A similar partitioned matrix arises from packing a cylindrical surface into a $4\times4$ array.) Notice that the partitioning becomes transparent for the helix, h=1.

With the partitioning thus invisible, the matrix simply represents one-dimensional convolution and we have an alternative analytical approach, Fourier Transform. We often need to solve sets of simultaneous equations with a matrix similar to (5). A costly method is to factor the matrix into upper and lower triangular form that can be ``backsolved'' which in this case amounts to recursive filtering.

The Fourier approach is similar but much faster. The (negative of the) Laplacian operator is regarded as an autocorrelation $(\cdots,-1,\cdots,-1,4,-1,\cdots,-1,\cdots)$.Using the Kolmogoroff spectral-factorization method, a ``minimum-phase'' wavelet is found which has this autocorrelation. This 2-D wavelet (rotated $90^\circ$ from (4)) is

 

$$\begin{array} {rrrrrrrrr} & & & & 1.000 & -.363 & -.024 & -... ...ts \\ \cdots &-.024 & -.048 & -.113 & -.312 & & & & \end{array}$$ (6)

The 2-D autocorrelation of (6) is the 2-D Laplace operator. The factorization (6) is nonunique in that causality can be defined along either x or y axes; thus the operator (6) can be rotated or reflected. Unfortunately, we see that the factored operator has an infinite number of nonzero terms, but fortunately they seem to be converging rapidly so truncation seems reasonable. We can use this feedback operator to solve Poisson's equation very rapidly. Using the seven coefficients shown, the cost is fourteen multiplications (because we need to run both ways) per mesh point. An example is shown in Figure 7.

 

lapfac
Figure 7 Deconvolution by a filter whose autocorrelation is the two-dimensional Laplacian operator. This amounts to solving the Poisson equation.

As a practical matter, this Poisson equation solver is a convenient isotropic smoothing filter.

Twenty years ago when migrations were two-dimensional, we developed powerful wave-equation methods for finite-difference downward-continuation of wave fields. These methods were frustrated in 3-D because we could not rapidly solve algebraic systems like $[1 - \alpha (\delta_{xx}+\delta_{yy})] P(z+\Delta z) = {\rm known}$.The helix should resuscitate these methods (although we will need to develop approximations for when $\alpha(x,y)$).