paper.tex

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

%then the filter is ``positive real,''
%hence ``minimum phase,'' hence stable in polynomial division.
%Ironically, it is the ``almost unstable'' filters
%that hold the greatest interest,
%because they are the ones that move energy long distances.
%To make long impulse-responses of deconvolution,
%I chose all the adjustable filter coefficients to be negative
%and to sum to $-1$,
%as does (\ref{eqn:filterfour}).
%For damping, I took the absolute value of the sum to be
%a little less than $1$.

\par
In seismology we often have occasion to steer summation along beams.
Such an impulse response is shown in Figure \ref{fig:dip}.
\plot{wrap-waves}{width=3in,height=1.5in}{ 
  A simple low-order 2-D filter whose inverse
  contains plane waves of two different dips.
  One of them is spatially aliased.
}

Of special interest are filters that destroy plane waves.  The inverse
of such a filter creates plane waves.  Such filters are like wave
equations.  A filter that creates two plane waves is illustrated in
figure \ref{fig:wrap-waves}.

\plot{dip}{width=6in,height=2.0in}{ 
  A simple low-order 2-D filter whose inverse times its inverse adjoint,
  is approximately a dipping seismic arrival.
}

\subsection{Coding multidimensional de/convolution}
Let us unroll the filter helix seen in Figure \ref{fig:sergey-helix}
and see what we have.
Start from the idea that a 2-D filter is generally made
from a cluster of values near one another in two dimensions
similar to the Laplacian operator in the figure.
We see that in the helical approach,
a 2-D filter is a 1-D filter containing some long intervals of zeros.
The intervals are about the length of a 1-D seismogram.

\par
Our program for 2-D convolution with a 1-D convolution program,
could convolve with the somewhat long 1-D strip,
but it is much more cost effective to ignore the many zeros,
which is what we do.
We do not multiply by the backside zeros, nor do we even store them in memory.
Whereas an ordinary convolution program would do time shifting
by a code line like {\tt iy=ix+lag},
Module
\texttt{helicon} \vpageref{lst:helicon}
ignores the many zero filter values on backside of the tube
by using the code {\tt iy=ix+lag[ia]}
where a counter {\tt ia} ranges over the nonzero filter coefficients.
Before operator {\tt helicon} is invoked,
we need to prepare two lists,
one list containing nonzero filter coefficients {\tt flt[ia]},
and the other list containing the corresponding lags {\tt lag[ia]}
measured to include multiple wraps around the helix.
For example, the 2-D Laplace operator
can be thought of as the 1-D filter
\begin{equation}
\label{eqn:2dlapfil}
\begin{array}{|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|} \hline
1&0&\cdots&0&1&-4&1&0& \cdots&0&1
\\ \hline
\end{array}
%
\quad
\xrightarrow{\rm helical \; boundaries}
\quad
\begin{array}{|r|r|r|}  \hline
& 1 & \\
\hline
1 & -4 & 1\\
\hline
& 1 & \\
\hline
\end{array}
\end{equation}
%
%
%\begin{equation}
%\label{eqn:2dlapfil}
%\left[ \begin{array}{ccc}
%& 1 & \\ 1 & -4 & 1\\ & 1 &
%\end{array} \right]
%\xrightarrow{\rm helical \; boundary \; conditions}
%(1,\,0, \; ... \; 0,\,1,\,-4,\,1,\,0, \; ... \; 0,\,1).
%\end{equation}
%
The first filter coefficient in equation (\ref{eqn:2dlapfil})
is $+1$ as implicit to module {\tt helicon}.
To apply the Laplacian on a $1000\times 1000$ mesh
requires the filter inputs:
\par\noindent
\footnotesize
\begin{verbatim}
                i    lag[i]   flt[i]
               ---   ------   -----
                0      999      1
                1     1000     -4
                2     1001      1
                3     2000      1
\end{verbatim}
\normalsize

