paper.tex

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

which is the in-plane information in Figure~\ref{fig:sneps}.

\par
To view genuinely 3-D information we must see a movie,
or we must compress the 3-D to 2-D.
There are only a small number of ways to compress 3-D to 2-D.
One is to select planes from the volume.
One is to sum the volume over one of its axes,
and the other is a compromise,
a filtering over the axis we wish to abandon before subsampling on it.
That filtering is a local smoothing.
If the local smoothing has motion
(out of plane dip) of various velocities (various dips),
then the desired process of smoothing the out of plane direction
is what we did in the in-plane direction
in Figure~\ref{fig:sneps}.
But Figure~\ref{fig:sneps} amounts to more than that.
It amounts to a kind of simultaneous smoothing in the {\it two}
most coherent directions
whereas in 3-D your eye can smooth in only {\it one} direction
when you turn your head along with the motion.

\par\noindent
\boxit{
        If the purpose of data processing
        is to collapse 3-D data volumes to 2-D
        where they are comprehensible to the human eye,
        then perhaps data-slope adaptive,
        low-pass filtering in the out-of-plane direction
        is the best process we can invent.
        }

\par

\par
My purpose in filtering the earthquake stacks
is to form a guiding ``pilot trace''
to the analysis of the traces {\it within} the bin.
Within each bin,
each trace needs small time shifts and perhaps a small temporal filter
to best compensate it to .\ .\ . to what? to the pilot trace,
which in these figures was simply a stack of traces in the bin.
Now that we have filtered in the range direction, however,
the next stack can be made with a better quality pilot.

\par
%Subroutine \texttt{losn2()} \vpageref{lst:losn2}
%does the ``patching,''
%invokes the signal partitioning operator
%\texttt{signoi} \vpageref{lst:signoi}
%on each patch,
%applies a tent-like weighting function,
%and reassembles the patches compensating for their overlap.
%The tents are easier to make here than they were in
%local prediction error 
%(\texttt{lopef2()} \vpageref{lst:lopef2})
%because they cover the space of filter {\it inputs} instead of {\it outputs.}
%(In early versions of this code,
%I used internal convolution
%instead of transient convolution,
%and it left me confused about the proper shape of the output data space.)
%\progdex{losn2}{local signoi}

%\begin{notforlecture}
%\subsection{Theory for noise along with missing data}
%Data $\bold d$ space can be decomposed into
%known plus missing parts,
%$\bold d = \bold k + \bold m$.
%We partition
%an identity operator $\bold I$ on the data space
%into parts that separate the known from the missing data
%$\bold I = \bold K + \bold M$.
%Thus data space can be written as
%\begin{equation}
%\bold d \eq \bold K \bold d + \bold M \bold m
%\end{equation}
%where
%$\bold K\bold d$ is zero-padded known data and
%all the components of $\bold m$ are freely adjustable.
%
%\par
%Data space $\bold d$
%can also be decomposed into
%signal plus noise,
%$\bold d = \bold s + \bold n$.
%Thus
%\begin{equation}
%\bold s + \bold n \eq \bold K \bold d + \bold M \bold m
%\label{eqn:snkm}
%\end{equation}
%
%\par
%Writing goals for signal and noise and then eliminating the noise
%by the constraint equation (\ref{eqn:snkm}) gives
%\begin{eqnarray}
%\bold 0&\approx&\bold N\bold n\eq \bold N(\bold K\bold d+\bold M\bold m-\bold s)\\
%\bold 0&\approx&\bold S\bold s\eq \bold S                               \bold s
%\end{eqnarray}
%Putting this in matrix form we have
%the operator needed in computation
%\begin{equation}
%\left[ 
%   \begin{array}{c}
%           \bold 0   \\
%           \bold 0 
%           \end{array}
%   \right] 
%%
%\quad\approx\quad
%%
%\left[ 
%   \begin{array}{cc}
%        -\bold N   & \bold N \bold M  \\
%         \bold S   & \bold 0
%          \end{array}
%   \right] \ 
%\left[ 
%        \begin{array}{c}
%                  \bold s \\
%                  \bold m
%          \end{array}
%\right] 
%\ + \ 
%\left[ 
%   \begin{array}{c}
%           \bold N \bold K \bold d   \\
%           \bold 0 
%           \end{array}
%   \right] 
%\end{equation}
%I have not had time to prepare an example.
%
%\end{notforlecture}



