Common-offset Depth Migrations
and Traveltime Tomography

William S. Harlan, Douglas W. Hanson, Mark Boyd,
Conoco, Inc.

submitted 1994, accepted 1995, withdrawn

abstract

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.

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.

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.

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.

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.

Describing a velocity function

We parameterize the transmission slowness $P({\svector x})$ (reciprocal velocity) as a smooth function of our spatial coordinates ${\svector x}$. 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 ({\svector x})$ centered at points ${\svector x}_i$. The widths of these basis functions are controlled by a scalar $w$.

$$P({\svector x}) \equiv \sum_{i} ~ P_{i} ~ w^{-1} ~ g [ ( {\svector x}- {\svector x}_i ) / w ] ,$$

$$\mbox{ where } \int g ({\svector x}) d{\svector x}= 1, \mbox{ and } \int g ({\svector x}) ~ \Vert {\svector x}\Vert^2 ~ d{\svector x}\approx 1 .$$ (1)

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

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:

  $$T_{bh} = \int_{0}^{a_{bh}} P [ {\svector x}_{bh} (a) ] d a.$$ (2)

For convenience, the raypath ${\svector x}_{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:

  $$e_{bh} = (t_{bh} - T_{bh} )/ \cos ({\theta}_{bh}/2) .$$ (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:

  $$\min_{P_{i}} = \sum_{b} \sum_{h} ( e_{bh} - {1 \over {N_b} } \sum_{h'} e_{bh'} )^2 .$$ (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.

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:

1. To avoid time-consuming hand optimization of prestack depth migration velocities, use tomographic velocity estimations whenever possible.

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

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

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

5. Pick data prior to migration when shallow lateral velocity anomalies are likely. (Migration will destroy evidence of “time sags.”)

6. When densely picked “stacking velocities” are available (twenty per cablelength), tomographically estimate an interval velocity routinely for depth migration or conversion.

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.

REFERENCES

Al-Yahya, K., 1989, Velocity analysis by iterative profile migration: Geophysics, 54, 718–729.

Asakawa, E, and Kawanaka, T, 1993, Seismic ray tracing using linear traveltime interpolation: Geophysical Prospecting, 41, 99–111.

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

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, 50, 903–923.

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., 90, 285–303.

Chiu., S.K.L., Kanasewich, E.R., and Phadke, S., 1986, Three-dimensional determination of structure and velocity by seismic tomography: Geophysics, 51, 1559–1571.

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, 1, 29–37.

Dyer, B.C., and Worthington, M.H., 1988, Seismic reflection tomography: a case study: First Break, 6, 354–366.

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

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

Harlan, W.S., 1989, Tomographic estimation of seismic velocities from reflected raypaths: 59th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 922–924.

Harlan, W.S., 1990, Tomographic estimation of shear velocities from shallow cross-well seismic data: 60th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 86–89.

Harlan, W.S., Hanson, D.W., and Boyd, M., 1991, Traveltime tomography with multioffset common-reflection points: 61st Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 974–976.

Jeannot, J.P., Faye, J.P., and Denelle, E., 1986, Prestack migration velocities from depth-focusing analysis: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1258–1261.

Liu, Z., and Bleistein, N., 1995, Migration velocity analysis: Theory and an iterative algorithm: Geophysics, 60, 142–153.

MacKay, S., and Abma, R., 1989, Refining prestack depth migration without remigration: 59th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1258–1261.

Matsuoka, T., and Ezaka, T., 1992, Ray tracing using reciprocity, Geophysics, 57, 326–333

Moser, T.J., 1991, Shortest path calculation of seismic rays: Geophysics, 56, 59–67.

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., 19, 229–240.

Sherwood, J.W.C., 1989, Depth sections and interval velocities from surface seismic data: Leading Edge, 8, 44–49.

Stork, C., and Clayton, R.W., 1991, Linear aspects of tomographic velocity analysis: Geophysics, 56, 483–495.

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

Toldi, J.L., 1989, Velocity analysis without picking: Geophysics, 54, 191–199.

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

Vidale, J.E., 1990, Finite-difference calculation of traveltimes in three dimensions: Geophysics, 55, 521–526.

.

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.

Prestack depth migration

Seismic amplitudes $u(t, {{\svector x}_s} , {{\svector x}_r} )$ (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 ${\svector x}$ are $(x,y,z)$, where $z$ increases with depth. For each surface source or receiver position ( ${{\svector x}_s}$ or ${{\svector x}_r}$) we extrapolate a table of traveltimes $T ( {\svector x}, {{\svector x}_s} )$ to many buried positions ${\svector x}$.

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

  $$1/v( {\svector x}) \equiv \Bigl \Vert {\svector {\nabla}}_{{\svector {\scriptstyle x}}} T [ {\svector x}, {\svector x}' ] \Bigr \Vert$$ (5)

The Eikonal equation is accompanied by transport equations, which specify the geometric changes in amplitude $R ({\svector x}, {{\svector x}_s} )$. 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 ${\svector x}'$ we also assume isotropy.

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

  $$u (t, {{\svector x}_s} , {{\svector x}_r} ) = \int d {\svector x}... ...}, {{\svector x}_s} ) ~ R( {\svector x}, {{\svector x}_r} ) ~ \mbox{gain}(t) .$$ (6)

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

  $$\hat m ( {{\svector x}_b} ) = \sum_{s,r} \int dt~ \dot u (t, {{\... ...\svector x}_s} ) ~ R( {{\svector x}_b} , {{\svector x}_r} ) ~ \mbox{gain}(t) .$$ (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 ${{\svector x}_c} \equiv ( {{\svector x}_r} + {{\svector x}_s} )/2$ rather than source position. An image at a constant “half offset” ${{\svector x}_h} \equiv ({{\svector x}_r} - {{\svector x}_s} )/2$ restricts the summation to source and receivers with a constant separation:

$$\hat m_h ( {{\svector x}_b} ) = \sum_{c} \int dt~ \dot u (t, {{\svector x}_c} - {{\svector x}_h} ... ...x}_h} ) - T ({{\svector x}_b} , {{\svector x}_c} + {{\svector x}_h} ) ] ~ \cdot$$

$$\cdot ~ R( {{\svector x}_b} , {{\svector x}_c} -{{\svector x}_h} ) ~ R( {{\svector x}_b} , {{\svector x}_c} + {{\svector x}_h} ) ~ \mbox{gain}(t) .$$ (8)

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

$$\hat u (t, {{\svector x}_c} - {{\svector x}_h} , {{\svector x}_c} + {{\svector x}_h} ) = \sum_{b} \hat m_h ( {{\svector x}_b} ) ~ \delta [ t - T ( {{\svec... ...x}_h} ) - T ({{\svector x}_b} , {{\svector x}_c} + {{\svector x}_h} ) ] ~ \cdot$$

$$\cdot ~ R( {{\svector x}_b} , {{\svector x}_c} -{{\svector x}_h} ) ~ R( {{\svector x}_b} , {{\svector x}_c} + {{\svector x}_h} ) ~ \mbox{gain}(t) .$$ (9)

When the traveltime table is consistent with the data, the constant-offset images $\hat m_h ({\svector x}_b)$ should not show changes in phase over different offsets ${{\svector x}_h}$. For geometric discussions of phase delays, we can ignore the smoothly varying gain and geometric scale factors.

Reconstructing raypaths from traveltimes

Usually, one constructs a traveltime table $T$ 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 ${\svector p}$ by treating the traveltime table $T$ as a scalar potential field:

  $${\svector p}[ {\svector x}, {\svector x}' ] \equiv {\svector {\nabla}}_{{\svector {\scriptstyle x}}} T[ {\svector x}, {\svector x}' ]$$ (10)

and

  $$T[ {\svector x}, {\svector x}' ] = \int_{{\svector {\scriptstyle ... ... x}}} {\svector p}[ {\svector x}'' , {\svector x}' ] \cdot d {\svector x}'' ,$$ (11)

where the line integral can follow any path. By construction, the vector slowness is irrotational (waves should not travel in a loop): ${\svector {\nabla}}\times {\svector p}= {\svector 0}$.

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

  $$P ({\svector x}) \equiv \Bigl \Vert {\svector p}[ {\svector x}, {\svector x}' ] \Bigr \Vert$$ (12)

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 ${\svector x}(a)$ be parameterized as a function of spatial distance $a$, so that $\bigl \Vert d {\svector x}(a) / da \bigr \Vert \equiv 1$. The raypath should also be be tangent to any slowness vector that originates from another point on the path:

  $${d \over {da}} {\svector x}(a) \equiv {\svector p}[ {\svector x}(a) , {\svector x}( a_0 ) ] / P [ {\svector x}(a) ] .$$ (13)

Thus,

$$T [ {\svector x}(a) , {\svector x}(a_0 ) ] = \int_{a_0}^{a} da' ~ { d \over {da'} } T [ {\svector x}( a' ) , {\svector x}( a_0 ) ]$$

$$= \int_{a_0}^{a} da' ~ {\svector {\nabla}}T [ {\svector x}(a') , {\svector x}( a_0 ) ] \cdot { d \over {da'} } {\svector x}( a' )$$

$$= \int_{a_0}^{a} da' ~ \Bigl \Vert {\svector p}[ {\svector x}(a')... ... ) ] \Bigr \Vert ~ \Bigl \Vert { d \over {da'} } {\svector x}( a' ) \Bigr \Vert$$

$$= \int_{a_0}^{a} da' ~ P [ {\svector x}(a') ]$$ (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 ${\svector x}$ (minimum traveltime). The calculus of variations allows us to reverse the derivation.

To extrapolate a raypath from a point ${\svector x}(a_0)$ in a known direction ${\svector p}(a_0)$, we can use equation (13) and the following, which derives from (13) and (12):

  $${d \over da} {\svector p}[ {\svector x}(a) , {\svector x}(a_0) ] = {\svector {\nabla}}_{{\svector {\scriptstyle x}}} P [ {\svector x}(a) ] .$$ (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).

Residual geometric modeling and migration

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 ${\svector x}_{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 ${\svector q}_{bh}$ that is normal to the migrated reflector. For convenience, assume a unit magnitude: $\Vert {\svector q}_{bh} \Vert \equiv 1$. Locally, the coherence of this reflection could be approximated to first order as a planar surface:

  $$\hat m_h ( {\svector x}_{b} ) \approx f [ ( {\svector x}_b - {\svector x}_{bh} ) \cdot {\svector q}_{bh} ]$$ (16)

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 $\{{\svector x}_{bh}$, ${\svector q}_{bh}\}$, for many $b$ and $h$. Migrated reflectors need only be continuous enough over ${\svector x}_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 (9) remodels the data. The migrated reflection point ${\svector x}_{bh}$ contributes to all source and receiver pairs with fixed “half offsets” ${{\svector x}_h} = ({\svector x}_r - {\svector x}_s)/2$. For each affected midpoint ${\svector x}_c = ({\svector x}_r + {\svector x}_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 ${\svector p}( {\svector x}_{bh} , {\svector x}_c + {{\svector x}_h} )$ and ${\svector p}( {\svector x}_{bh} , {\svector x}_c - {{\svector x}_h} )$.

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

  $$\max_{{\svector {\scriptstyle x}}_c} ~ \Bigl \Vert [{\svector p}(... ... , {{\svector x}_c} + {{\svector x}_h} )] \cdot {\svector q}_{bh} \Bigr \Vert$$ (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

$${\svector p}( {\svector x}_{bh} , {{\svector x}_c} - {{\svector x... ... ) = 2 ~ P ({\svector x}_{bh}) ~ \cos ( \theta_{cbh} / 2) ~ {\svector q}_{bh} ,$$

$$\mbox { where } \cos ( \theta_{cbh} ) \equiv [ {\svector p}( {\sv... ...x}_{bh} , {{\svector x}_c} + {{\svector x}_h} ) ] / P ( {\svector x}_{bh} )^2 .$$ (18)

$\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

  $$t_{ch} ({\svector x}_{bh}) = T({\svector x}_{bh} , {{\svector x}_... ...vector x}_h} ) + T({\svector x}_{bh} , {{\svector x}_c} + {{\svector x}_h} ) .$$ (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 ${{\svector x}_h}$ and midpoint ${{\svector x}_c}$ we pick a traveltime $t_{ch}$. According to the migration equation (8), this pick affects all migrated positions ${\svector x}_{bh}$ 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 ${\svector p}_{ch} = {\svector {\nabla}}_{{\svector {\scriptstyle x}}_c} t_{ch}$ where ${\svector {\nabla}}_{{\svector {\scriptstyle x}}_c} \equiv ({\svector {\nabla}}... ... {\scriptstyle x}}_s} + {\svector {\nabla}}_{{\svector {\scriptstyle x}}_r})/2$. 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

  $$\max_{{\svector {\scriptstyle x}}_{bh}} ~ \Bigl \Vert [ {\svecto... ... {{\svector x}_h} , {\svector x}_{bh} ) ] \cdot {\svector p}_{ch} \Bigr \Vert$$ (20)

Thus, a constant-offset time pick $\{t_{ch},{\svector p}_{ch}\}$ or migrated pick $\{{\svector x}_{bh}, {\svector q}_{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 (3) in the main text to abbreviate $t_{ch} ({\svector x}_{bh})$. The stationary phase approximations make the same high frequency assumptions as the Eikonal and ray equations, and all fail in similar situations.

Converting time errors to reflector errors

A perturbation of the reflection point will perturb the reflection traveltime (19) according to

  $$\Delta t_{ch} ({\svector x}_{bh}) \equiv t_{ch} ({\svector x}_{b... ...) \cos ( \theta_{cbh} / 2 ) {\svector q}_{bh} \cdot \Delta {\svector x}_{bh} .$$ (21)

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 modeled and measured traveltimes. Instead, we can convert traveltime errors into equivalent errors in reflector positions:

  $$\Delta \hat {\svector x}_{cbh} = { \Delta t_{ch} ({\svector x}_{b... ...er 2 P ({\svector x}_{bh} ) \cos ( \theta_{cbh} / 2 ) } ~ {\svector q}_{bh} .$$ (22)

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

  $$e_{cbh} \equiv { \Delta t_{ch} ({\svector x}_{bh}) \over \cos ( \theta_{cbh} / 2 ) } .$$ (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).

FIGURES

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.

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.

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.