paper.tex

国外免费地震资料处理软件包

%\def\Xactiveplot#1#2#3#4{}
%\def\Xactivesideplot#1#2#3#4{}

\title{The helical coordinate}
\author{Jon Claerbout}
\maketitle
\label{paper:hlx}

For many years it has been true that
our most powerful signal-analysis techniques
are in {\em  one}-dimensional space,
while our most important applications are in {\em  multi}-dimensional space.
The helical coordinate system makes a giant step
towards overcoming this difficulty.

\par
Many geophysical map estimation problems appear to be multidimensional,
but actually they are not.
To see the tip of the iceberg, consider this example:
On a two-dimensional cartesian mesh, the function
$
\begin{array}{|r|r|r|r|}
        \hline
        0 & 0 & 0 &0 \\
        \hline
        0 & 1 & 1 &0 \\
        \hline
        0 & 1 & 1 &0 \\
        \hline
        0 & 0 & 0 &0 \\
        \hline
\end{array}
$
\par\noindent
has the autocorrelation
$
\begin{array}{|r|r|r|} \hline
        1 & 2 & 1\\
        \hline
        2 & 4 & 2\\
        \hline
        1 & 2 & 1
        \\ \hline
\end{array}
$.

\par\noindent
Likewise, on a one-dimensional cartesian mesh,
\par\noindent
the function $
\mathcal{B} =
\begin{array}{|r|r|r|r|r|r|r|r|r|r|r|r|r|} \hline
 1&1&0&0& \cdots& 0& 1&1
        \\ \hline
\end{array}
$
\par\noindent
has the autocorrelation
$ \mathcal{R} =
\begin{array}{|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|} \hline
 1&2&1&0&\cdots&0&2&4&2&0&\cdots&1&2&1
        \\ \hline
\end{array}
$.


%\par\noindent
%the function $ \bold b = ( 1,1,0,0, \cdots, 0, 1,1) $
%\par\noindent
%has the autocorrelation
%$ \bold r = ( 1,2,1,0,\cdots,0,2,4,2,0,\cdots,1,2,1)$.

\par\noindent
Observe the numbers in the one-dimensional world are identical
with the numbers in the two-dimensional world.
This correspondence is no accident.


\section{FILTERING ON A HELIX}
\inputdir{helicon}
Figure \ref{fig:diamond} shows some two-dimensional shapes
that are convolved together.
The left panel shows an impulse response function,
the center shows some impulses,
and the right shows the superposition of responses.
\plot{diamond}{width=6in,height=2.0in}{
  Two-dimensional convolution
  as performed in one dimension
  by module
  \texttt{helicon} %\vpageref{lst:helicon}
}

\par
A surprising, indeed amazing, fact is that
Figure \ref{fig:diamond} was not computed with a two-dimensional
convolution program.
It was computed with a one-dimensional computer program.
It could have been done with anybody's one-dimensional convolution program,
either in the time domain or in the fourier domain.
This magical trick is done with the helical coordinate system.

\par
A basic idea of filtering, be it in one dimension, two dimensions, or more,
is that you have some filter coefficients and some sampled data;
you pass the filter over the data; 
at each location you find an output by crossmultiplying
the filter coefficients times the underlying data and summing the terms.

\par
The helical coordinate system is much simpler than you might imagine.
Ordinarily, a plane of data is thought of as a collection of columns,
side by side.
Instead, imagine the columns stored end-to-end,
and then coiled around a cylinder.
This is the helix.
Fortran programmers will recognize that fortran's way
of storing 2-D arrays in one-dimensional memory
is exactly what we need for this helical mapping.
Seismologists sometimes use the word ``supertrace''
to describe a collection of seismograms stored ``end-to-end''.

\inputdir{Math}

Figure \ref{fig:sergey-helix} shows a helical mesh for 2-D data on a cylinder.
Darkened squares depict a 2-D filter
shaped like the Laplacian operator $\partial_{xx}+\partial_{yy}$.
The input data, the filter, and the output data are all on helical meshes
all of which could be unrolled into linear strips.
A compact 2-D filter like a Laplacian,
on a helix is a sparse 1-D filter with long empty gaps.

%\activeplot{sergey-helix}{width=6.0in,bb=210 315 630 545}{} {
%\activeplot{sergey-helix}{width=6.0in,bb=186 134 657 380}{} {
\plot{sergey-helix}{width=6.0in,bb=210 155 630 390}{
  Filtering on a helix.
  The same filter coefficients overlay the same data values
  if the 2-D coils are unwound into 1-D strips.
  (\emph{Mathematica} drawing by Sergey Fomel)
}

\par
Since the values output from filtering can be computed
in any order, we can slide the filter coil
over the data coil in any direction.
The order that you produce the outputs is irrelevant.
You could compute the results in parallel.
We could, however, slide the filter over the data in the screwing order
that a nut passes over a bolt.
The screw order is the same order that would be used
if we were to unwind the coils into one-dimensional strips
and convolve them across one another.
The same filter coefficients overlay the same data values
if the 2-D coils are unwound into 1-D strips.
The helix idea allows us to obtain the same convolution output
in either of two ways, a one-dimensional way, or a two-dimensional way.
I used the one-dimensional way to compute the obviously two-dimensional
result in Figure \ref{fig:diamond}.



\subsection{Review of 1-D recursive filters}
Convolution is the operation we do on polynomial coefficients
when we multiply polynomials.
Deconvolution is likewise for polynomial division.
Often these ideas are described
as polynomials in the variable $Z$.
Take $X(Z)$ to denote the polynomial
whose coefficients are samples of input data,
and let $A(Z)$ likewise denote the filter.
The convention I adopt here is that the first coefficient
of the filter has the value +1, so the filter's polynomial
is $A(Z) = 1 + a_1Z + a_2Z^2 + \cdots$.
To see how to convolve, we now identify the coefficient
of $Z^k$ in the product $Y(Z)=A(Z)X(Z)$.
The usual case ($k$ larger than the number $N_a$ of filter coefficients) is
\begin{equation}
y_k \quad=\quad x_k + \sum_{i=1}^{N_a} a_i x_{k-i}
\label{eqn:convolution}
\end{equation}
Convolution computes $y_k$ from $x_k$ whereas deconvolution
(also called back substitution) does the reverse.
Rearranging (\ref{eqn:convolution}) we get
\begin{equation}
x_k \quad=\quad y_k - \sum_{i=1}^{N_a} a_i x_{k-i}
\label{eqn:deconvolution}
\end{equation}
where now we are finding the output $x_k$ from
its past outputs $x_{k-i}$ and from the present input $y_k$.
We see that the deconvolution process is essentially
the same as the convolution process,
except that the filter coefficients
are used with opposite polarity;
and they are applied to the past {\em  outputs}
instead of the past {\em  inputs}.
That is why deconvolution must be done sequentially
while convolution can be done in parallel.



\subsection{Multidimensional deconvolution breakthrough}
Deconvolution (polynomial division)
can undo convolution (polynomial multiplication).
A magical property of the helix is that we can consider
1-D convolution to be the same as 2-D convolution.
Hence is a second magical property:
We can use 1-D
 {\em  de}convolution to undo convolution,
whether that convolution was 1-D or 2-D.
Thus, we have discovered how to undo 2-D convolution.
We have discovered that 2-D deconvolution on a helix
is equivalent to        1-D deconvolution.
The helix enables us to do multidimensional deconvolution.
\par
Deconvolution is recursive filtering.
Recursive filter outputs cannot be computed in parallel,
but must be computed sequentially
as in one dimension, namely,
in the order that the nut
screws on the bolt.  
\par
Recursive filtering sometimes solves big problems with astonishing speed.
It can propagate energy rapidly for long distances.
Unfortunately, recursive filtering can also be unstable.
The most interesting case, near resonance, is also near instability.
There is a large literature and extensive technology
about recursive filtering in one dimension.
The helix allows us to apply that technology to two (and more) dimensions.
It is a huge technological breakthrough.
\par
In 3-D we simply append one plane after
another (like a 3-D fortran array).
It is easier to code than to explain or visualize
a spool or torus wrapped with string, etc.



\subsection{Examples of simple 2-D recursive filters}
Let us associate $x$- and $y$-derivatives with
a finite-difference stencil or template.
(For simplicity take $\Delta x=\Delta y=1$.)
\begin{equation}
 {\partial \over \partial x } \eq
 \begin{array}{|c|c|} \hline
   1    & -1
   \\ \hline
 \end{array}
 \label{eqn:partialx}
\end{equation}
\begin{equation}
  {\partial \over \partial y } \eq
  \begin{array}{|r|} \hline
    1  \\
    \hline
    -1
    \\ \hline
  \end{array}
  \label{eqn:partialy}
\end{equation}
Convolving a data plane with
the stencil (\ref{eqn:partialx})
forms the $x$-derivative of the plane.
Convolving a data plane with
the stencil (\ref{eqn:partialy})
forms the $y$-derivative of the plane.
On the other hand,
{\it deconvolving}
with (\ref{eqn:partialx}) integrates data along the $x$-axis for each $y$.
Likewise, deconvolving
with (\ref{eqn:partialy}) integrates data along the $y$-axis for each $x$.
Next we look at a fully two-dimensional operator
(like the cross derivative $\partial_{xy}$).

\inputdir{helicon}

\par
A nontrivial two-dimensional convolution stencil is
\par
\begin{equation}
        \begin{array}{|r|r|}   \hline
                0    & -1/4 \\
                \hline
                1    & -1/4 \\
                \hline
                -1/4 & -1/4
                \\ \hline
        \end{array}
\label{eqn:filterfour}
\end{equation}
We will convolve and deconvolve a data plane with this operator.
Although everything is shown on a plane,
the actual computations are done in one dimension
with equations
(\ref{eqn:convolution}) and
(\ref{eqn:deconvolution}).
Let us manufacture the simple data plane
shown on the left in Figure \ref{fig:wrap-four}.
Beginning with a zero-valued plane, we add
in a copy of the filter (\ref{eqn:filterfour})
near the top of the frame.
Nearby add another copy with opposite polarity.
Finally add some impulses near the bottom boundary.
The second frame in Figure \ref{fig:wrap-four} is the result
of deconvolution by the filter (\ref{eqn:filterfour})
using the one-dimensional equation (\ref{eqn:deconvolution}).
Notice that deconvolution
turns the filter itself into an impulse,
while it turns the impulses
into comet-like images.
The use of a helix is evident
by the comet images wrapping around the vertical axis.


%\activeplot{wrap}{width=6in,height=3in}{} { 
%\activeplot{wrap}{width=3in,height=1.5in}{} { 
\plot{wrap-four}{width=4in,height=2in} { 
  Illustration of 2-D deconvolution.
  Left is the input.
  Right is after deconvolution with
  the filter (\protect\ref{eqn:filterfour})
  as preformed by
  by module
  \texttt{polydiv} %\vpageref{lst:polydiv}
}

The filtering in Figure \ref{fig:wrap-four}
ran along a helix from left to right.
Figure \ref{fig:back-four}
shows a second filtering running from right to left.
Filtering in the reverse direction is the adjoint.
After deconvolving both ways, we have accomplished a symmetical smoothing.
The final frame undoes the smoothing to bring us exactly back
to where we started.
The smoothing was done with two passes of {\it deconvolution}
and it is undone by two passes of {\it convolution}.
No errors, no evidence remains of any of the boundaries
where we have wrapped and truncated.

%\activeplot{pdadj}{width=6in,height=6.0in}{} { 
\plot{back-four}{width=4in,height=2.0in}{ 
  Recursive filtering backwards (leftward on the space axis)
  is done by the {\em  adjoint} of 2-D deconvolution.
  Here we see that 2-D {\it deconvolution} compounded with its adjoint
  is exactly inverted by 2-D {\it convolution} and its adjoint.
}

%\begin{notforlecture}
\par
Chapter \ref{prc/paper:prc} explains the important practical role
to be played by a multidimensional operator for which
we know the exact inverse.  
Other than multidimensional Fourier transformation,
transforms based on polynomial multiplication and division
on a helix are the only known easily invertible linear operators.

%\par
%I found myself with a powerful convolution/deconvolution
%program, but I did not have a bag full of 2-D filters.
%Further, I knew that a great number of filters are unstable
%for deconvolution.
%I recalled a variant of Rouch$\acute{\rm e}$'s theorem
%that if the sum of the absolute values of the filter coefficients
%after the onset pulse is less than the pulse,