Review of traveltime computation methods

F. Audebert, D. Bevc, B. Biondi, D. Lumley, D. Nichols, G. Palacharla, T. Rekdal & H. Urdaneta

Author has no known email address

ABSTRACT

We overview a wide range of traveltime computation methods, most of them novelties from the last five years or so. We first classify the methods, then we give a standardized description. This unified form allows a convenient comparison of these methods. Each description first outlines the basics principles of the given method. It then identifies its main features: assumptions made on the velocity model, strong points, and weak points. Finally, we attempt a quantitative rating of all the methods. The rating is based on geophysical and efficiency criteria.

The last five years or so have witnessed a bloom of novel traveltime computation methods. These methods provide parameterized estimates of Green's function between a given source and all depth points in a 3-D volume. These Green's function estimates are in the form either of single arrival traveltime, or more sophisticated traveltime, amplitude and phase estimates. The Green's functions can be used as input to Kirchhoff-type migration methods. Thus, given the possible importance of kinematic Kirchhoff-type methods for 3D pre-stack processing, it seems necessary to have a clear view of the state of the art. Thanks to both novelty and diversity, it is dubious that any single individual has an up-to-date, clear perception of all that has been done recently in this domain. We think it is important, both for our future research work, and for the use of our sponsors, to create this overview.

TENTATIVE CLASSIFICATION

The next section presents a list of methods we have reviewed. One clear cut division is between infinite bandwidth family and finite bandwidth family. That is to say between methods explicitly working in the seismic frequency range (finite bandwidth) and methods assuming either implicitly or explicitly some infinite frequency or infinite bandwidth condition (ray theory, etc...).

FAMILY: Limited Bandwidth methods

These limited bandwidth methods are typically the extrapolation or modeling algorithms working in the time domain or the temporal frequency domain: phase-shift or finite-difference downward continuation, etc...  . These algorithms deal with frequency-dependent, complex-valued Green's functions. The traveltimes do not exist explicitly. The reference is a finite frequency-band wavelet propagated from the source to the depth points using a two-way wave equation. The modeled wavefield recorded at a given depth point contains the interference of many arrivals. They include direct arrivals, head waves, diffractions, and reflected energy. A simpler wavefield can be modeled using either a one-way wave equation () or a non-reflecting wave equation (). Again a full wavefield is modeled at the depth points, but now the events are all direct arrivals. This kind of modeling is commonly used in finite-difference shot profile migration, for the propagation of the model source wavelet.

The method of band-limited traveltimes is referred to here as Nichols' method. It is described in this report () and a full description is provided in his thesis (). It is a reduced version of a full wavefield modeling. The innovation is that only a decimated set of frequencies are processed. Decimating the frequencies means that the computational cost is decimated in proportion, but at the expense of the modeled wavefield, which becomes aliased. The wavefield is correctly modeled only in a small time window (around the maximum energy arrival). Nevertheless, by keeping track of this time window, it is possible to pick either first or most energetic arrivals within the window. The corresponding amplitudes and phases are also estimated. Thus, at a cost lower than that of full wavefield modeling (though still higher than the fastest ray-tracing methods) a reduced traveltime+amplitude+phase output is produced. This reduced (with respect to a full wave field) output is similar to the output of ray tracing, but it is expected to be more reliable. This method bridges the gap between full wavefield modeling and ray tracing. It is truly a traveltime computation method, but it takes into account the frequency-dependence of propagation.

FAMILY: Infinite frequency methods

The family of infinite frequency methods is best represented by the classic ray-theory methods. There are followed by more and more distant relatives, of which the only common feature is to be theoretically valid in the high frequency limit.

Genus: Ray based methods

The first `` genus '' of this family includes the classic ray tracing methods, paraxial and dynamic ray tracing (, , ), plus some recent methods which use ray tracing to mimic wavefront propagation, or, even better, frequency dependent wavefront propagation.

However classic, ray tracing is still evolving. For example, some recent improvements, described in this report (), try to take into account merging of multiple events that arrive within a wavelength.

Gaussian beam modeling and migration represent a first attempt to make ray tracing more robust. Those methods are classic and well described by various authors (, ). Therefore, we will not describe them in detail in this paper.