\section{SPACE-VARIABLE DECONVOLUTION}
%\begin{notforlecture}
\sx{time-variable deconvolution}
\sx{deconvolution ! time-variable}
\sx{filter ! time-variable}
Filters sometimes change with time and space.
We sometimes observe signals whose spectrum changes with position.
A filter that changes with position is called nonstationary.
We need an extension of our usual convolution operator
\texttt{hconest} \vpageref{lst:hconest}.
Conceptually,
little needs to be changed besides changing 
\texttt{aa(ia)} to
\texttt{aa(ia,iy)}.
But there is a practical problem.
Fomel and I have made the decision
to clutter up the code somewhat
to save a great deal of memory.
This should be important to people interested in
solving multidimensional problems with big data sets.
\par
Normally, the number of filter coefficients is many fewer
than the number of data points,
but here we have very many more.
Indeed, there are {\tt na} times more.
%The memory requirement for a filter was formerly less than for data.
Variable filters require {\tt na} times more memory than the data itself.
To make the nonstationary helix code more practical,
we now require the filters to be constant in patches.
The data type for nonstationary filters
(which are constant in patches)
is introduced in module
	\texttt{nhelix}, which is a simple modification of module
	\texttt{helix} \vpageref{lst:helix}.
	\moddex{nhelix}{non-stationary convolution}{12}{17}{user/gee}
What is new is the integer valued vector \texttt{pch(nd)}
the size of the one-dimensional (helix) output data space.
Every filter output point is to be assigned to a patch.
All filters of a given patch number will be the same filter.
Nonstationary helixes are created with
	\texttt{createnhelix}, which is a simple modification of module
	\texttt{createhelix} \vpageref{lst:createhelix}.
	\moddex{createnhelix}{create non-stationary helix}{29}{68}{user/gee}
Notice that the user must define the \texttt{pch(product(nd))}
vector before creating a nonstationary helix.
For a simple 1-D time-variable filter,
presumably \texttt{pch} would be something like
$(1,1,2,2,3,3,\cdots)$.
For multidimensional patching we need to think a little more.

\par
Finally, we are ready for the convolution operator.
The operator \texttt{nhconest} \vpageref{lst:nhconest}
allows for a different filter in each patch.
\opdex{nhconest}{non-stationary convolution}{51}{68}{user/gee}
A filter output \texttt{y[iy]}
has its filter from the patch \texttt{ip=aa->pch[iy]}.
\begin{comment}
The line
\texttt{t=a(ip,:)}
extracts the filter for the \texttt{ip}th patch.
If you are confused (as I am) about the difference
between \texttt{aa} and \texttt{a},
maybe now is the time to have a look at beyond Loptran
to the Fortran version.\footnote{
	http://sepwww.stanford.edu/sep/prof/gee/Lib/
	}
\end{comment}
%\end{notforlecture}
\par
Because of the massive increase in the number of filter coefficients,
allowing these many filters
takes us from overdetermined to very undetermined.
We can estimate all these filter coefficients
by the usual deconvolution fitting goal (\ref{mda/eqn:pefregression})
\begin{equation}
\bold 0 \quad\approx\quad
\bold r \eq \bold Y \bold K \bold a +\bold r_0
\end{equation}
but we need to supplement it with some damping goals, say
\begin{equation}\begin{array}{lll}
\bold 0  &\approx&      \bold Y  \bold K \bold a  +\bold r_0
\\
\bold 0  &\approx&      \epsilon\ \bold R \bold a 
\label{eqn:leakydecon}
\end{array}
\end{equation}
where $\bold R$ is a roughening operator to be chosen.

\par
Experience with missing data in Chapter \ref{iin/paper:iin}
shows that when the roughening operator $\bold R$ is a differential operator,
the number of iterations can be large.
We can speed the calculation immensely by ``preconditioning''.
Define a new variable $\bold m$ by
$\bold a=\bold R^{-1}\bold m$
and insert it into (\ref{eqn:leakydecon}) to get
the equivalent preconditioned system of goals.
\begin{eqnarray}
\label{eqn:goodleak}
\bold 0   &\approx &   \bold Y  \bold K \bold R^{-1}\bold m  \\
\bold 0   &\approx &   \epsilon \ \bold m
\label{eqn:dumbdamp}
\end{eqnarray}

