J. Geofisika Vol. 13 No.1/2012
Penentuan Hiposenter Gempa Mikro Menggunakan Metode Inversi Simulated Annealing pada Lapangan Geotermal “RR” Rexha Verdhora Ry, Andri Dian Nugraha Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan, Institut Teknologi Bandung, Jalan Ganesa 10 Bandung, 40132 Email.
[email protected]
Abstrak Penentuan lokasi hiposenter melibatkan proses pencarian solusi hiposenter yang memiki selisih antara waktu tempuh observasi dan kalkulasi minimum. Ketika memecahkan permasalahan inversi non-linier ini, teknk optimisasi lokal dapat dengan mudah menghasilkan solusi dengan meminimumkan fungsi error, namun fungsi tersebut bergantung pada model awal dan belum tentu mencapai minimum global. Metode lain seperti simulated annealing dapat diterapkan untuk masalah optimisasi global. Sebelumnya, lokasi hiposenter pada lapangan geotermal “RR” telah ditentukan menggunakan metode Geiger.Pada studi ini, metode simulated annealing diterapkan pada data gempa mikromenggunakan model kecepatan 1-D untuk merelokasi hiposenter dan meminimalisasi fungsi error. Waktu tempuh dihitung menggunakan ray tracing dengan metode shooting.Hasil penelitian ini menunjukkan bahwa lokasi hiposenter memiliki errorRMS yang lebih kecil dibandingkan dengan penelitian sebelumnya, yang secara statistik berasosiasi dengan solusi yang lebih baik. Kata kunci: penentuan hiposenter, inversi, metode Geiger, simulated annealing, metode shooting. Abstract Hypocenter determination involves finding a hypocenter location that has minimum difference between the observed and the calculated travel times. When solving this nonlinear inverse problem, a local optimization technique can easily produce a solution for which minimizes error function, but its function itself depends on initial model and does not necessarily reach its global minimum. Other methods such as simulated annealing can be applied to such global optimization problems. Unlike local methods, the convergence of the simulated annealing method is independent on the initial model. Previously, hypocenter location at “RR” Geothermal Field has been determined by using Geiger’s method. However, in this paper, simulated annealing method was applied on same data and 1-D velocity model to relocate hypocenter and minimize error function. The travel times were calculated using ray tracing with shooting method. Our results show, that hypocenter location has smaller Root Means Square (RMS) error compared to the previous study, can be statistically associated with better solution. Keywords:hypocenter determination, inversion, Geiger’s method, simulated annealing, shooting method.
1. Pendahuluan Pengamatan aktivitas gempa mikro dari suatu area geotermal dapat digunakan untuk mendeteksi permeabilitas struktur. Perubahan permeabilitas struktur dapat disebabkan oleh tekanan pori dan perubahan suhu akibat interaksi antara fluida reservoir yang bersirkulasi dengan hot rock, ataupun oleh proses stimulasi pada reservoir geotermal yang berhubungan dengan injeksi fluida. Agar dapat menginterpretasi proses perekahan tersebut dan menganalisa persebaran sumber gempa mikro, lokasi hiposenter gempa mikro perlu ditentukan terlebih dahulu. Penentuan lokasi hiposenter ini melibatkan proses inversi untuk mencari posisi hiposenter yang memiliki selisih (misfit) minimum antara waktu tempuh observasi dan kalkulasi.
Penelitian sebelumnya telah dilakukan untuk menentukan lokasi hiposenter pada lapangan geotermal “RR” menggunakan metode Geiger.Metode Geiger sendiri mengimplementasikan algoritma iterative least square untuk memperoleh misfit minimum. Teknik optimalisasi lokal seperti ini dapat menghasilkan solusi dengan meminimalkan fungsi misfit, namun proses tersebut sangat bergantung pada model awal.Penelitian tersebut mengidentifikasi 591 event gempa mikro.Data tersebut direkam oleh 19 stasiun dan terjadi selama satu tahun di lapangan geotermal “RR” dari periode Januari 2007 sampai Desember 2007.Gambar 1 menunjukkan persebaran hiposenter dan stasiun hasil penelitian tersebut.
23
J. Geofisika Vol. 13 No.1/2012
Menurut penelitian tersebut, gempa mikro pada lapangan geotermal “RR” terjadi pada rentang antara elevasi 0,75 km sampai -8,43 km (MSL = 0 km). Apabila dirincikan, persebaran event dapat dibagi menjadi 4 zona yaitu zona dangkal pada rentang antara elevasi 1 km sampai 0 km, zona menengah pada rentang antara elevasi 0 sampai -6,5 km, zona
dalam pada rentang antara elevasi -6,5 km sampai -9 km, dan zona jauh dimana event terjadi jauh dari persebaran event lainnya. Pada zona dangkal terdapat 67 event gempa mikro, sedangkan pada zona jauh terdapat 9 event gempa mikro, dari total 591 event. Persebaran gempa mikro ini ditunjukkan oleh Gambar 1.
U
MSL
Gambar 1. Peta persebaran eventgempa mikro yang dihasilkan oleh penelitian sebelumnya (bulat hijau) dan stasiun (segitiga biru terbalik) pada lapangan geotermal “RR”
24
hiposenter yang memiliki RMS error lebih kecil akan berasosiasi dengan solusi yang lebih baik secara statistik. 2 Vp Vs
1 0 -1 -2 Kedalam an (km )
Setiap metode inversi memiliki kelebihan dan kekurangannya masing-masing.Solusi hiposenter tersebut bisa saja terjebak pada minimum lokal karena model awal yang kurang baik.Metode inversi lain yaitu simulated annealing dapat diterapkan untuk kasus optimisasi global. Tidak seperti metode optimisasi lokal seperti iterative least square, kekonvergenan dari metode simulated annealing tidak bergantung pada model awal. Pada penelitian ini, metode tersebut digunakan pada data dan model kecepatan 1D (Gambar 2) yang sama dengan penelitian sebelumnya untuk menentukan lokasi hiposenter dan meminimumkan fungsi misfit. Kedua solusi akan dibandingkan dan lokasi
-3 -4 -5 -6 -7 -8 -9
0
2 4 6 8 Kecepatan Gelombang Seismik (km/s)
J. Geofisika Vol. 13 No.1/2012
Gambar 2. Plot model kecepatan untuk gelombang P (Vp) dan gelombang S (Vs) yang digunakan pada studi ini.
2. Metodologi
konfigurasi sistem dan fungsi objektif (misfit) sebagai energi. Faktor temperatur merupakan faktor pengontrol dengan satuan sama dengan fungsi objektif. Probabilitas perturbasi model dinyatakan oleh:
2.1 Ray Tracing Metode Shooting Untuk menentukan waktu tempuh gelombang P dan S yang merambat dari sumber gempa ke stasiun penerima melalui model kecepatan seismik 1-D (Gambar 2) dilakukan perhitungan menggunakan metode ray tracing metode shooting. Metode ini menerapkan hukum Snell untuk menentukan sudut datang dan transmisi pada setiap batas lapisan kecepatan gelombang seismik, diformulasikan oleh persamaan:
=
∅ ∅
(1)
2.2 Simulated Annealing Salah satu metode guided random search atau pencarian acak terarah adalah metode simulated annealing (SA). Metode simulated annealing (Grandis,2009) dalam inversi didasarkan pada proses termodinamika pembentukan kristal suatu substansi. Pada temperatur tinggi suatu substansi padat mencair, kemudian proses pendinginan secara perlahan-lahan menyebabkan terbentuknya kristal yang berasosiasi dengan energi sistem yang minimum. Pada proses simulated annealing ini, ruang model harus didefinisikan terlebih dahulu dengan menentukan secara “a priori” interval harga minimum dan maksimum parameter model, dalam penelitian kali ini parameter modelnya adalah posisi gempa. Pemilihan harga parameter model ditentukan secara acak sebagai bilangan sembarang dalam interval nilai minimum dan maksimum masing-masing. Caranya adalah mengambil bilangan acak R dengan probabilitas uniform antar 0 dan 1 yang dipetakan menjadi harga parameter model. Proses pembentukan kristal (annealing) dalam termodinamika diadopsi dalam penyelesaian masalah inversi dengan menggunakan parameter model untuk mendefinisikan
P(∆E) =exp −
∆
(2)
dimana ∆E adalah perubahan fungsi objektif atau perubahan misfit akibat perturbasi model tersebut. Jika ∆E < 0, maka perubahan model diterima. Namun jika ∆E≥0, maka penentuannya ditentukan secara probabilistik menggunakan bilangan acak R yang terdistribusi uniform pada interval [0,1]. Jika R < P(∆E) maka perubahan diterima, sebaliknya jika R≥P(∆E) perubahan ditolak dan kembali kemodel sebelumnya.Proses iterasi dimulai dengan temperatur cukup tinggi sehingga hampir semua perturbasi model diterima. Pada saat temperatur turun secara perlahan, perturbasi model yang diterima akan mengecil jika ∆E≥0. 3. Algoritma Ray tracing pada penelitian ini menggunakan kode program MATLAB yang dikembangkan oleh penulis. Cara kerja pemograman ray tracing tersebut adalah melakukan iterasi untuk menebak sudut datang pada batas lapisan yang kemudian menghasilkan sudut transmisi berdasarkan hukum Snell dan iterasi terus dilakukan sampai ray mengenai stasiun. Untuk mempermudah, jarak antara sumber gempa ke penerima dianggap sebagai suatu bidang lurus 1-D. Apabila ray telah berhasil mengenai stasiun, maka ray tersebut digunakan untuk menghitung panjang ray dan waktu tempuh gelombang di setiap lapisan kecepatan. Kemudian, waktu tempuh dari sumber gempa ke stasiun penerima adalah jumlah dari seluruh waktu tempuh gelombang tersebut. Penentuan lokasi hiposenter melibatkan proses inversi untuk mencari posisi hiposenter yang memiliki selisih (misfit) minimum antara waktu tempuh observasi dan kalkulasi. Metode inversi yang digunakan pada penelitian ini ialah metode simulated annealing. Proses inversi ini menggunakan kode program MATLAB yang dikembangkan 25
J. Geofisika Vol. 13 No.1/2012
oleh penulis. Algoritma dikembangkan agar bersifat cepat dan adaptif berdasarkan White (1984) dan Weber (2000). Untuk itu, beberapa faktor yang perlu diperhatikan ialah temperatur awal, formulasi penurunan temperatur, perturbasi model, jumlah simulasi pada suatu temperatur, dan kriteria pemberhentian iterasi. Temperatur awal harus dipilih sekecil mungkin agar proses inversi berlangsung cepat, namun dibuat sedemikian mungkin agar dapat tetap menerima dan mencoba banyak percobaan simulasi. Beberapa prosedur trial and error dilakukan untuk memperkirakan nilai optimum dari temperatur awal tersebut.Akhirnya diperoleh kesimpulan bahwa nilai optimum dari temperatur awal dapat ditentukan dari konfigurasi parameter yang ditentukan secara acak pada simulasi pertama. Nilai temperatur awal dipilih lebih besar atau sama dengan standar deviasi dari misfit simulasi pertama. Maka karena itu, nilai temperatur awal akan berubah (tidak tetap) setiap kali proses inversi dilakukan pada penelitian ini. Formulasi penurunan temperatur dipilih agar beberapa simulasi terus dilakukan pada setiap temperatur Tk sampai sistem stabil dan memperoleh solusi. Pada penelitian ini, penulis menggunakan formulasi penurunan sederhana yaitu Tk+1 = 0.99 Tk. Suatu simulasi menggambarkan suatu percobaan parameter model yang dilakukan. Pada setiap simulasi, suatu parameter model diperbarui dengan pemilihan secara acak kemudian dihitung waktu tempuh kalkulasinya, atau dapat disebut dengan perturbasi model. Berdasarkan Baghmishes dan Navarbaf (2008), model akan diperbarui sesuai dengan persamaan: = 0,1 − 0.5 .! 1 + # |& '()*,# +#| $
− 1]
(3)
sehingga pada simulasi ke-k dihasilkan dk yang bernilai [-1,1] dan besar pembaruan model di simulasi ke-k adalah:
26
∆, = , -, − , - (4) Pada persamaan 4 , interval [xmax , xmin] merupakan rentang nilai mungkin yang diizinkan pada simulasi ke-k. Pada proses inversi di penelitian ini, jumlah simulasi di temperatur Tk digambarkan oleh variabel Ak dan Lk. Nilai Akmenggambarkan banyaknya model baru yang diterima dan nilai Lk menggambarkan banyaknya perturbasi model yang dilakukan. Nilai Ak akan bertambah setiap ada model baru yang diterima, sampai pada batas Amin. Nilai Amin ditentukan sehingga terdapat jumlah minimum model baru yang diterima pada temperatur Tk, dimana Amin adalah angka bulat. Kemudian, pada temperatur yang sudah sangat kecil (mendekati nol), model baru yang diterima hampir tidak ada sehingga perlu adanya variabel lain untuk melanjutkan iterasi, yaitu variabel Lk. Nilai Lk akan bertambah setiap perturbasi model dilakukan, sampai pada batas Lmax. Nilai Lmax ditentukan sehingga terdapat jumlah maksimum perturbasi model yang dilakukan pada temperatur Tk, dimana Lmax adalah angka bulat. Faktor terakhir yang perlu diperhatikan ialah kriteria pemberhentian iterasi, digambarkan oleh variabel U. Nilai U adalah nol jika pada temperatur Tk terdapat pembaruan model (ada model baru yang diterima, Ak≥ 1). Nilai Uakan bertambah jika tidak ada pembaruan model di temperatur tersebut. Batas nilainya adalah Umax, dan kemudian prosedur berhenti. Dalam kata lain, semua prosedur akan berhenti jika konfigurasi sistem tidak berubah (tidak ada parameter model baru yang diterima, Ak == 0) untuk simulasi sebanyak UmaxxLmaxkali. Berdasarkan algoritma di atas, parameter annealing yang digunakan pada penelitian ini adalah: Amin = 10, Lmax = 15, dan Umax = 20. Gambar 3 menggambarkan algoritma keseluruhan prosedur inversi yang dilakukan sampai memperoleh solusi.
J. Geofisika Vol. 13 No.1/2012
Gambar 3. Diagram alir metode inversi simulatedannealing pada penentuan lokasi hiposenter pada studi ini.
model kecepatan gelombang seismik 1-D dan posisi stasiun yang sama dengan data riil. Data sintetis berupa waktu tempuh gelombang dibuat dengan menentukan suatu hiposenter yang telah diketahui posisinya, kemudian dihitung waktu tempuhnya ke setiap stasiun melalui model kecepatan 1-D. Data sintetis ini ditunjukkan pada Tabel 1.
4. Hasil dan Diskusi 4.1 Data Sintetis Tujuan dari pembuatan dan pengolahan data sintetis adalah untuk menguji kebenaran dan kehandalan dari algoritma yang telah dibuat. Pembuatan data sintetis ini menggunakan
Tabel 1. Katalog data sintetis waktu tempuh gelombang Posisi Hiposenter X (km)
Y (km)
11
11
Posisi Stasiun Z (km)
-3
No.
X (km)
Y (km)
Z (km)
Waktu Tempuh (s)
1
16,61666
9,3853
1,335
1,561148715
2
16,51366
11,5613
1,42
1,540462621
3
12,29366
11,4463
1,05
0,943305874
4
14,31366
9,6153
1,22
1,211803692
5
14,00566
8,2193
1,097
1,248522213
Algoritma inversi menggunakan metode simulated annealing kemudian diterapkan pada data sintetis tersebut. Perhitungan
dilakukan sebanyak 20 kali untuk melihat apakah solusi yang diperoleh stabil atau tidak.Hasil inversi ditunjukkan oleh Tabel 2.
Tabel 2. Katalog solusi hasil inversi pada data sintetis Percobaan ke-
Posisi Hiposenter X (km)
Y (km)
Z (km)
error RMS (s)
1
11,03191
10,99055
-2,99346
0,002495
2
11,02536
11,04477
-3,01915
0,003215
3
11,06267
11,06747
-3,0462
0,002582
27
J. Geofisika Vol. 13 No.1/2012
4
11,01087
11,03581
5
11,05361
11,01492
-3,0276
0,002157
6
11,05611
11,01796
-3,02093
0,001617
7
11,01033
11,00713
-3,00043
0,002506
8
11,02544
10,99621
-2,99265
0,002915
9
11,01741
11,00869
-3,02176
0,002446
10
11,03251
11,00522
-3,02504
0,001925
11
11,0046
10,98935
-2,98496
0,00305
12
11,00836
10,95518
-3,0037
0,002366
13
11,00018
10,98141
-2,98632
0,00247
14
11,04784
10,97016
-3,06595
0,003328
15
11,00541
11,00144
-2,97761
0,002881
16
11,0129
10,9679
-3,00386
0,002379
17
11,0188
10,97034
-3,01281
0,002822
18
10,97239
10,98823
-2,9635
0,002069
19
11,0188
11,02949
-2,99187
0,002202
20
11,04076
10,98443
-3,04058
0,002941
Dapat dilihat algoritma inversi simulatedannealing yang digunakan menghasilkan solusi yang stabil setelah 20 kali perhitungan. Solusi yang dihasilkan mendekati solusi sebenarnya dan memiliki nilai error RMS yang kecil. Maka, dapat disimpulkan bahwa solusi yang dihasilkan oleh algoritma tersebut memiliki presisi dan akurasi yang tinggi. 4.2 Data Riil Data 591 event gempa mikro dan model kecepatan 1-D dari penelitian sebelumnya akan digunakan pada penelitian ini. Berikut hasil dari penelitian ini menggunakan metode inversi simulated annealing. Gempa mikro pada lapangan geotermal “RR” (591 event) terjadi pada rentang antara elevasi 0.55 sampai -8.67 km (MSL = 0 km). Berdasarkan pembagian 4 zona di bagian pendahuluan, persebaran event dapat dirincikan sebagai berikut. Pada zona dangkal, terdapat 18 event gempa mikro.Pada zona menengah, 556 event gempa mikro.Pada zona dalam, terdapat 9
28
-2,98444
0,003117
event gempa mikro. Kemudian pada zona jauh, terdapat 9 event gempa mikro. Persebaran gempa mikro hasil penelitian ini ditunjukkan oleh Gambar 4. Pada penelitian sebelumnya, gempa mikro yang berada di zona dangkal diidentifikasi sebanyak 67 event.Sedangkan pada penelitian ini, gempa mikro yang berada di zona dangkal diidentifikasi sebanyak 18 event. Gempagempa yang sebelumnya diidentifikasi terjadi di kedalaman dangkal ternyata terjadi di kedalaman yang lebih dalam. Solusi yang dihasilkan dari penelitian ini dapat dianggap lebih baik. Pada zona jauh, gempa mikro yang teridentifikasi dari hasil inversi tidak jauh berbeda. Sebanyak 9 event yang berada di zona jauh dari hasil penelitian sebelumnya tidak berpindah posisi secara signifikan pada penelitian ini. Solusi yang dihasilkan pada zona tersebut belum dapat menghasilkan perbaikan yang signifikan.
J. Geofisika Vol. 13 No.1/2012
U
MSL
Gambar 4. Peta persebaran event yang dihasilkan oleh penelitian ini menggunakan simulated annealing (bulat merah) dan stasiun (segitiga biru terbalik).
Solusi hiposenter yang dihasilkan pada penelitian ini memiliki sebaran nilai RMS error yang lebih kecil dibandingkan solusi hiposenter pada penelitian sebelumnya.Hal ini ditunjukkan oleh Gambar 5. Dapat disimpulkan bahwa metode simulated annealing menghasilkan solusi yang lebih baik secara statistik.
permeabel ataupun rekahan. Dibutuhkan studi lanjut untuk menganalisa struktur-struktur tersebut, baik kajian geologi ataupun kajian geofisika lainnya. Salah satu studi lanjut yang bisa dilakukan untuk menganalisa gempa mikro tersebut ialah permodelan tomografi seismik 3-D untuk mencitrakan struktur bawah permukaan pada lapangan panas bumi “RR”.
Persebaran gempa mikro pada lapangan panas bumi “RR” akan berhubungan dengan struktur
29
J. Geofisika Vol. 13 No.1/2012
a)
b)
Gambar 5.Histogram (a) distribusi waktu residual dan (b) error RMS solusi hiposenter gempa mikro yang dihasilkan dari penelitian sebelumnya (kiri) dan simulated annealing pada penelitian ini (kanan).
5.
Kesimpulan
Dapat disimpulkan bahwa metode inversi simulated annealing merupakan metode inversi yang sangat baik dan berguna untuk diterapkan pada kasus penentuan lokasi hiposenter. Kekonvergenan dari metode ini tidak bergantung pada model awal. Metode simulated annealing menghasilkan solusi lokasi hiposenter dengan nilai RMS error yang lebih kecil dibandingkan dengan metode Geiger. Selain itu, metode Geiger menghasilkan banyak gempa mikro yang berada pada zona dangkal (elevasi antara 1 - 0 km). Sedangkan solusi metode simulated annealing pada event yang sama, gempa mikro ternyata terjadi lebih dalam dan dominan berada pada zona menengah (elevasi < 0 km). Hal ini bisa saja terjadi disebabkan oleh solusi yang dihasilkan metode Geiger 30
terjebak pada minimum lokal karena model awal yang kurang baik. Ke depannya, solusi hiposenter gempa mikro ini dapat digunakan untuk menganalisa struktur permeabel ataupun rekahan pada lapangan geotermal “RR”. Data gempa mikro ini dapat digunakan sebagai input pada permodelan tomografi seismik 3-D untuk mencitrakan struktur bawah permukaan. Ucapan Terima Kasih Penulis menghaturkan terimakasih ke PT Chevron Geothermal Indonesia atas penggunaan data dalam penelitian ini; ke Laboratorium Geofisika Dekat Permukaan, Teknik Geofisika, FTTM, ITB atas fasilitas yang mendukung penelitian ini.
J. Geofisika Vol. 13 No.1/2012
Daftar Pustaka Grandis, H., 2009. Pengantar Permodelan Inversi Geofisika. Penerbit Himpunan Ahli Geofisika Indonesia, pp. 93-100. Vakil-Baghmishes, M. dan Navarbaf, A., 2008. Modified Very Fast Simulated Annealing Algorithm. International Symposium on Telecommunications. Weber, Z., 2000. Seismic travel time tomography: a simulated annealing approach. Elsevier. Physics of Earth and Planetary Interiors, 199, pp 149159. White, S. R., 1984. Concepts of scale in simulated annealing. Proc. IEEE Int. Conference on Computer Design, Port Chester, pp. 646–651.
31
