SIMULASI PERHITUNGAN WAKTU TEMPUH GELOMBANG DENGAN METODA EIKONAL : SUATU CONTOH APLIKASI DALAM ESTIMASI KETELITIAN HIPOSENTER GEMPA

Simulasi Perhitungan Waktu Tempuh Gelombang dengan Metoda Eikonal : Suatu Contoh Aplikasi Dalam Estimasi Ketelitian Hiposenter Gempa (Yasa Suparman, Hendra Gunawan)

SIMULASI PERHITUNGAN WAKTU TEMPUH GELOMBANG DENGAN METODA EIKONAL : SUATU CONTOH APLIKASI DALAM ESTIMASI KETELITIAN HIPOSENTER GEMPA Yasa SUPARMAN dkk Pusat Vulkanologi dan Mitigasi Bencana Geologi Badan Geologi

Sari Persamaan eikonal merupakan bagian dari persamaan gelombang yang dapat dijadikan dasar perhitungan waktu tempuh gelombang pada medium akustik. Secara eksplisit persamaan ini menyatakan hubungan antara waktu tempuh dengan kecepatan lokal suatu medium yang dinyatakan dalam suatu persamaan diferensial parsial nonlinier orde pertama. Tulisan ini merupakan salah satu aplikasi dari suatu rangkaian pekerjaan dari beberapa pekerjaan sebelumnya. Salah satu aplikasi yang bisa dikembangkan dari simulasi ini untuk waktu ke depan adalah untuk penentuan lokasi gempa dekat dasar permukaan kawah aktif seperti gempa ”long-period” atau ”low frequency”. Pada tulisan ini hanya akan membahas hasil simulasi waktu tiba data gempa sintetik. Untuk membuat validasi hasil simulasi di atas metoda penentuan lokasi gempa ”Geiger Adaptive Damping” (GAD) akan digunakan sebagai pembanding. Kata kunci : persamaan eikonal, GAD

1. Pendahuluan Penentuan lokasi sumber gempa (hiposenter) merupakan tahap awal dalam melakukan analisa kegempaan. Berbagai software dapat digunakan untuk menentukan posisi hiposenter, misalnya: GAD, GrHypo, Hypo71, Hypoellipse dan sebagainya. Geiger’s Adaptive Damping (GAD), Nishi (2001), merupakan salah satu software yang umum digunakan untuk penentuan posisi hiposenter terutama dalam penentuan lokasi hiposenter gempa di daerah gunungapi atau pada daerah yang mempunyai jarak yang relatif dekat antara sumber gempa dan penerima (receiver). Data yang harus dipersiapkan untuk menjalankan program ini adalah data waktu tiba, posisi seismometer dan struktur kecepatan. Beberapa model sintetis dikembangkan untuk verifikasi program GAD sehingga dapat diketahui pengaruh input data, terutama data waktu tiba, terhadap output yang berupa posisi sumber gempa. Permasalahan dalam penentuan waktu tiba pada model terletak dalam penentuan waktu tempuh gelombang. Penentuan waktu tempuh gelombang ditentukan dari lintasan perambatan gelombang (raypath) dari sumber ke penerima. Lintasan perambatan gelombang (raypath) di GAD pada dasarnya menggunakan metode shooting. Pada metode shooting permasalahan dirumuskan dengan

mencari sudut atau arah tembak yang tepat sehingga lintasan gelombang berujung tepat pada koordinat stasion penerima. Pada tulisan ini penentuan waktu tempuh diperoleh dengan terlebih dahulu melakukan penentuan muka gelombang dan lintasan perambatan (raypath) gelombang. Muka gelombang yang dibangun berdasarkan penyelesaian persamaan eikonal dengan metoda beda hingga sedangkan lintasan perambatan gelombang ditentukan berdasarkan prinsip resiprositas (Vidale, 1990; Qin et al, 1992; Matsuoka dan Ezaka ,1992; Andri, 2006). Pemodelan penentuan muka gelombang dan lintasan perambatan (raypath) gelombang hanya dilakukan pada model lapisan medium yang sederhana sehingga penentuan waktu tempuh gelombang dengan menggunakan metoda shooting ataupun metoda eikonal akan menghasilkan nilai yang hampir sama. 2. Metoda Penentuan muka gelombang dibangun berdasarkan persamaan eikonal. Persamaan gelombang untuk medium akustik yang menggunakan asumsi frekuensi tinggi akan menghasilkan dua persamaan parsial nonlinier orde pertama, yaitu persamaan eikonal dan persamaan transport (Andri, 2006). Persamaan

Bulletin Vulkanologi dan Bencana Geologi, Volume 5 Nomor 1, Januari 2010: 1-6

Hal :1


gelombang akustik untuk medium dengan densitas homogen dinyatakan dengan:

1 ∂2 p ∇ p= 2 2 c ∂t 2

................................ (1)

p dan c masing-masing menyatakan tekanan dan kecepatan lokal suatu media. Jika diasumsikan persamaan (1) memiliki bentuk solusi:

p ( xi , t ) = P( xi ) exp[− iω (t − T ( xi ))] ...... (2) Persamaan (2) menggunakan asumsi nilai frekuensi ω tinggi, ω >> 0. Fungsi P(xi) dan T(xi) diasumsikan juga merupakan fungsi kontinyu. Dengan menggunakan identitas

∇.ab = b.∇a + a∇.b , bagian persamaan (1) menjadi:

kanan

pada

∇ 2 p = ∇.∇p = (iω (∇P + iωP∇T ).∇T ) + (∇ 2 P + iω∇ T .∇ P + iωP∇ 2 P ) exp [− iω (t − T ( x i )) ] Dengan menggunakan persamaan (3), maka persamaan (1) menjadi:

1⎤ ⎡ − ω 2 P ⎢(∇T ) 2 − 2 ⎥ + c ⎦ ⎣ iω 2∇P.∇T + P∇ 2T + ∇ 2 P = 0

[

]

..... (4)

Persamaan (4) berlaku juga untuk nilai frekuensi ω yang bukan nol, sehingga persamaan (4) mempunyai solusi yang memenuhi:

(∇T ) 2 =

1 c2

...............

(5)

Persamaan (5) disebut persamaan eikonal dan persamaan (6) disebut persamaan transport. Dalam prakteknya simulasi perhitungan waktu tempuh yang digunakan akan mengacu pada program yang sudah dikembangkan oleh beberapa penulis sebelumnya (Andri, 2006 dan Nishi, 2001). 3. Solusi dan Hasil Salah satu teknik untuk menyelesaikan persamaan eikonal secara numerik adalah dengan menggunakan metoda beda hingga (Vidale, 1988;1990). Metoda beda hingga dilakukan dengan terlebih dahulu membuat diskritisasi model menjadi elemen - elemen kecil (grid) dengan ukuran dan jumlah tertentu. Tiap grid memiliki empat titik sudut. Pada tiap tiap sudut, nilai fungsi (unknown) yang akan dicari nilainya dapat dijumpai. Vidale(1988) dan Qin et al. (1992) mengajukan skema di bawah ini untuk menyelesaikan persamaan eikonal. 1) Asumsikan titik A merupakan sumber rambatan gelombang (TA = 0). Waktu (3).... rambatan ................. gelombang dari titik A ke masing-masing titik B1, B2, B3 dan B4 (Gambar 1) dihitung dengan melakukan operasi perkalian antara jarak h dengan pelambatan rata - rata (s) antara kedua titik A dan Bi menggunakan persamaan,

⎛ sB + s A ⎞ ⎟⎟ TBi = h⎜⎜ i 2 ⎝ ⎠ Kemudian hitung pula waktu rambatan gelombang pada keempat titik C1, C2, C3 dan C4 dengan menggunakan persamaan,

Tci = T A + 2(h s i ) 2 − (TBi +1 − TBi ) 2 Ketika i = 4 maka

si =

dan

1 ( s A + sCi + s Bi + s Bi +1 ) 4

2∇P.∇T + P∇ 2T = 0 ......................... (6)

Hal :2

TBi +1 = TBi


dan


Gambar 1. Titik - titik nodal sekeliling titik sumber A. Sumber: Qin et al.(1992)

2) Mencari titik nodal dengan waktu rambatan minimum. Waktu rambatan yang didapat merupakan waktu rambatan minimum global karena data waktu dari perhitungan delapan titik sebelumnya disimpan pada suatu larik parameter. Dari titik nodal dengan waktu rambatan minimal ini muka gelombang berpropagasi hingga ke seluruh domain. Tiap kali berpropagasi elemen elemen larik parameter terus diperbarui dengan data yang baru.

V1

Gambar 2 memperlihatkan hasil penentuan muka gelombang dan lintasan perambatan gelombang (raypath) dari sumber (tanda bintang) ke penerima (koordinat (0,0)). Pada model tersebut digunakan medium dengan kecepatan V1 = 5800 m/s; V2 = 6700 m/s dan V3 = 7800 m/s. Perbedaan interval kontur muka gelombang (secara spasial) di tiap lapisan dapat terlihat jelas. Kontur muka gelombang yang lebih lebar berasosiasi dengan lapisan yang memiliki kecepatan relatif lebih besar, sedangkan spasi spasial yang relatif lebih rapat berasosiasi dengan lapisan yang memiliki kecepatan relatif lebih lambat. Lintasan perambatan gelombang (raypath) selain mengalami pembelokan pada bidang batas akibat perbedaan nilai kecepatan media dipengaruhi juga oleh nilai waktu rambatan minimal pada tiap nodal.

