Script Sederhana Penentuan Episenter Gempa Bumi

I'm sure we all have perfectly known what the earthquake is. Yes... when we feel the ground shakes that makes us feel swaying hither and thither, that's the time of earthquake. Recently, our beloved country has suffered many earthquakes, in Aceh, Yogyakarta, Pangandaran, and even Jakarta the capital city wasn't escaped from the earthquake. Some of these earthquake were followed by tsunami, the others were not. Why this could happen?


Beberapa minggu yang lalu saya ngeliat file-file dan arsip-arsip tugas yang pernah saya kerjakan semenjak awal kuliah. Dulu, banyak tentunya yang sudah saya kerjakan. Sekian lama hingga kini rasanya saya masih interest bangat dengan kajian-kajian studi tersebut. Teringat juga saat-saat ngerjain praktikum dan bikin laporannya 3 sampai 4 kali seminggu. Mulai dari practice, data processing, analyse, lalu buat resume alias reportnya. Wahh, kebayang gimana lelahnya. Blum lagi presentasi depan para senior tengah malem hingga pagi tuhh. Hihihii…waktu istirahatnya dikit bangat paling waktu kuliah (tidur..huusss). Kudu pas weekend doang bisa spend time to rest and relax.

So far, saya pikir-pikir apa sihh feedback nya, tugas-tugas yang sudah dikerjain selain dapetin nilai, mau dikemanain? Kalo di share n dipublish bagus nggak sihh? Lumayan kan jadi contoh-contoh journal gratisan…hehehe. Publish di blog, wordpress, atau apapunlah. Nahh, salah satunya yang pengen saya publish kali ini adalah kajian tentang penentuan episenter gempa bumi. Ini dasar bangat, mungkin masih banyak kajian yang lebih advanced lagi. Di sini, untuk menentukan posisi episenternya dengan memakai Metode Pendekatan Linier Bentuk Forward Modeling dan Pendekatan Global (Metoda Grid Search dan Ramdom Search). Kedua metode tersebut merupakan solusi dari Inversi Non Linier.

Kita tahu bahwa pusat gempa tektonik terdapat pada tempat-tempat terjadinya akumulasi energi. Gesekan yang terjadi pada batuan dapat mengakibatkan penumpukan energi di daerah tersebut sehingga pada akhirnya akan pecah dan terjadilah gempa. Akumulasi energi ini terjadi di daerah tiga batas lempeng yaitu di daerah konvergen (pertemuan dari dua buah lempeng sehingga salah satunya menunjam di bawah lempeng lainnya), daerah divergen (dua buah lempeng yang saling menjauh) , dan daerah sesar mendatar (dua lempeng tektonik yang saling berpapasan mengakibatkan terjadinya sesar secara mendatar, berbentuk memanjang).

Titik dalam perut bumi yang merupakan sumber gempa dinamakan hiposenter atau fokus. Proyeksi tegak lurus hiposenter ini ke permukaan bumi dinamakan episenter. Gelombang gempa merambat dari hiposenter ke patahan sesar fault rupture. Bila kedalaman fokus dari permukaan adalah 0 - 70 km, terjadilah gempa dangkal (shallow earthquake), sedangkan bila kedalamannya antara 70 - 700 km, terjadilah gempa dalam (deep earthquake). Gempa dangkal menimbulkan efek goncangan yang lebih dahsyat dibanding gempa dalam. Ini karena letak fokus lebih dekat ke permukaan, dimana batu-batuan bersifat lebih keras sehingga melepaskan lebih besar regangan (strain).

Berikut ini adalah contoh percobaan dengan menggunakan data sintetik untuk menentukan posisi episenter gempa bumi. Data sintetik ini diolah dengan kedua metode di atas. Data processing nya dengan memakai software programming Matlab, yang hasilnya diperoleh solusi inversi non linier yang akan diplot dalam bentuk kurva dan kontur fungsi obyektif hasil interpolasi yang tersebar secara random. Singkatnya, di bawah ini saya sertakan langkah-langkah dan script programming nya.

Langkah-Langkah Pengolahan Data dengan menggunakan Matlab.
Buat data sintetik waktu tempuh gelombang gempa dari episenter (Xo,Yo) atau hiposenter (Xo,Yo,Zo), yang Anda tentukan ke-4 atau 5 stasiun (XSi ,YSi) i= 1,2, …,4 atau (XSi ,YSi , ZSi) i= 1,2, …,4, pilih melingkupi sekitar episenter.
• Lakukan inverse
• Plot solusi tiap iterasi
• Plot ERMS


