fast marching in spherical coordinates

Details of solving the eikonal equation using the fast marching method in Cartesian coordinates can be found in Cao and Greenhalgh (1994) and Sethian (1996), while its application in spherical coordinates can be found in Alkhalifah and Fomel (1997). A summary of the approach, however, is revisited below.

The first-order nature of the fast marching method results in large errors for conventional sparse grid-point configurations. The largest of these errors occur when the curvature of the wavefront is large and the wavefront is traveling diagonally with respect to the orthogonal coordinate system. However, if the wavefront curvature is zero, or the wavefront is parallel to the coordinate system, the fast marching method becomes exact. For waves emanating from a point source, large errors appear in the Cartesian coordinate implementation of the method at angles of 45 degrees from the vertical, especially near the source, where the wavefront curvature is large. Plane waves, on the other hand, are calculated exactly in the Cartesian-coordinates, because in such coordinates plane waves have no curvature.

In spherical coordinates (Figure 1), waves emanating from a point source are effectively propagated as plane waves on a regular grid. For homogeneous media, these plane waves have fronts that are always parallel to the $\theta$-$\phi$ plane. As a result, traveltime calculation using the fast marching approach in spherical coordinates is always exact in homogeneous media Alkhalifah and Fomel (1997).

 

polar-dif Figure 1 A spherical coordinate system given by r, $\theta$, and $\phi$. The source, s(x0,y0,z0), is at the origin of the spherical coordinates where r=0. The parameter r (=$\sqrt{(x-x_0)^2+(y-y_0)^2+(z-z_0)^2}$) is the distance from the source to the point of interest along the wavefront, which is also the magnitude of the vector. The parameter $\phi$ is the angle between the x-axis and the projection of the vector onto the x-y plane. The parameter $\theta$ is the angle between the z-axis and the vector along a vertical plane.

 

march-pol Figure 2 The steps taken to implement the fast marching in polar coordinates. Black circles imply computed traveltimes that are set because of their minimum traveltime value along the front. Gray circles constitute the front of the wave, and they are stored in the wavefront heap array and sorted from minimum to maximum traveltime value. The traveltime is schematically given by the size of the circle; the larger the radius, the greater the traveltime. The minimum is always extracted first from the heap array at each step, its traveltime is set (given a black circle), and all surrounding grid points that are not set are computed and put into the heap array. We precede until all grid points are computed and set.

Figure 2 shows schematic plots of the progress of this method along a 2-D polar coordinate grid. The source is computed initially and set to zero for all angles $\theta$. When stretched on a regular grid, all points at the surface r=0 are set to zero. These points are inserted in the wavefront array and sorted from minimum traveltime to maximum. In the case of the source grid point, the sorting step is unnecessary because all traveltimes are equal to zero. The minimum is then extracted, and the traveltimes for neighboring grid points are computed using the following relation:

$$\begin{eqnarray} max(D_{ijk}^{-r}\,t,0)^2 & + & min(D_{ijk}^{+r}\,t,0)^2+ \nonum... ...ijk}^{-\phi}\,t,0)^2 & + & min(D_{ijk}^{+\phi}\,t,0)^2=s^2_{ijk}, \end{eqnarray}$$ (1)

where D_ijk^-r is the derivative of traveltime with respect to r at grid point i, j, k, given by

$$D_{ijk}^{-r}\,t = \frac{t_{i,j,k} - t_{i-1,j,k}}{\Delta r},$$

and

$$D_{ijk}^{+r}\,t = \frac{t_{i+1,j,k} - t_{i,j,k}}{\Delta r}.$$

The $D_{ijk}^{\theta}$ and $D_{ijk}^{\phi}$ derivatives are slightly different, given by

$$D_{ijk}^{-\theta}\,t = \frac{t_{i,j,k} - t_{i,j-1,k}}{r \Delta \theta},$$

$$D_{ijk}^{+\theta}\,t = \frac{t_{i,j+1,k} - t_{i,j,k}}{r \Delta \theta},$$

$$D_{ijk}^{-\phi}\,t = \frac{t_{i,j,k} - t_{i,j-1,k}}{r \sin \theta \Delta \phi},$$

and

$$D_{ijk}^{+\phi}\,t = \frac{t_{i,j+1,k} - t_{i,j,k}}{r \sin \theta \Delta \phi}.$$

The traveltime t_ijk and the slowness s_ijk correspond to the grid point that is being updated. Solving for t_ijk, the time at a new grid point, requires solving a quadratic equation that has two solutions. We choose the solution that reduces to $t_{ijk}=t_{i-1,j,k}+\frac{\Delta r}{v}$, or $t_{ijk}=t_{i+1,j,k}+\frac{\Delta r}{v}$, when the wavefront travels radially. We do the same for all other unset points surrounding the initial extracted point. Each newly computed grid point is added to the wavefront array, and by using a highly efficient heap method the sorting from minimum to maximum traveltime is done promptly. Then we extract the minimum once again, update all neighboring live grid points, and so on.

Using fast algorithms (heap sorting), the sorting part of the fast marching algorithm can be maintained at a computational cost proportional to $\log N$, where N is the number of grid points in the computational domain. As a result, the cost of the eikonal solver is roughly proportional to $N\log N$. This efficiency feature is maintained in spherical coordinates as well. The cost of transforming traveltime information from spherical coordinates to Cartesian ones is proportional to N, and thus is less significant.

 

vel-top
Figure 3 The salt structure displayed from the top. The vertical section (left) is the inline slice through the 3-D salt velocity model and horizontal section (right) is a depth slice through the same velocity model.

 

vel-bot
Figure 4 The salt structure displayed from the top (left) and from the bottom (right). The vertical section is the inline slice through the 3-D salt velocity model and horizontal section is a depth slice through the same velocity model.