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.

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:

The traveltime integral says that the traveltime for energy to propagate from source to receiver is computed by integrating the weighted slowness

There are 7 steps to inverting traveltime data:

- 1.
- Pick the traveltime
*t*_{i}from the*ith*raw record, and construct the*Mx*1 traveltime vector 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*Nx*1 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:

The traveltime equations in 2.2 can be recast in matrix-vector notation as , where 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
we have
which can be used to form the perturbed
set of traveltime equations:

where is the raypath matrix associated with the guessed model. Appendix A shows that can be approximated by to second-order in the perturbation parameter , i.e., . To first-order in the perturbation parameter equation (3) becomes:

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:

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

where*I*is the NxN Identity matrix. - 5.
- The slowness model is updated by

where is the step length. If the starting model is very close to the actual model so that the traveltime equations are truly linear, then 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 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 *NxN*matrix
.
This is too expensive in practice
(*N*^{3}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*^{2}) 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.