\par
The fitting (\ref{eqn:goodleak}) uses the operator $\bold Y\bold K\bold R^{-1}$.
For $\bold Y$ we can use subroutine
\texttt{nhconest()} \vpageref{lst:nhconest};
for the smoothing operator $\bold R^{-1}$ we can use nonstationary
polynomial division
with operator \texttt{npolydiv()}:
%\begin{notforlecture}
\opdex{npolydiv}{non-stationary polynomial division}{56}{90}{user/gee}

Now we have all the pieces we need.
As we previously estimated stationary filters with
the module \texttt{pef} \vpageref{lst:pef},
now we can estimate nonstationary PEFs with
the module \texttt{npef} \vpageref{lst:npef}.
The steps are hardly any different.
\moddex{npef}{non-stationary PEF}{29}{61}{user/gee}
Near the end
of module \texttt{npef}
is a filter \texttt{reshape}
from a 1-D array to a 2-D array.
\begin{comment}
If you find it troublesome that
\texttt{nhconest} \vpageref{lst:nhconest}
was using the filter
during the optimization as already multidimensional,
perhaps again,
it is time to examine the Fortran code.
The answer is that there has been a conversion
back and forth partially hidden by Loptran.
%\end{notforlecture}
\end{comment}

\inputdir{tvdecon}

\par
Figure~\ref{fig:tvdecon} shows a synthetic data example using these programs.
As we hope for deconvolution, events are compressed.
The compression is fairly good, even though each event has
a different spectrum.
What is especially pleasing is that satisfactory results
are obtained in truly small numbers of iterations (about three).
The example is for two free filter coefficients $(1,a_1,a_2)$
per output point.
The roughening operator $\bold R$ was taken to be $(1,-2,1)$
which was factored into
causal and anticausal finite difference.
\plot{tvdecon}{width=6.0in,height=3in}{
  Time variable deconvolution with
  two free filter coefficients and a gap of 6.
}
\par
I hope also to find a test case with field data,
but experience in seismology
is that spectral changes are slow,
which implies unexciting results.
Many interesting examples should exist
in two- and three-dimensional filtering, however,
because reflector dip is always changing
and that changes the {\it spatial} spectrum.
\par
In multidimensional space, the smoothing filter $\bold R^{-1}$
can be chosen with interesting directional properties.
Sergey, Bob, Sean and I have joked about this code being
the ``double helix'' program because there are two multidimensional
helixes in it, one the smoothing filter, the other the deconvolution filter.
Unlike the biological helixes, however, these two helixes
do not seem to form a symmetrical pair.

\begin{exer}
\item
Is 
\texttt{nhconest} \vpageref{lst:nhconest}
the inverse operator to
\texttt{npolydiv} \vpageref{lst:npolydiv}?
Do they commute?
\item
Sketch the matrix corresponding to operator
\texttt{nhconest} \vpageref{lst:nhconest}.  {\sc hints:}
Do not try to write all the matrix elements.
Instead draw short lines to indicate rows or columns.
As a ``warm up'' consider a simpler case where one filter
is used on the first half of the data and another filter
for the other half.
Then upgrade that solution from two to about ten filters.
\end{exer}


%\reference{
%       Shearer, P.M., 1991,
%       Imaging global body wave phases by stacking long-period seismograms.
%       J. Geophy. Res.,  v.~96, n.~B12,
%       pp. 20,535--20,324.
%       % November 10, 1991
%       }
%
%\reference{
%       Shearer, P.M., 1991,
%       Constraints on upper mantle discontinuities from observations
%       of long period reflected and converted phases.
%       J. Geophy. Res.,  v.~96, n.~B11,
%       pp. 18,147--18,182.
%       % October 10, 1991
%       }
%
%I need help from Diane here.
%%\putbib[MISC,SEP]

\bibliographystyle{sep}\bibliography{SEP,MISC}

%XXX putbib[SEP,MISC] busted