paper.tex

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

% Started 06/28/00

%\shortnote

\lefthead{Fomel}
\righthead{Plane-wave destructors}
\footer{SEP--105}

\title{Applications of plane-wave destruction filters}

\email{[e-mail]}
%\keywords{}

\author{Sergey Fomel}

\maketitle

\begin{abstract}
  Plane-wave destruction filters originate from a local plane-wave
  model for characterizing seismic data. These filters can be thought
  of as a $T$-$X$ analog of $F$-$X$ prediction-error filters and as an
  alternative to $T$-$X$ prediction-error filters. The filters are
  constructed with the help of an implicit finite-difference scheme
  for the local plane-wave equation.  On several synthetic and
  real-data examples, I demonstrate that finite-difference plane-wave
  destruction filters perform well in applications such as fault
  detection, data interpolation, and noise attenuation.
\end{abstract}

\section{Introduction}

Plane-wave destruction filters, introduced by
\cite{Claerbout.blackwell.92}, serve the purpose of characterizing
seismic images by a superposition of local plane waves.  They are
constructed as finite-difference stencils for the plane-wave
differential equation. In many cases, a local plane-wave model is a
very convenient representation of seismic data. Unfortunately, early
experiences with applying plane-wave destructors for interpolating
spatially aliased data
\cite[]{Nichols.sep.65.271,Claerbout.blackwell.92} demonstrated their
poor performance in comparison with that of industry-standard $F$-$X$
prediction-error filters \cite[]{GEO56-06-07850794}.
\par
For each given frequency, an $F$-$X$ prediction-error filter (PEF) can
be thought of as a $Z$-transform polynomial. The roots of the
polynomial correspond precisely to predicted plane waves
\cite[]{SEG-1984-S10.1}.  Therefore, $F$-$X$ PEFs simply represent a
spectral (frequency-domain) approach to plane-wave
destruction\footnote{The filters are designed to \emph{destruct} local
  plane waves. However, in applications such as data interpolation,
  they are often used to \emph{reconstruct} the missing parts of
  local waves. The choice of terminology should not confuse the
  reader.}  This powerful and efficient approach is, however, not
theoretically adequate when the plane-wave slopes or the boundary
conditions vary both spatially and temporally. In practice, this 
limitation is addressed by breaking the data into windows and assuming
that the slopes are stationary within each window.
\par
Multidimensional $T$-$X$ prediction-error filters
\cite[]{Claerbout.blackwell.92,gee} share the same purpose of predicting
local plane waves. They work well with spatially aliased data and
allow for both temporal and spatial variability of the slopes.  In
practice, however, $T$-$X$ filters appear as very mysterious objects,
because their construction involves many non-intuitive parameters. The
user needs to choose a raft of parameters, such as the number of
filter coefficients, the gap and the exact shape of the filter, the
size, number, and shape of local patches for filter estimation, the
number of iterations, and the amount of regularization. Recently
developed techniques for handling non-stationary PEFs
\cite[]{SEG-1999-11541157} performed well in
a variety of applications \cite[]{Crawley.sepphd.104,SEG-2001-13051308}, but the large
number of adjustable parameters still requires a significant level of
human interaction and remains the drawback of the method.
\par
\cite{SEG-1998-1851} have recently revived the original
plane-wave destructors for preconditioning tomographic problems with a
predefined dip field \cite[]{Clapp.sepphd.106}. The filters were named
\emph{steering filters} because of their ability to steer the solution
in the direction of the local dips. The name is also reminiscent of
\emph{steerable filters} used in medical image processing \cite[]{steer0,steer}.
\par
In this paper, I revisit Claerbout's original technique of
finite-difference plane-wave destruction. First, I develop an approach
for increasing the accuracy and dip bandwidth of the method.  Applying
the improved filter design to several data regularization problems, I
discover that the finite-difference filters often perform as well as,
or even better than, $T$-$X$ PEFs.  At the same time, they keep the
number of adjustable parameters to a minimum, and the only estimated quantity 
has a clear physical meaning of the local plane-wave slope. 
No local windows are required, because the slope is estimated as a 
smoothly variable continuous function of the data coordinates.
\par
Conventional methods for estimating plane-wave slopes are based on
picking maximum values of stacking semblance and other cumulative
coherency measures \cite[]{GEO36-03-04820497}. The differential approach
to slope estimation, employed by plane-wave destruction filters, is
related to the differential semblance method \cite[]{GEO56-05-06540663}.
Its theoretical superiority to conventional semblance measures for the
problem of local plane wave detection has been established by
\cite{symes} and \cite{kim}.

