Skip to content

Commit 4143d34

Browse files
DirackDirack
authored
Add program sfvfsacrsnh and sfnhcrssurf
The program _sfvfsacrsnh_ gets the zero offset CRS parameters (RN, RNIP and BETA) using Very Fast Simulated Aneeling (VFSA) global optimization. The program _sfnhcrssurf_ builds a 2D CRS approximated surface using [non hyperbolic CRS traveltime approximation](http://www.reproducibility.org/RSF/book/tccs/crs/paper_html/) and the CRS parameters. The VFSA algorithm fits the [non hyperbolic CRS traveltime approximation](http://www.reproducibility.org/RSF/book/tccs/crs/paper_html/) in a modeled data cube. The CRS parameters that produce the best fit will be the optimal.
1 parent 481123a commit 4143d34

File tree

5 files changed

+644
-0
lines changed

5 files changed

+644
-0
lines changed

‎Mnhcrssurf.c

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
/* Version 1.0 - Build Non-Hyperbolic CRS approximation surface giver RN, RNIP and BETA parameters.
2+
3+
Programer: Rodolfo A. C. Neves (Dirack) 19/09/2019
4+
5+
Email: [email protected]
6+
7+
License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
8+
9+
*/
10+
11+
#include <math.h>
12+
#include <stdio.h>
13+
#include <stdlib.h>
14+
#include <rsf.h>
15+
16+
int main(int argc, char* argv[])
17+
{
18+
19+
float m0; // central CMP
20+
float om; // CMP axis origin
21+
float dm; // CMP sampling
22+
int nm; // Number of CMP's
23+
float oh; // Offset axis origin
24+
float dh; // Offset sampling
25+
int nh; // Number of Offsets
26+
bool verb; // Key to turn On/Off active mode
27+
float v0; // Near surface velocity
28+
float t0; // Normal ray time travel
29+
float *c; // Temporary parameters vector - last iteration
30+
float **t; // CRS surface t(m,h)
31+
float RN, RNIP, BETA; // CRS parameters
32+
float Fd, Fd1, Fd2;
33+
float c1, a1, a2, b2;
34+
float m;
35+
float h;
36+
int nc, ih, im;
37+
38+
/* RSF files I/O */
39+
sf_file in, out, par;
40+
41+
/* RSF files axis */
42+
sf_axis ax,ay,az;
43+
44+
sf_init(argc,argv);
45+
46+
in = sf_input("in");
47+
par = sf_input("param");
48+
out = sf_output("out");
49+
50+
if (!sf_getfloat("m0",&m0)) m0=0;
51+
/* central CMP of the approximation (Km) */
52+
53+
if (!sf_getfloat("v0",&v0)) v0=1.5;
54+
/* Near surface velocity (Km/s) */
55+
56+
if (!sf_getfloat("t0",&t0)) t0=1.5;
57+
/* Normal ray traveltime (s) */
58+
59+
if (!sf_histint(in,"n1",&nh)) sf_error("No n1= in input");
60+
if (!sf_histfloat(in,"d1",&dh)) sf_error("No d1= in input");
61+
if (!sf_histfloat(in,"o1",&oh)) sf_error("No o1= in input");
62+
if (!sf_histint(in,"n2",&nm)) sf_error("No n2= in input");
63+
if (!sf_histfloat(in,"d2",&dm)) sf_error("No d2= in input");
64+
if (!sf_histfloat(in,"o2",&om)) sf_error("No o2= in input");
65+
66+
if(!sf_histint(par,"n1",&nc)) sf_error("No n1= in parameters input");
67+
68+
if(! sf_getbool("verb",&verb)) verb=0;
69+
/* 1: active mode; 0: quiet mode */
70+
71+
if (verb) {
72+
73+
sf_warning("Active mode on!!!");
74+
sf_warning("Command line parameters: ");
75+
sf_warning("m0=%f v0=%f t0=%f",m0,v0,t0);
76+
sf_warning("Input file parameters: ");
77+
sf_warning("n1=%i d1=%f o1=%f",nh,dh,oh);
78+
sf_warning("n2=%i d2=%f o2=%f",nm,dm,om);
79+
sf_warning("Param file parameters: ");
80+
sf_warning("n1=%i",nc);
81+
}
82+
83+
c = sf_floatalloc(nc);
84+
sf_floatread(c,nc,par);
85+
86+
RN = c[0];
87+
RNIP = c[1];
88+
BETA = c[2];
89+
90+
91+
t = sf_floatalloc2(nh,nm);
92+
93+
for (im=0; im < nm; im++){
94+
95+
for(ih=0;ih<nh;ih++){
96+
97+
m = om + (im * dm);
98+
m = m - m0;
99+
h = oh + (ih * dh);
100+
101+
a1=(2*sin(BETA))/(v0);
102+
a2=(2*cos(BETA)*cos(BETA)*t0)/(v0*RN);
103+
b2=(2*cos(BETA)*cos(BETA)*t0)/(v0*RNIP);
104+
c1=2*b2+a1*a1-a2;
105+
106+
Fd=(t0+a1*m)*(t0+a1*m)+a2*m*m;
107+
Fd2=(t0+a1*(m-h))*(t0+a1*(m-h))+a2*(m-h)*(m-h);
108+
Fd1=(t0+a1*(m+h))*(t0+a1*(m+h))+a2*(m+h)*(m+h);
109+
t[im][ih]=sqrt((Fd+c1*h*h+sqrt(Fd2*Fd1))*0.5);
110+
111+
}
112+
}
113+
114+
/* axis = sf_maxa(n,o,d)*/
115+
ax = sf_maxa(nh, oh, dh);
116+
ay = sf_maxa(nm, om, dm);
117+
az = sf_maxa(1, 0, 1);
118+
119+
/* sf_oaxa(file, axis, axis index) */
120+
sf_oaxa(out,ax,1);
121+
sf_oaxa(out,ay,2);
122+
sf_oaxa(out,az,3);
123+
sf_floatwrite(t[0],nh*nm,out);
124+
125+
exit(0);
126+
}

‎Mvfsacrsnh.c

Lines changed: 201 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,201 @@
1+
/* Version 1.0 - Zero offset CRS parameter inversion (RN, RNIP, BETA) with Very Fast Simulated Aneeling (VFSA) Global Optimization
2+
3+
This program the Non-Hyperbolic CRS approximation to fit data cube and get the parameters (Fomel, 2013).
4+
5+
Programer: Rodolfo A. C. Neves (Dirack) 19/09/2019
6+
7+
Email: [email protected]
8+
9+
License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
10+
11+
*/
12+
13+
#include "vfsacrsnh_lib.h"
14+
15+
int main(int argc, char* argv[])
16+
{
17+
18+
float m0; // central CMP
19+
float om; // CMP axis origin
20+
float dm; // CMP sampling
21+
int nm; // Number of CMP's
22+
float oh; // Offset axis origin
23+
float dh; // Offset sampling
24+
int nh; // Number of Offsets
25+
int nt; // Number of time samples
26+
float ot; // Time axis origin
27+
float dt; // Time sampling
28+
bool verb; // Key to turn On/Off active mode
29+
float v0; // Near surface velocity
30+
float t0; // Normal ray time travel
31+
float cnew[3]; // Temporary parameters vector - actual iteration
32+
float c[3]; // Temporary parameters vector - last iteration
33+
float *otm; // Optimazed parameters
34+
float otrn, otrnip, otbeta, otsemb; // Optimazed parameters - actual iteration
35+
float deltaE, PM; // Metrópolis criteria
36+
float Em0=0; // Major semblance
37+
float u; // Random number
38+
float ***t; // Data cube A(m,h,t)
39+
int q, i; // loop counter
40+
float semb; // Semblance - actual iteration
41+
float RN, RNIP, BETA; // CRS parameters
42+
float semb0; // Inicial semblance value
43+
float c0; // VFSA damping factor
44+
float temp0; // inicial VFSA temperature
45+
float temp; // VFSA temperature
46+
int repeat; // perform VFSA optimization more than once
47+
48+
/* RSF files I/O */
49+
sf_file in, out;
50+
51+
/* RSF files axis */
52+
sf_axis ax,ay,az;
53+
54+
sf_init(argc,argv);
55+
56+
in = sf_input("in");
57+
out = sf_output("out");
58+
59+
if (!sf_getfloat("m0",&m0)) m0=0;
60+
/* central CMP of the approximation (Km) */
61+
62+
if (!sf_getfloat("v0",&v0)) v0=1.5;
63+
/* Near surface velocity (Km/s) */
64+
65+
if (!sf_getfloat("t0",&t0)) t0=1.5;
66+
/* Normal ray traveltime (s) */
67+
68+
if (!sf_getfloat("c0",&c0)) c0=0.5;
69+
/* damping factor of VFSA */
70+
71+
if (!sf_getfloat("temp0",&temp0)) temp0=10;
72+
/* inicial VFSA temperature */
73+
74+
if(!sf_getint("repeat",&repeat)) repeat=1;
75+
/* How many times to perform VFSA global optimization */
76+
77+
if (!sf_histint(in,"n1",&nt)) sf_error("No n1= in input");
78+
if (!sf_histfloat(in,"d1",&dt)) sf_error("No d1= in input");
79+
if (!sf_histfloat(in,"o1",&ot)) sf_error("No o1= in input");
80+
if (!sf_histint(in,"n2",&nh)) sf_error("No n2= in input");
81+
if (!sf_histfloat(in,"d2",&dh)) sf_error("No d2= in input");
82+
if (!sf_histfloat(in,"o2",&oh)) sf_error("No o2= in input");
83+
if (!sf_histint(in,"n3",&nm)) sf_error("No n3= in input");
84+
if (!sf_histfloat(in,"d3",&dm)) sf_error("No d3= in input");
85+
if (!sf_histfloat(in,"o3",&om)) sf_error("No o3= in input");
86+
87+
if(! sf_getbool("verb",&verb)) verb=0;
88+
/* 1: active mode; 0: quiet mode */
89+
90+
if (verb) {
91+
92+
sf_warning("Active mode on!!!");
93+
sf_warning("Command line parameters: ");
94+
sf_warning("m0=%f v0=%f t0=%f c0=%f temp0=%f repeat=%i",m0,v0,t0,c0,temp0,repeat);
95+
sf_warning("Input file parameters: ");
96+
sf_warning("n1=%i d1=%f o1=%f",nt,dt,ot);
97+
sf_warning("n2=%i d2=%f o2=%f",nh,dh,oh);
98+
sf_warning("n3=%i d3=%f o3=%f",nm,dm,om);
99+
}
100+
101+
c[0] = 0;
102+
c[1] = 0;
103+
c[2] = 0;
104+
cnew[0] = 0;
105+
cnew[1] = 0;
106+
cnew[2] = 0;
107+
108+
srand(time(NULL));
109+
110+
/* Read seismic data cube */
111+
t=sf_floatalloc3(nt,nh,nm);
112+
sf_floatread(t[0][0],nm*nh*nt,in);
113+
114+
sf_fileclose(in);
115+
116+
semb0=0;
117+
118+
for(i=0;i<repeat;i++){
119+
120+
for (q=0; q <ITMAX; q++){
121+
122+
/* calculate VFSA temperature for this iteration */
123+
temp=getVfsaIterationTemperature(q,c0,temp0);
124+
125+
/* parameter disturbance */
126+
disturbParameters(temp,cnew,c);
127+
128+
RN = cnew[0];
129+
RNIP = cnew[1];
130+
BETA = cnew[2];
131+
132+
semb=0;
133+
134+
/* Calculate semblance: Non-hyperbolic CRS approximation with data */
135+
semb=semblance(m0,dm,om,oh,dh,dt,nt,t0,v0,RN,RNIP,BETA,t);
136+
137+
138+
/* VFSA parameters convergence condition */
139+
if(fabs(semb) > fabs(semb0) ){
140+
otsemb = semb;
141+
otrn = RN;
142+
otrnip = RNIP;
143+
otbeta = BETA;
144+
semb0 = semb;
145+
}
146+
147+
/* VFSA parameters update condition */
148+
deltaE = -semb - Em0;
149+
150+
/* Metrópolis criteria */
151+
PM = expf(-deltaE/temp);
152+
153+
if (deltaE<=0){
154+
c[0] = cnew[0];
155+
c[1] = cnew[1];
156+
c[2] = cnew[2];
157+
Em0 = -semb;
158+
} else {
159+
u=getRandomNumberBetween0and1();
160+
if (PM > u){
161+
c[0] = cnew[0];
162+
c[1] = cnew[1];
163+
c[2] = cnew[2];
164+
Em0 = -semb;
165+
}
166+
}
167+
168+
} /* loop over iterations */
169+
170+
} /* repeat VFSA global optimization */
171+
172+
free(t);
173+
174+
/* Save optimized parameters in 'param' file */
175+
otm=sf_floatalloc(8);
176+
otm[0] = otrn;
177+
otm[1] = otrnip;
178+
otm[2] = otbeta;
179+
otm[3] = otsemb;
180+
otm[4] = c0;
181+
otm[5] = temp0;
182+
otm[6] = t0;
183+
otm[7] = m0;
184+
185+
/* Show optimized parameters on screen before save them */
186+
sf_warning("Parâmetros otimizados:\n RN=%f, RNIP=%f, BETA=%f, SEMB=%f",otrn,otrnip,otbeta,otsemb);
187+
188+
/* axis = sf_maxa(n,o,d)*/
189+
ax = sf_maxa(8, 0, 1);
190+
ay = sf_maxa(1, 0, 1);
191+
az = sf_maxa(1, 0, 1);
192+
193+
/* sf_oaxa(file, axis, axis index) */
194+
sf_oaxa(out,ax,1);
195+
sf_oaxa(out,ay,2);
196+
sf_oaxa(out,az,3);
197+
sf_floatwrite(otm,8,out);
198+
199+
sf_close();
200+
exit(0);
201+
}

0 commit comments

Comments
 (0)