Variational principles on a grid

In this section, I derive a discrete traveltime computation procedure, based solely on Fermat's principle, and show that on a Cartesian rectangular grid it is precisely equivalent to the update formula (1) of the first-order eikonal solver.

 

triangle Figure 3 A geometrical scheme for the traveltime updating procedure in two dimensions.

For simplicity, let us focus on the two-dimensional case. Consider a line segment with the end points A and B, as shown in Figure 3. Let t_A and t_B denote the traveltimes from a fixed distant source to points A and B, respectively. Define a parameter $\xi$ such that $\xi=0$ at A, $\xi=1$ at B, and $\xi$ changes continuously on the line segment between A and B. Then for each point of the segment, we can approximate the traveltime by the linear interpolation formula  

$$t (\xi) = (1-\xi) t_A + \xi t_B\;.$$ (9)

Now let us consider an arbitrary point C in the vicinity of AB. If we know that the ray from the source to C passes through some point $\xi$ of the segment AB, then the total traveltime at C is approximately  

$$t_C = t (\xi) + s_C\,\sqrt{\vert AB\vert^2 (\xi-\xi_0)^2 + \rho_0^2}\;,$$ (10)

where s_C is the local slowness, $\xi_0$ corresponds to the projection of C to the line AB (normalized by the length |AB|), and $\rho_0$ is the length of the normal from C to $\xi_0$.

Fermat's principle states that the actual ray to C corresponds to a local minimum of the traveltime with respect to raypath perturbations. According to our parameterization, it is sufficient to find a local extreme of t_C with respect to the parameter $\xi$. Equating the $\xi$ derivative to zero, we arrive at the equation  

$$t_B - t_A + \frac{s_C\,\vert AB\vert^2\,(\xi-\xi_0)} {\sqrt{\vert AB\vert^2 (\xi-\xi_0)^2 + \rho_0^2}} = 0\;,$$ (11)

which has (as a quadratic equation) the explicit solution for $\xi$: 

$$\xi = \xi_0 \pm \frac{\rho_0\,(t_A - t_B)} {\vert AB\vert\,\sqrt{s_C^2\,\vert AB\vert^2 - (t_A - t_B)^2}}\;.$$ (12)

Finally, substituting the value of $\xi$ from (12) into equation (10) and selecting the appropriate branch of the square root, we obtain the formula

$$\begin{eqnarray} c\,t_C & = & \rho_0\,\sqrt{s_C^2 c^2 - (t_A - t_B)^2} + c\,t_... ... - (t_A - t_B)^2} + a\,t_A\,\cos{\beta} + b\,t_B\,\cos{\alpha}\;,\end{eqnarray}$$ (13)

where c = |AB|, a = |BC|, b = |AC|, angle $\alpha$corresponds to $\widehat{BAC}$, and angle $\beta$ corresponds to $\widehat{ABC}$ in the triangle ABC (Figure 3).

 

square Figure 4 A geometrical scheme for traveltime updating on a rectangular grid.

To see the connection of formula (13) with the eikonal difference equation (1), we need to consider the case of a rectangular computation cell with the edge AB being a diagonal segment, as illustrated in Figure 4. In this case, $\cos{\alpha} = \frac{a}{c}$, $\cos{\beta} = \frac{b}{c}$, $\rho_0 = \frac{ab}{c}$, and formula (13) reduces to  

$$t_C = \frac{ab\,\sqrt{s_C^2 (a^2 + b^2) - (t_A - t_B)^2} + a^2\,t_A + b^2\,t_B} {a^2+b^2}\;.$$ (14)

We can notice that (14) is precisely equivalent to the solution of the quadratic equation (13), which in our new notation takes the form  

$$\left(\frac{t_C - t_A}{b}\right)^2 + \left(\frac{t_C - t_B}{a}\right)^2 = s_C^2\;.$$ (15)

What have we accomplished by this analysis? First, we have derived a local traveltime computation formula for an arbitrary grid. The derivation is based solely on Fermat's principle and a local linear interpolation, which provides the first-order accuracy. Combined with the fast marching evaluation order, which is also based on Fermat's principle, this procedure defines a complete algorithm of first-arrival traveltime calculation. On a rectangular grid, this algorithm is exactly equivalent to the fast marching method of Sethian (1996a) and Sethian and Popovici (1997). Second, the derivation provides a general principle, which can be applied to derive analogous algorithms for other eikonal-type (Hamilton-Jacobi) equations and their corresponding variational principles.