% segskeleton.tex
\documentclass[12pt]{article}
%\documentstyle[seg,manuscript]{revtex} % This line gives a manuscript style 
%
%\documentstyle[seg]{revtex}  % Uncomment this line for a preprint style.
\usepackage{wsh}
% \usepackage[dvips]{epsfig}

\renewcommand{\baselinestretch}{1}
\def\listoffigures{}

\gdef\xv{{\svector x}} 
\gdef\xvs{{\svector {\scriptstyle x}}} 
\gdef\grad{{\svector {\nabla}}} 
\gdef\pv{{\svector p}}
\gdef\qv{{\svector q}}
\gdef\zerov{{\svector 0}} 
%%%%%%%%%%%%%%%
\gdef\xs{{\xv_s} }
\gdef\xr{{\xv_r} }
\gdef\xb{{\xv_b} }
\gdef\xh{{\xv_h} }
\gdef\xc{{\xv_c} }
\gdef\gain{\mbox{gain}}
%%%%%%%%%%%%%%%
\gdef\W{{\st{W}}}
%%%%%%%%%%%%%%%

\begin{document}
\bibliographystyle{plain}

\title{Common-offset Depth Migrations\\ and Traveltime Tomography}

\author{William S. Harlan, Douglas W. Hanson, Mark Boyd, \\
Conoco, Inc.}

\date{submitted 1994, accepted 1995, withdrawn}

\maketitle

\section{abstract}
We automatically improve an interval 
velocity model after picking
residual inconsistencies from con\-stant-off\-set depth migrations.
For generality, we
employ a reflection traveltime tomography algorithm,
which allows other applications and other sources of information.

Many methods of depth migration velocity analysis emphasize well-focused
images and use tools similar to semblance stacks.  
Others linearize and invert the effect of perturbed velocities on
migrated images.  We prefer to use developed methods of reflection traveltime
tomography by converting picked migrated reflections into equivalent
multi-offset traveltimes.

Reflection traveltime tomography finds interval velocities and reflection
geometries that best explain observed surface reflection times.
Reflection tomography has evolv\-ed away from layered models toward
independent parameters for velocities and reflectors.  
Interval velocities are parameterized as a smooth function of
spatial coordinates. Reflections are described by a collection of 
common-reflection points, which do not
assume more continuity than necessary to reconstruct picked segments of 
picked reflection times.

Migration facilitates prestack picking by simplifying diffracted reflections 
and dispersing noise.  The effective signal-to-noise
ratio improves.
Depth migration does not add information to reflections, however.  
In fact, the bias of a poor velocity model must be removed by
reconstructing the prestack traveltimes that produced the poor migration.  
To do so, we reconstruct the paths and surface geometries
for each of the picked migrated reflector positions.
Conventional dynamic ray methods or extrapolated traveltime 
tables suffice.  

Constant-offset sections of a North Sea line were independently migrated in
depth and viewed on a 3D interpretive workstation.  One reflection at the base
of chalk imaged at inconsistent depths over offset.  The migrated depths
of this and other reflections were picked over a range of offsets.  
Equivalent prestack traveltimes were
modeled through the migration velocity model.   The chosen method of traveltime
tomography implicitly encouraged consistency in common-reflection points for
raypaths at various offsets.  The final estimated velocity model showed an
increase in velocities near the base of the chalk, then a decrease in velocities
below. Remigration of the data with the revised velocities greatly increased
the visibility of the reflection at the base of the chalk.

Dynamic ray methods and explicit traveltime extrapolations identify
com\-mon-reflec\-tion points that best model prestack traveltimes.  The
error between a modeled and measured traveltime is 
converted into an equivalent positioning error for the reflection
point.  Velocities
are revised to minimize the {\sl variance} of these positioning 
errors for all offsets of each common-reflection point. 

\section{Introduction}
\label{sec:intro}

Velocity analysis of seismic data after prestack depth migration has largely
concentrated on better focused images of reflectors (e.g. Jeannot et al,
1986; Al-Yahya, 1989; and MacKay and Abma, 1989).  Others have formulated
tomographic methods that directly optimize the effect of velocities on migrated
depths (Fowler, 1988; Etgen, 1990; van Trier, 1990).  Velocity models are
expected to produce consistent images in depth from independently migrated
gathers: usually common-offset or common-shot.  Iteratively
linearized inversions can perturb velocity models to reduce these
inconsistencies.  Each of these methods requires an algorithm designed
specifically for depth migration, with no other obvious application.

Alternatively, we prefer to use prestack depth migrations as a source of
information for already existing methods of 
reflection traveltime tomography, such as
Sattlegger et al (1981), Bishop et al (1985), Bording et al (1987), Sword
(1987), Dyer and Worthington (1988), Sherwood (1989), Harlan et al (1989, 1991),
and Stork and Clayton (1991).  These methods usually require
lists of picked reflection times for many source and receiver combinations.
The estimated interval velocities are also used to detect the anomalous
velocities of gas and overpressure, and to correct the distortions
of structure by shallow velocity changes (``buried statics'').
Those interested only in applications to depth migration still 
benefit from simpler algorithms, with broader application,
and with better-understood properties.  Those interested most 
in the interpretation of velocities find that migration improves
the quality of prestack picking.

