suharlan.c

su 的源代码库

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

/* SUHARLAN: $Revision: 1.6 $ ; $Date: 2003/06/09 16:17:07 $	*/

#include "su.h"
#include "segy.h"
#include "header.h"
#include "taup.h"
#include <time.h>
#include <signal.h>

/**************************** self documentation *****************************/
char *sdoc[] = {
"									",
" SUHARLAN - signal-noise separation by the invertible linear		",
"	    transformation method of Harlan, 1984			",
"									",
"   susnlinsep <infile >outfile  [optional parameters]			",
"									",
" Required Parameters:						 	",
" <none>								",
"									",
" Optional Parameters:							",
" FLAGS:								",
" niter=1	number of requested iterations				",
" anenv=1	=1 for positive analytic envelopes			",
"		=0 for no analytic envelopes (not recommended)		",
" scl=0		=1 to scale output traces (not recommended)		",
" plot=3	=0 for no plots. =1 for 1-D plots only			",
"		=2 for 2-D plots only. =3 for all plots			",
" norm=1	=0 not to normalize reliability values			",
" verbose=1	=0 not to print processing information			",
" rgt=2		=1 for uniform random generator				",
"		=2 for gaussian random generator			",
" sts=1		=0 for no smoothing (not recommended)			",
"									",
" tmpdir= 	 if non-empty, use the value as a directory path	",
"		 prefix for storing temporary files; else if the	",
"	         the CWP_TMPDIR environment variable is set use		",
"	         its value for the path; else use tmpfile()		",
"									",
" General Parameters:							",
" dx=20		offset sampling interval (m)				",
" fx=0	  	offset on first trace (m)				",
" dt=0.004	time sampling interval (s)				",
"									",
" Tau-P Transform Parameters:						",
" gopt=1	=1 for parabolic transform. =2 for Foster/Mosher	",
"		=3 for linear. =4 for absolute value of linear		",
" pmin1=-400	minimum moveout at farthest offset for fwd transf(ms)	",
" pmax1=400	maximum moveout at farthest offset for fwd transf(ms)	",
" pmin2=pmin1	minimum moveout at farthest offset for inv transf(ms)	",
" pmax2=pmax1	maximum moveout at farthest offset for inv transf(ms)	",
" np=100	number of p-values for taup transform			",
" prewhite=0.01	prewhitening value (suggested between 0.1 and 0,01)	",
" offref=2000	reference offset for p-values (m)			",
" depthref=500	reference depth for Foster/Mosher taup (if gopt=4)	",
" pmula=pmax1	maximum p-value preserved in the data (ms)		",
" pmulb=pmax1	minimum p-value muted on the data (ms)			",
" ninterp=0	number of traces to interpolate in input data		",
"									",
" Extraction Parameters:						",
" nintlh=50	number of intervals (bins) in histograms		",
" sditer=5	number of steepest descent iterations to compute ps	",
" c=0.04	maximum noise allowed in a sample of signal(%)		",
" rel1=0.5	reliability value for first pass of the extraction	",
" rel2=0.75	reliability value for second pass of the extraction	",
"								   	",
" Smoothing Parameters:							", 
" r1=10		number of points for damped lsq vertical smoothing	",
" r2=2		number of points for damped lsq horizontal smoothing	",
"								   	",
"								   	",
" Output Files:								",
" signal=out_signal 	name of output file for extracted signal	",
" noise=out_noise 	name of output file for extracted noise		",
"									",
" 									",
" Notes:								",
" The signal-noise separation algorithm was developed by Dr. Bill Harlan",
" in 1984. It can be used to separate events that can be focused by a	",
" linear transformation (signal) from events that can't (noise). The	",
" linear transform is whatever is well siuted for the application at	",
" hand. Here, only the discrete Radon transform is used, so the program	",
" is capable of separating events focused by that transform (linear,	",
" parabolic or time-invariantly hyperbolic). Should other transform be	",
" required, the changes to the program will be relatively		",
" straightforward.							",
"									",
" The reliability parameter is the most critical one to determine what	",
" to extract as signal and what to reject as noise. It should be tested	",
" for every dataset. The way to test it is to start with a small value,	",
" say 0.1 or 0.01. If too much noise is present in the extracted noise,	",
" it is too low. If too much signal was extracted, that is, part of the	",
" signal was lost, it is too big. All other parameters have good default",
" values and should perhaps not be changed in a first encounter with the",
" program. The transform parameters are also critical. They should be	",
" chosen such that no aliasing is present and such that the range of	",
" interesting slopes is spanned by the transform but not much more. The ",
" program suradon.c has more documentation on the transform paramters.	",
"									",
NULL};

