Travel times and their computation
Ludek Klimes
The selection of the numerical method to calculate travel times
depends on the nature of the travel times, on the complexity
and computer representation of the seismic model of the geologic structure,
on the geometry of the set of points at which the
travel time is to be calculated, and on the accuracy required.
Type of travel times required by the application
When selecting the traveltime calculation method, we should consider
in the first place the kind of travel times to calculate.
In particular, the firstarrival travel times (often but not
always required in seismic traveltime tomography) should strictly
be distinguished from raytheory travel times (required, for
instance, in raytheory based migrations).
Model
Seismic model of the medium is, as a rule, specified by a finite set
of values (model parameters), and of the rules or procedures expressing
the dependence of spatial variations of material properties on these
values.
Since the exact material properties are fractal in the nature,
the material properties described by the model are, as a rule,
smoothed approximations of the exact ones. Thus the smoothness
is a natural property of a seismic model.
Moreover, probably all contemporary numerical methods of wavefield
or traveltime calculation require in some sense smooth models.
The calculation algorithm must be conformal with the selected model
specification, and the model specification should be selected
in such a way to enable application of most methods!
There is a variety of possible model specifications, for example
smooth model, gridded smooth model, cell model; smooth interfaces,
triangularized smooth interfaces, interfaces composed of triangles,
etc. Some of these representations may be converted into other ones
within required accuracy, whereas other conversions are very cumbersome.
It is thus advisable to carefully select primary model representation
in such a way to allow for conversion to other model representations
(secondary model representations) suitable for all important
contemporary traveltime calculation and wavefield modelling methods,
and also to maximize the possibility of convenient application of
unknown future methods.
After selecting the model parametrization, one should avoid to
calculate travel times for different model, although described
in terms of the same numerical values. For instance, gridded smooth
models should never be confused with cell models, but sometimes may be
substituted if the error of the substitution is correctly accounted for.
There are sevaral kinds of continuous model representations.
For instance, geological blocks may be described in terms of
volume modelling
(Gjoystdal, Reindhardsen & Astebol 1985,
Cerveny, Klimes & Psencik, 1988) or
surface modelling (CADlike) techniques.
Structural interfaces and material parameters (like propagation
velocities) may be represented by means of explicit,
implicit, or parametric representations of various degree
of smoothness.
For methods to calculate travel times (and also for methods to
calculate complete wavefields) on a rectangular grid of points,
the same grid slowness values may represent different models.
The most important kinds of a seismic model represented
in terms of grid values are
a cell model composed of rectangular homogeneous cells,
and gridded smooth model composed of discrete slowness values
at gridpoints, corresponding to a smooth model.
The smooth model may be a smooth model without interfaces,
or a smooth block model.
In the smooth model without interfaces, the slowness is
a smooth function of spatial coordinates.
The smooth block model is composed of blocks of arbitrary shapes
separated by structural interfaces covered by a finite number
of smooth surfaces. Inside each block, the slowness is
a smooth function of spatial coordinates.
Since the uncertainty of the interface position in
a gridded smooth block model is of the order of grid step h ,
the gridded smooth block models should be supplemented by gridded
block indices and by gridded information on the position of structural
interfaces if a better traveltime accuracy than of the order h
is required.
Calculation of the firstarrival travel times
In recent years, following the extension of randomaccess
computer memory (RAM), many authors proposed the methods of
calculating travel times on regular rectangular grids of points.
The grids represent the discretized forms of seismic models,
like in fullwave finite difference methods.
As a rule, the grid traveltime calculating methods have been
proposed for firstarrival travel times.
The accuracy of calculated travel times depends on grid interval h .
We refer the traveltime calculation or wavefieldmodelling method
to the method of the nth order if the error of
travel time is proportional to the nth power of h for sufficiently small grid
steps h . Because of memory limitations of present computers,
the methods of the first or smaller orders do not usually meet
the accuracy requirements of the contemporary geophysics
in inhomogeneous 3D structures.
Network shortestpath ray tracing
Network shortestpath ray tracing methods, based on Fermat's principle
and the graph theory, are designed to calculate firstarrival travel times.
They developed from the less accurate methods of the orders
between 0 and 2/3
by Nakanishi & Yamaguchi (1986) and Moser (1991),
to the selfadaptable algorithm
of the first order by Klimes & Kvasnicka (1994).
The method is very reliable but the accuracy of the order h
is not sufficient for 3D traveltime tomography and the method
should thus be supplemented with a bending method
(Moser, Nolet & Snieder 1992).
Another disadvantage is the computation time
increasing in 3D with minus fourth power of grid step h .
Let us emphasize that bending is not a complete method
of ray tracing. It may be especially valuable to postprocess
preliminary approximations of the rays estimated by means of
another method. The accuracy of preliminary ray has to be sufficient
in order to distinguish the required ray strictly from other extremal
paths.
Grid traveltime tracing of the first and second orders
Under grid traveltime tracing we understand here
the methods interpolating already calculated travel times between
gridpoints in order to evaluate the travel time in the closest next
gridpoint. These methods are often called "finitedifference"
eikonalequation solvers.
The accuracy of the grid traveltime tracing methods considerably
increases when replacing the plane wavefront approximation (firstorder
methods, traveltime error proportional to grid interval h )
by the secondorder methods (traveltime error proportional
to h¹) or even higherorder methods, if the updating scheme
is made stable for the interpolation method. It is also useful to locally
apply spherical wavefront approximation (but not to use global spherical
coordinates!).
Methods designed for gridded smooth models may be divided
into methods directly extrapolating travel times,
and to the methods calculating slowness vectors.
Noncausal secondorder (in smooth media) finitedifference
scheme for travel times was proposed by Vidale (1988),
and was made causal by several other authors, hopefully preserving
the accuracy. The firstorder method calculating slowness vectors
at gridpoints together with travel times
was proposed by Van Trier & Symes (1991), the corresponding secondorder
methods (in smooth media) were proposed by Pusey & Vidale (1991)
and by Klimes (1996). The secondorder methods
may likely be generalized for smooth block models using
traveltime matching at structural interfaces.
Very reliable methods applicable to general 2D and 3D cell models
were proposed by Podvin & Lecomte (1991) (firstorder method),
and by Ditmar (1995) (probably secondorder accuracy in cell models),
using direct extrapolation of travel times.
To enable quadratic interpolation of travel times, Ditmar calculates
travel times, in addition to gridpoints at the corners of cells,
also in additional points situated at the gridlines in the middle
between the pairs of neighbouring gridpoints. Thus the number of
calculated travel times is four times the number of gridpoints in 3D.
Since designed for cell models, Ditmar's (1995) algorithm is not applicable
to gridded smooth models, especially to models with smooth interfaces,
without a loss of accuracy.
The computation time of the grid traveltime tracing algorithms usually
increase in 3D with the minus third power of grid step h.
Calculation of raytheory travel times
Wavefront tracing
In the wavefront tracing method, complete wavefronts are discretized
(triangularized)
with a sufficient density and traced along short ray elements from
preceding to succeeding wavefront, see Vinje, Iversen & Gjoystdal (1993).
The density of rays is continuously adjusted by means of interpolation.
In this way the space is divided into the ray cells and travel time
may be very accurately interpolated to the receivers or gridpoints inside
individual ray cells.
Note that wavefront tracing may also be modified for first arrivals.
Shooting
Shooting method is based on application of some sort of
initialvalue ray tracing algorithm.
The principal problem of twopoint ray tracing in heterogeneous
block models consists in a sufficiently accurate determination
of boundaries between rays with different ray histories, i.e. passing
through different geological blocks, across different structural
interfaces, terminated in different areas for different reasons,
or passing through a different number of caustic points.
The demarcation belts between the different ray histories, corresponding
to numerical uncertainty of the boundaries, have to be kept reasonably
narrow, because all twopoint rays, situated inside the demarcation
belts, are lost (Bulant 1996). On the other hand,
determination of a twopoint ray within a single ray history is
usually a trivial task, although more frequently addressed in geophysical
literature than demarcating the ray histories.
Controlled initialvalue ray tracing
and traveltime interpolation within ray cells
Controlled initialvalue ray tracing is closely related to shooting.
Its aim is not to hit given receivers, but to cover the model with
a sufficiently dense system of rays (or to cover the given surface with
a sufficiently dense system of ray endpoints).
The demarcation of ray histories and the triangularization of takeoff
parameters of rays of the same history play the principal role in
controlled initialvalue ray tracing. In this way, controlled
initialvalue ray tracing should already be included in a good shooting
algorithm.
Once the system of rays is triangularized during controlled
initialvalue ray tracing,
the space is divided into the ray cells and travel time
may be very accurately interpolated to the receivers or gridpoints
inside individual ray cells, like during wavefront tracing
This method is in principle equivalent to wavefront tracing.
Weighting of paraxial ray approximations
Initialvalue ray tracing may yield the points along all traced rays,
stored with a given step. Unfortunately, the points are distributed
irregularly in space and the travel times are multivalued that prevents
application of routine interpolation techniques.
If the system of rays is sufficiently dense and the step along rays
is reasonable, the travel times and all other calculated quantities
may be interpolated by weighting paraxial ray approximations.
The method has been proposed by Klimes (1994) for nonlinear
hypocentre determination, and its accuracy and applicability has
been demonstrated by Klimes, Kvasnicka & Cerveny (1994).
Since the initialvalue ray tracing itself is not capable to guarantee
a sufficient coverage of the target zone, in which the grid is situated,
by rays, this method is unsafe unless controlled initialvalue
ray tracing or wavefront tracing is used. However, in such a case
the interpolation within individual ray cells is usually preferable.
Interpolation of travel times in rectangular grids of points
The interpolation in a rectangular grid of points may be much
faster than interpolation within ray cells or than weighting of
paraxial ray approximations.
That is why wavefront tracing, controlled initialvalue ray tracing,
or weighting of paraxial ray approximations should be supplemented
with several steps of consecutive interpolation in order to considerably
speed up the calculation of raytheory travel times at gridpoints
of very dense rectangular grids of points for, e.g., raytheory
based migrations. The interpolation is, of course, performed not
only between depth points, but also between surface points.
Examples of the firstarrival travel times and raytheory multivalued
travel times calculated in both the "smooth" and "hard"
INRIA benchmark versions of the Marmousi model.
Acknowledgements
The research has been supported
by the Grant Agency of the Czech Republic under Contract 205/95/1465,
and by the members of
the consortium "Seismic Waves in Complex 3D Structures".
References
 Bulant, P. (1996):
 Twopoint ray tracing in 3D.
