sudmofkcw.c

su 的源代码库

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

/* SUDMOFKCW: $Revision: 1.5 $ ; $Date: 2003/06/09 16:17:07 $	*/

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

/*********************** self documentation ******************************/
char *sdoc[] = {
" 									",
" SUDMOFKCW - converted-wave DMO via F-K domain (log-stretch) method for",
" 		common-offset gathers					",
" 									",
" sudmofkcw <stdin >stdout cdpmin= cdpmax= dxcdp= noffmix= [...]	",
" 									",
" Required Parameters:							",
" cdpmin		  minimum cdp (integer number) for which to apply DMO",
" cdpmax		  maximum cdp (integer number) for which to apply DMO",
" dxcdp		   distance between adjacent cdp bins (m)		",
" noffmix		 number of offsets to mix (see notes)		",
" 									",
" Optional Parameters:							",
" tdmo=0.0		times corresponding to rms velocities in vdmo (s)",
" vdmo=1500.0		rms velocities corresponding to times in tdmo (m/s)",
" gamma=0.5		 velocity ratio, upgoing/downgoing		",
" ntable=1000		 number of tabulated z/h and b/h (see notes)	",
" sdmo=1.0		DMO stretch factor; try 0.6 for typical v(z)	",
" flip=0		 =1 for negative shifts and exchanging s1 and s2",
" 			 (see notes below)				",
" fmax=0.5/dt		maximum frequency in input traces (Hz)		",
" verbose=0		=1 for diagnostic print				",
" 									",
" Notes:								",
" Input traces should be sorted into common-offset gathers.  One common-",
" offset gather ends and another begins when the offset field of the trace",
" headers changes.							",
" 									",
" The cdp field of the input trace headers must be the cdp bin NUMBER, NOT",
" the cdp location expressed in units of meters or feet.		",
" 									",
" The number of offsets to mix (noffmix) should typically equal the ratio of",
" the shotpoint spacing to the cdp spacing.  This choice ensures that every",
" cdp will be represented in each offset mix.  Traces in each mix will	",
" contribute through DMO to other traces in adjacent cdps within that mix.",
" 									",
" The tdmo and vdmo arrays specify a velocity function of time that is	",
" used to implement a first-order correction for depth-variable velocity.",
" The times in tdmo must be monotonically increasing. The velocity function",
" is assumed to have been gotten by traditional NMO. 			",
" 									",
" For each offset, the minimum time at which a non-zero sample exists is",
" used to determine a mute time.  Output samples for times earlier than this", 
" mute time will be zeroed.  Computation time may be significantly reduced",
" if the input traces are zeroed (muted) for early times at large offsets.",
" 									",
" z/h is horizontal-reflector depth normalized to half source-reciver offset",
" h.  Normalized shift of conversion point is b/h.  The code now does not",
" support signed offsets, therefore it is recommended that only end-on data,",
" not split-spread, be used as input (of course after being sorted into	",
" common-offset gathers).  z/h vs b/h depends on gamma (see Alfaraj's Ph.D.",
" thesis, 1993).							",
" 									",
" Flip factor = 1 implies positive shift of traces (in the increasing CDP",
" bin number direction).  When processing split-spread data, for example,",
" if one side of the spread is processed with flip=0, then the other side",
" of the spread should be processed with flip=1.  The flip factor also	",
" determines the actions of the factors s1 and s2, i.e., stretching or	",
" squeezing.								",
" 									",
" Trace header fields accessed:  nt, dt, delrt, offset, cdp.		",
"									",
NULL};

/* Credits:
 *	CWP: Mohamed Alfaraj
 *		Dave Hale
 *
 * Technical Reference:
 *	Transformation to zero offset for mode-converted waves
 *	Mohammed Alfaraj, Ph.D. thesis, 1993, Colorado School of Mines
 *
 *	Dip-Moveout Processing - SEG Course Notes
 *	Dave Hale, 1988
 */
/**************** end self doc *******************************************/

/* Prototypes of functions used internally */
static void mkvrms (int ntdmo, float *tdmo, float *vdmo,
	int nt, float dt, float ft, float *vrms);
static void dmooff (float offset, float fmax, int nx, float dx,
	int nt, float dt, float ft, float *vrms, float **ptx, float gamma,
	float *boh, float *zoh, int ntable, float s1, float s2, float sign);