More recently, we find wavefront propagation by (local) ray tracing: the best example is the wavefront reconstruction method (). In this paper we refer to it as the ``NORSAR'' method (only because the name of the research institute is more convenient than the long list of the authors). This method overcomes the problem of interpolating wildly diverging rays, simply by interpolating before the divergence of the rays reaches a certain level. Moreover, it allows the tracking of multiple arrivals, and corresponding amplitudes. Nevertheless, it simulates only infinite frequency wavefronts, since it is based on ray-tracing.

Lin Zhang's method, described later in this paper, belongs to this category, because it uses local ray tracing.

Hybrid Genus: Rays with finite-frequency behavior Ray tracing methods are based on the solution of the eikonal equation, that is a high-frequency approximation to the scalar wave-equation. When the velocity model is rough compared with the wavelengths of the seismic signal, the solution of the eikonal equation is not necessarily close to the solution of the wave equation with a band-limited source. Our numerical experiments with Marmousi () show that the differences between the ray-tracing solution of the eikonal equation and the band-limited wave-equation solution fall in three main categories. The most dramatic misbehavior is instability of the ray-tracing solution; when it goes undetected it can seriously deteriorate the final image. Fortunately, safeguards against numerical instability are available for most of the methods. The case when ray-tracing computes a solution that is a valid solution to the eikonal equation, but is very different from the desired band-limited solution (e. g. waves trapped in low velocity layer), is more insidious, since less easy to detect, than numerical instability. Finally, the most common situation is when the ray-tracing solution is close to the band-limited one, but the traveltimes, phases, and amplitudes are slightly different. In general, amplitude errors are larger and more frequent than traveltime errors. Obviously amplitude errors will seriously undermine the reliability of AVO analysis; but also the kinematics of migration results are affected, when only the most energetic arrivals are used during the summation step. Wrong amplitudes can often lead to the selection of arrivals that are apparently high-energy, but are actually low energy in the band-limited solution.

Smoothing the slowness function is the most common solution in dealing with the problems caused by the high-frequency assumption. There are few heuristic justifications for it. An obvious one is that a smoother model is closer to fulfill the validity condition for ray tracing (). Another justification is that waves are influenced by the velocity function in the whole Fresnel region (, ), and not only along the raypath. Unfortunately there are only few and limited attempts () to directly link velocity smoothing with the solution of the band-limited wave-equation. Very few guidelines exist to determine quantitatively the amount of smoothing required, which is even more important for applications. The simplest, and most commonly used smoothing is the convolution of the slowness function with a constant-width isotropic averaging function e.g. a Gaussian. A slightly more sophisticated strategy is to adapt the width of the smoothing function to the local wavelength at an assigned central frequency. Our numerical experiments show that this adaptive approach gives better results than constant smoothing. Finally, the most complex, but the most promising, approaches, are the ones in which smoothing is dependent on the raypaths. The slowness function is smoothed perpendicular to the rays with a smoothing function as wide as the local wavelength. The simplest of these methods is when the smoothing is applied a posteriori and does not influence the raypath (). On the contrary, in the Lomax method reviewed in this paper the velocity is smoothed during ray tracing, and rays bend according to the smoothed velocity function.

Genus: Expanding computation fronts

This genus is hard to define or even identify as a genus. Its members has the common characteristics that they cumulate traveltime from the source while computing incremental traveltime on a computational front. This expanding computation front somehow mimics an wavefront in a given slowness grid. All these methods invoke the necessary ingredients: eikonal equation, Huygens' and Fermat's principles.

Vidale's method for a finite difference solution to the eikonal equation () is historically the first member of the genus. The eikonal equation is solved by finite-differences on a Cartesian grid (constant slowness rectangular cells). The (incremental) traveltimes are calculated on an ever-expanding rectangular front. There is no relaxation, and no updating back into the inside of the computation front.

The upwind finite-difference method () is an improvement of Vidale's method. A modified eikonal, allowing for the direct computation of first arrival, is solved by finite differences in polar coordinates (constant slowness cells). The expanding front is a semi-circle, expected to better follow an actual wavefront than Vidale's square front.

