paper.tex

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

\footer{SEP--89}

\title{Antialiasing of Kirchhoff operators by reciprocal
  parameterization}

\author{Sergey Fomel}
\maketitle

\input{intro}

\section{Overview of existing methods}

I start with reviewing the existing approaches to operator
antialiasing and discussing their main principles and limitations. The
two reviewed approaches are temporal filtering, as suggested by
\cite{GPR40-05-05650572} and \cite{SEG-1994-1282}, and Hale's spatial
filtering technique, developed originally for an integral
implementation of the dip moveout operator \cite[]{GEO56-06-07950805}.

\subsection{Temporal filtering}
%%%%%%%%%%%%%%%%%%%%%%%%
The temporal filtering idea follows
from the well-known Nyquist sampling criterion. With application to
integral operators, the Nyquist criterion takes the form
\begin{equation}
\Delta x \leq {{\Delta t} \over {|\partial t / \partial x|}}\;,
\label{eqn:Nyquist}
\end{equation}
where $t(x)$ is the traveltime of the operator impulse response,
$\Delta x$ is the space sampling interval and $\Delta t$ is the time
sampling interval. In the steep parts of the traveltime curve, the
sampling criterion (\ref{eqn:Nyquist}) is not satisfied, which causes
aliasing artifacts in the output data. To overcome this problem, the
method of local triangle filtering
\cite[]{Claerbout.sep.73.371,SEG-1994-1282} suggests convolving the
traces of the generated impulse response with a triangle-shaped filter
of the length
\begin{equation}
\delta t = \Delta x\,|\partial t / \partial x|\;.
\label{eqn:dt}
\end{equation}

\inputdir{XFig}

Cascading operators of causal and anticausal numerical integration is
an efficient way to construct the desired filter shape.

Triangle filters approximate the ideal (sinc) low-pass time filters.
The idea of using low-pass filtering for antialiasing
\cite[]{GPR40-05-05650572} is illustrated in Figure \ref{fig:amolow}.
When a steeply dipping event is included in the operator, its
counterpart in the frequency domain wraps around to produce the
aliasing artifacts. Those are removed by a dip-dependent low-pass
filtering.

\plot{amolow}{width=4.5in,height=1.5in}{Schematic
illustration of low-pass antialiasing (triangle filters). The aliased
events are removed by low-pass filtration on the temporal frequency axis.
The width of the low-pass filter depends on dips of the aliased events.}

\inputdir{flt}

The method of low-pass filtering is less evident in the case of a
three-dimensional integral operator. We can take the length of a
triangle filter proportional to the absolute value of the time
gradient \cite[]{SEG-1994-1282}, the maximum of the gradient
components in the two directions of the operator space
\cite[]{GEO64-06-17831792}, or the sum of these components. The latter
follows from considering the 3-D operator as a double integration in
space. Decoupling the 3-D integral into a cascade of two 2-D operators
suggests convolving two triangle filters designed with respect to two
coordinates of the operator. In this case, the length of the resultant
filter is approximately equal to
\begin{equation}
\delta t =      \Delta x\,|\partial t / \partial x| + 
                \Delta y\,|\partial t / \partial y|\;,
\label{eqn:dt3}
\end{equation}
and its shape is smoother than that of a triangle filter (Figure
\ref{fig:amoflt}).

\plot{amoflt}{height=2.5in}{Building the smoothed filter for
3-D antialiasing by successive integration of a five-point wavelet. C
denotes the operator of causal integration, C' denotes its adjoint (the
anticausal integration). The result is equivalent to the convolution
of two equal triangle filters.}

The temporal filtering method was proven to be an efficient tool in
the design of stacking operators of different types. However, when the
operator introduces rapid changes in the length and direction of the
traveltime gradient, it leads to an inexact estimation of the filter
cutoff (triangle length for the method of triangle filtering) at the
curved parts of the operator. Consequently, the high-frequency part of
the output can be distorted, causing a loss in the image resolution.

