From 289151c2b091a4be35831144e4039f8c4968f6f3 Mon Sep 17 00:00:00 2001 From: JanThorbecke Date: Thu, 8 Aug 2019 15:14:04 +0200 Subject: [PATCH] taper in time window in primaries app --- marchenko/demo/oneD/iterations.scr | 24 +++++ marchenko/demo/oneD/primaries.scr | 8 +- marchenko/demo/oneD/primariesFrame.scr | 24 ++++- marchenko/marchenko.c | 4 +- marchenko/marchenko_primaries.c | 138 ++++++++++++++++++------- marchenko/writeDataIter.c | 23 +++-- 6 files changed, 165 insertions(+), 56 deletions(-) create mode 100755 marchenko/demo/oneD/iterations.scr diff --git a/marchenko/demo/oneD/iterations.scr b/marchenko/demo/oneD/iterations.scr new file mode 100755 index 0000000..07839ed --- /dev/null +++ b/marchenko/demo/oneD/iterations.scr @@ -0,0 +1,24 @@ +#!/bin/bash + +#second reflector at time: +# 800/1800+600/2300 +# .70531400966183574878 => sample 176 +#thrids reflector at time: +# 800/1800+600/2300+800/2000 +# 1.10531400966183574878 sample 276 + +select=451 + +#for istart in 176 200 276 400 +for istart in 176 +do +(( iend = istart + 1 )) + +../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \ + nshots=601 verbose=1 istart=$istart iend=$iend fmax=90 \ + pad=1024 niter=45 smooth=10 niterskip=600 shift=20 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su + +suximage < f1min_${istart}043.su clip=1 title="${istart}" & + +done + diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr index 69cbf38..f5a82f2 100755 --- a/marchenko/demo/oneD/primaries.scr +++ b/marchenko/demo/oneD/primaries.scr @@ -13,14 +13,14 @@ export OMP_NUM_THREADS=20 #shot record to remove internal multiples select=451 -makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1 +makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1 ../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \ - nshots=601 verbose=1 istart=40 fmax=90 \ - niter=15 niterskip=600 shift=20 file_rr=pred_rr.su T=0 + nshots=601 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \ + niter=25 smooth=10 niterskip=600 shift=20 file_rr=pred_rr.su T=0 #for reference original shot record from Reflection matrix -suwind key=fldr min=$select max=$select < shotsdx5_rp.su > shotsx0.su +suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su # for displaying results diff --git a/marchenko/demo/oneD/primariesFrame.scr b/marchenko/demo/oneD/primariesFrame.scr index c1a7bf4..5017c4b 100755 --- a/marchenko/demo/oneD/primariesFrame.scr +++ b/marchenko/demo/oneD/primariesFrame.scr @@ -23,13 +23,13 @@ fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su basop file_in=shotw.su choice=5 | sugain scale=-1 > DD.su nts=1024 -itime=300 # select time sample to work on +itime=300 # select time sample to work on 1.2 seconds shift=20 #mute time nts-ii+shift to compute G_d (( itmax = nts-itime+shift )) suzero itmax=$itmax < DD.su > G_d.su #f1min = -DD(-t) = shotw.su -cp shotw.su f1min.su +cp shotw.su f1min0.su #first iteration cp G_d.su N0.su @@ -39,6 +39,26 @@ fconv file_in1=shotsdx5_rp.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax= suwind key=fldr min=451 max=451 < fconvN1.su | suximage +# Ni(t) = -1*RN0(-t) +sustack < fconvN1.su key=fldr > RN0.su +basop file_in=RN0.su choice=5 | sugain scale=-1 > N1.su + +#apply mute window for samples above nts-ii +(( itmax = nts-itime+shift )) +suzero itmax=$itmax < N1.su > N1m.su + + +exit +#display +file=fconvN1.su +file_base=${file%.su} +sumax < $file mode=abs outpar=nep +clipref=`cat nep | awk '{print $1/2}'` +supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \ + label1="time (s)" label2="distance (m)" \ + n1tic=2 d2=5 x1beg=0 x1end=1.5 d1num=0.4 \ + f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps + exit; diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c index 3182ff4..69a8977 100644 --- a/marchenko/marchenko.c +++ b/marchenko/marchenko.c @@ -35,7 +35,7 @@ typedef struct _complexStruct { /* complex number */ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots, int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float tsq, float Q, float f0, int reci, int *nshots_r, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose); int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose); -int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter); +int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter); void name_ext(char *filename, char *extension); @@ -464,7 +464,7 @@ int main (int argc, char **argv) tsyn += t3 - t2; if (file_iter != NULL) { - writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, iter); + writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, iter); } /* N_k(x,t) = -N_(k-1)(x,-t) */ /* p0^-(x,t) += iRN = (R * T_d^inv)(t) */ diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c index 7942e75..80c2398 100644 --- a/marchenko/marchenko_primaries.c +++ b/marchenko/marchenko_primaries.c @@ -29,7 +29,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose); -int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter); +int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter); int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm); int readData(FILE *fp, float *data, segy *hdrs, int n1); int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2); @@ -102,7 +102,7 @@ int main (int argc, char **argv) { FILE *fp_out, *fp_rr, *fp_w; size_t nread; - int i, j, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath; + int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath; int size, n1, n2, ntap, tap, di, ntraces; int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn; int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft; @@ -114,9 +114,9 @@ int main (int argc, char **argv) double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii; float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc; float *G_d, *DD, *RR, dt, dx, dxs, scl, mem; - float *rtrace, *tmpdata, *f1min, *iRN, *Ni, *trace; + float *rtrace, *tmpdata, *f1min, *f1plus, *iRN, *Ni, *trace; float xmin, xmax, scale, tsq; - float Q, f0, *ixmask; + float Q, f0, *ixmask, *costaper; complex *Refl, *Fop, *ctrace, *cwave; char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter; segy *hdrs_out, hdr; @@ -157,9 +157,17 @@ int main (int argc, char **argv) if(!getparint("shift", &shift)) shift=20; if(!getparint("smooth", &smooth)) smooth = 5; if(!getparint("ishot", &ishot)) ishot=300; - if (reci && ntap) vwarn("tapering influences the reciprocal result"); + smooth = MIN(smooth, shift); + if (smooth) { + costaper = (float *)malloc(smooth*sizeof(float)); + scl = M_PI/((float)smooth); + for (i=0; i f_1^-(t) */ + if (iter % 2 == 0) { /* even iterations: => f_1^+(t) */ /* apply muting for the acausal part */ for (l = 0; l < Nfoc; l++) { for (i = 0; i < npos; i++) { + /* apply mute window for samples after ii */ for (j = ii-T*shift; j < nts; j++) { Ni[l*nxs*nts+i*nts+j] = 0.0; } - for (j = 0; j < shift; j++) { + for (j = ii-T*shift+smooth, k=0; j < ii-T*shift; j++, k++) { + Ni[l*nxs*nts+i*nts+j] *= costaper[k]; + } + /* apply mute window for delta function at t=0*/ + for (j = 0; j < shift-smooth; j++) { Ni[l*nxs*nts+i*nts+j] = 0.0; } + for (j = shift-smooth, k=1; j < shift; j++, k++) { + Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k]; + } + f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j]; + for (j = 1; j < nts; j++) { + f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j]; + } } } + if (file_iter != NULL) { + writeDataIter("f1plus.su", f1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter); + } } - else {/* odd iterations: => f_1^+(t) */ + else {/* odd iterations: => f_1^-(t) */ for (l = 0; l < Nfoc; l++) { for (i = 0; i < npos; i++) { ix = ixpos[i]; @@ -602,17 +647,30 @@ int main (int argc, char **argv) f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j]; } } - /* apply mute window */ - for (j = nts-shift; j < nts; j++) { - Ni[l*nxs*nts+i*nts+j] = 0.0; - } + /* apply mute window for delta function at t=0*/ + for (j = nts-shift+smooth; j < nts; j++) { + Ni[l*nxs*nts+i*nts+j] = 0.0; + } + for (j = nts-shift, k=0; j < nts-shift+smooth; j++, k++) { + Ni[l*nxs*nts+i*nts+j] *= costaper[k]; + } /* apply mute window for samples above nts-ii */ - for (j = 0; j < nts-ii+T*shift; j++) { - Ni[l*nxs*nts+i*nts+j] = 0.0; - } + /* important to apply this mute after updating f1min */ + //for (j = 0; j < nts-ii+T*shift; j++) { + // Ni[l*nxs*nts+i*nts+j] = 0.0; + //} + /* apply mute window for samples above nts-ii */ + for (j = 0; j < nts-ii+T*shift-smooth; j++) { + Ni[l*nxs*nts+i*nts+j] = 0.0; + } + for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) { + Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k]; + } } } - //fprintf(stderr,"f1min at %d = %f NI=%f\n", ii, f1min[(nxs/2)*nts+ii], Ni[(nxs/2)*nts+ii]); + if (file_iter != NULL) { + writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter); + } } /* end else (iter) branch */ t2 = wallclock_time(); @@ -627,13 +685,13 @@ int main (int argc, char **argv) } } - /* To Do optional write intermediate RR results to file */ + /* To Do? optional write intermediate RR results to file */ if (verbose) { if(!((iend-ii-istart)%perc)) fprintf(stderr,"\b\b\b\b%3d%%",(ii-istart)*100/(iend-istart)); if((ii-istart)==10)t4=wallclock_time(); - if((ii-istart)==50){ - t4=(wallclock_time()-t4)*((iend-istart)/40.0); + if((ii-istart)==20){ + t4=(wallclock_time()-t4)*((iend-istart)/10.0); fprintf(stderr,"\r %s: Estimated total compute time = %.2fs.\n %s: Progress: %3d%%",xargv[0],t4,xargv[0],(ii-istart)/((iend-istart)/100)); } //t4=wallclock_time(); @@ -645,6 +703,8 @@ int main (int argc, char **argv) free(Ni); free(G_d); + free(f1min); + free(f1plus); t2 = wallclock_time(); if (verbose) { diff --git a/marchenko/writeDataIter.c b/marchenko/writeDataIter.c index 5de0b5d..b1a51f3 100644 --- a/marchenko/writeDataIter.c +++ b/marchenko/writeDataIter.c @@ -17,7 +17,7 @@ void name_ext(char *filename, char *extension); int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2); -int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter) +int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter) { FILE *fp_iter; size_t nwrite; @@ -43,14 +43,19 @@ int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, floa hdrs[i].selev = NINT(zsyn[l]*1000); hdrs[i].sdepth = NINT(-zsyn[l]*1000); /* rotate to get t=0 in the middle */ - hdrs[i].f1 = -n1*0.5*hdrs[i].d1; - memcpy(&trace[0],&data[l*size+ix*n1],n1*sizeof(float)); - for (j = 0; j < n1/2; j++) { - trace[n1/2+j] = data[l*size+ix*n1+j]; - } - for (j = n1/2; j < n1; j++) { - trace[j-n1/2] = data[l*size+ix*n1+j]; - } + if (t0shift) { + hdrs[i].f1 = -n1*0.5*hdrs[i].d1; + for (j = 0; j < n1/2; j++) { + trace[n1/2+j] = data[l*size+ix*n1+j]; + } + for (j = n1/2; j < n1; j++) { + trace[j-n1/2] = data[l*size+ix*n1+j]; + } + } + else { + hdrs[i].f1 = 0.0; + memcpy(&trace[0],&data[l*size+ix*n1],n1*sizeof(float)); + } nwrite = fwrite(&hdrs[i], 1, TRCBYTES, fp_iter); assert(nwrite == TRCBYTES); nwrite = fwrite(trace, sizeof(float), n1, fp_iter);