paper.tex

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

  with that of the reconstructed image. Bottom left is the
  reconstruction error. Bottom right is the absolute error in the
  spectrum.}

\plot{vlcvel}{width=6in}{Result of modeling and
  migration with the finite-difference velocity continuation. 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.  }

\plot{vlcspv}{width=6in,height=5in}{Top:
modeling with Stolt method, migration with the Kirchhoff
method. Bottom: modeling with Stolt method, migration with the
finite-difference velocity continuation. Left plots show the
reconstructed models. Right plots show the reconstruction errors.}

These tests confirm that finite-difference velocity continuation is an
attractive migration method. It possesses remarkable invertability properties,
which may be useful in applications that \new{require inversion. While the
  traditional migration methods transform the data between two completely
  different domains (data-space and image-space), velocity continuation
  accomplishes the same trans\-for\-ma\-ti\-on by propagating the data in the
  extended domain along the velocity direction. Inverse propagation restores
  the original data.} According to \cite{Li.sep.48.85}, the computational
speed of this method compares favorably with that of Stolt migration. The
advantage is apparent for cascaded migration or migration with multiple
velocity models. In these cases, the cost of Stolt migration increases in
direct proportion to the number of velocity models, while the cost of velocity
continuation stays the same.

\par
\subsection{Fourier approach}
\par
The change of variable $\sigma = t^2$ transforms
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 do not 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\,{{d \hat{P}} \over {d v}} -
  v\,k^2\,\hat{P} = 0\;,
\end{equation}
where the ``frequency'' variable $\Omega$ corresponds to the stretched time
coordinate $\sigma$, and $k$ is the wavenumber in $x$:
\begin{equation}
\label{FFT}
\hat{P}(k,\Omega,v) = \int\!\!\int P(x,t,v)\,
e^{-i\,\Omega\,t^2 + i\,k\,x} d x\,dt
\end{equation}
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 leads to a very simple algorithm for the numerical velocity
continuation. The algorithm consists of the following steps:
\begin{enumerate}
\item \new{Input the zero-offset (post-stack) data migrated with velocity $v_0$ (or
  unmigrated if $v_0=0$).}
\item Transform the input from a regular grid in $t$ to a regular grid
  in $\sigma$.
\item Apply Fast Fourier Transform (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 generated from the model in Figure~\ref{fig:mod}
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 region, where the data
spectrum is zero. It 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.}

\new{The described algorithm}
 is very attractive from the practical point of view
because of its efficiency (based on the FFT algorithm). The operation count is
roughly the same as in the Stolt migration implemented with equation
(\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). The velocity continuation algorithm can be more efficient
than the Stolt method because of the simpler structure of the innermost loop
\new{(step 4 in the algorithm)}. 

\section{Numerical velocity continuation in the prestack domain}

To generalize the algorithm of the previous section to the prestack case, it is
first necessary to include the residual NMO term \cite[]{first}.  Residual normal
moveout can be formulated with the help of the differential equation:
\begin{equation}
{{\partial P} \over {\partial v}} + 
{{h^2} \over {v^3\,t}}\,{{\partial P} \over {\partial t}} = 0\;,
\label{eqn:ResNMOdyn} 
\end{equation}
where $h$ stands for the half-offset.  The analytical solution of
equation (\ref{eqn:ResNMOdyn}) has the form of the residual NMO
operator:
\begin{equation}
  P(t,h,v) = P_0\left(\sqrt{t^2 + h^2\,
    \left(\frac{1}{v_0^2} - \frac{1}{v^2}\right)},h\right)\;.
\label{eqn:ResNMOsol} 
\end{equation}
After transforming to the squared time $\sigma = t^2$ and the
corresponding Fourier frequency $\Omega$, equation
(\ref{eqn:ResNMOdyn}) takes the form of the ordinary differential
equation
\begin{equation}
  \frac{d \hat{P}}{d v} + 
  i \Omega\,\frac{2\,h^2}{v^3}\,\hat{P} = 0
\label{eqn:FFTResNMO} 
\end{equation}
with the analytical frequency-domain phase-shift solution
\begin{equation}
  \hat{P} (\Omega, h, v) = \hat{P_0} (\Omega,h) e^{i\,\Omega\,h^2\,
    \left(\frac{1}{v_0^2} - \frac{1}{v^2}\right)}\;.
