-
Notifications
You must be signed in to change notification settings - Fork 0
/
Fourier_PSO.m
251 lines (226 loc) · 9.74 KB
/
Fourier_PSO.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
% ***************************************************************
% *** Matlab function for inversion of sedimentary basin from observed
% gravity anomalies.
% *** 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
% ****************************************************************
%%Matlab function for inversion of sedimentary basin using PSO and adaptive
%%Fourier coefficients
function [depth_inverted,grav_inverted]=Fourier_PSO(x_obs_data,gravity_anomaly_data,density,max_depth)
%input:
% 1. x_obs = data file for horizontal location of
% observation points (single collumn data file in ascii format in m.)
% 2. gravity_anomaly = corresponding gravity anomaly data file at the
% observation points (single collumn data file in ascii format in mGal)
% 3. density = density contrast as a function of z in
% kg/m^3
% 4. max_depth = Maximum depth bound for inversion
%output:
% 1.depth_inverted = inverted depth profile of sedimentary basin
% at the locations of horizontal observation
% points in m.
% 1.grav_inverted = inverted gravity anomaly of sedimentary basin
% at the locations of horizontal observation
% points in mGal.
%% Gravity anomaly and its inversion for any sedimentary basin
%import data for observation points and gravity anomaly
x_obs=(importdata(x_obs_data))';
zz1=(importdata(gravity_anomaly_data))';
sz1=size(x_obs);
if sz1(1,1)~=1
x_obs=x_obs';
end
sz2=size(zz1);
if sz2(1,1)~=1
zz1=zz1';
end
%t and c are Legendre Gaussian quadrature points for numerical integration
[t_leg,c_leg]=lgwt(10,0,1);
z_obs=0; %height of observation point is in surface
%Fourier coefficients for Gravity anomaly
TT=1:2*length(zz1);
ll=(TT(end)-TT(1))/2;
%Gravity data for Fourier series
grav_data=[zz1 fliplr(zz1)];
%1st 100 Fourier Coefficients
for j=1:100
%coefficients for Gravity field
ss1=grav_data.*cos(j*pi*TT/ll);
aa_grv(j)=(1/ll)*trapz(TT,ss1);
ss2=grav_data.*sin(j*pi*TT/ll);
bb_grv(j)=(1/ll)*trapz(TT,ss2);
vv_grv(j+1)=sqrt((aa_grv(j))^2+(bb_grv(j))^2);
wn_num(j+1)=j*pi/ll;
end
%normalization of Fourier Coefficients of Gravity anomaly
vv_grv(1)=abs((1/ll)*trapz(TT,grav_data));
vv_grv=vv_grv./max(vv_grv);
wn_num(1)=0;
figure(1)
%Plotting the Fourier Spectrum of Gravity anomaly
semilogy(wn_num(1:51),vv_grv(1:51))
hold on
semilogy(wn_num(1:51),10^-2.*ones(51,1),'k--')
title('Fourier power spectrum for gravity anomaly')
xlabel('Wavenumbers')
ylabel('Amplitude')
%selecting number of Fourier Coefficients
n=find(vv_grv>=10^-2, 1, 'last' )-1;
%printing the number of parameters
fprintf('Number of parameters for Fourier coefficients of Depth profile = %d\n',n)
data_x=x_obs; %data for horizontal locations
data=zz1; %data for gravity anomaly
%% creating Fourier Transformation matrix for multiplication
T=1:2*length(data_x);
l=(T(end)-T(1))/2;
for i=1:length(T)
jj1=0;jj2=0;
for j=1:2*n+1
if j==1
A_mat(i,j)=(1/2);
elseif j>1 && j<=n+1
jj1=jj1+1;
A_mat(i,j)=cos(jj1*pi*T(i)/l);
else
jj2=jj2+1;
A_mat(i,j)=sin(jj2*pi*T(i)/l);
end
end
end
%% Problem Definition
%%%% Optimization of Sayula using PSO %%%%
c1=1.4; c2=1.7;
%Cost function with constraints for optimization
CostFunction =@(x) myCostFunction(x,x_obs,data,A_mat,t_leg,c_leg,density)+1000*(constrained1(x,A_mat)+constrained2(x,A_mat,max_depth));
nVar=(2*n+1); %Number of Unknown Variables
MaxIt = 300; %Maximum number of iterations
nPoP = 40; %Population size or swarm size
%Loop for 2 independent run
for ii=1:2
%Calling the function for optimization using PSO
[bst_var, best_cost,iter_count,error_energy] = WIPSO(CostFunction,nVar,MaxIt,nPoP,c1,c2);
%Saving all best cost and parameters for each independent run
all_best(ii)=best_cost;
all_var(:,ii)=bst_var;
err_enrgy(:,ii)=error_energy';
end
%finding best model
[val,id]=min(all_best);
%Parameters for best Model
best_var=(squeeze(all_var(:,id)))';
%error energy for best model
best_error=(squeeze(err_enrgy(:,id)))';
%% Plotting
%%%% Plotting the results for Optimization model %%%%
x=best_var'; % Fourier coefficients for optimized model
xx=x_obs; % observation points
%gravity field at 40 discrete data points
xx1=linspace(min(xx),max(xx),40);
%depth profile from optimized Fourier coefficients
yy_eval=(A_mat*x)';
yy=yy_eval(1:length(yy_eval)/2);
yy1=spline(xx,yy,xx1);
%close polygonal form of depth profile
x1(1:length(xx1)+2)=[xx1 max(xx1) min(xx1)];
y1(1:length(yy1)+2)=[yy1 0 0];
%gravity field for optimized depth profile
zz_eval=poly_gravityrho(xx,0,x1,y1,density,t_leg,c_leg);
depth_inverted=yy;
grav_inverted=zz_eval;
%true model
figure(2)
subplot(2,1,2)
hold on
xx_p=[xx1(1) xx1 xx1(end)];
depth_p=[0 yy1 0];
%patch plot for true model
zd_p=density(depth_p);
patch(xx_p,depth_p,zd_p)
colormap copper
%optimized model
plot(xx1,yy1,'linewidth',1.25,'color','r')
set(gca,'Ydir','reverse')
c = colorbar('location','southoutside');
c.Label.String = 'Density contrast (kg/m^3)';
xlabel('Distance (m)')
ylabel('Depth (m)')
title('Depth profile')
box on
subplot(2,1,1)
hold on
%Observed
plot(xx,data,'-o','linewidth',1.25)
%Inverted
plot(xx,zz_eval,'linewidth',1.25,'color','r')
%Axis lebeling
xlabel('Observation points(m)')
ylabel('Gravity anomaly (mGal)')
title('Gravity anomaly ')
legend('Observed Data','Inverted Data','location','southeast')
box on
%RMSE for gravity
N_g=length(data); %Number of Observation points
RMSE_g=sqrt((sum((data-zz_eval).^2))/N_g)/(max(data(:))-min(data(:)));
%Printing the RMSE error for depth and gravity profile
fprintf('RMSE in gravity field=%f\n',RMSE_g)
%%%% Objective functions and Constraints %%%%
function val=myCostFunction(x,x_obs,data,A_mat,t_leg,c_leg,density)
%%inputs are
%x= parameters for PSO algorithm
%data= observed gravity field
%A_mat= Matrix with Fourier coefficients
%t_leg and c_leg are Legendre Gaussian quadrature points for numerical integration
% subroutine for t_leg and c_leg evaluation is given in lgwt.m file
%observation points having 100 linearly spaced data points
xx_v=x_obs;
%40 linearly spaced data points in closed polygonal form
xx1_v=linspace(min(xx_v),max(xx_v),40);
yy1_v=(A_mat*x')';
yy_v=yy1_v(1:length(yy1_v)/2);
yy1_v=spline(xx_v,yy_v,xx1_v);
%close polygonal form of depth profile
x1(1:length(xx1_v)+2)=[xx1_v xx1_v(end) xx1_v(1)];
y1(1:length(yy1_v)+2)=[yy1_v 0 0];
%gravity field for given depth profile
zz1=poly_gravityrho(xx_v,0,x1,y1,density,t_leg,c_leg);
%zz1=poly_gravityrho(xx,0,x1,y1,density,t_leg,c_leg);
%misfit functional for observed and inverted gravity anomaly
val=norm(data-zz1);
end
%Constraints for upper bound
function val=constrained1(x,A_mat)
%%inputs are
%x= parameters for DE algorithm
%A_mat= Matrix with Fourier coefficients
%depth profile from Fourier domain to time domain
yy1=A_mat*x';
yy=yy1(1:length(yy1)/2);
%minimum of depth
m_yy=min(yy);
%penalty barrier method for minimum bound
gg=(-m_yy+0);
val=(max(0,gg))^2;
end
%Constraints for lower bound
function val=constrained2(x,A_mat,max_depth)
%%inputs are
%x= parameters for DE algorithm
%A_mat= Matrix with Fourier coefficients
%depth profile from Fourier domain to time domain
yy1=A_mat*x';
yy=yy1(1:length(yy1)/2);
%maximum of depth
m_yy=max(yy);
%penalty barrier method for minimum bound
gg=(m_yy-max_depth);
val=(max(0,gg))^2;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%