paper.tex

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

% Started 09/15/97

%\shortnote

\lefthead{Fomel \& Claerbout}
\righthead{Implicit extrapolation}

\title{Exploring three-dimensional implicit wavefield extrapolation
  with the helix transform}

\email{[e-mail], [e-mail]}
%\keywords{depth migration, post-stack, finite-difference, velocity continuation, helix}

\author{Sergey Fomel and Jon F. Claerbout}

\maketitle

\begin{abstract}
  Implicit extrapolation is an efficient and unconditionally stable
  method of wavefield continuation. Unfortunately, implicit wave
  extrapolation in three dimensions requires an expensive solution of
  a large system of linear equations. However, by mapping the
  computational domain into one dimension via the helix transform, we
  show that the matrix inversion problem can be recast in terms of an
  efficient recursive filtering. Apart from the boundary conditions,
  the solution is exact in the case of constant coefficients (that is,
  a laterally homogeneous velocity.) We illustrate this fact with an
  example of three-dimensional velocity continuation and discuss
  possible ways of attacking the problem of lateral variations.
\end{abstract}

\section{Introduction}

Implicit finite-difference wavefield extrapolation played an
exceptionally important role in the early development of seismic
migration methods. Using limited-degree approximations to the one-way
wave equation, implicit schemes have provided efficient and
unconditionally stable numerical wave extrapolation operators
\cite[]{Godfrey.sep.16.83,Claerbout.blackwell.85}. Unfortunately, the
advantages of \emph{implicit} methods were lost with the development
of three-dimensional seismic exploration. While the cost of 2-D
implicit extrapolation is linearly proportional to the mesh size, the
same approach, applied in the 3-D case, leads to a nonlinear
computational complexity. Primarily for this reason, implicit
extrapolators were replaced in practice by \emph{explicit} ones,
capable of maintaining linear complexity in all dimensions. A number
of computational tricks \cite[]{GEO56-11-17701777} allow the commonly
used explicit schemes to behave stably in practical cases.  However,
their stability is not unconditional and may break in unusual
situations \cite[]{SEG-1994-1266}.
\par
In this paper, we present an approach to three-dimensional
extrapolation, based on the helix transform of multidimensional
filters to one dimension \cite[]{Claerbout.gem.97}. The traditional
approach involves an inversion of a banded matrix (tridiagonal in the
2-D case and blocked-tridiagonal in the 3-D case). With the help of
the helix transform, we can recast this problem in terms of inverse
recursive filtering.  The coefficients of two-dimensional filters on a
helix are obtained by one-dimensional spectral factorization methods.
As a result, the complexity of three-dimensional implicit
extrapolation is reduced to a linear function of the computational
mesh size. This approach doesn't provide an exact solution in the
presence of lateral velocity variations. Nevertheless, it can be used
for preconditioning iterative methods, such as those described by
\cite{Nichols.sep.70.31}.  In this paper, we demonstrate the
feasibility of 3-D implicit extrapolation on the example of laterally
invariant velocity continuation and, in the final part, discuss
possible strategies for solving the problem of lateral variations.
\par
The main application of finite-difference wave extrapolation is
\emph{post-stack} depth migration. An application of similar methods
for \emph{prestack} common-shot migration is constrained by the
limited aperture of commonly used seismic acquisition patterns.
Recently developed acquisition methods, such as the vertical cable
technique \cite[]{SEG-1993-1376}, open up new possibilities for 3-D wave
extrapolation applications. An alternative approach is common-azimuth
migration \cite[]{Biondi.sep.80.109,Biondi.sep.93.1}. Other interesting
applications include finite-difference data extrapolation in offset
\cite[]{Fomel.sep.84.179}, migration velocity \cite[]{Fomel.sep.92.159},
and anisotropy \cite[]{Alkhalifah.sep.94.tariq3}.

\section{Implicit versus Explicit extrapolation}

The difference between implicit and explicit extrapolation is best
understood through an example. Following \cite{Claerbout.blackwell.85},
let us consider, for instance, the diffusion (heat conduction) equation
of the form
\begin{equation}
  {\frac{\partial T}{\partial t}} = {a (x)\,{\frac{\partial^2 T}{\partial x^2}}}\;.
