paper.tex

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

\texttt{mis2} \vpageref{lst:mis2}.


\section{MULTISCALE, SELF-SIMILAR FITTING}
\sx{multiscale fitting}
\sx{self-similar fitting}
\sx{goal ! multiscale self-similar}

Large objects often resemble small objects.
To express this idea we use {\it \bx{axis scaling}}
and we apply it to the basic theory
of prediction-error filter (PEF) fitting
and missing-data estimation.
\par
Equations (\ref{eqn:mspef}) and (\ref{eqn:msmis}) compute the same thing
by two different methods,
$ \bold r = \bold Y \bold a$ and
$ \bold r = \bold A \bold y$.
When it is viewed as fitting goals minimizing $||\bold r||$
and used along with suitable constraints,
(\ref{eqn:mspef}) leads to finding filters and \bx{spectra},
while
(\ref{eqn:msmis}) leads to finding \bx{missing data}.

\begin{equation}
\left[ 
\begin{array}{c}
  r_1 \\ 
  r_2 \\ 
  r_3 \\ 
  r_4 \\ 
  r_5 \\ 
  \hline
  r_6 \\
  r_7 \\
  r_8 \\
  r_9
  \end{array} \right] 
\eq
\left[ 
\begin{array}{ccc}
  y_2 & y_1 \\
  y_3 & y_2 \\
  y_4 & y_3 \\
  y_5 & y_4 \\
  y_6 & y_5 \\
  \hline
  y_3 & y_1 \\
  y_4 & y_2 \\
  y_5 & y_3 \\
  y_6 & y_4 
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  a_1 \\ 
  a_2 \end{array} \right]
\quad\quad\quad {\rm or}\quad\quad
\left[ 
\begin{array}{c}
  \bold r_1 \\ 
  \bold r_2
  \end{array} \right] 
\ =\ 
\left[ 
\begin{array}{c}
  \bold Y_1 \\
  \bold Y_2
  \end{array} \right] 
\;
\bold a
\label{eqn:mspef}
\end{equation}

\begin{equation}
\left[ 
\begin{array}{c}
  r_1 \\ 
  r_2 \\ 
  r_3 \\ 
  r_4 \\ 
  r_5 \\ 
  \hline
  r_6 \\ 
  r_7 \\ 
  r_8 \\ 
  r_9 
  \end{array} \right] 
\eq
\left[ 
\begin{array}{cccccc}
  a_2   & a_1   & \cdot & \cdot & \cdot  & \cdot \\
  \cdot & a_2   & a_1   & \cdot & \cdot  & \cdot \\
  \cdot & \cdot & a_2   & a_1   & \cdot  & \cdot \\
  \cdot & \cdot & \cdot & a_2   & a_1    & \cdot \\
  \cdot & \cdot & \cdot & \cdot & a_2    & a_1   \\
  \hline
  a_2   & \cdot & a_1   & \cdot & \cdot  & \cdot \\
  \cdot & a_2   & \cdot & a_1   & \cdot  & \cdot \\
  \cdot & \cdot & a_2   & \cdot & a_1    & \cdot \\
  \cdot & \cdot & \cdot & a_2   & \cdot  & a_1   
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6
  \end{array} \right]
\quad\quad {\rm or}\quad\quad
\left[ 
\begin{array}{c}
  \bold r_1 \\ 
  \bold r_2 
  \end{array} \right] 
\ =\ 
\left[ 
\begin{array}{c}
  \bold A_1 \\
  \bold A_2
\end{array}
\right]
\;
\bold y
\label{eqn:msmis}
\end{equation}

\par
A new concept embedded in (\ref{eqn:mspef}) and (\ref{eqn:msmis})
is that one filter can be applicable for different \bx{stretching}s
of the filter's time axis.
One wonders,
``Of all classes of filters,
what subset remains appropriate for stretchings of the axes?''

\subsection{Examples of scale-invariant filtering}
When we consider all functions with vanishing gradient,
we notice that the gradient vanishes whether it is represented as
$(1,-1)/\Delta x$ or as
$(1,0,-1)/2\Delta x$.
Likewise for the Laplacian, in one dimension or more.
Likewise for the wave equation, as long as there is no viscosity
and as long as the time axis and space axes are stretched
by the same amount.
The notion of ``dip filter'' seems to have no formal definition,
but the idea that the spectrum should depend mainly on slope
in Fourier space implies a filter that is scale-invariant.
I expect
the most fruitful applications
to be with \bx{dip filter}s.
\sx{filter ! dip}

