paper.tex

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

7/30 & -44/15 & 57/5 & -44/15 & 7/30 \\
\hline
2/5 & -14/15 & -44/15 & -14/15 & 2/5 \\
\hline
-1/60 & 2/5 & 7/30 & 2/5 & -1/60 \\
\hline
\end{tabular}
= 1/60
\begin{tabular}{|c|c|c|c|c|}
\hline 
-1 & 24 & 14 & 24 & -1 \\
\hline
24 & -56 & -176 & -56 & 24 \\
\hline
14 & -176 & 684 & -176 & 14 \\
\hline
24 & -56 & -176 & -56 & 24 \\
\hline
-1 & 24 & 14 & 24 & -1 \\
\hline
\end{tabular}
\end{center}
The Laplacian representation with the same order of accuracy has the
coefficients
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline 
-1/360 & 2/45 & 0 & 2/45 & -1/360 \\
\hline
2/45 & -14/45 & -4/5 & -14/45 & 2/45 \\
\hline
0 & -4/5 & 41/10 & -4/5 & 0 \\
\hline
2/45 & -14/45 & -4/5 & -14/45 & 2/45 \\
\hline
-1/360 & 2/45 & 0 & 2/45 & -1/360 \\
\hline
\end{tabular} =
1/360 
\begin{tabular}{|c|c|c|c|c|}
\hline 
-1 & 16 & 0 & 16 & -1 \\
\hline
16 & -112 & -288 & -112 & 16 \\
\hline
0 & -288 & 1476 & -288 & 0 \\
\hline
16 & -112 & -288 & -112 & 16 \\
\hline
-1 & 16 & 0 & 16 & -1 \\
\hline
\end{tabular}
\end{center}
For the sake of simplicity, we assumed equal spacing in the $x$ and $y$
direction. The coefficients can be easily adjusted for anisotropic
spacing. Figures~\ref{fig:specc} and~\ref{fig:specp} show the spectra
of the finite-difference representations of operator~(\ref{eqn:bit})
for different values of the tension parameter. The
finite-difference spectra appear to be fairly isotropic (independent of
angle in polar coordinates).  They match the exact expressions at small
frequencies.

\plot{specc}{width=4.5in}{Spectra of the finite-difference
  splines-in-tension schemes for different values of the tension
  parameter (contour plots).}

\plot{specp}{width=4.5in}{Spectra of the
  finite-difference splines-in-tension schemes for different values of
  the tension parameter (cross-section plots). The dashed lines show
  the exact spectra for continuous operators.}
\par
Regarding the finite-difference operators as two-dimensional
auto-correlations and applying the Wilson-Burg method of spectral
factorization, we obtain two-dimensional minimum-phase filters
suitable for inverse filtering.  The exact filters contain many
coefficients, which rapidly decrease in magnitude at a distance from
the first coefficient. For reasons of efficiency, it is advisable to
restrict the shape of the filter so that it contains only the
significant coefficients. Keeping all the coefficients that are $1000$
times smaller in magnitude than the leading coefficient creates a
53-point filter for $\lambda=0$ and a 35-point filter for $\lambda=1$,
with intermediate filter lengths for intermediate values of $\lambda$.
Keeping only the coefficients that are $200$ times smaller that the
leading coefficient, we obtain 25- and 16-point filters for
respectively $\lambda=0$ and $\lambda=1$.  The restricted filters do
not factor the autocorrelation exactly but provide an effective
approximation of the exact factors. As outputs of the Wilson-Burg
spectral factorization process, they obey the minimum-phase condition.

\inputdir{gtens}

\plot{splin}{width=6in,height=6in}{Inverse filtering with the tension
  filters. The left plots show the inputs composed of filters and
  spikes. Inverse filtering turns filters into impulses and turns
  spikes into inverse filter responses (middle plots). Adjoint
  filtering creates smooth isotropic shapes (right plots). The tension
  parameter takes on the values 0.3, 0.7, and 1 (from top to bottom).
  The case of zero tension corresponds to Figure~\ref{fig:thin42}.}
\par
Figure~\ref{fig:splin} shows the two-dimensional filters for different
values of $\lambda$ and illustrates inverse recursive filtering, which
is the essence of the helix method \cite[]{GEO63-05-15321541}. The case of
$\lambda=1$ leads to the filter known as \emph{helix derivative}
\cite[]{iee}.  The filter values are spread mostly in two columns. The
other boundary case ($\lambda=0$) leads to a three-column filter,
which serves as the minimum-phase version of the Laplacian. This
filter is similar to the one shown in Figure~\ref{fig:thin42}.  As
expected from the theory, the inverse impulse response of this filter
is noticeably smoother and wider than the inverse response of the
helix derivative.  Filters corresponding to intermediate values of
$\lambda$ exhibit intermediate properties.  Theoretically, the inverse
impulse response of the filter corresponds to the Green's function of
equation~(\ref{eqn:bit}). The theoretical Green's function for the
case of $\lambda=1$ is
\begin{equation}
  \label{eqn:g1}
  G = \frac{1}{2\pi}\ln{r}\;,
