kdmig3d.c

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

/* Copyright (c) Colorado School of Mines, 2005.*/
/* All rights reserved.                       */

/* SUKDMIG3D: $Revision: 1.4 $ ; $Date: 2003/06/09 16:17:07 $	*/

#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[]={
" 									",
"SUKDMIG3D - Kirchhoff Depth Migration of 3D poststack/prestack data	",
" 									",
" 	sukdmig3d datain= dataout= [parameters] 			",
"									",
" Required parameters:							",
" ttfile	 file for input tttables			       	",
"									",
" Optional Parameters:							",
"									",
" datain=stdin	 file for input seismic traces				",
" dataout=stdout file for common offset migration output		",
" crfile=NULL    file for cos theta and ray paths                       ",
"									",
"   The following 17 parameters describe tttables: (from ttfile header)	",
" fxgd= or from header (f1)		first x-sample in tttable	",
" nxt= or from header (ns)		number of x-samples in tttable	",
" dxgd= or from header (d1)		x-interval in tttable 		",
" fygd= or from header (f2)		first y-sample in tttable 	",
" nyt= or from header (ntr)		number of y-samples in tttable	",
" dygd= or from header (d2)		y-interval in tttable		",
" ixsf= or from header (sdel) 	        x in dxgd of first source	",
" nxs= or from header (nhs) 		number of sources in x		",
" ixsr= or from header (swdep) 	        ratio of source & gd spacing    ",
" iysf= or from header (gdel)	        y in dygd of first source 	",
" nys= or from header (nvs)		number of sources in y          ",
" iysr= or from header (gwdep)	        ratio of source & gd spacing    ",
" fzs= or from header (sdepth/1000)	first depth sample in tttable	",
" nzs= or from header (duse) 		number of depth samples in tttable",
" dzs= or from header (ep/1000)		depth interval in tttable	",
" nxgd= or from header (selev)		x size of the traveltime region ",
" nygd= or from header (gelev)		y size of the traveltime region ",
" multit= or from header (scalel)       number of multivalued traveltime",
"									",
"   The following two parameters are from data header			",
" dt= or from header (dt)	time sampling interval of input data	",
" ft= or from header (ft)	first time sample of input data		",
" dxm= or from header (d2)      mid point spacing of input data         ",
"									",
" Default: output is 5 times finer in depth and 2 times finer in lateral",
" fzo=fzs	z-coordinate of first point in output trace		",
" dzo=0.2*dzs	vertical spacing of output trace (5 times finer)	",
" nzo=5*(nzs-1)+1 number of points in output trace (5 times finer)	",	
" fxo=fxgd	x-coordinate of first output trace	 		",
" dxo=0.5*dxgd	horizontal spacing of output trace (2 times finer)	",
" nxo=2*(nxgd-1)+1 number of output traces (2 times finer)		",
" fyo=fygd	y-coordinate of first output trace			",
" dyo=0.5*dygd	horizontal spacing of output trace (2 times finer)	",
" nyo=2*(nygd-1)+1 number of output traces (2 times finer)		",
"									",
" Default: poststack migration						",	
" fxoffset=0		first offest in output in x			",
" fyoffset=0		first offest in output in y			",
" dxoffset=99999	offset increment in output in x	 		",
" dyoffset=99999	offset increment in output in y			",
" nxoffset=1		number of offsets in output in x		",
" nyoffset=1		number of offsets in output in y		",	
" xoffsetmax=99999	x-maximum absolute offset allowed in migration 	",
" yoffsetmax=99999	y-maximum absolute offset allowed in migration  ",
" xaper=nxt*dxgd/2.5	migration lateral aperature in x		",
" yaper=nyt*dygd/2.5	migration lateral aperature in y		",
" angmax=60             max angle to handle                             ",
" fmax=0.25/dt          max frequency in the data                       ",
" jpfile=stderr		job print file name 				",
" pptr=100		print verbal information at every pptr traces	",
" ntrmax=100000		maximum number of input traces to be migrated	",
" ls=0                  point =0 line source =1                         ",
"									",
" Notes:								",
" 1. Traveltime tables were generated by program SUTETRARAY (or other	",
"    ones) on very sparse tetrahedral model, with dimension nys*nxs*nzs ",
"    *nyt*nxt.                                                          ",
" 2. Input seismic traces must be SU format and can be any type of 	",
"    gathers (common shot, common offset, common CDP, and so on).	", 
" 3. Migrated traces are output in CDP gathers if velocity analysis	",
"    is required, with dimension nyoffset*nxoffset*nyo*nxo*nzo.		", 
" 4. If the offset value of an input trace is not in the offset array 	",
"    of output, the nearest one in the array is chosen. 		",
" 5. Memory requirement for this program is about			",
"    [nys*nxs*nzs*nyt*nxt+nyoffset*nxoffset*nxo*nyo*nzo+	       	",
"    nys*nxo*nzo*nyt*nxt]                                               ",
" 6. Input traces must specify source and receiver positions via header	",
"    fields tr.sx and tr.gx, as well as tr.sy and tr.gy. Offset is 	",
"    computed automatically.						",
"									",
" Disclaimer:								",
" This is a research code that will take considerable work to get into	",
" the form of a a production level 3D migration code. The code is	",
" offered as is, along with tetramod and sutetraray, to provide a starting",
" point for researchers who wish to write their own 3D migration codes.",
"									",
NULL};

/*
 * Author:  Zhaobo Meng, 01/10/97,  Colorado School of Mines 
 *
 * Trace header fields accessed: ns, dt, delrt, d2
 * Trace header fields modified: sx, gx
 */
/******************************* end self doc *******************************/

#define DEBUG

#define InfDistance 10
#define InfTime 7000 /*ms*/

#define InfByte1 255
#define InfByte2 65535

struct GD {
      float fxgd;
      float dxgd;
      int nxgd;	        /*number of samples in x in tttable*/
      float fygd;
      float dygd;
      int nygd;	        /*number of samples in tttable in y*/

      int nxt;	        /*number of lateral samples (in x) in tttable*/
      int nyt;	        /*number of lateral samples (in y) in tttable*/

      int nxo;	        /*number of samples in x in output*/
      float fxo;	/*first sample in x in output*/
      float dxo;	/*sample spacing in x in output*/
      int nyo;	        /*number of samples in y in output*/
      float fyo;	/*first sample in y in output*/
      float dyo;	/*sample spacing in y in output*/
      int nzo;	        /*number of samples in z in output*/
      float fzo;	/*first sample in z in output*/
      float dzo;	/*sample spacing in z in output*/

      float fxs;        /*first shot in x*/
      float dxs;        /*shot spacing in x*/
      int nxs;          /*number of sources in x*/
      float fys;        /*first shot in y*/
      float dys;        /*shot spacing in y*/
      int nys;          /*number of sources in y*/
      float fzs;        /*first shot in z*/
      float dzs;        /*shot spacing in z*/
      int nzs;          /*number of sources in z*/
      float xaper;	/*aperture in x*/
      float yaper;	/*aperture in y*/
      int multit;
      int ixsf;         /*first source in x in gd*/
      int ixsr;         /*source spacing/grid spacing*/
      int iysf;         /*first source in y in gd*/
      int iysr;         /*source spacing/grid spacing*/
      float sy;	        /*source position in y*/
      float gy;	        /*geophone position in y*/
      float sx;	        /*source position in x*/
      float gx;	        /*geophone position in y*/
      int nt;		/*number of samples of this trace*/
      float ft;	        /*first sample of this trace*/
      float dt;	        /*sampling rate in t*/
      float v0;
      float fmax;
      int ls;
      int ntrimax;
      FILE *crfp;
      FILE *jpfp;
};

struct TTAB {
      unsigned short t;  /*real t=t/InfByte2*InfTime*/
      unsigned char cs;  /*real cs=cs/InfByte1*/
      unsigned char r;   /*real r=r/InfByte1*InfDistance*/
};

void filt(float *trace, struct GD *gd, float *trf);

void get_bufs( struct GD *gd, float ******ttab, float *****ctab,
      float *****rtab, unsigned short ******buf, unsigned char *****cbuf,
      unsigned char  *****rbuf);

void mig3d( struct GD *gd, float *trace, float *trf, float ***mig,
      unsigned short ******buf, unsigned char  *****cbuf,
      unsigned char  *****rbuf);

void preproc( float *data, struct GD *gd);

/*segy trace*/
segy tr,tro,ttr;

