3-D dynamic ray tracing

Dynamic ray tracing will provide us dynamic information like amplitude and phase-shift, not only kinematic information like raypaths and traveltimes. More specifically, we can estimate the traveltime, amplitude, and phase-shift in the vicinity of the central ray. In order to do this, we need to introduce another two coordinate systems, i.e., ray-centered coordinate system (q₁, q₂, q₃) and ray parameter coordinate system ($\gamma_1$, $\gamma_2$, $\gamma_3$).

The ray-centered coordinate system is a curvilinear orthogonal coordinate system. Let us consider an arbitrarily selected ray $\Omega$, which is specified by ray parameter $\gamma_1$, $\gamma_2$. The ray-centered coordinates q₁, q₂, q₃ connected with the ray $\Omega$ are defined in the following way: q₃ corresponds to the traveltime along the ray $\Omega$. q₁ and q₂ form a 2D orthogonal Cartesian coordinate system in the plane perpendicular to $\Omega$, as shown in Figure 2.

 

ray-center Figure 2 The ray-centered coordinate system. q3 is tangent to the raypath. q1 and q2 forms a plane perpendicular to q3.

The physical meaning of the ray parameter coordinate system is very obvious. $\gamma_1$ and $\gamma_2$ are two ray parameters in the x₁ and x₂ direction, which are constant along the whole ray and variable from one ray to another. Their definitions are as follows

$$\gamma_1 = {\partial t \over \partial x_1}, \hspace*{0.3in} \gamma_2 = {\partial t \over \partial x_2}$$ (7)

$\gamma_3$ is chosen as q₃. It increases monotonously along the ray.

Dynamic ray tracing system can be expressed as  

$$\frac{d}{dt}{\bf Q}=s^{-2}{\bf P}, \hspace*{0.3in} \frac{d}{dt}{\bf P}=-s{\bf V}{\bf Q}$$ (8)

where Q, P, V are $2\times 2$ matrices.

$${\bf Q}= \left[\begin{array} {cc} Q_{11} & Q_{12} \\ Q_{21... ...} {cc} V_{11} & V_{12} \\ V_{21} & V_{22} \end{array} \right]$$ (9)

The definition of each component is

$$Q_{\rm IJ}={\partial q_{\rm I} \over \partial \gamma_{\rm J}... ...IJ}={\partial v \over {\partial q_{\rm I} \partial q_{\rm J}}}$$ (10)

where $q_{\rm I} (I=1,2)$ is the ray-centered coordinate. $\gamma_{\rm I} (I=1,2)$ is the ray parameter associated with each ray. $v=\frac{1}{s}$ is the medium velocity.

Mathematically, Q and P can be interpreted as transformation matrices: Q is a transformation matrix from the ray parameters $\gamma_1$, $\gamma_2$ to the ray-centered coordinate q₁, q₂, P is a transformation matrix from the ray parameters $\gamma_1$, $\gamma_2$ to the slowness vector component in the ray-centered coordinate system. Physically, matrix Q measures the derivations of paraxial rays from ray $\Omega$ and is also referred to as the matrix of geometrical spreading. Matrix P has no obvious physical meaning. But we can define a $2\times 2$ matrix M of second derivatives of the traveltime field with respect to the ray-centered coordinate q₁, q₂ by the relationship

$$M_{\rm IJ}={\partial^2 \tau \over {\partial q_{\rm I} \partial q_{\rm J}}}.$$ (11)

Then we can write

$${\bf M}(\tau)={\bf P}{\bf Q}^{-1}$$ (12)

In order to fit the requirement of Runge-Kutta package, equation (8) is expressed as

$$\frac{d}{d \tau} \left[\begin{array} {cc} Q_{11} & P_{11} \... ...s {\partial^2 v \over \partial q_1^2}Q_{22} \end{array} \right]$$ (13)

The initial condition is chosen as

$$\left[\begin{array} {cc} Q_{11}^0 & P_{11}^0 \\ Q_{12}^0 &... ...} {cc} 1 & 1 \\ 0 & 0 \\ 0 & 0 \\ 1 & 1 \end{array} \right]$$ (14)

Here we set matrix ${\bf Q^0}$ and ${\bf P^0}$ to be $2\times 2$ identity matrix. The reason for choosing an identity matrix is to make sure that ${\bf M}$ is a symmetric matrix.

The paraxial ray tracing solution of the pressure field can be written as

 

$$\psi ({\bf x})= C \biggl ( { \alpha \over \alpha_0 ~det{\bf ... ... e^{i \omega (\tau({\bf 0}) + {1 \over 2} {\bf q}^T {\bf M q})}$$ (15)

Sub-index denotes a value taken at the source. C is a factor dependent upon the source, and is constant along the ray.