Though traveltime computation is widely used in seismic modeling and routine data processing, attaining sufficient accuracy without compromising speed and robustness is problematic. Moreover, there is no easy way to obtain the traveltimes corresponding to the multiple arrivals that appear in complex velocity media.
The tradeoff between speed and accuracy becomes apparent in the choice between the two most commonly used methods, ray tracing and numerical solutions to the eikonal equation. Other methods reported in the literature (dynamic programming Moser (1991), wavefront construction Vinje et al. (1993), etc.) are less common in practice Audebert et al. (1994).
Eikonal solvers provide a relatively fast and robust method of traveltime computations Vidale (1990); van Trier and Symes (1991). They also avoid the problem of traveltime interpolation to a regular grid which imaging applications require. However, the eikonal solvers compute first-arrival traveltimes and lack the important ability to track multiple arrivals. In complex velocity structures, the first arrival does not necessarily correspond to the most energetic wave, and other arrivals can be crucially important for accurate modeling and imaging Geoltrain and Brac (1993); Gray and May (1994).
On the other hand, one-point ray tracing can compute multiple arrivals with great accuracy. Unfortunately, it lacks the robustness of eikonal solvers. Increasing the accuracy of ray tracing in the regions of complex velocity variations raises the cost of the method and makes it prohibitively expensive for routine large-scale applications. Mathematically, ray tracing amounts to a numerical solution of the initial value problem for a system of ordinary differential equations Cervený (1987). These ray equations describe characteristic lines of the eikonal partial differential equation.
Here, we propose a somewhat different approach to traveltime computation, that is both fast and accurate, and has the ability to find multiple arrival traveltimes. The theoretical construction is based on a system of differential equations, equivalent to the eikonal equation, but formulated in the ray coordinate system. Unlike eikonal solvers, our method produces the output in ray coordinates. Unlike ray tracing, it is computed by a numerical solution of partial differential equations. We show that the first-order discretization scheme has a remarkably simple interpretation in terms of the Huygens' principle and propose a Huygens wavefront tracing (from now on referred to as HWT ) scheme as a robust alternative to conventional ray tracing. Numerical examples demonstrate the following properties of the method: stability in media with strong and sharp lateral velocity variations, better coverage of the shadow zones, and greater speed than paraxial ray tracing (from now on referred to as PRT ).