paper.tex

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

surprisingly small number of iterations, the output closely resembles
the final output. The final output of interpolation with
preconditioning by recursive deconvolution is exactly
the same as that of the original method.

\plot{sall}{width=6in,height=4in}{Interpolation with
  preconditioning. Left plot: the top shows the input data; the
  middle, the result of interpolation; the bottom, the filter
  impulse response. The right plot shows the convergence process for
  the first four iterations.}

The next section extends the idea of preconditioning by inverse
recursive filtering to multiple dimensions.

\section{Multidimensional recursive filter preconditioning}

\inputdir{Math}

\cite{helix} proposed a \emph{helix} transform for mapping
multidimensional convolution operators to their one-dimensional
equivalents.  This transform proves the feasibility of
multidimensional deconvolution, an issue that has been in question for
more than 25 years \cite[]{Claerbout.blackwell.76}.  By mapping discrete
convolution operators to one-dimensional space, the inverse filtering
problem can be conveniently recast in terms of recursive filtering, a
well-known part of the digital filtering theory.

\plot{helix}{width=5in,bb=210 155 630 390}{Helix transform of
  two-dimensional filters to one dimension (a scheme).  The
  two-dimensional filter (plot \texttt{a}) is equivalent to the
  one-dimensional filter in (plot \texttt{d}), assuming that a shifted
  periodic condition is imposed on one of the axes (plots \texttt{b}
  and \texttt{c}.)}

The helix filtering idea is schematically illustrated in
Figure~\ref{fig:helix1}. The left plot (labeled ``\texttt{a}'' in the
figure) shows a two-dimensional digital filter overlayed on the
computational grid. A two-dimensional convolution computes its output
by sliding the filter over the plane. If we impose helical boundary
conditions on one of the axes, the filter will slide to the beginning
of the next trace after reaching the end of the previous one
(plot~``\texttt{b}''). As evident from plots~``\texttt{c}''
and~``\texttt{d}'', this is completely equivalent to one-dimensional
convolution with a long 1-D filter with internal gaps.  For
efficiency, the gaps are simply skipped in a helical convolution
algorithm. The gain is not in the convolution
itself, but in the ability to perform recursive inverse filtering
(deconvolution) in multiple dimensions. A multi-dimensional filter is
mapped to its 1-D analog by imposing helical boundary conditions on
the appropriate axes. After that, inverse filtering is applied
recursively in a one-dimensional manner.  Neglecting parallelization
and indexing issues, the cost of inverse filtering is equivalent to
the cost of convolution. It is proportional to the data size and to
the number of non-zero filter coefficients.

\inputdir{hlx}

\plot{waves}{width=4in,height=2in}{ Illustration of 2-D
  deconvolution with helix transform.  Left is the input: two spikes
  and two filters.  Right is the output of deconvolution.}

An example of two-dimensional recursive filtering is shown in
Figure~\ref{fig:waves}. The left plot contains two spikes and two
filter impulse responses with different polarity. After deconvolution
with the given filter, the filter responses turn into spikes, and the
initial spikes turn into long-tailed inverse impulse responses (right
plot in Figure~\ref{fig:waves}). Helical wrap-around, visible on the
top and bottom boundaries, indicates the direction of the helix.
\cite{gee} presents more examples and discusses all the issues of
multidimensional helical deconvolution in detail.

%As is known from the one-dimensional theory \cite[]{Claerbout.fgdp.76},
%a stable recursive filtering requires a minimum-phase filter, which
%can be constructed with a spectral factorization algorithm. 

\section{Multidimensional examples}

\inputdir{seab}

Our first multidimensional example is the SeaBeam dataset, a result of
water bottom measurements from a single day of acquisition. SeaBeam is
an apparatus for measuring water depth both directly under a ship and
somewhat off to the sides of the ship's track.  The dataset has been
used at the Stanford Exploration Project for benchmarking different
strategies of data interpolation.  The left plot in Figure
\ref{fig:seabdat} shows the original data. The right plot shows the
result of (unpreconditioned) missing data interpolation with the
Laplacian filter after 200 iterations.  The result is unsatisfactory,
because the Laplacian filter does not absorb the spatial frequency
distribution of the input dataset. We judge the quality of an
interpolation scheme by its ability to hide the footprints of the
acquisition geometry in the final result. The ship track from the
original acquisition pattern is clearly visible in the Laplacian
result, which is an indication of a poor interpolation
method.

\plot{seabdat}{width=6in,height=2.5in}{On the left, the
  SeaBeam data: the depth of the ocean under ship tracks; on the
  right, an interpolation with the Laplacian filter.}
\par
We can obtain a significantly better image (Figure \ref{fig:seabold})
by replacing the Laplacian filter with a two-dimensional
prediction-error filter estimated from the input data.  The
result in the left plot of Figure \ref{fig:seabold} was obtained after
200 conjugate-gradient iterations. If we stop after 20 iterations, the
output (the right plot in Figure \ref{fig:seabold}) shows only a small
deviation from the input data. Large areas of the image remain
unfilled. At each iteration, the interpolation process progresses only
to the length of the filter.

\plot{seabold}{width=6in,height=2.5in}{SeaBeam interpolation
  with the prediction-error filter. The left plot was taken after 200
  conjugate-gradient iterations; the right, after 20 iterations.}
\par
Inverting the PEF convolution with the help of the helix transform, we
can now apply the inverse filtering operator to precondition the
interpolation problem. As expected, the result after 200 iterations
(the left plot in Figure \ref{fig:seabnew}) is similar to the result
of the corresponding unpreconditioned (model-space) interpolation.
However, the output after just 20 iterations (the right plot in Figure
\ref{fig:seabnew}) is already fairly close to the solution.

