sustaticrrs.c

su 的源代码库

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

/* SUSTATICRRS: $Revision: 1.21 $ ; $Date: 2003/06/09 16:17:07 $	*/

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


/*********************** self documentation ******************************/
char *sdoc[] = {
"									",
" SUSTATICRRS - Elevation STATIC corrections, apply corrections from	",
"	      headers or from a source and receiver statics file,	",
"	      includes application of Residual Refraction Statics	",
"									",
"     sustaticrrs <stdin >stdout  [optional parameters]	 		",
"									",
" Required parameters:							",
"	none								",
" Optional Parameters:							",
"	v0=v1 or user-defined	or from header, weathering velocity	",
"	v1=user-defined		or from header, subweathering velocity	",
"	hdrs=0			=1 to read statics from headers		",
" 				=2 to read statics from files		",
"	sign=1			=-1 to subtract statics from traces(up shift)",
" Options when hdrs=2:							",
"	sou_file=		input file for source statics (ms) 	",
"	rec_file=		input file for receiver statics (ms) 	",
"	ns=240 			number of sources 			",
"	nr=335 			number of receivers 			",
"	no=96 			number of offsets			",
"                                                                       ",
" Options when hdrs=3:                                                  ",
"       blvl_file=              base of the near-surface model file (sampled",
"                                  at CMP locations)                    ",
"       refr_file=              horizontal reference datum file (sampled at",
"                                  CMP locations)                       ",
"       nsamp=                  number of midpoints on line             ",
"       fx=                     first x location in velocity model      ",
"       dx=                     midpoint interval                       ",
"       V_r=                    replacement velocity                    ",
"       mx=                     number of velocity model samples in     ",
"                                  lateral direction                    ",
"       mz=                     number of velocity model samples in     ",
"                                  vertical direction                   ",
"       dzv=                    velocity model depth interval           ",
"       vfile=                  near-surface velocity model             ",
"                                                                       ",
" Options when hdrs=4:                                                  ",
"       nsamp=                  number of midpoints on line             ",
"       fx=                     first x location in velocity model      ", 
"       dx=                     midpoint interval                       ", 
"                                                                       ",
" Options when hdrs=5:                                                  ",
"       none                                                            ",
"									",
" Notes:								",
" For hdrs=1, statics calculation is not performed, statics correction  ",
" is applied to the data by reading statics (in ms) from the header.	",
"									",
" For hdrs=0, field statics are calculated, and				",
" 	input field sut is assumed measured in ms.			",
" 	output field sstat = 10^scalel*(sdel - selev + sdepth)/swevel	",
" 	output field gstat = sstat - sut/1000.				",
" 	output field tstat = sstat + gstat + 10^scalel*(selev - gelev)/wevel",
"									",
" For hdrs=2, statics are surface consistently obtained from the 	",
" statics files. The geometry should be regular.			",
" The source- and receiver-statics files should be unformated C binary 	",
" floats and contain the statics (in ms) as a function of surface location.",
"									",
" For hdrs=3, residual refraction statics and average refraction statics",
" are computed.  For hdrs=4, residual refraction statics are applied,   ",
" and for hdrs=5, average refraction statics are applied (Cox, 1999).   ",
" These three options are coupled in many data processing sequences:    ",
" before stack residual and average refraction statics are computed but ",
" only residual refractions statics are applied, and after stack average",
" refraction statics are applied.  Refraction statics are often split   ",
" like this to avoid biasing stacking velocities.  The files blvl_file  ",
" and refr_file are the base of the velocity model defined in vfile and ",
" the final reference datum, as described by Cox (1999), respectively.  ",
" Residual refraction statics are stored in the header field gstat, and ",
" the average statics are stored in the header field tstat.  V_r is the ",
" replacement velocity as described by Cox (1999).  The velocity file,  ",
" vfile, is designed to work with a horizontal upper surface defined in ",
" refr_file.  If the survey has irregular topography, the horizontal    ",
" upper surface should be above the highest topographic point on the    ",
" line, and the velocity between this horizontal surface and topography ",
" should be some very large value, such as 999999999, so that the       ",
" traveltimes through that region are inconsequential.                  ",
NULL};

/* Credits:
 *	CWP: Jamie Burns
 *
 *	CWP: Modified by Mohammed Alfaraj, 11/10/1992, for reading
 *	     statics from headers and including sign (+-) option
 *
 *      CWP: Modified by Timo Tjan, 29 June 1995, to include input of
 *           source and receiver statics from files. 
 *
 *      CWP: Modified by Chris Robinson, 11/2000, to include the splitting
 *           of refraction statics into residuals and averages
 *
 * Trace header fields accessed:  ns, dt, delrt, gelev, selev,
 *	sdepth, gdel, sdel, swevel, sut, scalel
 * Trace header fields modified:  sstat, gstat, tstat
 *
 * References:
 *
 * Cox, M., 1999, Static corrections for seismic reflection surveys:
 *    Soc. Expl. Geophys.
 */

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


segy intrace, outtrace;

int
main(int argc, char **argv)
{
	int nt;		/* number of samples on output trace	*/
	float dt;	/* sample rate on outpu trace		*/
	int itime;	/* counter          			*/
	float tmin;	/* first time sample on output trace	*/
	float tsd=0.0;	/* time to move source to datum         */
	float trd=0.0;	/* time to move 0 offset receiver       */
	float v0;	/* weathering velocity			*/
	float v1;	/* subweathering velocity		*/
	int hdrs; 	/* flag to read statics from headers	*/ 
	float *t;	/* array of output times		*/
	float tstat=0.0; /* total (source and receiver) statics	*/
	int sign;	/* to add (+) or subtract (-) statics	*/
	int no;		/* number of offsets per shot 		*/
	int io;		/* offset counter 			*/
	int is;		/* source counter 			*/
	int ir;		/* receiver counter 			*/
	int ns;		/* number of sources = number of source statics */
	int nr;		/* number of receiver = number of rec. statics	*/
        int nsamp;      /* number of samples for each surface           */
        int i,j,res=0,bulk=0,mx,mz;
        float fx;      /* first sample position of surfaces            */
        float dx;
        int ic_c,is_c,ir_c;
        float z_blvl_s,z_refr_s;
        float z_blvl_r,z_refr_r;
	float *sou_statics=NULL;	/* array of source statics	*/
	float *rec_statics=NULL;	/* array of receiver statics	*/
        float *blvl_arr=NULL;           /* array of base LVL values     */
        float *refr_arr=NULL;           /* array of reference datum values */
        float *tcdp_stat=NULL, *tcdp_arr=NULL;
        float *num=NULL, *num_arr=NULL, **vel=NULL;
        float dzv;
        float V_r;
	FILE *fps, *fpr;       	/* file pointers for statics input 	*/
        FILE *fpb=NULL, *fpf=NULL;   /* file pointers for input surfaces     */
        FILE *g1=NULL, *g2=NULL, *g3=NULL, *g4=NULL;
        FILE *vfp;
	cwp_String sou_file, rec_file; /* statics filenames 		*/
        cwp_String blvl_file; /* surface filenames                      */
        cwp_String refr_file;            /*surface filename             */
        cwp_String vfile;

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

	/* Get information from first trace */
	if (!gettr(&intrace)) err("can't get first trace");
	nt   = intrace.ns;
	tmin = intrace.delrt/1000.0;
	dt   = ((double) intrace.dt)/1000000.0;
	
	/* Get parameters */
	if (!getparfloat("v1", &v1))          v1 = (float) intrace.swevel;
	if (!getparfloat("v0", &v0))
                v0 = (float) ((intrace.wevel) ? intrace.wevel : v1);
	if (!getparint("hdrs", &hdrs))        hdrs = 0;
	if (!getparint("sign", &sign))        sign = 1;

	/* Allocate vector of output times */
	t = ealloc1float(nt);

	/* reading source and receiver statics from files */
	if (hdrs == 2){

		/* getpar statics file related parameters */
		if (!getparint("ns", &ns))        ns = 240;
		if (!getparint("nr", &nr))        nr = 335;
		if (!getparint("no", &no))        no = 96;

		/* getpar statics file names */
        	getparstring("sou_file",&sou_file);
        	getparstring("rec_file",&rec_file);

		/* allocate space */
		rec_statics = alloc1float(nr);
        	sou_statics = alloc1float(ns);

		/* open and read from receiver statics file */
        	if((fpr=efopen(rec_file,"rb"))==NULL)
                	err("cannot open stat_file=%s\n",rec_file);
        	efread(rec_statics, sizeof(float),nr,fpr);