paper.tex

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

\sideplot{phase}{width=2.5in}{The phase error of the
  implicit depth extrapolation with the Crank-Nicolson method.}
\par
The unconditional stability property is not achievable with the
explicit approach, though it is possible to increase the stability of
explicit operators by using relatively long filters
\cite[]{GPR36-02-00990114,GEO56-11-17701777}.

\section{Spectral factorization and three-dimensional extrapolation}

In this section, we continue our review of extrapolation methods to
reveal the principal difficulties of three-dimensional extrapolation.
We then describe a new, helix-transform approach to this old and
fascinating problem.

\subsection{Inverse filter factorization}
The conventional way of applying implicit finite-difference schemes
reduces to solving a system of linear equations with a sparse matrix.
For example, to apply the scheme of equation (\ref{eqn:heatfk}), we
can put the filter denominator on the other side of the extrapolation
equation, writing it as
\begin{equation}
  \label{eqn:heattris}
  \left(\mathbf{I} - \frac{a-\beta}{2}\,\mathbf{D}_2\right) \mathbf{T}_{t+1} =
  \left(\mathbf{I} + \frac{a+\beta}{2}\,\mathbf{D}_2\right) \mathbf{T}_t\;,
\end{equation}
where $\mathbf{I}$ is the identity matrix, $\mathbf{D}$ is the convolution
matrix for filter (\ref{eqn:d2k}), and $\mathbf{T}_t$ is the vector of
temperature distribution at time level $t$. In the case of
two-dimensional extrapolation, the matrix on the left side of equation
(\ref{eqn:heattris}) takes the tridiagonal form
\begin{equation}\label{eqn:matrix}
  \mathbf{A} = \left(\mathbf{I} -c\,\mathbf{D}_2\right) =
  \left[\begin{array}{cccccc}
  1+2 c_{1}   & -c_{1}     &  0     & \cdots &        & 0      \\
  -c_{2}     & 1+2c_{2}   & -c_{2}     & 0      & \cdots &        \\
  0      & -c_{3}     & \cdots & \cdots &        & \cdots \\
  \cdots & 0      & \cdots &        &        &        \\
         & \cdots &        &        & \cdots & -c_{n-1}     \\
  0      &        &        &        & -c_{n}     & 1 + 2c_{n}
  \end{array}\right]\;,
\end{equation}
where $c = \frac{a-\beta}{2}$, and where, for simplicity, we assume
zero-slope boundary conditions. Like any positive-definite tridiagonal
matrix, matrix $\mathbf{A}$ can be inverted recursively by an $LU$
decomposition into two bidiagonal matrices. The cost of inversion is
directly proportional to the number of vector components.  The same
conclusion holds for the case of depth extrapolation [equation
(\ref{eqn:wave45})] with the substitution $c = \frac{\beta + i
  a}{1-4\,a^2}$.
\par
In the case of a laterally constant coefficient $a$, we can take a
different point of view on the tridiagonal matrix inversion. In this
case, the matrix $\mathbf{A}_2$ represents a convolution with a
symmetric three-point filter $1-c\,D_2(k)$. The $LU$
decomposition of such a matrix is precisely equivalent to filter
\emph{factorization} into the product of a causal minimum-phase filter
with its adjoint. This conclusion can be confirmed by the easily
verified equality
\begin{equation}\label{eqn:d2d1}
  1 + c (Z^{-1} - 2 + Z) = \frac{(1+b)^2}{4}\, \left(1 + \frac{1-b}{1+b} Z\right)
  \,\left(1 + \frac{1-b}{1+b} Z^{-1}\right)\;,
\end{equation}
where $b = \sqrt{1+ 4\,c}$. The inverse of the causal minimum-phase
filter $1 + \frac{1-b}{1+b} Z$ is a recursive inverse filter.
Correspondingly, the inverse of its adjoint pair, $1 + \frac{1-b}{1+b}
Z^{-1}$, is the same inverse filtering, performed in the adjoint mode
(backwards in space).  In the next subsection, we show how this
approach can be carried into three dimensions by applying the helix
transform.

\subsection{Helix and multidimensional deconvolution}

