stratinv.c.bac

su 的源代码库

#include "stratInvMain.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"                                                                        ",
" stratInv - Bayesian waveform inversion code for stratified (v(z))      ",
"            medium. This codes operating by maximizing a Gaussian       ",
"            a posteriori probability distribution along the lines       ",
"            described by the references below. A nonlinear version of   ",
"            conjugate gradient method is employed in the calculation.   ",
"            This code is implemented in a distributed fashion using the ",
"            PVM library, according to the master-slave topology. The    ",
"            input parameters for this code are as follows.              ",
"                                                                        ",
" model=[model]  Elastic model file. This is the initial model used in   ",
"                optimization. This is an ASCII file that should contain ",
"                one row per layer of the model in the following format: ",
"                                                                        ",
"                Thick(km)  Rho(g/cm^3)  Vp(km/s)  Qp  Vs(km/s)  Qs      ",
"                                                                        ",
"                Where:                                                  ",
"                Thick: Layer thickness                                  ",
"                Rho:   Layer density                                    ",
"                Vp:    P-wave velocity                                  ",
"                Qp:    P-wave quality factor                            ",
"                Vs:    S-wave velocity                                  ",
"                Qs:    S-wave quality factor                            ",
" data=[data]    Input seismograms for the inversion, in SU format.      ",
" scale=[1]      Scale factor for input data.                            ",
" datacovar=[0]  If set to 1 data covariance matrix required             ",
" corrData=[covD]Inverse of data covariance matrix.                      ",
" impedance=[0]  If set to 1, inversion is carried out in the impedance -",
"                density domain. Otherwise it is carried out in the      ",
"                velocity - density domain.                              ",
" p=[1]          P-wave (velocity or impedance) active.                  ",
" s=[1]          S-wave (velocity or impedance) active.                  ",
" r=[1]          Density active.                                         ",
" prior=[0]      If set to 1 prior information required.                 ",
" corrP=[covP]   Inverse of the P-wave model covariance matrix.          ",
" corrS=[covS]   Inverse of the S-wave model covariance matrix.          ",
" corrR=[covR]   Inverse of the density model covariance matrix.         ",
" targetbeg=[0]  Beggin of the target depth interval (km).               ",
" targetend=[1]  End of the target depth interval (km).                  ",
" t1=[0]         Beggin of the time misfit window (s).                   ",
" t2=[1]         End of the time misfit window (s).                      ",
" mutefile=[NULL]Muting information to mute SYNTHETIC data. This is as   ",
"                ASCII file with the following format for each row:      ",
"                time(s)     trace index                                 ",
"                ...                                                     ",
"                The program will linearly interpolate between times.    ",
" recfile=[NULL] File with receiver offsets. The format is one offset    ",
"                line. If none is specified (default), the acquisition   ",
"                is end on                                               ",
" directwave=1   Direct wave flag. If zero direct wave is not generated  ",
" r1=[0]         Minimum source-receiver offset (km)                     ",
" nr=[48]        Number of receivers                                     ",
" dr=[0.025]     Receiver spacing (km)                                   ",
" zs=[0.001]     Source depth (km)                                       ",
" u1=[0.0002]    Initial horizontal slowness (s/km)                      ",
" u2=[1]         Final horizontal slowness (s/km)                        ",
" nu=[1000]      Number of slownesses                                    ",
" f1=[2]         Initial frequency for modeling (Hz)                     ",
" f2=[50]        Final frequency for modeling (Hz)                       ",
" dt=[0.004]     Time sampling interval (0.004)                          ",
" tmax=[8]       Maximum modeling time (s)                               ",
" F1=[0]         EAST-WEST point force coordinate                        ",
" F2=[0]         NORTH-SOUTH point force coordinate                      ",
" F3=[1]         VERTICAL point force coordinate                         ",
" Hanning=[1]    Final trace is convolved with a hanning window          ",
" wu=[5]         Percentual of slowness interval that is windowed        ",
" ww=[5]         Percentual of frequency interval that is windowed       ",
" fr=[1]         Reference frequency used in absorption model            ",
" tau=[50]       Wrap-around attenuation parameter                       ",
" nProc=[1]      Number of processors used in the computation            ",
" nfreqproc=[0]  Number of frequencies allocated to each processor       ",
"                The default value means that the frequency              ",
"                components will be evenly split among the slave         ",
"                processors                                              ",
" verbose=[0]    =1 print useful information                             ",
"                                                                        ",
" stratInv outputs an ASCII file containing information regarding the    ",
" optimization work. This file also contains the subsurface model at     ",
" each iteration of the optimization.                                    ",
"                                                                        ", 
" Reminder: YOU SHOULD INSTALL PVM (version 3.3.9 or higher) BEFORE THE  ",
"           INSTALLATION OF THIS CODE.                                   ",
"                                                                        ",
"                                 Author: Wences Gouveia - 96-08-01      ",
"                                 Center For Wave Phenomena              ",
"                                 Colorado School of Mines               ",
NULL};   
/************************ end self doc ***********************************/ 

