Next: Computation of amplitudes Up: ALGORITHM Previous: Slowness model

## Computation of traveltimes

A good finite difference scheme should incorporate the relevant physical law into its representation. Extrapolating traveltimes is like propagating wavefronts; Huygens' principle and Fermat's principle should be respected. Extrapolating the gradients of traveltimes is like extending rays; Snell's law should be respected. Because these two extrapolations simulate the same physical phenomenon, the algorithm to do one can be translated and used to do the other. The algorithm described below does two extrapolations simultaneously.

 vfdtime1 Figure 1 The grid and notations.

The eikonal equation in polar coordinates is as follows:
 (17)
where subscripts r and denote the partial derivatives with respect to r and , respectively. The initial condition is

where r0 defines a circle centered at the source point, within which the slowness is constant. Starting from this circle, we extrapolate the traveltime field in the radial directions. As indicated in Figure , I use to denote the traveltime at grid point . Then, in the jth step of the traveltime extrapolation, we know the traveltime for all i and we want to find the traveltime for all i.

According to Snell's law, the tangential component of the traveltime gradient to a boundary is continuous across the boundary. I define ui,j to be the tangential component of the traveltime gradient to the constant radius boundary. If we approximate the wavefront within a cell with a local plane wave, then we have
 (18)
The orthogonal component of the traveltime gradient can be computed from the eikonal equation as follows:
 (19)
These two components determine the propagation direction of the local plane wave. If we define
 (20)
then we can use the principles of Podvin and Lecomte's scheme (1991) to calculate in the following four cases:

1.
and

2.
and ui,j > 0

3.
ui-1,j < 0 and

4.
ui-1,j < 0 and ui,j > 0

The appearance of slowness in these formulas indicates that the possible refractions along the constant boundary are taken into account. The min'' function ensures that case gives the minimum traveltime. The causality of the extrapolation requires that
 (21)
If we adapt the step size of the extrapolation to the propagation direction of the local plane wave and choose
 (22)
then the inequalities in equation (21) are guaranteed. Solving equation (20) for ti,j+1 gives
 (23)
which completes one step of the extrapolation.

The implementation of this algorithm is quite simple and the program is fully vectorizable. For example, the computation of for all i can be coded as follows:

        do i = 1,n
sign1 = 0.5*(1.+sign(1.,u(i-1)))
sign2 = 0.5*(1.+sign(1.,-u(i)))
hatw1 = sign1*min(w(i-1),s(i))+(1.-sign1)*s(i-1)
hatw2 = sign2*min(w(i),s(i-1))+(1.-sign2)*s(i)
hatw(i) = min(hatw1,hatw2)
enddo


One problem that has not been addressed yet is the possibility of refractions along the constant radius boundaries, which may initiate wave propagations toward the source. In practice, such situation seldom happen. To make the algorithm more complete, I handle such situation by adopting Podvin and Lecomte's method in polar coordinates and performing backward extrapolation as necessary. After obtaining the traveltimes in polar coordinates, I map them to Cartesian coordinates through interpolation.

Next: Computation of amplitudes Up: ALGORITHM Previous: Slowness model
Stanford Exploration Project
12/18/1997