sudatumk2ds.c

su 的源代码库

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

/* SUDATUMK2DS: $Revision: 1.5 $ ; $Date: 2006/11/07 22:58:42 $	*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
"SUDATUMK2DS - Kirchhoff datuming of sources for 2D prestack data	",
" 		(input data are receiver gathers) 			",
" 									",
"    sudatumk2ds  infile=  outfile=  [parameters] 			",
"									",
" Required parameters:							",
" infile=stdin		file for input seismic traces			",
" outfile=stdout	file for common offset migration output  	",
" ttfile=		file for input traveltime tables		",
"   The following 9 parameters describe traveltime tables:		",
" fzt=			first depth sample in traveltime table		",
" nzt= 			number of depth samples in traveltime table	",
" dzt=			depth interval in traveltime table		",
" fxt=			first lateral sample in traveltime table	",
" nxt=			number of lateral samples in traveltime table	",
" dxt=			lateral interval in traveltime table		",
" fs= 			x-coordinate of first source			",
" ns= 			number of sources				",
" ds= 			x-coordinate increment of sources		",
"                                                                       ",
" fxi=                   x-coordinate of the first surface location      ",
" dxi=                   horizontal spacing on surface                   ",
" nxi=                   number of input surface locations               ",
" sgn=                   Sign of the datuming process (up=-1 or down=1)  ",
"									",
" Optional Parameters:							",
" dt= or from header (dt) 	time sampling interval of input data	",
" ft= or from header (ft) 	first time sample of input data		",
" surf=\"0,0;99999,0\"  The first surface defined the recording surface ",
" surf=\"0,0;99999,0\"  and the second one, the new datum.              ",
"                       \"x1,z1;x2,z2;x3,z3;...\"                       ",
" fzo=fzt		z-coordinate of first point in output trace 	",
" dzo=0.2*dzt		vertical spacing of output trace 		",
" nzo=5*(nzt-1)+1 	number of points in output trace		",	
" fxso=fxt		x-coordinate of first shot	 		",
" dxso=0.5*dxt		shot horizontal spacing		 		",
" nxso=2*(nxt-1)+1  	number of shots 				",
" fxgo=fxt		x-coordinate of first receiver			",
" exgo=fxgo+(nxgo-1)*dxgo	x-coordinate of the last receiver	",
" dxgo=0.5*dxt		receiver horizontal spacing			",
" nxgo=nxso		number of receivers per shot			",
" fmax=0.25/dt		frequency-highcut for input traces		",
" offmax=99999		maximum absolute offset allowed in migration 	",
" aperx=nxt*dxt/2  	migration lateral aperature 			",
" angmax=60		migration angle aperature from vertical 	",
" v0=1500(m/s)		reference velocity value at surface		",
" dvz=0.0  		reference velocity vertical gradient		",
" antiali=1             Antialiase filter (no-filter = 0)               ",
" jpfile=stderr		job print file name 				",
" mtr=100  		print verbal information at every mtr traces	",
" ntr=100000		maximum number of input traces to be migrated	",
"									",
" Notes:								",
" 1. Traveltime tables were generated by program rayt2d (or other ones)	",
"    on relatively coarse grids, with dimension ns*nxt*nzt. In the	",
"    migration process, traveltimes are interpolated into shot/gephone 	",
"    positions and output grids.					",
" 2. Input traces must be SU format and organized in common shot gathers",
" 3. If the offset value of an input trace is not in the offset array 	",
"    of output, the nearest one in the array is chosen. 		",
" 4. Amplitudes are computed using the reference velocity profile, v(z),",
"    specified by the parameters v0= and dvz=.				",
" 5. Input traces must specify source and receiver positions via the header",
"    fields tr.sx and tr.gx. Offset is computed automatically.		",
"									",
NULL};
/*
 * Author:  Trino Salinas, 05/01/96,  Colorado School of Mines
 *
 * This code is based on sukzmig2d.c written by Zhenyue Liu, 03/01/95.
 * Subroutines from Dave Hale's modeling library were adapted in
 * this code to define topography using cubic splines.
 *
 * This code implements a Kirchhoff extraplolation operator that allows to
 * transfer data from one reference surface to another.  The formula used in
 * this application is a far field approximation of the Berryhill's original
 * formula (Berryhill, 1979).  This equation is the result of a stationary
 * phase analysis to get an analog asymptotic expansion for the two-and-one
 * half dimensional extrapolation formula (Bleistein, 1984).
 *
 * The extrapolation formula permits the downward continuation of upgoing
 * waves  and  upward  continuation  of  downgoing waves.  For upward conti-
 * nuation of upgoing waves and downward continuation of downgoing waves,
 * the conjugate transpose of the equation is used (Bevc, 1993).
 *
 * References :
 *
 * Berryhill, J.R., 1979, Wave equation datuming: Geophysics,
 *   44, 1329--1344.
 *
 * _______________, 1984, Wave equation datuming before stack
 *   (short note) : Geophysics, 49, 2064--2067.
 *
 * Bevc, D., 1993, Data parallel wave equation datuming with
 *   irregular acquisition topography :  63rd Ann. Internat.
 *   Mtg., SEG, Expanded Abstracts, 197--200.
 *
 * Bleistein, N., 1984, Mathematical methods for wave phenomena,
 *   Academic Press Inc. (Harcourt Brace Jovanovich Publishers),
 *   New York.
 *
**************** end self doc ***********************************/
/* TYPEDEFS */
typedef struct SurfaceSegmentStruct {
        float x;        /* x coordinate of segment midpoint */
        float z;        /* z coordinate of segment midpoint */
        float s;        /* x component of unit-normal-vector */
        float c;        /* z component of unit-normal-vector */
} SurfaceSegment;
typedef struct SurfaceStruct {
        int ns;               /* number of surface segments */
        float ds;             /* segment length */
        SurfaceSegment *ss;   /* array[ns] of surface segments */
} Surface;

