paper.tex

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

\lefthead{Fomel}
\righthead{Spectral velocity continuation}
\footer{SEP--97}

\title{Velocity continuation by spectral methods}

\email{[e-mail]} 
%\keywords{imaging, velocity continuation, spectral methods, Fourier,
%  Chebyshev}

\author{Sergey Fomel}

\maketitle

\begin{abstract}
I apply Fourier and Chebyshev spectral methods to derive accurate
  and efficient algorithms for velocity continuation.  As expected,
  the accuracy of the spectral methods is noticeably superior to that
  of the finite-difference approach. Both methods apply a
  transformation of the time axis to squared time. The Chebyshev
  method is slightly less efficient than the Fourier method, but has
  less problems with the time transformation and also handles
  accurately the non-periodic boundary conditions.
\end{abstract}

\section{Introduction}

\inputdir{impls}

\par
In a recent work \cite[]{me,Fomel.sep.92.159}, I
introduced the process of \emph{velocity continuation} to describe a
continuous transformation of seismic time-migrated images with a
change of the migration velocity. Velocity continuation generalizes
the ideas of residual migration
\cite[]{GEO50-01-01100126,Etgen.sepphd.68} and cascaded migrations
\cite[]{GEO52-05-06180643}. In the zero-offset (post-stack) case, the
velocity continuation process is governed by a partial differential
equation in midpoint, time, and velocity coordinates, first discovered
by \cite{Claerbout.sep.48.79}. \cite{hubral1} and
\cite{hubral2} describe this process in a broader context of
``image waves''. Generalizations are possible for the non-zero offset
(prestack) case \cite[]{Fomel.sep.92.159,SEG-1997-1762}.
\par
A numerical implementation of velocity continuation process provides
an efficient method of scanning the velocity dimension in the search
of an optimally focused image. The first implementations
\cite[]{Li.sep.48.85,Fomel.sep.92.159} used an analogy with Claerbout's
15-degree depth extrapolation equation to construct a
finite-difference scheme with an implicit unconditionally stable
advancement in velocity. \cite{Fomel.sep.95.sergey2} presented an
efficient three-dimensional generalization, applying the helix
transform \cite[]{Claerbout.sep.95.jon1}.
\par
\plot{fd-imp}{width=6in}{Impulse responses (Green's functions) of
  velocity continuation, computed by a second-order finite-difference
  method.  The left plots corresponds to continuation to a larger
  velocity ($+1$ km/sec); the right plot, smaller velocity, ($-1$
  km/sec).}
\par
A low-order finite-difference method is probably the most efficient
numerical approach to this method, requiring the least work per
velocity step. However, its accuracy is not optimal because of the
well-known numerical dispersion effect. Figure \ref{fig:fd-imp} shows
impulse responses of post-stack velocity continuation for three
impulses, computed by the second-order finite-difference method
\cite[]{Fomel.sep.92.159}. As expected from the residual migration
theory \cite[]{GEO50-01-01100126}, continuation to a higher velocity
(left plot) corresponds to migration with a residual velocity, and its
impulse responses have an elliptical shape. Continuation to a smaller
velocity (right plot in Figure \ref{fig:fd-imp}) corresponds to
demigration (modeling), and its impulse responses have a hyperbolic
shape. The dispersion artifacts are clearly visible in the figure.
\par
In this paper, I explore the possibility of implementing a numerical
velocity continuation by spectral methods. I adopted two different
methods, comparable in efficiency with finite differences.  The first
method is a direct application of the Fast Fourier Transform (FFT)
technique.  The second method transforms the time grid to Chebyshev
collocation points, which leads to an application of the
Chebyshev-$\tau$ method \cite[]{lanc,orszag,boyd}, combined with an
unconditionally stable implicit advancement in velocity.  Both methods
employ a transformation of the grid from time $t$ to the squared time
$\sigma = t^2$, which removes the dependence on $t$ from the
coefficients of the velocity continuation equation. Additionally, the
Fourier transform in the space (midpoint) variable $x$ takes care of
the spatial dependencies. This transform is a major source of
efficiency, because different wavenumber slices can be processed
independently on a parallel computer before transforming them back to
the physical space.
\par
\section{Problem formulation}
\par
The post-stack velocity continuation process is governed by a partial
differential equation in the domain, composed by the seismic image
coordinates (midpoint $x$ and vertical time $t$) and the additional
velocity coordinate $v$. Neglecting some amplitude-correcting terms
\cite[]{Fomel.sep.92.159}, the equation takes the form
\cite[]{Claerbout.sep.48.79}
\begin{equation}
  \frac{\partial^2 P}{\partial v\, \partial t} +
  v\,t\,\frac{\partial^2 P}{\partial x^2} = 0\;.
  \label{PDE}
\end{equation}
Equation (\ref{PDE}) is linear and belongs to the hyperbolic type. It
describes a wave-type process with the velocity $v$ acting as a
``time-like'' variable. Each constant-$v$ slice of the function
$P(x,t,v)$ corresponds to an image with the corresponding constant
velocity. The necessary boundary and initial conditions are 
\begin{equation}
  \label{BC}
  \left.P\right|_{t=T} = 0\;\quad \left.P\right|_{v=v_0} = P_0 (x,t)\;,
\end{equation}
where $v_0$ is the starting velocity, $T=0$ for continuation to a
smaller velocity and $T$ is the largest time on the image (completely
attenuated reflection energy) for continuation to a larger velocity.
The first case corresponds to ``modeling''; the latter case, to
seismic migration.
\par
Mathematically, equations (\ref{PDE}) and (\ref{BC}) define a
Goursat-type problem \cite[]{kurant}. Its analytical solution can be
constructed by a variation of the Riemann method in the form of an
integral operator \cite[]{me,Fomel.sep.92.159}:
\begin{equation}
P(t,x,v) = \frac{1}{(2\,\pi)^{m/2}}\,
{\int\, 
\frac{1}{(\sqrt{v^2-v_0^2}\,\rho)^{m/2}}\, 
\left(- \frac{\partial}{\partial t_0}\right)^{m/2}
P_0\left(\frac{\rho}{\sqrt{v^2-v_0^2}},x_0\right)\,dx_0}\;,
\label{kirch}
\end{equation}
where $\rho = \sqrt{(v^2-v_0^2)\,t^2 + (x - x_0)^2}$, $m=1$ in the 2-D
case, and $m=2$ in the 3-D case. In the case of continuation from zero
velocity $v_0=0$, operator (\ref{kirch}) is equivalent (up to the
amplitude weighting) to conventional Kirchoff time migration
\cite[]{GEO43-01-00490076}.  Similarly, in the frequency-wavenumber
domain, velocity continuation takes the form
\begin{equation}
  \label{stolt}
  \hat{P} (\omega,k,v) = \hat{P}_0 (\sqrt{\omega^2+k^2 (v^2-v_0^2)},k)\;,
\end{equation}
which is equivalent (up to scaling coefficients) to Stolt migration
\cite[]{GEO50-11-22192244}, regarded as the most efficient migration
method.
\par
If our task is to create many constant-velocity slices, there are
other ways to construct the solution of problem (\ref{PDE}-\ref{BC}).
Two alternative spectral approaches are discussed in the next two
sections.
\par
\section{Fourier approach}

\inputdir{sigvc}

\par
Introducing the change of variable $\sigma = t^2$, we can transform
equation (\ref{PDE}) to the form
\begin{equation}
  2\,\frac{\partial^2 P}{\partial v\, \partial \sigma} +
  v\,\frac{\partial^2 P}{\partial x^2} = 0\;,
  \label{PDE2}
\end{equation}
whose coefficients don't depend on the time variables.  Double Fourier
transform in $\sigma$ and $x$ further simplifies equation (\ref{PDE2})
to the ordinary differential equation
\begin{equation}
  \label{ODE}
  2\,i\Omega\,\frac{d^2 \hat{P}}{d v} -
  v\,k^2\,\hat{P} = 0\;,
\end{equation}
where the frequency $\Omega$ corresponds to the time coordinate
$\sigma$, and $k$ is the wavenumber in $x$. Equation (\ref{ODE}) has
an explicit analytical solution
\begin{equation}
  \label{ODEsol}
  \hat{P} (k,\Omega,v) = \hat{P}_0 (k,\Omega)\,
  e^{\frac{i k^2(v_0^2-v^2)}{4\Omega}}\;,
\end{equation}
which defines a very simple algorithm for the numerical velocity
continuation. The algorithms consists of the following steps:
\begin{enumerate}
\item Transform the input from a regular grid in $t$ to a regular grid
  in $\sigma$.
\item Apply FFT in $x$ and $\sigma$.
\item Multiply by the all-pass phase-shift filter $e^{\frac{i
      k^2(v_0^2-v^2)}{4\Omega}}$.
\item Inverse FFT in $x$ and $\sigma$.
\item Inverse transform to a regular grid in $t$.
\end{enumerate}

\plot{t2}{width=6in}{Synthetic seismic data before (left) and after
  (right) transformation to the $\sigma$ grid.}

Figure \ref{fig:t2} shows a simple synthetic model of seismic
reflection data from \cite[]{Claerbout.bei.95} before and after
transforming the grid, regularly spaced in $t$, to a grid, regular in
$\sigma$. The left plot of Figure \ref{fig:t2-fft} shows the Fourier
transform of the data. Except for the nearly vertical event, which
corresponds to a stack of parallel layers in the shallow part of the
data, the data frequency range is contained near the origin in the
$\Omega-k$ space.  The right plot of Figure \ref{fig:t2-fft} shows the
phase-shift filter for continuation from zero imaging velocity (which
corresponds to unprocessed data) to the velocity of 1 km/sec. The
rapidly oscillating part (small frequencies and large wavenumbers) is
exactly in the place, where the data spectrum is zero and corresponds
to physically impossible reflection events.

\plot{t2-fft}{width=6in}{Left: the real part of the data Fourier
  transform.  Right: the real part of the velocity continuation
  operator (continuation from 0 to 1 km/s) in the Fourier domain.}

Algorithm (\ref{ODEsol}) is very attractive from the practical point
of view because of its efficiency (based on the FFT algorithm). The
operations count is roughly the same as in Stolt migration
(\ref{stolt}): two forward and inverse FFTs and forward and inverse
grid transform with interpolation (one complex-number transform in the
case of Stolt migration). Algorithm (\ref{ODEsol}) can be even more
efficient than Stolt method because of the simpler structure of the
innermost loop.  However, its practical implementation faces two
difficult problems: artifacts of the $t^2$ grid transform and
wraparound artifacts
\par
\subsection{Improving the accuracy of the $t^2$ grid transform}
\par
The first problem is the loss of information in the transform to the
$t^2$ grid. As illustrated in Figure \ref{fig:t2}, the shallow part of
the data gets severely compressed in the $t^2$ grid. The amount of
compression can lead to inadequate sampling, and as a result, aliasing
artifacts in the frequency domain. Moreover, it can be difficult to
recover from the loss of information in the transformed domain when
transforming back into the original grid. A partial remedy for this
problem is to increase the grid size in the $t^2$ domain. The top plots
in Figure \ref{fig:fft-inv} show the result of back transformation to
the $t$ grid and the difference between this result and the original
model (plotted on the same scale). We can see a noticeable loss of
information in the upper (shallow) part of the data, caused by