sumiggbzo.c

su 的源代码库

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

/* SUMIGGBZO: $Revision: 1.10 $ ; $Date: 2006/09/28 19:24:23 $		*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SUMIGGBZO - MIGration via Gaussian Beams of Zero-Offset SU data	",
"									",
" sumiggbzo <infile >outfile vfile=  nz= [optional parameters]		",
"									",
" Required Parameters:							",
" vfile=                 name of file containing v(z,x)			",
" nz=                    number of depth samples			",
"									",
" Optional Parameters:							",
" dt=from header		time sampling interval			",
" dx=from header(d2) or 1.0	spatial sampling interval		",
" dz=1.0                 depth sampling interval			",
" fmin=0.025/dt          minimum frequency				",
" fmax=10*fmin           maximum frequency				",
" amin=-amax             minimum emergence angle; must be > -90 degrees	",
" amax=60                maximum emergence angle; must be < 90 degrees	",
" bwh=0.5*vavg/fmin      beam half-width; vavg denotes average velocity	",
" verbose=0		 =0 silent; =1 chatty				",
"									",
" Note: spatial units of v(z,x) must be the same as those on dx.	",
" v(z,x) is represented numerically in C-style binary floats v[x][z],	",
" where the depth direction is the fast direction in the data. Such	",
" models can be created with unif2 or makevel.				",
"									",
NULL};

/* Credits:
 *
 * CWP: Dave Hale (algorithm), Jack K. Cohen, and John Stockwell
 * (reformatting for SU)
 */

/**************** end self doc ***********************************/

/* Ray types */
/* one step along ray */
typedef struct RayStepStruct {
	float t;		/* time */
	float a;		/* angle */
	float x,z;		/* x,z coordinates */
	float q1,p1,q2,p2;	/* Cerveny's dynamic ray tracing solution */
	int kmah;		/* KMAH index */
	float c,s;		/* cos(angle) and sin(angle) */
	float v,dvdx,dvdz;	/* velocity and its derivatives */
} RayStep;

/* one ray */
typedef struct RayStruct {
	int nrs;		/* number of ray steps */
	RayStep *rs;		/* array[nrs] of ray steps */
	int nc;			/* number of circles */
	int ic;			/* index of circle containing nearest step */
	void *c;		/* array[nc] of circles */
} Ray;

/* Ray functions */
Ray *makeRay (float x0, float z0, float a0, int nt, float dt, float ft,
	int nx, float dx, float fx, int nz, float dz, float fz, float **v);
void freeRay (Ray *ray);
int nearestRayStep (Ray *ray, float x, float z);

/* Velocity functions */
void* vel2Alloc (int nx, float dx, float fx,
	int nz, float dz, float fz, float **v);
void vel2Free (void *vel2);
void vel2Interp (void *vel2, float x, float z,
	float *v, float *vx, float *vz, float *vxx, float *vxz, float *vzz);

/* Beam functions */
void formBeams (float bwh, float dxb, float fmin,
	int nt, float dt, float ft, 
	int nx, float dx, float fx, float **f,
	int ntau, float dtau, float ftau, 
	int npx, float dpx, float fpx, float **g);
void accBeam (Ray *ray, float fmin, float lmin,
	int nt, float dt, float ft, float *f,
	int nx, float dx, float fx, int nz, float dz, float fz, float **g);

/* functions defined and used internally */
static void miggbzo (float bwh, float fmin, float fmax, float amin, float amax,
	int nt, float dt, int nx, float dx, int nz, float dz, 
	float **f, float **v, float **g, int verbose);

segy tr;

int
main(int argc, char **argv)
{
	int nx;		/* number of horizontal samples (traces) in input */
	int nz;		/* number of vertical samples in output	          */	
	int nt;		/* number of time samples in input		  */
	int ix;		/* counter in x					  */
	int iz;		/* counter in z					  */

	float dx;	/* trace spacing in input			  */
	float dz;	/* vertical sampling interval in output		  */
	float dt;	/* time sampling interval in input	          */

	float fmin;	/* minimum frequency in migration		  */
	float fmax;     /* maximum frequency in migration 		  */
	float amin;	/* minimum ray angle				  */
	float amax;     /* maximum ray angle 			          */
	float vavg;	/* average velocity in input vfile		  */

	float bwh;	/* gaussian beam half-width			  */

	float **v=NULL;	/* pointer to background velocities		  */
	float **f=NULL; /* pointer to input data                          */
	float **g=NULL; /* pointer to migrated data 			  */

	char *vfile="";		/* velocity filename		*/
	FILE *vfp;		/* velocity file pointer	*/
	FILE *tracefp;		/* temp file to hold traces	*/
	int verbose;		/* verbose flag			*/

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


	/* get info from first trace */ 
	if (!gettr(&tr))  err("can't get first trace");
	nt = tr.ns;

	/* get required parameters */
	MUSTGETPARSTRING("vfile", &vfile);
	MUSTGETPARINT("nz", &nz);
	MUSTGETPARFLOAT("dz", &dz);
	if (!getparint("verbose",&verbose))	verbose = 0;


	/* let user give dt and/or dx from command line */
	if (!getparfloat("dt", &dt)) {
		if (tr.dt) { /* is dt field set? */
			dt = (float) tr.dt / 1000000.0;
		} else { /* dt not set, assume 4 ms */
			dt = 0.004;
			if (verbose) warn("tr.dt not set, assuming dt=0.004");
		}
	}
	if (!getparfloat("dx",&dx)) {
		if (tr.d2) { /* is d2 field set? */
			dx = tr.d2;
		} else {
			dx = 1.0;
			if (verbose) warn("tr.d2 not set, assuming d2=1.0");
		}
	}

	
	/* get optional parameters */
	if (!getparfloat("dz",&dz)) dz = 1.0;
	if (!getparfloat("fmin",&fmin)) fmin = 0.025/dt;
	if (!getparfloat("fmax",&fmax)) fmax = 10.0*fmin;
	if (!getparfloat("amax",&amax)) amax = 60.0;
	if (!getparfloat("amin",&amin)) amin = -amax;

	/* store traces in tmpfile while getting a count */
	tracefp = etmpfile();
	nx = 0;
	do { 
		++nx;
		efwrite(tr.data, FSIZE, nt, tracefp);
	} while (gettr(&tr));


	/* allocate workspace */
	f = ealloc2float(nt,nx);
	v = ealloc2float(nz,nx);
	g = ealloc2float(nz,nx);

	/* load traces into the zero-offset array and close tmpfile */
	rewind(tracefp);
	efread(*f, FSIZE, nt*nx, tracefp);
	efclose(tracefp);

	/* read and halve velocities; determine average */
	vfp = efopen(vfile,"r");
	if (efread(*v, FSIZE, nz*nx, vfp)!=nz*nx)
		err("error reading vfile=%s!\n",vfile);
	for (ix=0,vavg=0.0; ix<nx; ++ix) {
		for (iz=0; iz<nz; ++iz) {
			v[ix][iz] *= 0.5;
			vavg += v[ix][iz];
		}
	}
	vavg /= nx*nz;

	/* get beam half-width */
	if (!getparfloat("bwh",&bwh)) bwh = vavg/fmin;
	
	if (verbose) warn("bhw=%f",bwh);
	
	/* migrate */
	miggbzo(bwh,fmin,fmax,amin,amax,nt,dt,nx,dx,nz,dz,f,v,g,verbose);
	
	/* set header fields and write output */
	tr.ns = nz;
	tr.trid = TRID_DEPTH;
	tr.d1 = dz;
	tr.d2 = dx;
	
	for (ix=0; ix<nx; ++ix) {
		tr.tracl = ix + 1;
		tr.tracr = ix + 1;
		for (iz=0; iz<nz; ++iz) {
			tr.data[iz] = g[ix][iz];
		}
		puttr(&tr);
	}

	/* free workspace */
	free2float(f);
	free2float(v);
	free2float(g);
	
	return(CWP_Exit());
}

static void miggbzo (float bwh, float fmin, float fmax, float amin, float amax,
	int nt, float dt, int nx, float dx, int nz, float dz, 
	float **f, float **v, float **g, int verbose)
/*****************************************************************************
Migrate zero-offset data via accumulation of Gaussian beams.
******************************************************************************
Input:
bwh		horizontal beam half-width at surface z=0
fmin		minimum frequency (cycles per unit time)
fmax		maximum frequency (cycles per unit time)
amin		minimum emergence angle at surface z=0 (degrees)
amax		maximum emergence angle at surface z=0 (degrees)
nt		number of time samples
dt		time sampling interval (first time assumed to be zero)
nx		number of x samples
dx		x sampling interval
nz		number of z samples
dz		z sampling interval
f		array[nx][nt] containing zero-offset data f(t,x)
v		array[nx][nz] containing half-velocities v(x,z)

Output:
g		array[nx][nz] containing migrated image g(x,z)
*****************************************************************************/
{
	int nxb=0,npx=0,ntau=0,ipx,ix,ixb,ixlo,ixhi,nxw,iz;
	float ft,fx,fz,xwh,dxb=0,fxb,xb,vmin,dpx,fpx,px,
		taupad,dtau,ftau,fxw,pxmin,pxmax,
		a0,x0,z0,bwhc,**b;
	Ray *ray;

	/* first t, x, and z assumed to be zero */
	ft = fx = fz = 0.0;

	/* convert minimum and maximum angles to radians */
	amin *= PI/180.0;
	amax *= PI/180.0;
	if (amin>amax) {
		float atemp=amin;
		amin = amax;
		amax = atemp;
	}
	
	/* window half-width */
	xwh = 3.0*bwh;

	/* beam center sampling */
	dxb = NINT((bwh*sqrt(2.0*fmin/fmax))/dx)*dx;
	nxb = 1 + (nx-1)*dx/dxb;
	fxb = fx+0.5*((nx-1)*dx-(nxb-1)*dxb);

	if (verbose) warn("nxb=%d dxb=%d fxb=%d ",nxb,dxb,fxb);

	/* minimum velocity at surface z=0 */
	for (ix=1,vmin=v[0][0]; ix<nx; ++ix)
		if (v[ix][0]<vmin) vmin = v[ix][0];
	if (verbose) warn("vmin=%f fmin=%f fmax=%f",vmin,fmin,fmax);

	/* beam sampling */
	pxmin = sin(amin)/vmin;
	pxmax = sin(amax)/vmin;
	dpx = 1.0/(2.0*xwh*sqrt(fmin*fmax));
	npx = 1+(pxmax-pxmin)/dpx;
	fpx = pxmin+0.5*(pxmax-pxmin-(npx-1)*dpx);
	taupad = MAX(ABS(pxmin),ABS(pxmax))*xwh;
	taupad = NINT(taupad/dt)*dt;
	ftau = ft-taupad;
	dtau = dt;
	ntau = nt+2.0*taupad/dtau;

	if (verbose) warn("npx=%d dpx=%g fpx=%g",npx,dpx,fpx);
	if (verbose) warn("pxmin=%f pxmax=%f",pxmin,pxmax);
	if (verbose) warn("amin=%f amax=%f",amin,amax);

	/* zero migrated image */
	for (ix=0; ix<nx; ++ix)
		for (iz=0; iz<nz; ++iz)
			g[ix][iz] = 0.0;

	/* loop over beam centers */
	for (ixb=0,xb=fxb; ixb<nxb; ++ixb,xb+=dxb) {

		/* horizontal window */
		ix = NINT((xb-fx)/dx);
		ixlo = MAX(ix-NINT(xwh/dx),0);
		ixhi = MIN(ix+NINT(xwh/dx),nx-1);
		nxw = 1+ixhi-ixlo;
		fxw = fx+(ixlo-ix)*dx;

		if (verbose) warn("ixb/nxb = %d/%d  ix = %d",ixb,nxb,ix);

		/* allocate space for beams */
		b = alloc2float(ntau,npx);

		/* form beams at surface */
		formBeams(bwh,dxb,fmin,
			nt,dt,ft,nxw,dx,fxw,&f[ixlo],
			ntau,dtau,ftau,npx,dpx,fpx,b);
		
		/* loop over beams */
		for (ipx=0,px=fpx; ipx<npx; ++ipx,px+=dpx) {

			/* sine of emergence angle; skip if out of bounds */
			if (px*v[ix][0]>sin(amax)+0.01) continue;
			if (px*v[ix][0]<sin(amin)-0.01) continue;

			/* emergence angle and location */
			a0 = -asin(px*v[ix][0]);
			x0 = fx+ix*dx;
			z0 = fz;

			/* beam half-width adjusted for cosine of angle */
			bwhc = bwh*cos(a0);

			/* trace ray */
			ray = makeRay(x0,z0,a0,nt,dt,ft,nx,dx,fx,nz,dz,fz,v);

			/* accumulate contribution of beam in migrated image */
			accBeam(ray,fmin,bwhc,
				ntau,dtau,ftau,b[ipx],
				nx,dx,fx,nz,dz,fz,g);
			/*
			fwrite(g[0],sizeof(float),nx*nz,stdout);
			*/
			
			/* free ray */
			freeRay(ray);
		}
		/*
		fwrite(g[0],sizeof(float),nx*nz,stdout);
		*/

		/* free space for beams */
		free2float(b);
	}
}