Podvin's method () falls into the present genus mostly because it resembles some of the previous methods, but the expanding wavefront criterion is only a preferred implementation, not a basis of the method. In fact, this method can as well be used in a general relaxation mode. In many aspects, it resembles Vidale's method: it uses constant slowness Cartesian cells, and, in the preferred implementation, a rectangular expanding pseudo-wavefront. The incremental traveltimes are computed locally with a simple straight rays or plane wave approximation. It is applying plain Huygens' principle in constant slowness cell. A method described by Schneider et al. is very similar.

Local ray tracing () was envisioned by its author as an improvement over Podvin's method. It is again a method that cumulates traveltime upon an expanding computation front, but this time, the incremental traveltimes are computed by local ray tracing. The slowness grid is a triangulated slowness surface: it is expected to provide a better representation of a velocity model. This method is very rich in further possibilities (amplitude, multiple arrivals), but at a corresponding cost. Since Lin Zhang's method uses local ray tracing, it can also be placed in the ray methods category.

Genus: Graph theory + Fermat

Some authors have recently proposed to apply graph theory to the computation of first arrival traveltimes. Moser for instance suggested that finding first arrival traveltimes is exactly like finding the minimum length path in a network of local traveltimes. The graph in this case is just the local traveltime (=distance*slowness) mesh.f

LIST OF REVIEWED METHODS

Limited bandwidth methods

$\bullet$

Band-limited Green's functions (Nichols' method)

Infinite frequency methods

*

Ray based methods.

$\bullet$

Paraxial or dynamic ray tracing

$\bullet$

Gaussian beams

$\bullet$

Wavefront reconstruction (``NORSAR'' method)

$\bullet$

Waveray method (Lomax's method)

*

Expanding computation fronts.

• Finite differencing the eikonal equation

  $\bullet$

  Expanding rectangular wavefront (Vidale's method)

  $\bullet$

  Upwind finite-difference (van Trier & Symes' method)

• Local ray tracing + Fermat

  $\bullet$

  Local ray tracing (Zhang's method)

• Huygens + Fermat

  $\bullet$

  Podvin's method

  $\bullet$

  Schneider et al.'s method

*

Graph theory + Fermat

$\bullet$

Shortest path rays (Moser's method)

NICHOLS' METHOD

Principle

This method provides a traveltime solution in a user-defined frequency-band rather than in a high frequency approximation (). A small number of frequencies (8-16) in the seismic frequency band are extrapolated outwards from the source location using a paraxial one-wave equation in polar coordinates. At each radius a parametric approximation to the wavefield is estimated as follows.
1) Calculate the Green's functions for a sparse frequency sampling at the new radius.
2) Choose a time window centered around the traveltime from previous radius.
3) Calculate time domain data in the window by slow Fourier transform.
4) Pick the maximum energy sample.
5) Use a quadratic fit to find the traveltime of the local peak of the energy function.
6) Calculate the amplitude, and phase at this traveltime. Features

• Velocity model

  The velocity model is a polar coordinate grid of constant slowness cells. Any velocity model can be used as input, there are no restrictions on the smoothness of the model.

• Advantages
  • Provides traveltime, amplitude and phase, not just traveltime.
  • Initial condition at the source is easily specified since the computations are done in polar coordinates.
  • The traveltime of the global maximum amplitude energy is found.
  • The solution is found at every point in the subsurface (no shadow zones)
  • The solution is an estimate of the Green's function in the seismic frequency band not the solution at very high frequency.
• Weaknesses
  • Cost, it is from 2 to 10 times as expensive as a finite-difference solution to the eikonal equation.
  • The current implementation is 2-D only. A 3-D implementation is possible but it has not yet been written.
  • The maximum amplitude traveltime field is discontinuous. This makes it harder to interpolate.
  • No explicit raypaths are calculated. If the rays or takeoff angles are required, they must be inferred from the Green's function.

PARAXIAL RAY TRACING

Principle

The paraxial ray tracing is a high-frequency ray approximation of the elastodynamic equation. P- and S-waves are considered totally decoupled. The traveltime is assumed to be independent on amplitude factors. As in all ray tracing methods, the raypaths and the traveltimes are found by solving the eikonal equation. The geometrical spreading can be computed by solving the transport equation. The curvature of the wavefront is computed from the dynamic ray equations, and is used to estimate the traveltimes (and the displacement direction) in the vicinity of the rays: the wavefront curvature is approximated with a parabola. The dynamic ray tracing equations are valid in the far-field from a point source and far from singularities in weakly inhomogeneous media. Reflection- and transmission-coefficients, Snells law, phase- and ray parameter matching has been successfully introduced to use ray tracing in layered media (). Features

