usage.txt

射线追踪程序

Short guide for using the tomo3d package (S閎astien Judenherc, Mar 2002).

1/ install the package + dependencies
install libsismoutil *without* iasp support (see cvs on nilgiri.u-strasbg.fr)
install iasp
install libsismoutil *with* iasp ("make clean" and "./configure --with-iasp=...")
install libeve
install libolo
install tomo3d

For each of these packages, ./configure --help gives you the usefull
informations. Please read.

If you have a computer with more than one processor (or even if you have
not), you may want to benefit the use of multithreading implemented in
the longest part of the inversion: building the matrix and lsq
inversion. If so, consider using the *thread* options in the configure
script in tomo3d.

2/ tomo3d input files (only 3 files)
====================================

"name" stand for the name of your inversion (no blank).

Model file: name.lay
====================
example:

#----------------------------------cut here
48.0 -2.0
20  90.0    10  90.0    10.0    5.80 1
50  20.0    50  20.0    20.0    5.80 1
50  25.0    50  25.0    30.0    8.10 1
50  25.0    50  25.0    35.0    8.10 1
50  30.0    50  30.0    35.0    8.10 1
50  30.0    50  30.0    35.0    8.10 1
50  30.0    50  30.0    35.0    8.10 1
#----------------------------------cut here

First line: center of the grid lat/lon. The grid is built around this
location.
Next lines:
Nw-e    Dw-e    Ns-n    Ds-n    thick.  dummy   damp_fact

Na-b: number of blocks along direction a-b.
Da-b: block size along direction a-b.
thick: layer thickness
dummy: was the layer mean velocity.
damp_fact: multiplication factor for damping value in this layer. 1 is a
good value !

Remarks:
-There are N/2 blocks along each direction around the center of the grid.
-The grid must be large enough to contain all the rays. The size of the
initial grid does not affect the inversion time as only the inverted
blocks are considered.
-The first layer is a special one: one block for 1 station BUT the layer
must contain all the rays (ok, it is stupid but...). So use LARGE
blocks. You must give at least as many blocks as stations (giving more
is not a problem).

Configurations file: name.cnf
=============================

Defines env. variables. The format is the usual unix-like *rc format:
NAME="the value you want to use"

#----------------------------------cut here
NORM=-b
#type of normalization (b: baseline, E: average by event).
#b: only the total average (over ALL residuals) is removed, should be
#small. zero if the 1D model is correct :)
#e: event normalization. The average for each event is computed and the
#residuals for this event are corrected.

REBLK1=0
MAXHIT1=100
#REBLK1=1 if you want to vertically cut in 4 equal size blocks the
#blocks visited by more than MAXHIT1 rays.

REBLK2=0
MAXHIT2=200
#the same for the second pass. (the new block may be very small...)

MAXLAY=3
#at which layer do we stop the REBLK pass

damp="-D 80"
#damping value

MINH=10
#minimum number of ray per block

MINR=5
#minimum number of ray per event

SIG2=0.1
#a priori uncertainty on data (s)

INC=5.0
#incremental length for ray tracing in the 1D model

MAXDEPTH=700
#mximum depth for the 1D ray tracing (at least as deep as the model !)

#VREF="-c model-vs-iasp.dat fact"
#defines the correction applied to the iasp model to obtain your
#preferred 1D model for raytracing.
#model-vs-iasp.dat: "r Vpert(not in %!)" on each line
#fact indicates the multiplying factor applied:
#delta_Vfinal(Z)=Viasp*pert(6371-Z)*fact
#----------------------------------cut here


Data file: name.dat
===================

1 line for 1 residual
example (read 1 line):

S111   48.170   -2.923  10   0  -20.000 -177.765 587.0
151.59 349.76 7.27   921237399.490 1128.965 1128.453  -0.165
2.314 921237399.eve Fiji Islands reg

S111:   station code (4 chars)
48.170: station lat (deg, north positive).
-2.923: station long (deg, east positive)
10:     station elev (m)
0:      event number (must be the same for all the res. for the same
        phase of the same event). Must change at each new event. The
        value is not important: it must change i.e. you can go from 11
        to 20 without problem. You can also start at 10. Just try to
        increase the event number.
-20.000:    event lat (deg, north positive)
-177.765:   event long (deg, east positive)
587.0:  event depth (km)
151.59: great circle arc (deg)
349.76: back azimuth (deg)
7.27:   azimuth (deg)
921237399.490:  origin time (epoch)
1128.965:   theoretical travel time (s)
1128.453:   observed travel time (s)
-0.512: residual (s)
2.314:  slowness (s/deg)
921237399.eve:  event name 1 (no blank!)
Fiji Islands reg:   event name 2 (free format, 16 char).
                    Don't store any important information here.

3/ running the inverse script
=============================

run ./inverse without arg:
    first arg (name): the name of the inversion/dataset

    second arg: offset&average
        OA1: offset&average 1 (half block)
        OA2:                2 (1/3rd block)
        1: no offset&avg. (single inversion)

    third arg: ray tracing
        RT: ray tracing is to be done (always)
        NRT: no ray tracing, uses the on stored in name.ray

    fourth arg: param type
        B:  block
        N:  nodes (untested for a while... please don't use... VERY slow)

4/ running the inverse3d script
===============================
For 3D ray-tracing/iterative inversion (long, long, long!)
see the script for parameters!
Check that -W option for 3draytrace is correct for you parameter
geometry (gaussian filter width for the interpolation: using the average
inter-parameter distance is probably a good idea).

5/ output files
===============
name.blk:   block file (not usefull)
name.geo:   MAIN RESULT FILE
name.hit:   hitcount file (not usefull)
name.inv:   internal use
name.log:   log file
name.num:   internal use
name.out:   data_index res_before res_after
name.par:   internal use
name.ray:   binary ray tracing file, not usefull
name.res:   resolution matrix (binary file).
name.inf:   information matrix (binary file) if requested, same as res.
name.sta:   station file (for a map ?)
name.sys:   linear problem file (internal use)
name.var:   variance file (appended!)
name.xyz:   same as .geo but in km in stead of degrees.
name-blk.gmt:   block geometry as a multisegment file for psxy -M

Usefull files:
--------------
-name.geo:
    lon, lat, depth, pert(%), stderr(%), "layer-"[0-N], hit, res_diag

    pert:   velocity perturbation
    res_diag:   disgonal term of resolution matrix.

-name.res:
    binary file containing the complete resolution matrix:
    i(float),j(float),R(float), use with -b option in gmt. i&j start at 0!
    example: xyz2grd -R0/818/0/818 all.res -b -I1 -N0 -Gres.grd
        For an inversion with 819 blocks.

-name.inf:
    same as name.res for the information matrix (if requested, see lsq
    options and modify inverse script if needed).

-name-blk.gmt:
    block geometry for use with psxy -M.
    To extract layer X from this file, use:
        "grep ' X$' name-blk.gmt"
    then plot it with psxy -M

-name.log:
    to check computation time etc.

-name.var:
    variance file. Usefull to fix the damping value as this file is
    appended by lsq.
    1/ Run ./inverse once with the first damping
    2/ run lsq with various damping values
    3/ plot the resulting variance reduction, trace(R) etc.
    Format:
    data_var_init data_var_final model_variance trace(R) damping Ndata Nmodel

-name.gmt:
    if request by raytrace (with -G option). Ray path as 3D multisegment
    lines (huge file, long to dump!!!). X(km), Y(km), Z(km). Z positive
    downward.

-final.geo: this file is generated by the "inverse" script. It is in the
    same format as name.geo but contains the result of the
    offset&average inversions.

There are also output files from the 3D ray tracing. They are in a human
readable format. They are:

-name.har:
    amplitudes of the perturbations applied to the rays for each
    wavelength and starting/final travel time.

-name.r3t:
    same as name.gmt (see -g in 3draytrace)

-name.r3d:
    same as name.ray (binary file for internal use)

6/ plotting the results
=======================
Modify the scripts mkgrdpl.tri and trplall according to your data set or
write your own scripts !

7/ cross sections
=================
crossgeo is designed to caomput cross section in the file name.geo (in
case of use of offset&average, remeber to copy final.geo.sort to
name.geo).
invoking crossgeo without argument dump a short help.
The output file (name-c.geo) is:

offset_along_profile(km) 0.0 depth(km) pert(%) stddev "layer"(layer number) hit-count diag_term_of_res

8/ synthetics
=============
8.1/ isotropic model
--------------------
You have to change the source code to add your own velocity function.
Here is a steb-by-step explanation:

in synthe.c, in function parseVfunc, you first have to add the entry for
yout function (called MyVelocityIso):

  if (!strcmp(s1,"MyVelocityIso")) {
     tmp->fun=MyVelocityIso;
     return(tmp);
  }

There are 2 arguments for the velocity functions, they are usually used
as amplitude-length (name A and P) but of course, they can be used for
any other parameter.

Then , still in synthe.c, you have to add the body of MyVelocityIso
somewhere at the beginning of the file (before the parser):

float MyVelocityIso(float x, float y, float z, float A, float P)
{ if (z...) return(...);
  if ((z...)&&(x...)) etc;

  /* For harmonic test:
  return(1.0+A*sin(x*.../P)*sin(y*.../P)*sin(z*.../P));
  */

  return(1.0);
}

x, y and z are cartesian coordinates (km) relative to the ref point defined
in name.lay
The value returned by MyVelocityIso is a factor for the velocity
reference field e.g.  0.95 for -5%, 1.05 for +5% ...

Then, "make" rebuilds synthetrace and you can now use MyVelocityIso:

synthetrace -n name -g 0.1 -Z 300 -p name.dat -i 1 -F MyVelocityIso 5 30
5 and 30 are used as A and P in MyVelocityIso.

8.1/ ANisotropic model
----------------------
ugly...
in synthe.c, add 
extern float MyVelocityAniso(float x, float y, float z, float A, float P);
at the with the other similar definitions.

Add your function to the parser:

  if (!strcmp(s1,"MyVelocityAniso")) {
     tmp->fun=MyVelocityAniso;
     initaniso(); /* !!!!!!!!! THIS WAS NOT IN THE ISOTROPIC EXAMPLE !!! */
     return(tmp);
  }

In syntheaniso.c,

float MyVelocityAniso(float x, float y, float z, float A, float P)
{ struct aniso_t *aniso;
  A*=100.; /* because ANISOTX=100.0*(Vmax-Vmin)/Vavg; */

  if (...tests on x, y, z to define anisotropic volumes ...) {
     aniso=anisoNord;
  }

  if (...other tests ? ...) {
     aniso=anisoSud;
  }

  if (...des tests sur x, y, z...) {
     return(1.0); /* if you have an ISOtropic region ... */
  }

  /* propagation vector... */
  geoph2xyz(RAYAZ*180.0/M_PI,RAYINC*180./M_PI,
            &(prop[0]),&(prop[1]),&(prop[2]));

  /* Christoffel matrix */
  m=christoffel(aniso,prop,m);

  /* eigenvalues */
  eigen(m,3,veloc,pol,3);

  /* normalisation (remember, ANISOTX is in % but we have already done
   * A*=100 */
  return(1.0+A*((sqrt(veloc[0])/VISO)-1.0)/ANISOTX);
  /* VISO and ANISOTX are computed elsewhere */
}

anisoNord and anisoSud are the two different anisotropic structure thant
can be defined. The are set up in initaniso(), according to azNord,
dipNord, azSud, and dipSud which define azimuth and dipping a the
low-velocity symetry axis for domaines called Nord and Sud... which can
be anywhere of course !

If you need mode (combining isotropic and anisotropic velocity fields,
more than two domains), consider using synthetrace several times with th
"-k" option (Keep residual and correct them with new model):

synthetrace -p data.dat ... -F riri 1 2 > 1.dat
synthetrace -k -p 1.dat ... -F fifi 1 2 > 2.dat
synthetrace -k -p 2.dat ... -F loulou 1 2 > 3.dat

WARNING: -k REMOVES the new residual from the original one as it was
designed to correct real data from an apriori model.