paper.tex

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

\def\figdir{./Fig}
\def\cakedir{.}

\lefthead{Fomel}
\righthead{Velocity continuation}
\footer{SEP--92}

\title{Time migration velocity analysis by velocity continuation}

\email{[e-mail]}
\author{Sergey Fomel}
\maketitle

\begin{abstract}
Time migration velocity analysis can be performed by velocity
  continuation, an incremental process that transforms
  migrated seismic sections according to changes in the migration velocity.
  Velocity continuation enhances residual \new{normal moveout} correction by
  properly taking into account both vertical and lateral movements of
  \new{events on seismic images}. Finite-difference and spectral algorithms
  provide efficient practical implementations for velocity continuation.
  Synthetic and field data examples demonstrate the performance of the method
  and confirm theoretical expectations.
\end{abstract}

\section{Introduction}
%%%%%%%%%%%%%%%%%%
Migration velocity analysis is a routine part of prestack time
migration applications. It serves both as a tool for velocity
estimation \cite[]{FBR08-06-02240234} and as a tool for optimal
stacking of migrated seismic sections prior to modeling zero-offset
data for depth migration \cite[]{GEO62-02-05680576}. In the most
common form, migration velocity analysis amounts to residual moveout
correction on CIP (common image point) gathers. However, in the case
of dipping reflectors, this correction does not provide optimal
focusing of reflection energy, since it does not account for lateral
movement of reflectors caused by the change in migration velocity. In
other words, different points on a stacking hyperbola in a CIP gather
can correspond to different reflection points at the actual reflector.
The situation is similar to that of the conventional \new{normal
  moveout} (NMO) velocity analysis, where the reflection point
dispersal problem is usually overcome with the help of \new{dip
  moveout} \cite[]{FBR04-07-00070024,dmo}. An analogous correction is
required for optimal focusing in the post-migration domain. In this
paper, I propose and test velocity continuation as a method of
migration velocity analysis. The method enhances the conventional
residual moveout correction by taking into account lateral movements
of migrated reflection events.

Velocity continuation is a process of transforming time migrated images
according to the changes in migration velocity. This process has wave-like
properties, which have been described in earlier papers
\cite[]{me,SEG-1997-1762,first}.  \cite{hubral} and
\cite{GEO62-02-05890597} use the term \emph{image waves} to describe a
similar concept. \cite{adlerSEG,adler} generalizes the velocity
continuation approach for the case of variable background velocities, using
the term \emph{Kirchhoff image propagation}. \new{Although the velocity
continuation concept is tailored for time migration, it finds important
applications in depth migration velocity analysis by recursive methods} 
\cite[]{SEG-1999-17231726,SEG-2000-08740877}.

\par
Applying velocity continuation to migration velocity analysis involves
the following steps: 
\begin{enumerate}
\item prestack common-offset (and common-azimuth) migration - to
  generate the initial data for continuation,
\item velocity continuation with stacking and semblance analysis across
  different offsets - to transform the offset data dimension into the velocity
  dimension,
\item picking the optimal velocity and slicing through the migrated
  data volume - to generate an optimally focused image.
\end{enumerate}
The first step transforms the data to the image space. The regularity of this
space can be exploited for devising efficient algorithms for the next two
steps. The idea of slicing through the velocity space goes back to the work of
\cite{shurtleff}, \cite{SEG-1984-S1.8,Fowler.sepphd.58}, and
\cite{GEO57-01-00510059}. While the previous slicing methods constructed
the velocity space by repeated migration with different velocities, velocity
continuation navigates directly in the migration velocity space without
returning to the original data. This leads to both more efficient algorithms
and a better understanding of the theoretical continuation properties
\cite[]{first}.

In this paper, I demonstrate all three steps, using both synthetic data and a
North Sea dataset. I introduce and exemplify two methods for the efficient
practical implementation of velocity continuation: the finite-difference
method and the Fourier spectral method. The Fourier method is recommended as
optimal in terms of the accuracy versus efficiency trade-off. Although all the
examples in this paper are two-dimensional, the method easily extends to 3-D
under the assumption of common-azimuth geometry (one oriented offset). More
investigation may be required to extend the method to the multi-azimuth case.

It is also important to note that although the velocity continuation result
could be achieved in principle by using prestack residual migration in
Kirchhoff \cite[]{Etgen.sepphd.68} or frequency-wavenumber
\cite[]{GEO61-02-06050607} formulation, the first is inferior in efficiency, and
the second is not convenient for the conventional velocity analysis across
different offsets, because it mixes them in the Fourier domain
\cite[]{SEG-2000-09920995}.  \new{Fourier-domain angle-gather analysis}
\cite[]{SEG-2001-02960299,sandf} \new{opens new possibilities for the future
development of the Fourier-domain velocity continuation. New insights into the
possibility of extending the method to depth migration can follow from the
work of} \cite{adler}.

