3-D imaging and velocity analysis methods are computationally intensive. One reason is that they are iterative processes, which need a forward modeling algorithm as a kernel component. For example, 3-D Kirchhoff-type imaging usually needs a 3-D Green's function as input. There are many different ways of estimating Green's function. They can be roughly categorized into two families according to both accuracy and efficiency: limited bandwidth methods and infinite frequency methods Audebert et al. (1994).
Limited bandwidth methods usually deal with frequency-dependent, complex-valued Green's function. These methods are typically the extrapolation or modeling algorithms working in the time domain or the temporal frequency domain, e.g., phase-shift or finite-difference downward continuation. Nichols 1994 presented a method of calculating band-limited traveltime explicitly, which is essentially a compromise between full wavefield modeling and ray tracing. In terms of accuracy, it is a reduced version of a full wavefield modeling; in terms of efficiency, the cost of this new method is lower than that of full wavefield modeling, but higher than ray tracing method. One main weakness with this method is that no explicit raypaths are calculated and stored. In many applications, e.g., tomographic velocity analysis, raypaths are indeed necessary.
Ray tracing method is the main sub-family of the infinite frequency methods. The simplest and also the fastest ray tracing algorithm is geometric ray tracing. It transfers the eikonal equation to the ray equation using the characteristic method and then estimates the raypaths by solving the ray equation. Geometric ray tracing can only provide traveltime and cannot provide any amplitude and phase-shift information. It has been proved that first arrivals cannot be used to image complex structures correctly Rekdal and Biondi (1994). Dynamic ray tracing Cervený (1985), is an extension to geometric ray tracing. Dynamic ray tracing (including paraxial ray approximation and Gaussian Beam approximation) can provide not only multiple traveltimes but also amplitude and phase-shift information. Therefore it can also be regarded as a compromise between full wavefield modeling and ray tracing modeling.
In this paper, on the basis of Cervený's work, we develop an accurate and robust 3-D dynamic ray tracing system. First, our implementation is independent of the model gridding scheme, which can be layer, cube, or tetrahedron, etc. Therefore, it is easy to incorporate our implementation into different 3-D velocity model representations. Second, since dynamic ray tracing is a high-frequency approximation, it has some requirements over the smoothness of the model. In order to make the model more ``ray valid'', we apply a triangle smoothing filter to the model along the x, y, z axis. The smoothed model makes not only our numerical solver more stable, but also the results from dynamic ray tracing more like band-limited seismic data. We use a Fortran90 version Runge-Kutta solver in our application Brankin and Gladwell (1995). This solver has the ability to adapt its step-length in accordance with the local gradient of the slowness, which makes it more accurate than constant step-length solver and also enables us to deal with more complex structures without smoothing it dramatically. Another advantage is that Fortran90 is easily parallelized and we expect our implementation would be more efficient when implemented on a parallel computer.
In the next several sections, we first discuss our representation of the model and then present a brief version of the dynamic ray tracing theory. After that, we show three results from applying our dynamic ray tracing program to synthetic models. The linear model verifies that our implementation is correct when there is a theoretical solution. The layer model shows the importance of smoothing the model. The saltdome model demonstrates our method's robustness.