sumigtopo2d.c

su 的源代码库

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

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

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

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
" SUMIGTOPO2D - Kirchhoff Depth Migration of 2D postack/prestack data	",
"	     from the (variable topography) recording surface		",
" 									",
"    sumigtopo2d  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 input trace           ",
" dxi                   horizontal spacing of input data                ",
" nxi                   number of input trace locations in surface      ",
"									",
" Optional Parameters:							",
" dt= or from header (dt) 	time sampling interval of input data	",
" ft= or from header (ft) 	first time sample of input data		",
" dxm= or from header (d2) 	sampling interval of midpoints 		",
" surf=\"0,0;99999,0\"  Recording surface \"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		",	
" fxo=fxt               x-coordinate of first output trace 		",
" dxo=0.5*dxt		horizontal spacing of output trace 		",
" nxo=2*(nxt-1)+1  	number of output traces 			",	
" off0=0               	first offest in output 				",
" doff=99999		offset increment in output 			",
" noff=1       		number of offsets in output 			",	
" 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		",
" ls=1	                flag for line source				",
" 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 rayt2dtopo (or any 	",
"    other one that considers topography )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 seismic traces must be SU format and can be any type of 	",
"    gathers (common shot, common offset, common CDP, and so on).	", 
" 3. Migrated traces are output in CDP gathers if velocity analysis	",
"    is required, with dimension nxo*noff*nzo.  			", 
" 4. If the offset value of an input trace is not in the offset array 	",
"    of output, the nearest one in the array is chosen. 		",
" 5. Amplitudes are computed using the reference velocity profile, v(z),",
"    specified by the parameters v0= and dvz=.				",
" 6. Input traces must specify source and receiver positions via the header",
"    fields tr.sx and tr.gx. Offset is computed automatically.		",
"									",
NULL};
/*
 * Author:  Zhenyue Liu, 03/01/95,  Colorado School of Mines
 *
 *	    Trino Salinas, 07/01/96, Colorado School of Mines,
 *          extended the code to migrate data from a nonflat 
 *          recording surface.
 *
 * References :
 *
 * Bleistein, N., Cohen, J., and Hagin, F., 1987, Two and one-half
 *   dimensional Born inversion with arbitrary reference: Geophysics,
 *   52, 26-36.
 *
 * Liu,Z., 1993, A Kirchhoff approach to seismic modeling and 
 *   pre-stack depth migration: CWP Annual Report, CWP, Colorado
 *   School of Mines.
 *
 * Wiggins, J. W., 1984, Kirchhoff integral extrapolation and
 *   migration of nonplanar data: Geophysics, 49, 1239-
 *   1248.
 *
**************** end self doc ***********************************/
/* TYPEDEFS */
/* Geometry of the recording surface */
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;

/* Prototypes of functions used internally */
  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 mig2d(float *trace,int nt,float ft,float dt,
	float sx,float gx,float **mig,float aperx,
  	int nx,float fx,float dx,float nz,float fz,float dz,
	int ls,int mtmax,float dxm,float fmax,float angmax,
	float **tbs,float **tbg,float **pbs,
	float **pbg,float **angbs,float **angbg,int nr,float **tsum, 
	int nzt,float fzt,float dzt,int nxt,float fxt,float dxt,float *szif,
	float nangls,float nanglg,float **sigbs,float **sigbg);
  void decodeSurfaces(int *nrPtr,int **nxzPtr,float ***xPtr,float ***zPtr);
  int decodeSurface(char *string,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,nxo,ns,noff,nr,is,io,ixo,izo;
	int   ls,ntr,jtr,ktr,mtr,mtmax,nsrf,*nxzsrf,nxi;
	off_t nseek;
	float ft,fzt,fxt,fzo,fxo,fs,off0,dt,dzt,dxt,dzo,dxo,ds,doff,dxm,
	      ext,ezt,ezo,exo,es,s,scal,**xsrf,**zsrf,fxi,dxi;
	float v0,v00,dvz,fmax,angmax,offmax,rmax,aperx,sx,gx,nangls,
	      nanglg;
	float ***mig,***ttab,***tb,***pb,***angb,**tsum,**tt,*szif,
              *nangl,*sz,**tbs,**tbg,**pbs,**pbg,**angbs,**angbg,
              ***sigb,**sigbs,**sigbg;
        Surface *srf;
	
	char *infile="stdin",*outfile="stdout",*ttfile,*jpfile;
	FILE *infp,*outfp,*ttfp,*jpfp;


	/* 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.; 
 	if (!getparfloat("dxm",&dxm)) dxm = tr.d2;
	if  (dxm<0.0000001) err("dxm must be positive!\n");
	
	/* get traveltime tabel 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");

	/* optional parameters	*/
	if (!getparint("nxo",&nxo)) nxo = (nxt-1)*2+1;
	if (!getparfloat("fxo",&fxo)) fxo = fxt;
	if (!getparfloat("dxo",&dxo)) dxo = dxt*0.5;
	if (!getparint("nzo",&nzo)) nzo = (nzt-1)*5+1;
	if (!getparfloat("fzo",&fzo)) fzo = fzt;
	if (!getparfloat("dzo",&dzo)) dzo = dzt*0.2;
	exo = fxo+(nxo-1)*dxo;
	ezo = fzo+(nzo-1)*dzo;
	if(fxt>fxo || ext<exo || 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");
	mtmax = 2*dxm*sin(angmax*PI/180.)/(v0*dt);
	if(mtmax<1) mtmax = 1;
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;
	if (!getparfloat("offmax",&offmax)) offmax = 3000;
	if (!getparfloat("fmax",&fmax)) fmax = 0.25/dt;
	if (!getparint("noff",&noff))	noff = 1;
	if (!getparfloat("off0",&off0)) off0 = 0.;
	if (!getparfloat("doff",&doff)) doff = 99999;

	if (!getparint("ls",&ls))	ls = 1;
	if (!getparint("ntr",&ntr))	ntr = 100000;
	if (!getparint("mtr",&mtr))	mtr = 100;

        /* Estimation of the recording surface  */
        decodeSurfaces(&nsrf,&nxzsrf,&xsrf,&zsrf);
        makesurf(dxm*0.025,nsrf,nxzsrf,xsrf,zsrf,&srf);

	fprintf(jpfp,"\n");
	fprintf(jpfp," Migration 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\n",nxi,fxi,dxi);
        fprintf(jpfp," \n");
	fprintf(jpfp," nzo=%d fzo=%g dzo=%g\n",nzo,fzo,dzo);
	fprintf(jpfp," nxo=%d fxo=%g dxo=%g\n",nxo,fxo,dxo);
	fprintf(jpfp," \n");

	/* compute the z coordinate of the recording surface  */
        szif = ealloc1float(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/dxo);
        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,dxo,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," dxm=%g fmax=%g\n",dxm,fmax);
 	fprintf(jpfp," noff=%d off0=%g doff=%g\n",noff,off0,doff);
	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 ls=%d\n",ntr,mtr,ls);
	fprintf(jpfp," ================\n");
	fflush(jpfp);

	/* allocate space */
	mig = ealloc3float(nzo,nxo,noff);
	ttab = ealloc3float(nzt,nxt,ns);
	tt = ealloc2float(nzt,nxt);
	tsum = ealloc2float(nzt,nxt);
        tbs = alloc2float(nzt,nr);
        tbg = alloc2float(nzt,nr);
        pbs = alloc2float(nzt,nr);
        pbg = alloc2float(nzt,nr);
        sigbs = alloc2float(nzt,nr);
        sigbg = alloc2float(nzt,nr);
        angbs = alloc2float(nzt,nr);
        angbg = alloc2float(nzt,nr);


 	memset((void *) mig[0][0],0,noff*nxo*nzo*sizeof(float)); 

 	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);
		fread(ttab[is][0],sizeof(float),nxt*nzt,ttfp);
		s = fs+is*ds;
		resit(nxt,fxt,dxt,nzt,nr,dxo,tb[is],ttab[is],s);
		
	}

	fprintf(jpfp," start migration ... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
	

	jtr = 1;
	ktr = 0;
	do {
		/* determine offset index	*/
	    float as,res;
	    sx = tr.sx;
	    gx = tr.gx;
	    io = (int)((gx-sx-off0)/doff+0.5);
	    if(io<0) io = 0;
	    if(io>=noff) io = noff-1;

	    if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && 
	       MAX(gx-sx,sx-gx)<=offmax ){
		/*     migrate this trace	*/
	    	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);
                sum2(nr,nzt,1-res,res,sigb[is],sigb[is+1],sigbs);
                sum2(nr,nzt,1-res,res,angb[is],angb[is+1],angbs);

	    	as = (gx-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],tt);
		sum2(nxt,nzt,1,1,tt,tsum,tsum);
                sum2(nr,nzt,1-res,res,pb[is],pb[is+1],pbg);
                sum2(nr,nzt,1-res,res,tb[is],tb[is+1],tbg);
                sum2(nr,nzt,1-res,res,sigb[is],sigb[is+1],sigbg);
                sum2(nr,nzt,1-res,res,angb[is],angb[is+1],angbg);

                as = (sx-fxi)/dxi;
                is = (int)as;
                if(is==nxi-1) is=nxi-2;
                res = as-is;
                if(res<=0.01) res = 0.0;
                if(res>=0.99) res = 1.0;
                nangls = (1-res)*nangl[is] + res*nangl[is+1];

                as = (gx-fxi)/dxi;
                is = (int)as;
                if(is==nxi-1) is=nxi-2;
                res = as-is;
                if(res<=0.01) res = 0.0;
                if(res>=0.99) res = 1.0;
                nanglg = (1-res)*nangl[is] + res*nangl[is+1];

		mig2d(tr.data,nt,ft,dt,sx,gx,mig[io],aperx,
		  nxo,fxo,dxo,nzo,fzo,dzo,
		  ls,mtmax,dxm,fmax,angmax,tbs,tbg,pbs,pbg,angbs,angbg,