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

Workshop test problem: Travel times in the INRIA Marmousi models

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 .