/* circle for efficiently finding nearest ray step */
typedef struct CircleStruct {
	int irsf;               /* index of first raystep in circle */
	int irsl;               /* index of last raystep in circle */
	float x;                /* x coordinate of center of circle */
	float z;                /* z coordinate of center of circle */
	float r;                /* radius of circle */
} Circle;

/* functions defined and used internally */
Circle *makeCircles (int nc, int nrs, RayStep *rs);

Ray *makeRay (float x0, float z0, float a0, int nt, float dt, float ft,
	int nx, float dx, float fx, int nz, float dz, float fz, float **vxz)
/*****************************************************************************
Trace a ray for uniformly sampled v(x,z).
******************************************************************************
Input:
x0		x coordinate of takeoff point
z0		z coordinate of takeoff point
a0		takeoff angle (radians)
nt		number of time samples
dt		time sampling interval
ft		first time sample
nx		number of x samples
dx		x sampling interval
fx		first x sample
nz		number of z samples
dz		z sampling interval
fz		first z sample
vxz		array[nx][nz] of uniformly sampled velocities v(x,z)

Returned:	pointer to ray parameters sampled at discrete ray steps
******************************************************************************
Notes:
The ray ends when it runs out of time (after nt steps) or with the first 
step that is out of the (x,z) bounds of the velocity function v(x,z).
*****************************************************************************/
{
	int it,kmah;
	float t,x,z,a,c,s,p1,q1,p2,q2,
		lx,lz,cc,ss,
		v,dvdx,dvdz,ddvdxdx,ddvdxdz,ddvdzdz,
		vv,ddvdndn;
	Ray *ray;
	RayStep *rs;
	void *vel2;

	/* allocate and initialize velocity interpolator */
	vel2 = vel2Alloc(nx,dx,fx,nz,dz,fz,vxz);

	/* last x and z in velocity model */
	lx = fx+(nx-1)*dx;
	lz = fz+(nz-1)*dz;

	/* ensure takeoff point is within model */
	if (x0<fx || x0>lx || z0<fz || z0>lz) return NULL;

	/* allocate space for ray and raysteps */
	ray = (Ray*)alloc1(1,sizeof(Ray));
	rs = (RayStep*)alloc1(nt,sizeof(RayStep));

	/* cosine and sine of takeoff angle */
	c = cos(a0);
	s = sin(a0);
	cc = c*c;
	ss = s*s;

	/* velocity and derivatives at takeoff point */
	vel2Interp(vel2,x0,z0,&v,&dvdx,&dvdz,&ddvdxdx,&ddvdxdz,&ddvdzdz);
	ddvdndn = ddvdxdx*cc-2.0*ddvdxdz*s*c+ddvdzdz*ss;
	vv = v*v;

	/* first ray step */
	rs[0].t = t = ft;
	rs[0].a = a = a0;
	rs[0].x = x = x0;
	rs[0].z = z = z0;
	rs[0].q1 = q1 = 1.0;
	rs[0].p1 = p1 = 0.0;
	rs[0].q2 = q2 = 0.0;
	rs[0].p2 = p2 = 1.0;
	rs[0].kmah = kmah = 0;
	rs[0].c = c;
	rs[0].s = s;
	rs[0].v = v;
	rs[0].dvdx = dvdx;
	rs[0].dvdz = dvdz;

	/* loop over time steps */
	for (it=1; it<nt; ++it) {

		/* variables used for Runge-Kutta integration */
		float h=dt,hhalf=dt/2.0,hsixth=dt/6.0,
			q2old,xt,zt,at,p1t,q1t,p2t,q2t,
			dx,dz,da,dp1,dq1,dp2,dq2,
			dxt,dzt,dat,dp1t,dq1t,dp2t,dq2t,
			dxm,dzm,dam,dp1m,dq1m,dp2m,dq2m;

		/* if ray is out of bounds, break */
		if (x<fx || x>lx || z<fz || z>lz) break;

		/* remember old q2 */
		q2old = q2;
		
		/* step 1 of 4th-order Runge-Kutta */
		dx =  v*s;
		dz =  v*c;
		da =  dvdz*s-dvdx*c;
		dp1 = -ddvdndn*q1/v;
		dq1 = vv*p1;
		dp2 = -ddvdndn*q2/v;
		dq2 =  vv*p2;
		xt = x+hhalf*dx;