The major obstacle of applying an implicit extrapolation in three
dimensions is that the inverted matrix is no longer tridiagonal. If we
approximate the second derivative (Laplacian) on the 2-D plane with
the commonly used five-point filter $Z_x^{-1} + Z_y^{-1} -4 + Z_x +
Z_y$, then the matrix on the left side of equation
(\ref{eqn:heattris}), under the usual mapping of vectors from a
two-dimensional mesh to one dimension, takes the infamous
blocked-tridiagonal form \cite[]{birk}
\begin{equation}\label{eqn:matrix2}
  \tilde{\mathbf{A}} = \left(\mathbf{I} -c\,\mathbf{D}_2\right) =
  \left[\begin{array}{cccccc}
  \mathbf{A}_1   & -c_1 \,\mathbf{I}    &  0     & \cdots &        & 0      \\
  -c_2 \,\mathbf{I}    & \mathbf{A}_2   & -c_2 \,\mathbf{I}    & 0      & \cdots &        \\
  0      & -c_3 \,\mathbf{I}     & \cdots & \cdots &        & \cdots \\
  \cdots & 0      & \cdots &        &        &        \\
         & \cdots &        &        & \cdots & -c_{n-1} \,\mathbf{I}    \\
  0      &        &        &        & -c_{n} \,\mathbf{I}    &
  \mathbf{A}_n
  \end{array}\right]\;.
\end{equation}
Inspecting this form more closely, we see that the main diagonal of
$\tilde{\mathbf{A}}$, as well as the two offset diagonals formed by the
scaled identity matrices, remains continuous, while the second top and
bottom diagonals are broken. Therefore, even for constant $c$, the
inverted matrix does not have a simple convolutional structure, and
the cost of its inversion grows nonlinearly with the number of grid
points.
\par
A \emph{helix transform}, recently proposed by one of the authors
\cite[]{Claerbout.sep.95.jon1}, sheds new light on this old problem.
Let us assume that the extrapolation filter is applied by sliding it
along the $x$ direction in the $\{x,y\}$ plane. The diagonal
discontinuities in matrix $\tilde{\mathbf{A}}$ occur exactly in the
places where the forward leg of the filter slides outside the
computational domain. Let's imagine a situation, where the leg of the
filter that went to the end of the $x$ column, would immediately
appear at the beginning of the next column. This situation defines a
different mapping from two computational dimensions to the one
dimension of linear algebra. The mapping can be understood as the
helix transform, illustrated in Figure \ref{fig:helix1} and explained
in detail by \cite{Claerbout.sep.95.jon1}. According to this
transform, we replace the original two-dimensional filter with a long
one-dimensional filter. The new filter is partially filled with zero
values (corresponding to the back side of the helix), which can be
safely ignored in the convolutional computation.

\plot{helix}{width=5in,bb=210 155 630 390}{The helix transform of
  two-dimensional filters to one dimension. The two-dimensional filter
  in the left plot is equivalent to the one-dimensional filter in the
  right plot, assuming that a shifted periodic condition is imposed on
  one of the axes.}
\par
This is exactly the helix transform that is required to make all the
diagonals of matrix $\tilde{\mathbf{A}}$ continuous. In the case of
laterally invariant coefficients, the matrix becomes strictly Toeplitz
(having constant values along the diagonals) and represents a
one-dimensional convolution on the helix surface. Moreover, this
simplified matrix structure applies equally well to larger
second-derivative filters ( such as those described in Appendix B),
with the obvious increase of the number of Toeplitz diagonals.
Inverting matrix $\tilde{\mathbf{A}}$ becomes once again a simple
inverse filtering problem.  To decompose the 2-D filter into a pair
consisting of a causal minimum-phase filter and its adjoint, we can
apply spectral factorization methods from the 1-D filtering theory
\cite[]{Claerbout.blackwell.76,Claerbout.blackwell.92}, for example,
Kolmogorov's highly efficient method \cite[]{kolmog}. Thus, in the case
of a laterally invariant implicit extrapolation, matrix inversion
reduces to a simple and efficient recursive filtering, which we need
to run twice: first in the forward mode, and second in the adjoint
mode.

\inputdir{heat}

\plot{heat3d}{width=6in,height=2in}{Heat extrapolation in two
  dimensions, computed by an implicit scheme with helix recursive
  filtering. The left plot shows the input temperature distributions;
  the two other plots, the extrapolation result at different time
  steps.  The coefficient $a$ is 2.}
\par
Figure \ref{fig:heat3d} shows the result of applying the helix
transform to an implicit heat extrapolation of a two-dimensional
temperature distribution. The unconditional stability properties are
nicely preserved, which can be verified by examining the plot of
changes in the average temperature (Figure \ref{fig:heat-mean}).

\sideplot{heat-mean}{width=2.5in}{Demonstration of the stability of
  implicit extrapolation. The solid curve shows the normalized mean
  temperature, which remains nearly constant throughout the
  extrapolation time. The dashed curve shows the normalized maximum
  value, which exhibits the expected Gaussian shape.}