• Velocity model:

  The paraxial ray tracing can handle all velocity representations in both 2-D and 3-D. Analytical solutions exist for some sets of solutions where the velocity is assumed to have some linear properties within tetrahedrons or triangles. Otherwise a numerical method, e.g. a Runge Kutta solver, can be applied.

• Advantages:

  The paraxial ray tracing method is fairly cheap in memory, and is used in 3-D. The method handles complex model structures. The method computes traveltimes very accurately. It finds multiple arrivals, amplitude and phase-shifts from caustics. The amplitude information can be used, among others, to find the most energetic arrivals.

• Weaknesses:

  The method is a high-frequency approximation, and the ray validity conditions are not always fulfilled. It produces shadow zones in the wavefield. The method is slow compared to many first arrival traveltime techniques.

GAUSSIAN BEAM METHOD

Principle

The Gaussian beam can be considered as a bundle of generally complex rays traveling in the vicinity of a central ray. The amplitude variation is Gaussian (bell-shaped) with maximum at the central ray. The type of beams, i.e., the shape of the generally elliptical bell (the beam width) and the curvature of the ``wavefront'' of the beams can to some extent be chosen freely. The superposition is an integral over beam solutions at the receiver point taken over all values of takeoff angles. The weight-function in the integrand is defined such that the asymptotic solution to the integral, according to the method of stationary phase, is equal to the ray solution in regular regions. The beam width and the curvature of the beams wavefront, are expressed by the complex Hessian matrix of second derivatives of complex phase with respect to curvilinear coordinates perpendicular to the central ray. This matrix can be initiated at different points on the ray. The simplest is to initiate the matrix at the end point of the rays. For requirements on the choice of parameter values and for details on the Gaussian beam method, see Cervený and Psencík .

Features

• Velocity model:

  The method has the same features here as the paraxial ray tracing method.

• Advantages:

  The Gaussian beam method has the same advantages as the paraxial ray tracing method, but the method gives reasonable results in caustic regions. Gaussian beam can also give reasonable results in other regions such as the transition zone between an illuminated region and a shadow-zone, but the results here are highly dependent upon the choice of the ``beam-width'' parameters. The method gives a much smoother and frequency dependent solution than the paraxial ray tracing method.

• Weaknesses:

  The Gaussian beam method is based upon the high frequency, asymptotic ray approximation, and diffraction effects are not handled correctly.

NORSAR METHOD

Principle

The NORSAR (or wavefront construction) method calculates traveltimes and amplitude coefficients of both first and later arrivals, in a smooth and continuous velocity model (). The main idea behind NORSAR method is to calculate the ray parameters along wavefronts instead of tracing independent rays. At each time step, a new wavefront is constructed from a previous one, by ray tracing. New rays are interpolated between rays that are further apart than a predefined distance. In this way, areas of large geometrical spreading where conventional ray tracing may not give arrivals, are fully covered with NORSAR method.

As the wavefronts are constructed, ray cells are checked for the presence of receivers. Ray cells are defined by a pair of contiguous wavefronts and neighboring rays. The ray attributes (travel-times, amplitude coefficients, ray direction, etc.), which are known as a function of the location along the wavefronts, are interpolated to the actual location of the receivers. The interpolation is done from the four corners of a ray cell. In a fist arrival mode, later arrivals are removed from the wavefronts in order to save computational time and memory.

Features

• Velocity model

  The NORSAR-type code that is implemented at SEP, runs on a rectangular 2-D grid of constant velocities. The NORSAR group has also implemented a 3-D version of the method. The method requires a velocity field that has a continuous first derivative.

• Advantages

  The advantages of the NORSAR method are its flexibility, robustness and accuracy. First and later arrivals may be found at any points in the model. It can handle all types of turning rays. Another advantage with respect to conventional ray tracing, is that the estimation of travel-times, amplitudes etc. does not come from a posteriori interpolation between single, separate rays, but instead comes directly from an already reconstructed wavefront.