\section{High-order plane-wave destructors}

Following the physical model of local plane waves, we can define the
mathematical basis of the plane-wave destruction filters as the local
plane differential equation
\begin{equation}
  \frac{\partial P}{\partial x} + 
  \sigma\,\frac{\partial P}{\partial t} = 0\;,
  \label{eqn:pde}
\end{equation}
where $P(t,x)$ is the wave field, and $\sigma$ is the local slope, which may
also depend on $t$ and $x$. In the case of a constant slope,
equation~(\ref{eqn:pde}) has the simple general solution
\begin{equation}
  P(t,x) = f(t - \sigma x)\;,
  \label{eqn:plane}
\end{equation}
where $f(t)$ is an arbitrary waveform. Equation~(\ref{eqn:plane}) is
nothing more than a mathematical description of a plane wave.
\par
If we assume that the slope $\sigma$ does not depend on $t$, we can
transform equation~(\ref{eqn:pde}) to the frequency domain, where it
takes the form of the ordinary differential equation
\begin{equation}
  {\frac{d \hat{P}}{d x}} +
  i \omega\,\sigma\, \hat{P} = 0
  \label{eqn:ode}
\end{equation}
and has the general solution
\begin{equation}
  \hat{P} (x) = \hat{P} (0)\,e^{i \omega\,\sigma x}\;,
  \label{eqn:px}
\end{equation}
where $\hat{P}$ is the Fourier transform of $P$. The complex
exponential term in equation~(\ref{eqn:px}) simply represents a shift
of a $t$-trace according to the slope $\sigma$ and the trace separation
$x$. 

In the frequency domain, the operator for transforming the trace at
position $x-1$ to the neighboring trace\footnote{For simplicity, it is
  assumed that $x$ takes integer values that correspond to trace
  numbering.} and at position $x$ is a multiplication by $e^{i
  \omega\,\sigma}$. In other words, a plane wave can be perfectly
predicted by a two-term prediction-error filter in the $F$-$X$ domain:
\begin{equation}
  a_0 \, \hat{P} (x) + a_1\, \hat{P} (x-1) = 0\;,
  \label{eqn:pef}
\end{equation}
where $a_0 = 1$ and $a_1 = - e^{i \omega\,\sigma}$. The goal of
predicting several plane waves can be accomplished by cascading
several two-term filters. In fact, any $F$-$X$ prediction-error
filter represented in the $Z$-transform notation as
\begin{equation}
  A(Z_x) = 1 + a_1 Z_x + a_2 Z_x^2 + \cdots + a_N Z_x^N
  \label{eqn:pef2}
\end{equation}
can be factored into a product of two-term filters:
\begin{equation}
  A(Z_x) = \left(1 - \frac{Z_x}{Z_1}\right)\left(1 - \frac{Z_x}{Z_2}\right)
  \cdots\left(1 - \frac{Z_x}{Z_N}\right)\;,
  \label{eqn:pef3}
\end{equation}
where $Z_1,Z_2,\ldots,Z_N$ are the zeroes of
polynomial~(\ref{eqn:pef2}). According to equation~(\ref{eqn:pef}),
the phase of each zero corresponds to the slope of a local plane wave
multiplied by the frequency. Zeroes that are not on the unit circle
carry an additional amplitude gain not included in
equation~(\ref{eqn:ode}).
\par
In order to incorporate time-varying slopes, we need to return to
the time domain and look for an appropriate analog of the phase-shift
operator~(\ref{eqn:px}) and the plane-prediction
filter~(\ref{eqn:pef}). An important property of plane-wave
propagation across different traces is that the total energy of the
propagating wave stays invariant throughout the process: the energy of 
the wave at one trace is completely transmitted to the next trace.
 This property
is assured in the frequency-domain solution~(\ref{eqn:px}) by the fact
that the spectrum of the complex exponential $e^{i \omega\,\sigma}$ is
equal to one.  In the time domain, we can reach an equivalent effect
by using an all-pass digital filter. In the $Z$-transform notation,
convolution with an all-pass filter takes the form
\begin{equation}
\hat{P}_{x+1}(Z_t) = \hat{P}_{x} (Z_t) \frac{B(Z_t)}{B(1/Z_t)}\;,
\label{eqn:allpass}
\end{equation}
where $\hat{P}_x (Z_t)$ denotes the $Z$-transform of the corresponding
trace, and the ratio $B(Z_t)/B(1/Z_t)$ is an all-pass digital filter
approximating the time-shift operator~$e^{i \omega \sigma}$. In
finite-difference terms, equation~(\ref{eqn:allpass}) represents an
implicit finite-difference scheme for solving equation~(\ref{eqn:pde})
with the initial conditions at a constant $x$.  The coefficients of
filter $B(Z_t)$ can be determined, for example, by fitting the filter
frequency response at low frequencies to the response of the
phase-shift operator. The Taylor series technique (equating the
coefficients of the Taylor series expansion around zero frequency)
yields the expression
\begin{equation}
  B_3(Z_t) = 
  \frac{(1-\sigma)(2-\sigma)}{12}\,Z_t^{-1} + 
  \frac{(2+\sigma)(2-\sigma)}{6} +
  \frac{(1+\sigma)(2+\sigma)}{12}\,Z_t
  \label{eqn:b3}
\end{equation}
for a three-point centered filter $B_3(Z_t)$ and the expression
\begin{eqnarray}
  B_5(Z_t) & = &  
  \frac{(1-\sigma)(2-\sigma)(3-\sigma)(4-\sigma)}{1680}\,Z_t^{-2} +
  \frac{(4-\sigma)(2-\sigma)(3-\sigma)(4+\sigma)}{420}\,Z_t^{-1} + 
  \nonumber \\
  & & 
  \frac{(4-\sigma)(3-\sigma)(3+\sigma)(4+\sigma)}{280} + 
  \nonumber \\
  & & 
  \frac{(4-\sigma)(2+\sigma)(3+\sigma)(4+\sigma)}{420}\,Z_t +
  \frac{(1+\sigma)(2+\sigma)(3+\sigma)(4+\sigma)}{1680}\,Z_t^2
  \label{eqn:b5}
\end{eqnarray} 
for a five-point centered filter $B_5(Z_t)$. The derivation of
equations~(\ref{eqn:b3}-\ref{eqn:b5}) is detailed in the appendix. It
is easy to generalize these equations to longer filters.

\inputdir{Math}

Figure~\ref{fig:phase} shows the phase of the all-pass filters
$B_3(Z_t)/B_3(1/Z_t)$ and $B_5(Z_t)/B_5(1/Z_t)$ for two values of the
slope $\sigma$ in comparison with the exact linear function of
equation~(\ref{eqn:px}).  As expected, the phases match the exact line
at low frequencies, and the accuracy of the approximation increases
with the length of the filter.

\plot{phase}{width=6in}{Phase of the implicit
  finite-difference shift operators in comparison with the exact
  solution. The left plot corresponds to the slope of 
  $\sigma=0.5$, the right plot
  to $\sigma=0.8$.}
\par
Taking both dimensions into consideration,
equation~(\ref{eqn:allpass}) transforms to the prediction equation
analogous to~(\ref{eqn:pef}) with the 2-D prediction filter
\begin{equation}
  A(Z_t,Z_x) = 1 - Z_x \frac{B(Z_t)}{B(1/Z_t)}\;.
  \label{eqn:2dpef}
\end{equation}
In order to characterize several plane waves, we can cascade several
filters of the form~(\ref{eqn:2dpef}) in a manner similar to that of
equation~(\ref{eqn:pef3}). In the examples of this paper, I use a
modified version of the filter $A(Z_t,Z_x)$, namely the filter
\begin{equation}
  \label{eqn:2dpef2}
  C(Z_t,Z_x) = A(Z_t,Z_x) B(1/Z_t) = B(1/Z_t) - Z_x B(Z_t)\;,
\end{equation}
which avoids the need for polynomial division. In case of the 3-point
filter~(\ref{eqn:b3}), the 2-D filter~(\ref{eqn:2dpef2}) has exactly