\end{equation}
where $r$ is the distance from the impulse:
$r=\sqrt{\left(x-x_k\right)^2 + \left(y-y_k\right)}$. In the case of
$\lambda=0$, the Green's function is smoother at the origin:
\begin{equation}
  \label{eqn:g2}
  G = \frac{1}{8\pi}r^2\ln{r}\;.
\end{equation}
The theoretical Green's function expression for an arbitrary value of
$\lambda$ is unknown\footnote{\cite{mitas} derive an analytical Green's
  function for a different model of tension splines using special functions.}, 
but we can assume that its smoothness lies
between the two boundary conditions.
\par
In the next subsection, we illustrate an application of helical inverse
filtering to a two-dimensional interpolation problem.

\subsection{Regularization example}

We chose an environmental dataset \cite[]{iee} for a simple illustration
of smooth data regularization. The data were collected on a bottom
sounding survey of the Sea of Galilee in Israel \cite[]{zvi}.  The data
contain a number of noisy, erroneous and inconsistent measurements,
which present a challenge for the traditional estimation methods
\cite[]{galilee}.
%Addressing this challenge completely goes beyond the scope of this
%paper.
\par
Figure~\ref{fig:mesh} shows the data after a nearest-neighbor binning
to a regular grid. The data were then passed to an interpolation
program to fill the empty bins. The results (for different values of
$\lambda$) are shown in Figures~\ref{fig:gal} and~\ref{fig:cross}.
Interpolation with the minimum-phase Laplacian ($\lambda=0$) creates a
relatively smooth interpolation surface but plants artificial
``hills'' around the edge of the sea. This effect is caused by large
gradient changes and is similar to the sidelobe effect in the
one-dimensional example (Figure~\ref{fig:int}).  It is clearly seen in
the cross-section plots in Figure~\ref{fig:cross}.  The abrupt
gradient change is a typical case of a shelf break. It is caused by a
combination of sedimentation and active rifting. Interpolation with
the helix derivative ($\lambda=1$) is free from the sidelobe
artifacts, but it also produces an undesirable non-smooth behavior in
the middle part of the image. As in the one-dimensional example,
intermediate tension allows us to achieve a compromise: smooth
interpolation in the middle and constrained behavior at the sides of
the sea bottom.

\sideplot{mesh}{width=3in,height=4in}{The Sea of Galilee dataset after
  a nearest-neighbor binning. The binned data is used as an input for
  the missing data interpolation program.}

\plot{gal}{width=5.5in,height=7.33in}{The Sea of Galilee dataset after
  missing data interpolation with helical preconditioning. Different
  plots correspond to different values of the tension parameter. An
  east-west derivative filter was applied to illuminate the surface.}

\plot{cross}{width=6in,height=6in}{ Cross-sections of the Sea of
  Galilee dataset after missing-data interpolation with helical
  preconditioning.  Different plots correspond to different values of
  the tension parameter.}


\section{Conclusions}

\begin{comment}
We observe a significant (order-of-magnitude) speed-up in the
optimization convergence when preconditioning interpolation problems
with inverse recursive filtering. Since inverse filtering takes almost
the same time as forward convolution, this speed-up translates
straightforwardly into computational time savings.

The savings are hardly noticeable for simple test problems, but they
can have a direct impact on the mere feasibility of iterative
least-squares inversion for large-scale (seismic-exploration-size)
problems.
\end{comment}

The Wilson-Burg spectral factorization method, presented in this paper, allows
one to construct stable recursive filters. The method appears to have
attractive computational properties and can be significantly more efficient
than alternative spectral factorization algorithms. It is particularly
suitable for the multidimensional case, where recursive filtering is enabled
by the helix transform.

We have illustrated an application of the Wilson-Burg method for efficient
smooth data regularization. A constrained approach to smooth data
regularization leads to splines in tension. The constraint is embedded in a
user-specified tension parameter. The two boundary values of tension
correspond to cubic and linear interpolation.  By applying the method of
spectral factorization on a helix, we have been able to define a family of
two-dimensional minimum-phase filters, which correspond to the spline
interpolation problem with different values of tension.  We have used these
filters for accelerating data-regularization problems with smooth surfaces by
recursive preconditioning. In general, they are applicable for preconditioning
acceleration in various estimation problems with smooth models.

\section{Acknowledgments}
This paper owes a great deal to John Burg. We would like to thank him
and Francis Muir for many useful and stimulating discussions. The
first author also thanks Jim Berryman for explaining the variational
derivation of the biharmonic and tension-spline equations. Ralf Ferber and Ali
\"{O}zbek provided helpful reviews.

The financial support for this work was provided by the sponsors of
the Stanford Exploration Project.

\bibliographystyle{seg}
\bibliography{SEG,burg}

%\plot{name}{width=6in,height=}{caption}
%\sideplot{name}{height=1.5in,width=}{caption}
%
%\begin{equation}
%\label{eqn:}
%\end{equation}
%
%\ref{fig:}
%(\ref{eqn:})

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% TeX-master: t
%%% End: