sukdmig2d.c

su 的源代码库

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

/* SUKDMIG2D: $Revision: 1.24 $ ; $Date: 2006/11/07 22:58:42 $	*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
"SUKDMIG2D - Kirchhoff Depth Migration of 2D poststack/prestack data	",
" 									",
"    sukdmig2d  infile=  outfile=  ttfile=   [parameters] 		",
"									",
" Required parameters:							",
" infile=stdin		file for input seismic traces			",
" outfile=stdout	file for common offset migration output  	",
" ttfile=		file for input traveltime tables		",
"									",
" ...  The following 9 parameters describe traveltime tables:		",
" fzt= 			first depth sample in traveltime table		",
" nzt= 			number of depth samples in traveltime table	",
" dzt=			depth interval in traveltime table		",
" fxt=			first lateral sample in traveltime table	",
" nxt=			number of lateral samples in traveltime table	",
" dxt=			lateral interval in traveltime table		",
" fs= 			x-coordinate of first source			",
" ns= 			number of sources				",
" ds= 			x-coordinate increment of sources		",
"									",
" Optional Parameters:							",
" 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) 	sampling interval of midpoints 		",
" fzo=fzt		    z-coordinate of first point in output trace	",
" dzo=0.2*dzt		vertical spacing of output trace 		",
" nzo=5*(nzt-1)+1 	number of points in output trace		",	
" fxo=fxt		    x-coordinate of first output trace 		",
" dxo=0.5*dxt		horizontal spacing of output trace 		",
" nxo=2*(nxt-1)+1  	number of output traces 			",	
" off0=0		   	first offest in output 			",
" doff=99999		offset increment in output 			",
" noff=1	 	number of offsets in output 			",	
" absoff=0      flag for using absolute offsets of input traces		",
"               =0 means use offset=gx-sx                		",
"               =1 means use abs(gx-sx)                  		",
" limoff=0      flag for only using input traces that fall within the range",
"               of defined output offset bins (off0,doff,noff) 		",
"               =0 means use all input traces                 		",
"               =1 means limit traces used by offset           		",
" fmax=0.25/dt		frequency-highcut for input traces		",
" offmax=99999		maximum absolute offset allowed in migration 	",
" aperx=nxt*dxt/2  	migration lateral aperature 			",
" angmax=60		migration angle aperature from vertical 	",
" v0=1500(m/s)		reference velocity value at surface		",	
" dvz=0.0  		reference velocity vertical gradient		",
"									",
" ls=1			flag for line source				",
" jpfile=stderr		job print file name 				",
"									",
" mtr=100  		print verbal information at every mtr traces	",
" ntr=100000		maximum number of input traces to be migrated	",
" npv=0			flag of computing quantities for velocity analysis",
" rscale=1000.0 	scaling for roundoff error suppression		",
"									",
"   ...if npv>0 specify the following three files:			",
" tvfile=tvfile		input file of traveltime variation tables	",
"			tv[ns][nxt][nzt]				",
" csfile=csfile		input file of cosine tables cs[ns][nxt][nzt]	",
" outfile1=outfile1	file containning additional migration output   	",
"			with extra amplitude				",
"									",
" Notes:								",
" 1. Traveltime tables were generated by program rayt2d (or other ones)	",
"    on relatively coarse grids, with dimension ns*nxt*nzt. In the	",
"    migration process, traveltimes are interpolated into shot/gephone 	",
"    positions and output grids.					",
" 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 nxo*noff*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			",
"    	[ns*nxt*nzt+noff*nxo*nzo+4*nr*nzt+5*nxt*nzt+npa*(2*ns*nxt*nzt   ",
"	+noff*nxo*nzo+4*nxt*nzt)]*4 bytes				",
"    where nr = 1+min(nxt*dxt,0.5*offmax+aperx)/dxo. 			",
" 6. Amplitudes are computed using the reference velocity profile, v(z),",
"    specified by the parameters v0= and dvz=.				",
" 7. Input traces must specify source and receiver positions via the header",
"    fields tr.sx and tr.gx. Offset is computed automatically.		",
" 8. if limoff=0, input traces from outside the range defined by off0, doff, ",
"    noff, will get migrated into the extremal offset bins/planes.  E.g. if ",
"    absoff=0 and limoff=0, all traces with gx<sx will get migrated into the ",
"    off0 bin.",
"									",
NULL};
/*
 * Author:  Zhenyue Liu, 03/01/95,  Colorado School of Mines 
 * Modifcations:
 *    Gary Billings, Talisman Energy, Sept 2005:  added absoff, limoff.
 *
 * Trace header fields accessed: ns, dt, delrt, d2
 * Trace header fields modified: sx, gx
 */
 
