Pendekatan Inversi Linier dengan Matriks Jacobi pada Kasus Permodelan Gravity Sederhana
2 likes2,382 views
Script MATLAB ini melakukan inversi gravitasi untuk menentukan posisi pusat model bawah permukaan berdasarkan data observasi anomali gravitasi. Script ini melakukan iterasi dengan metode Jacobi untuk memperkecil standar deviasi misfit antara data kalkulasi dan observasi, hingga nilai epsilon mencapai 0,00001. Hasil akhir menunjukkan data kalkulasi sangat mendekati data observasi.
![SCRIPT/SYNTAX MATLAB
clc, clear all
%DEFINISI parameter model gravitasi
kons_G=6.67*10^(-11);
R=70;
rho_model=100;
%tentukan posisi PUSAT model
x0=500;
z0=250;
%tentukan lintasan pengukuran
xlin=0:50:950;
zlin=zeros(length(xlin),1);
%menghitung DATA g_obs
g_obs=zeros(length(xlin),1);
for i=1:length(xlin)
g_obs(i)=kons_G*(4/3*pi*(R^3)*z0*rho_model)*10^5 / ...
((xlin(i)-x0)^2+(zlin(i)-z0)^2)^(3/2); %sudah dalam mGal
End
%PLOT data observasi dan juga model bawah permukaan yang digunakan
subplot(2,2,1)
plot(xlin,g_obs,'or','MarkerFaceColor','r')
xlabel('Jarak (m)'); ylabel('Gravity Anomaly (mGal)');
title('Grafik g observasi')
subplot(2,2,3)
d = R*2; px = x0-R; py = z0-R;
rectangle('Position',[px py d d],'Curvature',[1,1],'FaceColor','r');
daspect([1,1,1]); ylim([0,500]); xlim([0,1000]); set(gca,'ydir','reverse');
xlabel('Jarak (m)'); ylabel('Kedalaman (m)'); title('Model Bawah
Permukaan')
%lakukan PREDIKSI MODEL awal
iterasi=1; eps=1;
while eps >= 0.00001
if iterasi==1
x0_model=100;
z0_model=100;
else
x0_model=x0_pertu;
z0_model=z0_pertu;
end
%menghitung DATA g_cal
g_cal=zeros(length(xlin),1);
for i=1:length(xlin)
g_cal(i)=kons_G*(4/3*pi*(R^3)*z0_model*rho_model)*10^5 / ...
((xlin(i)-x0_model)^2+(zlin(i)-z0_model)^2)^(3/2);
%sudah dalam mGal
end
%menghitung MISFIT
dg_misfit=g_obs - g_cal;
eps=std(abs(dg_misfit));
e_plot(iterasi)=eps;
%membuat kondisi IF untuk inversi jacobi
if (std(abs(dg_misfit))) >= 0.00001
%membuat matriks JACOBI
for i=1:length(xlin)
derivative_x(i)=kons_G*(4/3*pi*(R^3)*rho_model)* ...
(3*z0_model*(xlin(i)-x0_model))*10^5 / ...
((xlin(i)-x0_model)^2+(zlin(i)-z0_model)^2)^(5/2);
%sudah dalam mGal
derivative_z(i)=kons_G*(4/3*pi*(R^3)*rho_model)* ...
(zlin(i)^2 + zlin(i)*z0_model + xlin(i)^2 - 2*xlin(i)*x0_model ...](https://image.slidesharecdn.com/tugasinversipendekatangravitylinierdenganjacobi-121026065233-phpapp02/85/Pendekatan-Inversi-Linier-dengan-Matriks-Jacobi-pada-Kasus-Permodelan-Gravity-Sederhana-1-320.jpg)
![2
- 2*z0_model^2 + x0_model^2)*10^5 / ...
((xlin(i)-x0_model)^2+(zlin(i)-z0_model)^2)^(5/2);
%sudah dalam mGal
end
J=ones(length(xlin),2);
J(:,1)=derivative_x';
J(:,2)=derivative_z';
%menghitung PERTURBASI MODEL
dm_perturbasi=inv(J'*J)*J'*dg_misfit;
x0_pertu=x0_model + dm_perturbasi(1);
z0_pertu=z0_model + dm_perturbasi(2);
%PLOT model prediksi bawah permukaan ter-perturbasi
subplot(2,2,4)
d = R*2; px = x0_pertu-R; py = z0_pertu-R;
rectangle('Position',[px py d d],'Curvature',[1,1],'EdgeColor','b');
daspect([1,1,1]); ylim([0,500]); xlim([0,1000]); set(gca,'ydir','reverse');
xlabel('Jarak (m)'); ylabel('Kedalaman (m)'); title('Model Prediksi')
iterasi=iterasi+1;
end
end
%PLOT data kalkulasi dan std dari misfit
subplot(2,2,2)
plot(xlin,g_cal,'ob','MarkerFaceColor','b')
hold on
plot(xlin,g_obs,'-r')
xlabel('Jarak (m)'); ylabel('Gravity Anomaly (mGal)');
title('Grafik g kalkulasi (ForMod)')
figure(2)
plot([1:1:(length(e_plot))],e_plot,'Color','m','LineStyle','-
','LineWidth',2)
xlabel('Iterasi'); ylabel('Std Misfit');
title('Grafik Misfit')
maka dari script MATLAb diatas jika di-run akan memberikan output berupa grafik sebagai berikut,](https://image.slidesharecdn.com/tugasinversipendekatangravitylinierdenganjacobi-121026065233-phpapp02/85/Pendekatan-Inversi-Linier-dengan-Matriks-Jacobi-pada-Kasus-Permodelan-Gravity-Sederhana-2-320.jpg)
More Related Content
What's hot (20)
Viewers also liked (14)
Similar to Pendekatan Inversi Linier dengan Matriks Jacobi pada Kasus Permodelan Gravity Sederhana (20)
![5_Kalkulus_Turunan_(1)[1].pptx](https://cdn.slidesharecdn.com/ss_thumbnails/5kalkulusturunan11-230926092318-7e31b4aa-thumbnail.jpg?width=560&fit=bounds)
More from Fajar Perdana (6)
Recently uploaded (7)
Pendekatan Inversi Linier dengan Matriks Jacobi pada Kasus Permodelan Gravity Sederhana
- 1. SCRIPT/SYNTAX MATLAB clc, clear all %DEFINISI parameter model gravitasi kons_G=6.67*10^(-11); R=70; rho_model=100; %tentukan posisi PUSAT model x0=500; z0=250; %tentukan lintasan pengukuran xlin=0:50:950; zlin=zeros(length(xlin),1); %menghitung DATA g_obs g_obs=zeros(length(xlin),1); for i=1:length(xlin) g_obs(i)=kons_G*(4/3*pi*(R^3)*z0*rho_model)*10^5 / ... ((xlin(i)-x0)^2+(zlin(i)-z0)^2)^(3/2); %sudah dalam mGal End %PLOT data observasi dan juga model bawah permukaan yang digunakan subplot(2,2,1) plot(xlin,g_obs,'or','MarkerFaceColor','r') xlabel('Jarak (m)'); ylabel('Gravity Anomaly (mGal)'); title('Grafik g observasi') subplot(2,2,3) d = R*2; px = x0-R; py = z0-R; rectangle('Position',[px py d d],'Curvature',[1,1],'FaceColor','r'); daspect([1,1,1]); ylim([0,500]); xlim([0,1000]); set(gca,'ydir','reverse'); xlabel('Jarak (m)'); ylabel('Kedalaman (m)'); title('Model Bawah Permukaan') %lakukan PREDIKSI MODEL awal iterasi=1; eps=1; while eps >= 0.00001 if iterasi==1 x0_model=100; z0_model=100; else x0_model=x0_pertu; z0_model=z0_pertu; end %menghitung DATA g_cal g_cal=zeros(length(xlin),1); for i=1:length(xlin) g_cal(i)=kons_G*(4/3*pi*(R^3)*z0_model*rho_model)*10^5 / ... ((xlin(i)-x0_model)^2+(zlin(i)-z0_model)^2)^(3/2); %sudah dalam mGal end %menghitung MISFIT dg_misfit=g_obs - g_cal; eps=std(abs(dg_misfit)); e_plot(iterasi)=eps; %membuat kondisi IF untuk inversi jacobi if (std(abs(dg_misfit))) >= 0.00001 %membuat matriks JACOBI for i=1:length(xlin) derivative_x(i)=kons_G*(4/3*pi*(R^3)*rho_model)* ... (3*z0_model*(xlin(i)-x0_model))*10^5 / ... ((xlin(i)-x0_model)^2+(zlin(i)-z0_model)^2)^(5/2); %sudah dalam mGal derivative_z(i)=kons_G*(4/3*pi*(R^3)*rho_model)* ... (zlin(i)^2 + zlin(i)*z0_model + xlin(i)^2 - 2*xlin(i)*x0_model ...
- 2. 2 - 2*z0_model^2 + x0_model^2)*10^5 / ... ((xlin(i)-x0_model)^2+(zlin(i)-z0_model)^2)^(5/2); %sudah dalam mGal end J=ones(length(xlin),2); J(:,1)=derivative_x'; J(:,2)=derivative_z'; %menghitung PERTURBASI MODEL dm_perturbasi=inv(J'*J)*J'*dg_misfit; x0_pertu=x0_model + dm_perturbasi(1); z0_pertu=z0_model + dm_perturbasi(2); %PLOT model prediksi bawah permukaan ter-perturbasi subplot(2,2,4) d = R*2; px = x0_pertu-R; py = z0_pertu-R; rectangle('Position',[px py d d],'Curvature',[1,1],'EdgeColor','b'); daspect([1,1,1]); ylim([0,500]); xlim([0,1000]); set(gca,'ydir','reverse'); xlabel('Jarak (m)'); ylabel('Kedalaman (m)'); title('Model Prediksi') iterasi=iterasi+1; end end %PLOT data kalkulasi dan std dari misfit subplot(2,2,2) plot(xlin,g_cal,'ob','MarkerFaceColor','b') hold on plot(xlin,g_obs,'-r') xlabel('Jarak (m)'); ylabel('Gravity Anomaly (mGal)'); title('Grafik g kalkulasi (ForMod)') figure(2) plot([1:1:(length(e_plot))],e_plot,'Color','m','LineStyle','- ','LineWidth',2) xlabel('Iterasi'); ylabel('Std Misfit'); title('Grafik Misfit') maka dari script MATLAb diatas jika di-run akan memberikan output berupa grafik sebagai berikut,
- 3. 3 Terlihat dari hasil diatas jika anomali gravitasi hasil kalkulasi (titik biru) sangat mendekati data observasi model bawah permukaan kita (garis merah). Script diatas memakai nilai batas epsilon sebesar 0.00001, sehingga nilai standar deviasi misfit yang didapat sangat baik. Sehingga jika kita menurunkan nilai batas epsilon menjadi 0.0025 maka kita dapat melihat antara data kalkulasi forward modeling dan data observasi belum fit. Berikut gambarannya,