sudatumk2dr.c

su 的源代码库

        free2float(ssx);
        free2float(ssz);
        free1float(sss);
}

/* residual traveltime calculation based  on reference   time	*/
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,
		float **tb,float **t,float x0)
{
	int ix,iz,jr;
	float xi,ar,sr,sr0;

	for(ix=0; ix<nx; ++ix){
		xi = fx+ix*dx-x0;
		ar = abs(xi)/dr;
		jr = (int)ar;
		sr = ar-jr;
		sr0 = 1.0-sr;
		if(jr>nr-2) jr = nr-2;
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] -= sr0*tb[jr][iz]+sr*tb[jr+1][iz];
	}
} 

/* lateral interpolation	*/

/* sum of two tables	*/
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t)
{
	int ix,iz;

	for(ix=0; ix<nx; ++ix) 
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] = a1*t1[ix][iz]+a2*t2[ix][iz];
}
 
/* compute  reference traveltime and slowness	*/
      void timeb(int nr,int nz,float dr,float dz,float fz,float z0,float a,
	float v0,float **t,float **p,float **sig,float **ang)
{
	int  ir,iz;
	float r,z,v,rc,oa,temp,rou,zc;


	if( a==0.0) {
		for(ir=0,r=0;ir<nr;++ir,r+=dr)
			for(iz=0,z=fz-z0;iz<nz;++iz,z+=dz){
				rou = sqrt(r*r+z*z);
				if(rou<dz) rou = dz;
				t[ir][iz] = rou/v0;
				p[ir][iz] = r/(rou*v0);
				sig[ir][iz] = v0*rou;
				ang[ir][iz] = asin(r/rou);
			}
	} else {
		oa = 1.0/a; 	zc = v0*oa;
		for(ir=0,r=0;ir<nr;++ir,r+=dr)
			for(iz=0,z=fz+zc-z0;iz<nz;++iz,z+=dz){
				rou = sqrt(r*r+z*z);
				v = v0+a*(z-zc+z0);
				if(ir==0){ 
					t[ir][iz] = log(v/v0)*oa;
					p[ir][iz] = 0.0;
					ang[ir][iz] = 0.0;
					sig[ir][iz] = 0.5*(z+0.1*dz-zc+z0)
						*(v0+v);
				} else {
					rc = (r*r+z*z-zc*zc)/(2.0*r);
					rou = sqrt(zc*zc+rc*rc);
					t[ir][iz] = log((v*(rou+rc))
						/(v0*(rou+rc-r)))*oa;
					p[ir][iz] = sqrt(rou*rou-rc*rc)
						/(rou*v0);
                                        temp = v0*p[ir][iz];
                                        if(temp>1.0) temp = 1.0;
                                        ang[ir][iz] = asin(temp);
					sig[ir][iz] = a*rou*r;
				}
			}
	}
}

void filt(float *trace,int nt,float dt,float fmax,int m,float *trf,
	     int sgn);

  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,int verbose)
