\lefthead{Fomel}\righthead{Velocity continuation}\footer{SEP--92}\title{Velocity continuation and the anatomy of \\residual prestack time migration}\email{sergey@sep.stanford.edu}\author{Sergey Fomel}\maketitle\begin{abstract}   Velocity continuation is an imaginary continuous process of  seismic image transformation in the post-migration domain. It generalizes  the concepts of residual and cascaded migrations. Understanding the laws  of velocity continuation is crucially important for a successful application  of time migration velocity analysis.  These laws predict the changes  in the geometry and intensity of reflection events on migrated images with  the change of the migration velocity.  In this paper, I derive kinematic  and dynamic laws for the case of prestack residual migration from simple  geometric principles. The main theoretical result is a decomposition  of prestack velocity continuation into three different components  corresponding to residual normal moveout, residual dip moveout, and residual  zero-offset migration. I analyze the contribution and properties of each of  the three components separately. This theory forms the basis for  constructing efficient finite-difference and spectral algorithms for  time migration velocity analysis.\end{abstract}\section{Introduction}%%%%%%%%%%%%%%%%%%\begin{comment}Migration velocity analysis is a routine part of prestack timemigration applications. It serves both as a tool for velocityestimation \cite[]{FBR08-06-02240234} and as a tool for optimal stackingof migrated seismic sections and modeling zero-offset data for depthmigration \cite[]{GEO62-02-05680576}. In the most common form, migrationvelocity analysis amounts to residual moveout correction on CRP(common reflection point) gathers. However, in the case of dippingreflectors, this correction does not provide optimal focusing ofreflection energy, since it does not account for lateral movement ofreflectors caused by the change in migration velocity. In other words,different points on a stacking hyperbola in a CRP gather cancorrespond to different reflection points at the actual reflector. Thesituation is similar to that of the conventional NMO velocityanalysis, where the reflection point dispersal problem is usuallyovercome with the help of DMO \cite[]{FBR04-07-00070024,dmo}. Ananalogous correction is required for optimal focusing in thepost-migration domain. In this paper, I propose and test velocitycontinuation as a method of migration velocity analysis. The methodenhances the conventional residual moveout correction by taking intoaccount lateral movements of migrated reflection events.\end{comment}The conventional approach to seismic migration theory\cite[]{Claerbout.blackwell.85,Berkhout.mig.14A.1985} employs thedownward continuation concept. According to this concept, migrationextrapolates upgoing reflected waves, recorded on the surface, to theplace of their reflection to form an image of subsurfacestructures.%In order to understand the concept of velocity continuation, we need%to look at the fundamentals of seismic time migration. Post-stack time migration possesses peculiar properties, which canlead to a different viewpoint on migration.  One of the mostinteresting properties is an ability to decompose the time migrationprocedure into a cascade of two or more migrations with smallermigration velocities. This remarkable property is described by\cite{GEO50-01-01100126} as {\em residual migration}.\cite{GEO52-05-06180643} generalized the method of residualmigration to one of {\em cascaded migration.} Cascadingfinite-difference migrations overcomes the dip limitations ofconventional finite-difference algorithms \cite[]{GEO52-05-06180643};cascading Stolt-type {\em f-k} migrations expands their range ofvalidity to the case of a vertically varying velocity\cite[]{GEO53-07-08810893}. Further theoretical generalization setsthe number of migrations in a cascade to infinity, making each step inthe velocity space infinitesimally small. This leads to a partialdifferential equation in the time-midpoint-velocity space, discoveredby \cite{Claerbout.sep.48.79}. Claerbout's equation describes theprocess of {\em velocity continuation,} which fills the velocity spacein the same manner as a set of constant-velocity migrations. Slicingin the migration velocity space can serve as a method of velocityanalysis for migration with nonconstant velocity\cite[]{shurtleff,SEG-1984-S1.8,Fowler.sepphd.58, GEO57-01-00510059}.\new{The concept of velocity continuation was introduced in the earlierpublications} \cite[]{me,SEG-1997-1762}.  \cite{hubral} and\cite{GEO62-02-05890597} \new{use the term \emph{image waves} to describe asimilar idea.} \cite{adler} \new{generalizes it to the case of variablebackground velocities under the name \emph{Kirchhoff image propagation}.  Theimportance of this concept lies in its ability to predict changes in thegeometry and intensity of reflection events on seismic images with the changeof migration velocity. While conventional approaches to migration velocityanalysis methods take into account only vertical movement of reflectors}\cite[]{FBR08-06-02240234,GEO60-01-01420153}, \new{velocity continuation attempts to describe bothvertical and lateral movements, thus providing for optimal focusing invelocity analysis applications} \cite[]{SEG-2001-11071110,second}. In this paper, I describe the velocity continuation theory for the case ofprestack time migration, connecting it with the theory of prestack residualmigration \cite[]{Al-Yahya.sep.50.219,Etgen.sepphd.68,GEO61-02-06050607}. Byexploiting the mathematical theory of characteristics, a simplified kinematicderivation of the velocity continuation equation leads to a differentialequation with correct dynamic properties. %In the post-stack case, the%solution of the boundary-value problem, associated with this equation,%coincides precisely with the operators of Kirchoff migration, traditionally%derived from a completely different prospective.  In practice, one canaccomplish dynamic velocity continuation by integral, finite-difference, orspectral methods.\begin{comment}For practicalapplications, I chose the Fourier spectral method. The method has itslimitations \cite[]{Fomel.sep.97.sergey2}, but looks optimal in terms ofthe accuracy versus efficiency trade-off. Applying velocity continuation to migration velocity analysis involvesthe following steps: \begin{enumerate}\item prestack common-offset (and common-azimuth) migration - to  generate the initial data for continuation,\item velocity continuation with stacking 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 final part of this paper includes a demonstration of all threesteps on a simple two-dimensional dataset.\end{comment}The accompanying paper \cite[]{second} introduces one of the possiblenumerical implementations and demonstrates its application on a \new{field}data example.\new{The paper is organized into two main sections. First, I derive the kinematicsof velocity continuation from the first geometric principles. I identify threedistinctive terms, corresponding to zero-offset residual migration, residualnormal moveout, and residual dip moveout. Each term is analyzed separately toderive an analytical prediction for the changes in the geometry of traveltime curves (reflection events on migrated images) with the change ofmigration velocity. Second, the dynamic behavior of seismic images isdescribed with the help of partial differential equations and theirsolutions. Reconstruction of the dynamical counterparts for kinematicequations is not unique. However, I show that, with an appropriate selection ofadditional terms, the image waves corresponding to the velocity continuationprocess have the correct dynamic behavior. In particular, a special boundary value problem with the zero-offset velocity continuation equationproduces the solution identical to the conventional Kirchoff time migration.}\begin{comment}It is important to note that although the velocity continuation resultcould be achieved in principle by using prestack residual migration inKirchhoff \cite[]{Etgen.sepphd.68} or Stolt \cite[]{GEO61-02-06050607}formulation, the first is evidently inferior in efficiency, and thesecond is not convenient for velocity analysis across differentoffsets, because it mixes them in the Fourier domain.\end{comment}\section{KINEMATICS OF VELOCITY CONTINUATION}From the kinematic point of view, it is convenient to describe thereflector as a locally smooth surface $z = z(x)$, where $z$ is thedepth, and $x$ is the point on the surface ($x$ is a two-dimensionalvector in the 3-D problem). The image of the reflector obtained aftera common-offset prestack migration with a half-offset $h$ and aconstant velocity $v$ is the surface $z = z(x;h,v)$. Appendix Aprovides the derivations of the partial differential equationdescribing the image surface in the depth-midpoint-offset-velocityspace. The purpose of this section is to discuss the laws of kinematictransformations implied by the velocity continuation equation. Laterin this paper, I obtain dynamic analogs of the kinematicrelationships in order to describe the continuation of migratedsections in the velocity space.The kinematic equation for prestack velocity continuation, derived inAppendix A, takes the following form:\begin{equation}{{\partial \tau} \over {\partial v}} = v\,\tau\,\left({{\partial \tau} \over {\partial x}}\right)^2 +{{h^2} \over {v^3\,\tau}}\,- \frac{h^2 v}{\tau}\,\left({{\partial \tau} \over {\partial x}}\right)^2\,\left({{\partial \tau} \over {\partial h}}\right)^2\;.\label{eq:eikonal} \end{equation}Here $\tau$ denotes the one-way vertical traveltime $\left(\tau = {z\over v}\right)$. The right-hand side of equation (\ref{eq:eikonal})consists of three distinctive terms. Each has its own geophysical meaning. The first term is the only one remainingwhen the half-offset $h$ equals zero. This term corresponds to the procedure of{\em zero-offset residual migration}.  Setting the traveltime dip tozero eliminates the first and third terms, leaving the second,dip-independent one. One can associate the second term with the process of{\em residual normal moveout}. The third term is both dip- and offset-dependent. The process that it describes is {\em residual dipmoveout}. It is convenient to analyze each of the three processesseparately, evaluating their contributions to the cumulative processof prestack velocity continuation.\subsection{Kinematics of Zero-Offset Velocity Continuation}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%The kinematic equation for zero-offset velocity continuation is\begin{equation}{{\partial \tau} \over {\partial v}} = v\,\tau\,\left({{\partial \tau} \over {\partial x}}\right)^2\;.\label{eq:POMeikonal} \end{equation}The typical boundary-value problem associated with it is to find thetraveltime surface $\tau_2(x_2)$ for a constant velocity $v_2$, given thetraveltime surface $\tau_1(x_1)$ at some other velocity $v_1$. Both surfacescorrespond to the reflector images obtained by time migrations with thespecified velocities.  When the migration velocity approaches zero, post-stacktime migration approaches the identity operator. Therefore, the case of $v_1 =0$ corresponds kinematically to the zero-offset (post-stack) migration, andthe case of $v_2 = 0$ corresponds to the zero-offset modeling (demigration).\new{The variable $x$ in equation~(\ref{eq:POMeikonal}) describes both the surfacemidpoint coordinate and the subsurface image coordinate. One of them iscontinuously transformed into the other in the velocity continuation process.}The appropriate mathematical method of solving the kinematicproblem posed above is the method of characteristics \cite[]{kurant2}. Thecharacteristics of equation (\ref{eq:POMeikonal}) are the trajectoriesfollowed by individual points of the reflector image in the velocitycontinuation process. These trajectories are called {\em velocity rays} \cite[]{me,them,adler}. Velocity rays are defined by the system of ordinarydifferential equations derived from (\ref{eq:POMeikonal}) according to the\new{Hamilton-Jacobi theory}:\begin{eqnarray}{{{dx} \over {dv}} = - 2\,v\,\tau\,\tau_x} & , &{{{d\tau} \over {dv}} = - \tau_v}\;,\label{eq:velray1} \\{{{d\tau_x} \over {dv}} = v\,\tau_x^3} & , &{{{d\tau_v} \over {dv}} = \left(\tau + v\,\tau_v\right)\,\tau_x^2}\;,\label{eq:velray2} \end{eqnarray}\new{where $\tau_x$ and $\tau_v$ are the phase-space parameters}.An additional constraint for $\tau_x$ and $\tau_v$follows from equation (\ref{eq:POMeikonal}), rewritten in the form\begin{equation}\tau_v = v\,\tau\,\tau_x^2\;. \label{eq:equiveikonal} \end{equation}%One can easily solve the system of equations (\ref{eq:velray1}) and%(\ref{eq:velray2}) by the classic mathematical methods for ordinary%differential equations. The general solution of the system ofequations~(\ref{eq:velray1}-\ref{eq:velray2}) takes theparametric form\begin{eqnarray}x(v) & = & A - C v^2\;,\quad\tau^2(v) = B - C^2\,v^2\;,\label{eq:velrayg1} \\ \tau_x(v) & = & {C \over {\tau(v)}}\;,\quad\tau_v(v) = {{C^2\,v} \over {\tau(v)}}\;,\label{eq:velrayg2} \end{eqnarray}where $A$, $B$, and $C$ are constant along each individual velocityray. These three constants are determined from the boundary conditionsas\begin{equation}A = x_1 + v_1^2\,\tau_1\,{{\partial \tau_1} \over {\partial x_1}} = x_0\;,\label{eq:a} \end{equation}\begin{equation}B = \tau_1^2\,\left(1 + v_1^2\,\left({{\partial \tau_1} \over {\partial x_1}}\right)^2\right) = \tau_0^2\;,\label{eq:b} \end{equation}\begin{equation}C = \tau_1\,{{\partial \tau_1} \over {\partial x_1}} = \tau_0\,{{\partial \tau_0} \over {\partial x_0}}\;,\label{eq:c} \end{equation}where $\tau_0$ and $x_0$ correspond to the zero velocity (unmigratedsection), while $\tau_1$ and $x_1$ correspond to the velocity $v_1$.% Equations (\ref{eq:a}), (\ref{eq:b}), and (\ref{eq:c}) have a clear%geometric meaning illustrated in Figure \ref{fig:vlczor}. Thesimple relationship between the midpoint derivative of the verticaltraveltime and the local dip angle $\alpha$ (appendix A),\begin{equation}{{\partial \tau} \over {\partial x}} = {{\tan{\alpha}} \over v}\;,\label{eq:dtaudx} \end{equation}shows that equations (\ref{eq:a}) and (\ref{eq:b}) are precisely equivalentto the evident geometric relationships (Figure~\ref{fig:vlczor})\begin{equation}x_1 + v_1\,\tau_1\,\tan{\alpha} = x_0\;,\;{\tau_1 \over {\cos{\alpha}}} = \tau_0\;.\label{eq:evident}\end{equation}Equation (\ref{eq:c}) states that the points on a velocity ray correspondto a single reflection point, constrained by the values of $\tau_1$,$v_1$, and $\alpha$.  As follows from equations (\ref{eq:velrayg1}), theprojection of a velocity ray to the time-midpoint plane has theparabolic shape $x(\tau) = A + (\tau^2 - B) / C$, which has beennoticed by \cite{GEO46-05-07170733}. On the depth-midpoint plane, thevelocity rays have the circular shape $z^2(x) = (A - x)\,B / C - (A -x)^2$, described by \cite{them} as ``Thales circles.''\inputdir{XFig}\sideplot{vlczor}{width=0.9\textwidth}{Zero-offset reflection in a  constant velocity medium (a scheme).}For an example of kinematic continuation by velocity rays, let usconsider the case of a point diffractor. If the diffractor location inthe subsurface is the point ${x_d,z_d}$, then the reflection traveltime atzero offset is defined from Pythagoras's theorem as the hyperboliccurve\begin{equation}\tau_0(x_0) = {{\sqrt{z_d^2 + (x_0 - x_d)^2}} \over v_d}\;,\label{eq:dift0}\end{equation}where $v_d$ is half of the actual velocity. Applying equations(\ref{eq:velrayg1}) produces the following mathematical expressionsfor the velocity rays:\begin{eqnarray}x(v) & = & x_d\,{v^2 \over v_d^2} + x_0\,\left(1 -  {v^2 \over v_d^2}\right)\;,\label{eq:difrayg1} \\ \tau^2(v) & = & \tau_d^2 + {{(x_0 - x_d)^2} \over v_d^2}\,\left(1 -  {v^2 \over v_d^2}\right)\;,\label{eq:difrayg2} \end{eqnarray}where $\tau_d = {z_d \over v_d}$.Eliminating $x_0$ from the system of equations (\ref{eq:difrayg1}) and(\ref{eq:difrayg2}) leads to the expression for the velocity continuation``wavefront'': \begin{equation}\tau(x)=\sqrt{\tau_d^2 + {{(x - x_d)^2} \over {v_d^2 - v^2}}}\;.\label{eq:diffront}\end{equation}For the case of a point diffractor, the wavefront corresponds preciselyto the summation path of the residual migration operator\cite[]{GEO50-01-01100126}. It has a hyperbolic shape when $v_d > v$(undermigration) and an elliptic shape when $v_d < v$(overmigration). The wavefront collapses to a point when the velocity$v$ approaches the actual effective velocity $v_d$. At zerovelocity, $v=0$, the wavefront takes the familiar form of the post-stack migrationhyperbolic summation path. The form of the velocity rays and wavefrontsis illustrated in the left plot of Figure \ref{fig:vlcvrs}.\inputdir{Math}