Few independently developed methods of reflection traveltime tomography share
identical physical parameters, input data, or numerical methods.  This paper
attempts to isolate features that adapt to a variety of data with the fewest
physical constraints.  Sattlegger et al (1981) introduced the tomographic
optimization of layered models: continuous reflectors that vertically delimit
sharp changes in interval velocities, usually with smooth lateral changes. 
With few parameters, layer boundaries and velocities can be optimized
simultaneously.  Sherwood's survey (1989) shows the continuing popularity of
this model.  The first three-dimensional applications (Chiu et al, 1986)
extended the layered model.  

Bishop et al (1985), Bording et al (1987), Dyer and Worthington (1988), and 
Toldi (1989) preferred models that decouple velocities and reflector
geometries.  Velocities can vary continuously, with resolution dependent on
discretization and binning.  Sword (1987), Harlan et al (1989, 1991), 
Biondi (1990),
van Trier (1990), and others avoided continuous reflectors and estimated
common-reflection points.  The additional degrees of freedom raise concerns
about convergence.  Fowler (1988), Etgen (1990), and Stork and Clayton (1991)
carefully analyzed the effect of perturbed velocities on migrated reflection
points and concluded that both must be perturbed simultaneously.  We introduce
a simple method of doing so.  

These papers use a variety of input data: picked prestack traveltimes, picked
prestack depth migrations, constant velocity time migrations, picked ``stacking
velocities,'' 
semblance panels, local slant stacks, and beam stacks. We have been
able to optimize many of these alternative forms of data by treating them as
simple functions of traveltimes.   Although we pick migrated depths from our
example data, we optimize an equivalent set of prestack reflection times.  

\section {An example of depth migration errors}

Figure 1 displays a prestack (Kirchhoff) depth migration of a seismic
line from the Netherlands' North Sea, spanning 11.25~km of midpoints and 5
km depth. Constant-offset sections were migrated independently, then
stacked over offset to produce a single image.  The original velocities
were largely stratified and only increased with depth.  (500 traces
are spaced at 22.5 m---one for each original shot position.)

When the unstacked cube of migrated data was examined on a 3D
interpretative workstation, some reflections were seen to align poorly
over offset.   Figure 2 shows some ``common-image point'' (CIP) gathers. 
Each gather shows the image for a single horizontal position and a range
of depths and offsets (154~m to 2000~m offset).  Note that the reflector
at shot position 400 and 2750~m depth is very inconsistent
over offset.    (Constant-offset depth migrations do not have the
numerical artifacts from edge effects found in shot profile migrations. 
See other differences in Cox and Wapenaar, 1992.)

Figure 3 shows the picks of migrated reflections at various offsets.  
At least five offsets were picked for each reflector, always including a
near offset of 154 m.  The maximum pickable offset increased linearly
from 1300~m at  800~m depth to 3574~m at 4800~m depth.  The grey levels
in figure 3 show the transmission velocity model used to migrate the
data originally.  A few sample reflection raypaths are shown.
(This figure spans the same distances as figure 1.)

The picks of most reflections are almost indistinguishable.  The
reflector near 2500~m depth lies beneath a 1000~m thick interval of
chalk and shows considerable inconsistency over offsets 154~m to 2700~m.  
The chalk velocity cannot be adjusted to flatten this one reflection
without spoiling the images of deeper reflectors.  Although the chosen
velocity model may appear close to a solution, it is not.

\section {Migrating for signal enhancement}

After prestack depth migration, a cube of unstacked reflection seismic
data can become considerably easier to interpret and pick.  Migration
improves signal-to-noise ratios by averaging random noise over midpoint. 
Migration also simplifies reflections from structure with high curvature
(particularly diffractions), reduces overlapping of events, and allows
easier visual correlation over offset.

Depth migration does not add information to observed reflections,
however.  If anything, depth migration adds the bias of a particular
velocity model that, good or bad, describes only our previous
assumptions.  If the migration and ``true'' velocities differ by a
shallow velocity anomaly, for example, then migration will only diffuse
and weaken underlying reflections.

If we choose migration velocities only to improve the quality of picks,
then we may prefer to initialize our velocity optimization with other
models.  First, we must remove the bias of our migration velocities
from the picked migrated depths, so far as possible.  To do so, we 
reconstruct the prestack traveltimes that must
have imaged at the picked migrated depths.   

\section {Reflection times for tomography}

To reconstruct prestack traveltimes from the picked migrated depths in figure
3, we use geometric constant-offset modeling: that is, find surface midpoints
for reflections from picked reflectors with the proper locations, angles,
and offsets. The prestack traveltimes (and their spatial derivatives)
are given by the estimated raypaths through the reference velocity
model. See the appendix for details.

Conventional methods of dynamic ray shooting or relaxation suffice for
this modeling step.  Explicit extrapolation and tabulation of
traveltimes are recommended for their simplicity and speed (Vidale,
1990; van Trier, 1990; Moser, 1991; and Asakawa and Kawanaka, 1993).

Figure 4 shows the corresponding constant-offset time picks modeled from
the reflectors in figure 3.  These picks should be equivalent to the
prestack traveltimes and moveouts in the original unmigrated, unstacked
data.  We can now proceed with a conventional reflection traveltime
tomography, as if these picks were our original data.
The chosen method of reflection traveltime tomography will
implicitly encourage consistent images of common-reflection points.

\section {Describing a velocity function}

We parameterize the transmission slowness $P(\xv)$ 
(reciprocal velocity) as a {\sf
smooth} function of our spatial coordinates $\xv$.
Basis functions, splines, or smoothed grids serve equally well.  
We require only that the continuous slowness be a linear function of
its parameters.  The smoothness of the function should
also be adjustable so that resolution can be increased as an
inversion proceeds and as accuracy increases.

As a concrete example, let discrete coefficients $P_i$ scale
basis functions $g (\xv)$ centered at points $\xv_i$.  The widths of
these basis functions are controlled by a scalar $w$.

\begin{eqnarray} 
& & P(\xv)  \equiv \sum_{i} ~ P_{i} ~ w^{-1} ~ g [ ( \xv - \xv_i ) / w ] ,
\nonumber \\
& \mbox{ where } &  \int g (\xv) d\xv = 1, \mbox{ and }
\int g (\xv) ~ \| \xv \|^2 ~ d\xv \approx 1 . 
\label{eqn:myeqbasis}
\end{eqnarray} 

This basis function has a normalized area and width, so that 
the magnitudes of $P_i$ and $w$ are comparable to the slownesses
and spatial resolution respectively.  Multidimensional Gaussians
are convenient.  This continuous slowness model is 
a linear function of the coefficients, a
convenient property for optimization.  The resolution of this model can be
modified dynamically simply by adjusting the scalar $w$.

\section {Optimizing common-reflection points and velocities}

An unoptimized slowness model will not allow a fan of modeled rays to
share a common-reflection point and explain the measured traveltimes at
all offsets.  Dynamic ray tracing, shooting, and relaxation 
can find reflection paths that fit multi-offset reflection 
times as well as possible.  See the appendix for details.
We prefer the powerful combination of explicit traveltime 
extrapolation (e.g. Vidale, 1990; van Trier,
1900; Moser, 1991) with Fermat's principle 
to estimate representative raypaths (Harlan, 1990).  
Spatial derivatives of measured traveltimes constrain
the dips of reflectors. 

Assume that we have identified many different common-reflection points, 
indexed by $b$.  Each point reflects $N_b$ raypaths with measured traveltimes
$t_{bh}$ at offsets indexed by $h$.  If estimated raypaths are written as a
function of spatial distance $a$, then modeled 
traveltimes are line integrals of slowness along the paths:

\begin{equation}
T_{bh} = \int_{0}^{a_{bh}} P [ \xv_{bh} (a) ] d a.
\label{eqn:myeqtimes}
\end{equation}
For convenience, the 
raypath $\xv_{bh} (a) $ begins with $a=0$ at a source
position, increases along the raypath, through the reflection point, and
reaches the total length $a_{bh}$ of the ray at the receiver location.
This modeled traveltime is also a linear function of the slownesses
and of the parameters that describe these slownesses.

When raypaths do not include reflections, tomography iteratively
linearizes the modeling by holding raypaths constant and considering
only the effect of interval velocities on traveltimes.  Because of
Fermat's principle, perturbations of raypaths do not affect the 
traveltimes to first order.  The position of reflections, however, 
does affect traveltimes to first order.  By requiring perfect
agreement with picked times, we can measure the effect of perturbing
velocities on reflector positions.

In the vicinity of a reflection point, up- and down-going waves can be
approximated as plane waves.  Assume that a reflector has been displaced
perpendicular to its dip until the measured and modeled traveltimes 
($t_{bh}$ and $T_{bh}$) of a raypath agree.  If the up- and down-going 
rays meet at an angle $\theta_{bh}$, then the
following error measures the effect of such a displacement on the zero-offset
(normal-incidence) reflection time:

\begin{equation} 
e_{bh} = (t_{bh} - T_{bh} )/ \cos ({\theta}_{bh}/2) .
\label{eqn:myeqerror}
\end{equation}
See the appendix as well as 
Stork and Clayton (1992) for a justification of the cosine.
Notice that this positioning error increases as the angle 
of reflection increases.  

Since the velocity model is imperfect, we know that our original
positions for reflection points were incorrect.  We do not want
to discourage a new velocity model from moving the reflection points,
but we do want consistency from all offsets that share
a common-reflection point.

A revised velocity model need not drive the positioning errors 
(\ref{eqn:myeqerror}) to zero
but should make the errors depend on the reflection point b alone.  We
want to
find the slowness model that minimizes the {\sl variance} of these errors
over offset:

\begin{equation}\min_{P_{i}} = \sum_{b} \sum_{h} ( e_{bh} - {1 \over {N_b} }
\sum_{h'} e_{bh'} )^2 .
\label{eqn:myeqobj}
\end{equation}
Analogously, prestack depth migration 
must create consistent images from
different offsets, without constraining the depth of reflectors.  This
quadratic function of slowness lends itself to least-squares methods like
conjugate gradients or singular-value decomposition.

Figure 5 shows estimated transmission velocities and reflection geometries. 
These estimated depths vary much less over offset 
than do the original picks in figure 3.  Picks are discarded
if the range of offsets is inadequate to constraint a particular
reflection point.  (Single offsets and nearly identical offsets
do not harm the optimization, but do not help either.)
Figure 6 shows the subtraction of the original velocities in figure 2 from the
estimated velocities in figure 5.
A single reflector location is able to fit
modeled traveltimes to within a quarter wavelength. 
Note that velocity increases near the bottom of the chalk, then decreases again
below.  Well logs in the area show similar changes in chalk velocities. 

Figure 7 shows a remigration of the data with revised velocities.  This time,
the reflection at the bottom of the chalk appears very strong and coherent,
as it does before stack.  The common-image point gathers in figure 8 show
greater consistency over offset.  Although a few shallower reflections seem
slightly less coherent before stack, the residual inconsistencies are
distributed much more evenly.

