ABSTRACTWe 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.
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
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
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.
PARAXIAL RAY TRACING
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
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.
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.
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
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 .
The method has the same features here as the paraxial ray tracing method.
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.
The Gaussian beam method is based upon the high frequency, asymptotic ray approximation, and diffraction effects are not handled correctly.
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.
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.
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.
The only weaknesses we can see are the infinite frequency assumption (or rather, infinite bandwidth assumption) and the necessity of a C1 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.
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.
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.
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.
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 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 ().
VAN TRIER AND SYMES METHOD
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.
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.
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.
The velocity model has a continuous representation in triangular cells of constant velocity gradient.
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.
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.
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 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.
The velocity model is represented by a rectangular grid of constant slowness cells. The grid is either 2-D or 3-D.
SCHNEIDER ET AL'S METHOD
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.
The methods can handle any velocity model, but the slowness is constant within rectangular cells.
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.
The method can only compute absolute first arrival traveltimes.
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.
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.
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.
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.
|paraxial ray tracing||4||3||2|
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.
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
|3|c| 2-D to 3-D criteria|
|Cost in 2-D||Reasonable cost in 3-D|
|paraxial ray tracing||4||4|
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.
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.
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, ().