rayinvr.doc

利用fortran写的实现2D射线追踪的源代码！运行于linux下！

            DOCUMENTATION FOR RAYINVR AND RELATED PROGRAMS
            ----------------------------------------------

Date:       Dec 1993
Version:    1.3

Written by: C.A. Zelt
Address:    Geological Survey of Canada
            1 Observatory Crescent
            Ottawa, Ontario, Canada  K1A 0Y3
phone:      [phone]   
fax:        [phone]
e-mail:     [e-mail]


            Table of Contents
            -----------------


RAYINVR              
     Program Description            
     Files                         
     Input Parameters             
         Plotting parameters     
         Axes parameters        
         Ray tracing parameters
         Inversion parameters 
     Additional Notes        
     Array Sizes            
     Warning and Error Messages 
Other Programs                 
     VMODEL 
     DMPLSTSQR               
     RAYPLOT


RAYINVR: a program to trace rays in 2-D  media for rapid forward 
         modelling and  inversion of  refraction and  reflection  
         travel times

                        Program Description
  
  A 2-D (x,z) isotropic medium is assumed. The velocity  model is 
composed  of  a  sequence  of   layers  separated  by  boundaries  
consisting of  linked  linear segments  of  arbitrary dip.  Layer  
boundaries must  cross  the  model  from  left  to  right.  Layer  
thicknesses may be reduced to zero to model pinchouts or isolated 
bodies. The velocity within a layer is defined by velocity values 
specified at arbitrary x-coordinates along the  top and bottom of 
the layer. The x-coordinates at which  layer boundaries and upper 
and lower velocities are specified can  be completely general and 
independent within and  between layers.  Velocity discontinuities  
across layer  boundaries are  allowed but  not required.  For the  
purposes of ray tracing, the model is automatically broken up into 
an irregular network of  trapezoids, each with  dipping upper and 
lower  boundaries  and  vertical   left  and  right   sides.  The  
velocities at  the four  corners  of the  trapezoid  are used  to  
interpolate a  velocity field  within the  trapezoid so  that the  
velocity  varies  linearly  along  its   four  sides.  Therefore,  
horizontal as well as vertical velocity gradients may exist within 
a trapezoid. A simulation of smooth  layer boundaries is possible 
in which the incident and emergent ray angles are calculated using 
the slope of the smoothed boundary.
  The source(s) may be positioned anywhere in the model  and rays 
may be directed any angle. The receivers are always assumed to be 
at the top of  the model. Both  P- and S-wave  propagation can be 
considered including (multiple)  conversions. A  unique Poisson's  
ratio may be assigned to each trapezoid  of the model. Refracted, 
reflected and head waves may be  traced, each possibly containing 
multiple and/or surface reflections and conversions. Ray take-off 
angles are determined automatically by the  program for those ray 
groups specified by the user using an iterative shooting/bisection 
search mode. Ray tracing is performed  by numerically solving the 
ray tracing equations for 2-D media, a pair of first order ODE's, 
using a Runge Kutta method. The ray  step length is automatically 
adjusted at each  step to  maximize efficiency  while maintaining  
accuracy. Travel  times are  calculated by  numerical integration  
along ray paths using the  trapezoidal rule. A plot  of the model 
and all rays traced may be produced along  with a plot of reduced 
travel time versus distance for the observed and calculated data.
  The partial derivatives  of travel time  with respect to  those 
model  parameters   selected   for   adjustment  are   calculated   
analytically  during  ray   tracing;  these   parameters  include   
velocities and  the  vertical  position  of boundary  nodes.  The  
travel times  correspond to  any ray  paths which  can be  traced  
through the  model, being  either  first or  later arrivals.  The  
travel time residuals with respect to the  observed data are also 
calculated. The travel times and partial derivatives are linearly 
interpolated to the observed seismogram locations since two-point 
ray tracing is not required. The partial derivatives  and  travel
time residuals are output and used later as input to the  program
DMPLSTSQR which  updates  the model  parameters  by applying  the  
method of damped least-squares to the linearized inverse problem.
  The model parameterization is well  suited to the inversion  of 
refraction/reflection data  since realistic  earth models  can be  
represented by a minimum  number of model  parameters, i.e.,. the 
number and position  of parameters  specifying each layer  can be  
adapted to the data's subsurface  ray coverage. Layer boundaries, 
including the  surface,  may  be  horizontal (one  parameter)  or  
consist of numerous  straight line segments.  A layer may  have a 
constant velocity (one parameter) or the velocity structure may be 
defined by many upper and lower  layer velocity points. Different 
velocity points may be specified above and below a layer boundary 
if a velocity discontinuity is required across the boundary, or a 
single row of  velocity points may  be specified if  an interface 
with no discontinuity is  needed. The vertical  velocity gradient 
or layer thickness may be fixed in all or  part of a layer during 
the  inversion  if   there  is   insufficient  ray   coverage  to   
independently determine an upper and lower layer velocity or layer 
thickness.
                              Files

Executable file: RAYINVR

Input file: r.in contains program input parameters in five parts:
   (1)  the  PLTPAR  namelist  contains   parameters  related  to  
   plotting
   (2) the AXEPAR namelist contains parameters related to axes
   (3) the  TRAPAR namelist  contains parameters  related  to ray  
   tracing
   (4)  the  INVPAR  namelist  contains   parameters  related  to  
   inversion
   (5) after skipping three lines (column headings), the velocity 
   model is specified as follows:
     (a) the  layer number,  the x-coordinates  (km)  of a  layer 
        boundary entered  from  left to  right  (format: I2,  1X,  
        10F7.2)
     (b) the z-coordinates (km) of a layer boundary corresponding 
        to the x-coordinates listed above (format: 3X, 10F7.2)
     (c) a  0,  1  or  -1  for  each  z-coordinate  listed  above  
        depending on whether (1)  or not (0)  partial derivatives 
        are to be calculated for a particular boundary node or the 
        boundary depth is to be determined by fixing the thickness 
        of the  layer above  at that  x-coordinate  (-1) (format:  
        3X, 10I7)
     (d) the layer number, the  x-coordinates (km) of the  points 
        at which the  upper layer  velocity is  specified entered  
        from left to right (format: I2, 1X, 10F7.2)
     (e) the upper layer  P-wave velocities (km/s)  corresponding 
        to the x-coordinates listed above (format: 3X, 10F7.2)
     (f) a 0  or 1 for  each velocity listed  above depending  on 
        whether (1)  or not  (0)  partial derivatives  are  to be  
        calculated for a particular velocity (format: 3X, 10I7)
     (g) the layer number, the  x-coordinates (km) of the  points 
        at which the  lower layer  velocity is  specified entered  
        from left to right (format: I2, 1X, 10F7.2)
     (h) the lower layer  P-wave velocities (km/s)  corresponding 
        to the x-coordinates listed above (format: 3X, 10F7.2)
     (i) a 0, 1 or -1 for each velocity listed above depending on 
        whether (1)  or not  (0)  partial derivatives  are  to be  
        calculated for a particular velocity or the lower velocity 
        is to  be  determined  be  fixing the  vertical  velocity  
        gradient at that x-coordinate (-1) (format: 3X, 10I7)
  The above sequence  of nine lines  is repeated  for each model  
   layer, the  top-most layer  specified  first, the  bottom-most  
   last, and is ended by specifying the  bottom layer boundary of 
   the model as  in (a) and  (b) above.  If the number  of points  
   defining a boundary or the upper or  lower velocity of a layer 
   must exceed 10 then the points can be continued onto subsequent 
   lines of the  file as  follows: line  (b), (e)  or (h)  of the  
   particular parameter to be extended is modified to include a 1 
   in the second column so the complete format of the line becomes 
   I2, 1X, 10F7.2. The sequence of three lines (a)-(c), (d)-(f) or 
   (g)-(i) is then repeated  as many times as  is necessary using 
   the same format described above.

Input file: v.in contains  the velocity model in  the same format  
   as described in part (5) above for r.in.

Input file: f.in contains  the "floating" reflecting boundaries;
   these are interfaces with no associated velocity discontinuity
   (hence the word "floating") at which rays can reflect in 
   addition to the layer boundaries; the file format is similar to 
   that for v.in:
   (1) the number of nodes defining the floating reflector
       (format: I2)
   (2) the x-coordinates (km) of the nodes defining the floating
       reflector (format: 3X, <number of nodes>F7.2)
   (3) the z-coordinates (km) of the nodes defining the floating
       reflector (format: 3X, <number of nodes>F7.2)
   (4) a 0 or 1 for each node listed above depending on whether (1)  
       or not  (0) partial derivatives  are  to be  calculated for 
       this particular node (format: 3X, <number of nodes>I7)
   Lines (1) through (4) are repeated for each floating reflector 

Input file:  tx.in  contains  the  observed travel  time-distance  
    pairs in the following format:
   (1) the  x-coordinate  (km)  of  the  shot  point,  1  if  the 
      receivers are to the right of the  shot point or -1 if the  
      receivers are to the left, 0, and 0 (format: 3F10.3, I10)
   (2)  the  x-coordinate   (km)  of  the   observed  data,  the   
      corresponding unreduced travel time (s), the estimated 
      uncertainty of the travel time pick  (s), and a  non-zero 
      integer used to identify the type of arrival  to allow for 
      the appropriate comparison with the rays traced 
      (format: 3F10.3, I10)
   Line (2) is repeated for each  pick corresponding to the shot  
   point in line (1).  The sequence (1) and  (2) is repeated  for 
   each shot point of the data  set. The file is terminated with  
   the following line:
   (3) 0, 0, 0, -1 (format: 3F10.3, I10)

Output file: tx.out contains the  calculated travel time-distance 
    pairs in the same format as described above for tx.in.

Output file:  r1.out contains  summary information  for  each ray  
    traced including  the  shot  number,  take-off and  emergent  
    angle, range, reduced  travel time,  number of  points (step  
    lengths) defining the ray and the ray code.

Output  file:  r2.out  contains  detailed   parameters  for  each  
    trapezoid of the velocity model, a one-dimensional equivalent 
    (average) velocity  model and  summary information  for each  
    point of each ray traced.

Output file: i.out contains the  partial derivatives, travel time 
    residuals and travel time uncertainties used for input by the 
    program  DMPLSTSQR.  The  top  of   the  file  contains  the  
    following information:
       (1) the number of travel time picks used and the number of 
         model parameters  for  which  partial derivatives  were  
         calculated
       (2)  for  each  model   parameter  involved,  an   integer  
         identifying the type of model parameter (1 for boundary 
         or 2  for velocity),  the  current value  of  the model  
         parameter, and the  estimated uncertainty of  the model  
         parameter
       (3) the  partial derivatives,  travel  time residuals  and 
         corresponding uncertainties in matrix form
       (4) the number of travel time picks used
       (5) the RMS travel time residual
       (6) the normalized chi-squared

Output file: p.out  contains all plot  commands for the  run used  
    for input by the program RAYPLOT.

Output file: n.out contains the namelist parameter values.

                         Input Parameters

1) Plotting parameters (PLTPAR namelist):

a) Switches  (usually 0 = off, 1 = on):

iroute - equals 1 to plot to the screen, 2 to create a postscript 
   file, 3 to create a plot file for the VERSATEC plotter, or 4
   to create a colour postscript file; if iroute does not equal 1 
   there is no plotting to the screen  (default: 1)

iseg - create a Uniras segment(s) (default: 0)

iplot - generate the plot  during the run (1), or  write all plot 
   commands to the file p.out (0), or do both (2) (default: 1)

imod - plot model boundaries (default: 1)

ibnd - plot vertical model boundaries (default: 1)

idash - plot model boundaries as dashed lines (default: 1)

ivel - plot the P-wave (1) or  S-wave (-1) velocity values (km/s) 
   within each trapezoid of the model (default: 0)

iray - plot  all rays traced  (1) or  only those which  reach the  
   surface (or reflect off the floating reflector for ray groups for 
   which frbnd>0) (2); rays traced in the search mode are not 
   plotted unless irays=1 (see below) (default: 1)

irays - plot the rays traced in the search mode (default: 0)

irayps - plot  the P-wave  segments of  ray paths as  solid lines  
   and the S-wave segments as dashed lines (default: 0)

idot - plot a symbol at each point (step length) defining each ray 
   path (default: 0)

itx -  plot the calculated travel time-distance  values   as lines 
   (1), symbols of height symht (2),  symbols of height symht 
   interpolated  at  the  observed receiver locations contained in
   the  file  tx.in (invr=1 must be used, if not, itx=2 is used) (3); 
   or as residuals measured from the oberved data in tx.in (invr=1
   must be used) (4); no calculated travel times are plotted if itx=0 
   (default: 1)