sutifowler.c

su 的源代码库

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

/* SUTIFOWLER: $Revision: 1.3 $ ; $Date: 2006/08/02 17:44:28 $	*/

#include "su.h"
#include "segy.h"
#include "VND.h"
#define LTABLE 8	/* number of coef in phased sinc table */
#define NTABLE 513	/* number of table entries in phased sinc table */
#define NP 10003	/* number of entries in (velocity,slowness) tables */

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SUTIFOWLER   VTI constant velocity prestack time migration		",
"	      velocity analysis via Fowler's method			",
"									",
" sutifowler ncdps=250 vmin=1500 vmax=6000 dx=12.5			",
"									",
" Required Parameter:							",
" ncdps=		number of input cdp's				",
" Optional Parameters:							",
" choose=1	1 do full prestack time migration			",
"		2 do only DMO						",
"		3 do only post-stack migrations				",
"		4 do only stacking velocity analysis			",
" getcvstacks=0	flag to set to 1 if inputting precomputed cvstacks	",
"		(vmin, nvstack, and ncdps must match SUCVS4FOWLER job)	",
" vminstack=vmin	minimum velocity panel in m/s in input cvstacks	",
" etamin=0.		minimum eta (see paper by Tariq Alkhalifah)	",
" etamax=0.5	maximum eta (see paper by Tariq Alkhalifah)		",
" neta=1	number of eta values to image				",
" d=0.		Thomsen's delta						",
" vpvs=0.5	assumed vp/vs ratio (not critical -- default almost always ok)",
" dx=25.	cdp x increment						",
" vmin=1500.	minimum velocity panel in m/s to output			",
" vmax=8000.	maximum velocity panel in m/s to output			",
" nv=75	 number of velocity panels to output				",
" nvstack=180	number of stacking velocity panels to compute		",
"		     ( Let offmax be the maximum offset, fmax be	",
"		     the maximum freq to preserve, and tmute be		",
"		     the starting mute time in sec on offmax, then	",
"		     the recommended value for nvstack would be		",
"		     nvstack = 4 +(offmax*offmax*fmax)/(0.6*vmin*vmin*tmute)",
"		     ---you may want to make do with less---)		",
" nxpad=0	  number of traces to padd for spatial fft		",
"		     Ideally nxpad = (0.5*tmax*vmax+0.5*offmax)/dx	",
" lmute=24	 length of mute taper in ms				",
" nonhyp=1	  1 if do mute at 2*offset/vmin to avoid non-hyperbolic ",
"				moveout, 0 otherwise			",
" lbtaper=0	length of bottom taper in ms				",
" lstaper=0	length of side taper in traces				",
" dtout=1.5*dt	output sample rate in s,   note: typically		",
"				fmax=salias*0.5/dtout			",
" mxfold=120	maximum number of offsets/input cmp			",
" salias=0.8	fraction of output frequencies to force within sloth	",
"		     antialias limit.  This controls muting by offset of",
"		     the input data prior to computing the cv stacks	",
"		     for values of choose=1 or choose=2.		",
" file=sutifowler	root name for temporary files			",
" p=not		Enter a path name where to put temporary files.		",
"	  	specified  Can enter multiple times to use multiple disk",
"		systems.						",
"		     The default uses information from the .VND file	",
"		     in the current directory if it exists, or puts 	",
"		     unique temporary files in the current directory.	",
" ngroup=20	Number of cmps per velocity analysis group.		",
" printfile=stderr    The output file for printout from this program.	",
"									",
" Required trace header words on input are ns, dt, cdp, offset.		",
" On output, trace headers are rebuilt from scratch with		",
" ns - number of samples						",
" dt - sample rate in usec						",
" cdp - the output cmp number (0 based)					",
" offset - the output velocity						",
" tracf	- the output velocity index (0 based)				",
" fldr - index for velocity analysis group (0 based, groups of ngroup cdps)",
" ep - central cmp for velocity analysis group				",
" igc - index for choice of eta (0 based)				",
" igi - eta*100								",
" sx=gx	- x coordinate as icmp*dx					",
" tracl=tracr -sequential trace count (1 based)				",
"									",
" Note: Due to aliasing considerations, the small offset-to-depth	",
" ratio assumption inherent in the TI DMO derivation, and the		",
" poor stacking of some large-offset events associated with TI non-hyperbolic",
" moveout, a fairly stiff initial mute is recommended for the		",
" long offsets.  As a result, this method may not work well		",
" where you have multiple reflections to remove via stacking.		",
" 									",
" Note: The temporary files can be split over multiple disks by building",
" a .VND file in your working directory.  The .VND file is ascii text	",
" with the first line giving the number of directories followed by	",
" successive lines with one line per directory name.			",
"									",
" Note: The output data order has primary key equal to cdp, secondary	",
" key equal to eta, and tertiary key equal to velocity.			",
NULL};