void main (int argc, char **argv)
{
   /* declaration of variables */
   FILE *fp;                     /* file pointer */
   char *dataFile = " ";         /* data file */
   char *auxChar = " ";          /* auxiliar character */ 
   char *muteFile = " ";         /* muting file */
   char *orientation = " ";      /* orientation of recordings */
   char *modelFile = " ";        /* elastic model file */
   char *corrDataFile = " ";     /* data covariance file */
   char *corrModelFile[3];       /* model covariance file */
   int maxIter;                  /* maximum iterations at conjugate gradient*/ 
   int notDone;                  /* convergence flag */
   int i, j, k;                  /* counters */
   int nMutePoints;              /* number of traces used in misfit interp */
   float a, b;                   /* used in interpolation */
   float dZ;                     /* layer thickness within target zone */
   float depth;                  /* current depth used in defining limits */
                                 /* for Frechet derivatives */
   float fR;                     /* reference frequency */
   float limZ[2];                /* target interval (Km) */
   float tMod;                   /* maximum modeling time */
   float slope;                  /* gradient . search direction */
   float *grad0, *grad1;         /* gradient vectors */
   float *search;                /* search direction */   
   float *auxP;                  /* auxiliar variable */
   float *t1Aux, *offAux;        /* used in muting interpolation */
   float oF0, oF1;               /* objective function value */
   
   /* allocing for orientation */
   orientation = malloc(1);
   
   /* complex Zero */
   zeroC = cmplx(0, 0);

   /* getting input parameters */
   initargs(argc, argv);

   /* seismic data and model parameters */
   if (!getparstring("model", &modelFile)) modelFile = "model";
   if (!getparstring("data", &dataFile)) dataFile = "model.su";
   if (!getparfloat("scale", &scaleData)) 
      err("Specify scale for field data\n");
   if (!getparstring("corrData", &corrDataFile)) corrDataFile = "corrData";
   if (!getparint("prior", &PRIOR)) PRIOR = 1;
   if (!getparint("datacov", &DATACOV)) DATACOV = 1;
   if (!getparint("vp", &vpFrechet)) vpFrechet = 1;
   if (!getparint("vs", &vsFrechet)) vsFrechet = 1;
   if (!getparint("rho", &rhoFrechet)) rhoFrechet = 1;
   if (!getparint("verbose", &info->verbose)) info->verbose = 1;
   numberPar = vpFrechet + vsFrechet + rhoFrechet;
      
   if (PRIOR)
   {
      if (vpFrechet)
      {
	 if (!getparstring("corrVP", &corrModelFile[0])) 
	    corrModelFile[0] = "corrVP";
      }
      if (vsFrechet) 
      {
	 if (!getparstring("corrVS", &corrModelFile[1])) 
	    corrModelFile[1] = "corrVS";
      }
      if (rhoFrechet) 
      {
	 if (!getparstring("corrRHO", &corrModelFile[2])) 
	    corrModelFile[2] = "corrRHO";
      }
   }
   
   if (!getparstring("orientation", &orientation)) orientation[0] = 'Z';
   if (orientation[0] == 'z' || orientation[0] == 'Z')
   {
      info->vertical = 1; info->radial = 0;
   }
   else
   {
      info->vertical = 0; info->radial = 1;
   }
   
   if (!getparfloat("dz", &dZ)) dZ = .5;
   if (!getparfloat("targetbeg", &limZ[0])) limZ[0] = 0.5; 
   if (!getparfloat("targetend", &limZ[1])) limZ[1] = 1.0;

   /* geometry */
   if (!getparfloat("r1", &info->r1)) info->r1 = 0.25;
   if (!getparint("nr", &info->nR)) info->nR = 48;
   if (!getparfloat("dr", &info->dR)) info->dR = .025;
   if (!getparfloat("zs", &info->zs)) info->zs = .001;
   if (!getparfloat("F1", &info->F1)) info->F1 = 0;
   if (!getparfloat("F2", &info->F2)) info->F2 = 0;
   if (!getparfloat("F3", &info->F3)) info->F3 = 1;

   /* modeling */
   if (!getparstring("receiverfile", &auxChar)) auxChar = " ";
   sprintf(info->recFile, "%s", auxChar);    
   if (!getparfloat("u1", &info->u1)) info->u1 = 0.0;
   if (!getparfloat("u2", &info->u2)) info->u2 = 1.;
   if (!getparint("directwave", &info->directWave)) info->directWave = 1;
   if (!getparint("nu", &info->nU)) info->nU = 1000;
   if (!getparfloat("f1", &info->f1)) info->f1 = 2;
   if (!getparfloat("f2", &info->f2)) info->f2 = 50;
   if (!getparfloat("dt", &dt)) dt = 0.004;
   if (!getparfloat("tmod", &tMod)) tMod = 8;
   if (!getparfloat("tau", &info->tau)) info->tau = 50;

   /* misfit information */
   if (!getparfloat("t1", &t1)) t1 = 0;
   if (!getparfloat("t2", &t2)) t2 = tMod;   
   if (!getparstring("mutefile", &muteFile)) MUTE = 1; else MUTE = 0;
  
   if (MUTE) 
   {
      fp = fopen(muteFile, "r");
      if (fp == NULL) err("Can't open mute file\n");
      nMutePoints = 0;
      while (fscanf(fp, "%f %f\n", &aux, &aux) != EOF) nMutePoints++;
      t1Aux = alloc1float(nMutePoints);
      offAux = alloc1float(nMutePoints);
   
      rewind(fp);
      i = 0;
      while (fscanf(fp, "%f %f\n,", &t1Aux[i], &offAux[i]) != EOF) i++;
      fclose(fp);
   }

   /* DD
   for (i = 0; i < nMutePoints; i++)
      fprintf(stderr, "%f %f\n", t1Aux[i], offAux[i]);
   exit(-1);*/
   
   if (!getparint("hanning", &info->hanningFlag)) info->hanningFlag = 1;
   if (!getparfloat("wu", &info->percU)) info->percU = 10; info->percU /= 100;
   if (!getparfloat("ww", &info->percW)) info->percW = 25; info->percW /= 100;

   /* distributed computation */
   if (!getparint("nproc", &nProc)) nProc = 1;
   if (!getparint("nfreqproc", &info->nFreqProc) || nProc == 1) 
      info->nFreqProc = 0;

   /* checking number of receivers */
   fp = fopen(info->recFile, "r");
   if (fp != NULL)
   {
      info->nR = 0;
      while (fscanf(fp, "%f\n", &auxm1) != EOF) info->nR++;
   }
   fclose(fp);

   /* building muting area */   
   t1Mute = alloc1float(info->nR);
   recArray = alloc1float(info->nR);
   
   /* reading receiver configuration */
   fp = fopen(info->recFile, "r");
   if (fp == NULL)
   {
      /* standard end-on */
      if (info->verbose) fprintf(stderr, "No receiver file available\n");
      for (i = 0; i < info->nR; i++)
      {
	 recArray[i] = info->r1 + i * info->dR;
      }
   }
   else
   {
      if (info->verbose) fprintf(stderr, "Reading receiver file\n");
      for (i = 0; i < info->nR; i++)
      {
	 fscanf(fp, "%f\n", &recArray[i]);
      }
   }
   fclose(fp);

   if (MUTE)
   {
      /* correcting end points */
      /*offAux[0] = 1;
      offAux[nMutePoints - 1] = info->nR;
 
      t1Mute[0] = t1Aux[0];
      t1Mute[info->nR - 1] = t1Aux[nMutePoints - 1];*/

      for (k = 1, i = 0; i < nMutePoints - 1; k = j, i++)
      {