V1 V2

V2

V3

Gambar 2 Hasil penentuan muka gelombang dan lintasan perambatan gelombang (raypath) dari sumber (tanda bintang) ke penerima (koordinat (0,0)).


Hal :3


Jaringan seismometer pada model akan ditempatkan pada posisi yang sama (Tabel 1), koordinat (0,0) adalah puncak/kawah, perbedaan antara model hanya terletak pada posisi sumber (Gambar 3). Origin time pada model, yaitu 0 detik sehingga input data waktu tiba akan sama dengan waktu tempuh gelombang dari sumber ke seismometer/penerima. Input data waktu tiba pada masing-masing model diperlakukan beda, yaitu (1) hanya data waktu tiba gelombang P, (2) data waktu tiba gelombang P (Tp) dan gelombang S (Ts).

4. Verifikasi Geiger’s Adaptive Damping (GAD) Model sintetis dikembangkan untuk verifikasi program GAD (Nishi, 2001) sehingga dapat diketahui pengaruh input data terhadap output yang berupa posisi sumber gempa. Input data untuk menjalankan program ini adalah posisi seismometer, waktu tiba dan struktur kecepatan. Posisi sumber yang terletak di tengah dan terkepung jaringan seismometer/ penerima merupakan posisi ideal dalam penentuan koordinat hiposenter (Andri, 2006). Jaringan seismometer yang tersebar pada kecenderungan arah tertentu akan akurat pada arah tersebut pula. Misalkan apabila jaringan seismometer tersebar pada arah longitude maka hasilnya (koordinat hiposenter) akan akurat pada komponen longitudenya saja, begitu juga sebaliknya. Oleh karenanya maka jaringan seismometer pada model akan ditempatkan pada posisi yang seimbang baik arah latitude (y) maupun arah longitude (x).

Tabel 1. Posisi Seismometer X (km)

A B C D

2,2 0 -3 0

Y (km) 0 -3,4 0 2,6

Z (km) -1,8 -1,0 -1,2 -1,6

5

5

(a)

4

(b)

4

3

3

2

2

1

1

Y (km)

Y (km)

Nama Station

0

0

-1

-1

-2

-2

-3

-3

-4

-4 -5

-5 -5

-4

-3

-2

-1

0

X (km)

1

2

3

4

5

-5

-4

-3

-2

-1

0

1

2

3

4

5

X (km)

Gambar 3. Posisi sumber (bintang) dan jaringan seismometer (segitiga). (a) Posisi sumber pada koordinat (0;0;5). (b) Posisi sumber pada koordinat (0;4;5)

Hal :4



Pada model hanya diterapkan struktur kecepatan sederhana, yaitu homogen dan multilayer 2 lapisan. Perbandingan kecepatan gelombang P dan S (Vp/Vs) = 1,73. Hasil inversi GAD yang diterima adalah hasil jumlah perbedaan waktu tiba observasi dengan perhitungan (∑(Tobs – Tcalc)) adalah 0.00.

•

Medium satu lapisan (homogen) Medium lapisan menggunakan kecepatan homogen, yaitu Vp = 2,7 km/s. Kemudian selanjutnya dilakukan inversi dengan menggunakan GAD untuk mendapatkan lokasi hiposenter dengan menggunakan data waktu rambat empat stasion. Lokasi hiposenter kalkulasi (hasil inversi) ditampilkan pada Tabel 2.

Tabel 2. Hasil Inversi GAD untuk Media Homogen

Data Waktu Tiba 4 Tp 4 Tp dan 1 Ts Data Waktu Tiba 4 Tp 4 Tp dan 1 Ts

MODEL A (sumber di posisi (0;0;5)) Posisi Hiposenter Error Hiposenter di GAD X Y Z X Y Z time -0,011 -0,026 4,935 0,000 0,000 0,000 0,000 -0,008 -0,022 4,992 0,002 0,002 0,003 0,001 MODEL B (sumber di posisi (0;4;5)) Posisi Hiposenter Error Hiposenter di GAD X Y Z X Y Z time 0,316 3,470 4,668 0,000 0,000 0,000 0,000 0,388 3,582 4,868 0,003 0,003 0,002 0,001

