A second-order fast marching eikonal solver

Fast marching and the eikonal equation

Under a high frequency approximation, propagating wavefronts may be described by the eikonal equation,

$$\left( \frac{\partial t}{\partial x} \right)^2 + \left( \fra... ... + \left( \frac{\partial t}{\partial z} \right)^2 =s^2(x,y,z),$$ (1)

where $t$ is the traveltime, $s$ is the slowness, and $x$, $y$ and $z$ represent the spatial Cartesian coordinates.

The fast marching method solves equation (1) by directly mimicking the advancing wavefront. Every point on the computational grid is classified into three groups: points behind the wavefront, whose traveltimes are known and fixed; points on the wavefront, whose traveltimes have been calculated, but are not yet fixed; and points ahead of the wavefront. The algorithm then proceeds as follows:

1. Choose the point on the wavefront with the smallest traveltime.
2. Fix this traveltime.
3. Advance the wavefront, so that this point is behind it, and adjacent points are either on the wavefront or behind it.
4. Update traveltimes for adjacent points on the wavefront by solving equation (1) numerically.
5. Repeat until every point is behind the wavefront.

The update procedure (step 4.) requires the solution of the following quadratic equation for $t$,

$$\max(D_{ijk}^{-x} t, 0)^2+ \min(D_{ijk}^{+x} t, 0)^2 +$$

$$\max(D_{ijk}^{-y} t, 0)^2+ \min(D_{ijk}^{+y} t, 0)^2 +$$

$$\max(D_{ijk}^{-z} t, 0)^2+ \min(D_{ijk}^{+z} t, 0)^2 = s_{ijk}$$ (2)

where $D_{ijk}^{-x}$ is a backward $x$ difference operator at grid point, $ijk$, $D_{ijk}^{+x}$ is a forward $x$ operator, and finite-difference operators in $y$ and $z$ are defined similarly. The roots of the quadratic equation, $a t^2 + b t + c = 0$, can be calculated explicitly as

$$t = \frac{ -b \pm \sqrt{b^2 -4ac}}{2a}.$$ (3)

Solving equation (2) amounts to accumulating coefficients $a$, $b$ and $c$ from its non-zero terms, and evaluating $t$ with equation (3).

If we choose a two-point finite-difference operator, such as

$$D_{ijk}^{-x} t = \; \frac{t_{ijk}-t_{(i-1)jk}}{\Delta x}$$ (4)

$${\rm then} \; (D_{ijk}^{-x} t)^2 = \; \alpha t_{ijk}^2 + \beta t_{ijk} + \gamma$$ (5)

where $\alpha = \frac{1}{\Delta x^2}$, $\beta = -2 t_{(i-1)jk} \alpha$ and $\gamma = t_{(i-1)jk}^2 \alpha$. Coefficients $a$, $b$ and $c$ can now be calculated from $a = \Sigma_l \alpha_l$, $b = \Sigma_l \beta_l$, and $c = \Sigma_l \gamma_l - s^2$, where the summation index, $l$, refers to the six terms in equation (2) subject to the various min/max conditions.

This two-point stencil, however, is only accurate to first-order. If instead we choose a suitable three-point finite-difference stencil, we may expect the method to have second-order accuracy. For example, the second-order upwind stencil,

$$D_{ijk}^{-x} t = \; \frac{3t_{ijk}-4t_{(i-1)jk} + t_{(i-2)jk}}{2 \Delta x}$$ (6)

$${\rm gives} \; (D_{ijk}^{-x} t)^2 = \; \alpha' t_{ijk}^2 + \beta' t_{ijk} + \gamma'$$ (7)

$$\begin{eqnarray*} {\rm where \; this \; time \;} & \alpha' & = \; \frac{9}{4 \De... ...1)jk}-t_{(i-2)jk})^2}{4 \Delta x^2} = \alpha' {t'}_{(i-1)jk}^2, \end{eqnarray*}$$

$$\begin{eqnarray*} {\rm and} \; \; & {t'}_{(i-1)jk} & = \; \frac{1}{3}(4 t_{(i-1)jk}-t_{(i-2)jk}). \end{eqnarray*}$$

Coefficients, $a$, $b$ and $c$ can be accumulated from $\alpha'$, $\beta'$ and $\gamma'$ as before, and if the traveltime, $t_{(i-2)jk}$ is not available, first-order values may be substituted.

A second-order fast marching eikonal solver