Metode 1 : Solusi Inversi Non-Linier Dengan Pendekatan Linier Bentuk Forward Modeling

Script Program M-File:
-------------------------------------------------------------------------------------
clear all
clc
to=0; % origin time
vp=4; % kecepatan gelombang gempa
ti=[7.1;1.8;5;7.9]; % waktu tempuh masing2 stasiun
x=[20;50;40;10]; % posisi longitude
y=[10;25;50;40]; % posisi latitude
n=length(x);

xo=40; yo=30;
M=[40;30]; % posisi tebakan awal(estimasi model awal)
[X,Y]=meshgrid(0:10:80,0:10:80);

M1=[];
for n=1:5
gm=to+(1/vp)*(sqrt((x-M(1)).^2+(y-M(2)).^2));
dgm_dm1=(1/vp)*(-(x-M(1)))./(sqrt((x-M(1)).^2+(y-M(2)).^2));
dgm_dm2=(1/vp)*(-(x-M(2)))./(sqrt((x-M(1)).^2+(y-M(2)).^2));
J=[dgm_dm1 dgm_dm2];
Mo=M;
M=Mo+inv(J'*J)*J'*(ti-gm)
M1=[M1 M];
Mn=M1';
end;
plot(Mn(:,1),Mn(:,2),'o')
hold on
plot(Mn(:,1),Mn(:,2))
hold on
plot(x,y,'ko')
hold on
plot(xo,yo,'ks')
hold on
plot(X,Y,'.')
-------------------------------------------------------------------------------------

Hasil Plotting:
Plotting Kurva Penentuan Posisi Episenter Gempa Hasil Solusi Inversi Non-Linier Melalui Pendekatan Linier






Metode 2 : Solusi Inversi Non-Linier Dengan Pendekatan Global (Metoda Grid Search dan Ramdom Search)

Script Program M-File:
-------------------------------------------------------------------------------------
clear all
clc
to=0; % origin time
vp=4; % kecepatan gelombang gempa
ti=[7.1;1.8;5;7.9]; % waktu tempuh masing2 stasiun
x=[20;50;40;10];
y=[10;25;50;40]; % posisi stasiun
h=length(x);
plot(x,y,'h')
hold on

M=[40;30]; % posisi tebakan awal
[X,Y]=meshgrid(0:10:80,0:10:80);
plot(X,Y,'.')

M1=[];
for n=1:5
gm=to+(1/vp)*(sqrt((x-M(1)).^2+(y-M(2)).^2));
dgm_dx=(1/vp)*(-(x-M(1)))./(sqrt((x-M(1)).^2+(y-M(2)).^2));
dgm_dy=(1/vp)*(-(x-M(2)))./(sqrt((x-M(1)).^2+(y-M(2)).^2));
J=[dgm_dx dgm_dy];
Mo=M;
M=Mo+inv(J'*J)*J'*(ti-gm)
M1=[M1 M];
Mn=M1';
end;
plot(Mn(:,1),Mn(:,2),'o')
hold on
plot(Mn(:,1),Mn(:,2))

% oleh stasiun pertama
t_cal1=to+(1/vp)*sqrt((x(1)-X).^2+(y(1)-Y).^2);
Error1=(t_cal1-ti(1)).^2;

% oleh stasiun kedua
t_cal2=to+(1/vp)*sqrt((x(2)-X).^2+(y(2)-Y).^2);
Error2=(t_cal2-ti(2)).^2;

% oleh stasiun ketiga
t_cal3=to+(1/vp)*sqrt((x(3)-X).^2+(y(3)-Y).^2);
Error3=(t_cal3-ti(3)).^2;

% oleh stasiun keempat
t_cal4=to+(1/vp)*sqrt((x(4)-X).^2+(y(4)-Y).^2);
Error4=(t_cal4-ti(4)).^2;

Erms=sqrt((1/n)*(Error1+Error2+Error3+Error4));
[cs,h]=contour(X,Y,Erms,[0:0.5:15])
clabel(cs,h)
-------------------------------------------------------------------------------------


Hasil Plotting :
Kontur Fungsi Obyektif Hasil Interpolasi Yang Tersebar Secara Acak Beserta Posisi Episenter Gempa Yang Diplot Berdasarkan Estimasi




Try it out.

Posted in: