suinvco3d.c

su 的源代码库

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

/* SUINVCO3D: $Revision: 1.1 $ ; $Date: 2001/03/14 18:14:10 $	*/

#include "su.h"
#include "segy.h"
#include <string.h>
#include <stdio.h>
#include <stdlib.h>

/*********************** self documentation **********************/
char *sdoc[] = {
" SUINVCO3D - Seismic INVersion of Common Offset data with V(X,Y,Z) velocity	",
"	     function in 3D							",
" 										",
"     suinvco3d <infile >outfile [optional parameters] 				",
"										",
" Required Parameters:								",
" vfile=		  file containing velocity array v[nvy][nvx][nvz]	",
" nzv=		   number of z samples (1st dimension) in velocity		",
" nxm=			number of midpoints of input traces			",
" nym=			number of lines 					",
" geo_type=		geometry type						",
"			1 ---- general velocity distribution v(x,y,z)		",
"			2 ---- v(x,z) medium					",
"			3 ---- v(z) medium					",
" com_type=		computation type, determines what tables are needed	",
"			1 ---- only needs traveltime,	   weight=1.0		",
"			2 ---- traveltime, propagation angles,  weight=ctheta	",
"			3 ---- traveltime, angle and amplitude,			",
"						  weight=det/as/ag/(1+ctheta)	",
" nzt=		   number of z samples (1st dimension) in traveltime		",
" nxt=		   number of x samples (2nd dimension) in traveltime		",
" nyt=		   number of y samples (3rd dimension) in traveltime		",
" tfile		  file containing traveltime array t[nyt][nxt][nzt]		",
" ampfile		file containing amplitude array amp[nyt][nxt][nzt]	",
" d21file		file containing Beylkin determinant component array	",
" d22file		file containing Beylkin determinant component array	",
" d23file		file containing Beylkin determinant component array	",
" d31file		file containing Beylkin determinant component array	",
" d32file		file containing Beylkin determinant component array	",
" d33file		file containing Beylkin determinant component array	",
" a1file		 file containing ray propagation angle (polar) array	",
" b1file		 file containing ray propagation angle (azimuth) array	",
"										",
" Optional Parameters:								",
" dt= or from header (dt) 	time sampling interval of input data		",
" offs= or from header (offset) 	source-receiver offset 			",
" dxm= or from header (d2) 	x sampling interval of midpoints 		",
" fxm=0		  first midpoint in input trace					",
" dym=50.0		y sampling interval of midpoints 			",
" fym=0		  y-coordinate of first midpoint in input trace			",
" nxv=		   number of x samples (2nd dimension) in velocity		",
" nyv=		   number of y samples (3rd dimension) in velocity		",
" dxv=50.0		x sampling interval of velocity				",
" fxv=0.0		first x sample of velocity				",
" dyv=50.0		y sampling interval of velocity				",
" fyv=0.0		first y sample of velocity				",
" dzv=50.0		z sampling interval of velocity				",
" fzv=0.0		first z sample of velocity				",
" nxb=nx/2		band centered at midpoints (see note)			",
" fxo=0.0		x-coordinate of first output trace 			",
" dxo=15.0		horizontal spacing of output trace 			",
" nxo=101		number of output traces 				",	
" fyo=0.0		y-coordinate of first output trace			",
" dyo=15.0		y-coordinate spacing of output trace			",
" nyo=101		number of output traces in y-direction			",
" fzo=0.0		z-coordinate of first point in output trace 		",
" dzo=15.0		vertical spacing of output trace 			",
" nzo=101		number of points in output trace			",	
" dxt=100.0		x-coordinate spacing of input tables(traveltime, etc)	",
" dyt=100.0		y-coordinate spacing of input tables(traveltime, etc)	",
" dzt=100.0		z-coordinate spacing of input tables(traveltime, etc)	",
" xt0=0.0		x-coordinate of first input tables			",
" xt1=0.0		x-coordinate of last input tables			",
" yt0=0.0		y-coordinate of first input tables		 	",
" yt1=0.0		y-coordinate of last input tables			",
" fmax=0.25/dt		Maximum frequency set for operator antialiasing		",
" ang=180		Maximum dip angle allowed in the image			",
" apet=45		aperture open angle for summation			",
" alias=0		=1 to set the anti-aliasing filter			",
" verbose=1		=1 to print some useful information			",
"										",
" Notes:									",
"										",
" The information needed in the computation of the weighting factor		",
" in Kirchhoff inversion formula includes traveltime, amplitude, 		",
" and Beylkin determinant at each grid point for each source/receiver		",
" point. For a 3-D nonzero common-offset inversion, the Beylkin			",
" determinant is computed from a 3x3 matrix with each element 			",
" containing a sum of quantities from the source and the receiver.		",
" The nine elements in the Beylkin matrix for each source/receiver		",
" can be determined by eight quantities. These quantities can be		",
" computed by the dynamic ray tracer. Tables of traveltime, amplitude,		",
" and Beylkin matrix elements from each source and receiver are			",
" pre-computed and stored in files.						",
"										",
" For each trace, tables of traveltime, amplitude and Beylkin matrix		",
" at the source and receiver location are input or interpolated from		",
" neighboring tables. For the computation of weighting factor, linear		",
" interpolation is used to determine the weighting factor at each		",
" output grid point, and weighted diffraction summation is then 		",
" applied. For each midpoint, the traveltimes and weight factors are		",
" calculated in the horizontal range of (xm-nxb*dx-z*tan(apet),			",
" xm+nxb*dx+z*tan(apet)).							",
"										",
" Offsets are signed - may be positive or negative. 				", 
"										",
NULL};

/*
 * This algorithm is based on the inversion formulas in chaper 5 of
 * _Mathematics of Multimensional Seismic Migration, Imaging and Inversion_ 
 * (Springer-Verlag, 2000), by Bleistein, N., Cohen, J.K.
 * and Stockwell, Jr., J. W.
 */
/**************** end self doc ***********************************/

#define TINY 0.0000001 /*avoid devide by zero*/
#define LARGE 1000000

void GetFileNameType1(char *infile,float xs,float ys,float xs0,float xs1,float ys0,float ys1,char *xsfile); 
void GetFileNameType2(char *infille,float xs,float xs0,float xs1,char *xsfile);
short ReadTable(char* infile,int nxt,int nyt,int nzt,float*** outArray); 
short ReadVelocity(char *vfile,int nxv,int nyv,int nzv,float fxv,float dxv,
float fyv,float dyv,float fzv,float dzv,int nxo,int nyo,int nzo,float fxo,
float dxo, float fyo,float dyo,float fzo,float dzo,float ***vout); 

/*segy trace */ 
segy tr, tro;

int
main (int argc, char **argv)
{
	int 	ixm,nxm,iym,nt,nxo,nyo,nzo,ixo,iyo,izo,nsamp,
		nym,nxb,i,ix,ixg,iy,iyg,iz,
		verbose,alias,geo_type,com_type,
		nxt,nyt,nzt,*aperture,nxv,nyv,nzv;
	float   temp,xo,yo,zo,xi,xig,zi,sx,sz,sxg,tsd,tgd,yi,sy,yig,syg,
		xs_beg,xg_beg,ys_beg,fym,
		amps,ampg,weight,
		dym,samp,smax,fmax,ang,cosa,
		dxm,dt,fxm,fxo,dxo,fyo,dyo,fzo,dzo,offs,xm,
		odt, smaxx,smaxy,
		dxt,dyt,dzt,apet,dxv,dyv,dzv,fxv,fyv,fzv,
		pxs,pys,pzs,pxg,pyg,pzg, 
		a1s,a1g,b1s,b1g,
		d11,d12,d13,d21,d22,d23,d31,d32,d33,
		d21s,d22s,d23s,d31s,d32s,d33s,
		d21g,d22g,d23g,d31g,d32g,d33g,
		ctheta,det,aliasx,aliasy,
		xs,xg,ys,vs,xt0,xt1,yt0,yt1,***outtrace=NULL, 
		***txs=NULL,***txg=NULL,***ampxs=NULL,***ampxg=NULL,
		***a1xs=NULL,***a1xg=NULL,***b1xs=NULL,***b1xg=NULL,
		***d21xs=NULL,***d21xg=NULL,***d22xs=NULL,
		***d22xg=NULL,***d23xs=NULL,***d23xg=NULL,
		***d31xs=NULL,***d31xg=NULL,***d32xs=NULL,
		***d32xg=NULL,***d33xs=NULL,***d33xg=NULL,***vout=NULL;
	char *vfile="";
	char *tfile="";
	char *ampfile="";
	char *d21file="";
	char *d22file="";
	char *d23file="";
	char *d31file="";
	char *d32file="";
	char *d33file="";
	char *a1file="";
	char *b1file="";

	char xsfile[80];
	char xgfile[80];

	/* hook up getpar to handle the parameters */
	initargs(argc, argv);
	requestdoc(1);


	/* get information from the first header */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	temp = tr.dt;
	if (!getparfloat("dt",&dt)) dt = temp*0.000001;
	if (dt<0.0000001) err("dt must be positive!\n");
	if (!getparfloat("offs",&offs)) offs = tr.offset;
	if (!getparfloat("dxm",&dxm)) dxm = tr.d2;
	if  (dxm<0.0000001) err("dxm must be positive!\n");
	
	/* get parameters for input data*/
	if (!getparint("geo_type",&geo_type)) 
			err("must specify geometry type!\n");
	/* 1---general velocity distribution v(x,y,z)
	   2---v(x,z) medium
	   3---v(z) medium  */
	if (!getparint("com_type",&com_type)) 
			err("must specify computation type!\n");
	/* 1---only needs traveltime table, weight=1.0
	   2---traveltime, a1, b1 tables,   weight=ctheta
	   3---everything		   weight=det/as/ag/(1+ctheta) */  
	if (!getparint("nym",&nym)) err("must specify nym!\n");
	if (!getparint("nxm",&nxm)) err("must specify nxm!\n");
	if (!getparfloat("fxm",&fxm)) fxm = 0.0;
	if (!getparfloat("dxm",&dxm)) dxm = 50.;
	if (!getparfloat("fym",&fym)) fym = 0.0;
	if (!getparfloat("dym",&dym)) dym=50.;

	/* get parameters for velocity model*/
	if (!getparint("nyv",&nyv)) {
	    if (geo_type ==1 ) err("must specify nyv!\n");
	    else nyv =1;
	}
	if (!getparint("nzv",&nzv)) err("must specify nzv!\n");
	if (!getparint("nxv",&nxv)) {
	    if (geo_type == 1 || geo_type ==2) err("must specify nxv!\n");
	    else nxv =1;
	}
	if (!getparfloat("fxv",&fxv)) fxv = 0.0;
	if (!getparfloat("dxv",&dxv)) dxv = 50.;
	if (!getparfloat("fyv",&fyv)) fyv = 0.0;
	if (!getparfloat("dyv",&dyv)) dyv=50.;
	if (!getparfloat("fzv",&fzv)) fzv = 0.0;
	if (!getparfloat("dzv",&dzv)) dzv = 50.;
	if (!getparstring("vfile",&vfile)) err("must specify vfile!\n");
	
	/* get required parameters about wavefied tables*/
	if (!getparstring("tfile",&tfile)) err("must specify timefile!\n");
	if (!getparstring("ampfile",&ampfile) && com_type ==3) 
				err("must specify ampfile!\n");
	if (!getparstring("d21file",&d21file) && com_type ==3) 
				err("must specify d21file!\n");
	if (!getparstring("d22file",&d22file) && com_type ==3) 
				err("must specify d22file!\n");
	if (!getparstring("d23file",&d23file) && com_type ==3) 
				err("must specify d23file!\n");
	if (!getparstring("d31file",&d31file) && com_type ==3) 
				err("must specify d31file!\n");
	if (!getparstring("d32file",&d32file) && com_type ==3) 
				err("must specify d32file!\n");
	if (!getparstring("d33file",&d33file) && com_type ==3) 
				err("must specify d33file!\n");
	if (!getparstring("a1file",&a1file))  
	    if (com_type ==2 || com_type ==3  ) 
				err("must specify a1file!\n"); 
	if (!getparstring("b1file",&b1file)) 
	    if (com_type ==2 || com_type ==3  ) 
				err("must specify b1file!\n"); 

	if (!getparint("nyt",&nyt)) err("must specify nyt!\n");
	if (!getparint("nzt",&nzt)) err("must specify nzt!\n");
	if (!getparint("nxt",&nxt)) err("must specify nxt!\n");
	if (!getparfloat("dyt",&dyt)) dyt = 100.0;
	if (!getparfloat("dzt",&dzt)) dzt = 100.0;
	if (!getparfloat("dxt",&dxt)) dxt = 100.0;
	if (!getparfloat("xt0",&xt0)) xt0 = 0.0;
	if (!getparfloat("xt1",&xt1)) xt1 = 0.0;
	if (!getparfloat("yt0",&yt0)) yt0 = 0.0;
	if (!getparfloat("yt1",&yt1)) yt1 = 0.0;

	/* get optional parameters */
	if (!getparint("nxb",&nxb)) nxb = nxm/2;
	if (!getparfloat("apet",&apet)) apet=45.;
	if (!getparint("alias",&alias)) alias = 0;
	if (!getparint("verbose",&verbose)) verbose = 1;
	if (!getparfloat("fmax",&fmax)) fmax = 0.25/dt;
	if (!getparfloat("ang",&ang)) ang = 180;
	cosa = cos(ang*PI/180);

	/* get optional parameters for output image*/
	if (!getparfloat("fxo",&fxo)) fxo = 0.0;
	/* check the first output trace	*/
	if (fxo<fxm )
		err("Does it make sense that fxo<fxm? \n");
	if (!getparfloat("dxo",&dxo)) dxo = 15.;
	if (!getparint("nxo",&nxo)) nxo = 101;
	if (!getparfloat("fyo",&fyo)) fyo=0.0;
	if (!getparfloat("dyo",&dyo)) dyo = 15.;
	if (!getparint("nyo",&nyo)) nyo = 101;
	if (!getparfloat("fzo",&fzo)) fzo = 0.;
	if (!getparfloat("dzo",&dzo)) dzo = 15.;
	if (!getparint("nzo",&nzo)) nzo = 101;

	/* print out the read-in parameters */
	fprintf(stderr,"\n parameters have been read in...");
	fprintf(stderr,"\n geo_type=%d,com_type=%d",geo_type,com_type);
	fprintf(stderr,"\n dt=%f,nt=%d,offs=%f",dt,nt,offs);
	fprintf(stderr,"\n nxm=%d,nym=%d",nxm,nym);
	fprintf(stderr,"\n dxm=%f,dym=%f",dxm,dym);
	fprintf(stderr,"\n fxm=%f,fym=%f",fxm,fym);
	fprintf(stderr,"\n vfile=%s",vfile);
	fprintf(stderr,"\n nxv=%d,nyv=%d,nzv=%d",nxv,nyv,nzv);
	fprintf(stderr,"\n dxv=%f,dyv=%f,dzv=%f",dxv,dyv,dzv);
	fprintf(stderr,"\n fxv=%f,fyv=%f,fzv=%f",fxv,fyv,fzv);
	fprintf(stderr,"\n tfile=%s",tfile);
	fprintf(stderr,"\n ampfile=%s",ampfile);
	fprintf(stderr,"\n d21file=%s",d21file);
	fprintf(stderr,"\n d22file=%s",d22file);
	fprintf(stderr,"\n d23file=%s",d23file);
	fprintf(stderr,"\n d31file=%s",d31file);
	fprintf(stderr,"\n d32file=%s",d32file);
	fprintf(stderr,"\n d33file=%s",d33file);
	fprintf(stderr,"\n a1file=%s",a1file);
	fprintf(stderr,"\n b1file=%s",b1file);
	fprintf(stderr,"\n dxt=%f,nxt=%d",dxt,nxt);
	fprintf(stderr,"\n dyt=%f,nyt=%d",dyt,nyt);
	fprintf(stderr,"\n dzt=%f,nzt=%d",dzt,nzt);
	fprintf(stderr,"\n xt0=%f,xt1=%f",xt0,xt1);
	fprintf(stderr,"\n yt0=%f,yt1=%f",yt0,yt1);
	fprintf(stderr,"\n nxb=%d,ang=%f,apet=%f",nxb,ang,apet);
	fprintf(stderr,"\n dxo=%f,nxo=%d,fxo=%f",dxo,nxo,fxo);
	fprintf(stderr,"\n dyo=%f,nyo=%d,fyo=%f",dyo,nyo,fyo);
	fprintf(stderr,"\n dzo=%f,nzo=%d,fzo=%f",dzo,nzo,fzo);
	fprintf(stderr,"\n fmax=%f,alias=%d,verbose=%d\n",fmax,alias,verbose);
	
	odt = 1.0/dt;
	smax = 0.5/(MAX(dxm,dym)*fmax);
	aliasx = 0.5/dxm/fmax;
	aliasy = 0.5/dym/fmax;

	/* allocate space */
	outtrace = ealloc3float(nzo,nxo,nyo);

	
	/* initialize output traces	*/
	for(ixo=0; ixo<nxo; ++ixo)
	    for(iyo=0;iyo<nyo;++iyo)
		for(izo=0; izo<nzo; ++izo)
			outtrace[iyo][ixo][izo]=0.0;
	
	/* allocate space for traveltimes and other quantities	*/ 
	txs = ealloc3float(nzt,nxt,nyt);
	txg = ealloc3float(nzt,nxt,nyt);
	if (com_type == 3 ) {
	   ampxs = ealloc3float(nzt,nxt,nyt);
	   ampxg = ealloc3float(nzt,nxt,nyt);
	   d21xs = ealloc3float(nzt,nxt,nyt);
	   d21xg = ealloc3float(nzt,nxt,nyt);
	   d22xs = ealloc3float(nzt,nxt,nyt);
	   d22xg = ealloc3float(nzt,nxt,nyt);
	   d23xs = ealloc3float(nzt,nxt,nyt);
	   d23xg = ealloc3float(nzt,nxt,nyt);
	   d31xs = ealloc3float(nzt,nxt,nyt);
	   d31xg = ealloc3float(nzt,nxt,nyt);
	   d32xs = ealloc3float(nzt,nxt,nyt);
	   d32xg = ealloc3float(nzt,nxt,nyt);
	   d33xs = ealloc3float(nzt,nxt,nyt);
	   d33xg = ealloc3float(nzt,nxt,nyt);
	}
	if ( com_type==2 || com_type == 3 ) {
	   a1xs = ealloc3float(nzt,nxt,nyt);
	   a1xg = ealloc3float(nzt,nxt,nyt);
	   b1xs = ealloc3float(nzt,nxt,nyt);
	   b1xg = ealloc3float(nzt,nxt,nyt);
	}

	vout = ealloc3float(nzo,nxo,nyo);
	aperture = ealloc1int(nzo);

	/* initialize output traces	*/
	if (com_type == 3 ) {
		for(ix=0; ix<nxt; ++ix)
		for(iy=0;iy<nyt;++iy)
		for(iz=0; iz<nzt; ++iz) {
			ampxs[iy][ix][iz]=1.0;
			ampxg[iy][ix][iz]=1.0;
			}
	}

	/* Read in velocity file */
	ReadVelocity(vfile, nxv, nyv, nzv, fxv, dxv, fyv, dyv, fzv, dzv, 
		nxo, nyo, nzo, fxo, dxo, fyo, dyo, fzo, dzo, vout); 

	/* determine the aperture array */
	for(iz=0; iz<nzo; ++iz) {
		if (apet < 89.0 )
		aperture[iz] = nxb+(dzo*iz*tan(apet*PI/180)+0.5)/dxm;
		else
		aperture[iz] = nxb;
	}
	fprintf(stderr,"\nvfile have been read in...\n");
	
	if(ABS(offs)>2*nxb*dxm) err("\t band NXB is too small!");

	if((nxm-1)*dxm<(nxo-1)*dxo || (nym-1)*dym<(nyo-1)*dyo) err("\t not enough range in input  		for the imaging");			    				
	vs = vout[0][0][0];
	smaxx = 0.5*vs/dxm/fmax;
	smaxy = 0.5*vs/dym/fmax;