suradon.c

su 的源代码库

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

/* SURADON: $Revision: 1.16 $ ; $Date: 2005/12/07 17:11:15 $	*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SURADON - compute forward or reverse Radon transform or remove multiples",
"           by using the parabolic Radon transform to estimate multiples",
"           and subtract.						",
"									",
"     suradon <stdin >stdout [Optional Parameters]			",
"									",
" Optional Parameters:							",
" choose=0    0  Forward Radon transform				",
"             1  Compute data minus multiples				",
"             2  Compute estimate of multiples				",
"             3  Compute forward and reverse transform			",
"             4  Compute inverse Radon transform			",
" igopt=1     1  parabolic transform: g(x) = offset**2			",
"             2  Foster/Mosher psuedo hyperbolic transform		",
"                   g(x) = sqrt(depth**2 + offset**2)			",
"             3  Linear tau-p: g(x) = offset				",
"             4  abs linear tau-p: g(x) = abs(offset)			",
" offref=2000.    reference maximum offset to which maximum and minimum	",
"                 moveout times are associated				",
" interoff=0.     intercept offset to which tau-p times are associated	",
" pmin=-200       minimum moveout in ms on reference offset		",
" pmax=400        maximum moveout in ms on reference offset		",
" dp=16           moveout increment in ms on reference offset		",
" pmula=80        moveout in ms on reference offset where multiples begin",
"                     at maximum time					",
" pmulb=200       moveout in ms on reference offset where multiples begin",
"                     at zero time					",
" depthref=500.   Reference depth for Foster/Mosher hyperbolic transform",
" nwin=1          number of windows to use through the mute zone	",
" f1=60.          High-end frequency before taper off			",
" f2=80.          High-end frequency					",
" prewhite=0.1    Prewhitening factor in percent.			",
" cdpkey=cdp      name of header word for defining ensemble		",
" offkey=offset   name of header word with spatial information		",
" nxmax=120       maximum number of input traces per ensemble		",
" ltaper=7	  taper (integer) for mute tapering function		",
"									",
" Optimizing Parameters:						",
" The following parameters are occasionally used to avoid spatial aliasing",
" problems on the linear tau-p transform.  Not recommended for other	",
" transforms...								",
" ninterp=0      number of traces to interpolate between each input trace",
"                   prior to computing transform			",
" freq1=3.0      low-end frequency in Hz for picking (good default: 3 Hz)",
"                (Known bug: freq1 cannot be zero) ",
" freq2=20.0     high-end frequency in Hz for picking (good default: 20 Hz)",
" lagc=400       length of AGC operator for picking (good default: 400 ms)",
" lent=5         length of time smoother in samples for picker		",
"                     (good default: 5 samples)				",
" lenx=1         length of space smoother in samples for picker		",
"                     (good default: 1 sample)				",
" xopt=1         1 = use differences for spatial derivative		",
"                        (works with irregular spacing)			",
"                0 = use FFT derivative for spatial derivatives		",
"                      (more accurate but requires regular spacing and	",
"                      at least 16 input tracs--will switch to differences",
"                      automatically if have less than 16 input traces)	",
"									",
NULL};

/* Credits:
 *	CWP: John Anderson (visitor to CSM from Mobil) Spring 1993
 *
 * Multiple removal notes:
 *	Usually the input data are NMO corrected CMP gathers.  The
 *	first pass is to compute a parabolic Radon transform and
 * 	identify the multiples in the transform domain.  Then, the
 * 	module is run on all the data using "choose=1" to estimate
 * 	and subtract the multiples.  See the May, 1993 CWP Project 
 *	Review for more extensive documentation.
 *
 * NWIN notes:
 *	The parabolic transform runs with higher resolution if the
 * 	mute zone is honored.  When "nwin" is specified larger than
 *   	one (say 6), then multiple windows are used through the mute
 * 	zone.  It is assumed in this case that the input data are
 * 	sorted by the offkey header item from small offset to large
 * 	offset.  This causes the code to run 6 times longer.  The
 *      mute time is taken from the "muts" header word. 
 *      You may have to manually set this header field yourself, if
 *      it is not already set.
 *
 * References:
 * Anderson, J. E., 1993, Parabolic and linear 2-D, tau-p transforms
 *       using the generalized radon tranform, in May 11-14, 1993
 *       Project Review, Consortium Project on Seismic Inverse methods
 *       for Complex Structures, CWP-137, Center for Wave Phenomena
 *       internal report.
 * Other References cited in above paper:
 * Beylkin, G,.1987, The discrete Radon transform: IEEE Transactions
 *       of Acoustics, Speech, and Signal Processing, 35, 162-712.
 * Chapman, C.H.,1981, Generalized Radon transforms and slant stacks:
 *       Geophysical Journal of the Royal Astronomical Society, 66,
 *       445-453.
 * Foster, D. J. and Mosher, C. C., 1990, Multiple supression
 *       using curvilinear Radon transforms: SEG Expanded Abstracts 1990,
 *       1647-1650.
 * Foster, D. J. and Mosher, C. C., 1992, Suppression of multiples
 *       using the Radon transform: Geophysics, 57, No. 3, 386-395.
 * Gulunay, N., 1990, F-X domain least-squares Tau-P and Tau-Q: SEG
 *       Expanded Abstracts 1990, 1607-1610.
 * Hampson, D., 1986, Inverse velocity stacking for multiple elimination:
 *       J. Can. Soc. Expl. Geophs., 22, 44-55.
 * Hampson, D., 1987, The discrete Radon transform: a new tool for image
 *       enhancement and noise suppression: SEG Expanded Abstracts 1978,
 *       141-143.
 * Johnston, D.E., 1990, Which multiple suppression method should I use?
 *       SEG Expanded Abstracts 1990, 1750-1752.
 *
 * Trace header words accessed: ns, dt, cdpkey, offkey, muts
 */