\par
Resonance or \bx{viscosity} or damping easily spoils scale-invariance.
The resonant frequency of a filter shifts if we stretch the time axis.
The difference equations
\begin{eqnarray}
y_t - \alpha   y_{t-1} &=& 0  \label{eqn:sparsedecay1} \\
y_t - \alpha^2 y_{t-2} &=& 0  \label{eqn:sparsedecay}
\end{eqnarray}
both have the same solution $y_t = y_0 \alpha^{-t}$.
One difference equation has the filter
$(1,-\alpha)$,
while the other has the filter
$(1,0,-\alpha^2)$,
and $\alpha$ is not equal to $\alpha^2$.
Although these operators differ,
when $\alpha \approx 1$ they might provide the same general utility,
say as a roughening operator in a fitting goal.

\par
Another aspect to scale-invariance work is the presence
of ``parasitic'' solutions,
which exist but are not desired.
For example, another solution to $y_t -  y_{t-2}=0$
is the one that oscillates at the Nyquist frequency.

\par
(Viscosity does not necessarily introduce an inherent length
and thereby spoil scale-invariance.
The approximate frequency independence of sound absorption per wavelength
typical in real rocks is a consequence
of physical inhomogeneity at all scales.
See for example \bx{Kjartansson}'s \bx{constant Q} viscosity,
\sx{Q}
described in \bx{IEI}.
Kjartansson teaches that
the decaying solutions $t^{-\gamma}$ are scale-invariant.
There is no ``decay time'' for the function $t^{-\gamma}$.
Differential equations of finite order and
difference equations of finite order cannot produce
$t^{-\gamma}$ damping,
yet we know that such damping is important in observations.
It is easy to manufacture
$t^{-\gamma}$ damping
in Fourier space;
$\exp[(-i\omega)^{\gamma+1}]$ is used.
Presumably,
difference equations can make reasonable approximations
over a reasonable frequency range.)

\subsection{Scale-invariance introduces more fitting equations}
\inputdir{multiscale}
The fitting goals (\ref{eqn:mspef}) and (\ref{eqn:msmis}) have
about double the usual number of fitting equations.
Scale-invariance {\it introduces extra equations}.
If the range of scale-invariance is wide, there will be more equations.
Now we begin to see the big picture.
\begin{enumerate}
\item Refining a model mesh improves accuracy.
\item Refining a model mesh makes empty bins.
\item Empty bins spoil analysis.
\item If there are not too many empty bins we can find a PEF.
\item With a PEF we can fill the empty bins.
\item To get the PEF and to fill bins we need enough equations.
\item Scale-invariance introduces more equations.
\end{enumerate}
An example of these concepts is shown in Figure \ref{fig:mshole}.
\plot{mshole}{width=6.0in,height=3in}{
  Overcoming aliasing with multiscale fitting.
}

Additionally, when we have a PEF,
often we still cannot find missing data
because conjugate-direction iterations do not converge fast enough
(to fill large holes).
Multiscale convolutions should converge quicker
because they are like mesh-refinement, which is quick.
An example of these concepts is shown in Figure \ref{fig:msiter}.
\sideplot{msiter}{width=3.0in,height=3.9in}{
  Large holes are filled faster with
  multiscale operators.
}

%\begin{notforlecture}
\subsection{Coding the multiscale filter operator}
%Many applications could need slightly different equations, namely
%\begin{equation}
%\left[ 
%\begin{array}{c}
%  \bold 0 \\ 
%  \bold 0
%  \end{array} \right] 
%\quad = \quad
%\left[ 
%\begin{array}{cc}
%  w_1\bold I &    \bold 0\\ 
%     \bold 0 & w_2\bold I
%  \end{array} \right] 
%\
%\left[ 
%\begin{array}{c}
%  \bold Y_1 \\
%  \bold Y_2
%  \end{array} \right] 
%\;
%\bold a
%\ +\ 
%\left[ 
%\begin{array}{c}
%  \bold r_1 \\ 
%  \bold r_2
%  \end{array} \right] 
%\label{eqn:wtpef}
%\end{equation}
%and
%\begin{equation}
%\left[ 
%\begin{array}{c}
%  \bold 0 \\ 
%  \bold 0
%  \end{array} \right] 
%\quad = \quad
%\left[ 
%\begin{array}{cc}
%  w_1\bold I &    \bold 0\\ 
%     \bold 0 & w_2\bold I\\ 
%  \end{array} \right] 
%\
%\left[ 
%\begin{array}{c}
%  \bold A_1 \\
%  \bold A_2
%  \end{array} \right] 
%\;
%\bold y
%\ +\ 
%\left[ 
%\begin{array}{c}
%  \bold r_1 \\ 
%  \bold r_2
%  \end{array} \right] 
%\label{eqn:wtmis}
%\end{equation}
%where $w_1$ and $w_2$ are weighting scale factors.
%This added complication can be easily overcome when putting
%the scaling factor is put in the central loops in the program (which I did).
%It would be more efficient, however, to take the scale factor
%out of the inner loop and apply it to the residual directly.
%I didn't do this because it would clutter the program
%and an optimizing compiler
%might take care of it anyway.
\par
Equation~(\ref{eqn:mspef}) shows an example
where the first output signal is the ordinary one
and the second output signal used a filter interlaced with zeros.
We prepare subroutines that allow for more output signals,
each with its own filter interlace parameter
given in the table {\tt jump[ns]}.
Each entry in the jump table corresponds to
a particular scaling of a filter axis.
The number of output signals is {\tt ns} and the
number of zeros interlaced between filter points 
for the {\tt j}-th signal is {\tt jump[j]-1}.
%Also, there is the weighting scale factor for each scaling,
%{\tt wt(ns)}.
%\progdex{msicaf1}{multiscale conv 1-D}
\par
The multiscale helix filter is defined in module
\texttt{mshelix} \vpageref{lst:mshelix}, analogous to the
single-scale module \texttt{helix} \vpageref{lst:helix}.
A new function
\texttt{onescale}
extracts our usual helix filter
of one particular scale
from the multiscale filter.