/**************** end self doc ***********************************/
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,
	float **tb,float **t,float x0);
  void interpx(int nxt,float fxt,float dxt,int nx,float fx,float dx,
	int nzt,float **tt,float **t);
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t);
  void timeb(int nr,int nz,float dr,float dz,float fz,float a,
	float v0,float **t,float **p,float **sig,float **ang);
  void mig2d(float *trace,int nt,float ft,float dt,
	float sx,float gx,float **mig,float aperx,
  	int nx,float fx,float dx,float nz,float fz,float dz,
	int ls,int mtmax,float dxm,float fmax,float angmax,
	float **tb,float **pb,float **cs0b,float **angb,int nr,
	float **tsum,int nzt,float fzt,float dzt,int nxt,float fxt,float dxt,
	int npv,float **cssum,float **tvsum,float **mig1);

/* define */
#define RSCALE_KDMIG 1000.0

/* segy trace */
segy tr, tro;

int
main (int argc, char **argv)
{
	int nt;		/* number of time samples in input data		*/
	int nzt;	/* number of z-values in traveltime table	*/
	int nxt;	/* number of x-values in traveltime table	*/
	int nzo;	/* number of z-values in output data		*/
	int nxo;	/* number of x-values in output data		*/
	int ns;		/* number of sources				*/
	int noff;	/* number of offsets in output data		*/
	int nr;
	int is,io,ixo,izo; /* counters */
	int ls,ntr,jtr,ktr,mtr,npv,mtmax;
	int   offset,absoff,limoff;
	off_t nseek;
	float   ft,fzt,fxt,fzo,fxo,fs,off0,dt,dzt,dxt,dzo,dxo,ds,doff,dxm,
		ext,ezt,ezo,exo,es,s,scal;	
	float v0,dvz,fmax,angmax,offmax,rmax,aperx,sx=0,gx=0;
	float ***mig=NULL,***ttab=NULL,**tb=NULL,**pb=NULL;
	float **cs0b=NULL,**angb=NULL,**tsum=NULL,**tt=NULL;
	float **tvsum=NULL,***mig1=NULL,***cs=NULL,***tv=NULL,**cssum=NULL;

	float rscale;			/* scaling factor for roundoff */
	
	char *infile="stdin",*outfile="stdout",*ttfile,*jpfile,*tvfile,
	     *csfile,*outfile1;
	FILE *infp,*outfp,*ttfp,*jpfp,*tvfp=NULL,*out1fp=NULL,*csfp=NULL;


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

	/* open input and output files	*/
	if( !getparstring("infile",&infile)) {
		infp = stdin;
	} else  
		if ((infp=fopen(infile,"r"))==NULL)
			err("cannot open infile=%s\n",infile);
	if( !getparstring("outfile",&outfile)) {
		outfp = stdout;
	} else  
		outfp = fopen(outfile,"w");
	efseeko(infp,(off_t) 0,SEEK_CUR);
	efseeko(outfp,(off_t) 0,SEEK_CUR);
	if( !getparstring("ttfile",&ttfile))
		err("must specify ttfile!\n");
	if ((ttfp=fopen(ttfile,"r"))==NULL)
		err("cannot open ttfile=%s\n",ttfile);
	if( !getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	} else  
		jpfp = fopen(jpfile,"w");

	/* get information from the first header */
	if (!fgettr(infp,&tr)) err("can't get first trace");
	nt = tr.ns;
	if (!getparfloat("dt",&dt)) dt = ((double) tr.dt)/1000000.0; 
	if (dt<0.0000001) err("dt must be positive!\n");
	if (!getparfloat("ft",&ft)) ft = tr.delrt/1000.0; 
 	if (!getparfloat("dxm",&dxm)) dxm = tr.d2;
	if  (dxm<0.0000001) err("dxm must be positive!\n");
	
	/* get traveltime tabel parameters	*/
	if (!getparint("nxt",&nxt)) err("must specify nxt!\n");
	if (!getparfloat("fxt",&fxt)) err("must specify fxt!\n");
	if (!getparfloat("dxt",&dxt)) err("must specify dxt!\n");
	if (!getparint("nzt",&nzt)) err("must specify nzt!\n");
	if (!getparfloat("fzt",&fzt)) err("must specify fzt!\n");
	if (!getparfloat("dzt",&dzt)) err("must specify dzt!\n");
	if (!getparint("ns",&ns)) err("must specify ns!\n");
	if (!getparfloat("fs",&fs)) err("must specify fs!\n");
	if (!getparfloat("ds",&ds)) err("must specify ds!\n");
	if (!getparfloat("rscale",&rscale))   rscale = RSCALE_KDMIG;

	ext = NINT(rscale*(fxt+(nxt-1)*dxt))/rscale;
	ezt = NINT(rscale*(fzt+(nzt-1)*dzt))/rscale;
	es = NINT(rscale*(fs+(ns-1)*ds))/rscale;

	/* optional parameters	*/
	if (!getparint("nxo",&nxo)) nxo = (nxt-1)*2+1;
	if (!getparfloat("fxo",&fxo)) fxo = fxt;
	if (!getparfloat("dxo",&dxo)) dxo = dxt*0.5;
	if (!getparint("nzo",&nzo)) nzo = (nzt-1)*5+1;
	if (!getparfloat("fzo",&fzo)) fzo = fzt;
	if (!getparfloat("dzo",&dzo)) dzo = dzt*0.2;
	exo = NINT(rscale*(fxo+(nxo-1)*dxo))/rscale;
	ezo = NINT(rscale*(fzo+(nzo-1)*dzo))/rscale;

	fprintf(jpfp," fxt=%f fxo=%f \n",fxt,fxo);
	fprintf(jpfp," ext=%f exo=%f \n",ext,exo);
	fprintf(jpfp," fzt=%f fzo=%f \n",fzt,fzo);
	fprintf(jpfp," ezt=%f ezo=%f \n",ezt,ezo);
	fprintf(jpfp," es=%f \n",es);
	fprintf(jpfp," \n");
	if(fxt>fxo || ext<exo || fzt>fzo || ezt<ezo) {
		warn("This condition must NOT be satisfied: fxt>fxo || ext<exo || fzt>fzo || ezt<ezo");
		warn("fxt=%f fxo=%f ext=%f exo=%f fzt=%f fzo=%f ezt=%f ezo=%f",
			  fxt,fxo,ext,exo,fzt,fzo,ezt,ezo);
		err(" migration output range is out of traveltime table!\n");
	  }

	if (!getparfloat("v0",&v0)) v0 = 1500;
	if (!getparfloat("dvz",&dvz)) dvz = 0;
	if (!getparfloat("angmax",&angmax)) angmax = 60.;
	if  (angmax<0.00001) err("angmax must be positive!\n");
	mtmax = 2*dxm*sin(angmax*PI/180.)/(v0*dt);
	if(mtmax<1) mtmax = 1;
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;
	if (!getparfloat("offmax",&offmax)) offmax = 3000;
	if (!getparfloat("fmax",&fmax)) fmax = 0.25/dt;
	if (!getparint("noff",&noff))	noff = 1;
	if (!getparfloat("off0",&off0)) off0 = 0.;
	if (!getparfloat("doff",&doff)) doff = 99999;

	if (!getparint("ls",&ls))	ls = 1;
	if (!getparint("absoff",&absoff))   absoff=0;
	if (!getparint("limoff",&limoff))   limoff=0;
	if (!getparint("ntr",&ntr))	ntr = 100000;
	if (!getparint("mtr",&mtr))	mtr = 100;
	if (!getparint("npv",&npv))	npv = 0;
	if(npv){
		if( !getparstring("tvfile",&tvfile))
			err("must specify tvfile!\n");
 		tvfp = fopen(tvfile,"r");
		if( !getparstring("csfile",&csfile))
			err("must specify csfile!\n");
 		csfp = fopen(csfile,"r");
		if( !getparstring("outfile1",&outfile1))
			outfile1="outfile1";
 		out1fp = fopen(outfile1,"w");
	}

	fprintf(jpfp,"\n");
	fprintf(jpfp," Migration parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," infile=%s \n",infile);
	fprintf(jpfp," outfile=%s \n",outfile);
	fprintf(jpfp," ttfile=%s \n",ttfile);
	fprintf(jpfp," \n");
	fprintf(jpfp," nzt=%d fzt=%g dzt=%g\n",nzt,fzt,dzt);
	fprintf(jpfp," nxt=%d fxt=%g dxt=%g\n",nxt,fxt,dxt);
 	fprintf(jpfp," ns=%d fs=%g ds=%g\n",ns,fs,ds);
	fprintf(jpfp," \n");
	fprintf(jpfp," nzo=%d fzo=%g dzo=%g\n",nzo,fzo,dzo);
	fprintf(jpfp," nxo=%d fxo=%g dxo=%g\n",nxo,fxo,dxo);
	fprintf(jpfp," \n");
	
	/* compute reference traveltime and slowness  */
	rmax = MAX(es-fxt,ext-fs);
	rmax = MIN(rmax,0.5*offmax+aperx);
	nr = 2+(int)(rmax/dxo);
	tb = ealloc2float(nzt,nr);
	pb = ealloc2float(nzt,nr);
	cs0b = ealloc2float(nzt,nr);
	angb = ealloc2float(nzt,nr);
	timeb(nr,nzt,dxo,dzt,fzt,dvz,v0,tb,pb,cs0b,angb);

	fprintf(jpfp," nt=%d ft=%g dt=%g \n",nt,ft,dt);
 	fprintf(jpfp," dxm=%g fmax=%g\n",dxm,fmax);
 	fprintf(jpfp," noff=%d off0=%g doff=%g\n",noff,off0,doff);
	fprintf(jpfp," v0=%g dvz=%g \n",v0,dvz);
 	fprintf(jpfp," aperx=%g offmax=%g angmax=%g\n",aperx,offmax,angmax);
 	fprintf(jpfp," ntr=%d mtr=%d ls=%d npv=%d\n",ntr,mtr,ls,npv);
	fprintf(jpfp," absoff=%d limoff=%d\n",absoff,limoff);
	if(npv)
 	  fprintf(jpfp," tvfile=%s csfile=%s outfile1=%s\n",
		tvfile,csfile,outfile1);
	fprintf(jpfp," ================\n");
	fflush(jpfp);

	/* allocate space */
	mig = ealloc3float(nzo,nxo,noff);
	ttab = ealloc3float(nzt,nxt,ns);
	tt = ealloc2float(nzt,nxt);
	tsum = ealloc2float(nzt,nxt);
	if(npv){
		tv = ealloc3float(nzt,nxt,ns);
		tvsum = ealloc2float(nzt,nxt);
		cs = ealloc3float(nzt,nxt,ns);
		cssum = ealloc2float(nzt,nxt);
	}
	if(!npv) 
		mig1 = ealloc3float(1,1,noff);
	else
		mig1 = ealloc3float(nzo,nxo,noff);


 	memset((void *) mig[0][0],0,noff*nxo*nzo*sizeof(float)); 
	 if(npv)
		memset((void *) mig1[0][0], 0,
			noff*nxo*nzo*sizeof(float));

 	fprintf(jpfp," input traveltime tables \n");
			     
	/* compute traveltime residual	*/
	for(is=0; is<ns; ++is){
		nseek = (off_t) nxt*nzt*is;
		efseeko(ttfp,nseek*((off_t) sizeof(float)),SEEK_SET);
		fread(ttab[is][0],sizeof(float),nxt*nzt,ttfp);
		s = fs+is*ds;
		resit(nxt,fxt,dxt,nzt,nr,dxo,tb,ttab[is],s);
		if(npv) {
			efseeko(tvfp, nseek*((off_t) sizeof(float)),SEEK_SET);
			fread(tv[is][0],sizeof(float),nxt*nzt,tvfp);
			efseeko(csfp, nseek*((off_t) sizeof(float)),SEEK_SET);
			fread(cs[is][0],sizeof(float),nxt*nzt,csfp);
		}
		
	}

	fprintf(jpfp," start migration ... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
	

	jtr = 1;
	ktr = 0;

	  fprintf(jpfp," fs=%g es=%g offmax=%g\n",fs,es,offmax);

	do {
		/* determine offset index	*/
	    float as,res;

		if (tr.scalco) { /* if tr.scalco is set, apply value */
			if (tr.scalco>0) {
				sx = tr.sx*tr.scalco;
				gx = tr.gx*tr.scalco;
		   	} else { /* if tr.scalco is negative divide */
				sx = tr.sx/ABS(tr.scalco);
				gx = tr.gx/ABS(tr.scalco);
			}
			
		} else {
			     sx = tr.sx;
			     gx = tr.gx;
		}

		/* GWB 2005.09.22: */
		/* io = (int)((gx-sx-off0)/doff+0.5); */
		offset=gx-sx;
		if( absoff && offset<0 )offset=-offset;
		io = (int)((offset-off0)/doff+0.5);
		if( limoff && (io<0 || io>=noff) ) continue;
		/* end of GWB 2005.09.22 */

	    if(io<0) io = 0;
	    if(io>=noff) io = noff-1;

		/* fprintf(jpfp," read trace jtr=%d: sx=%g gx=%g io=%d ktr=%d\n",jtr,sx,gx,io,ktr); */

	    if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && 
		 MAX(gx-sx,sx-gx)<=offmax ){
		/*     migrate this trace	*/

		    /* fprintf(jpfp," Good! Condition NOT satisfied\n"); */

	    	as = (sx-fs)/ds;
	    	is = (int)as;
		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tsum);
		if(npv)  {
			sum2(nxt,nzt,1-res,res,tv[is],tv[is+1],tvsum);
			sum2(nxt,nzt,1-res,res,cs[is],cs[is+1],cssum);
		}
		
	    	as = (gx-fs)/ds;
	    	is = (int)as;
		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tt);
		sum2(nxt,nzt,1,1,tt,tsum,tsum);
		if(npv)  {
			sum2(nxt,nzt,1-res,res,tv[is],tv[is+1],tt);
			sum2(nxt,nzt,1,1,tt,tvsum,tvsum);
			sum2(nxt,nzt,1-res,res,cs[is],cs[is+1],tt);
			sum2(nxt,nzt,1,1,tt,cssum,cssum);
		}

		mig2d(tr.data,nt,ft,dt,sx,gx,mig[io],aperx,
		  nxo,fxo,dxo,nzo,fzo,dzo,
		  ls,mtmax,dxm,fmax,angmax,
		  tb,pb,cs0b,angb,nr,tsum,nzt,fzt,dzt,nxt,fxt,dxt,
		  npv,cssum,tvsum,mig1[io]);

		  ktr++;
		  if((jtr-1)%mtr ==0 ){
			fprintf(jpfp," migrated trace %d\n",jtr);
			fflush(jpfp);
	    	}
	    }
	    jtr++;
	} while (fgettr(infp,&tr) && jtr<ntr);

	fprintf(jpfp," migrated %d traces in total\n",ktr);

	memset((void *) &tro, 0, sizeof(segy));
	tro.ns = nzo;
	tro.d1 = dzo;
	tro.dt = 1000*(int)(1000*dt+0.5);
	tro.delrt = 0.0;
	tro.f1 = fzo;
	tro.f2 = fxo;
	tro.d2 = dxo;
	tro.trid = 200;

	scal = 4/sqrt(PI)*dxm/v0;
	for(ixo=0; ixo<nxo; ixo++) {
		for(io=0; io<noff; io++)  {
			memcpy((void *) tro.data,
				 (const void *) mig[io][ixo], nzo*sizeof(float));
			tro.offset = off0+io*doff;
			tro.tracr = 1+ixo;
			tro.tracl = 1+io+noff*ixo;
			tro.cdp = fxo+ixo*dxo;
			tro.cdpt = 1+io;

			for(izo=0; izo<nzo; ++izo)
				tro.data[izo] *=scal;

			/* write out */
			fputtr(outfp,&tro);