\plot{seabnew}{width=6in,height=2.5in}{SeaBeam interpolation
  with the inverse prediction-error filter. The left plot was taken
  after 200 conjugate-gradient iterations; the right, after 20
  iterations.}

\inputdir{cube}

For a more practical test, we chose the North Sea seismic reflection
dataset, previously used for testing azimuth moveout and
common-azimuth migration \cite[]{GEO63-02-05740588,SEG-1997-1375}.
Figure~\ref{fig:cmp-win} shows the highly irregular midpoint geometry
for a selected in-line and cross-line offset bin in the data. The data
irregularity is also evident in the bin fold map, shown in
Figure~\ref{fig:fold-win}.  The goal of data regularization is to
create a regular data cube at the specified bins from the irregular
input data, which have been preprocessed by normal moveout without
stacking.

\sideplot{cmp-win}{width=3.5in}{ Midpoint distribution for a 50 by 50 m
  offset bin in the 3-D North Sea dataset.}
 
\sideplot{fold-win}{width=3.5in}{Map of the fold distribution for the
  3-D data test.}

The data cube after normalized binning is shown in
Figure~\ref{fig:bin-win}.  Binning works reasonably well in the areas of
large fold but fails to fill the zero fold gaps and has an overall
limited accuracy.

\plot{bin-win}{width=6in}{3-D data after normalized binning.}

For efficiency, we perform regularization on individual time slices.
Figure~\ref{fig:smo2-win} shows the result of regularization using
bi-linear interpolation and smoothing preconditioning (data-space
regularization) with the minimum-phase Laplacian filter \cite[]{wilsonburg}.
The empty bins are filled in a consistent manner but the data quality
is distorted because simple smoothing fails to characterize the
complicated data structure.  Instead of continuous events, we see
smoothed blobs in the time slices. The events in the in-line and
cross-line sections are also not clearly pronounced.

\plot{smo2-win}{width=6in}{3-D data regularized with bi-linear
  interpolation and smoothing preconditioning.}

We can use the smoothing regularization result to estimate the local
dips in the data, design invertible local plane-wave destruction
filters \cite[]{Fomel.sepphd.107}, and repeat the regularization process.  Inverse
interpolation with plane-wave data-space regularization is shown in
Figure~\ref{fig:int4-win}. The result is noticeably
improved: the continuous reflection events become clearly visible in
the time slices.  Despite the irregularities in the input data, the
regularization result preserves both flat reflection events and
steeply-dipping diffractions. Preserving diffractions is important for
correct imaging of sharp edges in the subsurface structure
\cite[]{GEO63-02-05740588}.

For simplicity, we assumed only a single local dip component in the
data. This assumption degrades the result in the areas of multiple
conflicting dips, such as the intersections of plane reflections and
hyperbolic diffractions in Figure~\ref{fig:int4-win}. One could
improve the image by considering multiple local dips.
\cite{ofcon2} describes an alternative offset-continuation
approach, which uses a physical connection between neighboring offsets
instead of assuming local continuity in the midpoint domain.

\plot{int4-win}{width=6in}{3-D data regularized with cubic B-spline
  interpolation and local plane-wave preconditioning.}

The 3-D results of this paper were obtained with an efficient 2-D
regularization in time slices. This approach is computationally
attractive because of its easy parallelization: different slices can
be interpolated independently and in parallel.
Figure~\ref{fig:winslice} shows the interpolation result for four
selected time slices. Local plane waves, barely identifiable after
binning (left plots in Figure~\ref{fig:winslice}), appear clear and
continuous in the interpolation result (right plots in
Figure~\ref{fig:winslice}). The time slices are assembled
together to form the 3-D cube shown in Figure~\ref{fig:int4-win}.

\plot{winslice}{width=6in}{Selected time slices of the 3-D dataset.
  Left: after binning. Right: after plane-wave data regularization.
  The data regularization program identifies and continues local
  plane waves in the data.}

\section{Conclusions}

Regularization is often a necessary part of geophysical estimation.
Its goal is to impose additional constraints on the model and
to guide the estimation process towards the desired solution.

We have considered two different regularization methods. The first,
model-space approach involves a convolution operator that enhances the
undesired features in the model. The second, data-space, approach
involves inverse filtering (deconvolution) to precondition the model.
Although the two approaches lead to the theoretically equivalent
results, their behavior in iterative estimation methods is quite
different. Using several synthetic and real data examples, we have
demonstrated that the second, preconditioning approach is generally
preferable because it shows a significantly faster convergence at 
early iterations.

We suggest a constructive method for preconditioning multidimensional
estimation problems using the helix transform. Applying inverse
filtering operators constructed this way, we observe a significant
(order of magnitude) speed-up in the optimization convergence. Since
inverse recursive filtering takes almost the same time as forward
convolution, the acceleration translates straightforwardly into
computational time savings.
\par
For simple test problems, these savings are hardly noticeable. On the
other hand, for large-scale (seismic-exploration-size) problems, the
achieved acceleration can have a direct impact on the mere feasibility
of iterative least-squares estimation.

\section{Acknowledgments}
We are grateful to Jim Berryman, Bill Harlan, Dave Nichols, Gennady
Ryzhikov, and Bill Symes for insightful discussions about
preconditioning and optimization problems.

The financial support for this work was provided by the sponsors of
the Stanford Exploration Project (SEP). The 3-D North Sea dataset was
released to SEP by Conoco and its partners, BP and Mobil. The SeaBeam
dataset is courtesy of Alistair Harding of the Scripps Institution of
Oceanography.

%\newpage
\bibliographystyle{seglike}
\bibliography{SEG,SEP2,sergey}

%\APPENDIX{A}

%\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
%%% TeX-master: t
%%% End: