Common-offset Depth Migrations
and Traveltime Tomography

We automatically improve an interval
velocity model after picking
residual inconsistencies from constant-offset 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 evolved 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
common-reflection 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 variance of these positioning
errors for all offsets of each common-reflection point.

Introduction

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.

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.

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.

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.

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.

We parameterize the transmission slowness
(reciprocal velocity) as a smooth function of our spatial coordinates .
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 scale
basis functions
centered at points
. The widths of
these basis functions are controlled by a scalar .

(1)

This basis function has a normalized area and width, so that
the magnitudes of and 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 .

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 . Each point reflects raypaths with measured traveltimes
at offsets indexed by . If estimated raypaths are written as a
function of spatial distance , then modeled
traveltimes are line integrals of slowness along the paths:

(2)

For convenience, the
raypath
begins with at a source
position, increases along the raypath, through the reflection point, and
reaches the total length 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
( and ) of a raypath agree. If the up- and down-going
rays meet at an angle , then the
following error measures the effect of such a displacement on the zero-offset
(normal-incidence) reflection time:

(3)

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
(3) to zero
but should make the errors depend on the reflection point b alone. We
want to
find the slowness model that minimizes the variance of these errors
over offset:

(4)

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.

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:

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.

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.

Seismic amplitudes
(displacement or pressure)
are recorded as a function of time at the surface
source and receiver positions indexed by and .
The Cartesian elements of a coordinate vector are ,
where increases with depth.
For each surface source or receiver position (
or
)
we extrapolate a table of traveltimes
to many buried positions .

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

(5)

The Eikonal equation is accompanied by
transport equations, which
specify the geometric changes in amplitude
.
The arguments of and 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
we also assume isotropy.

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

(6)

The recorded
data are usually scaled by an increasing function of time, such as
, 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 .

(7)

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
rather than source position.
An image at a constant ``half offset''
restricts the summation
to source and receivers with a constant separation:

(8)

Similarly, we can remodel
data with different versions of the constant offset migrations:

(9)

When the traveltime table is consistent with the data,
the constant-offset images
should not show changes in phase over different offsets
.
For geometric discussions of phase delays,
we can ignore the smoothly varying
gain and geometric scale factors.

Usually, one constructs a traveltime table from
a particular velocity model. To study the properties of
the transforms (6) through (9), 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 by treating the traveltime table as
a scalar potential field:

(10)

and

(11)

where the line integral can follow any path.
By construction, the vector slowness is irrotational (waves should not
travel in a loop):
.

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

(12)

To derive traveltimes tables from local slownesses, we
need constants of integration. We can extrapolate a unique
traveltime table from if traveltimes are specified on
a point, curve, or surface, and if traveltimes
satisfy Laplace's equation
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
be parameterized as a function of spatial
distance , so that
.
The raypath should also be
be tangent to any slowness vector that originates from
another point on the path:

(13)

Thus,

(14)

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 (14) is stationary with
respect to the raypath (minimum traveltime).
The calculus of variations allows us to reverse the derivation.

To extrapolate a raypath from a point
in a known direction
, we can use equation (13) and the following,
which derives from (13) and (12):

(15)

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, (13)
and (15).
Other methods include shooting, relaxation, and the reciprocity
method (which we use), described in Harlan (1990) and
Matsuoka and Ezaka (1992).

After performing the constant offset migration in (8),
we identify the same continuous reflector at several constant
offsets.
We pick the migrated positions of this reflector
at a fixed lateral
position and allow the depth to change with
offset index .
Each coherent pick is indexed by .
Let us also pick the local dip
with a vector
that is normal to the migrated reflector.
For convenience, assume a unit magnitude:
.
Locally, the coherence of this reflection could be approximated to first
order as a planar surface:

(16)

where 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
,
, for many and
. Migrated
reflectors need only be continuous enough over
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 (9) remodels the data.
The migrated reflection point
contributes to all
source and receiver pairs with fixed ``half offsets''
.
For each affected midpoint
,
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
and
.

By looking for a stationary phase in the constant-offset modeling
integral (9) with the approximation (16),
we find
this reflection point contributes most to the midpoint which maximizes

(17)

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

(18)

gives the angle between
the two raypaths as they meet at the reflection point.

The total traveltime of a reflection is given by

(19)

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
and midpoint
we pick a traveltime
. According to the migration equation (8),
this pick affects all migrated positions
along
the arc described by equation (19).
To determine which of these midpoints contribute most, we
require more information.

We can also pick a dip of traveltime with respect to midpoint
where
.
We assume a corresponding coherence in the data and
look for stationary phase in equation (8).
The position along the arc in (19) that contributes
most to the picked reflection maximizes

(20)

Thus, a constant-offset time
pick
or migrated pick
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 in equation (3) in the main text
to abbreviate
.
The stationary phase
approximations make the same high frequency assumptions as the
Eikonal and ray equations, and all fail in similar situations.

More conveniently still, we can measure these errors
in reflector positions by the change in traveltime of a normal
reflection:

(23)

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 variance
in position errors, as in equation (4).

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.

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

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

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.

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.)

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