sumiggbzoan.c

su 的源代码库

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

/* SUMIGGBZOAN: $Revision: 1.2 $ ; $Date: 2006/04/14 20:49:11 $           */

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

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SUMIGGBZOAN - MIGration via Gaussian beams ANisotropic media (P-wave)	",
"									",
" sumiggbzoan <infile >outfile vfile= nt= nx= nz= [optional parameters]	",
"									",
" Required Parameters:							",
" a3333file=		name of file containing a3333(x,z)		",
" nx=                    number of inline samples (traces)		",
" nz=                    number of depth samples			",
" 									",
" Optional Parameters:							",
" dt=tr.dt               time sampling interval				",
" dx=tr.d2               inline sampling interval (trace spacing)	",
" 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	",
" 									",
" Files for general anisotropic parameters confined to a vertical plane:",
" a1111file=		name of file containing a1111(x,z)		",
" a1133file=          	name of file containing a1133(x,z)		",
" a1313file=          	name of file containing a1313(x,z)		",
" a1113file=          	name of file containing a1113(x,z)		",
" a3313file=          	name of file containing a3313(x,z)		",
"									",
" For transversely isotropic media Thomsen's parameters could be used:	",
" deltafile=		name of file containing delta(x,z)		",
" epsilonfile=		name of file containing epsilon(x,z)		",
" a1313file=          	name of file containing a1313(x,z)		",
"									",
" if anisotropy parameters are not given the program considers		", 
" the medium to be isotropic.						",
"									",
NULL};

/* Credits:
 *	CWP: Tariq Alkhalifah,  based on MIGGBZO by Dave Hale
 *      CWP: repackaged as an SU program by John Stockwell, April 2006
 *      
 *   Technical Reference:
 *
 *      Alkhailfah, T., 1993, Gaussian beam migration for
 *      anisotropic media: submitted to Geophysics.
 *
 *	Cerveny, V., 1972, Seismic rays and ray intensities 
 *	in inhomogeneous anisotropic media: 
 *	Geophys. J. R. Astr. Soc., 29, 1--13.
 *
 *	Hale, D., 1992, Migration by the Kirchhoff, 
 *	slant stack, and Gaussian beam methods:
 *      CWP,1992 Report 121, Colorado School of Mines.
 *
 *	Hale, D., 1992, Computational Aspects of Gaussian
 *      Beam migration:
 *     	CWP,1992 Report 121, Colorado School of Mines.
 *
 *
 */

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