\subsection{Hale's method}
%%%%%%%%%%%%%%%%%%%%

Considering the case of integral dip moveout, \cite{GEO56-06-07950805}
points out that the steep parts of the operator, while aliased in the
space (midpoint) coordinate, are not aliased with respect to the time
coordinate. He suggests replacing the conventional $t(x)$
parameterization of the DMO impulse response by $x(t)$
parameterization. Conventionally, the integral operators are
implemented by shifting the input traces in space and transforming
them in time.  According to Hale's method, the traces are shifted in
time and transformed along the $x(t)$ trajectories in space.
Interpolation in time, required in the conventional approach, is
replaced by interpolation in space. The idea of Hale's method is
related to the idea of the ``pixel-precise velocity transform''
\cite[]{Claerbout.blackwell.92}.

The steep parts of the operator satisfy the criterion
\begin{equation}
\Delta t \leq {{\Delta x} \over {|\partial x / \partial t|}}\;,
\label{eqn:Nyquist2}
\end{equation}
which is the the obvious reverse of inequality (\ref{eqn:Nyquist}).
Therefore, they are not aliased if defined on the time grid. In these
parts, one can implement the operator by constant time shifts equal to
the time sampling interval $\Delta t$. In the parts where the
criterion (\ref{eqn:Nyquist2}) is not valid (the flat part of the DMO
operator), Hale suggests reducing the length of the time shifts
according to equality (\ref{eqn:dt}), where $\delta t$ becomes less
than $\Delta t$. He formulates the following principle of operator
antialiasing:
\begin{quote}
  To eliminate spatial aliasing, simply never allow successive time
shifts applied to the input trace to differ by more than one time
sampling interval. Further restrict the difference between time shifts
so that the spacing between the corresponding output trajectories
never exceeds the CMP sampling interval.
\end{quote}

\inputdir{XFig}

The idea of Hale's method is illustrated in Figure \ref{fig:amosft}.
Increasing the density of spatial sampling by small successive time
shifts implies increasing the Nyquist boundaries of the spatial
wavenumber. Further interpolation is a low-pass spatial filtering that
removes the parts of the spectrum beyond the Nyquist frequency of the
output. If the dip of the operator does not vary between neighboring
traces (the operator is a straight line as in the slant stack case),
Hale's approach will produce essentially the same result as that of
temporal filtering. Triangle filters in this case approximately
correspond to linear interpolation in space between adjacent traces
\cite[]{Nichols.sep.77.283}. The difference between the two approaches
occurs if the local dip varies in space as in the case of a curved
operator, such as DMO. In this case, Hale's approach provides a more
accurate space interpolation of the operator and preserves the
high-frequency part of its spectrum from distortion.

\plot{amosft}{width=4.125in,height=1.5in}{Schematic
illustration of Hale's antialiasing. The aliased events are removed by
spatial interpolation. In the frequency domain, the interpolation
consists of widening and low-passing on the wavenumber axis. The
low-pass spatial filtering does not depend on dip.} 

Hale's method has proven to preserve the amplitude of flat reflectors
from aliasing distortions, which is the simplest antialiasing test on
a DMO operator. The most valuable advantage of this method in the fact
that the implied low-pass spatial filtering (interpolation) does not
depend on the operator dip and is controlled by the Nyquist boundary
of the spectrum only (compare Figures \ref{fig:amolow} and
\ref{fig:amosft}). This is especially important, when the local dip of
the operator changes rapidly and therefore cannot be estimated
precisely by finite-difference approximation at spatially separated
traces. Such a situation is common in dip moveout and azimuth moveout
integral operators, as well as in prestack Kirchhoff migration.

A weakness of the method is the necessity to switch from
interpolation in space to two-dimensional interpolation in both the
time and the space variables, when trying to construct the flat part
of the operator. 
%In the case of AMO, the 2-D spatial interpolation
%arises as a result of building the operator in the transformed
%coordinate system. 
In the next section, I show how to avoid the expense of the additional
time interpolation required by Hale's method of antialiasing.

\section{Proposed technique}
%%%%%%%%%%%%%%%%%%%%%%%%%%
We can use the reciprocity of the time parameterization and the space
parameterization of integral operators, discovered by Hale, to arrive
at the following antialiasing technique.

For simplicity, let us consider the two-dimensional case first.  The
linearity of a two-dimensional integral operator allows us to
decompose this operator into two parts. The first part corresponds to
the steep part of the travel-time function, which satisfies the
time-sampling criterion (\ref{eqn:Nyquist2}). The second term
corresponds to the flat part of the traveltime, which satisfies the
space-sampling criterion (\ref{eqn:Nyquist}). The first part is not
aliased with respect to the time sampling interval, while the second
one is not aliased with respect to the space sampling. We can apply
interpolation in time to construct the flat part.  Reciprocally,
interpolation in space is applied to construct the steep part of the
operator in the fashion of Hale's time-shifting method. 
%Linear
%interpolation in this case is a cheap substitution for the accurate
%but computationally expensive sinc interpolation. 
The amplitude
difference between the two integrals is simply the Jacobian term