/*****************************************************************************
Datum one trace 
******************************************************************************
Input:
*trace		one seismic trace 
nt		number of time samples in seismic trace
ft		first time sample of seismic trace
dt		time sampleing interval in seismic trace
sx,gx		lateral coordinates of source and geophone (ignored)
aperx		lateral aperture in the datuming process
nx,fx,dx,nz,fz,dz	dimension parameters of datuming region
mtmax		number of time samples in triangle filter
xm		datuming point ( source=sx or receiver=gx )
fmax		frequency-highcut for input trace	 
fxi             x-coordinate of the first surface location
dxi             horizontal spacingon surface
nxi             number of input surface locations
angmax		migration angle aperature from vertical 	 
(tb,pb,angb)	reference traveltime, lateral slowness, 
sigb		emergent angle and sigma
nr		number of lateral samples in reference quantities
tsum		sum of residual traveltimes from shot and receiver
nxt,fxt,dxt,nzt,fzt,dzt		dimension parameters of traveltime table
antiali         Anti-aliase filter flag
sgn             Sign of the datuming process(up or down ward)
szif            array[] of the z component of the reference surfaces
nangl		Angles that the normal unit vector of the
                current surface form with the vertical at source or geophone
                location

Output:
dat		Redatumed section
*****************************************************************************/
{
        int nxf,nxe,nxtf,nxte,it,ix,iz,iz0,izt0,nzp,jr,jz,jt,mt,jx,mr;
        float x,dis,rxz,ar,sr,sr0,z0,rdz,ampd,res0,am,am0,fxf,
	      ang,ax,ax0,pmin,odt=1.0/dt,pd,az,sz,sz0,
              at,td,res,temp,sig,sigma;
	float **tmt,**ampt,**ampti,*tm,*amp,*ampi,*tzt,*trf,*zpt;

	tmt = alloc2float(nzt,nxt);
	ampt = alloc2float(nzt,nxt);
	ampti = alloc2float(nzt,nxt);
	tm = alloc1float(nzt);
	tzt = alloc1float(nzt);
	amp = alloc1float(nzt);
	ampi = alloc1float(nzt);
	zpt = alloc1float(nxt);
	trf = alloc1float(nt+2*mtmax);

	z0 = (fz-fzt)/dzt + 0.0*sx + 0.0*gx;
	rdz = dz/dzt;
	pmin = 1.0/(2.0*dx*fmax);
	fxf = fxi + (nxi-2)*dxi;
	
	filt(trace,nt,dt,fmax,mtmax,trf,sgn);

	rxz = (angmax==90)?0.0:1.0/tan(angmax*PI/180.);
	nxtf = (xm-aperx-fxt)/dxt;
	if(nxtf<0) nxtf = 0;
	nxte = (xm+aperx-fxt)/dxt+1;
	if(nxte>=nxt) nxte = nxt-1;

	/* compute amplitudes and filter length	*/
	for(ix=nxtf; ix<=nxte; ++ix){
		x = fxt+ix*dxt;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = ((dis-dxt)*rxz-fzt)/dzt-1;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;

		ar = (xm>=x)?(xm-x)/dx:(x-xm)/dx;
		jr = (int)ar;
		if(jr>nr-2) jr = nr-2;
		sr = ar-jr;
		sr0 = 1.0-sr;
                sig = ((xm-x)<0)?1.0:-1.0;
		zpt[ix] = fzt+(nzt-1)*dzt;

		for(iz=izt0; iz<nzt; ++iz){
			sigma = sr0*sigb[jr][iz]+sr*sigb[jr+1][iz];
			ang = sr0*angb[jr][iz]+sr*angb[jr+1][iz];
                        ang = (sig*ang - nangl);
			ampd = cos(ang)/(cos(nangl)*sqrt(sigma));

	/* Filter of 90 degrees to the operator (Peels,1988)

           Peels, G. L., 1988, True amplitude wave field extrapolation
           with applications in  seismic shot record redatuming, PhD
           thesis, Delft University of Technology.              */

			if(ABS(ang)>=(PI/2)) {
				ampd=0.0;
			}

			if(ampd<0.0) ampd = -ampd;
			ampt[ix][iz] = ampd;

			pd = sr0*pb[jr][iz]+sr*pb[jr+1][iz];
			if(pd<0.0) pd = -pd;
			temp = pd*dx*odt;
			if(temp<1) temp = 1.0;
			if(temp>mtmax) temp = mtmax;
			ampti[ix][iz] = ampd/(temp*temp);
			tmt[ix][iz] = temp;
			if(pd<pmin && zpt[ix]>fzt+(nzt-1.1)*dzt) 
				zpt[ix] = fzt+iz*dzt;

		}
	}

	nxf = (xm-aperx-fx)/dx+0.5;
	if(nxf<0) nxf = 0;

	if((xm+aperx)>=fxf)
		nxe = (fxf-fx)/dx + 0.5;
	else
		nxe = (xm+aperx-fx)/dx + 0.5;

	if(nxe>=nx) nxe = nx-1;

	am = (fx-fxi)/dxi;

	mr = (int)am;
	am = am - mr;
	if(am<=0.01) am = 0.;
	if(am>=0.99) am = 1.0;
	am0 = 1.0-am;
	if(mr<0) mr = 0;
	if(mr+nxe>=nxi-1) 
		err("Topography definition is out of range!\n");
	
	/* interpolate amplitudes and filter length along lateral	*/
	for(ix=nxf; ix<=nxe; ++ix){
		x = fx+ix*dx;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = (dis*rxz-fzt)/dzt;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;
		iz0 = (dis*rxz-fz)/dz;
		if(iz0<0) iz0 = 0;
		if(iz0>=nz) iz0 = nz-1;

		ax = (x-fxt)/dxt;
		jx = (int)ax;
		ax = ax-jx;
		if(ax<=0.01) ax = 0.;
		if(ax>=0.99) ax = 1.0;
		ax0 = 1.0-ax;
		if(jx>nxte-1) jx = nxte-1;
		if(jx<nxtf) jx = nxtf;

		ar = (xm>=x)?(xm-x)/dx:(x-xm)/dx;
		jr = (int)ar;
		if(jr>nr-2) jr = nr-2;
		sr = ar-jr;
		sr0 = 1.0-sr;

		for(iz=izt0; iz<nzt; ++iz){
		    tzt[iz] = ax0*tsum[jx][iz]+ax*tsum[jx+1][iz]
				+sr0*tb[jr][iz]+sr*tb[jr+1][iz];

		    amp[iz] = ax0*ampt[jx][iz]+ax*ampt[jx+1][iz];
		    ampi[iz] = ax0*ampti[jx][iz]+ax*ampti[jx+1][iz];
		    tm[iz] = ax0*tmt[jx][iz]+ax*tmt[jx+1][iz];
		}

		nzp = (ax0*zpt[jx]+ax*zpt[jx+1]-fz)/dz+0.5;
		if(nzp<iz0) nzp = iz0;
                if(nzp>nz) nzp = nz;

                iz = (ABS(am0*szif[ix+mr][1]+am*szif[ix+mr+1][1])-fz)/dz;
                if(iz>=nz)
                        err("Datuming surface is out of output range!\n");
                az = z0+iz*rdz;
                jz = (int)az;
                if(jz>=nzt-1) jz = nzt-2;
                sz = az-jz;
                sz0 = 1.0-sz;
                td = sz0*tzt[jz]+sz*tzt[jz+1];
                at = (sgn*td-ft)*odt;
                if ((iz<nzp) && (antiali)) {
                /* interpolate along depth if operater aliasing   */
                        at = at + mtmax;
                        jt = (int)at;
                        ampd = sz0*ampi[jz]+sz*ampi[jz+1];
                        mt = (int)(0.5+sz0*tm[jz]+sz*tm[jz+1]);
                        res = ABS(at-jt);
                        res0 = 1.0-res;
                        for (it=0; it<nt; it++){
                          if(it+jt >= mtmax && jt < nt-it+mtmax-sgn){
                             temp = (res0*(-trf[it+jt-mt]+2.0*trf[it+jt]
                                 -trf[it+jt+mt])+res*(-trf[it+jt-mt+sgn]
                                 +2.0*trf[it+jt+sgn]-trf[it+jt+mt+sgn]))*ampd;
                             dat[ix][it] += temp;
                          }
                        }

                }
                /* interpolate along depth if not operater aliasing     */
                else{
                        jt = (int)at;
                        ampd = sz0*amp[jz]+sz*amp[jz+1];
                        res = ABS(at-jt);
                        res0 = 1.0-res;
                        for (it=0; it<nt; it++){
                          if(it+jt >= 0 && jt < nt-it-sgn){
                            temp=(res0*trace[it+jt]+res*trace[it+jt+sgn])*ampd;
                            dat[ix][it] += temp;
                          }
                        }

                }

	}

	free2float(ampt);
	free2float(ampti);
	free2float(tmt);
	free1float(amp);
	free1float(ampi);
	free1float(zpt);
	free1float(tm);
	free1float(tzt);
	free1float(trf);
}

