SU的PSTM只能做单个offset偏移，要求输入数据每道数据等间隔CDP，且要求速度CDP等间隔，稍作修改可以实现任意数据集的偏移。但实际处理的时候最好还是按offset输出结果用sustack叠加。注意输入速度文件需要连续CDP从1开始。

/* Copyright (c) Colorado School of Mines, 2011.*/

/* All rights reserved.       */

/* SUKTMIG2D: $Revision: 1.7 $ ; $Date: 2013/08/14 18:34:35 $*/

#include "su.h"

#include "segy.h"

#include "header.h"

#include "stdio.h"

#include "stdlib.h"

/*********************** self documentation **********************/

char *sdoc[] = {

" ",

" 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.cdpfirst cdp in data",

" firstcdp=fcdpdatafirst cdp number in velocity file",

" lastcdp=from headerlast cdp number in velocity file",

" dcdp=from headernumber of cdps between consecutive traces",

" angmax=40maximum aperture angle for migration (degrees)",

" hoffset=.5*tr.offsethalf offset (m)",

" nfc=16number 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)",

NULL };

#define LOOKFAC 2       /* Look ahead factor for npfaro   */

#define PFA_MAX 720720  /* Largest allowed nfft   */

/* Prototype of functions used internally */

void 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 cdp_all = 4200;

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 */

float **migcrp; = NULL;

float **crp = NULL;

float ofs1 = 0;

float ofs2 = 0;

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;

firstcdp = 1000000;

lastcdp = 1;

int *cdpnum;

float *offsetnum;

cdpnum=alloc1int(64028);

offsetnum = alloc1float(64028);

do {

++ntr;

if (intrace.cdp >= lastcdp) lastcdp = intrace.cdp;

if (intrace.cdp <= firstcdp) firstcdp = intrace.cdp;

cdpnum[ntr-1] = intrace.cdp;

offsetnum[ntr-1] = intrace.offset;

efwrite(&intrace, HDRBYTES, 1, hfp);

efwrite(intrace.data, FSIZE, nt, tracefp);

} while (gettr(&intrace));

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);

/* Set up FFT parameters */

nfft = npfaro(nt, LOOKFAC*nt);

if (nfft >= SU_NFLTS || nfft >= PFA_MAX)

err("Padded nt=%d -- 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= %d", nc);

/* Allocate space */

data = alloc2float(nt, ntr);

lowpass = alloc2float(nt, nc + 1);

mig = alloc2float(nt, cdp_all);

migcrp = alloc2float(nt, cdp_all);

crp = alloc2float(nt,42*72);

vel = alloc2float(nt, cdp_all);

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*cdp_all, vfp);

efclose(vfp);

/* Zero all arrays */

memset((void *)mig[0], 0, nt*cdp_all*FSIZE);

memset((void *)migcrp[0], 0, nt*cdp_all*FSIZE);

memset((void *)crp[0], 0, nt*42*72*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");

int itrace;

for (itrace = 0; itrace<ntr; ++itrace){

float perc;

imp = cdpnum[itrace];

mp = imp*dx;

perc = itrace*100.0 / (ntr - 1);

if (fmod(itrace * 100, ntr - 1) == 0 && verbose)

warn("migrated %gn ", perc);

warn("itrace= %d  cdp= %d  offset=  % ", itrace ,cdpnum[itrace],offsetnum[itrace]);

/* 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, data[itrace], 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][it];

xmax = tan((angmax + 10.0)*PI / 180.0)*v*t0;

lx = MAX(0, imp - ceil(xmax / dx));

ux = MIN(cdp_all, 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 = 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;

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, &amplo);

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, &amplo);

ints8r(8, dt, firstt, datahi, 0.0, 0.0, 1, &t, &amphi);

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 = 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, &amplo);

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, &amplo);

ints8r(8, dt, firstt, datahi, 0.0, 0.0, 1, &t, &amphi);

mig[iip][it] += geoms*obliq*angtaper*(wlo*amplo + whi*amphi);

}

}

}

}

/* Output migrated data */

erewind(hfp);

for (ix = 0; ix<cdp_all; ++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:

nfcnumber of Fourier coefficients to approximate ideal filter

nfftnumber of points in the fft

dttime sampling interval

fhicut-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);

}

}

}

二维偏移结果：

CRP道集：