\label{eqn:heat}
\end{equation}
Here $t$ denotes time, $x$ is the space coordinate, $T (x,t)$ is the
temperature, and $a$ is the heat conductivity coefficient.
Equation (\ref{eqn:heat}) forms a well-posed boundary-value problem if
supplied with the initial condition
\begin{equation}
  \label{eqn:heatinit}
  \left.T\right|_{t=0} = T_0 (x)
\end{equation}
and the appropriate boundary conditions. Our task is to build a
digital filter, which transforms a gridded temperature $T$ from one
time level to another.
\par
It helps to note that when the conductivity coefficient $a$ is
constant and the space domain of the problem is infinite (or periodic)
in $x$, the problem can be solved in the wavenumber domain. Indeed,
after the Fourier transform over the variable $x$, equation
(\ref{eqn:heat}) transforms to the ordinary differential equation
\begin{equation}
  {\frac{d \hat{T}}{d t}} = {- a k^2\, \hat{T}}\;,
  \label{eqn:heatk}	
\end{equation}
which has the explicit analytical solution
\begin{equation}
  \label{eqn:heatsol}
  \hat{T} (k,t) = \hat{T}_0 (k) e^{- a k^2 t}\;,
\end{equation}
where $\hat{T}$ denotes the Fourier transform of $T$, and $k$ stands
for the wavenumber. Therefore, the desired filter in the
wavenumber domain has the form
\begin{equation}
  \label{eqn:heatf}
  H (k) = e^{- a k^2}\;,
\end{equation}
where for simplicity the coefficient $a$ is normalized for the time
step $\triangle t$ equal to $1$.
\par
Returning now to the time-and-space domain, we can approach the filter
construction problem by approximating the space-domain response of
filter (\ref{eqn:heatf}) in terms of the differential operators
$\frac{\partial^2}{\partial x^2} = - k^2$, which can be approximated
by finite differences. An \emph{explicit} approach would amount to
constructing a series expansion of the form
\begin{equation}
  \label{eqn:heatexpl}
  H_{\mbox{ex}} (k) \approx a_0 + a_1 k^2 + a_2 k^4 + \ldots\;,
\end{equation}
and selecting the coefficients $a_j$ to approximate equation
(\ref{eqn:heatf}). For example, the three-term Taylor series expansion
around the zero wavenumber yields
\begin{equation}
  \label{eqn:heattayl}
  H_{\mbox{ex}} (k) = 1 - a\,{k^2}  +
  {\frac{{{a }^2}\,{k^4}}{2}} \;.
\end{equation}
The error of approximation (\ref{eqn:heattayl}) as a function of $k$
for two different values of $a$ is shown in the left plot of Figure
\ref{fig:error}.

\inputdir{Math}

\plot{error}{width=6in}{Errors of second-order explicit and implicit
  approximations for the heat extrapolation.}
\par
An \emph{implicit} approach also approximates the ideal filter
(\ref{eqn:heatf}), but with a rational approximation of the form
\begin{equation}
  \label{eqn:heatpade}
  H_{\mbox{im}} (k) \approx \frac{b_0 + b_1 k^2 + b_2 k^4 + \ldots}
  {1 + c_1 k^2 + c_2 k^4 + \ldots}
\;.
\end{equation}
One way of selecting the coefficients $b_i$ and $c_i$ is to apply an
appropriate Pad\'{e} approximation \cite[]{pade}\footnote{If the
  denominator and the numerator have the same order, Pad\'{e}
  approximants are equivalent to the corresponding continuous
  fraction expansions.}.  For example the $[2/2]$ Pad\'{e}
approximation is
\begin{equation}
  \label{eqn:heatcrank}
  H_{\mbox{im}} (k) =
  \frac{1 - \frac{a}{2}\,k^2}{1 + \frac{a}{2}\,k^2}
  \;.
\end{equation}
This approximation corresponds to the famous Crank-Nicolson implicit
method \cite[]{cn}. The error of approximation (\ref{eqn:heatcrank}) as
a function of $k$ for different values of $a$ is shown in the right
plot of Figure \ref{fig:error}. Not only is it significantly smaller
than the error of the same-order explicit approximation, but it also
has a negative sign. It means that the high-frequency numerical noise
gets suppressed rather than amplified. In practice, this property
translates into a stable numerical extrapolation.
\par
The second derivative operator $-k^2$ can be approximated in practice
by a digital filter. The most commonly used filter has the
$Z$-transform $D_2 (Z) = -Z^{-1} + 2 - Z$, and the Fourier transform
\begin{equation}
  \label{eqn:d2k}
  D_2 (k) = e^{-ik} - 2 + e^{-ik} = 2 (\cos{k} - 1) = -4
  \sin^2{\frac{k}{2}}\;.
\end{equation}
Formula (\ref{eqn:d2k}) approximates $-k^2$ well only for small
wavenumbers $k$. As shown in Appendix A, the implicit scheme allows
the accuracy of the second-derivative filter to be significantly
improved by a variation of the ``1/6-th trick''
\cite[]{Claerbout.blackwell.85}. The final form of the implicit
extrapolation filter is
\begin{equation}
  \label{eqn:heatfk}
   H_{\mbox{im}} (k) =
   \frac{1 + \frac{a+\beta}{2}\,D_2 (k)}{1 - \frac{a-\beta}{2}\,D_2 (k)}
   \;,
\end{equation}
where $\beta$ is a numerical constant, found in Appendix A.

\inputdir{heat}

\plot{heat}{width=6in,height=2.5in}{Heat extrapolation with explicit
  and implicit finite-different schemes. Explicit extrapolation
  appears stable for $a=2/3$ (left plot) and unstable for $a=4/3$
  (middle plot). Implicit interpolation is stable even for larger
  values of $a$ (right plot).}
\par
A numerical 1-D example is shown in Figure \ref{fig:heat}. The initial
temperature distribution is given by a step function. The
discontinuity at the step gets smoothed with time by the heat
diffusion. The left plot shows the result of an explicit extrapolation
with $a=2/3$, which appears stable. The middle plot is an explicit
extrapolation with $a=4/3$, which shows a terribly unstable behavior:
the high-frequency numerical noise is amplified and dominates the
solution. The right plot shows a stable (though not perfectly
accurate) extrapolation with the implicit scheme for the larger value of
$a=2$.
\par
The difference in stability between explicit and implicit schemes is
even more pronounced in the case of \emph{wave extrapolation}. For
example, let us consider the ideal depth extrapolation filter in the
form of the phase-shift operator
\cite[]{GEO43-07-13421351,Claerbout.blackwell.85}
\begin{equation}
  \label{eqn:gazdag}
  W (k) = e^{i \sqrt{a^2 - k^2}}\;,
\end{equation}
where $a = \omega / v$, $\omega$ is the time frequency, and $v$ is the
seismic velocity (which may vary spatially); we assume for simplicity
that both the depth step $\triangle z$ and the space sampling
$\triangle x$ are normalized to $1$.  A simple implicit approximation
to filter (\ref{eqn:gazdag}) is
\begin{equation}
  \label{eqn:wave45}
   W_{\mbox{im}} (k) = e^{i a}\,
   \frac{1 -4 a^2 + i a\,k^2}{1 - 4 a^2 - i a\,k^2} = e^{i \phi}\;,
\end{equation}
where $\phi = a - 2 \arctan{\frac{a\,k^2}{4 a^2-1}}$. We can see
that approximation (\ref{eqn:wave45}) is again a pure phase shift
operator, only with a slightly different phase. For that reason, the
operator is unconditionally stable for all values of $a$: the total
wave energy from one depth level to another is preserved. Operator
(\ref{eqn:gazdag}) corresponds to the Crank-Nicolson scheme for the
45-degree one-way wave equation \cite[]{Claerbout.blackwell.85}. Its
phase error as a function of the dip angle $\theta =
\arcsin{\frac{k}{a}}$ for different values of $a$ is shown in Figure
\ref{fig:phase}.

\inputdir{Math}