sunmo.c

su 的源代码库

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

/* SUNMO: $Revision: 1.24 $ ; $Date: 2006/10/31 22:25:06 $		*/
 
#include "su.h"
#include "segy.h"

/*********************** self documentation ******************************/
char *sdoc[] = {
"									     ",
" SUNMO - NMO for an arbitrary velocity function of time and CDP	     ",
"									     ",
"  sunmo <stdin >stdout [optional parameters]				     ",
"									     ",
" Optional Parameters:							     ",
" tnmo=0		NMO times corresponding to velocities in vnmo	     ",
" vnmo=1500		NMO velocities corresponding to times in tnmo	     ",
" anis1=0		two anisotropy coefficients making up quartic term   ",
" anis2=0		in traveltime curve, corresponding to times in tnmo  ",
" cdp=			CDPs for which vnmo & tnmo are specified (see Notes) ",
" smute=1.5		samples with NMO stretch exceeding smute are zeroed  ",
" lmute=25		length (in samples) of linear ramp for stretch mute  ",
" sscale=1		=1 to divide output samples by NMO stretch factor    ",
" invert=0		=1 to perform (approximate) inverse NMO	             ",
" ixoffset=0		do not consider cross-line offset                    ",
"                       =1 read cross-line offset from trace header          ",
" upward=0		=1 to scan upward to find first sample to kill       ",
"									     ",
" Notes:								     ",
"									     ",
" For constant-velocity NMO, specify only one vnmo=constant and omit tnmo.   ",
"									     ",
" The anisotropy coefficients anis1, anis2 permit non-hyperbolicity due	     ",
" to layering, mode conversion, or anisotropy. Default is isotropic NMO.     ",
"									     ",
" For NMO with a velocity function of time only, specify the arrays	     ",
"	   vnmo=v1,v2,... tnmo=t1,t2,...				     ",
" where v1 is the velocity at time t1, v2 is the velocity at time t2, ...    ",
" The times specified in the tnmo array must be monotonically increasing.    ",
" Linear interpolation and constant extrapolation of the specified velocities",
" is used to compute the velocities at times not specified.		     ",
" The same holds for the anisotropy coefficients as a function of time only. ",
"									     ",
" For NMO with a velocity function of time and CDP, specify the array	     ",
"	   cdp=cdp1,cdp2,...						     ",
" and, for each CDP specified, specify the vnmo and tnmo arrays as described ",
" above. The first (vnmo,tnmo) pair corresponds to the first cdp, and so on. ",
" Linear interpolation and constant extrapolation of 1/velocity^2 is used    ",
" to compute velocities at CDPs not specified.				     ",
" The same holds for the anisotropy coefficients as a function of time and   ",
" CDP.									     ",
"									     ",
" Moveout is defined by							     ",
"									     ",
"   1		 anis1							     ",
"  --- x^2 + ------------- x^4.						     ",
"  v^2	     1 + anis2 x^2						     ",
"									     ",
" Note: In general, the user should set the cdp parameter.  The default is   ",
"	to use tr.cdp from the first trace and assume only one cdp.          ",
" Caveat:								     ",
" Nmo cannot handle negative moveout as in triplication caused by	     ",
" anisotropy. But negative moveout happens necessarily for negative anis1 at ",
" sufficiently large offsets. Then the error-negative moveout- is printed.   ",
" Check anis1. An error (anis2 too small) is also printed if the	     ",
" denominator of the quartic term becomes negative. Check anis2. These errors",
" are prompted even if they occur in traces which would not survive the	     ",
" NMO-stretch threshold. Chop off enough far-offset traces (e.g. with suwind)",
" if anis1, anis2 are fine for near-offset traces.			     ",
"									     ",
" NMO interpolation error is less than 1% for frequencies less than 60% of   ",
" the Nyquist frequency.						     ",
"									     ",
" Exact inverse NMO is impossible, particularly for early times at large     ",
" offsets and for frequencies near Nyquist with large interpolation errors.  ",
NULL};

/* Credits:
 *	SEP: Shuki, Chuck Sword
 *	CWP: Shuki, Jack, Dave Hale, Bjoern Rommel
 *      Modified: 08/08/98 - Carlos E. Theodoro - option for lateral offset
 *      Modified: 07/11/02 - Sang-yong Suh -
 *          added "upward" option to handle decreasing velocity function.
 *
 * Technical Reference:
 *	The Common Depth Point Stack
 *	William A. Schneider
 *	Proc. IEEE, v. 72, n. 10, p. 1238-1254
 *	1984
 *
 * Trace header fields accessed: ns, dt, delrt, offset, cdp, sy
 */
/**************** end self doc *******************************************/

static void interpovv (int nt, int ncdp, float *cdp, 
	float **ovv, float **oa1, float **oa2, float cdpt, 
	float *ovvt, float *oa1t, float *oa2t);

segy tr;

