suinterp.c

su 的源代码库

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

/* SUINTERP: $Revision: 1.13 $ ; $Date: 2003/06/09 16:17:07 $                */

#include "su.h"
#include "segy.h"
#include "header.h"
#include "VND.h"
#include <signal.h>

/*********************** self documentation **********************/
char *sdoc[] = { 
"									",
" SUINTERP - interpolate traces using automatic event picking		",
"									",
"           suinterp < stdin > stdout					",
"									",
" ninterp=1    number of traces to output between each pair of input traces",
" nxmax=500    maximum number of input traces				",
" freq1=4.     starting corner frequency of unaliased range		",
" freq2=20.    ending corner frequency of unaliased range		",
" deriv=0      =1 means take vertical derivative on pick section        ",
"              (useful if interpolating velocities instead of seismic)  ",
" linear=0     =0 means use 8 point sinc temporal interpolation         ",
"              =1 means use linear temporal interpolation               ",
"              (useful if interpolating velocities instead of seismic)  ",
" lent=5       number of time samples to smooth for dip estimate	",
" lenx=1       number of traces to smooth for dip estimate		",
" lagc=400     number of ms agc for dip estimate			",
" xopt=0       0 compute spatial derivative via FFT			",
"                 (assumes input traces regularly spaced and relatively	",
"                  noise-free)						",
"              1 compute spatial derivative via differences		",
"                 (will work on irregulary spaced data)			",
" iopt=0     0 = interpolate",
"            1 = output low-pass model: useful for QC if interpolator failing",
"            2 = output dip picks in units of samples/trace		",
"									",
" verbose=0	verbose = 1 echoes information				",
"									",
" 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()		",
" 									",
" Notes:								",
" This program outputs 'ninterp' interpolated traces between each pair of",
" input traces.  The values for lagc, freq1, and freq2 are only used for",
" event tracking. The output data will be full bandwidth with no agc.  The",
" default parameters typically will do a satisfactory job of interpolation",
" for dips up to about 12 ms/trace.  Using a larger value for freq2 causes",
" the algorithm to do a better job on the shallow dips, but to fail on the",
" steep dips.  Only one dip is assumed at each time sample between each pair",
" of input traces.							",
" 									",
" The key assumption used here is that the low frequency data are unaliased",
" and can be used for event tracking. Those dip picks are used to interpolate",
" the original full-bandwidth data, giving some measure of interpolation",
" at higher frequencies which otherwise would be aliased.  Using iopt equal",
" to 1 allows you to visually check whether the low-pass picking model is",
" aliased.								",
" 									",
" Trace headers for interpolated traces are not updated correctly.	",
" The output header for an interpolated traces equals that for the preceding",
" trace in the original input data.  The original input traces are passed",
" through this module without modification.				",
"									",
" The place this code is most likely to fail is on the first breaks.	",
"									",
" Example run:    suplane | suinterp | suxwigb &			",
"									",
NULL}; 

/*
 * Credit: John Anderson (visiting scholar from Mobil) July 1994
 *
 * Trace header fields accessed: ns, dt
 */

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

void runav(int n,int len,float *a,float *b);
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,
		int deriv, int linear);

static void closefiles(void);

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

segy tr;	/* Input and output trace data of length nt */

int
main (int argc, char **argv)
{
	int	nx,nt,ix,nxmax,ninterp,k;
	int	lent,lenx,lagc,xopt,iopt,deriv,linear;
	float	dt,freq1,freq2;
	VND	*vndorig,*vndinterp;
	char	*file;
	int verbose;	/* flag for echoing info			*/
	char *tmpdir;	/* directory path for tmp files			*/
	cwp_Bool istmpdir=cwp_false;/* true for user-given path		*/

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

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

	if (!getparint("lagc",&lagc))lagc=400;
	if (!getparfloat("freq1",&freq1)) freq1=4.;
	if (!getparfloat("freq2",&freq2)) freq2=20.;
	if (!getparint("lent",&lent)) lent=5;
	if (!getparint("lenx",&lenx)) lenx=1;
	if (!getparint("nxmax",&nxmax)) nxmax=500;
	if (!getparint("xopt",&xopt)) xopt=0;
	if (!getparint("ninterp",&ninterp)) ninterp=1;
	if (!getparint("iopt",&iopt)) iopt=0;
	if (!getparint("deriv",&deriv)) deriv=0;
	if (!getparint("linear",&linear)) linear=0;
	if (!getparint("verbose", &verbose)) verbose = 0;

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

	file=VNDtempname("suinterp");
	vndorig = V2Dop(2,1000000,sizeof(float),
			file,nt,nxmax);
	VNDfree(file,"file");
	file=VNDtempname("suinterp");
	vndinterp = V2Dop(2,1000000,sizeof(float),
			file,
			nt,1+(nxmax-1)*(ninterp+1));
	VNDfree(file,"file");

	if (STREQ(tmpdir,"")) {
		headerfp = etmpfile();
		if (verbose) warn("using tmpfile() call");
	} else { /* user-supplied tmpdir */
		char directory[BUFSIZ];
		strcpy(directory, tmpdir);
		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);
		headerfp = efopen(headerfile, "w+");
      		istmpdir=cwp_true;		
		if (verbose) warn("putting temporary files in %s", directory);
	}

	/* Main loop for saving input traces */
	nx=0;
	do {
		V2Dw0(vndorig,nx,(char *)tr.data,1);
		efwrite(&tr,HDRBYTES,1,headerfp);
		nx++;
		if(nx>=nxmax) break;
	} while(gettr(&tr));

	jea_xinterpolate(vndorig,vndinterp,ninterp,nt,nx,freq1,freq2,
			lagc,lent,lenx,xopt,dt,iopt,deriv,linear);

	/* loop outputting results */
	if(iopt!=0) ninterp=0;
	efseeko(headerfp, (off_t) 0,SEEK_SET);
	for(ix=0;ix<nx-1;ix++) {
		efread(&tr,HDRBYTES,1,headerfp);
		for(k=0;k<=ninterp;k++) {
			V2Dr0(vndinterp,k+ix*(ninterp+1),(char *)tr.data,18);
			puttr(&tr);
		}
	} 
	efread(&tr,HDRBYTES,1,headerfp);
	V2Dr0(vndinterp,(nx-1)*(ninterp+1),(char *)tr.data,18);
	puttr(&tr);

	/* close files and return */
	VNDcl(vndorig,1);
	VNDcl(vndinterp,1);
	if(VNDtotalmem()!=0) {
		fprintf(stderr,"total VND memory at end = %ld\n",
		VNDtotalmem());
	}
	efclose(headerfp);
	if (istmpdir) eremove(headerfile);

	return(CWP_Exit());
}

void jea_xinterpolate(VND *vndorig, VND *vndinterp, int ninterp, 
		int nt, int n2, float freq1, float freq2, int lagc, 
		int lent, int len2, int xopt, float dt, int iopt,
		int deriv, int linear)
/*******************************************************************
interpolate input data in space placing "ninterp" synthetic traces 
between each pair of original input traces
******************************************************************
Function parameters:

VND *vndorig		VND file with input data
VND *vndinterp		VND file with output original plus interpolated data
int ninterp		number of traces to interpolate between each
			input trace
int nt			number of time samples
int n2			number of input traces
float freq1		low-end frequency in Hz for picking
						(good default: 3 Hz)
float freq2		high-end frequency in Hz for picking
						(good default: 20 Hz)
int lagc		length of AGC operator for picking
						(good default: 400 ms)
int lent		length of time smoother in samples for picker
                        (good default: 5 samples)
int len2		length of space smoother in samples for picker
                        (good default: 1 sample)
int xopt		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)
float dt		sample rate in sec
int iopt		0 = interpolate: output 1+(n2-1)*(1+ninterp) traces
                            with ninterp traces between each pair of
			    input traces
			1 = compute low-pass model: output n2 traces
                            on original trace locations -- This is typically
                            used for Quality Control if the interpolator
                            is failing for any reason
			2 = compute dip picks in units of samples/trace: 

                            output n2 traces on original trace locations
int deriv		0 = default
			1 = take vertical derivative for picking section
int linear		0 = default use 8 point sinc interpolators
			1 = use linear interpolation

This routine outputs 'ninterp' interpolated traces between each pair of 
input traces.  The values for lagc, freq1, and freq2 are only used for
event tracking. The output data will be full bandwidth with no agc.  The 
suggested default parameters typically will do a satisfactory job of 
interpolation for dips up to about 12 ms/trace.  Using a larger value for 
freq2 causes the algorithm to do a better job on the shallow dips, but to 
fail on the steep dips.  Only one dip is assumed at each time sample between 
each pair of input traces.  The original input traces are passed through
this routine without modification.

The key assumption used here is that the low frequency data are unaliased
and can be used for event tracking.  Those dip picks are used to interpolate
the original full-bandwidth data, giving some measure of interpolation
at higher frequencies which otherwise would be aliased.  Using iopt equal
to 1 allows you to visually check whether the low-pass picking model
is aliased.
If you can't visually pick dips correctly on the low-pass picking 
model, this computer routine will fail.

The place this code is most likely to fail is on the first breaks.

This routine assumes that the input and output files hav been allocated in
the calling routine as

vndorig = V2Dop(2,1000000,sizeof(float),VNDtempname("suinterp"),nt,n2max);
vndinterp = V2Dop(2,1000000,sizeof(float),VNDtempname("suinterp"),
			nt,1+(n2max-1)*(ninterp+1));

where n2max is the maximum number of input traces and nt is the number 
of time samples.
*******************************************************************