Tabel 2 memperlihatkan bahwa posisi sumber yang berada di tengah dan terkepung jaringan seismometer menghasilkan koordinat hiposenter yang jauh lebih baik, yakni posisi hiposenter hasil perhitungan lebih mendekati terhadap posisi hiposenter model. Hasil inversi Model A menghasilkan koordinat sumber yang mempunyai perbedaan sangat kecil dengan koordinat sumber model (error < 1%). Hasil inversi Model B memperlihatkan bahwa koordinat sumber mempunyai perbedaan yang lebih besar ( error > 3%). Input data waktu tiba dengan hanya data gelombang P menghasilkan error hasil inversi selalu 0,0 padahal koordinat yang dihasilkan mempunyai kesalahan/ tidak tepat. Data waktu tiba dengan input data gelombang P dan S

mempunyai error posisi lebih baik dibandingkan dengan tanpa menggunakan gelombang S. Input data waktu tiba pada media homogen dengan menggunakan hanya data gelombang P ataupun dengan data gelombang P dan S menghasilkan posisi hiposenter yang tidak terlalu berbeda. •

Medium multilayer dua lapisan Medium lapisan menggunakan kecepatan gelombang P (Vp) 2,7 km/s untuk lapisan atas dan 3.3 km/s untuk lapisan bawah. Batas lapisan atas dan bawah diletakkan pada kedalaman 0 meter (sama dengan permukaan laut rata-rata). Lokasi hiposenter kalkulasi (hasil inversi) ditampilkan pada Tabel 3.


Hal :5


Tabel 3. Hasil Inversi GAD untuk Media Multilayer Dua Lapisan

Data Waktu Tiba 4 Tp 4 Tp dan 1 Ts Data Waktu Tiba 4 Tp 4 Tp dan 1 Ts

MODEL A (sumber di posisi (0;0;5)) Posisi Hiposenter Error Hiposenter di GAD X Y Z X Y Z time -0,124 -0,255 3,325 0,000 0,000 0,000 0,000 0,004 -0,148 5,406 0,110 0,098 0,156 0,049 MODEL B (sumber di posisi (0;4;5)) Posisi Hiposenter Error Hiposenter di GAD X Y Z X Y Z time -0,251 1,514 0,701 0,000 0,000 0,000 0,000 0,127 3,976 5,027 0,079 0,101 0,075 0,033

Tabel 3 memperlihatkan bahwa, pada media multilayer dua lapisan, posisi hiposenter hasil inversi dengan input data waktu tiba gelombang P dan S menghasilkan posisi yang lebih baik dibandingkan hanya gelombang P. Posisi hiposenter hasil inversi pada Model A dengan input hanya gelombang P mempunyai error episenter kecil sedangkan pada komponen kedalaman (Z) errornya > 30%. Pada Model B, posisi hiposenter hasil inversi dengan input hanya gelombang P mempunyai error episenter, terutama pada komponen Y, lebih dari 20% dan pada komponen kedalaman (Z) errornya > 50%. Error episenter pada model B disebabkan karena posisi sumber tidak terletak di tengah dan terkepung jaringan seismometer. 5. Kesimpulan • Posisi sumber yang berada di tengah dan terkepung jaringan seismometer menghasilkan koordinat hiposenter yang jauh lebih baik. • Input data waktu tiba pada media homogen antara dengan menggunakan hanya data gelombang P ataupun dengan data gelombang P dan S menghasilkan posisi hiposenter yang tidak terlalu berbeda. • Input data waktu tiba pada media multilayer (minimal dua lapisan) harus menggunakan data waktu tiba gelombang P dan S. Posisi hiposenter hasil inversi pada media multilayer (minimal dua lapisan) dengan input hanya gelombang P menghasilkan Hal :6

error signifikan untuk posisi kedalaman (error > 30%). Daftar Pustaka Andri, 2006, Pemodelan Kedepan Struktur Kecepatan di Bawah Kompleks Krakatau dengan Menggunakan Solusi Persamaan Gelombang dan Persamaan Eikonal, Teknik Geofisika, ITB (tidak dipublikasikan). Nishi, K., 2001, A three-dimensional robust seismic ray tracer for volcanic regions, Earth Planets Space, 53, 101-109. Qin, F., Luo, Y., Olsen, K.B., Cai, W., Schuster, G.T., 1992, Finite-difference solution of the eikonal equation along expanding wavefronts, Geophysics, 57, p.478-487. Vidale, J. E., 1988, Finite-difference traveltime calculation: Bull., Seis. Soc. Am., 78, 2062-2076. Vidale, J.E., 1990, Finite-difference calculation traveltime in three diemnsions, Geophysics, 55, 512-526.


SIMULASI PERHITUNGAN WAKTU TEMPUH GELOMBANG DENGAN METODA EIKONAL : SUATU CONTOH APLIKASI DALAM ESTIMASI KETELITIAN HIPOSENTER GEMPA

Recommend Documents