int
main(int argc, char **argv)
{
	int nt;		/* number of time samples per trace */
	float dt;	/* time sampling interval */
	float ft;	/* time of first sample */
	int it;		/* time sample index */
	int ncdp;	/* number of cdps specified */
	float *cdp;	/* array[ncdp] of cdps */
	int icdp;	/* index into cdp array */
	int jcdp;	/* index into cdp array */
	int nvnmo;	/* number of vnmos specified */
	float *vnmo;	/* array[nvnmo] of vnmos */
	int ntnmo;	/* number of tnmos specified */
	float *tnmo;	/* array[ntnmo] of tnmos */
	float **ovv;	/* array[ncdp][nt] of sloth (1/velocity^2) functions */
	float *ovvt;	/* array[nt] of sloth for a particular trace */
	int nanis1;	/* number of anis1's specified */
	int nanis2;	/* number of anis2's specified */
	float *anis1;	/* array[nanis1] of anis1's */
	float *anis2;	/* array[nanis2] of anis2's */
	float **oa1;	/* array[ncdp][nt] of anis1 functions */
	float **oa2;	/* array[ncdp][nt] of anis2 functions */
	float *oa1t;	/* array[nt] of anis1 for a particular trace */
	float *oa2t;	/* array[nt] of anis2 for a particular trace */
	float smute;	/* zero samples with NMO stretch exceeding smute */
	float osmute;	/* 1/smute */
	int lmute;	/* length in samples of linear ramp for mute */
	int itmute=0;	/* zero samples with indices less than itmute */
	int sscale;	/* if non-zero, apply NMO stretch scaling */
	int invert;	/* if non-zero, do inverse NMO */
	float sy;	/* cross-line offset component */
	int ixoffset;	/* indes for cross-line offset component */
	long oldoffset;	/* offset of previous trace */
	long oldcdp;	/* cdp of previous trace */
	int newsloth;	/* if non-zero, new sloth function was computed */
	float tn;	/* NMO time (time after NMO correction) */
	float v;	/* velocity */
	float *qtn;	/* NMO-corrected trace q(tn) */
	float *ttn;	/* time t(tn) for NMO */
	float *atn;	/* amplitude a(tn) for NMO */
	float *qt;	/* inverse NMO-corrected trace q(t) */
	float *tnt;	/* time tn(t) for inverse NMO */
	float *at;	/* amplitude a(t) for inverse NMO */
	float acdp;	/* temporary used to sort cdp array */
	float *aovv;	/* temporary used to sort ovv array */
	float *aoa1;	/* temporary used to sort oa1 array */
	float *aoa2;	/* temporary used to sort oa2 array */
	float temp;	/* temporary float */
	float tsq;	/* temporary float */
	int i;		/* index used in loop */
	int upward;	/* scans upward if it's nonzero. */

	/* hook up getpar */
	initargs(argc, argv);
	requestdoc(1);

	/* get information from the first header */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = ((double) tr.dt)/1000000.0;
	ft = tr.delrt/1000.0;
	sy  = tr.sy;

	/* get velocity functions, linearly interpolated in time */
	ncdp = countparval("cdp");
	if (ncdp>0) {
		if (countparname("vnmo")!=ncdp)
			err("a vnmo array must be specified for each cdp");
		if (countparname("tnmo")!=ncdp)
			err("a tnmo array must be specified for each cdp");
		if (countparname("anis1")!=ncdp &&
		    countparname("anis1")!=0)
			err("an anis1 array must be specified for each cdp, "
			    "or omitted at all");
		if (countparname("anis2")!=ncdp &&
		    countparname("anis2")!=0)
			err("an anis2 array must be specified for each cdp, "
			    "or omitted at all");
	} else {
		ncdp = 1;
		if (countparname("vnmo")>1)
			err("only one (or no) vnmo array must be specified");
		if (countparname("tnmo")>1)
			err("only one (or no) tnmo array must be specified");
		if (countparname("anis1")>1)
			err("only one (or no) anis1 array must be specified");
		if (countparname("anis2")>1)
			err("only one (or no) anis2 array must be specified");
	}
	cdp = ealloc1float(ncdp);
	if (!getparfloat("cdp",cdp)) cdp[0] = tr.cdp;
	ovv = ealloc2float(nt,ncdp);
	oa1 = ealloc2float(nt,ncdp);
	oa2 = ealloc2float(nt,ncdp);
	for (icdp=0; icdp<ncdp; ++icdp) {
		nvnmo = countnparval(icdp+1,"vnmo");
		ntnmo = countnparval(icdp+1,"tnmo");
		nanis1 = countnparval(icdp+1,"anis1");
		nanis2 = countnparval(icdp+1,"anis2");
		if (nvnmo!=ntnmo && !(ncdp==1 && nvnmo==1 && ntnmo==0))
			err("number of vnmo and tnmo values must be equal");
		if (nanis1!=nvnmo && nanis1 != 0)
			err("number of vnmo and anis1 values must be equal");
		if (nanis2!=nvnmo && nanis2 != 0)
			err("number of vnmo and anis2 values must be equal");
		if (nvnmo==0) nvnmo = 1;
		if (ntnmo==0) ntnmo = nvnmo;
		if (nanis1==0) nanis1 = nvnmo;
		if (nanis2==0) nanis2 = nvnmo;
		/* equal numbers of parameters vnmo, tnmo, anis1, anis2 */
		vnmo = ealloc1float(nvnmo);
		tnmo = ealloc1float(nvnmo);
		anis1 = ealloc1float(nvnmo);
		anis2 = ealloc1float(nvnmo);
		if (!getnparfloat(icdp+1,"vnmo",vnmo)) vnmo[0] = 1500.0;
		if (!getnparfloat(icdp+1,"tnmo",tnmo)) tnmo[0] = 0.0;
		if (!getnparfloat(icdp+1,"anis1",anis1)) 
			for (i=0; i<nvnmo; i++) anis1[i] = 0.0;
		if (!getnparfloat(icdp+1,"anis2",anis2))
			for (i=0; i<nvnmo; i++) anis2[i] = 0.0;
		for (it=1; it<ntnmo; ++it)
			if (tnmo[it]<=tnmo[it-1])
				err("tnmo values must increase monotonically");
		for (it=0,tn=ft; it<nt; ++it,tn+=dt) {