kdmig3d.c

用于石油地震资料数字处理

      if (!getparint("nyt",&gd->nyt))	  gd->nyt=ttr.ntr;

      gd->nxgd=ttr.selev;
      gd->nygd=ttr.gelev;

      /*number of sources*/
      if (!getparint("nxs",&gd->nxs))	gd->nxs=ttr.nhs;
      if (!getparint("nys",&gd->nys))	gd->nys=ttr.nvs;
      if (!getparint("nzs",&gd->nzs))   gd->nzs=ttr.duse;

      /*first source*/
      if (!getparint("ixsf",&gd->ixsf))	gd->ixsf=ttr.sdel;
      if (!getparint("iysf",&gd->iysf))	gd->iysf=ttr.gdel;
      if (!getparfloat("fzs",&gd->fzs)) gd->fzs=ttr.sdepth/1000.0;

      /*source spacing*/
      if (!getparint("ixsr",&gd->ixsr))	gd->ixsr=ttr.swdep;
      if (!getparint("iysr",&gd->iysr))	gd->iysr=ttr.gwdep;
      if (!getparfloat("dzs",&gd->dzs)) gd->dzs=ttr.ep/1000.0;
     
      fclose(ttfp);
      if ((ttfp=fopen(ttfile,"r"))==NULL)
            err("Can not open ttfile=%s",ttfile);
      efseeko(ttfp,(off_t) 0,SEEK_CUR);

      if (!getparfloat("dyoffset",&dyoffset)) dyoffset=99999;

      if (!getparint("ntrmax",&ntrmax)) ntrmax=100000;
      if (!getparint("pptr",&pptr)) pptr=100;

      gd->fxs=gd->fxgd+gd->ixsf*gd->dxgd;
      gd->dxs=gd->ixsr*gd->dxgd;
      gd->fys=gd->fygd+gd->iysf*gd->dygd;
      gd->dys=gd->iysr*gd->dygd;

      exs=gd->fxs+(gd->nxs-1)*gd->dxs;	/*last source in x in tttable*/
      eys=gd->fys+(gd->nys-1)*gd->dys;	/*last source in y in tttable*/
      ezs=gd->fzs+(gd->nzs-1)*gd->dzs;	/*last source in z in tttable*/
      ext=gd->fxgd+(gd->nxgd-1)*gd->dxgd;	/*last sample in x in tttable*/
      eyt=gd->fygd+(gd->nygd-1)*gd->dygd;	/*last sample in y in tttable*/

      if (!getparint("nxo",&gd->nxo)) gd->nxo=(gd->nxgd-1)*2+1;
      if (!getparfloat("fxo",&gd->fxo)) gd->fxo=gd->fxgd;
      if (!getparfloat("dxo",&gd->dxo)) gd->dxo=gd->dxgd*0.5;

      if (!getparint("nyo",&gd->nyo)) gd->nyo=(gd->nygd-1)*2+1;
      if (!getparfloat("fyo",&gd->fyo)) gd->fyo=gd->fygd;
      if (!getparfloat("dyo",&gd->dyo)) gd->dyo=gd->dygd*0.5;

      if (!getparint("nzo",&gd->nzo)) gd->nzo=(gd->nzs-1)*5+1;
      if (!getparfloat("fzo",&gd->fzo)) gd->fzo=gd->fzs;
      if (!getparfloat("dzo",&gd->dzo)) gd->dzo=gd->dzs*0.2;

      exo=gd->fxo+(gd->nxo-1)*gd->dxo;	/*last sample in x in output*/
      eyo=gd->fyo+(gd->nyo-1)*gd->dyo;	/*last sample in y in output*/
      ezo=gd->fzo+(gd->nzo-1)*gd->dzo;	/*last sample in z in output*/

      /*****************************************************
      The migration area should lie within the traveltime
      region

	      (fxgd,fygd,fzs)
		    _________________________________
		    |(fxo,fyo,fzo)		     |
		    |       __________________       |
		    |      |	   	      |      |
		    |      |		      |      |
 		    |      |		      |      |
		    |      |		      |      |
		    |      |__________________|      |
		    |		   (exo,eyo,ezo)     |
		    |________________________________|
					       (ext,eyt,ezs)
      *****************************************************/

      /***************************************************
      The apertures for x and y is the influential region
      of the current trace
      ***************************************************/
      if (!getparfloat("xaper",&gd->xaper)) 
       	    gd->xaper=gd->nxt*gd->dxgd/2.5;
      if (!getparfloat("yaper",&gd->yaper)) 
	    gd->yaper=gd->nyt*gd->dygd/2.5;
      if (!getparfloat("xoffsetmax",&xoffsetmax)) xoffsetmax=3000;
      if (!getparfloat("yoffsetmax",&yoffsetmax)) yoffsetmax=3000;

      /**************************************************
      For effective MVA, we need some nxoffset, nyoffset 
      about 10
      **************************************************/
      if (!getparint("nxoffset",&nxoffset)) nxoffset=1;
      if (!getparint("nyoffset",&nyoffset)) nyoffset=1;
      if (!getparfloat("fxoffset",&fxoffset)) fxoffset=0.;
      if (!getparfloat("fyoffset",&fyoffset)) fyoffset=0.;
      if (!getparfloat("dxoffset",&dxoffset)) dxoffset=99999;

      /*********************************************************
      Get information from the first header of data
      *********************************************************/
      if (!fgettr(infp,&tr)) 
	    err("Can't get first trace");
      gd->nt=tr.ns;	/*number of samples*/
      if (!getparfloat("dt",&gd->dt)) 
	    gd->dt=((double)tr.dt)/1000.0;  /*ms*/
      if (gd->dt<0.0000001) 
	    err("dt must be positive!");
      if (!getparfloat("ft",&gd->ft)) 
	    gd->ft=tr.delrt/1000.0;
      if (!getparfloat("dxm",&dxm))
            dxm=tr.d2;
      if (dxm<0.0000001) err("dxm must be positive");
      if (!getparfloat("angmax",&angmax))
            angmax=60.0*PI/180.0;
      if (angmax<0.0000001) err("angmax must be positive");
      gd->ntrimax=2*dxm*sin(angmax)/(gd->v0*gd->dt);
      if (gd->ntrimax<1) gd->ntrimax=1;

      if (!getparfloat("fmax",&gd->fmax)) gd->fmax=0.25/gd->dt;

      /*********************************************************
      Preparing the job print
      *********************************************************/
      if ( !getparstring("jpfile",&jpfile)) {
	    gd->jpfp=stderr;
      } else
	    gd->jpfp=fopen(jpfile,"w");

      fprintf(gd->jpfp,"\n");
      fprintf(gd->jpfp," Migration parameters\n");
      fprintf(gd->jpfp," ================\n");
      fprintf(gd->jpfp," ttfile=%s \n",ttfile);
      fprintf(gd->jpfp," \n");
      fprintf(gd->jpfp," nxt=%d fxgd=%g dxgd=%g\n",gd->nxt,gd->fxgd,gd->dxgd);
      fprintf(gd->jpfp," nyt=%d fygd=%g dygd=%g\n",gd->nyt,gd->fygd,gd->dygd);
      fprintf(gd->jpfp," nxs=%d fxs=%g dxs=%g\n",gd->nxs,gd->fxs,gd->dxs);
      fprintf(gd->jpfp," nys=%d fys=%g dys=%g\n",gd->nys,gd->fys,gd->dys);
      fprintf(gd->jpfp," nzs=%d gd->fzs=%g dzs=%g\n",gd->nzs,gd->fzs,gd->dzs);
      fprintf(gd->jpfp," nxgd=%d nygd=%d\n\n",gd->nxgd,gd->nygd);
      fprintf(gd->jpfp," nzo=%d fzo=%g dzo=%g\n",gd->nzo,gd->fzo,gd->dzo);
      fprintf(gd->jpfp," nxo=%d fxo=%g dxo=%g\n",gd->nxo,gd->fxo,gd->dxo);
      fprintf(gd->jpfp," nyo=%d fyo=%g dyo=%g\n\n",gd->nyo,gd->fyo,gd->dyo);
      fprintf(gd->jpfp," nt=%d ft=%g dt=%g dxm=%g\n\n",gd->nt,gd->ft,gd->dt,dxm);
      fprintf(gd->jpfp," nxoffset=%d fxoffset=%g dxoffset=%g\n",
	    nxoffset,fxoffset,dxoffset);
      fprintf(gd->jpfp," nyoffset=%d fyoffset=%g dyoffset=%g\n",
	    nyoffset,fyoffset,dyoffset);
      fprintf(gd->jpfp," xaper=%g xoffsetmax=%g\n",
	    gd->xaper,xoffsetmax);
      fprintf(gd->jpfp," yaper=%g yoffsetmax=%g\n",
	    gd->yaper,yoffsetmax);
      fprintf(gd->jpfp," ntrmax=%d pptr=%d\n",ntrmax,pptr);
      fprintf(gd->jpfp," ext=%g exo=%g\n",ext,exo);
      fprintf(gd->jpfp," eyt=%g eyo=%g\n",eyt,eyo);
      fprintf(gd->jpfp," ezs=%g ezo=%g\n",ezs,ezo);
      fprintf(gd->jpfp," ====================================\n");
      fflush(gd->jpfp);

      if (gd->fxgd>gd->fxo || ext<exo ||
	  gd->fygd>gd->fyo || eyt<eyo) {
            fprintf(gd->jpfp,"fxgd=%g > fxo=%g ?\n",gd->fxgd,gd->fxo);
            fprintf(gd->jpfp,"ext=%g < exo=%g ?\n",ext,exo);
            fprintf(gd->jpfp,"fygd=%g > fyo=%g ?\n",gd->fygd,gd->fyo);
            fprintf(gd->jpfp,"eyt=%g < eyo=%g ?\n",eyt,eyo);
	    err(" migration output range is out of tttable");
      }

      /****************************************************
      Allocate space
      ****************************************************/
      fprintf(gd->jpfp,"allocating mig of size %g Mbytes\n",(4.0e-6)*
	    gd->nzo*gd->nxo*gd->nyo*nxoffset*nyoffset);

      mig =ealloc5float(gd->nzo,gd->nxo,gd->nyo,nxoffset,nyoffset);

      fprintf(gd->jpfp,"allocating ttab, ctab and rtab of size %g Mbytes\n",(12.0e-6)*
	    gd->nzs*gd->nxs*gd->nys*gd->nxt*gd->nyt);
      ttab=ealloc6float(gd->nxt,gd->nyt,gd->nzs,gd->nxs,gd->nys,gd->multit);
      ctab=ealloc5float(gd->nxt,gd->nyt,gd->nzs,gd->nxs,gd->nys);
      rtab=ealloc5float(gd->nxt,gd->nyt,gd->nzs,gd->nxs,gd->nys);

      fprintf(gd->jpfp,"allocating tbuf, cbuf and rbuf of size %g Mbytes\n",(4.0e-6)*
	    gd->nys*gd->nxo*gd->nzo*gd->nyt*gd->nxt);
      tbuf=ealloc6ushort(gd->multit,gd->nxt,gd->nyt,gd->nzo,gd->nxo,gd->nys);
      cbuf=ealloc5uchar(gd->nxt,gd->nyt,gd->nzo,gd->nxo,gd->nys);
      rbuf=ealloc5uchar(gd->nxt,gd->nyt,gd->nzo,gd->nxo,gd->nys);

      trf=ealloc1float(gd->nt + 2*gd->ntrimax);

      memset((void *)mig[0][0][0][0],(int)'\0',
	    nxoffset*nyoffset*gd->nxo*gd->nyo*gd->nzo*sizeof(float));

      fprintf(gd->jpfp,"\nReading ttables ...");
      fprintf(gd->jpfp,"\n");

      /*****************************************************		       
      Compute traveltime residual, which is the difference
      of the tttable and the traveltimes calculated
      from the reference model	
      *****************************************************/
      for (iys=0;iys<gd->nys;iys++) {
	    for(ixs=0;ixs<gd->nxs;ixs++){			
		  for (izs=0;izs<gd->nzs;izs++) {
			fprintf(gd->jpfp,
			      "Reading source No. (%d,%d,%d)\n",ixs,iys,izs);
			for (iyt=0;iyt<gd->nyt;iyt++) {
			      if (!fgettr(ttfp,&ttr))
				    err("Can't read ttfile");
			      for (imul=0;imul<gd->multit;imul++) {
			            for (ixt=0;ixt<gd->nxt;ixt++) {
				          ttab[imul][iys][ixs][izs][iyt][ixt]
					        =ttr.data[ixt+imul*gd->nxt];
			            }
			      }
                              if (gd->crfp!=NULL) {
                                    if (fread(ctab[iys][ixs][izs][iyt],
					  sizeof(float),gd->nxt,gd->crfp)!=gd->nxt) 
				          err("Can not read in crfile");
				    for (ixt=0;ixt<gd->nxt;ixt++) 
				          ctab[iys][ixs][izs][iyt][ixt]=1.0;

				    if (fread(rtab[iys][ixs][izs][iyt],
					  sizeof(float),gd->nxt,gd->crfp)!=gd->nxt) 
				          err("Can not read in crfile");
				    for (ixt=0;ixt<gd->nxt;ixt++) 
				          rtab[iys][ixs][izs][iyt][ixt]=1.0;
			      }
			}
		  }
	    }
      }

      fclose(ttfp);
      fprintf(gd->jpfp,"\nTable reading done\n");	

      fprintf(gd->jpfp,"\nGet traveltime buffer\n");
      get_bufs(
            gd, 
            ttab,       /*tttable, of size multit*nys*nxs*nzs*nyt*gd->nxt*/
            ctab,       /*tttable, of size nys*nxs*nzs*nyt*nxt*/
            rtab,       /*tttable, of size nys*nxs*nzs*nyt*nxt*/
            tbuf,       /*shot at nys*nxo*nzo*nyt*nxt*/
	    cbuf,       /*shot at nys*nxo*nzo*nyt*nxt*/
            rbuf);      /*shot at nys*nxo*nzo*nyt*nxt*/

      fprintf(gd->jpfp,"Start migration ...\n");
      fflush(gd->jpfp);
	
      imigtr=1;
      itotmigtr=0; 

      do { 
	    /*****************************************************	
	    To determine offset index
	    *****************************************************/
	    gd->sx=tr.sx;    /*source position in x*/
	    gd->gx=tr.gx;    /*geophone position in x*/
	    gd->sy=tr.sy;    /*source position in y*/
	    gd->gy=tr.gy;    /*geophone position in y*/ 

            #ifdef DEBUG
	    gd->sy=3;	
	    gd->gy=3;
	    fprintf(gd->jpfp,"trace %d: sx,gx,sy,gy=%g %g %g %g\n",
		  imigtr,gd->sx,gd->gx,gd->sy,gd->gy);
	    fprintf(gd->jpfp,"trace in line: %d; trace in reel: %d\n",
		  tr.tracl,tr.tracr);
            #endif

            xm=0.5*(gd->sx+gd->gx);
            ym=0.5*(gd->sy+gd->gy);

	    if (xm<gd->fxs || xm>exs || 
	        MAX(gd->gx-gd->sx,gd->sx-gd->gx)>xoffsetmax  || 
		ym<gd->fys || ym>eys || 
		MAX(gd->gy-gd->sy,gd->sy-gd->gy)>yoffsetmax )
		  continue; 

	    /****************************************************
	    Determine the index for offset in x, fxoffset is the 
	    starting offset
	    ****************************************************/
	    ixoffset=(int)((gd->gx-gd->sx-fxoffset)/dxoffset+0.5);
	    if (ixoffset<0) ixoffset=0;	
	    if (ixoffset>=nxoffset) ixoffset=nxoffset-1;

	    /****************************************************
	    Determine the index for offset in y, fxoffset is the 
	    starting offset
	    ****************************************************/
	    iyoffset=(int)((gd->gy-gd->sy-fyoffset)/dyoffset+0.5);
	    if (iyoffset<0) iyoffset=0;	
	    if (iyoffset>=nyoffset) iyoffset=nyoffset-1;

	    /*******************************************
	    Take the t-derivative of the data
	    *******************************************/
	    /*For the 1D synthetic data, turn this off*/

            filt(tr.data,gd,trf); 

	    mig3d(gd,
	       	  tr.data,  /*the data of this trace*/
                  trf,      /*integrated trace*/
		  mig[iyoffset][ixoffset],/*of nyo*nxo*nzo*/
		  tbuf,	    
		  cbuf,
		  rbuf); 

	    itotmigtr++;
	    if ((imigtr-1)%pptr ==0 ){
		  fprintf(gd->jpfp,"Migrated trace %d\n",imigtr);
		  fflush(gd->jpfp);
    	    }
	    imigtr++;
      } while (fgettr(infp,&tr) && imigtr<ntrmax);

      fprintf(gd->jpfp,"Migrated %d traces in total\n",itotmigtr);

      memset((void *)&tro,(int)'\0',sizeof(segy));
      tro.ns=gd->nzo;
      tro.dt=1000*(int)(1000*gd->dt+0.5);
      tro.delrt=0.0;
      tro.f1=gd->fzo;
      tro.d1=gd->dzo;
      tro.f2=gd->fxo;
      tro.d2=gd->dxo;
      for (iyo=0;iyo<gd->nyo;iyo++) {
	    for (ixo=0;ixo<gd->nxo;ixo++) {
		  for(iyoffset=0;iyoffset<nyoffset;iyoffset++) {
			for(ixoffset=0;ixoffset<nxoffset;ixoffset++)  {
			      memcpy((void *)tro.data,(const void*)
			      mig[iyoffset][ixoffset][iyo][ixo],
			      gd->nzo*sizeof(float));
			      tro.offset=fxoffset+ixoffset*dxoffset;