\moddex{mshelix}{multiscale helix filter definition}{27}{45}{user/gee}

We create a multscale helix with module
\texttt{createmshelix} \vpageref{lst:createmshelix}.
An expanded scale helix filter is like an ordinary helix filter
except that the lags are scaled according to a \texttt{jump}.

\moddex{createmshelix}{create multiscale helix}{41}{81}{user/gee}

\par
First we examine code for estimating a prediction-error filter
that is applicable at many scales.
We simply invoke the usual filter operator
\texttt{hconest} \vpageref{lst:hconest}
for each scale.

\opdex{mshconest}{multiscale convolution, adjoint is the filter}{48}{51}{user/gee}


\sx{filter ! multiscale prediction-error}

\par
The \bx{multiscale prediction-error filter} finding subroutine
is nearly identical to the usual subroutine
\texttt{find\_pef()} \vpageref{lst:pef}.
(That routine cleverly ignores missing data while estimating a PEF.)
To easily extend {\tt pef} to multiscale filters
we replace its call to the ordinary helix
filter module
\texttt{hconest} \vpageref{lst:hconest}
by a call to \texttt{mshconest}.
\moddex{mspef}{multiscale PEF}{26}{49}{user/gee}
The purpose of \texttt{pack(dd,.true.)}
is to produce the one-dimensional array expected by
our solver routines.


\par
%Likewise, I coded up
Similar code applies to
the operator in
(\ref{eqn:msmis})
which is needed for missing data problems.
This is like
\texttt{mshconest} \vpageref{lst:mshconest}
except the adjoint is not the filter but the input.

\opdex{mshelicon}{multiscale convolution, adjoint is the input}{27}{31}{user/gee}
The multiscale missing-data module \texttt{msmis2}
is just like the usual missing-data module 
\texttt{mis2} \vpageref{lst:mis2}
except that
the filtering is done with the multiscale filter
\texttt{mshelicon}. % \vpageref{lst:mshelicon}.

\moddex{msmis2}{multiscale missing data}{29}{53}{user/gee}

%\putbib[MISC]

%\section{ONE WAY EQUATIONS}
%An interesting question is how to construct a one-way equation
%like the downward-continuation equation
%from two-way equations given us by physics.
%Even the existence of one-way equations is a fundamental theoretical question.
%Although I have nothing to say about existence,
%if they exist, here are some ideas how to construct them.
%These ideas might be useful for exotic wavetypes.
%I include them here mainly to draw together ideas of filters
%and partial differential equations.
%\par
%Take the differencing molecule of any wave-like physical process
%like the scalar wave equation.
%This will need to be a physical process that damps out as it goes to infinity.
%Suppose the thickness of this molecule in some direction is $m$ points.
%Place on the plane random impulses separated by $m$ or more points.
%Solve the missing data problem to fill in the values
%between the random points.
%This should be the solution to the given difference equation
%in the source free regions of space between pulses.
%Given this solution plane, solve for the 2-D prediction-error filter (PEF).
%The PEF should be a one-way equation for the original physical problem.
%We could dispense with the random numbers
%and solve two stages of least squares on the wave molecule itself.
%Perhaps we can find a filter whose autocorrelation is the Laplace operator.
%The filter's phase deserves investigation.

\section{References}
\reference{\bx{Canales}, L.L., 1984,
        Random noise reduction:
        54th Ann. Internat. Mtg.,
        Soc. Explor. Geophys.,
        Expanded Abstracts, 525-527.}

\reference{\bx{Rothman}, D., 1985,
        Nonlinear inversion, statistical mechanics,
        and residual statics estimation: Geophysics, {\bf 50}, 2784-2798
        }

\reference{\bx{Spitz}, S., 1991,
        Seismic trace interpolation in the F-X domain:
        Geophysics, {\bf 56}, 785-794.
        }
\newpage

%\end{notforlecture}