sudmotivz.c

su 的源代码库

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

/* SUDMOTIVZ: $$Revision: 1.2 $ ; $Date: 2003/06/09 16:17:07 $		*/


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


/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SUDMOTIVZ - DMO for Transeversely Isotropic V(Z) media for common-offset",
"            gathers							",
"									",
" sudmotivz <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:							",
" vnfile=        binary (non-ascii) file containing NMO interval	",
"                  velocities (m/s)					",
" vfile=         binary (non-ascii) file containing interval velocities	(m/s)",
" etafile=       binary (non-ascii) file containing eta interval values (m/s)",
" tdmo=0.0       times corresponding to interval velocities in vdmo (s)	",
" vndmo=1500.0   NMO interval velocities corresponding to times in tdmo (m/s)",
" vdmo=vndmo    interval velocities corresponding to times in tdmo (m/s)",
" etadmo=1500.0  eta interval values corresponding to times in tdmo (m/s)",
" fmax=0.5/dt    maximum frequency in input traces (Hz)			",
" smute=1.5      stretch mute used for NMO correction			",
" speed=1.0      twist this knob for speed (and aliasing)		",
" 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.",
"									",
" vnfile, vfile and etafile should contain the regularly sampled interval",
" values of NMO velocity, velocity and eta respectivily as a		",
" function of time.  If, for example, vfile is not supplied, the interval",
" velocity function is defined by linear interpolation of the values in the",
" tdmo and vdmo arrays.  The times in tdmo must be monotonically increasing.",
" If vfile or vdmo are not given it will be equated to vnfile or vndmo. ",
"									",
" For each offset, the minimum time to process is determined using the	",
" smute parameter.  The DMO correction is not computed for samples that	",
" have experienced greater stretch during NMO.				",
"									",
" Trace header fields accessed:  nt, dt, delrt, offset, cdp.		",
NULL};
/**************** end self doc ********************************/

/* Credits:
 *	CWP: Tariq Alkhalifah, Craig Artley
 *
 * Technical Reference:
 *
 *	Alkhalifah, T., 1994,
 *	Transformation to zero offset in transversely isotropic media,
 *	Center for Wave Phenomena, CWP-162.
 *
 *	Alkhalifah, T., Tsvankin, Ilya, 1994,
 *	Velocity analysis for transversely isotropic media,
 *	published at Center for Wave Phenomena report CWP-145.
 *
 *	Artley, C., 1992,
 *	Dip-moveout processing for depth-variable velocity,
 *	M.S. Thesis, Colorado School of Mines;
 *	also published at Center for Wave Phenomena report CWP-115.
 *	
 */

#define NPTAB 401
#define NITER 10
#define TOL 1.0e-4

/* function prototypes for functions defined below */
void dmooff (float **x, float **a, float **s, float **tau,
	 int nttab, float dt, float ft,
	int nptab, float dptab, float fp, float lp,
	int nx, float dx, float offset, int nt,
	float *vrms, float fmax, float smute, float speed, int verbose,
	float **q);
void dmomap (float **x, float **a, float **tau, float **s, int nt, float dt, float ft,
	int nptab, float dptab, int np, float dp2, float fp, float lp,
	int itn, float tn, float h, float *vrms,
	int verbose, float **t0table);
void dmoapply (int nx, float dx, int nt, float dt, int np, float dp2,
	float h, float **tntable, int *it0min, int *nt0,
	float tmin, float fmax, float **qx);
void jakub1k (float k, float h, int np, float dp2, int nt, float dt, float ft,
	float **tntable, int *it0min, int *nt0, float fmax, complex *qq);
float tmute(float *vrms, int nt, float dt, float ft, float h, float smute);
void getX (float tsg, float h, float p0, float v, float *vars);
void getFJ (float **x, float **a, float **tau, float **s,
	int nt, int np, float dt, float dp, float tsg, float h, float p0,
	float *vars, float *eqns, float **jacob);
void intder (float **f, int nx, float dx, float fx,
	int ny, float dy, float fy, float x, float y,
	float *fxy, float *dfdx, float *dfdy);