No further iteration was necessary.  If substantial inconsistencies had
remained over offset, then repicking would not have helped unless new
reflections became visible before stack.  In this case, revised velocities
affected only the migrated depths of reflectors before stack, not their
coherence or strength.


\section{Recommendations}

The example in this paper was chosen to demonstrate the 
equivalence of depth migration velocity analysis and reflection tomography.
Most of our applications of reflection tomography begin with
densely picked stacking functions that best describe the
unmigrated prestack moveouts of reflections over offset.
The following guidelines are appropriate:

\begin{enumerate}
\item   To avoid time-consuming hand optimization of prestack depth 
migration velocities, use tomographic velocity estimations whenever 
possible.

\item   Use post-migration picks when unmigrated data 
are too noisy for prestack interpretation, or when 
complex structure overlaps considerably in time.

\item   Use post-migration picks to improve an already existing
interval velocity model that requires some minor improvement.

\item  Use unmigrated prestack traveltimes to estimate an interval
velocity model from scratch, when data quality allows.  

\item  Pick data prior to migration when shallow 
lateral velocity anomalies are likely.  
(Migration will destroy evidence of ``time sags.'')

\item  When densely picked ``stacking velocities'' are available (twenty per
cablelength), tomographically estimate an interval velocity routinely for 
depth migration or conversion.

\end{enumerate}

\section{Conclusions}

Already existing tools for reflection traveltime tomography are easily
adapted to prestack migrated data.  Migration eases picking by improving
signal-to-noise ratios and by simplifying the appearance of reflections. 
Those interested only in migrated images will benefit from using a more
general algorithm, capable of incorporating traveltime information from other
sources.  When the initial velocity model is poor, some reflections
may be easier to pick without migration.
Post-migration picks can be
converted and combined with pre-migration picks, and even 
with picks from ``stacking
velocity'' analyses.  One tomographic algorithm can serve for many varieties
of data.  

No repicking of data appears to be necessary, except to 
eliminate multiples, cycle skipping, and other mistakes.  
Traveltime tomography is sufficiently
iterative to allow for the non-linearities of ray-bending, constrained
velocities, and so on.  If tomographically estimated velocities and reflectors
do not fit the picked data, then the picks may not be consistent with the
physical assumptions.  Tomography provides the best estimate
of migrated depths from surface information alone.  Focusing analysis can
remove any remaining unexplained inconsistencies.
Tools also exist for interpretive modification of the best tomographic
model, particularly to add or adjust sharp velocity
contrasts, such as salt interfaces.

Identifying common-reflection points improves the robustness and convergence of
estimated interval velocities.  Errors in modeled traveltimes can be converted
into equivalent displacements of the reflection point for each raypath.  An
optimum velocity model encourages these displacements to be as consistent as
possible, without attempting to preserve the original positions.

%\section{Acknowledgments}

\section{REFERENCES}
\def\reference#1\par{
  \hangafter=1
  \hangindent2em
  #1
}
\parindent=0pt

\def\segmeeting{{Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded
Abstracts,~}}
\reference 
Al-Yahya, K., 1989, Velocity analysis by iterative profile migration:
Geophysics, {\bf 54}, \mbox{718--729.}
\par

\reference
Asakawa, E, and Kawanaka, T, 1993, Seismic ray tracing using linear
traveltime interpolation: Geophysical Prospecting, {\bf 41}, 99--111.
\par

\reference
Biondi, B., 1990, Seismic velocity estimation by beam stack: Ph.D.
thesis, Stanford Univ.
\par

\reference
Bishop, T.N., Bube, K.P., Cutler, R.T., Langan, R.T., Love, P.L.,
Resnick, J.R., Shuey, R.T., Spindler, D.A., and Wyld, H.W., 1985, Tomographic
determination of velocity and depth in laterally varying media:
Geophysics, {\bf 50}, 903--923.
\par

\reference
Bording, R.P., Gersztenkorn, A., Lines, L.R., Scales, J.A., and Treitel,
S., 1987, Applications of seismic traveltime tomography: Geophys. J.
Roy. Astr. Soc., {\bf 90}, 285--303.
\par

\reference
Chiu., S.K.L., Kanasewich, E.R., and Phadke, S., 1986, Three-dimensional
determination of structure and velocity by seismic tomography:
Geophysics, {\bf 51}, 1559--1571.
\par

\reference
Cox, H.L.H., and Wapenaar, C.P.A., 1992, Macro model estimation by common
offset migration and by shot record migration: Journal of Seismic Exploration,
{\bf 1}, 29--37.
\par

\reference
Dyer, B.C., and Worthington, M.H., 1988, Seismic reflection tomography:
a case study: First Break, {\bf 6}, 354--366.
\par

\reference
Etgen, 1990, Residual prestack migration and interval velocity
estimation: Ph.D. thesis, Stanford Univ.
\par

\reference
Fowler, P.J, 1988, Seismic velocity estimation using prestack time
migration, Ph.D. thesis, Stanford Univ.
\par

\reference
Harlan, W.S., 1989, Tomographic estimation of seismic velocities from
reflected raypaths: 59th \segmeeting 922--924.
\par

\reference
Harlan, W.S., 1990, Tomographic estimation of shear velocities from
shallow cross-well seismic data: 60th \segmeeting 86--89. 
\par

\reference
Harlan, W.S., Hanson, D.W., and Boyd, M., 
1991, Traveltime tomography with multioffset
common-reflection points: 61st \segmeeting 974--976.
\par

\reference
Jeannot, J.P., Faye, J.P., and Denelle, E., 1986, Prestack migration
velocities from depth-focusing analysis: 56th \segmeeting 1258--1261.
\par