void filt(float *trace,int nt,float dt,float fmax,int m,
	  float *trf,int sgn)
/*******************************************************************
  Low-pass filter, integration and phase shift for input data	 
  input: 
    trace(nt)	single seismic trace
    fmax	high cut frequency
    sgn		Sign of the Datuming process
  output:
    trace(nt) 	filtered and phase-shifted seismic trace 
    tracei(nt) 	filtered, integrated and phase-shifted seismic trace 
********************************************************************/
{
	static int nfft=0, itaper, nw, nwf;
	static float *taper, *amp, *ampi, dw;
	int  it,iw,itemp;
	float temp, ftaper, const2, *rt;
	complex *ct;

	fmax *= 2.0*PI;
	ftaper = 0.1*fmax;
	const2 = 0.5*sqrt(2.0);

	if(nfft==0) {
        	/* Set up FFT parameters */
        	nfft = npfaro(2*(nt+2*m), 4*(nt+2*m));
        	if (nfft >= SU_NFLTS || nfft >= 720720)
                	err("Padded nt=%d -- too big", nfft);

        	nw = nfft/2 + 1;
		dw = 2.0*PI/(nfft*dt);

		itaper = 0.5+ftaper/dw;
		taper = ealloc1float(2*itaper+1);
		for(iw=-itaper; iw<=itaper; ++iw){
			temp = (float)iw/(1.0+itaper); 
			taper[iw+itaper] = (1-temp)*(1-temp)*(temp+2)/4;
		}

		nwf = 0.5+fmax/dw;
		if(nwf>nw-itaper-1) nwf = nw-itaper-1;
		amp = ealloc1float(nwf+itaper+1);
		ampi = ealloc1float(nwf+itaper+1);
		amp[0] = ampi[0] = 0.;
		for(iw=1; iw<=nwf+itaper; ++iw){
			amp[iw] = sqrt(dw*iw)/nfft;
			ampi[iw] = 0.5/(1-cos(iw*dw*dt));
		}
	}

        /* Allocate fft arrays */
        rt   = ealloc1float(nfft);
        ct   = ealloc1complex(nw);

        memcpy(rt, trace, nt*FSIZE);
        memset((void *) (rt + nt), 0, (nfft-nt)*FSIZE); 
        pfarc(1, nfft, rt, ct);

	for(iw=nwf-itaper;iw<=nwf+itaper;++iw){
		itemp = iw-(nwf-itaper);
		ct[iw].r = taper[itemp]*ct[iw].r; 
		ct[iw].i = taper[itemp]*ct[iw].i; 
	}
	for(iw=nwf+itaper+1;iw<nw;++iw){
		ct[iw].r = 0.; 
		ct[iw].i = 0.; 
	}

	for(iw=0; iw<=nwf+itaper; ++iw){
		/* phase shifts PI/4 - Half derivative	*/
		temp = (ct[iw].r-sgn*ct[iw].i)*amp[iw]*const2;
		ct[iw].i = (sgn*ct[iw].r + ct[iw].i)*amp[iw]*const2;
		ct[iw].r = temp;
	}
              
        pfacr(-1, nfft, ct, rt);
		
        /* Load traces back in */
	for (it=0; it<nt; ++it) trace[it] = rt[it];

        /* Integrate traces   */
	for(iw=0; iw<=nwf+itaper; ++iw){
		ct[iw].i = ct[iw].i*ampi[iw];
		ct[iw].r = ct[iw].r*ampi[iw];
	}
        pfacr(-1, nfft, ct, rt);
        for (it=0; it<m; ++it)  trf[it] = rt[nfft-m+it];
        for (it=0; it<nt+m; ++it)  trf[it+m] = rt[it];

	free1float(rt);
	free1complex(ct);
}