PAGEOPH, 148, 421447.
 Cerveny, V., Klimes, L. & Psencik, I. (1988):

Complete seismicray tracing in threedimensional structures.
In: Doornbos, D.J. (ed.), Seismological Algorithms, pp. 89168.
Academic Press, New York.
 Ditmar, P.G. (1995):
 An approach to 2D ray tracing intended for nonlinear seismic tomography.
In: Lavrentiev, M.M., ed.: Computerized Tomography, pp. 138154.
VSP, Utrecht.
 Gjoystdal, H., Reindhardsen, J.E. & Astebol, K. (1985):
 Computer representation of complex 3D geological structures using
a new "solid modeling" technique.
Geophys. Prospecting, 33, 11951211.
 Klimes, L. (1994):

Kinematic hypocentre determination using the paraxial ray approximation.
Stud. geophys. geod., 38, 140156.
 Klimes, L. & Kvasnicka, M. (1994):

3D network ray tracing.
Geophys. J. int., 116, 726738.
 Klimes, L., Kvasnicka, M. & Cerveny, V. (1994):

Grid computations of rays and travel times.
In: Seismic Waves in Complex 3D Structures. Report 1., pp. 85114,
Department of Geophysics, Charles University, Prague.
 Klimes, L. (1996):

Grid traveltime tracing:
secondorder method for the first arrivals in smooth media.
PAGEOPH, 148, 539563.
 Moser, TJ. (1991):
 Shortest path calculation of seismic rays.