\label{eqn:ResNMOshift} 
\end{equation}
To obtain a Fourier-domain prestack velocity continuation algorithm, one just
needs to combine the phase-shift operators in equations (\ref{ODEsol}) and
(\ref{eqn:ResNMOshift}) and to include stacking across different offsets. 
\new{The
exact velocity continuation theory also includes the residual DMO term}
\cite[]{first}, \new{which has a second-order effect, pronounced only at small
depths. It is neglected here for simplicity.}  The algorithm takes the
following form:
\begin{enumerate}
  \item Input a set of common-offset images, migrated with velocity $v_0$.
  \item Transform the time axis $t$ to the squared time coordinate:
    $\sigma=t^2$.
  \item Apply a fast Fourier transform (FFT) on both the squared time
    and the midpoint axis. The squared time $\sigma$ transforms to the
    frequency $\Omega$, and the midpoint coordinate $x$ transforms to
    the wavenumber $k$. 
  \item Apply a phase-shift operator to transform to different velocities $v$:
    \begin{equation}
      \label{eqn:zophase2}
      \hat{P}(\Omega,k,v) = \sum_{h} \hat{P}_0 (\Omega,k,h)\,
      e^{i\,\frac{k^2\left(v_0^2 - v^2\right)}{4\,\Omega} +
        i\,\Omega\,h^2\, \left(\frac{1}{v_0^2} -
          \frac{1}{v^2}\right)}\;.
    \end{equation}
    To save memory, the continuation step is immediately followed
    by stacking. \new{For velocity analysis purposes, a semblance measure}
    \cite[]{GEO36-03-04820497} \new{is computed in addition to the simple stack
    analogously to the standard practice of stacking velocity analysis.}
    
  \new{Implementing the residual moveout correction in the Fourier
    domain allows one to package it conveniently with the phase-shift
    operator without the need to transform the continuation result
    back to the time domain. The offset dimension in
    equation~(\ref{eqn:zophase2}) is replaced by the velocity
    dimension similarly to the velocity transform of the conventional
    stacking velocity analysis} \cite[]{yilmaz}.

  \item Apply an inverse FFT to transform from $\Omega$ and $k$ to
    $\sigma$ and $x$.
  \item Apply an inverse time stretch to transform from $\sigma$ to $t$.
  \end{enumerate}  
  One can design similar algorithms by using the finite difference method.
  \new{Although the finite-difference approach offers a faster continuation speed,
  the spectral algorithm has a higher accuracy while maintaining an acceptable
  cost.}

\inputdir{vcspk}

%The complete theory of prestack velocity continuation also requires a
%residual DMO operator
%\cite[]{Etgen.sepphd.68,Fomel.sep.92.159,Fomel.segab.97}. However, the
%difficulty of implementing this operator is not fully compensated by
%its contribution to the full velocity continuation. For simplicity, I
%decided not to include residual DMO in the current implementation.
\par  
Figure~\ref{fig:velimp} shows impulse responses of prestack velocity
continuation. The input for producing this figure was a time-migrated
constant-offset section, corresponding to an offset of 1 km and a constant
migration velocity of 1 km/s. In full accordance with the theory \cite[]{first},
three spikes in the input section transformed into shifted ellipsoids after
continuation to a higher velocity and into shifted hyperbolas after
continuation to a smaller velocity. \new{
Padding of the time axis helps to avoid
the wrap-around artifacts of the Fourier method. Alternatively, one could use
the artifact-free but more expensive Chebyshev spectral method}
\cite[]{Fomel.sep.97.sergey2}.

\plot{velimp}{width=6in}{Impulse responses of prestack velocity
  continuation. Left plot: continuation from 1 km/s to 1.5 km/s.
  Right plot: continuation from 1 km/s to 0.7 km/s. Both plots
  correspond to the offset of 1 km.}
\par

Velocity continuation creates a time-midpoint-velocity cube (four-dimensional
for 3-D data), which is convenient for picking imaging velocities in the same
way as the result of common-midpoint or common-reflection-point velocity
analysis. The important difference is that velocity continuation provides an
optimal focusing of the reflection energy by properly taking into account both
vertical and lateral movements of reflector images with changing migration
velocity. \new{An experimental evidence for this conclusion is provided in
  the examples section of this paper.}

\begin{comment}
Figure \ref{fig:consmb} compares velocity spectra
(semblance panels) at a CIP location of about 11.5 km after residual
NMO and after prestack velocity continuation. Although the overall
difference between the two panels is small, the velocity continuation
panel shows a noticeably better focusing, especially in the region of
conflicting dips between 1 and 2 seconds. 
\end{comment}
The next subsection discusses
the velocity picking step in more detail.

%\plot{consmb}{width=6in}{Velocity spectra around 11.5 km CRP after
%  residual NMO (left) and after prestack velocity continuation
%  (right). The right plot shows improved focusing in the region between
%  1 and 2 seconds.}

\subsection{Velocity picking and slicing}

After the velocity continuation process has created a time-midpoint-velocity
cube, one can pick the best focusing velocity from that cube and create an
optimally focused image by slicing through the cube. \new{This step is common in
other methods that involve velocity slicing}
\cite[]{shurtleff,SEG-1984-S1.8,GEO57-01-00510059}. \new{The algorithm described
below has been also adopted by} \cite{SEG-2000-09920995} \new{for velocity
analysis in wave-equation migration.}

A simple automatic velocity picking algorithm follows from solving the
following regularized least-squares system: