sukdsyn2d.c

su 的源代码库

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

/* SUKDSYN2D: $Revision: 1.10 $ ; $Date: 2006/11/07 22:58:42 $	*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
" SUKDSYN2D - Kirchhoff Depth SYNthesis of 2D seismic data from a	",
"             migrated seismic section					",
" 									",
"   sukdsyn2d  infile  outfile [parameters] 				",
"									",
" Required parameters: 							",
" infile=stdin		input migrated section				",
" outfile=stdout	file for output seismic traces  		",
" ttfile		file for input traveltime tables		",
"									",
" The following 9 parameters describe traveltime tables:		",
" fzt= 			first depth sample 				",
" nzt= 			number of depth samples 			",
" dzt=			depth interval 					",
" fxt=			first lateral sample 				",
" nxt=			number of lateral samples 			",
" dxt=			lateral interval 				",
" fs=			x-coordinate of first source			",
" ns=			number of sources				",
" ds=			x-coordinate increment of sources		",
"									",
" The following 6 parameters describe the migration section:		",
" fz=                   first z-coordinate in migrated section 		",
" dz=     		vertical spacing of migrated section 		",
" nz=           	number of depth points in migrated section	",
" fx=                   first x-coordinate of migrated section 		",
" dx=     		horizontal spacing of migrated section 		",
" nx=           	number of lateral points in migrated section  	",
"									",
" Optional Parameters:							",
" nt=501        	number of time samples				",
" dt=0.004      	time sampling interval (sec)			",
" ft=0.0        	first time (sec)				",
" nxo=1                 number of source-receiver offsets		",
" dxo=25                offset sampling interval  			",
" fxo=0.0               first offset  					",
" nxs=101               number of shotpoints  				",
" dxs=25                shotpoint sampling interval  			",
" fxs=0.0               first shotpoint 				",
" fmax=1/(4*dt)         maximum frequency in migration section (Hz)	",
" aperx=nxt*dxt/2  	modeling lateral aperature 			",
" angmax=60		modeling 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	",
"									",
" Notes:								",
" This program takes a migrated seismic section and a set of travel time",
" tables generated using rayt2d for a specific background velocity model",
" and generates synthetic seismic data in the form of common shot gathers.",
" (Common offset gathers may be generated by using nxo=1.) (Demigration.)",
"									",
" This program is a tool which may be used for the migration residual	",
" statics estimation technique of Tjan, Audebert, and Larner 1994.	",
"									",
"1. The traveltime tables are generated by program rayt2d (or other ones)",
"   on relatively coarse grids, with dimension ns*nxt*nzt. In the	",
"   modeling process, traveltimes are interpolated into shot/geophone 	",
"   positions and migration section grids. 				",
"2. The input migration section must be an array of binary floats (no SU",
"   headers).								", 
"3. The synthesized traces are output in common-shot gathers in SU format.",
"4. The memory requirement for this program is about			",
"    	(ns*nxt*nzt+nx*nz+4*nr*nzt+3*nxt*nzt)*4 bytes 			",
"    where nr = 1+min(nxt*dxt,0.5*offmax+aperx)/dxo. 			",
"									",
NULL};
/*
 * Author:  CWP: Zhenyue Liu, 07/24/95,  Colorado School of Mines 
 *
 * References: 
 *
 * Tjan, T., F. Audebert, and K. Larner, 1994,
 *    Prestack migration for residual statics estimation in complex media
 *    (Appeared in 1994 Project Review, CWP-153.)
 *
 * Tjan, T., 1995, Residual statics estimation for data from structurally
 *    complex areas using prestack depth migration: M.Sc. thesis, Colorado
 *    School of Mines. (In progress.)
 *
 * Larner, K., and Tjan, T., 1995, Simultaneous statics and velocity
 *    estimation for data from structurally complex areas.
 *    (Appeared in 1995 Project Review, CWP-185.)
 *
 *
 * Trace header fields set: ns, dt, delrt, tracl, tracr, fldr, tracf
 *                          offset, sx, gx, trid, counit
 */

/**************** end self doc ***********************************/
 
/* prototypes of functions used internally */
void integ(float **mig,int nz,float dz,int nx,int m,float **migi);
void resit(int nx,float fx,float dx,int nz,int nr,float dr,
	float **tb,float **t,float x0);
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 **sigb,float **cosb);
void maketrace(float *trace,int nt,float ft,float dt,
	float xs,float xg,float **mig,float **migi,float aperx,
  	int nx,float fx,float dx,int nz,float fz,float dz,
	int mzmax,int ls,float angmax,float v0,float fpeak,Wavelet *w,
	float **tb,float **pb,float **sigb,float **cosb,int nr,float **tsum,
	int nxt,float fxt,float dxt,int nzt,float fzt,float dzt);

/* parameters for half-derivative filter */
#define LHD 20
#define NHD 1+2*LHD

/* segy trace */
segy tr;

int
main (int argc, char **argv)
{
	int 	nt,nzt,nxt,nz,nxo,nxs,ns,nx,nr,is,ixo,ixs;
	int 	ls,jtr,mtr,tracl,mzmax;
	off_t nseek;
	float   ft,fzt,fxt,fz,fx,fxo,fs,fxs,dt,dzt,dxt,dz,dx,dxo,ds,dxs,
		ext,ezt,es,ex,ez,xo,s,xs,xg,fpeak;	
	float 	v0,dvz;
	float 	fmax,angmax,offmax,rmax,aperx;
	float **mig,**migi,***ttab,**tb,**pb,**cosb,**sigb,**tsum,**tt;
	
	char *infile="stdin",*outfile="stdout",*ttfile,*jpfile;
	FILE *infp,*outfp,*ttfp,*jpfp;
	Wavelet *w;


	/* 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);
	warn("Got A");
	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 for seismogram traces */
	if (!getparint("nt",&nt)) nt = 501;
	if (!getparfloat("dt",&dt)) dt = 0.004;
	if (!getparfloat("ft",&ft)) ft = 0.0;
	if (!getparfloat("fmax",&fmax)) fmax = 1.0/(4*dt);
	fpeak = 0.2/dt;
	if (!getparint("nxs",&nxs)) nxs = 101;
	if (!getparfloat("dxs",&dxs)) dxs = 15;
	if (!getparfloat("fxs",&fxs)) fxs = 0.0;
	if (!getparint("nxo",&nxo)) nxo = 1;
	if (!getparfloat("dxo",&dxo)) dxo = 50;
	if (!getparfloat("fxo",&fxo)) fxo = 0.0;
	
	/* get traveltime table 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;

	/* check source and receiver coordinates */
 	for (ixs=0; ixs<nxs; ++ixs) {
		xs = fxs+ixs*dxs;
		for (ixo=0; ixo<nxo; ++ixo) {
			xg = xs+fxo+ixo*dxo;
			if (fs>xs || es<xs || fs>xg || es<xg) 
		err("shot or receiver lie outside of specified (x,z) grid\n");
		}
	} 

	/* get migration section parameters	*/
	if (!getparint("nx",&nx)) err("must specify nx!\n");
	if (!getparfloat("fx",&fx)) err("must specify fx!\n");
	if (!getparfloat("dx",&dx)) err("must specify dx!\n");
	if (!getparint("nz",&nz)) err("must specify nz!\n");
	if (!getparfloat("fz",&fz)) err("must specify fz!\n");
	if (!getparfloat("dz",&dz)) err("must specify dz!\n");
	ex = fx+(nx-1)*dx;
	ez = fz+(nz-1)*dz;
	if(fxt>fx || ext<ex || fzt>fz || ezt<ez) 
		err(" migration section 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");
	mzmax = dx*sin(angmax*PI/180.)/dz;
	if(mzmax<1) mzmax = 1;
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;

	if (!getparint("ls",&ls))	ls = 1;
	if (!getparint("mtr",&mtr))	mtr = 100;

	fprintf(jpfp,"\n");
	fprintf(jpfp," Modeling 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," nz=%d fz=%g dz=%g\n",nz,fz,dz);
	fprintf(jpfp," nx=%d fx=%g dx=%g\n",nx,fx,dx);
	fprintf(jpfp," \n");
	fprintf(jpfp," nxs=%d fxs=%g dxs=%g\n",nxs,fxs,dxs);
	fprintf(jpfp," nxo=%d fxo=%g dxo=%g\n",nxo,fxo,dxo);
	fprintf(jpfp," \n");
	
	/* compute reference traveltime and slowness  */
	offmax = MAX(ABS(fxo),ABS(fxo+(nxo-1)*dxo));
	rmax = MAX(es-fxt,ext-fs);
	rmax = MIN(rmax,0.5*offmax+aperx);
	nr = 2+(int)(rmax/dx);
	tb = ealloc2float(nzt,nr);
	pb = ealloc2float(nzt,nr);
	sigb = ealloc2float(nzt,nr);
	cosb = ealloc2float(nzt,nr);
	timeb(nr,nzt,dx,dzt,fzt,dvz,v0,tb,pb,sigb,cosb);

	fprintf(jpfp," nt=%d ft=%g dt=%g fpeak=%g \n",nt,ft,dt,fpeak);
	fprintf(jpfp," v0=%g dvz=%g \n",v0,dvz);
 	fprintf(jpfp," aperx=%g angmax=%g offmax=%g\n",aperx,angmax,offmax);
 	fprintf(jpfp," mtr=%d ls=%d\n",mtr,ls);
	fprintf(jpfp," ================\n");
	fflush(jpfp);

	/* allocate space */
	mig = ealloc2float(nz,nx);
	migi = ealloc2float(nz+2*mzmax,nx);
	ttab = ealloc3float(nzt,nxt,ns);
	tt = ealloc2float(nzt,nxt);
	tsum = ealloc2float(nzt,nxt);


	/* make wavelet */
	makericker(fpeak,dt,&w);

	/* set constant segy trace header parameters */
	memset((void *) &tr, 0, sizeof(segy));
	tr.trid = 1;
	tr.counit = 1;
	tr.ns = nt;
	tr.dt = 1.0e6*dt;
	tr.delrt = 1.0e3*ft;

	fprintf(jpfp," read traveltime tables \n");
	fflush(jpfp);

	/* 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,dx,tb,ttab[is],s);		
	}

	fprintf(jpfp," read migration section \n");
	fflush(jpfp);

	/* read migration section	*/
	fread(mig[0],sizeof(float),nx*nz,infp);

	/* integrate migration section for constructing anti-aliasing 
		filter	*/
	integ(mig,nz,dz,nx,mzmax,migi);
                       
	fprintf(jpfp," start synthesis ... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
	
	jtr = 0;
	/* loop over shots  */
	for (ixs=0,xs=fxs,tracl=0; ixs<nxs; ++ixs,xs+=dxs) {
		/* loop over offsets */
		for (ixo=0,xo=fxo; ixo<nxo; ++ixo,xo+=dxo) {
	    		float as,res;
	    		int is;
			xg = xs+xo; 
			/* set segy trace header parameters */
			tr.tracl = tr.tracr = ++tracl;
			tr.fldr = 1+ixs;
			tr.tracf = 1+ixo;
			tr.offset = NINT(xo);
			tr.sx = NINT(xs);
			tr.gx = NINT(xg);

	    		as = (xs-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);

	    		as = (xg-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);
				fflush(jpfp);

			/* make one trace */
			maketrace(tr.data,nt,ft,dt,xs,xg,mig,migi,aperx,
			  nx,fx,dx,nz,fz,dz,mzmax,ls,angmax,v0,fmax,w,
			  tb,pb,sigb,cosb,nr,tsum,nxt,fxt,dxt,nzt,fzt,dzt);
			
			/* write trace */
			fputtr(outfp,&tr);

	        	jtr++;
	        	if((jtr-1)%mtr ==0 ){
				fprintf(jpfp," generated trace %d\n",jtr);
				fflush(jpfp);
	    		}
		}
	}


	fprintf(jpfp," generated %d traces in total\n",jtr);


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

	efclose(jpfp);
	efclose(outfp);

	    
	free2float(tsum);
	free2float(tt);