/**************** end self doc ********************************/

static void forward_p_transform(VND *vnda,VND *vndb,int nx, int nt, float *g,
	float dt, int ntfft, int np, float pmin, float dp,
	float *mutetime, float *offset, int nk,float f1,
	float f2,float prewhite);

static void inverse_p_transform(VND *vnda,int nx, float *g,
	float dt, int ntfft, int np, float pmin, float dp, int ip1,
	float f1,float f2);
static void compute_r(float w, int nx, float *g, int np, float dp, complex *r);
static void compute_rhs(float w, int nx, float *g, complex *data, int np, 
		float pmin, float dp, complex *rhs);
static int ctoep(int n, complex *r, complex *a, complex *b,
		complex *f, complex *wrk);
static int ctoephcg(int niter, int n, complex *r, complex *a, complex *b,
		complex *wrk1, complex *wrk2, complex *wrk3, complex *wrk4 );
static float rcdot(int n, complex *a, complex *b);
static void htmul(int n, complex *a, complex *x, complex *y);
static float freqweight(int j, float df, float f1, float f2);
static float gofx(int igopt, float offset, float intercept_off,
		  float refdepth);
static void taupmute(int ip,int ipa,int ipb,int nt, int ltap, float *rt);
static void jea_xinterpolate(VND *vndorig, VND *vndinterp, int ninterp, 
		int nt, int nx, float freq1, float freq2, int lagc, 
		int lent, int lenx, int xopt, float dt, int iopt);
static void runav(int n,int len,float *a,float *b);

