BOMCRATR


BOMCRATR (Bureau of Mines Curved Ray Tomographic Reconstruction) is a curved ray tomographic program for seismic travel times that uses the SIRT method of reconstruction. It contains a ray tracing algorithm adapted from one used by Michael Weber. Input can be prepared with the program PREBMCRA. Manual entry of the pixel grid is tedious.

Written by Daryl Tweeton, Bureau of Mines, 5629 Minnehaha Ave. South, Minneapolis, MN 55417-3099. Phone (612) 725-4670. Was revised Feb 2, 1993. (Increased MAXTTM to 1850.) Fixed velocity nodes are listed in matrix form. Additional revisions Dec 26, 1993. Corrected calling of SORT3 for receivers on top border and changed STOP criteria in RAYTRC. Changed F1.NE.F2 comparisons for guarding against division by zero with ABS(F1-F2).GT..00002.

The following Bureau of Mines standard disclaimer applies. The Bureau of Mines expressly declares that there are no warrenties expressed or implied which apply to the software contained herein. By acceptance and use of said software, which is conveyed to the user without consideration by the Bureau of Mines, the user expressly waives any and all claims for damage and/or suits for or by reason of personal injury, or property damage, including special, conseqential, or other similar damages arising out of or in any way connected with the use of the software contained herein.


Operation of Program
Operation of Program Label numbers
Read or create input file 100- 999
Calculate tomographic corrections and velocities 1000-3660
Perform optional averaging or smoothing 3670-3900
Write output files 3910-6999
End program, possible error messages 9000-9999


Description of input file

Units for length, travel times, and velocities must be self-consistent. Usual units are length in meters, time in milliseconds, velocity in m/ms. For labeling in the program, sources are assumed to be on left and receivers on right side. Pixel 1 is at top left. Pixel 2 is next pixel in top row.

All data are in list-directed format. Input numbers must be separated by a space or comma. They do not need to be in particular columns. For readability, some variable names are longer than six characters, but the first six characters are distinct for portability.

HEADER
One line of comments for the top of output file, first column blank.

ITRMAX MSMTH MSTRT
ITRMAX = maximum number of iterations allowed. Set MSMTH=1 for smoothing after each iteration. Weighting is 20, 4, 1 for pixel, its nearest neighbors, and corner neighbors, respectively. Set MSTRT=1 to force all rays to be straight.

IPDAT IPITR IPSEG IRAYGRF IRAYTBL
IPDAT=1 causes printing of source, rec positions and travel times. IPITR=1 causes printing of ray information in *.RAY for all iterations instead of just the final one. IPSEG=1 causes printing of segment lengths for each path and pixel. (It makes a very large output file, useful mostly for debugging.) IRAYGRF=1 creates output file of simple graphical display of ray paths. IRAYTBL=1 creates output file of x,z values at pixel sides along ray path. Any value other than 1 suppresses the corresponding printing or display.

XTOPL ZTOPL XBOTL ZBOTL
XTOPL and ZTOPL are the x and z coordinates of the top left node. XBOTL AND ZBOTL are the x and z coordinates of the bottom left node. They should be set so that the sources are not exactly on a grid border or on pixel sides or corners. Grid borders are assumed to be straight for this calculation, but border node positions can be adjusted to account for curved boreholes.

XTOPR ZTOPR XBOTR ZBOTR
XTOPR and ZTOPR are the x and z coordinates of the top right node. XBOTR AND ZBOTR are the x and z coordinates of the bottom right node.

NNODH NNODV NLATOP NLABOT NFXLFT NFXRHT
NNODH=number of columns of nodes (counting horizontally). NNODV=number of rows of nodes (counting vertically). If desired pixel height is EPIXV, set NNODV=1+(ZBOTR-ZTOPL)/EPIXV. NLATOP and NLABOT are the number of laterally invariant rows at top and bottom of the grid. The velocity in each of those rows is set equal to the average for that row after each iteration. NFXLFT and NFXRHT are the number of columns at the left and right sides of the grid in which the velocities are fixed at the initial values.

IDSOR(ISOR) XSOR(ISOR) ZSOR(ISOR) (One line for each ISOR = 1 to NSOR)
End with -99 -99. -99. IDSOR(ISOR) is the identification number of the source position. They do not need to be sequential. XSOR AND ZSOR are X and Z coordinates of source positions. Z increases with depth. Source positions should not be exactly on a border or the side of pixel.

IDREC(IREC) XREC(IREC) ZREC(IREC) (One line for each IREC = 1 to NREC)
End with -99 -99. -99. IDREC(IREC) is the identification number of the receiver position. XREC and ZREC are X and Z coordinates of receiver positions. Receiver positions must be on one of the four borders, not in a corner.

ADJTTM
ADJTTM is added to all experimental travel times.

IDSORT(IEXPTM) IDRECT(IEXPTM) EXPTT(IEXPTM) (One line for each IEXPTM = 1 to NEXPTM)
End with -99 -99 -99. IDSORT and IDRECT are the identification numbers of the source and receiver positions for the travel time EXPTM(IEXPTM)

SMFPLL VELMIN VELMAX INVOPT
SMFPLL is the lower limit for the sum of fractions of path lengths in a pixel (segment length in pixel/length of path from source to rec). If the sum is less than or equal SMFPLL, then the velocity for the pixel is displayed as ####, and is not calculated.
VELMIN and VELMAX are the low and high limits for the calculated velocities. If a velocity is out of that range, it is set = the limit. If no limit is desired, set the limit to -99.
INVOPT specifies the option for inputting the initial guesses of velocities.

If INVOPT = 1, input is uniform and the format is: VEL(1) BOMCRATR sets initial velocity = VEL(1) for every node.
If INVOPT = 2, input is row by row and the format is: IROWV VEL(IROWV) (One line for each row, NNODV lines) IROWV is the row number. Rows can be input in any order. BOMCRATR sets initial velocities = VEL(IROWV) for all nodes in that row. This option facilitates varying the pattern of initial guesses to test the uniqueness of a solution.
If INVOPT = 3, input is node by node and the format is: IROWV (VEL(I),I=1,NNODH) (One line for each row, NNODV lines) IROWV is the row number. VEL(I) is the initial velocity for row IROWV, column I. This option is used when resuming a run using a data file as input.

I XNOD(I) ZNOD(I) (One line for each node.)
I - number of node, starting with I=1,monotonically increasing, maximum set to MAXNOD (see parameter definition) XNOD(I),ZNOD(I) - x,z-coordinates of node I
Terminate this set with -99 -99. -99.

XNOD(J) every second position, at top of grid ZNOD(INODV) LGLVAR(J + (INODV-1)*NNODH) J=1,number of node columns One line for each row of velocity node points. XNOD(J) every second position, at bottom of grid
This matrix of logical variables is for conveniently specifying whether the velocity should be fixed to the initial value (F) or allowed to vary (T). XNOD and ZNOD are displayed around the matrix to facilitate selecting the velocity points to fix. They are read into an array that is temporarily not in use, and not used as new data. Fixing velocities with NFXLFT or NFXRHT overrides T in this matrix.

JPIX KNOD(JPIX,1) KNOD(JPIX,2) KNOD(JPIX,3) LJ(JPIX,1) LJ(JPIX,2) LJ(JPIX,3) (One line for each triangular pixel.)
JPIX - number of triangles,starting with JPIX=1,monotonically increasing, maximum set to MAXTRI (see parameter definition).
KNOD(JPIX,1), KNOD(JPIX,2), and KNOD(JPIX,3) are the nodepoints (I) of this triangle, in mathematically positive direction, i.e. counterclockwise.
LJ(JPIX,1), LJ(JPIX,2), and LJ(JPIX,3) are indices of neighboring triangles of triangle JPIX. LJ = -32111,-32222,-32333,-32444 indicate left, bottom, right, top borders of grid, respectively. LJ is given in counter- clockwise direction. LJ(JPIX,1) is the index of the neighboring triangle across the triangle side with the nodes KNOD(JPIX,1) and KNOD(JPIX,2), etc.
Terminate with -99 -99 -99 -99 -99 -99 -99


Files
Files Default name Unit
Input (output from PREBMCRA). Output files all have the same root filename supplied by the user and these extensions supplied by the program. BOMCRATR.IN LUIN
Time and path length for each ray going to a receiver for each iteration. BOMCRATR.RAY LURAY
Summary output from last iteration, can be used as input for contouring program TOMOPLOT, supplied with BOMCRATR. BOMCRATR.SUM LUSUM
Output file in same format as input file, with new velocities, for resuming run. BOMCRATR.RES LURES
Optional text mode graphs of ray paths at final iteration. BOMCRATR.GRF LUGRF
Optional table of x,z values at edge each pixel at final iteration. BOMCRATR.TBL LUTBL

Description of some program variables

There are NANG equally spaced take-off angles from the source between ANG1TRY and ANG2TRY. Angles are taken from positive z-axis. They range from -180 degrees (up) to -90.degrees (left) to 0 degrees (down) to +90 degrees (right) to +180 degrees (up). ANG1FAN and ANG2FAN are the angles that give sufficient coverage, and are used for the final ray tracing.

ANG - takeoff angle at source. Q - ray parameter at ray endpoint (horizontal slowness) i.e., at LJ=negative (if JNXTPIX=0 => Q=0.) TQX,TQZ - x,z-component of unit-tangent to ray at ray endpoint.

JNXTPIX = 0 bad ray , not used; = negative, ray ended at border.

XQ,ZQ - x,z-coordinates of ray endpoint. VQ - velocity at ray endpoint. T - travel time to ray endpoint.

XTOLER and ZTOLER are the distances within which receivers must be of the border to be counted along that border. They are set to XTOLMUL and ZTOLMUL times the distance between receivers 1 and 2.

parameter:
MAXNOD - max. # of nodes allowed.
MAXNDH - max. # of horizontal nodes allowed
MAXTRI - max. # of triangles allowed.
MAXRAY - max. # of rays allowed per source.
MAXSOR - max. # of source positions.
MAXREC - max. # of receiver positions.
MAXTTM - max. # of experimental travel times.
MAXSEG - max. # of ray path segments in each ray.

The program attempts to trace rays to both sides of all receivers, so receivers should not be on the corners of the model. The program adjusts the input fan of rays and interpolates travel times for specified receiver positions on the receiver plane.

An 'active' receiver has at least one ray on each side. For each active receiver a travel time is calculated with straight-line interpolation on the two closest enclosing rays. This interpolated travel time is output. Plotting data for the ray closest to each active receiver is output optionally to PATHGRF.GRF and PATHTBL.TBL.
© William P. Clement 2000