• Weaknesses

  The only weaknesses we can see are the infinite frequency assumption (or rather, infinite bandwidth assumption) and the necessity of a C¹ slowness medium. Those weaknesses can be overcome by combining NORSAR and Lomax approaches, or by using NORSAR method in a medium smoothed with over a wavelength.

LOMAX'S METHOD

Principle

The Lomax or waveray method is an approximate method for kinematic modeling of band-limited wave propagation in heterogeneous velocity structures (). It combines Huygen's principle (to track the motion of narrow band wavefronts at a number of center frequencies) and wavelength-dependent velocity smoothing. The motion through time of the narrow band wavefronts determines wavepaths, which are frequency-dependent. The narrow band wavefronts at any time are approximated by a straight line normal to the wavepath. The wavelength-dependent velocity is obtained by using a Gaussian type weighted average along the wavefront. The wavefield then is constructed by summing over the contributions of all narrow band ``wavefields'' at all frequencies arriving at a given receiver. Traveltimes and amplitudes, for a given central frequency, are estimated as in conventional ray methods, from the traveltime and geometrical spreading of adjacent wavepaths passing near the receiver location.

Features

• Velocity model

  There are no restrictions on the model roughness, except that the medium must be smoothly varying relative to a given wavelength. It has only been implemented for 2-D models.

• Advantages

  Refracted direct waves are accurately reproduced in homogeneous or smoothly varying regions. Transmitted refractions, wide-angle reflections and headwaves are reproduced approximately at discontinuities. Frequency-dependent scattering of some wave types is reproduced as a function of the ratio of wavelength to characteristic size of scattering region. A portion of the diffracted energy is produced in geometrical shadow regions.

• Weaknesses

  It has not being formally derived from basic equations, but is the result of the application of physical principles. Just like in conventional ray tracing methods its deficiencies include the following: (1) In regions with strong velocity gradients there are no high angle reflections. (2) In critical regions where the geometrical spreading function is small or singular there may be instability in the amplitude estimates. (3) There is incomplete modeling of diffracted waves in geometrical shadow zones. (4) The use of a finite number of wavepaths can lead to poor sampling of parts of the structure and inaccurate synthesis of corresponding parts of the wavefield.

VIDALE'S METHOD

Principle

Vidale's method computes traveltimes by solving the eikonal equation on a rectangular grid using expanding wavefronts (). The calculation begins by assuming a seismic source at a grid point and calculating four adjacent traveltimes by averaging local slowness. The corner traveltimes are then calculated by finite-differencing the eikonal equation and assuming either planar or circular wavefronts. The assumption of planar wavefronts is sufficiently accurate for most smoothly varying velocity models at grid points far away from the source. Once a square grid of traveltimes is found, a non-centered finite-difference of the eikonal equation is used to expand the computation out from the edge of the square. For stability, the solution must follow causality, therefore, relative minima are found along the edges of the square and the computational front is expanded from these points. Since there can be multiple relative minima along a computational front, some care is required in expanding the traveltime front in such a way as to guarantee stability and retain first arrivals.

The method is easily extended to three dimensions (). Geometric amplitudes can be efficiently calculated using the traveltimes from four closely spaced source points, and without explicit knowledge of raypaths ().

Features

• Velocity model
  • The velocity model is represented by a rectangular grid.
  • The algorithm is not accurate for rapidly varying velocity models.
• Advantages
  • The method is stable for smooth velocity models,
  • and is straight-forward to understand and to implement.
• Weaknesses
  • Involves sorting of points on the computational front. This makes the scheme non-vectorizable and non-parallelizable.
  • Computes only the first arrival field.

VAN TRIER AND SYMES METHOD

Principle

This method computes the traveltimes by solving the eikonal equation on a polar grid using an upwind finite difference scheme (). It is an explicit scheme. The method solves a conservation law derived from the eikonal equation that describes changes in gradient components of the traveltime field. The solution gives the first arrival traveltime field. The traveltimes are extrapolated using a plane wave approximation. The scheme treats all points equally, unlike Vidale's scheme where the minimum traveltime point has to be found at each extrapolation step. The computations proceed along an expanding front, starting from an initial condition specified at the source. The method may become unstable if the radial step is too large. Many authors have implemented the CFL stability condition to perform adaptive refinement of the radial step-size. This guarantees a stable algorithm in 2-D, but the cost may rise rapidly in a complex velocity model.

Features