/* Credits:
 *	CWP: John Anderson (visitor to CSM from Mobil) Spring 1993
 *
 */
/**************** end self doc ********************************/


static void cvstack(VND *vnda, VND *vnd, int icmp, int noff, float *off,
		float *mute, int lmute, int nv, float *p2,
		float dt, float dtout);
static void vget( float a, float b, float e, float d,
		float theta, float *vel);
VND *ptabledmo(int nv, float *v, float etamin, float deta, int neta,
		float d, float vsvp, int np, float dp, float dp2, char *file);
VND *ptablemig(int nv, float *v, float etamin, float deta,
		int neta, float d, float vsvp, int np, char *file);
static void taper (int lxtaper, int lbtaper,
		int nx, int ix, int nt, float *trace);

segy tr;	/* input and output SEGY data */
FILE *fpl;	/* file pointer for print listing */
int main(int argc, char **argv)
{
	VND *vnd=NULL;	/* big file holding data, all cmps, all etas, all velocities */
	VND *vnda=NULL;	/* holds one input cmp gather */
	VND *vndb=NULL;	/* holds (w,v) for one k component */
	VND *vndvnmo=NULL;	/* holds (vnmo,p) table for ti dmo */
	VND *vndvphase=NULL;	/* holds (vphase,p) table for ti Stolt migration */
	long N[2];	/* holds number of values in each dimension for VND opens */
	long key[2];	/* holds key in each dimension for VND i/o */
	char **dir=NULL; /* could hold list of directories where to put VND temp files */
	char *file;	/* root name for temporary files */
	char *printfile; /* name of file for printout */
	complex *crt;
	complex *ctemp;
	complex czero;
	float *rt;
	char *ccrt;
	char *fname;
	
	float etamin;	/* minimum eta scan to compute */
	float etamax;	/* maximum eta scan to compute */
	float deta;	/* increment in eta to compute for eta scan */
	float dx;	/* cmp spatial sampling interval */
	float dk;	/* wavenumber increment */
	float dv;	/* velocity increment */
	float vmin;	/* minimum output velocity */
	float vmax;	/* maximum output velocity */
	float dt;	/* input sample rate in seconds */
	float dtout;	/* output sample rate in seconds */
	float *mute;	/* array of mute times for this cmp */
	float *off;	/* array of offsets for this cmp */
	float *v;	/* array of output velocities */
	float *p2stack;	/* array of stacking 1/(v*v) */
	float *rindex;	/* array of interpolation indices */
	float dp2=0.0;	/* increment in slowness squared for input cvstacks */
	float scale;	/* used for trace scale factor */
	float p;	/* horizontal slowness */
	float p2;	/* p*p */
	float v2;	/* velocity squared */
	float ak;	/* horizontal wavenumber */
	float dw;	/* angular frequency increment */
	float *w;	/* array holding w values for Fowler */
	float factor;	/* scale factor */
	float d;	/* Thomsen's delta */
	float vsvp;	/* vs/vp ratio */
	float dp;	/* increment of slowness values in vndvnmo table */
	float rp;	/* real valued index in p */
	float wgt;	/* weight for linear interpolation */
	float fmax;	/* maximum frequency to use for antialias mute */
	float salias;	/* fraction of frequencies to be within sloth antialias limit */
	float dpm;	/* slowness increment in TI migration table */
	float fw;	/* first w in Stolt data table */
	float vminstack;/* only used if reading precomputed cvstacks, minimum stacking vel */

	int neta;	/* number of eta scans to compute */
	int ichoose;	/* defines type of processing to do */
	int ncmps;	/* number of input and output cmps */
	int nv;		/* number of output velocity panels to generate */
	int nvstack;	/* number of cvstack panels to generate */
	int ntpad;	/* number of time samples to padd to avoid wraparound */
	int nxpad;	/* number of traces to padd to avoid wraparound */
	int lmute;	/* number of samples to taper mute */
	int lbtaper;	/* length of bottom time taper in ms */
	int lstaper;	/* length of side taper in traces */
	int mxfold;	/* maximum allowed number of input offsets/cmp */
	int icmp;	/* cmp index */
	int ntfft;	/* length of temporal fft for Fowler */
	int ntffts;	/* length of temporal fft for Stolt */
	int nxfft;	/* length of spatial fft */
	int ntfftny;	/* count of freq to nyquist */
	int nxfftny;	/* count of wavenumbers to nyquist */
	int nmax;	/* used to compute max number of samples for array allocation */
	int oldcmp;	/* current cdp header value */
	int noff;	/* number of offsets */
	int k;		/* wavenumber index */
	int iwmin;	/* minimum freq index */
	int TI;		/* 0 for isotropic, 1 for transversely isotropic */
	long it;	/* time index */
	long iw;	/* index for angular frequency */
	long nt;	/* number of input time samples */
	long ntout;	/* number of output time samples */
	long iv;	/* velocity index */
	long ip;	/* slowness index */
	long ieta;
	int nonhyp;	/* flag equals 1 if do mute to avoid non-hyperbolic events */
	int getcvstacks;/* flag equals 1 if input cvstacks precomputed */
	int ngroup;	/* number of traces per vel anal group */
	int ndir;	/* number of user specified directories for temp files */

/******************************************************************************/
/* 	input parameters, allocate buffers, and define reusable constants     */
/******************************************************************************/
	initargs(argc, argv);
	requestdoc(1);

	/* get first trace and extract critical info from header */
	if(!gettr(&tr)) err("Can't get first trace \n");
	nt=tr.ns;
	dt=0.000001*tr.dt;
	oldcmp=tr.cdp;

	if (!getparstring("printfile",&printfile)) printfile=NULL;
	if (printfile==NULL) {
		fpl=stderr;
	}else{
		fpl=fopen(printfile,"w");
	}

	if (!getparfloat("salias",&salias)) salias=0.8;
 	if(salias>1.0) salias=1.0;
	if (!getparfloat("dtout",&dtout)) dtout=1.5*dt;
	ntout=1+nt*dt/dtout;
	if (!getparint("getcvstacks",&getcvstacks)) getcvstacks=0;
	if(getcvstacks) {
		dtout=dt;
		ntout=nt;
	}
	fmax=salias*0.5/dtout;
	fprintf(fpl,"sutifowler: ntin=%ld dtin=%f\n",nt,dt);
	fprintf(fpl,"sutifowler: ntout=%ld dtout=%f\n",ntout,dtout);
	if (!getparstring("file",&file)) file="sutifowler";
	if (!getparfloat("dx",&dx)) dx=25.;
	if (!getparfloat("vmin",&vmin)) vmin=1500.;
	if (!getparfloat("vmax",&vmax)) vmax=8000.;
	if (!getparfloat("vminstack",&vminstack)) vminstack=vmin;
	if (!getparfloat("d",&d)) d=0.0;
	if (!getparfloat("etamin",&etamin)) etamin=0.0;
	if (!getparfloat("etamax",&etamax)) etamax=0.5;
	if (!getparfloat("vsvp",&vsvp)) vsvp=0.5;
	if (!getparint("neta", &neta)) neta = 1;
	if (fabs(etamax-etamin)<1.0e-7) neta = 1;
	if (neta < 1) neta = 1;
	if (!getparint("choose", &ichoose)) ichoose = 1;
	if (!getparint("ncdps", &ncmps)) err("sutifowler: must enter ncdps");
	if (!getparint("nv", &nv)) nv = 75;
	if (!getparint("nvstack", &nvstack)) nvstack = 180;
	if (!getparint("ntpad", &ntpad)) ntpad = 0.1*ntout;
	if (!getparint("nxpad", &nxpad)) nxpad = 0;
	if (!getparint("lmute", &lmute)) lmute = 24;
	lmute=1 + 0.001*lmute/dtout;
	if (!getparint("lbtaper", &lbtaper)) lbtaper = 0;
	if (!getparint("lstaper", &lstaper)) lstaper = 0;
	if (!getparint("mxfold", &mxfold)) mxfold = 120;
	if (!getparint("nonhyp",&nonhyp)) nonhyp=1.;
	if (!getparint("ngroup", &ngroup)) ngroup = 20;
	ndir = countparname("p");
	if(ndir==0) {
		ndir=-1;
	}else{
		dir = (char **)VNDemalloc(ndir*sizeof(char *),"dir");
		for(k=0;k<ndir;k++) {
			it=getnparstring(k+1,"p",&dir[k]);
		}
	}
	lbtaper=lbtaper/(1000.*dt);
	TI=0;
	if(fabs(d)>0. || fabs(etamin)>0 || neta>1 ) TI=1;
	if(TI) fprintf(fpl,"sutifowler: operation in TI mode\n");
	deta = 0.;
	if(neta>1) deta=(etamax-etamin)/(neta-1);
	dp=1./(vmin*(NP-5));
	if(TI) dp=dp*sqrt(1.+2.*fabs(etamin));
	if(ichoose>2) nvstack=nv;
	if(ichoose==1 || ichoose==2 || ichoose==3) {
		ntfft=ntout+ntpad;
	}else{
		ntfft=1;
	}
	if(ichoose==1 || ichoose==3) {
		ntffts=2*ntout/0.6;
	}else{
		ntffts=1;
	}
	ntfft=npfao(ntfft,2*ntfft);
	ntffts=npfao(ntffts,2*ntffts);
	dw=2.*PI/(ntfft*dtout);
	nxfft=npfar(ncmps+nxpad);
	dk=2.*PI/(nxfft*dx);
	fprintf(fpl,"sutifowler: ntfft=%d ntffts=%d nxfft=%d\n",ntfft,ntffts,nxfft);
	czero.r=czero.i=0.;
	scale=1.;
	if(ichoose<5) scale=1./(nxfft);
	if(ichoose==1 || ichoose==2 ) scale*=1./ntfft;
	if(ichoose==1 || ichoose==3 ) scale*=1./ntffts;
	nxfftny = nxfft/2 + 1;
	ntfftny = ntfft/2 + 1;
	nmax = nxfftny;
	if(ntfft > nmax) nmax=ntfft;
	if((NP/2+1)>nmax) nmax=(NP/2+1);
	if(nvstack>nmax) nmax=nvstack;
	if(nv*neta>nmax) nmax=nv*neta;
	ctemp = (complex *)VNDemalloc(nmax*sizeof(complex),"allocating ctemp");
	rindex=(float *)VNDemalloc(nmax*sizeof(float),"allocating rindex");
	if(ntffts > nmax) nmax=ntffts;
	crt = (complex *)VNDemalloc(nmax*sizeof(complex),"allocating crt");
	rt = (float *)crt;
	ccrt = (char *)crt;
	fprintf(fpl,"sutifowler: nv=%d nvstack=%d\n",nv,nvstack);
	v=(float *)VNDemalloc(nv*sizeof(float),"allocating v");
	p2stack=(float *)VNDemalloc(nvstack*sizeof(float),"allocating p2stack");
	mute=(float *)VNDemalloc(mxfold*sizeof(float),"allocating mute");
	off=(float *)VNDemalloc(mxfold*sizeof(float),"allocating off");
	fprintf(fpl,"sutifowler: allocating and filling w array\n");
	w=(float *)VNDemalloc(ntfft*sizeof(float),"allocating w");
	for(iw=0;iw<ntfft;iw++) {
		if(iw<ntfftny){
			w[iw]=iw*dw;
		}else{
			w[iw]=(iw-ntfft)*dw;
		}
		if(iw==0) w[0]=0.1*dw;  	/* fudge for dc component */
	}

/******************************************************************************/
	fprintf(fpl,"sutifowler: building function for stacking velocity analysis\n");
/******************************************************************************/
	dv=(vmax-vmin)/MAX((nv-1),1);
	for(iv=0;iv<nv;iv++) v[iv]=vmin+iv*dv;
	if(ichoose>=3){
	  	for(iv=0;iv<nvstack;iv++) {
			p2stack[iv]=1./(v[iv]*v[iv]);
			fprintf(fpl,"	    stacking velocity %ld %f\n",iv,v[iv]);
		}
	}else{
		if(nvstack<6) err("sutifowler: nvstack must be 6 or more");
		dp2 = 1./(vminstack*vminstack*(nvstack-5));
		for(iv=0;iv<nvstack;iv++) {
			p2stack[iv]=iv*dp2;
			if(iv>0) {
				factor=1./sqrt(p2stack[iv]);
				fprintf(fpl,"	    stacking velocity %ld %f\n",iv,factor);
			}else{