From 5643c815311563487e6d94705f5067b381839fa3 Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Thu, 12 Nov 2020 16:00:23 +0100
Subject: [PATCH] preparing for GPY MME paper publication
---
marchenko/demo/mme/clean | 4 ++
marchenko/demo/mme/epsModel.scr | 8 +--
marchenko/demo/mme/epsPrimaries.scr | 20 +++++---
marchenko/marchenko_primaries.c | 78 +++++++++++++++++++++++------
4 files changed, 83 insertions(+), 27 deletions(-)
create mode 100755 marchenko/demo/mme/clean
diff --git a/marchenko/demo/mme/clean b/marchenko/demo/mme/clean
new file mode 100755
index 0000000..3890128
--- /dev/null
+++ b/marchenko/demo/mme/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci
+
diff --git a/marchenko/demo/mme/epsModel.scr b/marchenko/demo/mme/epsModel.scr
index d71acdc..82d2785 100755
--- a/marchenko/demo/mme/epsModel.scr
+++ b/marchenko/demo/mme/epsModel.scr
@@ -21,13 +21,13 @@ cat << EOF3 > line3
EOF3
#model
-supsimage hbox=4 wbox=6 labelsize=12 < model10_cp.su \
+supsimage hbox=4 wbox=6 labelsize=16 < model10_cp.su \
x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \
label1="depth (m)" label2="lateral distance (m)" \
curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_cp_line.eps
-supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \
+supsimage hbox=4 wbox=6 labelsize=16 < model10_ro.su \
x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \
label1="depth (m)" label2="lateral distance (m)" \
curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
@@ -36,12 +36,12 @@ supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \
#wavelet
dt=0.0005
supsgraph < wavefw.su \
- labelsize=12 d1=$dt style=normal \
+ labelsize=16 d1=$dt style=normal \
label1="time (s)" label2="amplitude" \
d1num=0.15 wbox=6 hbox=3 x1end=0.9 > wavefw.eps
sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \
- labelsize=12 style=normal \
+ labelsize=16 style=normal \
label1="frequency (1/s)" label2="amplitude" \
d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps
diff --git a/marchenko/demo/mme/epsPrimaries.scr b/marchenko/demo/mme/epsPrimaries.scr
index 67341c2..e74cf81 100755
--- a/marchenko/demo/mme/epsPrimaries.scr
+++ b/marchenko/demo/mme/epsPrimaries.scr
@@ -102,6 +102,10 @@ fi
# second iteration
if [[ "$1" == "Figure4" ]];
then
+if [[ ! -f "shotsx0.su" ]]; then
+echo "ERR: File shotsx0.su is not yet created, please run ./epsPrimaries.scr Figure3 first."
+exit
+fi
suwind itmax=1023 < Mi_276001.su > N0.su
#compute R*N0
fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=90
@@ -182,7 +186,7 @@ echo $piter
file=Mi_${istart}${piter}.su
file_base=${file%.su}
#possibly add suflip flip=3 to flip the time-axis
-supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+supsimage < $file hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
@@ -197,7 +201,7 @@ piter=$(printf %03d $iter)
echo $piter
file=Mi_${istart}${piter}.su
file_base=${file%.su}
-supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+supsimage < $file hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
@@ -218,14 +222,14 @@ file=Mi_${istart}${piter}.su
#shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
#basop choice=shift shift=$shift file_in=$file | \
file_base=${file%.su}
-suflip flip=3 < $file | supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+suflip flip=3 < $file | supsimage hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > ${file_base}flip.eps
file=k1min_${istart}${piter}.su
file_base=${file%.su}
-supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
+supsimage hbox=6 wbox=4 labelsize=16 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
@@ -242,7 +246,7 @@ for (( istart=246; istart<=316; istart+=10 ))
do
file=k1min_${istart}${piter}.su
file_base=${file%.su}
-supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
+supsimage hbox=6 wbox=4 labelsize=16 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
@@ -259,7 +263,7 @@ echo $piter
file=Mi_${istart}${piter}.su
file_base=${file%.su}
#possibly add suflip flip=3 to flip the time-axis
-supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+supsimage < $file hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}T.eps
@@ -296,12 +300,12 @@ file=WindowOdd276.su
file_base=${file%.su}
suwind key=tracl min=1 max=1 < $file | suwind itmin=1024 | \
supsgraph d1=1 style=normal f1=0 \
- labelsize=12 label1="time sample number" label2="amplitude" \
+ labelsize=14 label1="time sample number" label2="amplitude" \
d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps
file=WindowEven276.su
file_base=${file%.su}
suwind key=tracl min=1 max=1 < $file | suwind itmax=1024 | \
supsgraph d1=1 style=normal f1=0 \
- labelsize=12 label1="time sample number" label2="amplitude" \
+ labelsize=14 label1="time sample number" label2="amplitude" \
d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index 12ef798..0e4843c 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -109,6 +109,13 @@ char *sdoc[] = {
NULL};
/**************** end self doc ***********************************/
+/*
+ * The equation and Figure numbers mentioned in the comments refer to the paper:
+ * "Implementation of the Marchenko Multiple Elimination algorithm",
+ * Thorbecke, Jan; Zhang, Lele; Wapenaar, Kees; Slob, Evert,
+ * Accepted for publication, todo add full reference after publication
+ */
+
int main (int argc, char **argv)
{
FILE *fp_out, *fp_rr, *fp_w, *fp_ud;
@@ -197,7 +204,7 @@ int main (int argc, char **argv)
if (reci && ntap) vwarn("tapering influences the reciprocal result");
-/* defines the smooth transition zone of the time-window */
+/* defines the smooth transition zone of the time-window (Figure 5) */
smooth = MIN(smooth, shift);
if (smooth) {
@@ -431,7 +438,7 @@ int main (int argc, char **argv)
}
/*================ Defining focusing operator(s) from R ================*/
-/* M0 = -R(ishot,-t) */
+/* M0 = -R(ishot,-t) equation (3) */
/* use ishot from Refl, complex-conjugate(time reverse), scale with -1 and convolve with wavelet */
if (file_tinv == NULL) {
@@ -702,7 +709,7 @@ int main (int argc, char **argv)
for (j = 0; j < nts; j++) {
M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j];
}
- /* apply mute window for samples above nts-ii */
+ /* apply mute window for samples above nts-ii : the Heaviside function in equation (3) */
for (j = 0; j < MIN(nts, nts-iw+isms); j++) {
M0[l*nxs*nts+i*nts+j] = 0.0;
}
@@ -710,6 +717,7 @@ int main (int argc, char **argv)
M0[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
}
+ /* Initilisation of k1min: the first term (R) in the right hand side of equation(6) */
for (i = 0; i < npos; i++) {
ix = ixpos[i];
j = 0;
@@ -743,6 +751,7 @@ int main (int argc, char **argv)
}
}
}
+ /* If enabled this write the result of equation(3) to disk */
if (file_iter != NULL) {
writeDataIter("M0.su", M0, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii);
}
@@ -760,15 +769,17 @@ int main (int argc, char **argv)
t2 = wallclock_time();
/*================ construction of Mi(-t) = - \int R(x,t) Mi(t) ================*/
- /* synthesis process is the compute kernel in equations (4) and (5) in the Geophysics implementation paper */
+ /* synthesis process is the compute kernel in equations (4) and (5) (rewritten from equation (2)) */
synthesisp(Refl, Fop, Mi, RMi, nx, nt, nacq, nts, dt, xsyn, Nfoc,
xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
+ /* If enbled the result of equation(4) or (5) (depending on odd/even) is written to disk */
if (file_iter != NULL) {
writeDataIter(file_iter, RMi, hdrs_out, ntfft, nacq, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
}
+ /* The Energy is computed to check the convergence of the scheme see Figure 6c and Figure 8 */
if (verbose >=2) {
for (l = 0; l < Nfoc; l++) {
if (iter % 2 == 0) {
@@ -790,7 +801,8 @@ int main (int argc, char **argv)
t3 = wallclock_time();
tsyn += t3 - t2;
- /* N_k(x,t) = -N_(k-1)(x,-t) */
+ /* M_i(x,t) = RMi_(i-1)(x,-t) : time-reversal of the outcome of the synthesis process */
+ /* see equation(3) where the time reversal of R is taken as input */
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
j = 0;
@@ -808,18 +820,19 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) {
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
- /* store results in kplus */
+ /* store results in kplus the 'plus-part' of equation (6) */
for (j = 0; j < nts; j++) {
k1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
}
- /* apply mute window for samples after ii */
+ /* Apply mute window for samples after ii : the Heaviside function in equation (5) */
+ /* This also defines v1_plus in equation (8) */
for (j = MAX(0,iw-isme); j < nts; j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = MAX(0,iw-isme), k=0; j < iw-isms; j++, k++) {
Mi[l*nxs*nts+i*nts+j] *= costaper[k];
}
- /* apply mute window for delta function at t=0*/
+ /* Apply mute window for delta function at t=0*/
iw = NINT((twplane[ix])/dt);
for (j = 0; j < MAX(0,iw+shift-smooth); j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0;
@@ -827,6 +840,7 @@ int main (int argc, char **argv)
for (j = MAX(0,iw+shift-smooth), k=1; j < MAX(0,iw+shift); j++, k++) {
Mi[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
+ /* This defines v1_plus in equation (8) */
for (j = 0; j < nts; j++) {
v1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
}
@@ -844,11 +858,13 @@ int main (int argc, char **argv)
for (j = 0; j < nts; j++) {
Mi[l*nxs*nts+i*nts+j] += -DD[l*nxs*nts+ix*nts+j];
}
+ /* store results in kmin as defined in equation (6) */
j = 0;
k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j];
}
+ /* This collects the eliminated internal multiples and is the sum in right hand-side in equation (1) */
if (file_update != NULL) {
j=0;
Mup[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j] - k1min[l*nxs*nts+i*nts+j];
@@ -858,11 +874,13 @@ int main (int argc, char **argv)
}
}
else {
+ /* store results in kmin as defined in equation (6) */
j = 0;
k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
}
+ /* Mup collects the eliminated internal multiples and is the sum in right hand-side in equation (1) */
if (file_update != NULL) {
j=0;
Mup[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
@@ -871,7 +889,8 @@ int main (int argc, char **argv)
}
}
}
- /* apply mute window for samples above nts-ii */
+ /* Apply mute window for samples above nts-ii : the Heaviside function in equation (4)
+ * This also defines u1_min, for t>=t2-epsilon in equation (7) */
iw = NINT((ii*dt+twplane[ix])/dt);
for (j = 0; j < MIN(nts,nts-iw+isms); j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0;
@@ -907,7 +926,8 @@ int main (int argc, char **argv)
}
}
- /* write optional u- u+ v- v+ to disk */
+ /* Write optional u- u+ v- v+ to disk */
+ /* Based on equation(6) and (7) k1min is split in vmin and umin */
if (file_vmin != NULL) {
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
@@ -952,16 +972,44 @@ int main (int argc, char **argv)
}
writeDataIter(file_umin, uv, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, ii);
}
+ /*
+ * \begin{alignat}{2}
+ * \label{k1p}
+ * k_{1,i}^+({\bf x}'_0,{\bf x}''_0,t,t_2) = & \int_{t'=0}^{+\infty} \int_{\partial \setD_0} R({\bf x}'_0,{\bf x}_0,-t')
+ * v_{1,i}^-({\bf x}_0,{\bf x}''_0,t-t',t_2) {\rm d}t' {\rm d} {\bf x}_0.
+ * \end{alignat}
+ *
+ * Equation (\ref{k1p}) can be further split in the time domain as follows;
+ * %
+ * \begin{equation} \label{uv+}
+ * k_{1,i}^+({\bf x}'''_0,{\bf x}''_0,t,t_2) = \begin{cases}
+ * v_{1,i}^+({\bf x}'''_0,{\bf x}''_0,t,t_2) & t < t_2-\varepsilon \\
+ * u_{1,i}^+({\bf x}'''_0,{\bf x}''_0,t,t_2) & t \geq t_2-\varepsilon
+ * \end{cases},
+ * \end{equation}
+ * %
+ *
+ */
+ /* The equation (\ref{k1p}) above is the k1plus equivalent of k1min of equation(6)
+ * and the time window in equation (\ref{uv+}) defines vplus and uplus */
+ /* Based on these equations and k1plus is split in vplus and uplus */
+
if (file_vplus != NULL) {
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
/* apply mute window for delta function at t=0*/
iw0 = NINT((twplane[ix])/dt);
- for (j = 0; j < MAX(0,iw0+shift-smooth); j++) {
- uv[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = MAX(0,iw0+shift-smooth), k=1; j < MAX(0,iw0+shift); j++, k++) {
- uv[l*nxs*nts+i*nts+j] = k1plus[l*nxs*nts+i*nts+j]*costaper[smooth-k];
+ /* to apply muting window around epsilon below t=0 uncomment this part */
+ /* This (partly) removes the delta function at t=0 */
+ // for (j = 0; j < MAX(0,iw0+shift-smooth); j++) {
+ // uv[l*nxs*nts+i*nts+j] = 0.0;
+ // }
+ // for (j = MAX(0,iw0+shift-smooth), k=1; j < MAX(0,iw0+shift); j++, k++) {
+ // uv[l*nxs*nts+i*nts+j] = k1plus[l*nxs*nts+i*nts+j]*costaper[smooth-k];
+ // }
+ /* this keeps the delta pulse below at at t=0 */
+ for (j = 0; j < iw0+shift; j++) {
+ uv[l*nxs*nts+i*nts+j] = k1plus[l*nxs*nts+i*nts+j];
}
iw = NINT((ii*dt+twplane[ix])/dt);
for (j = MAX(0,iw0+shift); j < iw-isme; j++, k++) {
--
GitLab