Theory

Traveltime tomography is the procedure for reconstructing the earth's velocity model from picked traveltimes. These traveltimes can be measured from seismic data associated with a variety of source-receiver configurations, and can be used to extract velocity models of varying spatial resolution.

Figure 2.1 shows 3 types of traveltime data that can be used in exploration geophysics: refraction traveltimes measured as first breaks in CSP gathers, transmission traveltimes in crosswell data (Ivansson, 1986; Bishop et al., 1985), and reflection arrivals in CDP data. In this lecture we will confine ourselves to refraction tomography and crosswell tomography where we pick first arrival traveltimes from a suite of Common Shot Point gathers.

  

Figure 2.1: Three types of traveltime data: CMP refraction, CMP reflection and crosswell transmission traveltimes.

Classical traveltime tomography uses raypaths to model traveltimes and so assumes that the data are of high frequency, i.e., the wavelength is at least several times smaller than the spatial variations of the velocity model. In this case, the traveltime measurements are connected to the model via the traveltime integral:

 

$$\int_{raypath} \frac{1}{c(x,z)} ds.$$ t(x,z) = (2.1)

The traveltime integral says that the traveltime for energy to propagate from source to receiver is computed by integrating the weighted slowness s(x,z)=1/c(x,z) along the raypath that connects the source and receiver.

There are 7 steps to inverting traveltime data:

1.

Pick the traveltime t_i from the ith raw record, and construct the Mx1 traveltime vector ${\bf t}$ comprised of M traveltime picks. A finite-difference solution to the eikonal equation can be used to quickly compute these traveltimes (Vidale, 1988; Qin et al., 1992). Discretize the earth model into N slowness cells, where each ith cell has an unknown constant slowness s_i.

  

Figure 2.2: Cell model of an earth discretized into 16 cells. These slowness values form the Nx1 slowness vector s.

2.

Denote the segment length of the ith ray in the jth cell by l_ij. In this case the traveltime integral reduces to a summation of weighted slownesses, where the weight is the segment length of the rays:

 

$$\sum_{j}l_{ij} s_j.$$ t_i = (2.2)

The traveltime equations in 2.2 can be recast in matrix-vector notation as ${\bf L' s' = t'}$, where ${\bf L'}$ is the MxN raypath segment matrix, with element values equal to l_ij and s' is the actual slowness model.

3.

For some initial guess slowness model $\bf s$ we have ${\bf t = L s,}$ which can be used to form the perturbed set of traveltime equations:

 

$${\bf t - t'} {\bf L s - L' s' },$$ = (2.3)

where $\bf L$ is the raypath matrix associated with the guessed $\bf s$ model. Appendix A shows that ${\bf L'}$ can be approximated by $\bf L$ to second-order in the perturbation parameter ${\bf ds = s-s'}$, i.e., $\bf {L' = L + O(ds^2)}$. To first-order in the perturbation parameter equation (3) becomes:

 

$${\bf dt} {\bf L ds }.$$ = (2.4)

This is known as the linearization step where the relationship between the model and data is linearized by equation (4); this implicitly assumes that the segment lengths do not change when the slowness model is slightly perturbed from the actual slowness model.

4.

There are usually many more equations than unknowns, so that equation (4) can be solved by minimizing the sum of the squared errors subject to a length constraint:

 

$$\epsilon ({\bf Lds-dt})^T ({\bf Lds-dt}) + \alpha {\bf ds^T ds },$$ = (2.5)

where $\alpha$ is a small positive scalar. Minimizing equation (5) w/r to the components of $\bf s$ yields the damped least squares solution:

 

$${\bf ds} ({\bf [L^T L]} + \alpha {\bf I})^{-1} {\bf L^T dt},$$ = (2.6)

where I is the NxN Identity matrix.

5.

The slowness model is updated by

 

$${\bf s} { \bf s} + \beta {\bf ds} ,$$ = (2.7)

where $\beta$ is the step length. If the starting model is very close to the actual model so that the traveltime equations are truly linear, then $\beta = 1$ and equation (7) is the final answer. We are usually not so lucky so we use equation (7) as our new starting model and repeat steps 3-5 to get a new updated model. This procedure, known as a non-linear Gauss-Newton method, is repeated until convergence in the model and the traveltime residual.

Figure 3 depicts the contours of a linear and a non-linear error surface.

  

Figure: Contours of a linear and a non-linear error surface. If the relationship between the data and model are approximately linear and the system of equations is somewhat well-conditioned then the damped Gauss-Newton direction given by equation (6) with $\epsilon =0$ points to the global minimum of the error surface, as shown on the RHS figure. Contrarily, if the error surface is highly non-linear (i.e., the traveltime equation that connects the data and model are highly non-linear) then the Gauss-Newton direction might point to a local minimum, as shown on the LHS figure. Hence further iterations (combined with smoothing to jump out of the local minimum) are needed to get to the global minimum.

Equation (6) demands the direct calculation of the inverse matrix for the NxNmatrix ${\bf L^TL}$. This is too expensive in practice (N³algebraic operations/iteration) because there are a large number of N unknowns in a typical tomography problem.

To avoid this direct calculation we resort to an avoid this indirect or iterative solution method which only requires O(N²) algebraic operations/iteration. We will examine the steepest descent method with the understanding that a more effective method such as a Conjugate Gradient-type method is used in practice.