\section{Numerical velocity continuation in the post-stack domain}
\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[]{first}, 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
propagation 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'' (demigration); the latter case, to seismic
migration.
\par
Mathematically, equations (\ref{PDE}) and (\ref{BC}) define a
Goursat-type problem \cite[]{kurant2}. Its analytical solution can be
constructed by a variation of the Riemann method in the form of an
integral operator \cite[]{me,first}:
\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[]{GEO43-01-00230048}, regarded as the most efficient constant-velocity
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 approaches are discussed in the next two
subsections.
\subsection{Finite-difference approach}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The differential equation~(\ref{PDE}) has a mathematical form analogous to
that of the 15-degree wave extrapolation equation
\cite[]{Claerbout.blackwell.76}. Its finite-difference implementation, first
described by Claerbout \shortcite{Claerbout.sep.48.79} and Li
\shortcite{Li.sep.48.85}, is also analogous to that of the 15-degree equation,
except for the variable coefficients. One can write the implicit
unconditionally stable finite-difference scheme for the velocity continuation
equation in the form
\begin{equation}
({\bf I} + a_{j+1}^{i+1}\,{\bf T})\,{\bf P}_{j+1}^{i+1} - 
({\bf I} - a_{j}^{i+1}\,{\bf T})\,{\bf P}_{j}^{i+1} -
({\bf I} - a_{j+1}^{i}\,{\bf T})\,{\bf P}_{j+1}^{i} + 
({\bf I} + a_{j}^{i}\,{\bf T})\,{\bf P}_{j}^{i} = 0\;, 
\label{eqn:fds} 
\end{equation}
where index $i$ corresponds to the time dimension, index $j$
corresponds to the velocity dimension, ${\bf P}$ is a vector along the
midpoint direction, ${\bf I}$ is the identity matrix, ${\bf T}$
represents the finite-difference approximation to the 
second-derivative operator in midpoint,
and $a_{j}^{i} = v_j\,t_i\,{\Delta v\,\Delta t}$.

\inputdir{XFig}

In the two-dimensional case, equation \ref{eqn:fds} reduces to a
tridiagonal system of linear equations, which can be easily inverted.
In 3-D, a straightforward extension can be obtained by using either
directional splitting or helical schemes \cite[]{SEG-1998-1124}. The
direction of stable propagation is either forward in velocity and
backward in time or backward in velocity and forward in time as shown
in Figure \ref{fig:vlcfds}.

\sideplot{vlcfds}{height=2in}{Finite-difference scheme for the
velocity continuation equation. A stable propagation is either forward
in velocity and backward in time (a) or backward in velocity
and forward in time (b).}

\inputdir{vcmod}

In order to test the performance of the finite-difference velocity
continuation method, I use a simple synthetic model from
\cite{Claerbout.bei.95}. The reflectivity model is shown in
Figure \ref{fig:mod}. It contains several features that challenge the
migration performance: dipping beds, unconformity, syncline,
anticline, and fault. \new{The velocity is taken to be constant 
$v=1.5\,\mbox{km/s}$.}

\sideplot{mod}{height=2.5in}{Synthetic model for testing
finite-difference migration by velocity continuation.}

\new{Figures~\ref{fig:vlckaa}--\ref{fig:vlcspv}} compare invertability of different
migration methods. In all cases, constant-velocity modeling (demigration) 
%by the adjoint operator \cite[]{Claerbout.blackwell.92}
was followed by migration with the correct velocity.  Figures
\ref{fig:vlckaa} and \ref{fig:vlcsto} show the results of modeling and
migration with the Kirchhoff \cite[]{GEO43-01-00490076} and $f$-$k$
\cite[]{GEO43-01-00230048} methods, respectively. These figures should
be compared with Figure \ref{fig:vlcvel}, showing the analogous result
of the finite-difference velocity continuation.  The comparison
reveals a remarkable invertability of velocity continuation, which
reconstructs accurately the main features and frequency content of the
model. Since the forward operators were different for different
migrations, this comparison did not test the migration properties
themselves. For such a test, I compare the results of the Kirchhoff
and velocity-continuation migrations after Stolt modeling.  The result
of velocity continuation, shown in Figure \ref{fig:vlcspv}, is
noticeably more accurate than that of the Kirchhoff method.

%\activeplot{vlckir}{width=6in,height=2.5in}{ER}{Result of modeling and
%migration with the fast Kirchhoff method. Left plot shows the
%reconstructed model. Right plot compares the average amplitude
%spectrum of the true model with that of the reconstructed image.}

\plot{vlckaa}{width=6in}{Result of modeling and migration with the
Kirchhoff method. Top left plot shows the reconstructed model. Top
right plot compares the average amplitude spectrum of the true model
with that of the reconstructed image.  Bottom left is the
reconstruction error. Bottom right is the absolute error in the
spectrum.}

%\activeplot{vlcpha}{width=6in,height=2.5in}{ER}{Result of modeling and
%migration with the phase-shift method. Left plot shows the
%reconstructed model. Right plot compares the average amplitude
%spectrum of the true model with that of the reconstructed image.}

\plot{vlcsto}{width=6in}{Result of modeling and migration with the
  Stolt method. Top left plot shows the reconstructed model.  Top
  right plot compares the average amplitude spectrum of the true model