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:  

$$\tau^2_r+{1 \over r^2}\tau^2_\theta=m^2(\theta,r),$$ (17)

where subscripts r and $\theta$ denote the partial derivatives with respect to r and $\theta$, respectively. The initial condition is

$$\tau(\theta,r)=m(\theta,r)r \ \ \ \ \ \ \hbox{for}\ r \le r_0,$$

where r₀ 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 $\tau_{i,j}$ to denote the traveltime at grid point $(\theta_i, r_j)$. Then, in the jth step of the traveltime extrapolation, we know the traveltime $\tau_{i,j}$ for all i and we want to find the traveltime $\tau_{i,j+1}$ 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  

$$u_{i,j}={1 \over r_j}{\partial \tau \over \partial \theta} ={\tau_{i+1,j}-\tau_{i,j} \over r_j\Delta \theta}.$$ (18)

The orthogonal component of the traveltime gradient can be computed from the eikonal equation as follows:  

$$w_{i,j}=\sqrt{s^2_{i,j}-u^2_{i,j}}.$$ (19)

These two components determine the propagation direction of the local plane wave. If we define  

$$\hat{w}_{i,j}={\tau_{i,j+1}-\tau_{i,j} \over \Delta r}$$ (20)

then we can use the principles of Podvin and Lecomte's scheme (1991) to calculate $\hat{w}_{i,j}$ in the following four cases:

1.

$u_{i-1,j} \ge 0$ and $u_{i,j} \le 0$

$$\hat{w}_{i,j}=\min(w_{i-1,j},w_{i,j},s_{i-1,j},s_{i,j})$$

2.

$u_{i-1,j} \ge 0$ and u_i,j > 0

$$\hat{w}_{i,j}=\min(w_{i-1,j},s_{i,j})$$

3.

u_i-1,j < 0 and $u_{i,j} \le 0$

$$\hat{w}_{i,j}=\min(w_{i,j},s_{i-1,j})$$

4.

u_i-1,j < 0 and u_i,j > 0

$$\hat{w}_{i,j}=\min(s_{i-1,j},s_{i,j}).$$

The appearance of slowness in these formulas indicates that the possible refractions along the constant $\theta$ boundary are taken into account. The ``min'' function ensures that case gives the minimum traveltime. The causality of the extrapolation requires that  

$$w_{i-1,j} \gt s_{i-1,j}{\Delta r \over \sqrt{\Delta r^2+(r_j... ...j} \gt s_{i,j}{\Delta r \over \sqrt{\Delta r^2+(r_j\theta)^2}}.$$ (21)

If we adapt the step size $\Delta r$ of the extrapolation to the propagation direction of the local plane wave and choose

$$\Delta r = \min_{\hbox{all}\ i}[ r_j\Delta \theta \left\vert{w_{i,j} \over u_{i,j}}\right\vert],$$ (22)

then the inequalities in equation (21) are guaranteed. Solving equation (20) for t_i,j+1 gives

$$t_{i,j+1}=t_{i,j}+\hat{w}_{i,j}\Delta r \ \ \ \ \ \ \ \ \hbox{for all}\ i,$$ (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 $\hat{w}_{i,j}$ 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.