\reference
Liu, Z., and Bleistein, N., 1995, Migration velocity analysis:
Theory and an iterative algorithm:
Geophysics, {\bf 60}, 142--153.
\par

\reference
MacKay, S., and Abma, R., 1989, Refining prestack depth migration
without remigration: 59th \segmeeting 1258--1261.
\par

\reference
Matsuoka, T., and Ezaka, T., 1992, Ray tracing using reciprocity,  Geophysics,
{\bf 57}, 326--333
\par

\reference
Moser, T.J., 1991, Shortest path calculation of seismic rays:
Geophysics, {\bf 56}, 59--67.
\par

\reference
Sattlegger, J.W., Rohde, J., Egbers, H., Dohr, G.P. Stiller, P.K., and
Echterhoff, J.A., 1981, INMOD---Two dimensional inverse modeling based
on ray theory: Geophys. Prosp., {\bf 19}, 229--240.
\par

\reference
Sherwood, J.W.C., 1989, Depth sections and interval velocities from
surface seismic data: Leading Edge, {\bf 8}, 44--49.
\par

\reference
Stork, C., and Clayton, R.W., 1991, Linear aspects of tomographic
velocity analysis: Geophysics, {\bf 56}, 483--495.
\par

\reference
Sword, C.H., 1987, Tomographic determination of interval velocities from
reflection seismic data: the method of controlled directional reception:
Ph.D. thesis, Stanford Univ.
\par

\reference
Toldi, J.L., 1989, Velocity analysis without picking: Geophysics, {\bf
54}, 191--199.
\par

\reference
van Trier, J.A., 1990, Tomographic determination of structural velocities 
from depth migrated seismic data: Ph.D. thesis, Stanford Univ.
\par

\reference
Vidale, J.E., 1990, Finite-difference calculation of traveltimes in
three dimensions: Geophysics, {\bf 55}, 521--526.
\par

%\reference
%%% the following line should be used when using bibtex, then citations will 
%%% be pulled in from the file called: segmaster.bib
%%%
%\bibliography{segmaster}

%%% this block should be commented out when using bibtex
%%%
%\begin{thebibliography}{}
%\bibitem[\protect\citeauthoryear{Lindsey}{1993}]{lindsey93}
%Lindsey, J. P, Jr., 1993, Instructions to authors: Geophysics, {\bf 58}, 2--9.
%\footnote{We encourage the use of BIB\TeX\ to create References. See
%the paper {\tt easybib.tex} for more details.}
%
%\bibitem[\protect\citeauthoryear{Claerbout}{1991}]{clr91}
%Claerbout, J. F., 1991, A scrutiny of the introduction: The Leading Edge, 
%{\bf 10}, 39--40.
%
%\bibitem[\protect\citeauthoryear{Landes}{1966}]{lan66}
%Landes, K. K., 1966, A scrutiny of the abstract, {II}: AAPG Bull., {\bf50},
 %1992.
%
%\bibitem[\protect\citeauthoryear{Smith and Wesson}{1911}]{smith11}
%Smith, A. B., and Wesson, W. E, 1911, Better than four aces: Texas Review, 
%{\bf 1}, 1--10.
%
%\end{thebibliography}
%%%

\appendix

\section{}
\label{app:first}

%
%\begin{equation}
%\label{eqn:myeqxxx}
%\end{equation}
%
%
%\tracingcommands=1
%\gdef\marginal#1{}

In this appendix, we fill in algorithmic details
with a notation chosen to minimize ambiguity.  
First, we define how partial prestack depth migrations transform
data with a summation (``Kirchhoff'') formulation.  Then 
we relate the effect of migration on coherent reflections
to the raypath approximations used by traveltime tomography.
Picks of constant-offset migrated depths are used to find equivalent
picked prestack reflection times, and vice-versa.
Finally, we examine how perturbations of 
reflector locations affect the modeled traveltimes, so that
tomography can simultaneously optimize reflection points and interval
velocities.

\subsection {Prestack depth migration}

Seismic amplitudes $u(t, \xs, \xr)$  (displacement or pressure)
are recorded as a function of time $t$ at the surface 
source and receiver positions indexed by $s$ and $r$.  
The Cartesian elements of a coordinate vector $\xv$ are $(x,y,z)$,
where $z$ increases with depth.
For each surface source or receiver position ($\xs$ or $\xr$)
we extrapolate a table of traveltimes 
$T ( \xv , \xs )$ 
to many buried positions $\xv$.

Traveltimes are understood to satisfy an Eikonal equation.  
The gradient of traveltime has a magnitude equal to 
reciprocal velocity, or slowness:

\begin{equation}
1/v( \xv ) \equiv \Bigl \| \grad _{\xvs} T [ \xv , \xv' ] \Bigr \| 
\label{eqn:myeqeik}
\end{equation}
The Eikonal equation is accompanied by 
transport equations, which
specify the geometric changes in amplitude $R (\xv , \xs )$.
The arguments of $T$ and $R$ both can be reversed symmetrically 
(a result of reciprocity).
Single-valued functions such as these do not allow caustics
or multiple arrivals.
By making the slowness and velocity independent of $\xv'$
we also assume {\sl isotropy}.

The data are assumed to be a linear function of the migrated image

\begin{equation}
u (t, \xs, \xr) = \int d \xv  ~ m( \xv ) ~ \delta [t - T ( \xv ,
\xs ) - T ( \xv , \xr ) ] ~ R( \xv , \xs ) ~ R( \xv , \xr ) ~ \gain(t) .
\label{eqn:myeqmodel}
\end{equation}
The recorded
data are usually scaled by an increasing function of time, such as 
$\gain(t) = t^2$, to reduce the dynamic range.  Intentionally,
this gain cancels some of the scaling by geometric factors.

A generalized inverse of this linear equation would be preferable,
but for efficiency, a modified adjoint operation gives an approximate inverse.
This summation method is often called a 
``Kirchhoff'' method although it need not
use the integral and approximation by that name.  
The image locations will be indexed by $b$.

\begin{equation}
\hat m ( \xb ) =  
\sum_{s,r} \int dt~ \dot u (t, \xs, \xr)
~ \delta [ t - T ( \xb , \xs ) -  T (\xb , \xr) ] 
~ R( \xb , \xs ) ~ R( \xb , \xr ) ~ \gain (t)
.
\label{eqn:myeqmigration}
\end{equation}
The summation is over 
recorded source and receiver positions.
A time differentiation of the data (a ``rho'' filter)
partially corrects the phase distortion of the model.

For our purposes, a partial migration will be more useful.
We find it useful to perform the summation
over the midpoint coordinate $\xc \equiv ( \xr + \xs)/2$
rather than source position.  
An image at a constant ``half offset'' 
$ \xh \equiv (\xr - \xs)/2 $  restricts the summation
to source and receivers with a constant separation:

\begin{eqnarray} 
& \hat m_h ( \xb )  = &
 \sum_{c} \int dt~ \dot u (t, \xc - \xh , \xc + \xh)
~ \delta [ t - T ( \xb , \xc -\xh) -  T (\xb, \xc + \xh)  ] 
~ \cdot \nonumber \\
& & \cdot ~ R( \xb , \xc-\xh ) ~ R( \xb , \xc + \xh ) ~ \gain (t) .
\label{eqn:myeqoffmig}
\end{eqnarray} 
Similarly, we can remodel 
data with different versions of the constant offset migrations:

\begin{eqnarray} 
& \hat u (t, \xc - \xh , \xc + \xh) = &
 \sum_{b} \hat m_h ( \xb )  
~ \delta [ t - T ( \xb , \xc -\xh) -  T (\xb, \xc + \xh)  ] 
~ \cdot 
\nonumber \\
& & \cdot ~ R( \xb , \xc-\xh ) ~ R( \xb , \xc + \xh ) ~ \gain (t) .
\label{eqn:myeqoffmod}
\end{eqnarray} 
When the traveltime table is consistent with the data,
the constant-offset images $\hat m_h (\xv_b)$ 
should not show changes in phase over different offsets $\xh$.
For geometric discussions of phase delays, 
we can ignore the smoothly varying
gain and geometric scale factors.

\subsection {Reconstructing raypaths from traveltimes}

Usually, one constructs a traveltime table $T$ from
a particular velocity model.  To study the properties of
the transforms (\ref{eqn:myeqmodel}) through (\ref{eqn:myeqoffmod}), we will
find it useful to take the traveltime table as given and 
deduce other properties from it.  We will then find it easier
to improve the velocity model and traveltime table.  

Define a slowness vector $\pv$ by treating the traveltime table $T$ as
a scalar potential field:

\begin{equation}
\pv [ \xv , \xv' ] \equiv \grad_{\xvs} T[ \xv , \xv' ]
\label{eqn:myeqpdef}
\end{equation}
and 

\begin{equation}
T[ \xv , \xv' ] =
\int_{\xvs'}^{\xvs} 
\pv [ \xv'' , \xv' ] \cdot d \xv'' ,
\label{eqn:myeqpint}
\end{equation}
where the line integral can follow any path.
By construction, the vector slowness is irrotational (waves should not
travel in a loop): $\grad \times \pv = \zerov$.  

The magnitude of the slowness vector
is the slowness $P$, the reciprocal of
the local velocity---a restatement of the Eikonal equation:

\begin{equation}
P (\xv ) \equiv \Bigl \| \pv [ \xv , \xv' ] \Bigr \| 
\label{eqn:myeqeikp}
\end{equation}
To derive traveltimes tables from local slownesses, we 
need constants of integration.  We can extrapolate a unique
traveltime table $T$ from $P$ if traveltimes are specified on
a point, curve, or surface, and if traveltimes
satisfy Laplace's equation $\nabla^2 T = 0$ elsewhere (sourceless).
Unfortunately, caustics of crossing slowness vectors 
easily form during extrapolation, producing 
multivalued traveltimes.  In practice, single-valued extrapolations
select either minimum traveltimes or those with the strongest
geometric scale factors.

Let a raypath $\xv (a)$ be parameterized as a function of spatial
distance $a$, so that $\bigl \| d \xv (a) / da \bigr \| \equiv 1$.  
The raypath should also be 
be tangent to any slowness vector that originates from
another point on the path:

\begin{equation}
{d \over {da}}  \xv (a)   \equiv 
\pv [ \xv (a) , \xv ( a_0 ) ] / 
  P [ \xv (a) ]  .
\label{eqn:myeqray}
\end{equation}
Thus,

\begin{eqnarray} 
&T [  \xv (a) , \xv (a_0 ) ] 
&= \int_{a_0}^{a} da' 
~ { d \over {da'} } T [ \xv ( a' ) , \xv ( a_0 ) ]  
\nonumber \\
& &= \int_{a_0}^{a} da' 
~ \grad T [ \xv (a') , \xv ( a_0 ) ] \cdot { d \over {da'} } \xv ( a' ) 
\nonumber \\
& &= \int_{a_0}^{a} da' 
~ \Bigl \| \pv [ \xv (a') , \xv ( a_0 ) ] \Bigr \| 
~ \Bigl \| { d \over {da'} } \xv ( a' ) \Bigr \| 
\nonumber \\
& &= \int_{a_0}^{a} da' 
~  P [ \xv (a') ]  
\label{eqn:myeqintt}
\end{eqnarray} 
We have the conventional result 
that the traveltime is the integration
of slowness along a raypath.  

The raypath was defined as tangent
to the slowness vector, but we could make the equivalent assumption
that the final integral in equation (\ref{eqn:myeqintt}) is stationary with
respect to the raypath $\xv$ (minimum traveltime).  
The calculus of variations allows us to reverse the derivation.

To extrapolate a raypath from a point $\xv(a_0)$ in a known direction
$\pv(a_0)$, we can use equation (\ref{eqn:myeqray}) and the following,
which derives from (\ref{eqn:myeqray}) and (\ref{eqn:myeqeikp}):

\begin{equation}
{d \over da} \pv [ \xv (a) , \xv (a_0) ] = \grad_{\xvs} P [ \xv (a) ]
.
\label{eqn:myeqrayp}
\end{equation}
This equation describes how a ray is bent by local changes
in slowness (Snell's Law).  Dynamic ray tracing uses finite differences
to extrapolate the ray differential equations, (\ref{eqn:myeqray}) 
and (\ref{eqn:myeqrayp}).
Other methods include shooting, relaxation, and the reciprocity
method (which we use), described in Harlan (1990) and
Matsuoka and Ezaka (1992).

\subsection {Residual geometric modeling and migration}

After performing the constant offset migration in (\ref{eqn:myeqoffmig}),
we identify the same continuous reflector at several constant
offsets.  
We pick the migrated positions of this reflector 
$\xv_{bh} = [x_b, y_b, z_{bh}]$ at a fixed lateral 
position $(x_b, y_b)$ and allow the depth $z_{bh}$ to change with
offset index $h$. 
Each coherent pick is indexed by $b$.
Let us also pick the local dip 
with a vector $\qv_{bh}$ that is normal to the migrated reflector.
For convenience, assume a unit magnitude: $\| \qv_{bh} \| \equiv 1$.
Locally, the coherence of this reflection could be approximated to first
order as a planar surface:

\begin{equation}
\hat m_h ( \xv_{b} ) \approx f [ ( \xv_b - \xv_{bh} ) \cdot \qv_{bh} ]
\label{eqn:myeqmoddip}
\end{equation}
where $f(\cdot)$ is a simple wavelet describing the
local coherence perpendicular to the surface.

All our picked data, such as found in figure 3, will be
summarized as a list of \mbox{$\{\xv_{bh}$, $\qv_{bh}\}$,} for many $b$ and
$h$.  Migrated 
reflectors need only be continuous enough over $\xv_b$ to allow
the picking of a local dip.
What coherence in the original unmigrated data would have produced
these picks? Can we derive a set of equivalent unmigrated traveltime
picks?  

We will find it easiest to answer these questions by
seeing how equation (\ref{eqn:myeqoffmod}) remodels the data.
The migrated reflection point $\xv_{bh}$ contributes to all 
source and receiver pairs with fixed ``half offsets''
$\xh = (\xv_r - \xv_s)/2$.
For each affected midpoint $\xv_c = (\xv_r + \xv_s)/2$,
we can draw a raypath from the source and receiver to the
reflection point.  The two rays reach the reflection point
with known slowness vector directions
$\pv ( \xv_{bh} , \xv_c + \xh )$ and $\pv ( \xv_{bh} , \xv_c - \xh )$.

By looking for a stationary phase in the constant-offset modeling 
integral (\ref{eqn:myeqoffmod}) with the approximation (\ref{eqn:myeqmoddip}), 
we find
this reflection point contributes most to the midpoint which maximizes

\begin{equation}
\max_{\xvs_c}
~ \Bigl \|
[\pv ( \xv_{bh} , \xc - \xh ) 
+\pv ( \xv_{bh} , \xc + \xh )] 
\cdot \qv_{bh}
\Bigr \|
\label{eqn:myeqmaxc}
\end{equation}
In other words, the rays should reflect
symmetrically about the normal to the reflector.
Compare this argument to that of Liu and Bleistein (1995).
When this dot product is maximized we find that

\begin{eqnarray} 
 \pv ( \xv_{bh} , \xc - \xh ) 
+\pv ( \xv_{bh} , \xc + \xh ) 
= 2   ~ P (\xv_{bh}) ~ \cos ( \theta_{cbh} / 2) ~
 \qv_{bh} ,  \nonumber \\
 \mbox { where } \cos ( \theta_{cbh} ) \equiv 
[ \pv ( \xv_{bh} , \xc - \xh ) \cdot
\pv ( \xv_{bh} , \xc + \xh ) ] / P ( \xv_{bh} )^2 .
\label{eqn:myeqcosine}
\end{eqnarray}
$\theta_{cbh}$ gives the angle between
the two raypaths as they meet at the reflection point.  

The total traveltime of a reflection is given by

\begin{equation}
t_{ch} (\xv_{bh}) = T(\xv_{bh} , \xc - \xh ) + T(\xv_{bh} , \xc + \xh ) .
\label{eqn:myeqtimefit}
\end{equation}
We see how to reconstruct raypaths, 
traveltimes, and surface positions from picks of migrated reflectors.  
For completeness, we outline how to reverse this procedure.

Let us define an equivalent set of traveltime picks from the original 
unmigrated data.
For each offset $\xh$ and midpoint $\xc$ we pick a traveltime
$t_{ch}$.  According to the migration equation (\ref{eqn:myeqoffmig}), 
this pick affects all migrated positions $\xv_{bh}$ along
the arc described by equation (\ref{eqn:myeqtimefit}).
To determine which of these midpoints contribute most, we
require more information.

We can also pick a dip of traveltime with respect to midpoint 
$\pv_{ch} = \grad_{\xvs_c} t_{ch} $ where
$\grad_{\xvs_c} \equiv (\grad_{\xvs_s} +  \grad_{\xvs_r})/2 $.
We assume a corresponding coherence in the data and 
look for stationary phase in equation (\ref{eqn:myeqoffmig}).
The position along the arc in (\ref{eqn:myeqtimefit}) that contributes 
most to the picked reflection maximizes

\begin{equation}
\max_{\xvs_{bh}}
~ \Bigl \| 
[ \pv( \xc - \xh, \xv_{bh} ) 
+ \pv( \xc + \xh, \xv_{bh} ) ]
\cdot \pv_{ch} 
\Bigr \|
\label{eqn:myeqdipfit}
\end{equation}
Thus, a constant-offset time
pick $\{t_{ch},\pv_{ch}\}$
or migrated pick
$\{\xv_{bh}, \qv_{bh}\}$ are interchangeable, and can be used
to derive each other and 
construct the same set of raypaths.  
To distinguish traveltimes that are reconstructed from migrated picks,
we use the index $t_{bh}$ in equation (\ref{eqn:myeqerror}) in the main text
to abbreviate $t_{ch} (\xv_{bh})$.
The stationary phase 
approximations make the same high frequency assumptions as the 
Eikonal and ray equations, and all fail in similar situations.

\subsection{Converting time errors to reflector errors}

A perturbation of the reflection point
will perturb the reflection traveltime (\ref{eqn:myeqtimefit}) according to

\begin{equation}
\Delta t_{ch} (\xv_{bh}) \equiv 
t_{ch} (\xv_{bh} + \Delta \xv_{bh} ) -
t_{ch} (\xv_{bh} ) 
= 2 P (\xv_{bh} ) \cos ( \theta_{cbh} / 2 ) \qv_{bh} \cdot \Delta
\xv_{bh} .
\label{eqn:myeqpert}
\end{equation}
Only a perturbation perpendicular to the reflector will
affect the traveltime.  
We have effectively assumed that the wavefront is planar in the vicinity
of the reflection point.

Traditionally, tomography minimizes
errors between mod\-el\-ed and measured traveltimes.  Instead, we
can convert traveltime errors into equivalent errors in reflector
positions:

\begin{equation}
\Delta \hat \xv_{cbh} =
{ \Delta t_{ch} (\xv_{bh}) \over
 2 P (\xv_{bh} ) \cos ( \theta_{cbh} / 2 ) } ~ \qv_{bh} .
\label{eqn:myeqerrx}
\end{equation}
More conveniently still, we can measure these errors
in reflector positions by the change in traveltime of a normal
reflection:

\begin{equation}
e_{cbh} \equiv 
{ \Delta t_{ch} (\xv_{bh}) \over
 \cos ( \theta_{cbh} / 2 ) }  .
\label{eqn:myeqerrt}
\end{equation}
As we optimize the slowness model, 
we do not wish to minimize these errors in reflector positions absolutely
because we do not know the correct absolute location. Rather we wish
the locations to be consistent over offset, with a minimum {\sl variance}
in position errors, as in equation (\ref{eqn:myeqobj}).

\newpage \subsection{\mbox{\hfill FIGURES \hfill}}
\def\figurex#1. #2 {\vfil \eject ~  \vfill FIG. #1. #2 }
\def\figure#1. #2 {\par FIG. #1. #2 }

\figure 1. Prestack depth migration of Netherlands North Sea data with a
simple stratified velocity model.  Shotpoint units index
the locations of seismic sources, which are spaced at 22.5~m.

\figure 2. Common-image-point gathers.  
Constant-offset migrated images were sorted
over offset at selected common-midpoint locations.  Note inconsistent imaging
of reflector at position 400 and 2750~m depth.

\figure 3. A simple stratified velocity model with picked constant-offset
migrated depths.  At least five offsets were picked for each reflection,
including a near offset of 154~m (distance between source and
receiver).  The maximum pickable offset increased
from 1300~m at 800~m depth to 3575~m at 4800~m depth.  Note the
inconsistency of depths at different 
offsets for the reflection near 2300~m depth.

\figure 4. A reconstruction of constant-offset traveltimes from the 
constant-offset picks and velocities in figure 3.  These are sufficient
data for traveltime tomography.

\figure 5. An iteratively optimized model for the interval velocities
and migrated reflection depths.  The consistency of reflectors over
offset has improved.

\figure 6. A subtraction of the original velocity model in figure 3
from the final estimated velocity model in figure 5.  Note that the
velocity has increased above 2500~m depth and decreased below.

\figure 7. A revised prestack depth migration with the optimized 
interval velocity model in figure 5.  The previously weak reflector
near 2500~m depth is now very strong.  (Local gain weakens some
neighboring reflectors.)

\figure 8. Revised common-image point gathers.  Errors in residual
moveout are much better distributed. 

\end{document}