• Velocity model

  The velocity model is represented by a rectangular grid of constant slowness cells which are then mapped to a polar grid. The grid is 2-D. A simple bilinear interpolation is used to do the mapping. The algorithm has stability problems in three dimensions.

• Advantages
  • Does not involve any sorting of points on the computational front. This makes the scheme fully vectorizable.
  • The initial condition at the source is easily specified since the computations are done in polar coordinates.
  • The boundary conditions are easily implemented since they use the same stencil as the interior of the model.
• Weaknesses
  • Local plane wave assumption (not a limitation, only affect the precision in the vicinity of the source).
  • Computes only the first arrival traveltime field.
  • No amplitude information.
  • Fails when the time field to be computed does not have an outward-pointing gradient at each front.
  • Stability problems in three dimensions.

ZHANG'S METHOD

Principle The method of wavefront propagation by local ray tracing (Zhang, 1993) was developed to take advantage of the strengths of finite-difference eikonal solvers: computational speed and evaluation on a regular grid without shadow zones; and to overcome their inherent weaknesses: discontinuous velocity field representation, low order finite-difference accuracy, and large traveltime gradient, amplitude and ray-angle errors. The velocity field is parameterized by regular triangular cells on a hexagonal computational mesh. Each triangular cell contains a constant velocity gradient and the velocity field is continuous along cell boundaries, which allows all ray tracing calculations to be evaluated analytically. Zhang propagates five ray properties from one hexagonal ring to the next: global traveltime, local traveltime gradient, ray angle, local wavefront radius of curvature, and global geometric spreading factor. The following wavefront attributes can be computed on the full mesh independently: first arrival time (global), most energetic amplitude (local), upgoing or downgoing wavefronts (local), and left-going or right-going wavefronts (local). The method allows for turning rays which do not turn back across the expanding hexagonal computational ring.

Features

• Velocity model:

  The velocity model has a continuous representation in triangular cells of constant velocity gradient.

• Advantages:

  The accuracy is an order of magnitude better than finite-difference eikonal schemes, and the method is robust in the models of Marmousi-type complexity. The method computes several wavefront attributes: traveltime, geometric spreading amplitude, local plane-wave angle and takeoff angle, and traveltime gradient. The method is capable of selecting different traveltime branches based on first arrivals, most energetic arrivals, or propagation direction. A 3-D version exists.

• Weaknesses:

  The method is two to three times slower than vanilla finite-differencing of the eikonal (e.g., Vidale). Some extreme turning rays are not accounted for, and the traveltime branch selection option is based on local rather than global attributes. The method is not well-suited for vectorizing.

PODVIN'S METHOD

Principle

The keywords of this method are finite-difference (discretization in space of the traveltime), Huygen's principle (to compute local traveltime) and Fermat's principle (to obtain first arrival) (). Given a source location on a rectangular slowness grid, a traveltime table is initialized at t = 0 at the source location, and $t = \infty$ elsewhere. The idea is to find, through successive relaxations from neighboring nodes to neighboring nodes, the minimum traveltime from the source to every depth location. The increment of traveltime are estimated locally between the nodes of a cubic slowness grid. The incremental traveltimes between a given node and its neighbors is estimated by applying Huygen's principle along the faces of the cube surrounding this node. The cumulated traveltime is computed from the source, on an expanding computational front. At every step, increments in traveltimes are computed. Fermat's principle is then invoked, at any node and at any computational step, to keep only the absolute minimum traveltime computed.

Features

• Velocity model

  The velocity model is represented by a rectangular grid of constant slowness cells. The grid is either 2-D or 3-D.

• Advantages
  • Pondvin's method is unconditionally stable and absolutely robust.
  • There is no constraint or assumption on the velocity distribution, other than the rectangular, constant slowness cells representation.
  • The accuracy of the traveltimes depends essentially on the velocity grid resolution.
  • It always finds first-arrival traveltime. It can handle all kind of turning rays.
  • There is the possibility to compute multiple-arrivals, which is more difficult to program and also more computer expensive.
  • A 3-D version exists.
• Weaknesses
  • Local plane wave assumption (not really a limitation, it only affects the precision in the vicinity of the source).
  • No amplitude information available, excepted a posteriori estimates.
  • Not at all or little vectorizable. Only coarse grain parallelizable.
  • Very slow.

