3-D geometric ray tracing

Assume that a ray is specified by the parametric equation

$$x_{\rm i} = x_{\rm i}(\tau), \hspace*{0.3in}i=1,2,3$$ (2)

where $\tau$ is the traveltime along the ray, with $\tau=0$ at the location of the source. The independent variable is not necessarily $\tau$. It can be any other variable which increases monotonously along the ray, e.g., arclength is a common choice in the ray tracing algorithm. Here we choose traveltime because seismic data are sampled along time axis. Therefore, it is more convenient to compare the ray tracing results with seismic data.

The three coordinates are the solution of geometric ray tracing system  

$$\frac{dx_{\rm i}}{d\tau}=s^2 p_{\rm i}, \hspace*{0.3in} \f... ...{\rm i}}{d\tau}=\frac{1}{s}{\partial s\over \partial x_{\rm i}}$$ (3)

where s is the slowness along the ray. The reason we use slowness, not velocity, is to follow the convention of tomography. Here, $p_{\rm i}$ are the Cartesian components of the slowness vector $\vec{p}$ and

$$p_{\rm i}={\partial t \over \partial x_{\rm i}}, \hspace*{0.3in} i=1,2,3$$ (4)

Equation (3) is first-order ordinary differential equations. Usually, we cannot find an analytical solution for a general velocity model. It is more common to use some numerical algorithms, e.g., the Runge-Kutta method. In this paper, we use a Fortran90 version Runge-Kutta solver. In order to fit our expression to the convention of rksuite90 package, the geometric ray tracing system is expressed as

$$\frac{d}{d \tau} \left[\begin{array} {cc} x_1 & p_1 \\ x_2... ...\frac{1}{s} {\partial s \over \partial x_3} \end{array} \right]$$ (5)

The initial condition is

$$\left[\begin{array} {cc} x_1^0 & p_1^0 \\ x_2^0 & p_2^0 \\... ...n \beta_0 \\ x_{30} & s_0 \, \cos \alpha_0 \end{array} \right]$$ (6)

where $\alpha_{0}$ and $\beta_{0}$ are the two initial take-off angles, as shown in Figure 1.

 

take-off Figure 1 The initial condition of each central ray. (x, y, z) forms a Cartesian coordinate system. O is the source location.