void rayvt (float p, int nt, float dt,
	float a1111[], float a3333[], float a1313[], float a1133[],
	float tau[], float x[], float a[], float s[]);

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,*etadmo,*vndmo;/* rms velocities at times specified in tdmo */
	int ntdmo;	/* number tnmo values specified */
	int itdmo;	/* index into tnmo array */
	int nvdmo,netadmo,nvndmo;	/* number vnmo values specified */
	float fmax;	/* maximum frequency */
	float *vt,*etat,*vnt;	/* uniformly sampled v(t) */
	float *vrms;	/* uniformly sampled vrms(t) */
	float **q;	/* traces for one offset - common-offset gather */
	float **qmix;	/* 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 */
	char *vfile=""; /* name of interval velocity file */
	char *vnfile="";/* name of interval NMO velocity file */
	char *etafile="";/* name of interval eta velocity file */
	FILE *hfp;	/* file pointer for temporary header file */
	int nttab;	/* number of times in ray tables */
	int nptab=NPTAB;/* number of ray parameters in ray tables */
	int ip;		/* ray parameter counter */
	float smute;	/* NMO stretch mute */
	float t;	/* time */
	float sum;	/* sum accumulator */
	float lt;	/* last time */
	float lttab;	/* last time in ray tables */
	float p;	/* ray parameter */
	float fp;	/* first ray parameter */
	float lp;	/* last ray parameter */
	float dptab;	/* ray parameter sampling interval in ray tables */
	float speed;	/* speed knob --- twist for speed and aliasing */
	float **x;	/* table [nptab][nttab] of horizontal distance ray */
	float **tau;	/* table [nptab][nttab] of time depth of ray */
	float **a;	/* table [nptab][nttab] of propagation angle of ray */
	float *a1111,*a3333,*a1313,*a1133,delta,vnt2,vt2,**s;

	/* 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;
	if (ft!=0.0) err("ft=%f --- nonzero ft not yet implemented!");
	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");
	if (!getparfloat("fmax",&fmax)) fmax = 0.5/dt;
	if (!getparint("verbose",&verbose)) verbose = 0;
	if (!getparfloat("speed",&speed)) speed = 1.0;
	if (!getparfloat("smute",&smute)) smute = 1.5;

	/* allocate space for velocity arrays */
	nttab = 2*nt-1;
        vt = ealloc1float(nttab);
	vnt = ealloc1float(nttab);
	etat = ealloc1float(nttab);
	a1111 = ealloc1float(nttab);
	a3333 = ealloc1float(nttab);
	a1313 = ealloc1float(nttab);
	a1133 = ealloc1float(nttab);
	vrms = ealloc1float(nttab);

        /* determine NMO velocity function vnm(t) */
        if (!getparstring("vnfile",&vnfile)) {
                ntdmo = countparval("tdmo");
                if (ntdmo==0) ntdmo = 1;
                tdmo = ealloc1float(ntdmo);
                if (!getparfloat("tdmo",tdmo)) tdmo[0] = 0.0;
                nvndmo = countparval("vndmo");
                if (nvndmo==0) nvndmo = 1;
                if (nvndmo!=ntdmo) err("number of tdmo and vndmo must be equal");
                vndmo = ealloc1float(nvndmo);
                if (!getparfloat("vndmo",vndmo)) vndmo[0] = 1500.0;
                for (itdmo=1; itdmo<ntdmo; ++itdmo)
                        if (tdmo[itdmo]<=tdmo[itdmo-1])
                                err("tdmo must increase monotonically");
                for (it=0,t=0.0; it<nt; ++it,t+=dt)
                        intlin(ntdmo,tdmo,vndmo,vndmo[0],vndmo[ntdmo-1],
                                1,&t,&vnt[it]);
        } else {
                if (fread(vnt,sizeof(float),nt,fopen(vnfile,"r"))!=nt)
                        err("cannot read %d velocities from file %s",nt,vnfile);
        }

	/* constant extrapolation of interval velocities */
	for (it=nt; it<nttab; ++it)
		vnt[it] = vnt[nt-1];

	/* determine eta function eta(t) */
        if (!getparstring("etafile",&etafile)) {
                ntdmo = countparval("tdmo");
                if (ntdmo==0) ntdmo = 1;
                tdmo = ealloc1float(ntdmo);
                if (!getparfloat("tdmo",tdmo)) tdmo[0] = 0.0;
                netadmo = countparval("etadmo");
                if (netadmo==0) netadmo = 1;
                if (netadmo!=ntdmo) err("number of tdmo and etadmo must be equal");
                etadmo = ealloc1float(netadmo);
                if (!getparfloat("etadmo",etadmo)) etadmo[0] = 0.0;
                for (itdmo=1; itdmo<ntdmo; ++itdmo)
                        if (tdmo[itdmo]<=tdmo[itdmo-1])
                                err("tdmo must increase monotonically");
                for (it=0,t=0.0; it<nt; ++it,t+=dt)
                        intlin(ntdmo,tdmo,etadmo,etadmo[0],etadmo[ntdmo-1],
                                1,&t,&etat[it]);
        } else {
                if (fread(etat,sizeof(float),nt,fopen(etafile,"r"))!=nt)
                        err("cannot read %d velocities from file %s",nt,etafile);
        }

	/* constant extrapolation of interval velocities */
	for (it=nt; it<nttab; ++it)
		etat[it] = etat[nt-1];

	/* determine velocity function v(t) */
        if (!getparstring("vfile",&vfile)) {
                ntdmo = countparval("tdmo");
                if (ntdmo==0) ntdmo = 1;
                tdmo = ealloc1float(ntdmo);
                if (!getparfloat("tdmo",tdmo)) tdmo[0] = 0.0;
                nvdmo = countparval("vdmo");
                if (nvdmo!=ntdmo){ 
			for (it=0; it<nt; ++it)
				vt[it] = vnt[it];
		} else {
                	vdmo = ealloc1float(nvdmo);
                	if (!getparfloat("vdmo",vdmo)) vdmo[0] = vnt[0];
                	for (itdmo=1; itdmo<ntdmo; ++itdmo)
                        	if (tdmo[itdmo]<=tdmo[itdmo-1])
                                	err("tdmo must increase monotonically");
                	for (it=0,t=0.0; it<nt; ++it,t+=dt)
                        	intlin(ntdmo,tdmo,vdmo,vdmo[0],vdmo[ntdmo-1],
                                	1,&t,&vt[it]);
		}
        } else {
                if (fread(vt,sizeof(float),nt,fopen(vfile,"r"))!=nt)
                        err("cannot read %d velocities from file %s",nt,vfile);
        }

	/* constant extrapolation of interval velocities */
	for (it=nt; it<nttab; ++it)
		vt[it] = vt[nt-1];

	/* compute vrms */
	vrms[0] = vnt[0];
	for (it=1,sum=0.0; it<nttab; ++it) {
		sum += vnt[it-1]*vnt[it-1];
			vrms[it] = sqrt(sum/it);
	}
        
	/*convert eta and NMO velocity values to elastic coeffecients*/
	for(it=0; it<nttab; ++it){
		vnt2   =vnt[it]*vnt[it];
		vt2    =vt[it]*vt[it];
		a1111[it]=vnt2*(1+2*etat[it]);
		a3333[it]=vt2;
		a1313[it]=0.25*vt2;
		delta= 0.5*(vnt2/vt2-1);
		a1133[it]=sqrt(2*delta*a3333[it]*
			(a3333[it]-a1313[it])+(a3333[it]-
			a1313[it])*(a3333[it]-a1313[it]))-a1313[it];
	}



	/* compute half-velocities */
        for (it=0; it<nttab; ++it){
                a1111[it] *= 0.25;
		a3333[it] *= 0.25;
		a1313[it] *= 0.25;
		a1133[it] *= 0.25;
		vt[it]    *= 0.5;
	}

	/* compute extreme values */
	ft = 0.0;
	lt = ft+(nt-1)*dt;
	lttab = ft+(nttab-1)*dt;
	fp = 0.0;
	lp = 1.0/sqrt(a1111[0]);

	/* compute p sampling for ray tables */
	dptab = lp/(nptab-1);

	/* allocate space for ray tables */
	tau = ealloc2float(nttab,nptab);
	x = ealloc2float(nttab,nptab);
	a = ealloc2float(nttab,nptab);
	s = ealloc2float(nttab,nptab);

	/* compute ray tables tau(p,t), x(p,t), a(p,t) */
	if (verbose) fprintf(stderr,"Computing ray tables...\n");
	for (ip=0,p=fp; ip<nptab; ++ip,p+=dptab)
		rayvt(p,nttab,dt,a1111,a3333,a1313,
			a1133,tau[ip],x[ip],a[ip],s[ip]);
		
	if (verbose) fprintf(stderr,"Ray tables ready...\n");
	free1float(vt);
	free1float(vnt);
	free1float(etat);
	free1float(a1111);
	free1float(a3333);
	free1float(a1133);
	free1float(a1313);
	

	/* determine number of cdps to process */
	nx = cdpmax-cdpmin+1;
	
	/* allocate and zero common-offset gather q(t,x) */
	nxpad = 0.5*ABS(offset/dx);
	nxfft = npfar(nx+nxpad);
	q = ealloc2float(nt,nxfft+2);
	for (ix=0; ix<nxfft; ++ix)
		for (it=0; it<nt; ++it)
			q[ix][it] = 0.0;

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

	/* 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(x,a,tau,s,nttab,dt,ft,nptab,dptab,fp,lp,
				nx,dx,oldoffset,nt,vrms,fmax,smute,speed,
				verbose,q);
			
			/* add dmo-corrected traces to mix */
			for (ix=0; ix<nx; ++ix)
				for (it=0; it<nt; ++it)
					qmix[ix][it] += q[ix][it];
			
			/* count offsets in mix */
			++noff;
							
			/* free space for old common-offset gather */
			free2float(q);
			
			/* if beginning a new offset */
			if (offset!=oldoffset) {
				
				/* allocate space for new offset */
				nxpad = 0.5*ABS(offset/dx);
				nxfft = npfar(nx+nxpad);
				q = ealloc2float(nt,nxfft+2);
				for (ix=0; ix<nxfft; ++ix)
					for (it=0; it<nt; ++it)
						q[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 *) qmix[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)
					qmix[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;