SYNOPSIS
PERL PROGRAM NAME:
AUTHOR: Juan Lorenzo (Perl module only)
DATE:
DESCRIPTION:
Version:
USE
NOTES
Examples
SYNOPSIS
SEISMIC UNIX NOTES SUKTMIG2D - prestack time migration of a common-offset
section with the double-square root (DSR) operator
suktmig2d < infile vfile= [parameters] > outfile
Required Parameters:
vfile= rms velocity file (units/s) v(t,x) as a function
of time
dx= distance (units) between consecutive traces
Optional parameters:
fcdpdata=tr.cdp first cdp in data
firstcdp=fcdpdata first cdp number in velocity file
lastcdp=from header last cdp number in velocity file
dcdp=from header number of cdps between consecutive traces
angmax=40 maximum aperture angle for migration (degrees)
hoffset=.5*tr.offset half offset (m)
nfc=16 number of Fourier-coefficients to approximate
low-pass
filters. The larger nfc the narrower the filter
fwidth=5 high-end frequency increment for the low-pass
filters
in Hz. The lower this number the more the number
of lowpass filters to be calculated for each
input trace.
Caveat: this code may need some work
Notes:
Data must be preprocessed with sufrac to correct for the
wave-shaping factor using phasefac=.25 for 2D migration.
Input traces must be sorted into offset and cdp number.
The velocity file consists of rms velocities for all CMPs as a
function of vertical time and horizontal position v(t,x)
in C-style binary floating point numbers. It's easiest to
supply v(t,x) that has the same dimensions as the input data to
be migrated. Note that time t is the fast dimension in these
the input velocity file.
The units may be feet or meters, as long as these are
consistent.
Antialias filter is performed using (Gray,1992, Geoph. Prosp),
using nc low- pass filtered copies of the data. The cutoff
frequencies are calculated as fractions of the Nyquist
frequency.
The maximum allowed angle is 80 degrees(a 10 degree taper is
applied to the end of the aperture) define LOOKFAC 2 /* Look ahead factor for npfaro
define PFA_MAX 720720 /* Largest allowed nfft
Prototype of functions used internallyvoid lpfilt(int nfc, int nfft, float dt, float fhi, float *filter);
segy intrace; /* input traces
segy outtrace; /* migrated output traces
int
main(int argc, char **argv)
{
int i,k,imp,iip,it,ix,ifc; /* counters
int ntr,nt; /* x,t
int verbose; /* is verbose? */
int nc; /* number of low-pass filtered versions */
/* of the data for antialiasing */
int nfft,nf; /* number of frequencies */
int nfc; /* number of Fourier coefficients for low-pass filter
int fwidth; /* high-end frequency increment for the low-pass
/* filters */
int firstcdp=0; /* first cdp in velocity file */
int lastcdp=0; /* last cdp in velocity file */
int oldcdp=0; /* temporary storage */
int fcdpdata=0; /* first cdp in the data */
int olddeltacdp=0;
int deltacdp;
int ncdp=0; /* number of cdps in the velocity file */
int dcdp=0; /* number of cdps between consecutive traces
float dx=0.0; /* cdp sample interval
float hoffset=0.0; /* half receiver-source
float p=0.0; /* horizontal slowness of the migration operator
float pmin=0.0; /* maximum horizontal slowness for which there's
/* no aliasing of the operator
float dt; /* t sample interval
float h; /* offset
float x; /* aperture distance
float xmax=0.0; /* maximum aperture distance
float obliq; /* obliquity factor
float geoms; /* geometrical spreading factor
float angmax; /* maximum aperture angle
float mp,ip; /* mid-point and image-point coordinates
float t; /* time
float t0; /* vertical traveltime
float tmax; /* maximum time
float fnyq; /* Nyquist frequency
float ang; /* aperture angle
float angtaper=0.0; /* aperture-angle taper
float v; /* velocity
float *fc=NULL; /* cut-frequencies for low-pass filters
float *filter=NULL; /* array of low-pass filter values
float **vel=NULL; /* array of velocity values from vfile
float **data=NULL; /* input data array*/
float **lowpass=NULL; /* low-pass filtered version of the trace
float **mig=NULL; /* output migrated data array
register float *rtin=NULL,*rtout=NULL;/* real traces
register complex *ct=NULL; /* complex trace
/* file names
char *vfile=""; /* name of velocity file
FILE *vfp=NULL;
FILE *tracefp=NULL; /* temp file to hold traces*/
FILE *hfp=NULL; /* temp file to hold trace headers
float datalo[8], datahi[8];
int itb, ite;
float firstt, amplo, amphi;
cwp_Bool check_cdp=cwp_false; /* check cdp in velocity file */
/* Hook up getpar to handle the parameters
initargs(argc,argv);
requestdoc(0);
/* Get info from first trace
if (!gettr(&intrace)) err("can't get first trace");
nt=intrace.ns;
dt=(float)intrace.dt/1000000;
tmax=(nt-1)*dt;
MUSTGETPARFLOAT("dx",&dx);
MUSTGETPARSTRING("vfile",&vfile);
if (!getparfloat("angmax",&angmax)) angmax=40;
if (!getparint("firstcdp",&firstcdp)) firstcdp=intrace.cdp;
if (!getparint("fcdpdata",&fcdpdata)) fcdpdata=intrace.cdp;
if (!getparfloat("hoffset",&hoffset)) hoffset=.5*intrace.offset;
if (!getparint("nfc",&nfc)) nfc=16;
if (!getparint("fwidth",&fwidth)) fwidth=5;
if (!getparint("verbose",&verbose)) verbose=0;
h=hoffset;
/* Store traces in tmpfile while getting a count of number of traces
tracefp = etmpfile();
hfp = etmpfile();
ntr = 0;
do {
++ntr;
/* get new deltacdp value
deltacdp=intrace.cdp-oldcdp;
/* read headers and data
efwrite(&intrace,HDRBYTES, 1, hfp);
efwrite(intrace.data, FSIZE, nt, tracefp);
/* error trappings.
/* ...did cdp value interval change?
if ((ntr>3) && (olddeltacdp!=deltacdp)) {
if (verbose) {
warn("cdp interval changed in data");
warn("ntr=0 olddeltacdp=0 deltacdp=0"
,ntr,olddeltacdp,deltacdp);
check_cdp=cwp_true;
}
}
/* save cdp and deltacdp values
oldcdp=intrace.cdp;
olddeltacdp=deltacdp;
} while (gettr(&intrace));
/* get last cdp and dcdp
if (!getparint("lastcdp",&lastcdp)) lastcdp=intrace.cdp;
if (!getparint("dcdp",&dcdp)) dcdp=deltacdp - 1;
checkpars();
/* error trappings
if ( (firstcdp==lastcdp)
|| (dcdp==0)
|| (check_cdp==cwp_true) ) warn("Check cdp values in data!");
/* rewind trace file pointer and header file pointer
erewind(tracefp);
erewind(hfp);
/* total number of cdp's in data
ncdp=lastcdp-firstcdp+1;
/* Set up FFT parameters
nfft = npfaro(nt, LOOKFAC*nt);
if(nfft>= SU_NFLTS || nfft >= PFA_MAX)
err("Padded nt=0 -- too big",nfft);
nf = nfft/2 + 1;
/* Determine number of filters for antialiasing
fnyq= 1.0/(2*dt);
nc=ceil(fnyq/fwidth);
if (verbose)
warn(" The number of filters for antialiasing is nc= 0",nc);
/* Allocate space
data = alloc2float(nt,ntr);
lowpass=alloc2float(nt,nc+1);
mig= alloc2float(nt,ntr);
vel= alloc2float(nt,ncdp);
fc = alloc1float(nc+1);
rtin= ealloc1float(nfft);
rtout= ealloc1float(nfft);
ct= ealloc1complex(nf);
filter= alloc1float(nf);
/* Read data from temporal array
for (ix=0; ix<ntr; ++ix){
efread(data[ix],FSIZE,nt,tracefp);
}
/* read velocities
vfp=efopen(vfile,"r");
efread(vel[0],FSIZE,nt*ncdp,vfp);
efclose(vfp);
/* Zero all arrays
memset((void *) mig[0], 0,nt*ntr*FSIZE);
memset((void *) rtin, 0, nfft*FSIZE);
memset((void *) filter, 0, nf*FSIZE);
memset((void *) lowpass[0], 0,nt*(nc+1)*FSIZE);
/* Calculate cut frequencies for low-pass filters
for(i=1; i<nc+1; ++i){
fc[i]= fnyq*i/nc;
}
/* Start the migration process
/* Loop over input mid-points first
if (verbose) warn("Starting migration process...\n");
for(imp=0; imp<ntr; ++imp){
float perc;
mp=imp*dx;
perc=imp*100.0/(ntr-1);
if(fmod(imp*100,ntr-1)==0 && verbose)
warn("migrated 0\n ",perc);
/* Calculate low-pass filtered versions
/* of the data to be used for antialiasing
for(it=0; it<nt; ++it){
rtin[it]=data[imp][it];
}
for(ifc=1; ifc<nc+1; ++ifc){
memset((void *) rtout, 0, nfft*FSIZE);
memset((void *) ct, 0, nf*FSIZE);
lpfilt(nfc,nfft,dt,fc[ifc],filter);
pfarc(1,nfft,rtin,ct);
for(it=0; it<nf; ++it){
ct[it] = crmul(ct[it],filter[it]);
}
pfacr(-1,nfft,ct,rtout);
for(it=0; it<nt; ++it){
lowpass[ifc][it]= rtout[it];
}
}
/* Loop over vertical traveltimes
for(it=0; it<nt; ++it){
int lx,ux;
t0=it*dt;
v=vel[imp*dcdp+fcdpdata-1][it];
xmax=tan((angmax+10.0)*PI/180.0)*v*t0;
lx=MAX(0,imp - ceil(xmax/dx));
ux=MIN(ntr,imp + ceil(xmax/dx));
/* loop over output image-points to the left of the midpoint
for(iip=imp; iip>lx; --iip){
float ts,tr;
int fplo=0, fphi=0;
float ref,wlo,whi;
ip=iip*dx;
x=ip-mp;
ts=sqrt( pow(t0/2,2) + pow((x+h)/v,2) );
tr=sqrt( pow(t0/2,2) + pow((h-x)/v,2) );
t= ts + tr;
if(t>=tmax) break;
geoms=sqrt(1/(t*v));
obliq=sqrt(.5*(1 + (t0*t0/(4*ts*tr))
- (1/(ts*tr))*sqrt(ts*ts - t0*t0/4)*sqrt(tr*tr - t0*t0/4)));
ang=180.0*fabs(acos(t0/t))/PI;
if(ang<=angmax) angtaper=1.0;
if(ang>angmax) angtaper=cos((ang-angmax)*PI/20);
/* Evaluate migration operator slowness p to determine
/* the low-pass filtered trace for antialiasing
pmin=1/(2*dx*fnyq);
p=fabs((x+h)/(pow(v,2)*ts) + (x-h)/(pow(v,2)*tr));
if(p>0){fplo=floor(nc*pmin/p);}
if(p==0){fplo=nc;}
ref=fmod(nc*pmin,p);
wlo=1-ref;
fphi=++fplo;
whi=ref;
itb=MAX(ceil(t/dt)-3,0);
ite=MIN(itb+8,nt);
firstt=(itb-1)*dt;
/* Move energy from CMP to CIP
if(fplo>=nc){
for(k=itb; k<ite; ++k){
datalo[k-itb]=lowpass[nc][k];
}
ints8r(8,dt,firstt,datalo,0.0,0.0,1,&t,&lo);
mig[iip][it] +=geoms*obliq*angtaper*amplo;
} else if(fplo<nc){
for(k=itb; k<ite; ++k){
datalo[k-itb]=lowpass[fplo][k];
datahi[k-itb]=lowpass[fphi][k];
}
ints8r(8,dt,firstt,datalo,0.0,0.0,1,&t,&lo);
ints8r(8,dt,firstt,datahi,0.0,0.0,1,&t,&hi);
mig[iip][it] += geoms*obliq*angtaper*(wlo*amplo + whi*amphi);
}
}
/* loop over output image-points to the right of the midpoint
for(iip=imp+1; iip<ux; ++iip){
float ts,tr;
int fplo=0, fphi;
float ref,wlo,whi;
ip=iip*dx;
x=ip-mp;
t0=it*dt;
ts=sqrt( pow(t0/2,2) + pow((x+h)/v,2) );
tr=sqrt( pow(t0/2,2) + pow((h-x)/v,2) );
t= ts + tr;
if(t>=tmax) break;
geoms=sqrt(1/(t*v));
obliq=sqrt(.5*(1 + (t0*t0/(4*ts*tr))
- (1/(ts*tr))*sqrt(ts*ts
- t0*t0/4)*sqrt(tr*tr
- t0*t0/4)));
ang=180.0*fabs(acos(t0/t))/PI;
if(ang<=angmax) angtaper=1.0;
if(ang>angmax) angtaper=cos((ang-angmax)*PI/20.0);
/* Evaluate migration operator slowness p to determine the
/* low-pass filtered trace for antialiasing
pmin=1/(2*dx*fnyq);
p=fabs((x+h)/(pow(v,2)*ts) + (x-h)/(pow(v,2)*tr));
if(p>0){
fplo=floor(nc*pmin/p);
}
if(p==0){
fplo=nc;
}
ref=fmod(nc*pmin,p);
wlo=1-ref;
fphi=fplo+1;
whi=ref;
itb=MAX(ceil(t/dt)-3,0);
ite=MIN(itb+8,nt);
firstt=(itb-1)*dt;
/* Move energy from CMP to CIP
if(fplo>=nc){
for(k=itb; k<ite; ++k){
datalo[k-itb]=lowpass[nc][k];
}
ints8r(8,dt,firstt,datalo,0.0,0.0,1,&t,&lo);
mig[iip][it] +=geoms*obliq*angtaper*amplo;
} else if(fplo<nc){
for(k=itb; k<ite; ++k){
datalo[k-itb]=lowpass[fplo][k];
datahi[k-itb]=lowpass[fphi][k];
}
ints8r(8,dt,firstt,datalo,0.0,0.0,1,&t,&lo);
ints8r(8,dt,firstt,datahi,0.0,0.0,1,&t,&hi);
mig[iip][it] += geoms*obliq*angtaper*(wlo*amplo + whi*amphi);
}
}
}
}
/* Output migrated data
erewind(hfp);
for (ix=0; ix<ntr; ++ix) {
efread(&outtrace, HDRBYTES, 1, hfp);
for (it=0; it<nt; ++it) {
outtrace.data[it] = mig[ix][it];
}
puttr(&outtrace);
}
efclose(hfp);
return(CWP_Exit());}
void
lpfilt(int nfc, int nfft, float dt, float fhi, float *filter)
lpfilt -- low-pass filter using Lanczos Smoothing
(R.W. Hamming:"Digital Filtering",1977)Input:
nfc number of Fourier coefficients to approximate ideal filter
nfft number of points in the fft
dt time sampling interval
fhi cut-frequency
Output:
filter array[nf] of filter values
Notes: Filter is to be applied in the frequency domain
Author: CWP: Carlos Pacheco 2006
{
int i,j; /* counters
int nf; /* Number of frequencies (including Nyquist)
float onfft; /* reciprocal of nfft
float fn; /* Nyquist frequency
float df; /* frequency interval
float dw; /* frequency interval in radians
float whi;/* cut-frequency in radians
float w; /* radian frequency
nf= nfft/2 + 1;
onfft=1.0/nfft;
fn=1.0/(2*dt);
df=onfft/dt;
whi=fhi*PI/fn;
dw=df*PI/fn;
for(i=0; i<nf; ++i){
filter[i]= whi/PI;
w=i*dw;
for(j=1; j<nfc; ++j){
float c= sin(whi*j)*sin(PI*j/nfc)*2*nfc/(PI*PI*j*j);
filter[i] +=c*cos(j*w);
}
}}
User's notes (Juan Lorenzo) untested
CHANGES and their DATES
Import packages
instantiation of packages
Encapsulated hash of private variables
sub Step
collects switches and assembles bash instructions by adding the program name
sub note
collects switches and assembles bash instructions by adding the program name
sub clear
sub ang
sub angmax
sub angtaper
sub c
sub check_cdp
sub ct
sub data
sub dcdp
sub deltacdp
sub df
sub dt
sub dw
sub dx
sub fc
sub fcdpdata
sub filter
sub firstcdp
sub firstt
sub fn
sub fnyq
sub fphi
sub fplo
sub fwidth
sub geoms
sub h
sub hfp
sub hoffset
sub i
sub ifc
sub iip
sub imp
sub ip
sub it
sub itb
sub ite
sub ix
sub j
sub k
sub lastcdp
sub lowpass
sub lx
sub mig
sub mp
sub nc
sub ncdp
sub nf
sub nfc
sub nfft
sub nt
sub ntr
sub obliq
sub oldcdp
sub olddeltacdp
sub onfft
sub p
sub perc
sub phasefac
sub pmin
sub ref
sub rtin
sub rtout
sub t
sub t0
sub tmax
sub tr
sub tracefp
sub ts
sub ux
sub v
sub vel
sub verbose
sub vfile
sub vfp
sub w
sub whi
sub wlo
sub x
sub xmax
sub get_max_index
max index = number of input variables -1