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 travel-time calculation method, we should consider
in the first place the kind of travel times to calculate.
In particular, the first-arrival travel times (often but not
always required in seismic travel-time tomography) should strictly
be distinguished from ray-theory travel times (required, for
instance, in ray-theory 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 travel-time 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 travel-time 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 (CAD-like) 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 travel-time accuracy than of the order h
is required.
Calculation of the first-arrival travel times
In recent years, following the extension of random-access
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 full-wave finite difference methods.
As a rule, the grid travel-time calculating methods have been
proposed for first-arrival travel times.
The accuracy of calculated travel times depends on grid interval h .
We refer the travel-time calculation or wavefield-modelling method
to the method of the n-th order if the error of
travel time is proportional to the n-th 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 3-D structures.
Network shortest-path ray tracing
Network shortest-path ray tracing methods, based on Fermat's principle
and the graph theory, are designed to calculate first-arrival 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 self-adaptable 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 3-D travel-time tomography and the method
should thus be supplemented with a bending method
(Moser, Nolet & Snieder 1992).
Another disadvantage is the computation time
increasing in 3-D 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 travel-time tracing of the first and second orders
Under grid travel-time 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 "finite-difference"
eikonal-equation solvers.
The accuracy of the grid travel-time tracing methods considerably
increases when replacing the plane wavefront approximation (first-order
methods, travel-time error proportional to grid interval h )
by the second-order methods (travel-time error proportional
to h¹) or even higher-order 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.
Non-causal second-order (in smooth media) finite-difference
scheme for travel times was proposed by Vidale (1988),
and was made causal by several other authors, hopefully preserving
the accuracy. The first-order method calculating slowness vectors
at gridpoints together with travel times
was proposed by Van Trier & Symes (1991), the corresponding second-order
methods (in smooth media) were proposed by Pusey & Vidale (1991)
and by Klimes (1996). The second-order methods
may likely be generalized for smooth block models using
travel-time matching at structural interfaces.
Very reliable methods applicable to general 2-D and 3-D cell models
were proposed by Podvin & Lecomte (1991) (first-order method),
and by Ditmar (1995) (probably second-order 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 3-D.
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 travel-time tracing algorithms usually
increase in 3-D with the minus third power of grid step h.
Calculation of ray-theory 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
initial-value ray tracing algorithm.
The principal problem of two-point 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 two-point rays, situated inside the demarcation
belts, are lost (Bulant 1996). On the other hand,
determination of a two-point ray within a single ray history is
usually a trivial task, although more frequently addressed in geophysical
literature than demarcating the ray histories.
Controlled initial-value ray tracing
and travel-time interpolation within ray cells
Controlled initial-value 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 take-off
parameters of rays of the same history play the principal role in
controlled initial-value ray tracing. In this way, controlled
initial-value ray tracing should already be included in a good shooting
algorithm.
Once the system of rays is triangularized during controlled
initial-value 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
Initial-value 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 initial-value 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 initial-value
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 initial-value 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 ray-theory travel times at gridpoints
of very dense rectangular grids of points for, e.g., ray-theory
based migrations. The interpolation is, of course, performed not
only between depth points, but also between surface points.
Examples of the first-arrival travel times and ray-theory multi-valued
travel times calculated in both the "smooth" and "hard"
INRIA bench-mark 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 3-D Structures".
References
- Bulant, P. (1996):
- Two-point ray tracing in 3-D.
PAGEOPH, 148, 421-447.
- Cerveny, V., Klimes, L. & Psencik, I. (1988):
-
Complete seismic-ray tracing in three-dimensional structures.
In: Doornbos, D.J. (ed.), Seismological Algorithms, pp. 89-168.
Academic Press, New York.
- Ditmar, P.G. (1995):
- An approach to 2-D ray tracing intended for non-linear seismic tomography.
In: Lavrentiev, M.M., ed.: Computerized Tomography, pp. 138-154.
VSP, Utrecht.
- Gjoystdal, H., Reindhardsen, J.E. & Astebol, K. (1985):
- Computer representation of complex 3-D geological structures using
a new "solid modeling" technique.
Geophys. Prospecting, 33, 1195-1211.
- Klimes, L. (1994):
-
Kinematic hypocentre determination using the paraxial ray approximation.
Stud. geophys. geod., 38, 140-156.
- Klimes, L. & Kvasnicka, M. (1994):
-
3-D network ray tracing.
Geophys. J. int., 116, 726-738.
- Klimes, L., Kvasnicka, M. & Cerveny, V. (1994):
-
Grid computations of rays and travel times.
In: Seismic Waves in Complex 3-D Structures. Report 1., pp. 85-114,
Department of Geophysics, Charles University, Prague.
- Klimes, L. (1996):
-
Grid travel-time tracing:
second-order method for the first arrivals in smooth media.
PAGEOPH, 148, 539-563.
- Moser, T-J. (1991):
- Shortest path calculation of seismic rays.
Geophysics, 56, 59-67.
- Moser, T-J., Nolet, G. & Snieder, R. (1992):
- Ray bending revisited.
Bull. Seism. Soc. Am., 82, 259-288.
- Nakanishi, I. & Yamaguchi, K. (1986):
- A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure.
J. Phys. Earth, 34, 195-201.
- 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, 271-284.
- Pusey, L.C. & Vidale, J.E. (1991):
- Accurate finite-difference calculation of WKBJ traveltimes and amplitudes.
Expanded Abstracts, 61st Annual Int. SEG Meeting, pp. 1513-1516, Soc. Exploration Geophysicists, Tulsa.
- Van Trier, J. & Symes, W.W. (1991):
- Upwind finite-difference calculation of traveltimes.
Geophysics, 56, 812-821.
- Vidale, J.E. (1988):
- Finite-difference calculation of travel times.
Bull. seismol. Soc. Am., 78, 2062-2076.
- Vinje, V., Iversen, E. & Gjoystdal, H. (1993):
- Traveltime and amplitude estimation using wavefront construction.
Geophysics, 58, 1157-1166.
Presented at the workshop
"Computation of Multi-Valued Traveltimes"
held in INRIA Rocquencourt, France, on September 16-18, 1996.
SW3D
- main page of consortium Seismic Waves in Complex 3-D Structures .