Review of HWT theory

Given an isotropic heterogeneous medium, wavefronts are represented by surfaces of equal traveltime, constrained by the eikonal equation  

$$\left(\frac{\partial \tau}{\partial x}\right)^2 + \left(\f... ...c{\partial \tau}{\partial z}\right)^2 = \frac{1}{v^2(x,y,z)}\;$$ (1)

and appropriate boundary conditions.

Each point on a wavefront can be parametrized by either its Cartesian coordinates x, y, and z, or its ray coordinates, which consist of the traveltime $\tau$, and the two shooting angles at the source, $\gamma$ and $\phi$.

For complex velocity fields, the ray coordinates as a function of the Cartesian coordinates become multi-valued, in other words, there is more than one ray going through a given point in the subsurface. In contrast, the Cartesian coordinates as a function of the ray coordinates remain single-valued, that is, there is one unique position in the subsurface where a ray, shot with two particular shooting angles, arrives at a given time. Figure 1 illustrates the difference between the two representations of the wavefronts.

 

coord
Figure 1 The ray coordinates as a function of the Cartesian coordinates are multi-valued (left). The Cartesian coordinates as a function of the ray coordinates are single-valued (right).

Since $x (\tau,\gamma,\phi)$, $y (\tau,\gamma,\phi)$, and are uniquely defined for arbitrarily complex velocity fields, the eikonal equation (Equation 1) can be transformed to another form that is better suited for analysis in ray coordinates Sava and Fomel (1998):  

$$\left(\frac{\partial x}{\partial \tau}\right)^2 + \left(\f... ...\partial z}{\partial \tau}\right)^2 = v^2 \left(x,y,z\right).$$ (2)

Converting equation (2) to a finite-difference equation using a first-order discretization scheme, we obtain  

$$\left(x_{j+1}^{i,k}-x_j^{i,k}\right)^2 + \left(y_{j+1}^{i,... ...z_{j+1}^{i,k}-z_j^{i,k}\right)^2 = \left(r_j^{i,k}\right)^2\;,$$ (3)

where j is the index of the current wavefront, j+1 is the index of the new wavefront to be computed, and i and k are the indices of the shooting angles. This equation represents a sphere, the wavefront of a secondary Huygens source placed at (x_j^i,k,y_j^i,k,z_j^i,k) on the current wavefront.

According to the Huygens principle, the new wavefront is the envelope of all the secondary wavefronts. Mathematically, the position of the new wavefront is described by a system of three equations composed of Equation (3) and the following two equations Sava and Fomel (1998):

$$\begin{eqnarray} \left(x_j^{i,k} - x_{j+1}^{i,k} \right)\,\left(x_j^{i+1,k} - ... ...i-1,k}\right) = r_j^{i,k}\left(r_j^{i+1,k} - r_j^{i-1,k}\right)\;\end{eqnarray}$$ (4)

and

$$\begin{eqnarray} \left(x_j^{i,k} - x_{j+1}^{i,k} \right)\,\left(x_j^{i,k+1} - ... ...,k-1}\right) = r_j^{i,k}\left(r_j^{i,k+1} - r_j^{i,k-1}\right)\;.\end{eqnarray}$$ (5)

Figure 2 contains a simple geometrical interpretation of the system described by Equations (3), (4), and (5).

 

huygens3d Figure 2 A geometrical updating scheme for 3-D HWT in the physical domain. Five points on the current wavefront, represented by the five spheres, not all visible, with radii defined by the velocities at the corresponding points of the wavefront, are used to compute a point on the next wavefront. The sphere in the middle represents equation (3), and the planes represent equations (4) and (5).

Five points on the current wavefront, represented by the five spheres, not all visible, with radii defined by the velocities at the corresponding points of the wavefront, are used to compute a point on the next wavefront. The sphere in the middle represents equation (3), while the planes represent equations (4) and (5).

Huygens wavefront tracing, based on the system of equations (3), (4), and (5), is nothing but an explicit finite-difference method in the ray coordinate system. The coordinates of the new wavefronts are computed according to those of the current wavefront.

 

scheme
Figure 3 Finite-difference traveltime computation scheme. A 3-point stencil is needed in 2-D to compute the centered finite-difference representation of the derivative with respect to the shooting angle (left). A 5-point stencil is needed in 3-D to compute the centered finite-difference representation of the derivatives with respect to the shooting angles (right)

A three-point stencil is needed in two dimensions to compute the centered finite-difference representation of the derivative with respect to the shooting angle, while a five-point stencil is needed in three dimensions to compute the centered finite-difference representation of the derivatives with respect to the shooting angles. Figure 3 is a graphical illustration of the finite-difference stencils.