Next: Computation of amplitudes
Up: ALGORITHM
Previous: Slowness model
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 r_{0} 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 u_{i,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 u_{i,j} > 0
 3.
 u_{i1,j} < 0 and
 4.
 u_{i1,j} < 0 and u_{i,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 t_{i,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(i1)))
sign2 = 0.5*(1.+sign(1.,u(i)))
hatw1 = sign1*min(w(i1),s(i))+(1.sign1)*s(i1)
hatw2 = sign2*min(w(i),s(i1))+(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