suharlan.c

su 的源代码库

				fread(&tro1, 1, HDRBYTES, headerfp);
				for (it=0; it<nt; it++)
					tro1.data[it]=out_signal[itr][it];
				fputtr(signal_file, &tro1);
			}
			efclose (signal_file);
		}
	}

	/* write extracted diffractions */
	if (*noisefile!='\0') {
		if ((noise_file=efopen(noisefile,"w"))==NULL)
			err("cannot open noise file=%s\n",noisefile);
		erewind(headerfp);
		{	register int itr;
			for (itr=0; itr<ntr; itr++) {
				efread(&tro2, 1, HDRBYTES, headerfp);
				for (it=0; it<nt; it++)
					tro2.data[it]=out_noise[itr][it];
				fputtr(noise_file, &tro2);
			}
			efclose (noise_file);
		}
	}

	/* free workspace */
	free2float(out_signal);
	free2float(out_noise);
	free2float(traces);
	efclose(headerfp);
	if (istmpdir) eremove(headerfile);
	if (istmpdir) eremove(tracefile);

	return(CWP_Exit());
}

/******************************************************************************

	Subroutine to separate (extract) focused energy (signal)
	from unfocused energy (noise), after applying an invertible
			linear transformation.

******************************************************************************/
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)
/******************************************************************************
Input parameters:
Flags:
verbose		=1 to print processing information 
		=0 not to
rgt		=1 for uniform random generator
		=2 for gaussian random generator
anenv		=0 not to use analytic envelopes
		=1 to use signed analytic envelopes for extractions
		=2 to use positive analytic envelopes for extractions
		=3 to use negative analytic envelopes for extractions
scl		=0 do not apply trace scaling
		=1 apply trace by trace scaling to extracted signal
		=2 apply scaling by time slices
		=3 apply both scalings
plot		=0 not to produce any plots
		=1 to produce plots of histograms, Esd and reliability only
		=2 to also produce 2-D plots (in t-x and taup domains) 
		=3 to produce only 2-D plots (in t-x and taup domains) 
gopt		=1 for parabolic tau-p transform: g(x)=offset**2
		=2 for Foster/Mosher pseudo-hyperbolic tau-p transform
			g(x)=sqrt(ref_depth**2+offset**2)-ref_depth
		=3 linear tau-p: g(x)=offset
		=4 abs linear tau-p: g(x)=abs(offset)
		=5 for new pseudo-hyperbolic tau-p transform
			g(x)=1/ref_vel*sqrt(ref_vel**2+offset**2)


General Parameters:
seed		seed for random number generator
niter		number of iterations for extraction process
nintlh		number of intervals (bins) for making histograms
dt		time sampling interval
nt		number of time samples
nx		number of horizontal samples (traces)
fx		offset on first trace
c		maximum allowed error for reliability computation (in percent)
rel1		minimum reliability value to accept a sample as reliable 
		(first pass)
rel2		minimum reliability value to accept a sample as reliable 
		(second pass)

Taup transform parameters:
prewhite	prewhitening factor to stabilize inverse tau-p transform
		(in percent)
ninterp		number of traces to interpolate between each pair of input 
		traces prior to taup computation
offref		reference maximum offset to which tau-p times are associated 
pmin	    minimum moveout (slope) in ms on reference offset
pmin	    maximum moveout (slope) in ms on reference offset
np		number of slopes (traces) in tau-p domain
traces	  input traces

Smoothing parameters:
r1		vertical (time) soothing factor. Usually 1<r1<=20
r2		horizontal soothing factor. Usually 1<r2<=20

Input Data:
in_traces	2-D array[nx][nt] of input traces

Output Data:
signal		2-D array[nx][nt] of extracted signal
noise		2-D array[nx][nt] of noise
*******************************************************************************
Note:
In this subroutine, the different iterations are basically the same except that
an updated noise model is used. The original data array is used in every 
iteration to compute the data histogram.
******************************************************************************/
{
	int ix,ip,it,iter;		/* loop counters */
	float dp;			/* moveout sampling interva (ms) */
	float rand=0.;			/* random number */
	float rel=0.;			/* minimum reliability for extraction */
	float **traces_taup;		/* array to store tau_p traces */
	float **temp_traces_taup;
	float **noise_taup;
	char *plotname="";

	if (plot<0||plot>3) err ("plot flag has to be 0<=plot<=3");
	
	/* number of required dips to avoid aliasing in tau-p domain */
	dp=(pmax1-pmin1)/(np-1);

	/* if requested, print processing information */
	if (verbose==1) warn("nx=%d nt=%d\n",nx,nt);

	/* allocate working space */
	temp_traces_taup=alloc2float(nt,np);
	traces_taup = alloc2float(nt,np);
	noise_taup = alloc2float(nt,np);

	/* compute forward global slant stack of input traces */
		forward_p_transform (nx, nt, dt, pmax1, pmin1, dp, depthref,
		60., 80., 4., 20., 400, 5, 1, 0, ninterp, 1, prewhite,
		interoff, offref, gopt, dx, fx, pmula, pmulb, in_traces,
		traces_taup);

	/* if requested, plot input data and its taup transform */
	if (plot==2||plot==3) { 
		plotname="input";
		plot_two_d (nx*nt, in_traces[0], plotname);
		plotname="input_taup";
		plot_two_d (nx*nt, traces_taup[0], plotname);
	}

	/* save input traces */
	for (ix=0; ix<nx; ix++) 
		for (it=0; it<nt; it++) {
			noise[ix][it]=in_traces[ix][it];
	}

	/* save taup traces */
	for (ip=0; ip<np; ip++) 
		for (it=0; it<nt; it++) {
			temp_traces_taup[ip][it]=traces_taup[ip][it];
	}

	/* do bed reflections extraction (one or two iterations) */
	for (iter=0; iter<niter; iter++) {

		/* compute seed for random generator */
		if (rgt==1) sranuni(seed);	
		else if (rgt==2) srannor(seed);

			/* compute array of noise (randomly reversed) traces */
			for (ix=0; ix<nx; ix++) {
			
			if (rgt==1) rand=2.0*franuni()-1.0;
			else if (rgt==2) rand=frannor();
			
				if (rand < 0.0) {	
					for (it=0; it<nt; it++) 
						noise[ix][it] *=-1.0;
			}
		}


		/* compute global slant stack of noise traces */
		forward_p_transform (nx, nt, dt, pmax1, pmin1, dp, depthref,
			60., 80., 4., 20., 400, 5, 1, 0, ninterp, 1, prewhite,
			interoff, offref, gopt, dx, fx, pmula, pmulb,
			noise, noise_taup);

		/* if requested, plot RRT and their tau-p transform */
		if (plot==2||plot==3) {		
			if (iter==0) plotname="rand1";
			else if (iter==1) plotname="rand2";
			plot_two_d (nx*nt, noise[0], plotname);
			if (iter==0) plotname="rtaup1";
			else if (iter==1) plotname="rtaup2";
			plot_two_d (nx*nt, noise_taup[0], plotname);
		}

		/* extract signal from diffractions and noise in taup domain */
		if (iter==0) rel=rel1;
		else if (iter==1) rel=rel2;
		extract_signal (verbose, norm, anenv, sts, plot, sditer, 
			nintlh, c, rel, np, nt, r1, r2, traces_taup,
			noise_taup);

		/* inverse slant stack to get extracted reflections in t-x dom*/
		inverse_p_transform (nx, nt, dt, pmax2, pmin2, dp, depthref,
			60., 80., interoff, offref, gopt, dx, fx, traces_taup,
			signal);

		/* if requested, output plot of extracted signal ad its taup*/
		if (plot==2||plot==3) { 
			if (iter==0) plotname="signal1";
			else if (iter==1) plotname="signal2";
			plot_two_d (nx*nt, signal[0], plotname);
			if (iter==0) plotname="staup1";
			else if (iter==1) plotname="staup2";
			plot_two_d (nx*nt, traces_taup[0], plotname);
		}

		/* if requested, scale filtered traces with respect to input */
		if (scl !=0) {
			scale_traces (scl, nx, nt, in_traces, signal); 

			/* if requested, output plot of scaled signal */
			if (plot==2) {
				if (iter==0) plotname="s1scaled";
				else if (iter==1) plotname="s2scaled";
				plot_two_d (nx*nt, signal[0], plotname);
			}
		}
	
		/* subtract signal from data to estimate noise */
			for (ix=0; ix<nx; ix++) 
					for (it=0; it<nt; it++) {
				if (scl !=-1) {
					noise[ix][it]=in_traces[ix][it]-
							signal[ix][it];
				}
			}

		/* recover the taup-transformed data */
			for (ip=0; ip<np; ip++) 
					for (it=0; it<nt; it++) {
				traces_taup[ip][it]=temp_traces_taup[ip][it];
		}

		/* if requested, output plot of extracted noise */
		if (plot==2||plot==3) {  
			if (iter==0) plotname="noise1";
			else if (iter==1) plotname="noise2";
			plot_two_d (nx*nt, noise[0], plotname);
		}
	}

	/* clean up */
	free2float(temp_traces_taup);
	free2float(traces_taup);
	free2float(noise_taup);
}