segy tr;
segy tro;
int main(int argc, char **argv)
{
	char *cdpkey;		/* key denoting the ensemble */
	char *offkey;		/* key denoting trace labeling in an ensemble */
	char *headerfile;	/* temporary file containing trace headers */
	char *fname;	
	int j;
	int it;
	int icount;
	int oldcmp;
	int icmp;
	int nxmax;
	int nx;
	int nxinterp;
	int nxout=0;
	int np;
	int ipa;
	int ipb;
	int ix;
	int ninterp;
	int k;
	int nt;
	int lent;
	int lenx;
	int lagc;
	int xopt=0;
	int ntfft;
	int nxm;
	int nmax;
	int choose;
	int nk;
	int igopt;
	int ltaper;
	int cdpindex;
	int offindex;
	int iend;
	int ieod;
	float offref;
	float depthref;
	float intercept_off;
	float pmin;
	float pmax;
	float pmula;
	float pmulb;
	float dp;
	float dt;
	float freq1;
	float freq2;
	float f1;
	float f2;
	float prewhite;
	float fac;
	float d;
	float *rt;
	float *xin;
	float *offset;
	float *mutetime;
	float *g;
	float *gg;
	float *trace;
	Value hdrwd;
	FILE *headerfp;
	VND *vndorig;
	VND *vndinterp;
	VND *vndresult;

	initargs(argc, argv);
	requestdoc(1);

	if (!getparstring("cdpkey",&cdpkey)) cdpkey="cdp";
	if (!getparstring("offkey",&offkey)) offkey="offset";
	if (!getparint("choose",&choose)) choose=0;
	if (!getparint("ninterp",&ninterp)) ninterp=0;
	if (!getparint("nwin",&nk)) nk=1;
	if (!getparint("igopt",&igopt)) igopt=1;
	if (!getparint("nxmax",&nxmax)) nxmax=240;
	if (!getparint("ninterp",&ninterp)) ninterp=1;
	if (!getparint("lagc",&lagc))lagc=400;
	if (!getparint("lent",&lent)) lent=5;
	if (!getparint("lenx",&lenx)) lenx=1;
	if (!getparfloat("freq1",&freq1)) freq1=4.;
	if (!getparfloat("freq2",&freq2)) freq2=20.;
	if (!getparfloat("offref",&offref)) offref=2000.;
	if (!getparfloat("f1",&f1)) f1=60.;
	if (!getparfloat("f2",&f2)) f2=80.;
	if (!getparfloat("pmin",&pmin)) pmin=-200.;
	if (!getparfloat("pmax",&pmax)) pmax=400.;
	if (!getparfloat("dp",&dp)) dp=16.;
	if (!getparfloat("pmula",&pmula)) pmula=80.;
	if (!getparfloat("pmulb",&pmulb)) pmulb=200.;
	if (!getparfloat("depthref",&depthref)) depthref=500.;
	if (!getparfloat("interoff",&intercept_off)) intercept_off=0.;
	if (!getparfloat("prewhite",&prewhite)) prewhite=0.1;
	if (!getparint("ltaper",&ltaper)) ltaper=7;

	cdpindex=getindex(cdpkey);
	offindex=getindex(offkey);
	if(nk<0) nk=1;
	fac=1000.*gofx(igopt,offref,intercept_off,depthref);
	pmin/=fac;
	pmax/=fac;
	dp/=fac;
	pmula/=fac;
	pmulb/=fac;
	ipa=( pmula -  pmin)/ dp;
	ipb=( pmulb -  pmin)/ dp;
	np=1+( pmax -  pmin)/ dp;
	if(np<1)err("Range of PMIN and PMAX invalid");
	if(choose==0) ipa=0;
	if(choose==3) ipa=0;
	if(choose==4) ipa=0;
	if(ipa<0) ipa=0;

	if(!gettr(&tr)) err("Can't get first trace \n");
	nt=tr.ns;
	dt = ((double) tr.dt)/1000000.0;
	gethval(&tr,cdpindex,&hdrwd);
	oldcmp=hdrwd.i;

	ntfft=npfar(nt);
	nxinterp = (1+ ninterp)*(nxmax-1)+1;
	nmax=MAX(nxinterp, np);
	nmax=MAX(nmax, ntfft+4);

	offset   = (float *)VNDemalloc(nxmax*sizeof(float),
				       "suradon_main: offset");
	rt       = (float *)VNDemalloc(nmax*sizeof(float),
				       "suradon_main: rt");
	trace    = (float *)VNDemalloc(nt*sizeof(float),
				       "suradon_main: trace");
	xin      = (float *)VNDemalloc(nxmax*sizeof(float),
				       "suradon_main: xin");
	g        = (float *)VNDemalloc(nxmax*sizeof(float),
				       "suradon_main: g");
	gg       = (float *)VNDemalloc(nxinterp*sizeof(float),
				       "suradon_main: gg");
	mutetime = (float *)VNDemalloc(nxmax*sizeof(float),
				       "suradon_main: mutetime");

	headerfile=VNDtempname("radontmp");
	if( (headerfp = fopen(headerfile,"w+"))==NULL) 
		err("couldn't open temp file for trace headers");
	fname = VNDtempname("radontmp");
	vndorig = V2Dop(2,1000000,sizeof(float),fname,nt, nxmax);
	VNDfree(fname,"suradon_main: fname 1");
	fname = VNDtempname("radontmp");
	vndinterp = V2Dop(2,2000000,sizeof(float),fname,nt, nxinterp);
	VNDfree(fname,"suradon_main: fname 2");
	fname = VNDtempname("radontmp");
	nmax=MAX( nxinterp, np);
	vndresult = V2Dop(2,1000000,sizeof(float),fname,ntfft+2, nmax);
	VNDfree(fname,"suradon_main: fname 3");

	icount=0;
	icmp=0;
	nx=0;
	iend=0;
	ieod=0;
	do {
		gethval(&tr,cdpindex,&hdrwd);
		if(hdrwd.i!=oldcmp || ieod) { /* process this ensemble */
		   icmp++;
		   nxm= nx-1;
		   if( nx>1) {
			nxinterp=1 + nxm*(ninterp+1);
			if( choose==4) { /* count number of original offsets */
				k=1;
				for(j=1;j<nx;j++) {
					if(fabs(offset[j]-offset[j-1])>0.001)k++;
				}
				np=nx;
				nx=k;
				for(j=0;j< np;j++){
					V2Dr0(vndorig,j,(char *)rt,1001);
					V2Dw0(vndresult,j,(char *)rt,1002);
				}
			}else{
				jea_xinterpolate( vndorig, vndinterp,
					 ninterp, nt, nx,
					 freq1, freq2, lagc,
					 lent, lenx, xopt, dt,0);
				d=1./(1.+ ninterp);
				for(j=0;j<nxinterp;j++) rt[j]=j*d;
				for(j=0;j<nx;j++) xin[j]=j;
				intlin(nx,xin,offset,offset[0],
					 offset[ nxm],nxinterp,rt,gg);
				for(j=0;j<nxinterp;j++)
					 gg[j]=gofx(igopt,gg[j],
						 intercept_off,depthref);
				forward_p_transform( vndinterp,		
					 vndresult, nxinterp,nt,gg,dt,
					 ntfft,np,pmin,dp,mutetime, 
					 offset,nk,f1,f2,prewhite);  
				nxout= np;
			} 
			if(choose>=1) {
				if( choose==1){
					/* do tau-p mute here */
					for(j=0;j< ipb;j++) {
						V2Dr0( vndresult,j,
						      (char *) rt,1003);
						taupmute(j,ipa,ipb,
							 ntfft,ltaper,rt);
						V2Dw0( vndresult,j,
						      (char *) rt,1004);
					}
				}
				inverse_p_transform( vndresult, nx,
					 g, dt, ntfft, np,
					 pmin, dp, ipa, f1, f2);
				if(choose==1){
				    for(ix=0;ix< nx;ix++) {
					V2Dr0( vndorig,ix,(char *)trace,1005);
					V2Dr0( vndresult,ix,(char *) rt,1006);
					for(it=0;it< nt;it++) 
						 rt[it]=trace[it]- rt[it];
					V2Dw0( vndresult,ix,(char *) rt,1007);
				    }
				}