WHY THE RECURSIVE MULTIDIMENSIONAL PEF?

Geophysicists estimate multidimensional functions, especially the earth's impedance and velocity. We are often required to produce such functions to fill maps and 3-D volumes while our measurements provide us with only part of the information we need. Thus we are left with unavoidable problems such as ``null space'' and/or ``empty bins''. We should fill this void with ``something that looks like everything else.'' We should produce a gentle transition from the known to the unknown, preserving the statistical properties (spectrum) while we do so. This is also called ``hiding the data acquisition footprint.'' In one dimension we have the frequency spectrum and in higher dimensions we have also the dip spectrum and the velocity spectrum. These ideas are gathered in my book PVI. They have not yet reached industrial-scale practice, perhaps because of the slow rate of convergence of optimization methods such as conjugate gradients (CG) when applied to massive problems.

When I tried optimization methods (such as CG) to achieve ``minimum wiggliness'' I found I was solving differential equations, not efficiently by using integration (recursion), but inefficiently by applying differential operators millions of times. This suggests that unadorned CG is far from optimum and that huge gains should be possible. This raised my interest in using recursive methods. But how can we use recursion in multidimensional space?

Signal analysis and filter theory provide us with a wealth of knowledge of one-dimensional transformations, but the basic ideas, poles and zeros, appear to fail us in two dimensions, because mathematics teaches us little about polynomials in two variables. In one dimension, it is easy to write a convolution program and its inverse, and we have adequate theory to assure stability of the inverse. Furthermore, we can easily find such operators and their inverses for most spectra of interest. In two dimensions, I found I could not name a single roughening operator for which I knew a reasonably efficient inverse operator. For example, it is easy to speak of the inverse of combinations of one-dimensional derivatives, but hopeless in two dimensions where the gradient operator has no obvious inverse operator and where finding the inverse of the Laplacian operator amounts to solving Laplace's equation, a significant chore.

It is easy to code a two-dimensional convolution program. It amounts manipulating the coefficients of a two-dimensional polynomial. This suggests that we prepare a program for polynomial division by a two-dimensional polynomial. It is easy, easy only so long as we limit ourselves to 2-D filters with the ``1'' on the corner. A 2-D PEF, however, has its ``1'' along a side. I was unable to write a polynomial-division for such a filter until I fell upon the idea of helical boundary conditions. A few short hours after I came upon the helical-boundary trick, I understood that a longstanding obstacle, the stability of 2-D deconvolution, would be overcome too. Thus my hopes grew to their present size, hopes for bringing many theoretical estimation techniques up to industrial-scale projects.

Several pilot projects are well underway. The heart of the matter is using multidimensional recursive filters as preconditioners to solve geophysical data fitting problems far more rapidly. Early results are very encouraging. I plan to publish them with my coworkers, Sergey Fomel and Bob Clapp.