/* Prototype */
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,float **tb,
	     float **t,float x0);
  void interpx(int nxt,float fxt,float dxt,int nx,float fx,float dx,int nzt,
	     float **tt,float **t);
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t);
  void timeb(int nr,int nz,float dr,float dz,float fz,float sz,float a,
             float v0, float **t,float **p,float **sig,float **ang);
  void dat2d(float *trace,int nt,float ft,float dt,float sx,float gx,
             float **dat,float aperx,int nx,float fx,float dx,float nz,
             float fz,float dz,int mtmax,float xm,float fmax,
             int nxi,float fxi,float dxi,float angmax,float **tb,float **pb,
             float **angb,int nr,float **tsum,int nzt,
             float fzt,float dzt,int nxt,float fxt,float dxt,int antiali,
             int sgn,float **szif,float nangl,float **sigb);
  void decodeSurfaces(int sgn, int *nrPtr, int **nxzPtr, float ***xPtr,
	     float ***zPtr);
  int decodeSurface(char *string, int sgn, int *nxzPtr, float **xPtr,
	     float **zPtr);
  void makesurf(float dsmax,int nr,int *nu,float **xu,float **zu,
             Surface **r);
  void zcoorSurfaces(float fx,float dx,int nx,float fxs,float dxs,int nxs,
             Surface *srf,float **szif,float *sz,float *nangl);

/* segy trace */
segy tr, tro;

int
main (int argc, char **argv)
{
	int   nt,nzt,nxt,nzo,nxso,nxgo,ns,nr,is,io,ixo,it,antiali,
	      ntr,jtr,ktr,mtr,mtmax,nxgt,*ng,nsrf,*nxzsrf,nxi,
	      sgn;
	off_t  nseek;
	float ft,fzt,fxt,fzo,fxso,fxgo,fs,dt,dzt,dxt,dzo,dxso,
	      dxgo,ds,ext,ezt,ezo,exso,exgo,es,s,scal,fxin,as,res;
	float v0,dvz,fmax,angmax,offmax,rmax,aperx,sx,gx,**xsrf,**zsrf,
	      v00,nangls,dxi,fxi;
        float ***datg,***ttab,***tb,***pb,***angb,**tsum,**tt,**szif,
                *nangl,*sz,**tbs,**pbs,**angbs,***sigb,**sigbs;
        Surface *srf;
	
	char  *infile="stdin",*outfile="stdout",*ttfile,*jpfile;
	FILE  *infp,*outfp,*ttfp,*jpfp,*hdrfp;


	/* hook up getpar to handle the parameters */
	initargs(argc, argv);
	requestdoc(1);

	/* open input and output files	*/
	if(!getparstring("infile",&infile)) {
		infp = stdin;
	} else  
		if ((infp=fopen(infile,"r"))==NULL)
			err("cannot open infile=%s\n",infile);
	if( !getparstring("outfile",&outfile)) {
		outfp = stdout;
	} else  {
		outfp = fopen(outfile,"w");
	}

	efseeko(infp,(off_t) 0,SEEK_CUR);
	efseeko(outfp,(off_t) 0,SEEK_CUR);

	if( !getparstring("ttfile",&ttfile))
		err("must specify ttfile!\n");
	if ((ttfp=fopen(ttfile,"r"))==NULL)
		err("cannot open ttfile=%s\n",ttfile);
	if( !getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	} else  
		jpfp = fopen(jpfile,"w");

	/* get information from the first header */
	if (!fgettr(infp,&tr)) err("can't get first trace");
	nt = tr.ns;
	if (!getparfloat("dt",&dt)) dt = tr.dt/1000000.0; 
	if (dt<0.0000001) err("dt must be positive!\n");
	if (!getparfloat("ft",&ft)) ft = (float)tr.delrt/1000.; 
	
	/* get traveltime table parameters	*/
	if (!getparint("nxt",&nxt)) err("must specify nxt!\n");
	if (!getparfloat("fxt",&fxt)) err("must specify fxt!\n");
	if (!getparfloat("dxt",&dxt)) err("must specify dxt!\n");
	if (!getparint("nzt",&nzt)) err("must specify nzt!\n");
	if (!getparfloat("fzt",&fzt)) err("must specify fzt!\n");
	if (!getparfloat("dzt",&dzt)) err("must specify dzt!\n");
	if (!getparint("ns",&ns)) err("must specify ns!\n");
	if (!getparfloat("fs",&fs)) err("must specify fs!\n");
	if (!getparfloat("ds",&ds)) err("must specify ds!\n");
	ext = fxt+(nxt-1)*dxt;
	ezt = fzt+(nzt-1)*dzt;
	es = fs+(ns-1)*ds;

        if (!getparint("nxi",&nxi)) err("must specify nxi!\n");
        if (!getparfloat("fxi",&fxi)) err("must specify fxi!\n");
        if (!getparfloat("dxi",&dxi)) err("must specify dxi!\n");
        if (!getparint("sgn",&sgn)) err("must specify sgn!\n");

	/* optional parameters	*/
	if (!getparint("nxso",&nxso)) nxso = (nxt-1)*2+1;
	if (!getparfloat("fxso",&fxso)) fxso = fxt;
	if (!getparfloat("dxso",&dxso)) dxso = dxt*0.5;
        if (!getparint("nxgo",&nxgo)) nxgo = (nxt-1)*2+1;
        if (!getparfloat("fxgo",&fxgo)) fxgo = fxt;
        if (!getparfloat("dxgo",&dxgo)) dxgo = dxt*0.5;
        if (!getparfloat("exgo",&exgo)) exgo = fxgo + (nxgo-1)*dxgo;
        if (!getparint("antiali",&antiali)) antiali=1;
	if (!getparint("nzo",&nzo)) nzo = (nzt-1)*5+1;
	if (!getparfloat("fzo",&fzo)) fzo = fzt;
	if (!getparfloat("dzo",&dzo)) dzo = dzt*0.2;

        fzo = fzo*sgn;
        dzo = dzo*sgn;

	exso = fxso+(nxso-1)*dxso;
	ezo = fzo+(nzo-1)*dzo;
	if(fxt>fxso||fxt>fxgo||ext<exso||ext<exgo||fzt>fzo||ezt<ezo) 
		err(" migration output range is out of traveltime table!\n");

	if (!getparfloat("v0",&v0)) v0 = 1500;
	if (!getparfloat("dvz",&dvz)) dvz = 0;
	if (!getparfloat("angmax",&angmax)) angmax = 60.;
	if  (angmax<0.00001) err("angmax must be positive!\n");
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;
	if (!getparfloat("offmax",&offmax)) offmax = 999999;
	if (!getparfloat("fmax",&fmax)) fmax = 0.25/dt;

	if (!getparint("ntr",&ntr))	ntr = 100000;
	if (!getparint("mtr",&mtr))	mtr = 100;
        decodeSurfaces(sgn,&nsrf,&nxzsrf,&xsrf,&zsrf);
        makesurf(0.025*dxso,nsrf,nxzsrf,xsrf,zsrf,&srf);

	fprintf(jpfp,"\n");
	fprintf(jpfp," Datuming parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," infile=%s \n",infile);
	fprintf(jpfp," outfile=%s \n",outfile);
	fprintf(jpfp," ttfile=%s \n",ttfile);
        fprintf(jpfp," \n");

        fprintf(jpfp," nzt=%d fzt=%g dzt=%g\n",nzt,fzt,dzt);
	fprintf(jpfp," nxt=%d fxt=%g dxt=%g\n",nxt,fxt,dxt);
 	fprintf(jpfp," ns=%d fs=%g ds=%g\n",ns,fs,ds);
	fprintf(jpfp," \n");
        fprintf(jpfp," nxi=%d fxi=%g dxi=%g sgn=%d\n",nxi,fxi,dxi,sgn);
        fprintf(jpfp," \n");
	fprintf(jpfp," nzo=%d fzo=%g dzo=%g\n",nzo,fzo,dzo);
	fprintf(jpfp," nxso=%d fxso=%g dxso=%g\n",nxso,fxso,dxso);
        fprintf(jpfp," nxgo=%d fxgo=%g dxgo=%g\n",nxgo,fxgo,dxgo);
	fprintf(jpfp," \n");

        /* compute Z coordinate for the surfaces  */
        szif = ealloc2float(2,nxi);
        sz = ealloc1float(ns);
        nangl = ealloc1float(nxi);
        zcoorSurfaces(fxi,dxi,nxi,fs,ds,ns,srf,szif,sz,nangl);
	
	/* compute reference traveltime and slowness  */
	rmax = MAX(es-fxt,ext-fs);
	rmax = MIN(rmax,0.5*offmax+aperx);
	nr = 2+(int)(rmax/dxgo);
        tb = ealloc3float(nzt,nr,ns);
        pb = ealloc3float(nzt,nr,ns);
        sigb = ealloc3float(nzt,nr,ns);
        angb = ealloc3float(nzt,nr,ns);
        for (is=0;is<ns;is++){
            v00 = v0 + sz[is]*dvz;
            timeb(nr,nzt,dxso,dzt,fzt,sz[is],dvz,v00,tb[is],pb[is],sigb[is],
                angb[is]);
        }

	fprintf(jpfp," nt=%d ft=%g dt=%g \n",nt,ft,dt);
 	fprintf(jpfp," fmax=%g\n",fmax);
	fprintf(jpfp," v0=%g dvz=%g \n",v0,dvz);
 	fprintf(jpfp," aperx=%g offmax=%g angmax=%g\n",aperx,offmax,angmax);
 	fprintf(jpfp," ntr=%d mtr=%d \n",ntr,mtr);
	fprintf(jpfp," ================\n");
	fflush(jpfp);

	/* allocate space */
	nxgt = (int)((exgo-fxgo)/dxgo+1);
	ng = ealloc1int(nxgt);
	ttab = ealloc3float(nzt,nxt,ns);
	tt = alloc2float(nzt,nxt);
	tsum = alloc2float(nzt,nxt);
        tbs = alloc2float(nzt,nr);
        pbs = alloc2float(nzt,nr);
        sigbs = alloc2float(nzt,nr);
        angbs = alloc2float(nzt,nr);

	if (nxso<nxgo) {
		datg = ealloc3float(nt,nxso,nxgt);
	 	memset((void *) datg[0][0],0,nxgt*nxso*nt*sizeof(float));
	} else {
		datg = ealloc3float(nt,nxgo,nxgt);
		memset((void *) datg[0][0],0,nxgt*nxgo*nt*sizeof(float));
	}

	memset((void *) ng,0,nxgt*sizeof(int));

 	fprintf(jpfp," input traveltime tables \n");
                       
	/* compute traveltime residual	*/
	for(is=0; is<ns; ++is){
		nseek = (off_t) nxt*nzt*is;
		efseeko(ttfp,nseek*((off_t) sizeof(float)),SEEK_SET);
		efread(ttab[is][0],sizeof(float),nxt*nzt,ttfp);
		s = fs+is*ds;
		resit(nxt,fxt,dxt,nzt,nr,dxgo,tb[is],ttab[is],s);
	}

	fprintf(jpfp," start datuming sources... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
	
        mtmax = 2*dxso*sin(angmax*PI/180.)/(v0*dt);
        if(mtmax<1) mtmax = 1;

	jtr = 1;
	ktr = 0;
	fxin = fxso;
	io = 0;
	s = 0.;

        for(is=0; is<nxgt; ++is)
                ng[is]=0;

	/* Store headers in tmpfile while getting a count */
	hdrfp = etmpfile(); 

	do {

	    sx = tr.sx;
	    gx = tr.gx;

	    if(gx!=s){
		fxin = sx;
		s = gx;
	    }

	    io = (int)((gx-fxgo)/dxgo);

	    if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && 
	        MAX(gx-sx,sx-gx)<=offmax ){

		/* Number of traces in receiver's gathers  */ 
                ng[io]++;

                efwrite(&tr, 1, HDRBYTES, hdrfp);

		/*  Down or up ward continuation of sources  */
	    	as = (sx-fs)/ds;
	    	is = (int)as;
		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tsum);
                sum2(nr,nzt,1-res,res,pb[is],pb[is+1],pbs);
                sum2(nr,nzt,1-res,res,tb[is],tb[is+1],tbs);