\par
Here we choose to use
``declaration of a type'',
a modern computer language feature that is absent from Fortran 77.
Fortran 77 has the built in complex arithmetic type.
In module \texttt{helix}
we define a type \texttt{filter}, actually, a helix filter.
After making this definition, it will be used by many programs.
The helix filter consists of three vectors,
a real valued vector of filter coefficients,
an integer valued vector of filter lags,
and an optional vector
that has logical values ``\texttt{true}''~for
output locations that will not be computed
(either because of boundary conditions or because of missing inputs).
The filter vectors are the size of the nonzero filter coefficents
(excluding the leading 1.) while the logical vector is long
and relates to the data size.
The \texttt{helix} module allocates and frees memory for a helix filter.
By default, the logical vector is not allocated but
is set to \texttt{null}
with the \texttt{nullify} operator and ignored.
\moddex{helix}{definition for helix-type filters}{27}{32}{user/gee}
\par
For those of you with no C  experience,
the ``\verb#->#'' appearing in the helix module denotes a pointer.
Fortran 77 has no pointers (or everything is a pointer).
%The C, C++, and Java languages use ``\texttt{.}'' to denote pointers.
%C and C++ also have a second type of pointer denoted by ``\texttt{->}''.
The behavior of pointers is somewhat different in each language.
Never-the-less, the idea is simple.
In module \texttt{helicon} \vpageref{lst:helicon}
you see the expression
\verb#aa->flt[ia]#.
It refers to the filter named \texttt{aa}.
Any filter defined by the \texttt{helix} module
contains three vectors, one of which is named \texttt{flt}.
The second component of the \texttt{flt} vector
in the \texttt{aa} filter
is referred to as
\verb#aa->flt[1]# which
in the example above refers to the value 4.0
in the center of the laplacian operator.
For data sets like above with 1000 points on the 1-axis,
this value 4.0 occurs after 1000 lags,
thus \verb#aa->lag[1]=1000#.

\par
Our first convolution operator
\texttt{tcai1}
\vpageref{lst:tcai1}
was limited to one dimension and a particular choice of end conditions.
With the helix and C pointers,
the operator
\texttt{helicon} \vpageref{lst:helicon}
is a {\it multidimensional} filter
with considerable flexibility (because of the \texttt{mis} vector)
to work around boundaries and missing data.
\opdex{helicon}{helical convolution}{36}{55}{user/gee}
The code fragment
\verb#aa->lag[ia]#
corresponds to 
\texttt{b-1}
in \texttt{tcai1} \vpageref{lst:tcai1}.


\par
Operator {\tt helicon} did the convolution job for Figure \ref{fig:diamond}.
As with
\texttt{tcai1} \vpageref{lst:tcai1}
the adjoint of filtering is filtering backwards
which means unscrewing the helix.

\par
The companion to convolution is deconvolution.
The module \texttt{polydiv} \vpageref{lst:polydiv}
is essentially the same as
\texttt{polydiv1} \vpageref{lst:polydiv1},
but here it was coded using
our new \texttt{filter} type in
module \texttt{helix} \vpageref{lst:polydiv1}
which will simplify our many future uses of
convolution and deconvolution.
Although convolution allows us to work around missing input values,
deconvolution does not
(any input affects all subsequent outputs),
so \texttt{polydiv} never references \verb#aa->mis[ia]#.
\opdex{polydiv}{helical deconvolution}{39}{70}{user/gee}

\begin{exer}
\item
Observe the matrix (\ref{ajt/eqn:contran1})
which corresponds to
subroutine
\texttt{tcai1} \vpageref{lst:tcai1}.
What is the matrix corresponding to
\texttt{helicon} \vpageref{lst:helicon}?
\end{exer}

\subsection{Causality in two-dimensions}
In one dimension, most filters of interest have a short memory.
Significant filter coefficients are concentrated shortly after $t=0$.
The favorite example in Physics is the
damped harmonic oscillator,
all of which is packed into a two-lag filter
(second order differential equation).
The complete story is rich in mathematics and in concepts,
but to sum up, filters fall into two categories according to the
numerical values of their coefficients.
There are filters for which equations
(\ref{eqn:convolution}) and
(\ref{eqn:deconvolution})
work as desired and expected.
These filters are called ``minimum phase''.
There are also filters for which
(\ref{eqn:deconvolution}) is a disaster numerically,
the feedback process diverging to infinity.

\par
Divergent cases correspond to physical processes
that require boundary conditions.
Equation (\ref{eqn:deconvolution})
only allows for initial conditions.
I oversimplify by trying to collapse an entire book (FGDP)
into a few sentences by saying here that
for any fixed spectrum
there exist many filters.
Of these, only one has
stable polynomial division.
That filter has its energy compacted
as soon as possible after the ``1.0'' at zero lag.