\par
In principle, we could also treat the case of a laterally invariant
coefficient with the help of the Fourier transform. Under what
circumstances does the helix approach have an advantage over Fourier
methods? One possible situation corresponds to a very large input data
size with a relatively small extrapolation filter. In this case, the
$O(N log N)$ cost of the fast Fourier transform is comparable with the
$O(N_f N)$ cost of the space-domain deconvolution (where $N$
corresponds to the data size, and $N_f$ is the filter size). Another
situation is when the boundary conditions of the problem have an
essential lateral variation. The latter case may occur in applications
of velocity continuation, which we discuss in the next section.  Later
in this paper, we return to the discussion of problems associated with
lateral variations.

\section{Three-dimensional implicit velocity continuation}
\inputdir{vc3}

Velocity continuation is a process of navigating in the migration
velocity space, applicable for time migration, residual migration, and
migration velocity analysis \cite[]{Fomel.sep.92.159}. In the
zero-offset (post-stack) case, the velocity continuation process is
described by the simple partial differential equation
\cite[]{Claerbout.sep.48.79,me}
\begin{equation}
  {\frac{\partial^2 P}{\partial v\,\partial t}} +
  {v\,t\,\left({\frac{\partial^2 P}{\partial x^2}} + {\frac{\partial^2
      P}{\partial y^2}}\right)} = 0\;,
\label{eqn:velcon}
\end{equation}
where $t$ is the vertical time coordinate of the migrated image, $x$
and $y$ are spatial (midpoint) coordinates, and $v$ is the migration
velocity. Slightly different versions of two-dimensional implicit
extrapolation with equation (\ref{eqn:velcon}) have been described by
\cite{Li.sep.48.85} and \cite[]{Fomel.sep.92.159}. 

\plot{velcon}{width=6in}{Impulse responses of the velocity
  continuation operator, computed by an implicit, unconditionally
  stable extrapolation via the helix transform. The left plot
  corresponds to continuation towards higher velocities (migration
  mode); the right plot, smaller velocities (modeling mode).}
\par
The helix approach has allowed us to modify the old code for three
dimensions. Figure \ref{fig:velcon} shows impulse responses of an
implicit helix-based three-dimensional velocity continuation.

\sideplot{qdome}{width=3in}{Qdome synthetic model, used for testing
  the 3-D velocity continuation program.}
\par
Figure \ref{fig:modmig} illustrates the velocity continuation process
on the Qdome synthetic model \cite[]{Claerbout.gem.97}, shown in Figure
\ref{fig:qdome}. Continuation backward in velocity corresponds to the
``modeling'' mode, while forward continuation corresponds to the
``migration'' mode. It is possible to balance the amplitudes of the
two processes so that the finite-difference velocity continuation
behaves as a unitary operator
\cite[]{Fomel.sep.92.159,Fomel.sep.92.267}.

\plot{modmig}{width=6in}{Modeling (left) and migration (right) with
  the Qdome synthetic model, obtained by running the 3-D velocity
  continuation backward and forward in velocity.}

\section{Depth extrapolation and the v(x) challenge}

Can the constant-velocity result help us achieve the challenging goal
of a stable implicit depth extrapolation through media with lateral
velocity variations?
\par
The first idea that comes to mind is to replace the space-invariant
helix filters with a precomputed set of spatially varying filters,
which reflect local changes in the velocity fields. This approach
would merely reproduce the conventional practice of explicit depth
extrapolators, popularized by \cite{GPR36-02-00990114} and
\cite{GEO56-11-17701777}. However, it hides the danger of losing
the property of unconditional stability, which is obviously the major
asset of implicit extrapolators.
\par
Another route, partially explored by \cite{Nichols.sep.70.31}, is
to implement the matrix inversion in the three-dimensional implicit
scheme by an iterative method. In this case, the helix inversion may
serve as a powerful preconditioner, providing an immediate answer in
constant velocity layers and speeding up the convergence in the case
of velocity variations. To see why this might be true, one can write
the variable-coefficient matrix $\tilde{\mathbf{A}}$ in the form
\begin{equation}
  \label{eqn:golub}
  \tilde{\mathbf{A}} = \mathbf{B} + \mathbf{D}\;,
\end{equation}
where matrix $\mathbf{B}$ corresponds to some constant average velocity,
and $\mathbf{D}$ is the matrix of velocity perturbations. The system of
linear equations that we need to solve is then