SCHNEIDER ET AL'S METHOD

Principle

This method, or family of techniques, computes the first arrival traveltimes in complex velocity model (). The velocity is constant within each cell. The source is positioned at one of the edges in the model. In the first approach to obtain the minimum traveltime, traveltimes (preliminary first arrivals) are first computed at the nodes on this edge. This is done by assuming straight ray paths along the edge. Traveltimes at the lower corner in next column is computed from the two node values at the last updated column. Either a plane wave or a circular wave assumption can be made. After all the nodes in the column has a traveltime value, a new pass is done in the other direction. This time the values from the two lower grid points are used to compute a traveltime at the upper corner in the column. The cells are passed four times in different directions, and the minimum traveltime obtained at each node is kept. The method has many similarities with Podvin's relaxation method (see section about Podvin's method). The second approach uses a mapping similar to Vidale's approach (see section about Vidale's method). The traveltimes are computed on the edge of an expanding rectangle. The rays are assumed to pass through the edge of a specific rectangle only once.

Features

• Velocity model:
  

  The methods can handle any velocity model, but the slowness is constant within rectangular cells.

• Advantages:
  

  The method computes a traveltime in every grid point. The method is efficient, and can be extended to three dimensions. The code is vectorized. The method can handle arbitrarily rough velocity models. The methods are stable and computes accurate traveltimes.

• Weaknesses:
  

  The method can only compute absolute first arrival traveltimes.

MOSER'S METHOD

Principle

The idea underlying Moser's method is to ``discretize'' the problem of finding seismic raypaths, and then, based on Fermat's principle, approximate the problem into calculating the shortest traveltime paths through a network, using graph theory. The network, which represents the model over which rays are traced, consists of nodes that are connected to neighboring nodes by arcs with a weight equal to the traveltime of a seismic wave along it. The problem of obtaining the single source shortest-paths on a network, is solved using Dijkstra's algorithm. Better performances are obtained by using a modified Dijkstra in which traveltimes are arranged into a binary heap. Two optimizations to Moser's method were presented by Fisher and Lees (1993) to obtain a faster and more accurate algorithm by perturbing shortest traveltime ``rays'' at interfaces according to Snell's law and by tracing straight rays in regions of low velocity contrast. Moser's method uses a very efficient way to calculate the raypaths and traveltimes of absolute first arrivals to all points in the model. Later arrivals, caused by reflections on interfaces or by multiples, can also be computed by posing constraints to the shortest paths. The computation time for Moser's method using the modified Dijkstra is almost linearly dependent on the number of nodes. The accuracy is quadratically dependent upon the number of points per coordinate direction and the number of arcs per node.

Features

• Velocity model

  There are no restrictions on the complexity or the dimensionality of the velocity model. The velocity model is represented by a rectangular grid of constant slowness cells.

• Advantages

  There are no restrictions of classical ray theory: diffracted raypaths and paths to shadow zones are found correctly. There are no problems with convergence of trial raypaths toward a specified receiver nor with raypaths with only a local minimal traveltime. All shortest paths from one source are constructed simultaneously. There is no notion of dimensionality of the space, the same algorithm that is used for the two-dimensional ray tracing can be used for three dimensions, just by adding more nodes. The problem of finding first arrival raypaths is framed inside graph theory, where its powerful tools can be applied.

• Weaknesses

  No amplitudes are obtained. It only computes absolute first arrival traveltimes.

Comparative rating: Geophysical aspects

The approximate rating ranges from 0 to 5. A `` 0 '' mark means a total incapacity to handle the case, a `` 1 '' means that the point is either disregarded or very roughly taken into account. `` 2 '' and `` 3 '' marks mean quite a good job, a `` 4 '' is very good if not best in the category, while a `` 5 '' means more or less perfect. Of course, the rating is quite subjective, it is very difficult to estimate precisely all aspects of all methods. Most probably, many authors would disagree with the rating of their method, and sometimes with the rating of the others too. It is not our claim to give an absolute nor definitive value to this approximate classification.

multiple arrivals amplitude robustness

Nichols 3 5 4

paraxial ray tracing 4 3 2

Gaussian beams 5 4 3

NORSAR 4 3 3

Lomax 4 3 2

Vidale 1 1

