kdmig2d.c

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

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

/* SUKDMIG2D: $Revision: 1.18 $ ; $Date: 2004/12/09 17:47:22 $	*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
"SUKDMIG2D - Kirchhoff Depth Migration of 2D poststack/prestack data	",
" 									",
"    sukdmig2d  infile=  outfile=  [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 			",	
" 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",
"   ...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]	",
" dataout1=dataout1	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.		",
"									",
NULL};
/*
 * Author:  Zhenyue Liu, 03/01/95,  Colorado School of Mines 
 *
 * 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);

/* segy trace */
segy tr, tro;

int
main (int argc, char **argv)
{
	int 	nt,nzt,nxt,nzo,nxo,ns,noff,nr,is,io,ixo,izo;
	int 	ls,ntr,jtr,ktr,mtr,npv,mtmax;
	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,gx;
	float ***mig,***ttab,**tb,**pb,**cs0b,**angb,**tsum,**tt;
	float **tvsum=NULL,***mig1=NULL,***cs=NULL,***tv=NULL,
		**cssum=NULL;
	
	char *datain="stdin",*dataout="stdout",*ttfile,*jpfile,*tvfile,
	     *csfile,*dataout1;
	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("datain",&datain)) {
		infp = stdin;
	} else  
		if ((infp=fopen(datain,"r"))==NULL)
			err("cannot open datain=%s\n",datain);
	if( !getparstring("dataout",&dataout)) {
		outfp = stdout;
	} else  
		outfp = fopen(dataout,"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");
	ext = fxt+(nxt-1)*dxt;
	ezt = fzt+(nzt-1)*dzt;
	es = fs+(ns-1)*ds;

	/* 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 = fxo+(nxo-1)*dxo;
	ezo = fzo+(nzo-1)*dzo;
	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=%g fxo=%g ext=%g exo=%g fzt=%g fzo=%g ezt=%g ezo=%g",
                    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("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("dataout1",&dataout1))
			dataout1="dataout1";
 		out1fp = fopen(dataout1,"w");
	}

	fprintf(jpfp,"\n");
	fprintf(jpfp," Migration parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," datain=%s \n",datain);
	fprintf(jpfp," dataout=%s \n",dataout);
	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);
	if(npv)
 	  fprintf(jpfp," tvfile=%s csfile=%s dataout1=%s\n",
		tvfile,csfile,dataout1);
	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],(int) '\0',noff*nxo*nzo*sizeof(float)); 
	 if(npv)
		memset((void *) mig1[0][0], (int)'\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);
        fprintf(jpfp," Must NOT satisfy this: if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && MAX(gx-sx,sx-gx)<=offmax )\n");

	do {
		/* determine offset index	*/
	    float as,res;
	    sx = tr.sx/1000;
	    gx = tr.gx/1000;
	    io = (int)((gx-sx-off0)/doff+0.5);
	    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, (int) '\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); 

		   if(npv){
			memcpy((void *) tro.data,
			       (const void *) mig1[io][ixo],nzo*sizeof(float));
			for(izo=0; izo<nzo; ++izo)
				tro.data[izo] *=scal;
			/* write out */
			fputtr(out1fp,&tro); 
		    }		 
		}
	}

	fprintf(jpfp," \n");
	fprintf(jpfp," output done\n");
	fflush(jpfp);

	efclose(jpfp);
	efclose(outfp);

	    
	free2float(tsum);
	free2float(tt);
	free2float(pb);
	free2float(tb);
	free2float(cs0b);
	free2float(angb);
	free3float(ttab);
	free3float(mig);
	free3float(mig1);
	if(npv){
		free3float(tv);