Program PRONTO implements an iterative tomographic inversion procedure for reconstructing a two-dimensional velocity model from measured first arrival times. Two key features of the technique are (i) use of a finite difference algorithm for rapid forward modeling of traveltimes, and (ii) inclusion of constraint information in order to restrict the nonuniqueness inherent in this large scale, nonlinear inverse problem.
The tomographic inversion procedure consists of four basic steps:
This four-step process is initiated with an estimate of the true slowness function, and is repeated until an acceptable match is obtained between observed and calculated traveltimes. A final velocity model is obtained by reciprocating the computed slowness model.
The inversion algorithm is not restricted to any particular source-receiver configuration. For example, it may be used with crosswell, VSP/RVSP, or surface-to-surface recording geometries. Moreover, undulating surface topography is easily accommodated.
A more complete description of the tomographic inversion procedure, including synthetic and field data examples, is given by Aldridge and Oldenburg (1993a, b).
The inversion algorithm works with a two-dimensional gridded slowness (i.e., reciprocal velocity) model defined at the points
x(i) = xmin + (i-1)*h, i=1,2,3,...,nx
z(j) = zmin + (j-1)*h, j=1,2,3,...,nz
where h is the grid interval. Note that this interval is the same in both spatial dimensions. In most applications, the x-axis is a horizontal axis and the z-axis is a depth axis (with depth increas- ing vertically downwards):
O | +x
--------------------------------->
|
|
|
|
|
|
|
+z \|/
Since the inversion algorithm is two-dimensional, all sources and receivers should be confined (approximately) to a plane. This is often the case in conventional line surveys, crosswell surveys, and VSP surveys.
The minimum horizontal grid coordinate (xmin) and the total number of grid points in the horizontal direction (nx) are automatically determined by the program so that all sources and receivers reside within the grid bounds. The analogous quantities for the vertical dimension of the grid (zmin and nz) are calculated in a similar manner. The user is responsible for specifying the grid interval (h).
The inversion algorithm is not sensitive to the unit of measure chosen for length. Either English or metric (or any other) unit may be used, as long as consistency is maintained.
The observed traveltimes used for the tomographic inversion are obtained from an external file. The file is read with the following simple FORTRAN code:
read (11,*) ns
do 1 i=1,ns
read (11,*) xs,zs,nr(i)
do 1 j=1,nr(i)
1 read (11,*) xr,zr,tobs,wait
where (xs, zs) are the source coordinates, (xr, zr) are the receiver coordinates, tobs is the observed traveltime, and wait is a (dimensionless) data weight value. The data in the input file should be organized in a `shot ordered' manner, with the total number of shot gathers equal to ns. Note that the number of first break pick times in each shot gather need not be the same, and is specified by the array values nr(i).
Each observed traveltime is tested against user specified lower and upper bounds. If it falls outside of these bounds, then the traveltime is discarded. This is a convenient mechanism for deleting bad picks from a dataset prior to performing a tomographic inversion. Bad picks may arise from noisy traces, dead traces, tailspreads, unpicked shot gathers, etc.
After the tomographic inversion is complete, either predicted (i.e., computed) traveltimes or traveltime residuals (i.e., observed minus predicted traveltimes) can be output. The format of the output file is identical to that of an input data file, described in 1) above. Hence, this file may be used to restart the algorithm.
The total number of output traveltimes (or residuals) can be less than the total number of input observed times. This will occur if certain input traveltimes are excluded from the tomographic inversion procedure via the limits described in point 2) above. Only those traveltimes (or residuals) that are actually used by the inversion algorithm will be output.
If traveltime residuals are output, then the minimum, maximum, mean, and rms values of the residuals are computed and written to standard output. Traveltime residuals provide useful diagnostic information for checking the quality of a tomographic inversion. Furthermore, the mean (or median) of the residuals associated with a particular source or receiver site can be used as a residual static correction for that location.
Finally, in order to conserve file storage space, the user may request that no traveltimes or residuals be written out.
The initial velocity model is obtained either from (i) the linear velocity formula, or (ii) an external file. If option (i) is selected, then an initial velocity model is internally generated by evaluating the linear velocity formula:
v(x,z) = v0 + ax*(x-x0) + az*(z-z0)
at the points
x(i) = xmin + (i-1)*h,
z(j) = zmin + (j-1)*h,
where h is the grid interval. In the above expression (x0,z0) are the coordinates of a reference point where the velocity equals the value v0, and (ax,az) are the horizontal and vertical components of the velocity gradient vector. These components can be expressed as
ax = a*cos(phi),
az = a*sin(phi),
where a is the magnitude of the gradient vector, and the angle phi defines the direction of the gradient relative to the +x direction. The angle phi is considered positive when it opens in a clockwise sense from the +x-axis. For example, the common situation where the velocity gradient vector points vertically downward is specified by phi=+90 degrees. The units of the gradient are m/s per m (v/m).
If option (ii) is selected, a gridded velocity model consisting of nx columns and nz rows is read into the program with the FORTRAN code:
do 1 j=1,nz
do 1 i=1,nx
1 read (12,*) x,z,v(i,j)
Velocity sample v(i,j) represents the velocity at location
x = xmin + (i-1)*h,
z = zmin + (j-1)*h,
where h is the grid interval.
After an initial velocity model is obtained, it is reciprocated to provide an initial slowness model for the tomographic inversion algorithm.
After the tomographic inversion is complete, the computed gridded slowness model is reciprocated to obtain a final velocity model. This model is written to an external file with the same format as an input velocity file described in point 4) above. Velocity sample v(i,j) represents the velocity at location
x = xmin + (i-1)*h,
z = zmin + (j-1)*h,
where h is the grid interval.
The raypath density associated with the final velocity model is written to an external file. The ray density measure used is a dimensionless quantity defined as the total raypath length per slowness cell, divided by the grid interval. A slowness cell is the square area bounded by four adjacent slowness grid points. The computed ray density is assigned to the center of the cell for plot purposes. A raypath density map consisting of (nx-1) columns and (nz-1) rows is written to the output file with the FORTRAN code:
do 1 j=1,nz-1
z=zmin+(j-0.5)*h
do 1 i=1,nx-1
x=xmin+(i-0.5)*h
1 write (16,*) x,z,rayden(i,j)
where rayden(i,j) represents the raypath density at location
x = xmin + (i-0.5)*h,
z = zmin + (j-0.5)*h,
and h is the grid interval.
Ray density provides useful diagnostic information for judging the reliability of a tomographic inversion. Portions of the velocity model that are not covered by any raypaths (typically near the corners) have a raypath density equal to zero (0.0).
In the Vidale (1988) scheme for solving the eikonal equation via finite differences, the horizontal and vertical grid intervals are equal. Thus, only one value for the grid interval needs to be specified. Typically, the grid interval is approximately equal to the source or receiver spacing for the survey.
As indicated in Section II above, the spatial limits of the grid are automatically determined from the minimum and maximum source and receiver coordinates contained in the input data file. This grid may be expanded, in any or all of the four coordinate directions (+x, -x, +z, -z), via user-specified extension distances.
There are two main reasons for expanding the grid. First, the raytracing algorithm performs better if sources and/or receivers are not located on edges of the velocity model. Second, in a surface seismic survey, the expansion distance in the +z direction controls the total depth to which a tomographic velocity model is determined. For such a survey, the depth should be chosen large enough so that raypaths can be successfully traced through the model to the far offset receivers.
The principal disadvantage of expanding the grid in the horizontal directions is that numerous cells may be added to the slowness model that have little or no raypath coverage. Hence, horizontal extensions should be kept small.
Each tomographic iteration consists of a forward modeling step (i.e., computation of traveltimes and residuals) and an inversion step (i.e., computation of a slowness model update). The maximum number of iterations must be specified.
Algorithm LSQR (Paige and Saunders, 1982; Nolet, 1987) is used to solve the large and sparse system of linear algebraic equations for an update to the slowness model. This algorithm is designed to seek the minimum norm, least squares solution of a set of linear equations. LSQR is iterative, and thus the number of iterations must be specified.
After traveltimes and raypaths are computed on each tomographic iteration, the root-mean-square value of the traveltime residuals is calculated. A traveltime residual is defined as an observed (i.e., picked) traveltime minus the corresponding predicted (i.e., computed) traveltime. Iterations terminate when the rms residual decreases below a user specified value. This value is usually chosen to be approximately equal to the estimated rms error in the picked times.
On every tomographic iteration, the traveltime residual associated with each source-receiver pair is tested against a user specified limit. If the absolute value of the residual exceeds this limit, then the equation corresponding to that residual is weighted by zero (0.0) in the system to be solved by LSQR. Thus, this equation has no influence on the least squares solution for the slowness model update vector.
During the course of the tomographic iterations, constraints can be applied to the slowness model. Constraints are used to stabilize the inversion and to introduce a priori geological or geophysical knowledge into the inversion procedure. The constraints may be applied:
The basic philosophy of the tomographic inversion procedure is to construct a slowness model that simultaneously satisfies the observed traveltime data and the specified constraints. The degree to which the constraint equations are fit is determined by user defined weight values. If a weight value equals zero, the corresponding constraint equations are ignored. Increasing the weight value implies that the constraint equations will assume a more important role in the inversion, at the expense of the data equations. Finally, in the limit of very large weight values, the data equations will be ignored and the inversion algorithm will construct a slowness model that satisfies the constraints only.
Unfortunately, there are no general rules for choosing numerical values for the constraint equation weights. Rather, values have to be determined by experimentation with the program. Experience with a few datasets indicates that values in the range 5 to 100 produce reasonable results.
The constraints are designed to produce a velocity model with features similar to those of a prescribed reference velocity model. Commonly, the reference model is identical to (i) the initial velocity model. However, this restriction is not mandatory. Thus, the reference model can also be obtained either from (ii) the linear velocity formula, or (iii) an external file. See the pertinent remarks under point 4) above.
A fuller description of the use of constraints in the inversion procedure is given by Aldridge and Oldenburg (1993a,b).
A rectangular smoothing filter can be convolved with the gridded slowness model between tomographic iterations. The horizontal and vertical dimensions of the filter are specified by the user. The weight distribution of the filter can be described as `tent-shaped': the edge weights are assumed to be 1.0, and the center weight (i.e., at the apex of the tent) is prescribed by the user. All weights are then normalized so that the filter response at wavenumbers kx = kz = 0 is unity. Hence, a constant slowness is passed unaltered by the filter. Note that a simple moving average filter can be obtained by setting the center weight equal to 1.0.
Filter edge effects are handled by conceptually extending the slowness model beyond the defined grid with the local gridpoint slowness values.
The filter response is zero-phase. The dimensions of the filter need not be integer multiples of the grid interval h. If the horizontal (vertical) dimension of the filter is smaller than two grid intervals (2*h), the the 2D filter degenerates into a 1D vertical (horizontal) line filter. If both filter dimensions are less than 2*h, then the filter request is ignored.
The seismic raypaths linking source-receiver pairs should be located beneath the earth's surface. This will usually occur if the velocity associated with grid points above the surface is sufficiently low. Ideally, the velocity at these points should equal sound speed in air (approximately 350 m/s or 1150 ft/s). However, a strong velocity contrast at the surface can create problems for the finite difference algorithm used to compute arrival times. Thus, the velocity at these grid points should only be `somewhat lower' than the velocity immediately below the ground surface.
If the `topograpy option' is activated, then the velocity at grid nodes above the earth's surface will be reduced. At horizontal coordinate x(i)=xmin+(i-1)*h, the velocity at these points will be set equal to a fraction of the velocity at the grid point just below the surface with the same x-coordinate. The fractional value is specified by the user and must be between 0.0 and 1.0.
The reduced velocity values are used only for forward modeling purposes (i.e., computation of traveltimes and raypaths) in order to ensure that raypaths are located below the earth's surface. When the final velocity model is written to an external file, the velocity at gridpoints above the topographic surface will be set equal to the (user specified) speed of sound in air.
If the tomographic survey is far from the earth's surface (as in many crosswell surveys) or there is no significant variation in the elevation of the surface (as in marine seismic surveys), then the topography option should not be activated. Use this option only if there is substantial elevation variation along a land seismic line.
Aldridge, D.F., and Oldenburg, D.W., 1993b, Two-dimensional tomographic inversion with finite-difference traveltimes: Journal of Seismic Exploration, 2, 257-274.
Nolet, G., 1987, Seismic wave propagation and seismic tomography: in Seismic Tomography, ed. by G. Nolet, D. Reidel Publishing Company, 1-23.
Paige, C.C., and Saunders, M.A., 1982, LSQR: an algorithm for sparse linear equations and sparse least squares: ACM Transactions on Mathematical Software, 8, 43-71.
Vidale, J., 1988, Finite-difference calculation of travel times: Bulletin of the Seismological Society of America, 78, 2062-2076.
FORTRAN code's header documentation
Sample run parameter file
Tips on running PRONTO
Sample reference velocity file
How to get PRONTO
| Last updated on December 22, 1999 by William P. Clement |