van Trier-Symes 1 1

Lin Zhang 1 2 3

Podvin 3 1 5

Moser 5

$\bullet$

Multiple arrivals: Nichols' and the ray tracing family methods are able to track the first or the most energetic arrivals, at user's choice. Lin Zhang's might do so, Podvin can do as well, but both at very high computational expenses.

$\bullet$

Handling of amplitude:

Nichols' band-limited gives a very good estimate of the amplitude, while Gaussian beams and other ray tracing family methods would be perfect but only for infinite frequency. Lin Zhang's method is able to compute amplitude information, but it will be very expensive to track all the ``possible'' most energetic arrivals. All the ``first arrival, no amplitude'' methods are given a 1 mark. They can yield a minimalist estimate of amplitude as a simple a posteriori function of traveltime.

$\bullet$

Robustness:

This criterion covers robustness and stability of the method, and the reliability of the results. Podvin and Moser get maximal mark for their unconditional stability and because they always find first arrival traveltime, for a given velocity model. In its multiple arrival implementation, Podvin's is probably robust as well. Nichols' is stable, and is reasonably reliable for both first and most energetic arrival traveltimes: hence a 4 mark. Ray methods get medium marks, for the good stability of ray tracing, balanced by the uncertain reliability of the results in complex media. Lin Zhang is stable, robust, even precise but may not always find first or more energetic arrivals. Vidale and van Trier-Symes get a 1 mark, because they are neither stable nor able to always find first arrivals.

Comparative rating: En route to 3-D

Cost in 2-D Reasonable cost in 3-D

Nichols 3 3

paraxial ray tracing 4 4

gaussian beams 3 3

NORSAR 4 4

Lomax 2 2

Vidale 5 4

van Trier-Symes 4 4

Lin Zhang 3 2

Podvin 3 2

Moser 5 5

$\bullet$

Cost in 2-D: The cost in 2-D is rather easy to assess, and might give a first idea of the eventual cost in 3-D. At the top of the list, we have the first arrival methods: finite-differencing of the eikonal (Vidale, van Trier) and graph theory. Ray based methods are reasonably cheap. At the bottom of the list we have Podvin and Lin Zhang methods, not very convenient for computer optimization. We also find Nichols' method expensive in 3D, because it is a rich and thus expensive method.

$\bullet$

Feasibility of a 3-D version:

This criterion was not rated. All methods reviewed are theoretically speaking portable to 3-D. A few of them already exist in 3-D: ray tracing, Podvin's, Norsar, Vidale's. A 3-D version of Moser's (graph theory) is unproblematic. A 3-D version of Zhang's code requires a lot of programming effort, and it is the same for Nichols' method. Both authors are working on a 3-D implementation. 3-D Gaussian beams are not yet current practice, and finally van Trier-Symes has some stability problems in 3-D.

$\bullet$

Reasonable cost of the 3-D version:

The cost in 3-D, can be extrapolated from the cost in 2-D (previous table) with some scaling factor not easy to estimate and depending on the method: for instance the unit cost of 3-D Podvin is around 9 times the unit cost of the 2-D version. On the field of 3-D cost, the ray tracing family keeps a serious advantage, for the paucity of data it handles (rays instead of wavefield), and for the intrinsic coarse grain parallelism of rays. For reason of vectorization and parallelization, Vidale and van Trier-Symes (finite differencing the eikonal), should be a top performer, once stability problems are solved. Graph methods (Moser) should be very efficient too. For its resistance to both vectorization and parallelization, Podvin's receives a mere 2. Nichols' and Lin Zhang's pay the due price of richness.

CONCLUSION

At the end of this survey, we hope that the reader will have a clearer view of the principles and the characteristics of the analyzed method. These traveltimes methods were designed for a specific goal: delivering approximate Green's functions for Kirchhoff migration and tomographic inversion or velocity analysis. From the point of view of those processes, the critical aspects are the cost efficiency, the reliability, the kinematical accuracy and the dynamical accuracy, of course all of it in a 3-D world. There is no winner on all these aspects, and hardly any on one aspect at a time. The reader is thus invited to select from our review the method more appropriate to his specific purpose. Some aspects more specifically related to imaging and migration are further analyzed in the next paper, ().

[GEOTLE,MISC,SEP,geophysics,rest]