Geophysics, 56, 5967.
 Moser, TJ., Nolet, G. & Snieder, R. (1992):
 Ray bending revisited.
Bull. Seism. Soc. Am., 82, 259288.
 Nakanishi, I. & Yamaguchi, K. (1986):
 A numerical experiment on nonlinear image reconstruction from firstarrival times for twodimensional island arc structure.
J. Phys. Earth, 34, 195201.
 Podvin, P. & Lecomte, I. (1991):
 Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools.
Geophys. J. int., 105, 271284.
 Pusey, L.C. & Vidale, J.E. (1991):
 Accurate finitedifference calculation of WKBJ traveltimes and amplitudes.
Expanded Abstracts, 61st Annual Int. SEG Meeting, pp. 15131516, Soc. Exploration Geophysicists, Tulsa.
 Van Trier, J. & Symes, W.W. (1991):
 Upwind finitedifference calculation of traveltimes.
Geophysics, 56, 812821.
 Vidale, J.E. (1988):
 Finitedifference calculation of travel times.
Bull. seismol. Soc. Am., 78, 20622076.
 Vinje, V., Iversen, E. & Gjoystdal, H. (1993):
 Traveltime and amplitude estimation using wavefront construction.
Geophysics, 58, 11571166.
Presented at the workshop
"Computation of MultiValued Traveltimes"
held in INRIA Rocquencourt, France, on September 1618, 1996.
SW3D
 main page of consortium Seismic Waves in Complex 3D Structures .