Introduction

The largest portion of the computational effort in migrating data in heterogeneous media with an integral method is the calculation of the required Green's functions. With asymptotic ray theory, the approximate Green's functions for compressional waves and shear waves can be decoupled (Cervený et al., 1977). Each Greens' function is characterized by a traveltime function and an amplitude function. Traditionally, these traveltime and amplitude functions are obtained by performing dynamic ray tracing to compute the traveltimes and amplitudes along rays, and then interpolating them onto a regular grid (Cervený,1987; Beydoun and Keho, 1987). However, for complicated slowness models, the interpolation process may be theoretically invalid and computationally expensive.

Vidale (1988) developed a finite difference method that computes the traveltimes directly on a regular grid in Cartesian coordinates. The method can be interpreted as a wavefront extrapolation process that starts from the source location and proceeds on an expanding square box. Van Trier and Symes (1991) described a slightly different algorithm in which the gradients of the traveltimes, instead of the traveltimes themselves, are extrapolated. Their method conducts the calculation in polar coordinates and transforms the results to the Cartesian coordinates afterwards. Van Trier and Symes showed that their method often achieves better efficiency and accuracy than Vidale's method. However, the results of many numerical experiments show that both algorithms encounter serious stability problems when applied to models containing large slowness discontinuities.

Podvin and Lecomte (1991) reported several significant improvements on Vidale's method. First, their algorithm systematically applies Huygens' principle in the finite difference approximation so that the method provides accurate traveltimes in the presence of extremely severe, arbitrarily shaped slowness contrasts. Second, they developed a massively parallel implementation of the algorithm. In this implementation, each grid point is associated with a processor. All the processors simultaneously update the traveltimes at the associated grid points, which is like a relaxation process. Third, to ensure the results of the sequential implementation are minimum traveltimes, Podvin and Lecomte introduced backward extrapolation. This procedure captures the first arrivals that propagate into the square computational box.

While acknowledging the promising results given by Podvin and Lecomte, I find three weaknesses in their algorithm, for which better solutions may be explored. First the massive parallel implementation does not effectively utilize computational power because, during the computation, most of the processors away from the first arrival wavefront perform dummy operations. In 3-D, the computation involved in each processor is quite complicated. Second, similar to Vidale's algorithm, the sequential implementation can not be vectorized. Third, the backward extrapolation introduces extra cost and should be avoided whenever possible.

I have developed a new traveltime calculation method that combines the advantages of the existing algorithms. I use a simplified finite difference representation of Podvin and Lecomte's scheme so that the new algorithm works for large contrast slowness models. I choose to work in polar coordinates to make the algorithm fully vectorizable. This choice also avoids performing backward extrapolation if one assumes that waves seldom propagate back toward the source. The stability of the algorithm is ensured by adapting the step size of the extrapolation to the directions of the wave propagation. The algorithm can also be used to solve the partial differential equations associated with the amplitude calculation. The resulted algorithm is relatively more efficient than Vidale's method (1990). The cost of adding the amplitude calculation to the algorithm is about twice the cost of calculating the traveltimes only.

This paper is organized into three parts. The first part explains the required theory. I begin with a brief review of some well-known results in the asymptotic ray theory. Then, I address the theoretical aspects of the traveltime and amplitude calculations. In the second part, I describe the implementational details of the finite difference algorithm. The third part of the paper demonstrates the method with several examples. Finally, I conclude the paper by discussing the applications of the method.