A discretization scheme and the Huygens' principle

A natural first-order discretization scheme for equation (5) leads to the difference equation  

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

where the index i corresponds to the ray parameter $\gamma$, j corresponds to the traveltime $\tau$, $r_j^i = \triangle \tau\,v_j^i$,$\triangle \tau$ is the increment in time, and v_jⁱ is the velocity at the $\{i,j\}$ grid point. It is easy to notice that equation (7) simply describes a sphere (or a circle in two dimensions) with the center at $\{ x_j^i, z_j^i \}$ and the radius r_jⁱ. This sphere is, of course, the wavefront of a secondary Huygens source.

This observation suggests that we apply the Huygens' principle directly to find an appropriate discretization for equation (6). Let us consider a family of Huygens spheres, centered at the points along the current wavefront. Mathematically, this family is described by an equation analogous to (7), as follows:  

$$\left(x-x (\gamma)\right)^2 + \left(z-z (\gamma)\right)^2 = r^2 (\gamma)\;.$$ (8)

Here the ray parameter $\gamma$ serves as the parameter that distinguishes a particular Huygens source. According to the Huygens' principle, the next wavefront corresponds to the envelope of the wavefront family. To find the envelop condition, we can simply differentiate both sides of equation (8) with respect to the family parameter $\gamma$. The result takes the form  

$$\left(x (\gamma) - x \right)\, x'(\gamma) + \left(z (\gamma) - z \right)\, z'(\gamma) = r (\gamma) r'(\gamma)\;,$$ (9)

which is clearly a semidiscrete analog of equation (6). To complete the discretization, we can represent the $\gamma$-derivatives in (9) by a centered finite-difference approximation. This representation yields the scheme  

$$\left(x_j^i - x_{j+1}^i \right)\, \left(x_j^{i+1} - x_j^{i-... ..._j^{i-1}\right) = r_j^i \left(r_j^{i+1} - r_j^{i-1}\right)\;,$$ (10)

which supplements the previously found scheme (7) for a unique determination of the point $\{x_{j+1}^i,z_{j+1}^i\}$ on the i-th ray and the (j+1)-th wavefront. Formulas (7) and (10) define an update scheme, depicted in Figure 1. To fill the $\{\tau,\gamma\}$ plane, the scheme needs to be initialized with one complete wavefront (around the wave source) and two boundary rays.

The solution of system (7-10) has the explicit form

$$\begin{eqnarray} x_{j+1}^i & = & x_j^i - r_j^i \left( \alpha\,\left(x_j^{i+1} ... ...}\right) \mp \beta\,\left(x_j^{i+1} - x_j^{i-1}\right)\right)\;, \end{eqnarray}$$ (11) (12)

where

$$\begin{eqnarray} \alpha & = & \frac{r_j^{i+1} - r_j^{i-1}} {\left(x_j^{i+1} - x_j^{i-1}\right)^2 + \left(z_j^{i+1} - z_j^{i-1}\right)^2} \;, \ \end{eqnarray}$$ (13)

and

$$\begin{eqnarray} \beta & = & \mbox{sign}\left( x_j^{i+1} - x_j^{i-1} \right) ... ...j^{i-1}\right)^2 + \left(z_j^{i+1} - z_j^{i-1}\right)^2 } \;.\end{eqnarray}$$ (14)

 

scheme Figure 1 An updating scheme for HWT. Three points on the current wavefront (A, B, and C) are used to advance in the $\tau$ direction.

Figure 2 shows a geometric interpretation of formulas (7) and (10). Formula (10) is clearly a line equation. Thus, the new point D in Figure 2 is defined as one of the two intersections of this line with the B sphere, defined by formula (7). It is easy to show geometrically that the newly created ray segment BD is orthogonal to the common tangent of spheres A and C. Within the finite-difference approximation, the common tangent reflects local wavefront behavior.

 

huygens Figure 2 A geometrical updating scheme for HWT in the physical domain. Three points on the current wavefront (A, B, and C) are used to compute the position of the D point. The bold lines represent equations (7) and (10). The tangent to circle B at point D is parallel to the common tangent of circles A and C.