/******************************************************************************

	Subroutine to extract linear reflections from behind diffractions
				and 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)
/******************************************************************************
Input:

Flags:
anenv		=0 not to use analytic envelopes
		=1 to use signed analytic envelopes
		=2 to use positive analytic envelopes	
		=3 to use negative analytic envelopes
sts		=0 not to smooth estimaded signal amplitudes
		=1 do damped least squares temporal and spatial smoothing
plot		=0 don't make any plots
		=1 make plots of histograms, reliability and Esd

General Parameters:
sditer		number of steepest descent iterations for ps computation
nintlh		number of intervals (bins) per histogram
nx		number of horizontal samples (traces)
nt		number of vertical samples

Smoothing Parameters:
r1		vertical smoothing factor (if (sts=1). 1<=r1<=20
r2		horizontal smoothing factor (if (sts=1). 1<=r2<=20

Data:
traces		2-D array[nx][nt] of tau-p transformed input traces
rand_traces	2-D array[nx][nt] of randomly reversed tau-p transformed traces

Output:
traces		2-D array[nx][nt] of tau-p transformed extracted signal traces
rand_traces	2-D array[nx][nt] of tau-p transformed extracted noise traces
******************************************************************************/
{
	int ix;			/* loop counters */
	int rindex=0;		/* reference index for histograms */
	float spd,sps,spn;	/* auxiliary variables to test unit area */
	float spspn;		/* auxilliary variable			*/
	float *pd;	   	/* array of data probability values */
	float *ps;		/* array of signal probability values */
	float *pn;		/* array of noise probability values */
	float **anenv_traces;	/* 2-D array of analytic envelopes */
	float *reliability;	/* 2-D array of computed reliability indicator*/
	float *Esd;		/* 2-D array of expected values */
	float *xamps;		/* 1-d array of reference amplitudes */
	float *amps;		/* 1-d array of amplitudes in histograms */
	float *pspn;		/* 1-d array of ps*pn */
	char *plotname="";	/* pointer to name of output plot files */
	static int count;

	/* if requested, print processing information */
	if (verbose==1) {
	    warn("PARAMETERS: c=%g rel=%g r1=%g r2=%g\n",c,r,r1,r2); 
	    warn("FLAGS: anenv=%d plot=%d sts=%d\n",anenv,plot,sts);
	}

	/* allocate space */
	pd = alloc1float(nintlh);
	ps = alloc1float(nintlh);
	pn = alloc1float(nintlh);
	amps=alloc1float(nintlh);
	xamps=alloc1float(nintlh);
	pspn=alloc1float(nintlh);
	Esd = alloc1float(nintlh);
	reliability = alloc1float(nintlh);
	anenv_traces=alloc2float(nt,nx);

	/* if requested, compute analytic envelopes */
	if (anenv !=0) 
		compute_analytic_envelopes (anenv,nx,nt,traces,anenv_traces); 

	/* compute histogram parameters */
	compute_histogram_stuff (verbose, nx, nt, nintlh, traces, rand_traces,
		&rindex, amps, xamps);

	/* compute pd(x) and pn(x) via local histograms */
	make_histogram (nintlh, nt, nx, traces, amps, pd);
	make_histogram (nintlh, nt, nx, rand_traces, amps, pn);

	/* smooth noise histogram */
	smooth_histogram (nintlh, pn);

	/* compute ps(x) via minimization of cross-entropy */
	deconvolve_histograms (nintlh, rindex, sditer, pd, pn, ps);

	/* compute ps(x) convolved with pn(x) */
	conv (nintlh, -rindex, ps, nintlh, -rindex, pn, nintlh, -rindex, pspn);

	/* coumpute array of expected values */
		expected_value (nintlh, nx, nt, rindex, pspn, pn, ps, xamps, traces,
		Esd);

	/* compute reliability as a bounded convolution */
	compute_reliability (norm, nx, nt, nintlh, rindex, c, ps, pn, pspn,
		traces, reliability);

	/* zero out noisy samples depending on their reliability */
	zero_noisy_samples (anenv, sts, nx, nt, r1, r2, r, nintlh, amps,