%\end{notforlecture}
\par
Now let us turn to two dimensions.
Filters of interest will correspond to energy
concentrated near the end of a helix.
Let us examine the end of a helix.
At the very end, as in the 1-D case,
is a coefficient with the numerical value 1.0.
Keeping only coefficients within two mesh points
in any direction from the 1.0,
we copy the coefficients from
near the end of the helix to a cartesian mesh like this:
%A causal filter in one dimension has
%a curious shape on the two-dimensional helix.
%This strange shape plays a central role in coming chapters.
%I use the convention that the zero-lag
%response of the 1-D filter has the value ``1''.
%In one dimension,
%the causal filter has zeros before the ``1'' and various values after it.
%Let us consider a compact 2-D filter,
%one whose nonzero filter coefficients lie within
%a short distance in two dimensions (two lags) from the ``1''.
%Consider this compact filter to be compressed up against
%one end of a helix coil.
%We pluck it from the coil and
%display it as a two-dimensional array
\begin{equation}
\label{eqn:2dpef}
\begin{array}{ccccc}
        \begin{array}{ccc}
                h  & c &   0 \\
                p  & d &   0 \\
                q  & e &  \bold 1    \\
                s  & f &   a \\
                u  & g &   b
        \end{array}
        &\quad=\quad&
        \begin{array}{ccc}
                h  & c &  \cdot   \\
                p  & d &  \cdot   \\
                q  & e &  \cdot    \\
                s  & f &   a \\
                u  & g &   b
        \end{array}
         &\quad +\quad&
        \begin{array}{ccc}
                \cdot &\cdot &   0 \\
                \cdot &\cdot &   0 \\
                \cdot &\cdot &  \bold 1    \\
                \cdot &\cdot &  \cdot    \\
                \cdot &\cdot &  \cdot  
        \end{array}
  \\
  \\
  \\
        {\rm 2-D \ filter}
           &=&
        {\rm variable}
          &\quad +\quad&
        {\rm constrained}
\end{array}
\end{equation}
where $a,b,c,...,u$ are adjustable coefficients.

\par
Which side of the little rectangular patch of coefficients
we choose to place the 1.0 is rather arbitrary.
The important matter is that as a matter of principle,
the 1.0 is expected to lie along one side of the little patch.
It is rarely (if ever) found at a corner of the patch.
It is important that beyond the 1.0 (in whatever direction that may be)
the filter coefficients must be zero because in one dimension,
these coefficients lie before zero lag.
Our foundations,
the basic convolution-deconvolution pair
(\ref{eqn:convolution}) and
(\ref{eqn:deconvolution})
are applicable only to filters with all coefficients {\it after} zero lag.


\par
Time-series analysis is rich with concepts that
the helix now allows us to apply to many dimensions.
First is the notion of an impulse function.
Observe that an impulse function on the 2-D surface
of the helical cylinder maps to an impulse function
on the 1-D line of the unwound coil.
An autocorrelation function that is an impulse
corresponds both to a white (constant) spectrum in 1-D and
to a white (constant) spectrum in 2-D.
Next we look at a particularly important autocorrelation function
and see how 2-D is the same as 1-D.



\section{FINITE DIFFERENCES ON A HELIX}
The function
\begin{equation}
\label{eqn:2dlapin1d}
\mathcal{R}
\eq
\begin{array}{|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|} \hline
-1&0&\cdots&0&-1&4&-1&0& \cdots&0&-1
\\ \hline
\end{array}
\end{equation}
is an autocorrelation function.
It is symmetrical about the ``4'' and its Fourier transform
is positive for all frequencies.
Digging out our old textbooks\footnote{
        PVI or FGDP, for example,
        explain spectral factorization.
        More concisely in PVI, more richly in FGDP.
        }
we discover
how to compute a causal wavelet with this autocorrelation.
I used the ``Kolmogoroff spectral-factorization method''
to find this wavelet $\mathcal{H}$:
\begin{equation}
{\mathcal{H}}
\eq
\begin{array}{|r|r|r|r|r|r|r|r|r|r|r|r|r|r|r|}
\hline