static void maketu (float offset, int itmin, float fmax,
	int nt, float dt, float ft, float *vrms, float **uoftp,
	int *nup, float *dup, float *fup, float **tofup, float *tconp);
static float getbbohitn (float x, int itmin, int nt, float dt, float *v,
	int ntable, float *boh, float *zoh, float gamma,
	float **bbohp, int **t0p);
static void table (int ntable, float gamma, float *zoh, float *boh);
static void stretchfactor (float sdmo, float gamma, float *s1, float *s2);

/* Globals */
segy tr,tro;

int
main(int argc, char **argv)
{
	int nt;		/* number of time samples per trace */
	float dt;	/* time sampling interval */
	float ft;	/* time of first sample */
	int it;		/* time sample index */
	int cdpmin;	/* minimum cdp to process */
	int cdpmax;	/* maximum cdp to process */
	float dx;	/* cdp sampling interval */
	int nx;		/* number of cdps to process */
	int nxfft;	/* number of cdps after zero padding for fft */
	int nxpad;	/* minimum number of cdps for zero padding */
	int ix;		/* cdp index, starting with ix=0 */
	int noffmix;	/* number of offsets to mix */
	float *tdmo;	/* times at which rms velocities are specified */
	float *vdmo;	/* rms velocities at times specified in tdmo */
	float gamma;	/* upgoing to downging velocity ratio */
	float *zoh=NULL;/* tabulated z/h */
	float *boh=NULL;/* tabulated b/h */
	int ntable;	/* number of tabulated zoh and boh */
	float sdmo;	/* DMO stretch factor */
	float s1;	/* DMO stretch factor */
	float s2;	/* DMO stretch factor */
	float temps;	/* temp value used in excahnging s1 and s2 */
	int flip;	/* apply negative shifts and exchange s1 and s2 */
	float sign; 	/* + if flip=0, negative if flip=1 */
	int ntdmo;	/* number tnmo values specified */
	int itdmo;	/* index into tnmo array */
	int nvdmo;	/* number vnmo values specified */
	float fmax;	/* maximum frequency */
	float *vrms;	/* uniformly sampled vrms(t) */
	float **p;	/* traces for one offset - common-offset gather */
	float **q;	/* DMO-corrected and mixed traces to be output */
	float offset;	/* source-receiver offset of current trace */
	float oldoffset;/* offset of previous trace */
	int noff;	/* number of offsets processed in current mix */
	int ntrace;	/* number of traces processed in current mix */
	int itrace;	/* trace index */
	int gottrace;	/* non-zero if an input trace was read */
	int done;	/* non-zero if done */
	int verbose;	/* =1 for diagnostic print */
	FILE *hfp;	/* file pointer for temporary header file */

	/* hook up getpar */
	initargs(argc, argv);
	requestdoc(1);

	/* get information from the first header */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = tr.dt/1000000.0;
	ft = tr.delrt/1000.0;
	offset = tr.offset;

	/* get parameters */
	if (!getparint("cdpmin",&cdpmin)) err("must specify cdpmin");
	if (!getparint("cdpmax",&cdpmax)) err("must specify cdpmax");
	if (cdpmin>cdpmax) err("cdpmin must not be greater than cdpmax");
	if (!getparfloat("dxcdp",&dx)) err("must specify dxcdp");
	if (!getparint("noffmix",&noffmix)) err("must specify noffmix");
	ntdmo = countparval("tdmo");
	if (ntdmo==0) ntdmo = 1;
	tdmo = ealloc1float(ntdmo);
	if (!getparfloat("tdmo",tdmo)) tdmo[0] = 0.0;
	nvdmo = countparval("vdmo");
	if (nvdmo==0) nvdmo = 1;
	if (nvdmo!=ntdmo) err("number of tdmo and vdmo must be equal");
	vdmo = ealloc1float(nvdmo);
	if (!getparfloat("vdmo",vdmo)) vdmo[0] = 1500.0;
	for (itdmo=1; itdmo<ntdmo; ++itdmo)
		if (tdmo[itdmo]<=tdmo[itdmo-1])
			err("tdmo must increase monotonically");
	if (!getparfloat("gamma",&gamma)) gamma = 0.5;
	if (!getparint("ntable",&ntable)) ntable = 1000;
	if (!getparfloat("sdmo",&sdmo)) sdmo = 1.0;
	if (!getparint("flip",&flip)) flip=0;
	if (flip) 
		sign = -1.0;
	else
		sign = 1.0;
	if (!getparfloat("fmax",&fmax)) fmax = 0.5/dt;
	if (!getparint("verbose",&verbose)) verbose=0;

	/* allocate and generate tables of b/h and z/h if gamma not equal 1 */
	if(gamma != 1.0){
		zoh=alloc1float(ntable);
		boh=alloc1float(ntable);
		table(ntable, gamma, zoh, boh);
	}
	
	/* make uniformly sampled rms velocity function of time */
	vrms = ealloc1float(nt);
	mkvrms(ntdmo,tdmo,vdmo,nt,dt,ft,vrms);
	
	/* determine number of cdps to process */
	nx = cdpmax-cdpmin+1;
	
	/* allocate and zero common-offset gather p(t,x) */
	nxpad = 0.5*ABS(offset/dx);
	nxfft = npfar(nx+nxpad);
	p = ealloc2float(nt,nxfft+2);
	for (ix=0; ix<nxfft; ++ix)
		for (it=0; it<nt; ++it)
			p[ix][it] = 0.0;

	/* allocate and zero offset mix accumulator q(t,x) */
	q = ealloc2float(nt,nx);
	for (ix=0; ix<nx; ++ix)
		for (it=0; it<nt; ++it)
			q[ix][it] = 0.0;
		
	/* open temporary file for headers */
	hfp = tmpfile();
	
	/* initialize */
	oldoffset = offset;
	gottrace = 1;
	done = 0;
	ntrace = 0;
	noff = 0;

	/* get DMO stretch/squeeze factors s1 and s2 */
	stretchfactor (sdmo,gamma,&s1,&s2);
	if(flip) {
		temps = s1;
		s1 = s2;
		s2 = temps; 
	}

	/* print useful information if requested */
	if (verbose)fprintf(stderr,"stretching factors: s1=%f s2=%f\n",s1,s2);

	/* loop over traces */
	do {
		/* if got a trace, determine offset */
		if (gottrace) offset = tr.offset;
		
		/* if an offset is complete */
		if ((gottrace && offset!=oldoffset) || !gottrace) {
		
			/* do dmo for old common-offset gather */
			dmooff(oldoffset,fmax,nx,dx,nt,dt,ft,vrms,p,
			gamma,boh,zoh,ntable,s1,s2,sign);
			
			/* add dmo-corrected traces to mix */
			for (ix=0; ix<nx; ++ix)
				for (it=0; it<nt; ++it)
					q[ix][it] += p[ix][it];
			
			/* count offsets in mix */
			noff++;
							
			/* free space for old common-offset gather */
			free2float(p);
			
			/* if beginning a new offset */
			if (offset!=oldoffset) {
				
				/* allocate space for new offset */
				nxpad = 0.5*ABS(offset/dx);
				nxfft = npfar(nx+nxpad);
				p = ealloc2float(nt,nxfft+2);
				for (ix=0; ix<nxfft; ++ix)
					for (it=0; it<nt; ++it)
						p[ix][it] = 0.0;
			}
		}
		
		/* if a mix of offsets is complete */
		if (noff==noffmix || !gottrace) {
			
			/* rewind trace header file */
			efseeko(hfp, (off_t) 0,SEEK_SET);
			
			/* loop over all output traces */
			for (itrace=0; itrace<ntrace; ++itrace) {
			
				/* read trace header and determine cdp index */
				efread(&tro,HDRBYTES,1,hfp);
				
				/* get dmo-corrected data */
				memcpy((void *) tro.data,
					(const void *) q[tro.cdp-cdpmin],
					nt*sizeof(float));
				
				/* write output trace */
				puttr(&tro);
			}
			
			/* report */
			if (verbose) 
				fprintf(stderr,"\tCompleted mix of "
					"%d offsets with %d traces\n",
					noff,ntrace);
			
			/* if no more traces, break */
			if (!gottrace) break;
			
			/* rewind trace header file */
			efseeko(hfp, (off_t) 0,SEEK_SET);
			
			/* reset number of offsets and traces in mix */
			noff = 0;
			ntrace = 0;
			
			/* zero offset mix accumulator */
			for (ix=0; ix<nx; ++ix)
				for (it=0; it<nt; ++it)
					q[ix][it] = 0.0;
		}
			
		/* if cdp is within range to process */
		if (tr.cdp>=cdpmin && tr.cdp<=cdpmax) {
	
			/* save trace header and update number of traces */
			efwrite(&tr,HDRBYTES,1,hfp);
			ntrace++;

			/* remember offset */
			oldoffset = offset;

			/* get trace samples */
			memcpy((void *) p[tr.cdp-cdpmin],
				(const void *) tr.data, nt*sizeof(float));
		}

		/* get next trace (if there is one) */
		if (!gettr(&tr)) gottrace = 0;
		
	} while (!done);

	return(CWP_Exit());
}
	
static void mkvrms (int ndmo, float *tdmo, float *vdmo,
	int nt, float dt, float ft, float *vrms)
/*****************************************************************************
make uniformly sampled vrms(t) for DMO
******************************************************************************
Input:
ndmo		number of tdmo,vdmo pairs
tdmo		array[ndmo] of times
vdmo		array[ndmo] of rms velocities
nt		number of time samples
dt		time sampling interval
ft		first time sample

Output:
vrms		array[nt] of rms velocities
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 10/03/91
*****************************************************************************/
{
	int it;
	float t,(*vdmod)[4];
	
		vdmod = (float(*)[4])ealloc1float(ndmo*4);
	cmonot(ndmo,tdmo,vdmo,vdmod);
	for (it=0,t=ft; it<nt; ++it,t+=dt)
		intcub(0,ndmo,tdmo,vdmod,1,&t,&vrms[it]);
	free1float((float*)vdmod);
}

static void dmooff (float offset, float fmax, int nx, float dx,
	int nt, float dt, float ft, float *vrms, float **ptx, float gamma,
	float *boh, float *zoh, int ntable, float s1, float s2, float sign)
/*****************************************************************************
perform dmo in f-k domain for one offset
******************************************************************************
Input:
offset		source receiver offset
fmax		maximum frequency
s1		DMO stretch factor
s2		DMO stretch factor
sign		sign of shift
nx		number of midpoints
dx		midpoint sampling interval
nt		number of time samples
dt		time sampling interval
ft		first time
vrms		array[nt] of rms velocities 
ptx		array[nx][nt] for p(t,x), zero-padded for fft (see notes)

Output:
ptx		array[nx][nt] for dmo-corrected p(t,x)
******************************************************************************
Notes:
To avoid having to allocate a separate work space that is larger than the
array ptx[nx][nt], ptx must be zero-padded to enable Fourier transform from x
to k via prime factor FFT.  nxpad (nx after zero-padding) can be estimated by
	nxpad = 2+npfar(nx+(int)(0.5*ABS(offset/dx)));
where npfar() is a function that returns a valid length for real-to-complex 
prime factor FFT.  ptx[nx] to ptx[nxfft-1] must point to zeros.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 08/08/91
*****************************************************************************/
{
	int nxfft,itmin,nu,nufft,nw,nk,ix,iu,iw,ik,it,iwn,
		iwmin,iwmax,nupad,ikmax,*itn,lnt;
	float dw,dk,tcon,wwscl,wwscl2,scale,scales,kmax,lt,
		amp,phase,fr,fi,pwr,pwi,cphase,sphase,os1,os2,
		wmin,wmax,fftscl,du,fu,w,k,*uoft,*tofu,g1,h,vmin,hk,t,*bboh; 
	complex czero=cmplx(0.0,0.0),**ptk,*pu,*pw;

	/* number of cdps after padding for fft */
	nxfft = npfar(nx+(int)(0.5*ABS(offset/dx)));

	/* get minimum time of first non-zero sample */
	for (ix=0,itmin=nt; ix<nx; ++ix) {
		for (it=0; it<itmin && ptx[ix][it]==0.0; ++it);
		itmin = it;
	}
	
	/* if all zeros, simply return */
	if (itmin>=nt) return;