-
Notifications
You must be signed in to change notification settings - Fork 1
/
final_real_Chai1988_2a.m
300 lines (269 loc) · 11.2 KB
/
final_real_Chai1988_2a.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
% ***************************************************************
% *** Matlab function for real model-2 with varying density is a part of SPoDEA programe that includes a set of *.m files to compute basement depth of the complex sedimentary basin.
% *** Source Code is mainly written for research purposes. The codes are
% *** having copyrights and required proper citations whenever it is used.
% *** Originated by:
% *** Mr. Arka Roy (email: [email protected])
% *** Dr. Chandra Prakash Dubey (email:[email protected])
% *** Mr. M. Prasad (email:[email protected])
% *** Crustal Processes Group, National Centre for Earth Science Studies,
% *** Ministry of Earth Sciences, Government of India
% *** Thiruvanthapuram, Kerala, India
% ****************************************************************
%Real example of San Jacinto Graben for gravity field and inversion
clear all
close all
%%%%%%%%%%%%%%%%%%%%%%%% Part 1. %%%%%%%%%%%%%%%%%%%%%%%%
%%%% Gravity anomaly of San Jacinto Graben and depth profile from Chai Hinze, 1988 %%%%
%Estimated depth profile of San Jacinto Graben from Chai Hinze, 1988
depth=(importdata('depth_Chai1988_2a.dat'))';
x_obs=(importdata('x_obs_Chai1988_2a.dat'))';
%Observed gravity anomaly of San Jacinto Graben from Chai Hinze, 1988
zz1=(importdata('gravity_anomaly_Chai1988_2a.dat'))';
%t and c are Legendre Gaussian quadrature integral points.
[t_leg,c_leg]=lgwt(10,0,1);
%plotting of depth profile estimated by Chai Hinze, 1988
figure(1)
subplot(2,1,2)
plot(x_obs,depth)
set(gca,'Ydir','reverse')
xlabel('Distance in meter')
ylabel('Depth in meter')
title('Depth profile of the San Jacinto Graben')
%Plotting the observed Gravity field anomaly
figure(1)
subplot(2,1,1)
plot(x_obs,zz1)
xlabel('Observation points in meter')
ylabel('Gravity anomaly in mGal')
title('Gravity anomaly for San Jacinto Graben ')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% Part 2. %%%%%%%%%%%%%%%%%%%%%%%%
%%%% Inversion of Gravity anomaly using Differential Evolution %%%%
data=zz1; %data is gravity field for San Jacinto Graben
%Finding inversion of the same gravity profile for varying density
%Parameters for DE algorithm
Cross_rate=0.8; %Crossover Rate
nVar = 14; %Number of Genes
K= 0.5; F=.8; %Scalling and Combinition Factor
%Parameters for b-Spline
n=13; % number of knots
k=4; % order of polynomial
%Population size
nPoP = 30; %Number of Chromosome
%Total number of Generations
MaxIt = 1000; %Number of Generations
%%Problem Definition
%Objective Function
CostFunction =@(x,data) myCostFunction(x,data,n,k,t_leg,c_leg)+1000*(constrained1(x,n,k)+constrained2(x,n,k));
VarSize = [1 nVar]; %Matrix size of Decision variables
VarMin= -500*ones(1,nVar); %Lower Bound of Unknown variable
VarMax= 500*ones(1,nVar); %Upper Bound of Unknown variable
%%Initialization
Empty.Particle.Position =[]; %creating Empty vectors for each Cromosome
Empty.Particle.Cost =[]; %creating Empty vectors for cost of Cromosome
Particle=repmat(Empty.Particle, nPoP,1); %Empty Matrix for all cromosome
Vector=Particle; %Cromosome
Next_Particle=Particle;
%Loop for initialization of Population
for i=1:nPoP
%initialize position with random number from VarMin and VarMax
for j=1:nVar
Particle(i).Position(j) =unifrnd(VarMin(j),VarMax(j));
end
%initial cost function for each Chromosome
Particle(i).Cost = CostFunction(Particle(i).Position,data);
end
%Loop for generations
for jj=1:MaxIt
%% Mutation of DE
%loop for all Chromosome
for it=1:nPoP
%Any random Chromosome which is not same as i
ss=rand_spl(nPoP,i);
%Choosing 3 random population
for jt=1:length(ss)
temp_Particle(jt)=Particle(ss(jt));
end
%Step for Mutation
Vector(it).Position=Particle(i).Position+K.*([temp_Particle(1).Position]-[Particle(i).Position])+...
F.*([temp_Particle(2).Position]-[temp_Particle(3).Position]);
end
%% Crossover of DE
%loop for all Chromosome
for i=1:nPoP
%loop for all genes
for j=1:nVar
%Parent Chromosome is mixing with mutated Chromosome
if rand(1)>Cross_rate && j~=randperm(nVar,1)
New_Particle(i).Position(j)=Particle(i).Position(j);
else
New_Particle(i).Position(j)=Vector(i).Position(j);
end
end
%After mixing the objective function value
New_Particle(i).Cost= CostFunction(New_Particle(i).Position,data);
end
%% Selection of DE
for i=1:nPoP
%Offspring process for next generation
if New_Particle(i).Cost<=Particle(i).Cost
New_Offspring(i)= New_Particle(i);
else
New_Offspring(i)= Particle(i);
end
end
%Selection of best offsprings for next generation
[x,idx]=sort([New_Offspring.Cost],'ascend'); %Sorting in ascending order
fmin=New_Offspring(idx(1)).Cost; %minimum value of objective function
best_var=New_Offspring(idx(1)).Position; %best parameters from optimization
Particle=New_Offspring; %all Chromosomes after optimization
%printing the best solution for each generation
fprintf('Cost Function=%.7f\n',fmin)
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% Part 3. %%%%%%%%%%%%%%%%%%%%%%%%
%%%% Plotting the results obtained form Differential Evolution %%%%
xx=x_obs; %observation points
x=best_var;%optimized parameters of gravity inversion
%bSpline method for finding depth from parameters
t=linspace(0,1,101);
%All basis functions for b-Spline
N_final=b_spline_basis(t,n,k);
%loop for finding depth profile from optimized parameters
yy_eval=0;
for i=1:n+1
%Depth profile
yy_eval=yy_eval+x(i)*N_final(:,i);
end
%Horizontal position of depth profile
xx1=linspace(min(xx),max(xx),30);
%Vertical position of depth profile
yy1=spline(xx,yy_eval,xx1);
%End point adjustment for closed form
x1=[xx1 x_obs(end) 0];
y1=[yy1 0 0];
%Gravity anomaly for optimized depth profile
zz1=poly_gravityrho(xx,0,x1,y1,@(z)(-0.08-0.42*exp(-0.522*z*(10^-3)))*1000,t_leg,c_leg);
% Plotting of observed and inverted Gravity anomaly
figure(2)
subplot(2,1,1)
hold on
%Observed Anomaly
plot(xx,data,'-o','linewidth',1.25)
%Inverted anomaly
plot(xx,zz1,'linewidth',1.25)
%Axis lebeling
xlabel('Observation points in meter')
ylabel('Gravity anomaly in mGal')
title('Gravity anomaly for San Jacinto Graben ')
legend('Observed','Inverted','Location','southeast')
box on
figure(2)
subplot(2,1,2)
hold on
%plotting of true model and optimized model
zd=(-0.08-0.42*exp(-0.522*y1*10^-3))*1000; %Depth varying density
%Optimized depth
patch(x1,y1,zd)
colormap winter
%Model depth from Chai Hinze, 1988
plot(x_obs,depth,'linewidth',1.25,'color','r')
box on
set(gca,'Ydir','reverse')
c = colorbar('location','southoutside');
c.Label.String = 'Density in kg/m^3';
xlabel('Distance in meter')
ylabel('Depth in meter')
legend('Best Model','Chai1988','Location','southeast')
title('Depth profile for San Jacinto Graben, California')
%rmse for gravity
N_g=length(data);
RMSE_g=sqrt((sum((data-zz1).^2))/N_g)/(max(data(:))-min(data(:)));
fprintf('RMSE in gravity field=%f\n',RMSE_g)
%plotting of residuals
figure(3)
subplot(2,1,1)
%plotting for gravity residual
plot(xx,data-zz1)
xlabel('Observation points in meter')
ylabel('Gravity anomaly in mGal')
title('Residual of Gravity anomaly for Chai Hinze Theoritical profile')
figure(3)
subplot(2,1,2)
%plotting for depth residual as per Chai Hinze, 1988
plot(xx,yy_eval'-depth)
set(gca,'Ydir','reverse')
xlabel('Distance in meter')
ylabel('Depth in meter')
title('Residual in depth profile of Chai Hinze Theoritical profile')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% Part 4. %%%%%%%%%%%%%%%%%%%%%%%%
%%%% Objective functions and Constraints %%%%
function val=myCostFunction(x,data,n,k,t_leg,c_leg)
%%inputs are
%x= parameters for DE algorithm
%n= number of knots
%k= order of polynomial for bspline
%100 linearly spaced data points from 0 to 1 for bspline basis
%observation points
xx=linspace(0,10.16*10^3,101);
t=linspace(0,1,101);
%all bspline basis for given knots and polynomial order
N_final=b_spline_basis(t,n,k);
%loop for depth profile reconstruction using bspline basis
yy=0;
for i=1:n+1
yy=yy+x(i)*N_final(:,i);
end
%close polygonal form of depth profile
xx1=linspace(min(xx),max(xx),20);
yy1=spline(xx,yy,xx1);
x1(1:22)=[xx1 10.16*10^3 0];
y1(1:22)=[yy1 0 0];
%gravity field for given depth profile
zz1=poly_gravityrho(xx,0,x1,y1,@(z) (-0.08-0.42*exp(-0.522*z*10^-3))*1000,t_leg,c_leg);
%misfit functional for observed and inverted gravity anomaly
val=norm(data-zz1);
end
function val=constrained1(x,n,k)
%%inputs are
%x= parameters for DE algorithm
%n= number of knots
%k= order of polynomial for bspline
%100 linearly spaced data points from 0 to 1 for bspline basis
t=linspace(0,1,101);
%all bspline basis for given knots and polynomial order
N_final=b_spline_basis(t,n,k);
%loop for depth profile reconstruction using bspline basis
yy=0;
for i=1:n+1
yy=yy+x(i)*N_final(:,i);
end
%minimum of depth
m_yy=min(yy);
%penalty barrier method for minimum bound
gg=(-m_yy+0);
val=(max(0,gg))^2;
end
function val=constrained2(x,n,k)
%%inputs are
%x= parameters for DE algorithm
%n= number of knots
%k= order of polynomial for bspline
%100 linearly spaced data points from 0 to 1 for bspline basis
t=linspace(0,1,101);
%all bspline basis for given knots and polynomial order
N_final=b_spline_basis(t,n,k);
%loop for depth profile reconstruction using bspline basis
yy=0;
for i=1:n+1
yy=yy+x(i)*N_final(:,i);
end
%maximum of depth
m_yy=max(yy);
%penalty barrier method for maximum bound
gg=(m_yy-10*1000);
val=(max(0,gg))^2;
end