Cross-borehole geometry: Shots are assumed to be in the left well, receivers in the right well. Currently set up with wells on the left and right sides of a rectangular grid of boxes.
VSP geometry: Shots are assumed to be at the top or bottom of the model, receivers in a well on the right edge of the model.
Cross-borehole and VSP data can be used together.
In each Gauss-Newton iteration, traces rays through current model and generates the derivative matrix of travel times with respect to slowness parameters using R.T. Langan's program RAYTRC. Incorporates a priori information using what is called here the "omega model". Regularizes ill-conditioned linear systems using penalty terms (penalizing horizontal and/or vertical derivatives). Uses IRLS (iteratively reweighted least squares) to solve the linearized problem in a hybrid L1/L2 sense. Allows for different hybrid "norms" on data, omega, horizontal derivative, and vertical derivative parts of vector. Solves each linear least squares sub-problem using the Paige and Saunders LSQR algorithm, a variant of the conjugate gradient algorithm, with preconditioning. Uses continuation method to decrease some penalty terms to improve accuracy in a robust fashion. Stopping criteria for both the continuation steps (several nonlinear Gauss-Newton iterations) and the overall continuation procedure are built in; the diagnostics to make these decisions have been corrected to be consistent with the hybrid "norm".
The important subroutines are:
Version of program: September 15, 1989. K.P.B.Efficient May 1988 version of R.T.L. raytrace program; adjustments to shot/receiver logic in DATSET, RAYS, ONESHT, FAN, ONERAY, DIVIDE, and INTERP, December 1988.
R.T.L. adjustments to logic concerning IDXRAY in subroutines RAYS, ONESHT incorporated, March 1989.
To run program: submit xtomo. The input files must be in certain places. The slowness model must be in file fort.10. The shot/receiver geometry file, if used, must be file fort.16. The travel time data are in file fort.12.
Once these files are created and properly named, execute xtomo. Remember that many changes to the parameters are necessary to successfully run the program. xtomo outputs several files which are discussed below.
The model used as the true solution for both the error estimates and to construct the omega model is read in from the latest version of slowmod.dat (fort.10 -- see subroutine MODEL.
The weighting terms (omhdi(m)) used to indicate the spatial variation of the penalty term for deviating from the omega model are either constructed or read in from the latest version of omweight.dat, according to the flag 'iomhdi' -- see subroutine MODEL
The model used as the starting (gray) model is either constructed from the omega model or read in from the latest version of graymod.dat, according to the flag 'igrflg' -- see subroutine MODEL
The shot/receiver geometry is either constructed in DATSET or read in from the latest version of geom.dat (fort.16, according to the flag 'igeom' -- see subroutine DATSET.
The traveltime data are obtained from one of two sources
according to the parameter itrvtm:
Generally, for the flags iomhdi, igeom, itrvtm, etc.:
At least two output files are generated:
See cbtvsp.jcl (xtomo.f?) for the disposition of other files: trvtm_out.dat (travel times), geom.dat (shot/receiver geometry), omweight.dat (weights), cbtvsp.tb1, and cbtvsp.tb2. (These might be fort.30 and fort.32).
fort.30 and fort.32 contain a list of parameters
and the true model. Sample output files for fort.30
and fort.32 are available. The difference between
the files is after the listing of the parameters.
In fort.30, the headings are:
STEP DATA RMS NRAY NCOMP CPURT FRBNRM ICNT OMPEN HDPEN VDPEN DAT L2 OMEGL2 HDL2 VDL2 TRESID L2 ALPHA
In fort.32, the headings are:
STEP ICNT NSTEP NITN ISTOP ACOND CPUPS DP L2 B*DP L2 THETA DP RMS MODEL RMS DPFRAC LEVCNT
To change models:
Statements in the source code to change are enclosed by $$$$$; they include:
The specific changes to the FORTRAN source code are documented in another file.
raytrc SUBROUTINE FLOW DIAGRAM ***** means commented
out or removed for now
---------------------------
I I----MODEL
I MAIN/RAYTRC I----DATA/DATSET
I I----INIT
I I----RESULT
---------------------------
I I
I *****I**************
----- * -------+ *
RAYS * |PLTINF|--SHTPLT *
----- * -------+--FANPLT *
I ********************
I
-------------------------------------------------
I ONESHT I
-------------------------------------------------
I I
I I
I----------I I-------I--FANINF*****
I I--SHRINK I I +------+-STAT
STAT--I ONERCV I--PATHFO I FAN I--|DIVIDE|-FANINF*****
I I--BIGMAT I I +------+
I----------I--STORE***** I-------I--STAT I
I I I
I I_____ I
I I I
I I I
STAT----INTERP I I
I I ___I
I I I
ONERAY ONERAY
I I
I I
------------- ------------
DERIV---I ONEPXL I DERIV---I ONEPXL I
------------- ------------
I I I I I I
I I I I I I
HYPTN SIDE TOPBOT HYPTN SIDE TOPBOT
The main program was changed substantially, but still includes calls to MODEL, DATSET, INIT, and RAYS. In addition, INIT and RAYS are called in the loop which performs the successive Gauss-Newton iterations.
MODEL now defines the true model, the omega model (best a priori information), and the gray model (initial guess). Parameters like NWB, etc., are now defined in the main program.
DATSET now sets up shots and receivers so that their spacing is aspect/nshpb and aspect/nrcpb, respectively, so that box boundaries come half way between shot locations and half way between receiver locations; it also sets up a flag array ifsr, currently used to eliminate rays whose offset is out of range.
RAYS now also passes offmax and offmin to ONESHT.
ONESHT now eliminates rays whose (absolute value of) offset does not satisfy offmin <= offset <= offmax (using ifsr).
Actually the solution of the system of equations is more than just those obtained from ray tracing. In addition to these traveltime derivatives, there are several ways in which the problem is regularized (sometimes called Tikhonov regularization). One can put in an a priori model, and weighting to this a priori model (e.g., one may know velocities at the well bore, but not between the wells), vertical smoothing constraints, and horizontal smoothing constraints. Because the ray path coverage is limited, usualyy 3 or 4 times as much horizontal smooth is used, as compared to the vertical smoothing.
In addition, one can solve the problem in an L2 sense, and L1 sense, or a hybrid L1/L2 sense, with the user choosing the transition between L1 and L2 (usually an estimate of the noise level suffices for setting the L1/L2 transition). The L1/L2 problem is solved by using a method called 'iterative reweighted least squares'. The way the equation solver is set up(Paige-Saunders variant of conjugate gradients), there is not much of a penalty in time for doing this. In addition, the method takes advantage of the sparseness of the system.
Please feel free to ask questions.
Bob Langan