/*
 * Credits:
 * 	Gabriel Alvarez CWP (1995) 
 *	Some subroutines are direct translations to C from Fortran versions
 * 	written by Dr. Bill Harlan (1984)
 *
 * References:
 *
 * 	Harlan, S., Claerbout, J., and Roca, F. (1984), Signal/noise
 *	separation and velocity estimation, Geophysics, v. 49, no. 11,
 *	p 1869-1880. 
 *
 * 	Harlan, S. (1988), Separation of signal and noise applied to
 *	vertical seismic profiles, Geophysics, v. 53, no. 7,
 *	p 932-946. 
 *
 *	Alvarez, G. (1995), Comparison of moveout-based approaches to
 *	ground roll and multiple suppression, MSc., Department of 
 *	Geophysics, Colorado School of Mines, (Chapter 3 deals
 *	exclusively with this method).
 *
 */
/************************** end self doc *************************************/

/* Prototypes of functions used internally */
void separate_signal_noise(int verbose, int rgt, int sts, int norm, int anenv,
	int scl, int plot, int seed, int sditer, int niter, int nintlh,
	float dt, int nt, int nx, int gopt, float prewhite, int ninterp,
	float offref, float depthref, float interoff, float pmin1, float pmax1,
	float pmin2, float pmax2, int np, float dx, float c, float rel1,
	float rel2, float r1, float r2,float fx, float pmula, float pmulb,
	float **in_traces, float **signal, float **noise);

void extract_signal (int verbose, int norm, int anenv, int sts, int plot,
	int sditer, int nintlh, float c, float r, int nx, int nt, float r1,
	float r2, float **traces, float **rand_traces);

void expected_value(int nintlh, int nx, int nt, int rindex, float *pspn,
	float *pn, float *ps, float *xamps, float **traces, float *Esd);

void compute_reliability(int norm, int nx, int nt, int nintlh, int rindex,
	float c, float *ps, float *pn, float *pspn, float **traces,
	float *reliability);

void zero_noisy_samples(int anenv, int sts, int nx, int nt, float r1,
	float r2, float r, int nintlh, float *amps, float **traces,
	float **anenv_traces, float *rel);

void compute_histogram_stuff(int verbose, int nx, int nt, int nintlh,
	float **traces, float **rand_traces, int *rindex, float *amps,
	float *xamps);

void make_histogram(int nintlh, int nt, int nx, float **traces, float *amps, 
	float *pdf);

void compute_analytic_envelopes(int sgn, int nx, int nt, float **traces, 
	float **analytic_envelopes);

void scale_traces(int scl, int nx, int nt, float **in_traces, 
	float **out_traces);

void scale_one_trace(int ns, float *in_trace, float *out_trace);

void matrix_transpose(int n1, int n2, float **matrix, float **tr_matrix);

float dot_product(int ns, float *vector1, float *vector2);

void compute_max_min_sum(int nx, int nt, float *min, float *max, float *sum,
	float **data);

void plot_one_d(int npoints, float *xamps, float *data, char *plotname);

void plot_two_d(int npoints, float *data, char *plotname);

void conv1(int nx, int fx, float *x, int ny, int fy, float *y, int nz,
	int fz, float *z, int flag, float perc);

void deconvolve_histograms(int nsamples, int mean_index, int niter,
	float *pd, float *pn, float *ps);

void gradient(int ns, int si, float rmax, float *pd, float *pn,
	float *ps, float *grad);

void cross_entropy(int ns, int si, float rmax, float *pd, float *pn,
	float *ps, float *fvalue);

float divide(float rmax, float a, float b);

void golden_search(float fvalue, int *iter, float *xvalue, float *alpha);

static void closefiles(void);

/* Globals (so can trap signal) defining temporary disk files */
char tracefile[BUFSIZ];	/* filename for the file of traces	*/
char headerfile[BUFSIZ];/* filename for the file of headers	*/
FILE *tracefp;		/* fp for trace storage file		*/
FILE *headerfp;		/* fp for header storage file		*/


segy tro1,tro2,tro3;

int
main(int argc, char **argv)
{
	int ix, it;		/* loop counters */
	int anenv;		/* =1 for analytic envelopes */
	int scl;		/* =1 to apply trace scaling */
	int plot;		/* flag for producing plots */
	int seed;		/* seed for random number generator */
	int verbose;		/* flag to print processing information */
	int rgt;		/* random number generator type */
	int sts;		/* vertical and horizontal smoothing */
	int sditer;		/* # of steepest descent iterations for ps */
	int nintlh;		/* number of intervals per local histogram */
	int norm;		/* option to normalize reliability indicator */
	float c;		/* maximum error allowed for reliability */
	float rel1;		/* minimum allowed reliability (first pass) */
	float rel2;		/* minimum allowed reliability (second pass) */
	int nt;			/* number of time samples */
	int niter;		/* number of iterations */
	int ntr;		/* number of traces */
	float dx;		/* horizontal sampling interval */
	float fx;		/* offset of first trace */
	float dt;		/* time sampling interval (ms) */
	float r1;		/* vertical smoothing factor for dlsq method */
	float r2;		/*horizontal smoothing factor for dlsq method*/
	int ninterp;		/* traces to interpolate for tau-p transform*/
	int gopt;		/* options for offset function g(x) in taup */
	float offref;		/* offset reference for tau-p transform */
	float pmin1;		/* min moveout in ms at ref offset for fwd tr*/
	float pmax1;		/* max moveout in ms at ref offset for fwd tr*/
	float pmin2;		/* min moveout in ms at ref offset for inv tr*/
	float pmax2;		/* max moveout in ms at ref offset for inv tr*/
	float prewhite;		/* prewhitening factor for tau-p transform */
	float depthref;		/* reference depth if gopt=2 */
	float pmula;		/* maximum slope to keep in the data */
	float pmulb;		/* minimum slope rejected from the data */
	float interoff;		/* intercept offset for tau-p times */
	int np;			/* number of slopes (traces) in Tau-p domain */
	float **traces; 	/* Array for input traces */
	float **out_signal;	/* Array of extracted signal */
	float **out_noise; 	/* Array of extracted noise */
	char *signalfile="";	/* file name for output signal */
	char *noisefile="";	/* file name for output noise */
	FILE *signal_file;	/* File pointer to output signal */
	FILE *noise_file;	/* File pointer to output noise */
	char *tmpdir;		/* directory path for tmp files	*/
	cwp_Bool istmpdir=cwp_false;/* true for user-given path */
	
	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(1);

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

/* get general flags and set their defaults */
	if (!getparint("anenv",&anenv))			anenv	= 1;
	if (!getparint("scl",&scl))			scl	= 0;
	if (!getparint("plot",&plot))			plot	= 3;
	if (!getparint("norm",&norm))			norm	= 1;
	if (!getparint("verbose",&verbose))		verbose	= 1;
	if (!getparint("rgt",&rgt))			rgt	= 2;
	if (!getparint("sts",&sts))			sts	= 1;
	if (!getparint("gopt",&gopt))			gopt	= 1;

	/* get general parameters and set their defaults */
	if (!getparfloat("dx",&dx)) 			dx	= 20.0;
	if (!getparfloat("fx",&fx)) 			fx	= 0.0;
	if (!getparint("nintlh",&nintlh))		nintlh	= 50;
	if (!getparint("sditer",&sditer))		sditer	= 5;
	if (!getparfloat("c",&c)) 			c	= 0.05;
	if (!getparfloat("rel1",&rel1)) 		rel1	= 0.5;
	if (!getparfloat("rel2",&rel2)) 		rel2	= 0.75;
	if (!getparint("niter",&niter))			niter	= 1;
	if (!getparfloat("dt",&dt))			dt	= 0.004;

	/* get smoothing parameters and set their defaults */
	if (!getparfloat("r1",&r1))			r1	= 10.;
	if (!getparfloat("r2",&r2))			r2	= 2.;

	/* get suradon taup parameters and set their defaults */
	if (!getparfloat("pmin1",&pmin1))		pmin1	= -400.0;
	if (!getparfloat("pmax1",&pmax1))		pmax1	= 400.0;
	if (!getparfloat("pmin2",&pmin2))		pmin2	= pmin1;
	if (!getparfloat("pmax2",&pmax2))		pmax2	= pmax1;
	if (!getparint("np",&np)) 			np	= 100;
	if (!getparfloat("prewhite",&prewhite))		prewhite= 0.01;
	if (!getparfloat("offref",&offref))		offref	= 2000.;
	if (!getparfloat("interoff",&interoff))		interoff= 0.;
	if (!getparfloat("depthref",&depthref))		depthref= 0.;
	if (!getparfloat("pmula",&pmula))		pmula	= pmax1;
	if (!getparfloat("pmulb",&pmulb))		pmulb	= pmax1;
		if (!getparint("ninterp",&ninterp))	ninterp	= 0;

	/* get names of output files */
	if (!getparstring("signal",&signalfile))    signalfile="out_signal";
	if (!getparstring("noise",&noisefile))	    noisefile="out_noise";

	/* get random generator seed */
	if (!getparint("seed", &seed)) { /* if not supplied, use clock */
		if (-1 == (seed = (unsigned int) time((time_t *) NULL))) {
			seed=1;
			warn("time() failed to set seed, setting it to one");
		}
	}

	/* Look for user-supplied tmpdir */
	if (!getparstring("tmpdir",&tmpdir) &&
	    !(tmpdir = getenv("CWP_TMPDIR"))) tmpdir="";
	if (!STREQ(tmpdir, "") && access(tmpdir, WRITE_OK))
		err("you can't write in %s (or it doesn't exist)", tmpdir);

	/* Store traces in tmpfile while getting a count */
	if (STREQ(tmpdir,"")) {
		tracefp = etmpfile();
		headerfp = etmpfile();
		if (verbose) warn("using tmpfile() call");
	} else { /* user-supplied tmpdir */
		char directory[BUFSIZ];
		strcpy(directory, tmpdir);
		strcpy(tracefile, temporary_filename(directory));
		strcpy(headerfile, temporary_filename(directory));
		/* Trap signals so can remove temp files */
		signal(SIGINT,  (void (*) (int)) closefiles);
		signal(SIGQUIT, (void (*) (int)) closefiles);
		signal(SIGHUP,  (void (*) (int)) closefiles);
		signal(SIGTERM, (void (*) (int)) closefiles);
		tracefp = efopen(tracefile, "w+");
		headerfp = efopen(headerfile, "w+");
      		istmpdir=cwp_true;		
		if (verbose) warn("putting temporary files in %s", directory);
	}

	ntr = 0;
	do {
		++ntr;
		fwrite(&tro1, 1, HDRBYTES, headerfp);
		fwrite(tro1.data, FSIZE, nt, tracefp);
	} while (gettr(&tro1));

	/* allocate space */
	traces = alloc2float(nt,ntr);
	out_signal = alloc2float(nt,ntr);
	out_noise = alloc2float(nt,ntr);
	
	/* load traces into an array and close temp file */
	rewind(tracefp);
	for (ix=0; ix<ntr; ix++)
		fread (traces[ix], FSIZE, nt, tracefp);
	efclose (tracefp);

	/* separate bed reflections from diffractions and noise */
	separate_signal_noise (verbose, rgt, sts, norm, anenv, scl, plot, seed,
		sditer, niter, nintlh, dt, nt, ntr, gopt, prewhite, ninterp,
		offref, depthref, interoff, pmin1, pmax1, pmin2, pmax2, np,
		dx, c, rel1, rel2, r1, r2, fx, pmula, pmulb, traces,
		out_signal, out_noise);

	/* write extracted reflections */
	if (*signalfile !='\0') {
		if ((signal_file=fopen(signalfile,"w"))==NULL)
			err("cannot open signal file=%s\n",signalfile);
		erewind(headerfp);
		{	register int itr;
			for (itr=0; itr<ntr; itr++) {