From c85c18fd4dcae921b320f9b3469af02c1fc9307f Mon Sep 17 00:00:00 2001
From: Jan Thorbecke <janth@xs4all.nl>
Date: Thu, 5 Jul 2018 16:12:40 +0200
Subject: [PATCH] added smoothing on 4x oversampled grid to reduce straircase
effect in FD
---
utils/makemod.c | 95 ++++++++++++++++++++++++++++++++++++++++++++++---
1 file changed, 91 insertions(+), 4 deletions(-)
diff --git a/utils/makemod.c b/utils/makemod.c
index 83e85f7..9e3a93d 100644
--- a/utils/makemod.c
+++ b/utils/makemod.c
@@ -86,6 +86,7 @@ char *sdoc[] = {
" dtype=0 .................. diffractor type for diffr and randdf",
" OUTPUT ",
" writeint=0 ............... interfaces as function of x (ext: _int)",
+ " supersmooth=0 ............ samples at dx/4 and apply 5x5 smoothing operator(1)",
" reflectivity=0 ........... compute reflectivity (ext: _rfl)",
" rayfile=0 ................ interfaces as function of x in ASCII file.mod",
" skip=5 ................... number of x position to skip in file_ray",
@@ -121,6 +122,7 @@ char *sdoc[] = {
" Option for gradunit, gradient unit per layer:",
" - 0 = gradient unit per layer is m/s per dz (default)",
" - 1 = gradient unit per layer is m/s per m",
+ " (1) Reference: Zeng and West: Geophysics 1996 ",
" ",
" makemod builds a gridded subsurface file which can be used in a migration",
" or finite difference program. The gridded model is stored in files with",
@@ -139,13 +141,13 @@ int main(int argc, char **argv)
{
FILE *fpint, *fpcp, *fpcs, *fpro;
int example, verbose, writeint, nb;
- int above, diffrwidth, dtype, reflectivity;
+ int above, diffrwidth, dtype, reflectivity, supersmooth;
int Ngp, Ngs, Ngr, Np, Ns, Nr, Ng, Ni, Nv, Nvi, No, Noi;
int jint, jcount, j, ix, iz, nx, nz, nxp, nzp, *zp, nmaxx, nminx, optgrad, poly, gradt;
int ncp, nro, ncs, nvel, skip, rayfile, store_int;
- long lseed;
+ long lseed;
size_t nwrite;
- float *data_out, orig[2], cp0, cs0, ro0, Z1, Z2;
+ float *data_out, orig[2], cp0, cs0, ro0, Z1, Z2, sigma;
float *x, *z, *var, *interface, **inter;
float back[3], sizex, sizez, dx, dz;
float **cp ,**cs, **ro, aver, gradlen, gradcp, gradcs, gradro;
@@ -212,12 +214,18 @@ int main(int argc, char **argv)
if(!getparint("gradunit", &gradunit)) gradunit = 0;
if(!getparint("writeint", &writeint)) writeint = 0;
if(!getparint("reflectivity", &reflectivity)) reflectivity = 0;
+ if(!getparint("supersmooth", &supersmooth)) supersmooth = 0;
+ if(!getparfloat("sigma", &sigma)) sigma = 0.0;
if(!getparint("rayfile", &rayfile)) rayfile = 0;
if(!getparint("skip", &skip)) skip = 5;
if(!getparint("above", &above)) above=0;
if(!getparint("verbose", &verbose)) verbose=0;
if(!getparint("dtype", &dtype)) dtype = 0;
+ if (supersmooth==1) {
+ dx = dx/4;
+ dz = dz/4;
+ }
if ((writeint == 1) || (rayfile == 1)) store_int = 1;
else store_int = 0;
@@ -631,6 +639,85 @@ int main(int argc, char **argv)
}
} /* end of loop over interfaces */
+/* apply supersampled data smoothing operator according to Zeng and West (1996) */
+
+ if (supersmooth==1) {
+ float Wsmooth[5][5], C, iC, xx, xz, *dataS, smooth;
+ int ixi, izi, nxout, nzout;
+
+ sigma = -1.0*log(sigma)/(dx*(powf(0.25*2.0,2.0)));
+ for(ixi = -2; ixi < 3; ixi++) {
+ for(izi = -2; izi < 3; izi++) {
+ xx = 0.25*ixi;
+ xz = 0.25*izi;
+ Wsmooth[ixi+2][izi+2] = exp(-sigma*dx*(xx*xx+xz*xz) );
+ //fprintf(stderr,"Wsmooth[%d][%d] = %f\n", ixi, izi, Wsmooth[ixi+2][izi+2]);
+ C += Wsmooth[ixi+2][izi+2];
+ }
+ }
+ iC = 1.0/C;
+ //fprintf(stderr,"sigma=%f %f C=%f\n",sigma, exp(-sigma*dx*powf(0.25*2.0,2.0)),C);
+
+ nxout = (nx-1)/4+1;
+ nzout = (nz-1)/4+1;
+ fprintf(stderr,"nxout=%d nzout=%d nx=%d nz=%d \n",nxout,nzout, nx, nz);
+ dataS = malloc(nxout*nzout*sizeof(float));
+ for(ix = 0; ix < nxout-1; ix++) {
+ for(iz = 0; iz < nzout-1; iz++) {
+ smooth = 0.0;
+ for(ixi = -2; ixi < 3; ixi++) {
+ for(izi = -2; izi < 3; izi++) {
+ smooth += Wsmooth[ixi+2][izi+2]/gridcp[4*ix+2+ixi][4*iz+2+izi];
+ }
+ }
+ dataS[ix*nzout+iz] = 1.0/(iC*smooth);
+ }
+ iz = nzout-1;
+ dataS[ix*nzout+iz] = dataS[ix*nzout+iz-1];
+ }
+ ix = nxout-1;
+ for(iz = 0; iz < nzout; iz++) {
+ dataS[ix*nzout+iz] = dataS[(ix-1)*nzout+iz-1];
+ }
+
+ for(ix = 0; ix < nxout; ix++) {
+ for(iz = 0; iz < nzout; iz++) {
+ gridcp[ix][iz] = dataS[ix*nzout+iz];
+ }
+ }
+
+ /* smoothin densities */
+ if (Nr > 0 || getparfloat("ro0", &ro0)) {
+ for(ix = 0; ix < nxout-1; ix++) {
+ for(iz = 0; iz < nzout-1; iz++) {
+ smooth = 0.0;
+ for(ixi = -2; ixi < 3; ixi++) {
+ for(izi = -2; izi < 3; izi++) {
+ smooth += Wsmooth[ixi+2][izi+2]*gridro[4*ix+2+ixi][4*iz+2+izi];
+ }
+ }
+ dataS[ix*nzout+iz] = iC*smooth;
+ }
+ iz = nzout-1;
+ dataS[ix*nzout+iz] = dataS[ix*nzout+iz-1];
+ }
+ ix = nxout-1;
+ for(iz = 0; iz < nzout; iz++) {
+ dataS[ix*nzout+iz] = dataS[(ix-1)*nzout+iz-1];
+ }
+ for(ix = 0; ix < nxout; ix++) {
+ for(iz = 0; iz < nzout; iz++) {
+ gridro[ix][iz] = dataS[ix*nzout+iz];
+ }
+ }
+ }
+ nx = nxout;
+ nz = nzout;
+ dx = dx*4;
+ dz = dz*4;
+ free(dataS);
+ }
+ fprintf(stderr,"nx=%d nz=%d \n", nx, nz);
if (verbose) vmess("Writing data to disk.");
hdrs = (segy *) calloc(nx,sizeof(segy));
@@ -657,7 +744,7 @@ int main(int argc, char **argv)
assert(fpint != NULL);
}
-
+
/* write P-velocities in file */
strcpy(filename, file_base);
name_ext(filename, "_cp");
--
GitLab