/* functions defined and used internally */
/* one step along ray */
typedef struct RayStepStruct {
	float t;		/* time */
	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;

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

/* 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);
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);

static void miggbzoaniso(float bwh, float fmin, float fmax, 
				float amin, float amax, int nt, 
				float dt, int nx, float dx, int nz, 
				float dz, float **a1111, float **a3333,
				float **a1133, float **a1313, float **a1113,
				float **a3313, float **f, float **g);

segy tr;

int
main(int argc, char **argv)
{
	int nx;		/* number of traces in input		*/
	int nz; 	/* number of depth values in output	*/
	int nt;		/* number of time samples in input	*/
	int ix,iz;	/* counters				*/
	int ia1111=1,ia1133=1,ia1313=1,ia1113=1,ia3313=1,
		idelta=1,iepsilon=1;
	float dx,dz,dt,fmin,fmax,amin,amax,vavg,vavg1,vavg3,bwh,**a1111,**a3333,
		**a1133,**a1313,**a1113,**a3313,**delta,**epsilon,**f,**g;
	char *a1111file="",*a3333file="",*a1133file="",*a1313file="",
		*a1113file="",*a3313file="",*deltafile="",*epsilonfile="";
	FILE *a1111fp,*a3333fp,*a1133fp,
		*a1313fp,*a1113fp,*a3313fp,*deltafp,*epsilonfp;
	FILE *tracefp;	/* temporary file to hold traces */

	/* 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("a3333file",&a3333file);
	MUSTGETPARINT("nt",&nt);
	MUSTGETPARINT("nx",&nx);
	MUSTGETPARINT("nz",&nz);

	/* 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;
                        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;
                        warn("tr.d2 not set, assuming d2=1.0");
                }
	}
	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 = alloc2float(nt,nx);
	a1111 = alloc2float(nz,nx);
	a3333 = alloc2float(nz,nx);
	a1133 = alloc2float(nz,nx);
	a1313 = alloc2float(nz,nx);
	a1113 = alloc2float(nz,nx);
	a3313 = alloc2float(nz,nx);
	delta = alloc2float(nz,nx);
	epsilon = alloc2float(nz,nx);
	g = alloc2float(nz,nx);

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

	/*anisotropic parameters*/
	if (!getparstring("a1111file",&a1111file)) 
			ia1111=0;
	if (!getparstring("a1133file",&a1133file)) 
			ia1133=0;
	if (!getparstring("a1313file",&a1313file)) 
			ia1313=0;
	if (!getparstring("a1113file",&a1113file)) 
			ia1113=0;
	if (!getparstring("a3313file",&a3313file)) 
			ia3313=0;
	if (!getparstring("deltafile",&deltafile)) 
			idelta=0;
	if (!getparstring("epsilonfile",&epsilonfile)) 
			iepsilon=0;

	
	/* read and halve required velocities*/
	if ((a3333fp=fopen(a3333file,"r"))==NULL)
		err("error opening a3333file=%s!\n",a3333file);
	if (fread(a3333[0],sizeof(float),nz*nx,a3333fp)!=nz*nx)
		err("error reading a3333file=%s!\n",a3333file);

	/* read and halve velocities*/
	if(ia1111){
		if ((a1111fp=fopen(a1111file,"r"))==NULL)
			err("error opening a1111file=%s!\n",a1111file);
		if (fread(a1111[0],sizeof(float),nz*nx,a1111fp)!=nz*nx)
			err("error reading a1111file=%s!\n",a1111file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a1111[ix][iz] = a3333[ix][iz];
	}		
	
	if(ia1133){
		if ((a1133fp=fopen(a1133file,"r"))==NULL)
			err("error opening a1133file=%s!\n",a1133file);
		if (fread(a1133[0],sizeof(float),nz*nx,a1133fp)!=nz*nx)
			err("error reading a1133file=%s!\n",a1133file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a1133[ix][iz] = .4*a3333[ix][iz];
	}

	if(ia1313){
		if ((a1313fp=fopen(a1313file,"r"))==NULL)
			err("error opening a1313file=%s!\n",a1313file);
		if (fread(a1313[0],sizeof(float),nz*nx,a1313fp)!=nz*nx)
			err("error reading a1313file=%s!\n",a1313file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a1313[ix][iz] = .3*a3333[ix][iz];
	}

	if(ia1113){
		if ((a1113fp=fopen(a1113file,"r"))==NULL)
			err("error opening a1113file=%s!\n",a1113file);
		if (fread(a1113[0],sizeof(float),nz*nx,a1113fp)!=nz*nx)
			err("error reading a1113file=%s!\n",a1113file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a1113[ix][iz] = 0;
	}

	if(ia3313){
		if ((a3313fp=fopen(a3313file,"r"))==NULL)
			err("error opening a3313file=%s!\n",a3313file);
		if (fread(a3313[0],sizeof(float),nz*nx,a3313fp)!=nz*nx)
			err("error reading a3313file=%s!\n",a3313file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a3313[ix][iz] = 0;
	}

	/* if specified read Thomsen's parameters*/
	if(idelta==1 || iepsilon==1) {
		if(idelta){
			if ((deltafp=fopen(deltafile,"r"))==NULL)
				err("error opening deltafile=%s!\n",deltafile);
			if (fread(delta[0],sizeof(float),nz*nx,deltafp)!=nz*nx)
				err("error reading deltafile=%s!\n",deltafile);
		}else{
			for (ix=0; ix<nx; ++ix)
				for (iz=0; iz<nz; ++iz)
					delta[ix][iz] = 0;
		}
		if(iepsilon){
			if ((epsilonfp=fopen(epsilonfile,"r"))==NULL)
				err("error opening epsilonfile=%s!\n",epsilonfile);
			if (fread(epsilon[0],sizeof(float),nz*nx,epsilonfp)!=nz*nx)
				err("error reading epsilonfile=%s!\n",epsilonfile);
		}else{
			for (ix=0; ix<nx; ++ix)
				for (iz=0; iz<nz; ++iz)
					epsilon[ix][iz] = 0;
		}
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz){
				a1111[ix][iz] = a3333[ix][iz]+2*epsilon[ix][iz]
					*a3333[ix][iz];
				a1133[ix][iz] = sqrt(2*delta[ix][iz]*a3333[ix][iz]*
					(a3333[ix][iz]-a1313[ix][iz])+(a3333[ix][iz]-
					a1313[ix][iz])*(a3333[ix][iz]-a1313[ix][iz]))
					-a1313[ix][iz];
			}
	}
	
	/* free workspace */
	free2float(delta);
	free2float(epsilon);
		

	/*determine average velocity*/
	vavg3=0.0; 
	for (ix=0,vavg1=0.0; ix<nx; ++ix) {
		for (iz=0; iz<nz; ++iz) {
			a1111[ix][iz] *= 0.25;
			vavg1 += sqrt(a1111[ix][iz]);
			a3333[ix][iz] *= 0.25;
			vavg3 += sqrt(a3333[ix][iz]);
			a1133[ix][iz] *= 0.25;
			a1313[ix][iz] *= 0.25;
			a1113[ix][iz] *= 0.25;
			a3313[ix][iz] *= 0.25;
		}
	}
	vavg = (vavg1+vavg3)/2;
	vavg /= nx*nz;

	/* get beam half-width */
	if (!getparfloat("bwh",&bwh)) bwh = vavg/fmin;
	
	/* migrate */
	miggbzoaniso(bwh,fmin,fmax,amin,amax,nt,dt,nx,dx,nz,dz,a1111,a3333,
		a1133,a1313,a1113,a3313,f,g);

        /* 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;
               for (iz=0; iz<nz; ++iz) {
                        tr.data[iz] = g[ix][iz];
                }
                puttr(&tr);
        }

	/* free workspace */
	free2float(f);
	free2float(a1111);
	free2float(a3333);
	free2float(a1133);
	free2float(a1313);
	free2float(a1113);
	free2float(a3313);
	free2float(g);
	
	return(CWP_Exit());
}


static void miggbzoaniso (float bwh, float fmin, float fmax, float amin, float amax,
	int nt, float dt, int nx, float dx, int nz, float dz, 
	float **a1111, float **a3333, float **a1133, float **a1313,
	float **a1113, float **a3313,
	float **f, float **g)
/*****************************************************************************
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,npx,ntau,ipx,ix,ixb,ixlo,ixhi,nxw,iz,ixm,i;
	float ft,fx,fz,xwh,dxb,fxb,xb,vmin,dpx,fpx,px,
		taupad,dtau,ftau,fxw,pxmin,pxmax,