int
main (int argc, char **argv)
{
      struct GD gd0;
      struct GD *gd=&gd0;

      /*************************************************************
	Memory needed for the computed area (for nxoff=nyoff=1):
		nxo * nyo * nzo * 4; (=100*100*100*4=8Mbytes)	

	Memory needed for the ttab:
		nys * nxo *nzo * nxt * nyt * 2; 
		(=20*100*200*40*40*2=1280Mbytes)

		(0,nygd)
			 ______________________________________ 
			|       |	    |		       |
			|       |	    |		       |
			|       |  ttable   |		       |
			|       |  for 1 s  |	     	       |
			|       |	    |		       |
			|_______|___________|__________________|
		(0,0) (ixsf+ixsr*ixs-nxt/2,0) (ixsf+ixsr*ixs+nxt/2,0) (nxgd,0)	

      **************************************************************/
      int nxoffset;     /*number of offsets in output in x*/
      int nyoffset;     /*number of offsets in output in y*/
      int ixoffset;     /*index for nxoffset*/
      int iyoffset;     /*index for nyoffset*/
      int ixo;	        /*index for nxo*/
      int iyo;	        /*index for nyo*/	
      int ixt;	        /*index for nxt*/
      int iyt;	        /*index for nyt*/
      int izs;	        /*index for nzs*/
      int ntrmax;	/*maximum number of input traces to be migrated*/
      int imigtr;	/*index for migrated traces*/
      int itotmigtr;    /*index for migrated traces*/
      int pptr;	        /*print verbal information for every pptr traces*/

      float fxoffset;   /*first offset in output in x*/
      float fyoffset;   /*first offset in output in y*/
      float dxm;        /*mid point spacing*/

      float angmax;     /*max angle to be able to handle*/
      float dxoffset;   /*offset increment in output in x*/
      float dyoffset;   /*offset increment in output in y*/
      float ext;	/*ending t in x=fxgd+(nxt-1)*dxgd*/
      float eyt;	/*ending t in y=fygd+(nyt-1)*dygd*/
      float ezs;	/*last time in z=fzs+(nzs-1)*dzs*/
      float ezo;	/*last point in output ezo=fzo+(nzo-1)*dzo*/
      float exo;	/*last point in output ex0=fxo+(nxo-1)*dxo*/
      float eyo;	/*last point in output ey0=fyo+(nyo-1)*dyo*/
      float exs;	/*last source exs=fxs+(nxs-1)*dxs*/
      float eys;	/*last source eys=fys+(nys-1)*dys*/
      float xoffsetmax; /*maximum absolute x offset allowed in migration*/
      float yoffsetmax; /*maximum absolute y offset allowed in migration*/
      float *****mig;   /*migrated data of size nyoffset*nxoffset*nyo*nxo*nzo*/
      float ******ttab; /*travetime table of size multit*nys*nxs*nzs*nyt*nxt*/
      float *****ctab;  /*travetime table of size nys*nxs*nzs*nyt*nxt*/
      float *****rtab;  /*travetime table of size nys*nxs*nzs*nyt*nxt*/

      unsigned short ******tbuf; /*traveltime buffer of size nys*nxo*nzo*nyt*nxt*/
      unsigned char  *****cbuf;  /*cos theta buffer of size nys*nxo*nzo*nyt*nxt*/
      unsigned char  *****rbuf;  /*cos theta buffer of size nys*nxo*nzo*nyt*nxt*/

      int iys,ixs;      /*indices for nys, nxs*/

      float *trf;       /*integrated traced*/
      int imul;

      char *datain="stdin";/*input data*/
      char *dataout="stdout";/*output data*/
      char *ttfile;     /*tttable to be read in*/
      char *crfile;     /*cos and r to be read in*/
      char *jpfile;     /*job print file name*/
      FILE *infp;	/*input file pointer*/
      FILE *outfp;	/*output file pointer*/
      FILE *ttfp;	/*tttable file pointer*/

      float xm,ym;      /*mid point*/

      /*hook up getpar to handle the parameters*/
      initargs(argc,argv);
      requestdoc(1);

      /*open input and output files*/
      if (!getparstring("datain",&datain)) {
	    infp=stdin;
      } else
	    if ((infp=fopen(datain,"r"))==NULL)
		  err("cannot open datain=%s",datain);

      if ( !getparstring("dataout",&dataout)) {
	    outfp=stdout;
      } else
	    outfp=fopen(dataout,"w");

      /***********************************************************
      Repositon a file on a logical unit: e.g. infp;
      efseeko(infp,offset,from), offset is an offset in bytes
      relative to the position specified by offset, offset=0 
      means beginning;=1 means current position;=2 means the end
      return 0 when successful, otherwise a system error code.
      ***********************************************************/
      efseeko(infp,(off_t) 0,SEEK_CUR);
      efseeko(outfp,(off_t) 0,SEEK_CUR);

      /*****************************************************
      Open tttable file, get parameters from the table
      *****************************************************/
      if ( !getparstring("ttfile",&ttfile))
	    err("must specify ttfile!");
      if ((ttfp=fopen(ttfile,"r"))==NULL)
	    err("Can not open ttfile=%s",ttfile);
      if (!fgettr(ttfp,&ttr))
	    err("Can't get first trace from tttable");

      if (getparstring("crfile",&crfile)) {
            if ((gd->crfp=fopen(crfile,"r"))==NULL)
	          err("Can not open crfile=%s",crfile);
      }
	else
	gd->crfp=NULL;
	

      if (!getparfloat("fxgd",&gd->fxgd)) gd->fxgd=ttr.f1;
      if (!getparfloat("dxgd",&gd->dxgd)) gd->dxgd=ttr.d1;
      if (!getparint("nxt",&gd->nxt)) 	  gd->nxt=ttr.ns;

      if (gd->dxgd<0.1) /* km used I guess */
            gd->v0=2.5;
      else
            gd->v0=2500.0;

      if (!getparint("multit",&gd->multit)) gd->multit=ttr.scalel;
      if (gd->multit==0) gd->multit=1;
      gd->nxt/=gd->multit;

      if (!getparint("ls",&gd->ls)) gd->ls=0;

      if (!getparfloat("fygd",&gd->fygd)) gd->fygd=ttr.f2;
      if (!getparfloat("dygd",&gd->dygd)) gd->dygd=ttr.d2;