THE BASIC FINITE-DIFFERENCE SCHEME

The Engquist-Osher scheme used by Van Trier and Symes is an interesting balance between computational speed and precision. Although the scheme introduces approximations in the calculation of the travel time field that can be adverse for accurate calculations of differential operators (which involve the local derivative of the traveltime field), the scheme is more stable than a straightforward finite-difference scheme. In spite of all the instability problems one may encounter using this algorithm, all of the simpler alternatives I have tried are stable only for much smaller grid sizes.

The fundamental algorithm described by Van Trier and Symes (1989) is based on the eikonal equation

 

u²+v²=s² (1)

where

$$\left \{ \begin{array} {l} u={\partial{t(x,z)} \over \partia... ...\\ \\ v={\partial{t(x,z)} \over \partial z} \end{array} \right.$$

and where s(x,z) is the 2-dimensional slowness model and t(x,z) is the travel time field. The second equation used is the equality of the partial derivatives of the fields u(x,z) and v(x,z),  

$${\partial u \over \partial z}={\partial^2 t \over {\partial x \partial z}} ={\partial v \over \partial x}.$$ (2)

In cylindrical coordinates the eikonal equation becomes  

$${u^2 \over r^2}+v^2=s^2$$ (3)

where

$$\left \{ \begin{array} {l} u={\partial{t(r,\theta)} \over \p... ...{\partial{t(r,\theta)} \over \partial r} \\ \end{array} \right.$$

and the mixed partial derivatives are  

$${\partial u \over \partial r}={\partial^2 t \over {\partial \theta \partial r}} ={\partial v \over \partial \theta}.$$ (4)

The finite difference implementation of equations (1), (2) and in cylindrical coordinates (3), (4) is based on advancing the computational front for the functions u(x,z) and v(x,z). The traveltime field t(x,z) is found subsequently by integrating the function $v(r,\theta)$ with respect to r. In cylindrical coordinates the scheme is based on using equation (4) to compute the values of $u(r,\theta)$ on a new computational front of constant radius, by using the values of $u(r,\theta)$ and $v(r,\theta)$ from the previous computational front. Equation (4) becomes  

$$u(r+\Delta r,\theta)=u(r,\theta)+{{\Delta r} \over {\Delta \theta}}\Delta v(r,\theta).$$ (5)

Starting with a constant velocity condition in the immediate vicinity of the source location $(u(r,\theta)=0, v(r,\theta)=s(r,\theta))$,we can advance the computational front using the finite-difference equation (5). Once the values of $u(r,\theta)$ are known, the values of $v(r,\theta)$ can be computed using the eikonal equation:  

$$v(r,\theta)=\sqrt{s^2(r,\theta)-{{u^2(r,\theta)} \over r^2}}.$$ (6)