JURNAL ISSN: 0854-4352 Vol. 13 No.1/2012 GE FISIKA Publikasi Ilmiah Sains dan Teknologi Geofisika, Meteorologi, Oseanografi, Geologi, Geodesi Pemograman Ray Tracing Metode Pseudo-Bending Medium 3-D Untuk Menghitung Waktu Tempuh Antara Sumber Dan Penerima Interpretasi Data Anomali Medan Magnetik Total Untuk Pemodelan Struktur Bawah Permukaan Daerah Manifestasi Mud Vulcano (Studi Kasus Bledug Kuwu Grobongan) Penentuan Hiposenter Menggunakan Simulated Annealing dan Guided Error Serta Penentuan Model Kecepatan Gelombang Seismik 1-D pada Lapangan “Geothermal” Penentuan Hiposenter Gempa Mikro Menggunakan Metode Inversi Simulated Annealing pada Lapangan Geotermal “RR” Pencitraan Struktur 3-D Vp, Vs, Rasio Vp/Vs Menggunakan Tomografi Double Difference di Willayah Bali Himpunan Ahli Geofisika Indonesia Indonesian Association of Geophysicists Jurnal GEOFISIKA Jurnal GEOFISIKA adalah majalah ilmiah yang diterbitkan oleh Himpunan Ahli Geofisika Indonesia (HAGI). Jurnal ini merupakan jurnal yang bereferensi, yaitu setiap judul makalah dinilai dan ditelaah oleh minimal 2 (dua) orang penelaah yang ahli dalam bidang yang bersangkutan dan berasal dari institusi perguruan tinggi dan industri baik dari dalam maupun luar negeri. Jurnal GEOFISIKA diterbitkan secara berkala setiap tahun dengan dua edisi penerbitan. Tim Editor Penasehat 1. Prof. Sri Widiyantoro (Presiden HAGI) 2. Mailendra (Sekjen HAGI) Pimpinan Editor Irwan Meilano (ITB) Sekretaris Zulfakriza Zulhan (ITB) Anggota Abdul Haris (UI) Hamzah Latief (ITB) Wiwit Suryanto (UGM) Prof. Phil Cummins (ANU) Erdinc Saygin (ANU) Makky Jaya (GFZ Postdam) Eddy Arus Sentani (PND) Cecep Rudiana (Hess) Alamat Sekretariat: Gedung Patra Office Tower Lt. 20 Suite 2045 Jl. Jend. Gatot Subroto Kav 32-34 Jakarta 12950 Phone/Fax. 021-52500040 Email: secretariat@hagi.or.id Website: hub.hagi.or.id Daftar Isi Vol. 13 No. 1/2012  Daftar Isi ..…...……………………………….i  Pengantar Editor …………………………….ii  Pemograman Ray Tracing Metode PseudoBending Medium 3-D Untuk Menghitung Waktu Tempuh Antara Sumber Dan Penerma Ahmad Syahputra dan Andri Dian Nugraha …………………………...1  Interpretasi Data Anomali Medan Magnetik Total Untuk Pemodelan Struktur Bawah Permukaan Daerah Manifestasi Mud Vulcano (Studi KasusBledug Kuwu Grobongan) Sigit Darmawan, Hernowo Danusaputro, dan Toni Yulianto ….……………………….7  Penentuan Hiposenter Menggunakan Simulated Annealing dan Guided Error Serta Penentuan Model Kecepatan Gelombang Seismik 1-D pada Lapangan “Geothermal” Akhmad Fanani Akbar, Andri Dian Nugraha, M. Rachmad Sule, dan Aditya A. Juanda ……………...…………...16  Penentuan Hiposenter Gempa Mikro Menggunakan Metode Inversi Simulated Annealing pada Lapangan Geotermal “RR”. Rexha Verdhora Ry dan Andri Dian Nugraha ………………….……23  Pencitraan Struktur 3-D Vp, Vs, Rasio Vp/Vs Menggunakan Tomografi Double Difference di Willayah Bali Putu Kusuma Yadnya, Andri Dian Nugraha dan Supriyanto Rohadi ….………………....32 Keterangan gambar sampul depan: Diagram skematik interpretasi terhadap hasil relokasi dan hasil tomogram Vp di bawah gunung Agung Bali i PENGANTAR EDITOR Para pembaca yang terhormat, Tim editor Jurnal GEOFISIKA mamanjatkan puji dan syukur pada Tuhan Yang Maha Kuasa atas terbitnya Jurnal GEOFISIKA Vol. 13 No.1 Tahun 2012. Pada edisi kali ini, Jurnal GEOFISIKA menyajikan makalah-makalah ilmiah yang telah melewati proses penilaian dan pengeditan sehingga trepilih lima makalah dengan bidang kajian dan metodologi yang beragam. Lima makalah yang dimuat menampilkan kajian yang menarik untuk dibaca dan dicermati. Kelima makalah tersebut adalah :  Pemograman Ray Tracing Metode Pseudo-Bending Medium 3-D Untuk Menghitung Waktu Tempuh Antara Sumber Dan Penerima, yang ditulis oleh Ahmad Syahputra dan Andri Dian Nugraha.  Interpretasi Data Anomali Medan Magnetik Total Untuk Pemodelan Struktur Bawah Permukaan Daerah Manifestasi Mud Vulcano (Studi KasusBledug Kuwu Grobongan), yang ditulis oleh Sigit Darmawan, Hernowo Danusaputro, dan Toni Yulianto  Penentuan Hiposenter Menggunakan Simulated Annealing dan Guided Error Serta Penentuan Model Kecepatan Gelombang Seismik 1-D pada Lapangan “Geothermal”, yang ditulis oleh Akhmad Fanani Akbar, Andri Dian Nugraha, M. Rachmad Sule, dan Aditya A. Juanda  Penentuan Hiposenter Gempa Mikro Menggunakan Metode Inversi Simulated Annealing pada Lapangan Geotermal “RR”, yang ditulis oleh Rexha Verdhora Ry dan Andri Dian Nugraha  Pencitraan Struktur 3-D Vp, Vs, Rasio Vp/Vs Menggunakan Tomografi Double Difference di Willayah Bali, yang ditulis oleh Putu Kusuma Yadnya, Andri Dian Nugraha dan Supriyanto Rohadi Terakhir Tim Editor berharap kelima makalah yang dimuat pada Jurnal GEOFISIKA edisi kali ini dapat memperkaya khasanah keilmuan geofisika di Indonesia dan bermanfaat untuk menambah wawasan baru bagi para pembaca. Tim Editor J. GEOFISIKA ii J. Geofisika Vol. 13 No.1/2012 Pemograman Ray Tracing Metode Pseudo-Bending Medium 3-D Untuk Menghitung Waktu Tempuh Antara Sumber Dan Penerima Ahmad Syahputra dan Andri Dian Nugraha Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan Institut Teknologi Bandung, Jalan Ganesha 10 Bandung, 40132 Email: ahmadsy_07@student.itb.ac.id Abstrak Perhitungan waktu tempuh gelombang melewati suatu medium dari sumber ke stasiun penerima menggunakan rekonstruksi jejak sinar gelombang dikenal dengan ray tracing. Salah satu prinsip dari perambatan gelombang adalah prinsip Fermat yang menyatakan bahwa gelombang merambat melewati suatu medium dengan waktu tempuh tercepat. Kami telah mengembangkan dan menguji pemograman aplikasi metode pseudo-bending melewati beberapa model kecepatan 3-D. Diantaranya adalah model anomali kecepatan rendah, model anomali kecepatan tinggi, dan model dengan perubahan kecepatan secara bergradasi serta model acak untuk menguji kestabilan metode ini. Lintasan perambatan gelombang selalu berusaha melewati medium dengan kecepatan lebih tinggi dengan waktu tempuh minimum. Metode ini membutuhkan waktu yang singkat dalam perhitungan sehingga mendukung studi-studi inversi tomografi 3-D seperti tomografi gempa lokal dan tomografi gempa mikro untuk aplikasi di bidang geotermal. Kata kunci: Jejak Sinar Gelombang, PrinsipFermat, Model Kecepatan 3-D. Abstract The calculation of seismic wave travel time through a medium from source to receiver using the wave ray reconstuction is known as ray tracing. One of the wave propagation principle is Fermat's principle which stated that the wave propagates through a medium using the fastest travel time path.We have developed and tested the application programming of pseudo-bending through several 3-D velocity models. There are low velocity anomaly model, high velocity anomaly model, and gradation change velocity model and random velocity model for the stability test of this method. Raypath of wave propagation is always trying to through the high velocity medium with minimum travel time. This method requires a short time in the calculation that support the 3-D tomographic inversion studies, such as local earthquake tomography and micro-earthquake tomography in the geothermal field. Keywords: Ray tracing, Fermat’s Principle, 3-D velocity model, pseudo-bending 1. Pendahuluan Perhitungan waktu tempuh gelombang dan rute dari jejak sinar gelombang dibutuhkan dalam tahapan pemodelan kedepan dalam berbagai studi geofisika, seperti inversi tomografi seismik dan relokasi hiposenter gempa bumi atau gempa mikro. Banyak metode yang telah dikembangkan untuk menghitung waktu tempuh dan jejak sinar gelombang, pseudo-bending (Um dan Thurber, 1987) merupakan algoritma yang cepat dalam waktu komputasinya dibandingkan metode lainnya. Dalam studi ini telah dilakukan pemograman dalam bahasa Matlab untuk proses ray tracing metode pseudo-bending (Um dan Thurber, 1987) dalam medium 3-D yang merupakan pengembangan dari studi ray-tracing 2-D (Syahputra, A., 2011; Nugraha, A. D. dkk., 1 J. Geofisika Vol. 13 No.1/2012 2011) dengan beberapa modifikasi dalam proses komputasinya. Tujuan dari studi ini yaitu untuk membuat pemograman ray-tracing yang dapat diaplikasikan pada inversi tomografi skala lokal 3-D untuk keperluan tomografi gempa lokal ataupun gempa mikro untuk memperoleh gambaran struktur kecepatan gelombang seimik 3-D di lapangan geotermal dan gunungapi. 2. Metode Persamaan integral sebagai ekspresi waktu tempuh (T) sepanjang lintasan gelombang dari sumber hingga penerima (Thurber, 1993), sebagai berikut: (1) dimana dl = segmen panjang lintasan gelombang dan V = kecepatan medium pada titik lintasan yang dilewati sinar gelombang. Lintasan sinar gelombang dapat didiskritasi dalam sebanyak n, jumlah titik bending + 2. seperti yang ditunjukkan pada , , ... , pada Gambar 1. Setelah ditekuk sepanjang Rc arah dengan tanpa mengubah posisi dan didapat titik lintasan yang baru . minimum (Fermat’s Principle). Hasil pertubasi pertama menjadi lintasan awal dan pada pertubasi ini , yang kemudian arah tekukan dan sejauh Rc dihitung kembali. Proses pertubasi ini diulang terus hingga mencapai konvergensi dan waktu tempuh minimum. merupakan vektor anti normal dari vektor titik ke titik . Vektor ini paralel dengan arah gradient kecepatan ( ) pada diturunkan dari medium 3 dimensi. hubungan persamaan sebagai berikut : ! "# $! "%&'()* '(+* ,-&'()* '(+*, .'()* '(+* ./ (2) Dan jarak Rc dihitung dengan rumus sebagai berikut : 1 23 01 # 451 % ! " 6 2 3" 8 !1 27 451 % ! " 6 2 dimana L = . (+* "8: # :1 9 ; 8 (3) . dan c = ! ()* 2 Sehingga didapat titik lintasan sinar gelombang yang baru, sebagai berikut : 2 01 (4) dimana . . Gambar 1. Ilustrasi dari skema tiga titik pertubasi ( , , ). Proses ray tracing berawal dari sinar dan adalah gelombang antara titik lintasan garis lurus. Kemudian titik tengah antara kedua titik ini, (pada pertubasi pertama ) ditekuk ke arah sejauh Rc. Kemudian skema tiga titik pertubasi ini (Gambar 2) diaplikasi ke sepanjang lintasan sinar gelombang yang mengalami gangguan tetapi belum mencapai waktu tempuh 2 Gambar 2. Skema urutan titik pertubasi dari kiri ke kanan yang digunakan dalam pemograman ray tracing pada studi ini. 2.1 Algoritma Berangkat dari posisi sumber dan penerima, kemudian parameterisasi model dilakukan. Parameterisasi model yang digunakan pada J. Geofisika Vol. 13 No.1/2012 penelitian ini berupa blok 3-D. Ukuran blok ditentukan tergantung dari tingkat heterogenitas dalam arah vertikal dan arah horizontal dari model kecepatan. Setelah membuat model kecepatan sesuai dengan parameterisasi model kemudian dihitung gradien kecepatan yang merupakan turunan kecepatan terhadap jarak spasial arah X, Y dan Z untuk kasus 3-D. < >? <= ? 2 <@ A 2 CD <B < < (5) Salah satu tahapan penting dalam algortima ini adalah penentuan jumlah titik tekuk dan banyaknya pertubasi. Tahapan ini sangat berpengaruh terhadap optimasi waktu komputasi dalam mencapai nilai konvergensinya. Penentuan kedua nilai tersebut dipengaruhi oleh parameterisasi model yang berujung dengan tingkat heterogenitas model awal. Untuk memudahkan dalam pemograman ray tracing metode pseudo-bending ini, kami membuat diagram alir yang ditunjukkan pada Gambar 3. Pada pemograman ini, ray tracing berawal dengan lintasan ray lurus. Kemudian lintasan ray yang lurus ini diberi gangguan arah sejauh Rc pada setiap titik tekuknya. Lintasan ray diperbaharui sebanyak jumlah pertubasi. Masing-masing ray hasil setiap pertubasi dihitung pajangnya pada setiap blok dengan cara membagi ray tersebut menjadi segmensegmen yang lebih kecil. Semakin kecil segmennya semakin tinggi tingkat ketelitian dalam menghitung ray pada setiap blok. Waktu tempuh gelombang merambat dihitung dengan mengalikan panjang ray setiap blok dengan nilai slowness (1/kecepatan) pada setiap blok. EFCGH IJKHL M NO 9O (6) Dimana Sf adalah slowness pada blok ke-f yang dilewati oleh ray. dLf merupakan segmen panjang ray pada blok ke-f yang dilewati ray. Kemudian dari waktu tempuh masing-masing pertubasi pada ray tracing dipilih waktu minimumnya. Pertubasi ke-i dengan waktu minimum ini menjadi ray tracing akhir yang memenuhi prinsip Fermat. Panjang ray dalam setiap blok ini merupakan komponen baris matriks tomografi dalam inversi tomografi. 3. Hasil Uji Ray-Tracing Dalam menguji pemograman ray tracing 3-D yang telah dibuat pada penelitian ini, kami membuat beberapa model kecepatan. Diantaranya model kecepatan anomali positif (Gambar 4) dan anomali negatif (Gambar 6), model kecepatan gradasi terhadap kedalaman (Gambar 8, 10 dan 12), dan model kecepatan acak (Gambar 14). Gambar 3. Diagram alir pemograman ray tracing metode pseudo-bending pada penelitian ini (Nugraha, A. D. dkk., 2011). Volume dari parameterisasi model 20 x 20 x 20 km3 dengan dimensi blok 1 x 1 x 1 km3. Dimana 20 stasiun penerima dengan jarak pisah 4 km arah X dan Y dipasangkan dipermukaan yang merekam 1 sumber di bawah permukaan (Gambar 4). Selain menampilkan jejak sinar yang dilalui gelombang dengan waktu minimum, pada penelitian ini akan memperlihat dampak waktu tempuh untuk mencapai ke semua titik 3 J. Geofisika Vol. 13 No.1/2012 di permukaan dengan interval titik 1 km arah X dan Y. 3.1 Model kecepatan beranomali positif Gambar 4.Jejak sinar gelombang pertubasi (garis hitam) dan jejak sinar gelombang waktu tempuh minimum (garis merah) melewati medium homogen dengan kecepatan 2 km/s dengan anomali kecepatan positif 6 km/s melengkung ke arah anomali kecepatan tinggi (biru). Gambar 6. Jejak sinar gelombang pertubasi (garis hitam) dan jejak sinar gelombang waktu tempuh minimum (garis merah) melewati medium homogen berkecepatan 6 km/s dengan anomali kecepatan negatif 2 km/s melengkung menjauhi anomali kecepatan rendah (merah). Gambar 7. Waktu tempuh gelombang dari sumber ke masing-masing penerima (segitiga hitam). Terdapat zona waktu tempuh tinggi (merah) sebagai dampak anomali kecepatan rendah. 3.3 Model kecepatan gradasi Gambar 5. Waktu tempuh gelombang dari sumber ke masing-masing penerima (segitiga hitam). Terdapat zona waktu tempuh rendah (biru) sebagai dampak anomali kecepatan tinggi. 3.2 Model kecepatan beranomali negatif 4 Gambar 8.Jejak sinar gelombang pertubasi (garis hitam) berawal dari garis lurus kemudian ditekuk hingga jejak sinar gelombang waktu minimum (garis merah). Jejak sinar gelombang final (garis merah) ditekuk ke arah kecepatan yang lebih tinggi pada medium kecepatan gradasi(Vz = 2 + 0.2z) ini sehingga didapat waktu tempuh minimum. J. Geofisika Vol. 13 No.1/2012 Gambar 9.Waktu tempuh gelombang dari sumber ke masing-masing penerima (segitiga hitam) yang diterima di permukaan menunjukkan kontur kesamaan waktu yang melingkar sempurna karena tidak terdapat anomali kecepatan di bawah permukaan. Gambar 10. Waktu tempuh gelombang dari sumber ke masing-masing penerima (segitiga hitam). Terdapat zona waktu tempuh rendah (biru) sebagai dampak anomali kecepatan tinggi. 3.5 Model kecepatan gradasi beranomali 3.4 Model kecepatan gradasi beranomali negatif positif Gambar 10. Jejak sinar gelombang pertubasi (garis hitam) berawal dari garis lurus kemudian ditekuk hingga jejak sinar gelombang waktu minimum (garis merah). Jejak sinar gelombang final (garis merah) ditekuk ke arah daerah anomalikecepatan 6 km/s (biru) pada medium kecepatan bergradasi(Vz = 2 + 0.2z). Gambar 11.Jejak sinar gelombang (garis merah) melewati medium kecepatan gradasi (Vz = 2 + 0.2z) dengan anomali kecepatan negatif 2 km/s melengkung menjauhi anomali kecepatan rendah (merah). Gambar 12. Waktu tempuh gelombang dari sumber ke masing-masing penerima (segitiga hitam). Terdapat zona waktu tempuh tinggi 5 J. Geofisika Vol. 13 No.1/2012 (merah) sebagai dampak anomali kecepatan rendah. 3.6 Model kecepatan acak Gambar 13.Jejak sinar gelombang (garis merah) melewati medium kecepatan acak dengan interval nilai > 0 dan <6 km/s melengkung secara acak dengan arah menuju kecepatan yang lebih tinggi dibandingkan lingkungan sekitarnya. Pada model kecepatan acak ini, metode pseudo-bending tetap stabil. bending ini sangat baik diterapkan dalam rekonstruksi penjejakan sinar gelombang yang memenuhi prinsip fermat dengan waktu tempuh tercepat karena membutuhkan waktu yang lebih singkat dalam perhitungannya dibandingkan metode lain dan stabil dalam menghadapi berbagai medium kecepatan 3-D yang memiliki tingkat heterogenitas bervariasi. Perhitungan nilai Rc kadang dapat memiliki nilai imajiner atau bernilai sangat besar karena hal ini sangat berhubungan dengan gradien kecepatan pada titik jejak lintasan sinar gelombang tersebut. Dalam menjaga kestabilan tekukkan (gangguan) nilai Rc yang dapat diterima jika bernilai 0 – 1 dan jika Rc ditemukan bernilai imajiner maka Rc dianggap bernilai 0 pada pertubasi tersebut (Nugraha, A.D. dkk., 2011). Pemograman ray tracing pseudo bending dalam penelitian ini, dapat diaplikasikan pada inversi tomografi waktu tempuh pada gempa lokal ataupun mikro di bawah gunungapi dan geotermal untuk mendapatkan struktur kecepatan gelombang seismik bawah permukaan. Daftar Pustaka Gambar 14.Waktu tempuh gelombang dari sumber ke masing-masing penerima (segitiga hitam) yang melewati medium kecepatan acak. Waktu tempuh yang dibutuhkan dari sumber mencapai titik-titik di permukaan tidak berpola baik (acak). 4. Kesimpulan Jumlah titik tekuk dan jumlah pertubasi yang berdampak linier pada lamanya waktu komputasi ditentukan diawal serta double paths segment dalam Um dan Thurber (1987) tidak dilakukan, merupakan modifikasi metode ray tracing pseudo-bending untuk menjaga kestabilan dalam menentukan titik jejak sinar gelombang. Dari pemograman dan pengujian beberapa model kecepatan 3-D, ray tracing pseudo6 Nugraha, A.D., Syahputra, A., dan Fatkhan., 2011. Pemograman ray tracing metode pseudo-bending mediun 2-D untuk menghitung waktu tempuh antara sumber dan penerima. Jurnal Geofisika, No. 1/2. Syahputra, A., 2011. Pengembangan perangkat lunak tomografi 2-D dan 3-D: Aplikasi tomografi lubang bor dan gunungapi. Tugas Akhir. Program Studi Teknik Geofisika, ITB, Bandung. Thurber, C. H., 1993. Local earthquake tomography velocities and Vp/Vs theory, in Seismic Tomography: Theory and Practice, pp. 563-583, edited by H. M. Iyer and K. Hirahara, CRC Press, Boca Raton, Fla. Um, J.dan Thurber, C., 1987. A fast algorithm for two point seismic ray tracing. Bull. Seismol. Soc. Am., Vol.77, No.33, pp. 972-986. J. Geofisika Vol. 13 No.1/2012 Interpretasi Data Anomali Medan Magnetik Total Untuk Permodelan Struktur Bawah Permukaan Daerah Manifestasi Mud Vulcano (Studi Kasus Bledug Kuwu Grobogan) Sigit Darmawan, Hernowo Danusaputro, Tony Yulianto Jurusan Fisika Fakultas Matematika dan Ilmu Pengetahuan Alam Universitas Diponegoro, Semarang Email sigit_x@yahoo.com Abstrak Telah dilakukan transformasi reduksi ke kutub data anomali medan magnetik total pada daerah Bledug Kuwu, Grobgan untuk interpretasi struktur bawah permukaan. Data yang ditransformasi hasil pengukuran 8-10 April 2006 dengan menggunakan Proton Precession Magnetometer (PPM) untuk mengukur medan magnet total, dengan luas daerah penelitian ± 7,5 ha, yang menghasilkan 135 titik pengukuran. Penentuan posisi menggunakan Global Positioning System (GPS) dan kompas geologi. Pengolahan data medan magnetik total dimulai dari koreksi variasi harian dan koreksi IGRF sehingga diperoleh anomali medan magnetik total pada topografi. Efek anomali lokal pada lokasi penelitian dieliminasi dengan metode upward continuation. Data hasil upward continuation pada ketinggian 3000 m di atas referensi spheroid direduksi ke kutub utara magnet bumi. Anomali medan magnetik total reduksi ke kutub pada ketinggian 3000 m diinterpretasi dengan permodelan Talwani 2,5 dimensi menggunakan software Mag2DC for Window. Hasil pemodelan dua dimensi menghasilkan benda penyebab anomali dengan suseptibilitas yaitu: untuk benda pertama 0,003 cgs, dan benda kedua (-0,001) cgs, sedangkan benda di bawah anomali memiliki nilai suseptibilitas yang kecil karena temperatur yang tinggi. Benda anomali berada pada kedalaman ± (270-350) m dari permukaan dan diidentifikasi berupa garam yang bercampur shale terjebak dalam cekungan sedimentasi. Kata kunci: anomali magnetik, reduksi ke kutub, suseptibilitas. Abstract Total magnetic field anomaly data at Bledug Kuwu, Grobogan area has been reduced to the pole in order to interpret underground structure. Data collecting has been done on April 8-10, 2006 by means of Proton Precession Magnetometer (PPM) to measure total magnetic field data, with research area of 7,5 ha, produces 135 points of measurement. Global Positioning System (GPS) was used to determine the position and a geological compass to find the geographic north. The first total magnetic field data processing is diurnal and International Geomagnetic Reference Field (IGRF) correction. The total magnetic field anomaly on the irregular surface was transformed to horizontal surface 3000 m above spheroid reference. Before reduction to the pole, local effect was eliminated by upward continuation as high as 3000 m above spheroid reference. Mag2DC for Window was carried out for interpretation of total magnetic anomaly at 3000 m above reference spheroid which has been reduced to the pole. The modeling software based on the 2.5 D Talwani’s method. The result of 2-D modeling produces anomaly objects with susceptibilities: the first object: 0,003 cgs and the second object: -0,001 cgs, whereas the object under anomaly with small susceptibilities because high temperature. The anomaly objects are in the depth of ± (270-350) meter below the surface and are interpreted as salt and shale mixture in a sedimentary dome. Keyword: magnetic anomaly, reduction to pole, susceptibilities. 7 J. Geofisika Vol. 13 No.1/2012 1. Pendahuluan Bledug Kuwu merupakan lokasi wisata dengan keajaiban alam yaitu adanya fenomena gunung api lumpur atau mud volcanoes (Bemmelen, 1949). Lokasi wisata ini luasnya ± 45 hektar, terletak di desa Kuwu, kec. Kradenan, kab. Grobogan. Daerah ini mempunyai posisi geografis terletak 111007’ BT dan 07007’ LS dan terletak di dataran rendah bersuhu 28-360C. Fenomena yang dapat dilihat yaitu berupa letupan gelembung lumpur raksasa yang mengandung garam, beserta gas yang mengandung unsur belerang dan hidrokarbon. Menurut Manurung (1989), erupsi lumpur yang terjadi di daerah Bledug Kuwu terbentuk di atas zona patahan (fault zone). Suhu dan tekanan lebih besar di bagian dalam dari daerah cekungan ini menyebabkan larutan dan gas mengalir melalui rekahanrekahan pada zona patahan tersebut dan mendorong lumpur naik ke atas. Dalam penelitiannya, Manurung mengambil daerah penelitian di sekitar daerah Kuwu dengan luas (10 x10) km2 dengan jarak tiap titiknya (100-300) meter. Tujuan penelitiannya adalah menampilkan penampang bawah permukaan yang bersifat regional. Metode magnetik merupakan salah satu metode geofisika yang sering digunakan untuk survai pendahuluan pada eksplorasi minyak dan gas bumi, penyelidikan batuan mineral dan penyelidikan tentang panas bumi. Di Jawa, telah banyak dilakukan penelitian dengan metode ini diantaranya: dalam penyelidikan panas bumi misalnya di Gunung Ungaran (Haryono, 2002), (Nurdiyanto, 2004), di Gunung Tangkuban Perahu (Yulianto, 2000) dan untuk pemodelan sesar regional Gunung Merapi-Merbabu (Ismail, 2001). Metode ini mempunyai akurasi pengukuran yang relatif tinggi, pengoperasian di lapangan relatif sederhana, mudah dan cepat jika dibandingkan dengan metode 8 geofisika lainnya. Metode magnetik bekerja berdasarkan sifat-sifat megnetik batuan yang terdapat di bawah permukaan bumi. Diharapkan dari hasil interpretasi akan diketahui struktur bawah permukaan di daerah Bledug Kuwu. Dari hasil tersebut dapat digunakan untuk menentukan penyebaran daerah yang masih berpotensi terjadi letupan lumpur sehingga dapat digunakan untuk pengembangan fasilitas lokasi wisata Bledug Kuwu. 2. Teori 2.1 Kontinuasi ke Atas Konsep dasar pengangkatan ke atas berasal dari identifikasi tiga teorema Green. Teorema ini menjelaskan bahwa apabila suatu fungsi U adalah harmonik, kontinu dan mempunyai turunan yang kontinu di sepanjang daerah R, maka nilai U pada suatu titik P di dalam daerah R dapat dinyatakan (Blakely, 1995): U(P) = ∂1 1 1 ∂U dS −U 4π S r ∂n ∂n r (1) dengan S menunjukkan permukaan daerah R, n menunjukkan arah normal keluar dan r adalah jarak dari titik P ke suatu titik pada permukaan S. Persamaan (1) menggambarkan secara dasar prinsip dari pengangkatan ke atas, dimana suatu medan potensial dapat dihitung pada setiap titik di dalam suatu daerah berdasarkan sifat medan pada permukaan yang melingkupi daerah tersebut. 2.3 Reduksi ke Kutub Baranov dan Naudy (1964) telah menggambarkan metode transformasi ke kutub untuk menyederhanakan interpretasi data magnetik pada daerah-daerah berlintang rendah dan menengah. Metode reduksi ke kutub magnetik bumi dapat mengurangi salah tahap yang rumit dari proses interpretasi, dengan anomali medan magnetik menunjukkan langsung posisi bendanya. J. Geofisika Vol. 13 No.1/2012 n̂ S Region R Q( x' , y' , z ' ) r . z0 P ( x , y , z 0 − ∆z ) ρ x x P ' ( x , y , z 0 + ∆z ) z y n̂ Sumber Gambar 1 Pengangkatan ke atas dari permukaan horizontal (Blakely, 1995) Medan Medan Magnetik Observasi Reduksi Ke kutub Pseudo Gravity Medan Magnet Jarak Magnetisasi Depth Gambar 2. Hubungan antara medan magnet observasi, reduksi ke kutub dan pseudogravity (Tchernychev, 2001) Formulasi yang umum sebagai hubungan antara medan potensial (f ) dengan distribusi material sumber (s): f (P ) = s(Q )Ψ (P, Q )dv (2) R Fungsi f (P ) adalah medan potensial atau anomali total medan magnetik pada P, sedangkan s(Q) kuantitas fisis magnetisasi pada Q dan Ψ (P, Q ) suatu fungsi Green berupa anomali total medan magnetik dipole tunggal yang bergantung pada geometris tempat titik observasi P dan titik distribusi sumber Q. Proses transformasi reduksi ke kutub dilakukan dengan mengubah arah magnetisasi dan medan utama dalam arah vertikal. Reduksi ke kutub dilakukan dengan cara membuat sudut inklinasi benda menjadi 900 dan deklinasinya 00. Karena pada kutub magnetik arah dari medan magnet bumi ke bawah dan arah dari induksi magnetisasinya ke bawah juga. 2.4 Geologi Daerah Penelitian Keadaan geologi regional menunjukkan bahwa mulai dari Semarang kearah timur hingga daerah Kuwu merupakan endapan alluvial yang termasuk zona Randublatung (Cipluk beds serta Lower Kalibeng beds). Daerah penelitian mempunyai kenampakan morfologi datar. Di bagian Utara terdapat perbukitan bergelombang lemah dan sedang. 9 J. Geofisika Vol. 13 No.1/2012 Sedangkan di bagian Selatan dibatasi oleh bagian darat formasi Kendeng Ridge (Bemmelen, 1949). Di sebelah timur daerah penelitian terdapat jalur patahan yang berarah Barat-Timur, yang merupakan patahan normal. Juga di sebelah selatan terdapat jalur patahan yang berarah barat-timur yang merupakan patahan naik. Tegak lurus patahan tersebut terdapat patahan normal. Geologi secara jelas dapat dilihat pada gambar 3. A B 45 ’ 111 0 15 ’ 30 ’ (BT) 70 15 ’ Alat yang dipergunakan meliputi: satu buah Proton Precession magnetometer (PPM) model G-856 Geometrics untuk merekam waktu dan medan magnet total (dalam satuan nT), satu buah Global Positioning System (GPS) model Trimble 4600TM LS frekuensi tunggal untuk menentukan posisi penelitian dengan ketelitian 0,1 dan sebuah kompas geologi untuk menentukan arah geografi lokasi pengukuran serta alat komunikasi. Proton Precession magnetometer (PPM) yang digunakan hanya satu maka pengambilan data dilakukan dengan cara Loopping, artinya setelah melakukan pengambilan data pada titik-titik pengukuran medan magnet yang telah ditentukan, maka harus kembali ke base untuk mengukur medan magnetnya lagi. Setelah itu pengukuran dilanjutkan pada titiktitik berikutnya dan kembali ke base lagi. Selang waktu pengukuran antar base harus kurang dari satu jam atau waktunya singkat agar variasi hariannya masih terpantau dengan baik. Setiap titik pengukuran diambil lima kali data yang berbeda dalam jarak ± 1 meter, diambil nilai terbaik atau nilai rata-rata. Posisi titik pengukuran dilakukan menggunakan Global Positioning System (GPS) yaitu dalam satuan derajat untuk lintang dan bujurnya. 3.2 Perhitungan Anomali Medan Magnet 30 ’ (LS) Gambar 3. A. Lokasi daerah penelitian, B. Peta Geologi Daerah Penelitian dan Sekitarnya (Direktorat Geologi, 1963) 3. Metode Penelitian 3.1 Pengambilan Data Lokasi penelitian mencakup daerah kawasan wisata Mud Vulcano Bledug Kuwu, yang terletak di desa Kuwu, kec. Kradenan, kab. Grobogan yang secara geografis terletak 111007’ BT dan 07007’ LS, dengan luas daerah penelitian ± (300 x 250) meter atau ± 7.5 ha. Pengambilan data dimulai tanggal 8 sampai 10 April 2006 sebanyak ± 135 titik yang terletak di sekitar letupan lumpur Bledug Kuwu. Pada penelitian ini, setiap titik pengukuran berjarak ± (20-25) meter. 10 Di dalam survai dengan metode magnetik digunakan satu set magnetometer yang pengambilannya dilakukan dengan cara loopping. Maka dalam survai, setelah pengukuran di tiap titik-titik pengukuran harus kembali ke base (dalam beberapa menit atau kurang dari satu jam). Pengukuran base ini diulang-ulang terus untuk mendapatkan variasi harian yang diakibatkan efek medan magnet luar bumi, dan untuk mengoreksi titiktitik pengukuran. Sedangkan medan magnet bumi dihitung berdasarkan pada persamaan International Geomagnetic Reference Field (IGRF), sehingga anomali magnetiknya diberikan oleh persamaan sebagai berikut: ∆T = Tobs − TIGRF ± Tvh dengan Tobs adalah medan (3) magnetik komponen total terukur, TIGRF adalah medan magnet teoritis berdasarkan IGRF pada J. Geofisika Vol. 13 No.1/2012 stasion dan Tvh adalah koreksi medan magnet akibat variasi harian. TIGRF dihitung pada titik pengukuran dengan memasukkan nilai posisi dan tanggal pengukuran dengan paket program IGRF yang telah terdapat pada beberapa software misalnya Magpick, WMM, dll. Sedangkan Tobs terukur pada saat magnetometer merekam data pada titik pengukuran. Hasil pengolahannya dengan Microsoft office Excel didapatkan data anomali medan magnet total. 3.3 Peta Anomali Medan Magnet Berdasarkan hasil pengolahan data yang diperoleh, dibuat peta anomali medan magnet dengan menggunakan software paket surfer version 8 yang menunjukkan hubungan antara posisi pengukuran dan nilai anomali medan magnet total ditunjukkan dalam gambar 5A. Penghalusan data pengamatan untuk mengeliminasi efek lokal dilakukan dengan kontinuasi ke atas (upward continuation) sebesar 3000m. Hasil kontinuasi ini terlihat pada gambar 5B yang memperlihatkan anomali yang muncul semakin jelas. Setelah dilakukan kontinuasi ke atas, data anomali medan magnetik total ini direduksi ke kutub. Kedua tahap ini dilakukan dengan menggunakan program MagPick atau reduksi ke kutubnya dengan sofware Signpro. Gambar 4. Diagram blok pengolahan data magnetik total 11 J. Geofisika Vol. 13 No.1/2012 mendorong keluar. Dengan adanya aktivitas ini, maka batuan akan mengalami penurunan sifat kemagnetannya, sesuai dengan aktivitas bledug yang mengeluarkan erupsi lumpur yang mengandung garam dan gas belerang serta gas metan lainnya. Gambar 5. A. Anomali Medan magnet total, B. Anomali medan magnet total (upward continuation 3000m). 4. Hasil dan Diskusi 4.1 Interpretasi Kualitatif Secara kualitatif peta anomali diperoleh menunjukkan penyebaran pasangan pola kontur tertutup (besar-kecil) yang terdapat pada masing-masing bledug. Penentuan pasangan ini didasarkan pada kecenderungan arah grid setiap pasangan kontur tertutup. Oleh karena itu, dapat terlihat anomali berarah utara-selatan untuk bledug pertama dan anomali berarah barat-timur untuk bledug ke dua, dengan pusat benda anomali ditafsirkan berada di tengah pasangan pola kontur tertutup itu. Dari pola-pola anomali yang terlihat mempunyai gradien anomali horisontal yang tinggi (gradiennya tajam) dari pada daerah sekitarnya. Di daerah dekat pusat bledug terlihat adanya anomali yang menunjukkan bahwa pada daerah inilah yang mengakibatkan terjadinya letupan lumpur. Ditafsirkan adanya aktifitas panas dari dalam yang berupa gas yang 12 Gambar 6. A. Sayatan anomali medan magnet total, B. Pemodelan pada sayatan profil A-A’, C. Pemodelan pada sayatan profil B-B’ Sayatan pertama dibuat dari pasangan kontur tertutup yang berarah utara-selatan yaitu A-A’ yang melewati bledug pertama (sebelah timur). Sayatan ke dua dibuat dari pasangan kontur tertutup berarah barat timur yaitu B-B’ yang melewati daerah bledug ke dua (sebelah barat). Dari kedua sayatan ini, akan digunakan untuk permodelan struktur bawah permukaan daerah Bledug Kuwu. Sayatan ini diambil dari data peta anomali yang telah dilakukan upward continuation 3000m. Dengan metode J. Geofisika Vol. 13 No.1/2012 ini, akan mempertajam anomali pasangan kontur dari data peta anomali medan magnet total. Hasil sayatan ini (upward continuation), kemudian dilakukan reduksi ke kutub untuk mengubah arah magnetisasi benda dalam arah vertikal sehingga anomali medan magnetik dapat menunjukkan langsung posisi benda penyebabnya. 4.2 Interpretasi Kuantitatif Pada pemodelan profil A-A’ dan B-B’, bentuk kurvanya hampir sama Hasil pemodelan profil ini didapatkan 2 benda penyebab anomali dengan nilai -0,012cgs untuk benda pertama (warna merah) dan (-0,016)cgs untuk benda ke dua (warna biru), dengan arah barat-timur. Benda pertama berada di sebelah barat dengan kedalaman (270-330) meter dari permukaan dan benda ke dua dengan kedalaman (270350) meter dengan sisi yang berbatasan dengan benda pertama (dekat batuan pertama) mempunyai lapisan lebih tipis. Di bawah kedua batuan anomali ini, terdapat batuan sedimen yang sangat dipengaruhi panas (warna hitam) dengan kontras suseptibilitas (0,014)cgs sebagai sumber tekanan. Sumber tekanan (gas) ini mencari tempat pada benda anomali yang lemah untuk dilaluinya sampai ke permukaan bumi dan daerah ini dinamakan zona lemah. Dari harga nilai suseptibilitas batuan sekitar (k0), kontras suseptibilitas batuan ( k) dan suseptibilitas batuan target (k1), maka dengan persamaan 4 berikut dapat dicari benda penyebab anomalinya (berdasar nilai suseptibilitas) sebagai ∆k = k1 − k0 (4) Tabel 2. Hubungan suseptibilitas batuan sekitar (k0), kontras suseptibilitas ( k) dan suseptibilitas batuan target (k1) Pemodelan Nilai k0 (cgs) Nilai k Nilai k1 (cgs) Kemagnetan Batuan 1 A-A’ 0,015 -0,012 0,003 Paramagnet Batuan 2 A-A’ 0,015 -0,016 -0,001 Diamagnet Batuan 1 B-B’ 0,015 -0,012 0,003 Paramagnet Batuan 2 B-B’ 0,015 -0,016 -0,001 Diamagnet Batuan 3 0,015 -0,014 0,001 Paramagnet Dari pengolahan data dan pemodelan perhitungan nilai suseptibilas ini, dapat dilihat bahwa batuan 1 A-A’ dan 1 B-B’ merupakan batuan dengan suseptibilitas kecil dan positif 0,003cgs merupakan batuan paramagnet, kemudian batuan 2 A-A’ dan 2 B-B’ adalah batuan dengan suseptibilitas kecil dan negatif (-0,001)cgs merupakan batuan diamagnet. Batuan yang suseptibilitasnya negatif ini diidentifikasi sebagai batuan garam (rocksalt) dapat berupa padat, lumpur maupun cairan. Dari tabel ini menunjukkan batuan pertama dan ke dua untuk kedua sayatan adalah batuan yang sama dengan nilai suseptibilitas (cgs) 0,003cgs dan (-0,001) cgs. Sedangkan batuan 3 di bawahnya merupakan batuan yang sangat dipengaruhi suhu dan tekanan sehingga suseptibilitasnya kecil (-0,001) cgs. Tekanan dan suhu yang tinggi menyebabkan batuan yang dilaluinya menjadi kehilangan sifat kemagnetannya. Nilai kontras suseptibilitas batuan yang cenderung lebih kecil dari batuan sekitarya menunjukkan aktifitas panas telah banyak mempengaruhi batuan tersebut. Dengan kata lain, batuan dengan kontras suseptibilitas lebih kecil (lebih negatif) menunjukkan tekanan dari bawah 13 J. Geofisika Vol. 13 No.1/2012 akibat aktifitas panas lebih besar dilakukan kepadanya dari pada batuan sekitarnya. Sehingga batuan ini akan menghasilkan letupan yang lebih besar dan periodenya cepat. Hal ini sesuai dengan fenomena di daerah penelitian, bahwa bledug pertama lebih aktif menghasilkan letupan dan periode letupannya lebih cepat dibandingkan bledug ke dua. Terjadinya letupan dikarenakan adanya tekanan dari bawah mampu mendorong batuan yang dilaluinya terangkat naik. Oleh karena itu, batuan ini harus bersifat lemah terhadap tekanan atau mudah dilalui gas (sumber tekanan). Selanjutnya, harus ada pula sumber tekanan dari bawah yang besar dan keluar melewati batuan ini. Pada prinsipnya benda di dalam bumi akan keluar ke permukaan karena di dalam bumi suhu dan tekanannya besar. Bila batuan dasarnya sangat keras maka benda dengan tekanan besar ini seperti terperangkap dan tidak bisa keluar. Benda di dalam bumi ini dapat keluar jika terdapat rekahan, patahan, ataupun karena adanya aktifitas pemboran. Sehingga syarat terjadinya letupan pada daerah Bledug Kuwu harus ada patahan yang terjadi di bawah batuan hasil pemodelan yang telah disebutkan di atas. Hasil interpretasi, kemudian dibuat struktur bawah permukaan mengenai terjadinya letupan di daerah penelitan yang ditunjukkan gambar 7. Terjadinya letupan hanya terjadi di atas batuan yang suseptibilitasnya kecil dan negatif sebagai deretan bledug dari besar sampai kecil. Tekanan yang melalui batas perlapisan menyebabkan tekanan memusat pada batas kontak batuan pertama dan ke dua menghasilkan tekanan paling besar. 5. Kesimpulan Dan Saran 5.1 Kesimpulan 1. Struktur bawah permukaan di daerah Bledug Kuwu terdiri dari: • Batuan penyebab anomali, ada dua jenis yaitu dengan suseptibilitas 0,003cgs, dan suseptibilitas 0,001cgs. • Batuan di atas anomali (batuan sekitar) adalah shale. • Batuan yang berada di bawah anomali berkurang sifat kemagnetannya yaitu dengan suseptibilitas 0,001 cgs. 2. Kedalaman benda anomali rata-rata adalah (270-350) meter. 3. Daerah di atas batuan penyebab anomali dengan suseptibilitas (-0,001)cgs, merupakan daerah potensial terjadi letupan. 4. Dari interpretasi menunjukkan bahwa batuan daerah penelitian adalah sedimen yaitu shale yang telah berkurang sifat kemagnetannya dan mengandung salt, water sebagai anomali. 5.2 Saran 1. Penelitian harus ditambah lagi atau diperluas di daerah di sekitar Bledug Kuwu dengan jarak tiap titik-titiknya pendek. 2. Dapat dilakukan mengolahan dan interpretsi magnetik dengan cara atau metode yang berbeda. 3. Hasil penelitian magnetik ini harus dicocokkan dengan data lubang bor, data seismik, dan data lainnya. 4. Perlu dilakukan penelitian lebih lanjut mengenai pengaruh perubahan temperatur terhadap nilai suseptibilitas batuan. Daftar Pustaka Baranov, V. and Naudy, H., 1964, Numeric Calculation of the Formula of reduction to pole, Geophysics, 29, 6769. Bemmelen, R.W., van, 1949, The Geology of Indonesia, V.I.A, Martinus Nijhoff, The hague. Gambar 7. Pemodelan dan interpretasi struktur bawah permukaan Bledug Kuwu 14 Bhaskara, R.D., Ramesh, B.N., 1991, A Rapid Method for Three-dimentional Modelling of Magnetik anomalies: Geophysics. 56,1729-1737. J. Geofisika Vol. 13 No.1/2012 Blakely, R.J., 1995, Potential theory in gravity and magnetic applications, Cambridge Univ Press, New York. Parkinson,W.D.,1983, Geomagnetism, Press London. Breiner, S., 1973, Applications Manual for Portable Magnetometers, Geometrics, USA. Reid, A.B., Allsop, J.M., Granser, H., Millet, A.J., and Somerton, I.W., 1990, Magnetic interpretation in three dimensions using Euler deconvolution: Geophysics, 55, 80– 91. Grant, F.S.,West, 1965, Interpretation Theory in Applied Geophysics, McGraw Hill Corporation. Haryono, A., 2002, Pemodelan Sesar Regional di Daerah Gunungapi Ungaran menggunakan Data Anomali Medan Magnetik Reduksi ke Kutub, Tesis, FMIPA, UGM. IAGA Working Group V-8, 1995, International Geomagnetic Reference Field, 1995 revision. Submitted to EOS Trans. Am. Geophys. Un., Geophysics, Geophys. J. Int., J. Geomag. Geoelectr.,Phys. Earth Planet.Int., and others. Ismail, N.,2001, Interpretasi Data Anomali medan Magnetik Total Reduksi ke Kutub Untuk Pemodelan sesar Regional di Daerah Gunung MerapiMerbabu, Tesis, FMIPA, UGM. Manurung, P., 1989, Penyelidikan Anomali Medan Magnet Total di Daerah Kuwu, Grobodan, Jawa Tengah, Skripsi UGM. McLean, S., S. Macmillan, S. Maus, V. Lesur, A.Thomson, and D. Dater, 2004, TheUS/UK World Magnetic Model for 2005-2010, NOAA Technical Report NESDIS/NGDC-1. Nurdiyanto, B., 2004, Analisis Data Anomali Medan Magnet Total Untuk Menafsirkan Struktur Bawah Permukaan Daerah Manifestasi Air Panas di Lereng Utara Gunungapi Ungaran, Skripsi, FMIPA, UGM. Introduction to Scottish academic Robinson, E. S., Coruh, C., 1998, Basic Exploration Geophysics, John Willey & Sons. Sharma, P.V., 1997, Environmental and Engineering Geophysics, Cambridge University Press. Sulindra , N., 2005, Interpretasi Data Anomali medan Magnet Total reduksi ke Kutub untuk Pemodelan Sesar, Skripsi, FMIPA, UGM. Talwani, M. and Heirtzler, J.R.,1964, Computation of Magnetic anomalies Caused by two Dimentional Structures of Arbitary Shapein The Mineral Industries, Stanford University Publications Geological Sciences Vol. 9, No.1. Tchernychev, M.,2001, Magpick-magnetic map & profile processing, user guide. Telford, W.M., Geldart, R.E., Sheriff, D.A., and Keys, 1979, Applied Geophysics, Cambridge University Press. University of Birmingham, 2004, Classification of magnetic material, Applied Alloy Chemistry Group. U.S. Geological Survey Information Service, 2005, World IGRF Magnetic Chart, web page:www.ngdc.noaa.gov. Yulianto, T., 2000, Pengukuran dan interpretasi anomali magnetik daerah Gunung Tangkuban Perahu, Tesis Pasca Sarjana ITB. 6 2 15 J. Geofisika Vol. 13 No.1/2012 Penentuan Hiposenter Menggunakan Simulated Annealing Dan Guided Error Search Serta Penentuan Model Kecepatan Gelombang Seismik 1-D Pada Lapangan “Geothermal” Akhmad Fanani Akbar1, Andri Dian Nugraha1, M. Rachmat Sule1, Aditya Abdurrahman Juanda2 1) Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan, Institut Teknologi Bandung, Jalan Ganesha 10 Bandung, 40132 2) PT. Pertamina Geothermal Energy, Indonesia Email : fanani_akbar@yahoo.co.uk Abstrak Penentuan hiposenter gempa mikro untuk lapangan geotermal gunung “Tinggi” telah dilakukan dengan menggunakan metoda simulated annealing dan guided error search dan model kecepatan gelombang seismik satu dimensi (1D). Untuk mempercepat proses penentuan lokasi hiposenter metoda perpotongan tiga lingkaran digunakan untuk memfokuskan area dimana penentuan lokasi hiposenter dilakukan. Data yang digunakan adalah waktu tiba gelombang P dan S. Dalam proses simulated annealing dan guided error search, penghitungan waktu tempuh minimum dari source ke receiver dilakukan dengan menggunakan ray tracing dengan metoda shooting. Hasil dari simulated annealing dan guided error search menunjukkan bahwa gempa terjadi pada kedalaman 3-4 km di bawah permukaan laut. Hal ini sesuai dengan studi sebelumnya yang menyatakan bahwa pada kedalaman tersebut terdapat area paling aktif dimana tempat terjadinya banyak fracture. Hasil posisi hiposenter gempa tersebut digunakan sebagai salah satu data awal pada program VELEST yang berfungsi untuk melakukan perbaikan model kecepatan gelombang seismik 1D. Hasil dari VELEST menunjukkan bahwa terdapat Vp/Vs yang rendah pada kedalaman 3-4km. Kemungkinan berhubungan dengan batuan yang tersaturasi uap (gas). Kata kunci: Hiposenter, Simulated Annealing, Guided Grid Search, Model Kecepatan 1-D. Abstract In this study, hypocenter determination of micro-earthquakes of Mount “Tinggi” has been conducted by employing simulated annealing and guided error search method using a 1D velocity model. In order to speed up the hypocenter determination process a three-circle intersection method has been used to guide the simulated annealing and guided error search process. We have used P and S arrival time data. In the simulated annealing and guided error search process, the minimum travel time from a source to a receiver has been calculated by using ray tracing with shooting method. The results show hypocenters of microearthquakes occurred at depths of about 3-4 km below mean sea level. The location of microearthquake in this study are correlated with high fractured zone that were inferred from previous study. We then determine 1-D velocity model by applying VELEST method. The results of VELEST indicate there are low Vp/Vs ratio values at depths of 3-4km. We can interpret this feature as a rock layer which is saturated by vapor (gas). Keywords: Hypocenter, Simulated Annealing, Guidded Grid Search, 1-D Velocity Model 16 J. Geofisika Vol. 13 No.1/2012 1. Pendahuluan Aktivitas gunung api dan pergerakan lempeng dapat diketahui dengan pemetaan lokasi gempa. Selain itu pemetaan lokasi gempa juga dapat digunakan untuk memantau dan menganalisis reservoir geotermal khususnya gempa lokal dan gempa mikro. Namun penentuan lokasi absolut gempa tersebut dipengaruhi oleh faktor-faktor sebagai berikut yaitu geometri stasiun pengamat, akurasi pembacaaan waktu tiba, fasa gelombang yang tersedia, dan pengetahuan tentang struktur geologi pada daerah studi (Gomberg dkk., 1990). Untuk skala lokal, faktor geometri dan ketersediaan data bisa diperbaiki dengan menambah jumlah stasiun pengamat pada daerah penelitian. Faktor kesalahan model kecepatan gelombang seismik dan penentuan waktu tiba dapat diminimalisir dengan analisis yang lebih lanjut. Model kecepatan gelombang seismik bawah permukaaan tidak bisa ditentukan secara pasti karena keterbatasan data dan kompleksitas struktur bawah permukaan. Oleh karena itu, diperlukan model sederhana bawah permukaan untuk menentukan posisi gempa dengan baik. Untuk menentukan posisi gempa dengan baik diperlukan adanya proses pemodelan ke belakang. Metode pemodelan ke belakang ini bertujuan untuk mencari posisi yang memiliki nilai selisih antara data observasi dan data lapangan yang paling kecil (minimum global). Teknik pemodelan ke belakang ini pada dasarnya adalah teknik pemodelan dengan cara mencoba-coba dan memodifikasi parameter model sehingga diperoleh kecocokan antara data perhitungan atau data estimasi dengan data lapangan. Namun karena heterogenitas dari kondisi geologi bawah permukaan bumi dan beberapa alasan lainnya seperti kekurangan data dalam membatasi atau mendefiniskan solusi, adanya pengaruh bising pada data lapangan, maka solusi dari pemodelan kebelakang ini menjadi tidak unik, artinya satu model bawah permukaan yang berbeda dapat memberikan respon dipermukaan yang sama. Oleh karena itu dibutuhkan suatu mekanisme pemodelan ke belakang yang dapat menghasilkan banyak realisasi atau kemungkinan model yang sesuai. Salah satu cara penyelesaiannya yaitu dengan menggunakan metode Simulated Annealing. Pada kenyataannya pencarian solusi optimum dengan metode Simulated Annealing masih menimbulkan masalah, yaitu karena pada aturan ini semua kemungkinan yang ada disertakan dalam pencarian solusi yang menyebabkan proses perhitungan menjadi sangat lama dan kisaran nilai dari hasil pencarian solusi yang didapat relatif masih terlalu luas atau belum mendekati unik. Oleh karena itu pada penelitian ini dilakukan modifikasi pada proses perturbasi model yaitu dengan menggunakan metode Guided Error Search. Latar belakang digunakan metode ini adalah agar batas ruang eksplorasi model semakin lama semakin dipersempit sehingga proses iterasi lebih cepat dan memiliki error yang minimum. 2. Metode Untuk menentukan lokasi suatu gempa dibutuhkan waktu terjadinya gempa yang nantinya akan digunakan untuk menentukan travel time. Pada studi ini, waktu tiba gelombang P dan S untuk gempa mikro di-pick secara manual (Gambar 1). Pada penelitian kali ini digunakan diagram Wadati untuk menentukan waktu terjadinya gempa mikro (Gambar 2). Beda waktu tiba gelombang P dan S (tS - tP) diplot terhadap waktu tiba gelombang P. Karena di hiposenter nilai tS - tP bernilai nol, maka titik potong garis lurus hasil dari regresi merupakan pendekatan waktu terjadinya gempa. 17 J. Geofisika Vol. 13 No.1/2012 Search dibutuhkan lokasi gempa awal menggunakan metode perpotongan tiga lingkaran. Persamaan umum lingkaran adalah sebagai berikut: (6) Gambar 1. Contoh waveform gempa mikro yang terekam oleh empat stasiun pada komponen vertikal dalam studi ini. Dimana nilai a dan b merupakan pusat lingkaran dengan jari-jari r. Waktu tiba gelombang P dan S diperlukan untuk menentukan jari-jari lingkaran tersebut. Setelah didapatkan nilai episenter dari gempa, kita bisa menentukan posisi kedalaman gempa terhadap salah satu stasiun. 2.1. Simulated Annealing Gambar 2. Contoh diagram Wadati dalam penentuan waktu terjadinya gempa mikro dalam studi ini. Dari Gambar di atas didapatkan persamaan sebagai berikut dan (1) Dimana merupakan waktu terjadinya gempa dan dengan Vp > Vs dan ts > tp maka diperoleh persamaan: (2) (3) (4) Dari persamaan 3 dan 4 didapat persamaan sebagai berikut: (5) Untuk perhitungan mempercepat proses metode Simulated Annealing dan Guided Error 18 Pencarian solusi secara acak kurang efisien mengingat jumlah model yang harus dievaluasi masih cukup besar. Pemilihan model secara acak murni menyebabkan semua model dalam ruang model memiliki probabilitas yang sama untuk dipilih sebagai sampel dalam perhitungan fungsi obyektif. Model pada daerah yang jauh dari solusi juga harus dihitung responsnya. Untuk meningkatkan efisiensi pencarian acak, pemilihan model dimodifikasi sehingga model pada daerah tertentu yang mengarah atau dekat dengan solusi memiliki probabilitas lebih besar untuk dipilih. Metode semacam ini disebut sebagai metode guided random search. 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 energy sistem yang minimum (Gambar 3). Pada temperatur tinggi sistem mengalami perturbasi dan perturbasi yang menghasilkan konfigurasi dengan energi rendah memiliki kemungkinan lebih besar. Pada saat temperatur menurun, perturbasi yang menghasilkan J. Geofisika Vol. 13 No.1/2012 konfigurasi dengan energi lebih rendah memiliki probabilitas makin besar, sedangkan perturbasi yang menghasilkan konfigurasi dengan energi lebih tinggi probabilitasnya makin kecil. 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 masingmasing. Caranya adalah mengambil bilangan acak R dengan probabilitas uniform antar 0 dan 1 yang dipetakan menjadi harga parameter model menggunakan persamaan berikut: Modeli=Modelimin+R(Modelimax-Modelimin) (7) Proses pembentukan kristal (annealing) dalam termodinamika diadopsi dalam penyelesaian masalah inversi dengan menggunakan parameter model untuk mendefinisikan konfigurasi sistem dan fungsi objektif (misfit) sebagai energi. Faktor temperatur merupakan faktor pengontrol dengan satuan sama dengan fungsi objektif. Probabilitas perturbasi model dinyatakan oleh: P( E) =exp (8) Dimana E adalah perubahan fungsi objektif atau perubahan misfit akibat 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. model diterima. Pada saat temperatur turun secara perlahan perturbasi model diterima mengecil jika E 0. Hal ini bertujuan untuk menghindari konvergensi solusi keminimum lokal. Gambar 3. Algoritma Simulated Annealing sederhana untuk inversi non-linier (Grandis, 2009) (7) Mekanisme penurunan temperatur merupakan faktor yang perlu dilakukan secara coba-coba agar sesuai dengan permasalahan yang ditinjau. Proses SA berlangsung sebanyak jumlah iterasi tertentu dengan penurunan temperatur secara bertahap, pada penelitian ini skema penurunan temperatur dirumuskan: Tn = T n-1 (9) Tn adalah temperatur ke-n dan adalah konstanta penurunan temperatur. Penurunan temperatur tidak boleh terlalu cepat ataupun terlalu lambat. Penurunan yang terlalu cepat akan menyebabkan proses terjebak pada minimum lokal,sedangkan penurunan yang terlalu lambat akan membutuhkan waktu yang lama menuju fungsi objektif minimum. Dengan pemilihan skema temperatur yang tepat maka diharapkan akan memperoleh solusi minimum global. Skema penurunan temperatur dapat menggunakan fungsi linier ataupun logaritmik. Selain mempengaruhi waktu pengolahan, fungsi tersebut juga sangat mempengaruhi tingkat akurasi hasil akhir. Proses iterasi dimulai dengan temperature cukup tinggi sehingga hampir semua perturbasi 19 J. Geofisika Vol. 13 No.1/2012 2.2. Guided Error Search Pencarian menggunakan guided error search (Gambar 4) merupakan salah satu pencarian yang pada dasarnya memperkecil ruang eksplorasi model. Berbeda dengan simulated annealing yang menggunakan ruang eksplorasi model yang tetap dari awal hingga akhir sehingga proses untuk mencari nilai minimum lokal sangat lama. merupakan perubahan fungsi objektif akibat perubahan model. E = E(Mn+1) – E(Mn) (10) Jika E < 0 maka perubahan model diterima, dan nilai error dari model tersenut dikalikan suatu konstanta (dalam penelitian ini digunakan kecepatan gelombang seismik) yang nantinya hasil perkalian tersebut digunakan untuk jarak antara model dengan batas eksplorasi ruang model. Setiap kali iterasi berlangsung, algoritma akan mencari satu titik referensi yang nantinya dibutuhkan untuk melakukan perhitungan untuk iterasi berikutnya. 2.3. Penentuan Model Kecepatan 1-D Gambar 4. Algoritma guided error search Sama halnya dengan simulated annealing, Pada proses guided error search 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 menggunakan persamaan 7. Pada model awal (Mn) dan model yang dipilih secara acak (Mn+1) tersebut ditentukan nilai error RMSnya, kemuudian dicari selisih diantara keduanya E dimana nilai tersebut 20 Pada penelitian ini model kecepatan gelombang seismik bawah permukaan yang digunakan adalah model 1D, hal ini dilakukan sebagai bentuk penyederhanaan masalah terhadap model bawah permukaan sebenarnya. Model kecepatan gelombang seismik 1D biasa digunakan sebagai prosedur dalam penentuan lokasi gempa dan sebagai inisial model untuk seismik tomografi (Kissling, 1995; Kissling et al., 1994). Salah satu metode penentuan model kecepatan gelombang seismik 1-D, adalah metode coupled velocity-hypocenter menggunakan program VELEST versi 3.1 (10.4.95) (Kissling, 1995). Metode coupled velocity-hypocenter merupakan metode relokasi gempa, penentuan model kecepatan gelombang seismik bawah permukaan 1-D, dan koreksi stasiun secara bersamaan menggunakan prinsip metode Geiger. Jumlah parameter modelnya (m) adalah 5 + N, (x, y, z, t0, koreksi stasiun, dan N adalah jumlah lapisan model kecepatan gelombang seismik 1D). Sebagai tahap pertama, didefinisikan m0 yaitu parameter hiposenter (x, y, z, t0), model kecepatan gelombang seismik (1D) dan koreksi stasiun. Selanjutnya forward modeling dilakukan dengan ray tracing dari gempa ke J. Geofisika Vol. 13 No.1/2012 stasiun sehingga memperoleh Tcal (Waktu tempuh perhitungan). Inverse modeling dilakukan dengan menyelesaikan Matriks Damped Least Square [At A + L] (A = Matriks Jacobi, At = Tranpos Matriks Jacobi; L = Matriks damping). Penggunaan nilai damping akan mempengaruhi nilai perturbasi parameter model ( m), dengan hubungan antara besarnya damping dan nilai m adalah berkebalikan. Hasil dari inverse modeling adalah vektor perbaikan parameter model ( m) yang selanjutnya diperoleh nilai parameter hiposenter, model kecepatan gelombang seismik 1D, dan koreksi stasiun. Dalam tahapan berikutnya nilai-nilai tersebut di forward modeling untuk memperoleh nilai Tcal baru yang akan dibandingkan misfitnya dengan Tcal sebelumnya dan demikianlah tahapan dalam VELEST untuk satu iterasi. Dalam setiap iterasinya, tercantum nilai RMS antara data waktu tempuh observasi dan waktu tempuh perhitungan, sehingga jumlah iterasi dapat diatur hingga memenuhi kriteria RMS yang diharapkan. 3. Hasil dan Diskusi Gambar 5. Persebaran hiposenter gempa menggunakan simulated annealing, guided error search, dan tiga lingkaran. Setelah didapatkan nilai hiposenter, kami melakukan penentuan model kecepatan 1-D yang salah satu input-nya merupakan lokasi hiposenter. Seperti yang terdapat pada Gambar 6, bahwa pada kedalaman 3-4 km terjadi penurunan nilai rasio Vp/Vs hasil dari studi ini. Dimana menurut Boitnott dkk (1997), hal ini mengindikasikan adanya vapour reservoir. Turunnya nilai rasio Vp/Vs ini, kemungkinan dikarenakan nilai Vp yang rendah karena meningkatnya compressibility dan meningkatnya Vs karena terjadi reduksi pada pore pressure yang mengakibatkan naiknya nilai shear modulus. Pada hasil penentuan posisi gempa dengan menggunakan simulated annealing dan guided error search, persebaran hiposenter berada pada daerah gunung “Tinggi” dan sesar pada kedalaman sekitar 3 km sampai dengan 4.4 km di bawah permukaan laut (Gambar 5). Hal ini sesuai dengan apa yang diinterpretasikan oleh ELC-Electroconsukt (1983), dimana pada kedalaman tersebut merupakan area paling aktif dimana tempat terjadinya banyak fracture yang juga merupakan tempat sirkulasi vertical (Gambar 6). Gambar 6. Sirkuasi hydrothermal gunung “Tinggi” (ELC-Electroconsult, 1983) beserta kurva Vp, Vs, dan Vp/Vs yang diperoleh dari studi ini. 21 J. Geofisika Vol. 13 No.1/2012 Tabel 1. Stasiun koreksi gelombang P dan S Stasiun ST1 ST2 ST3 ST4 ST5 ST6 Gel. P 0.02 0 0.06 0 0.02 0.06 Gel. S 0 0.07 0.15 -0.01 0.01 0.05 Hasil output dari software VELEST tidak hanya berupa nilai kecepatan model 1D saja, namun juga terdapat koreksi stasiun yang mencerminkan kondisi batuan di daerah sekitar stasiun (Tabel 1). Stasiun yang memiliki nilai koreksi stasiun yang kecil menandakan kondisi batuan di sekitar stasiun tersebut bersifat relatif kompak sedangkan jika nilai koreksi stasiunnya besar maka menandakan kondisi batuan di sekitar stasiun tersebut bersifat relatif lunak. 4. Kesimpulan Dari penelitian penentuan lokasi gempa dengan menggunakan simulated annealing dan guided error search serta penentuan kecepatan gelombang seismik 1-D dengan menggunakan VELEST pada lapangan “Geothermal” dapat disimpulkan sebagai berikut: 1. Gempa lokal yang terekam tersebar pada daerah sesar dan gunung “Tinggi” pada kedalaman 3 km sampai dengan 4.4 km di bawah permukaan laut. 2. Pada kedalaman 3 km sampai dengan 4 km terdapat nilai rasio Vp/Vs yang rendah yang dapat diinterpretasikan sebagai vapour reservoir. 3. Dengan menggunakan data koreksi stasiun dapat di interpretasikan bahwa pada daerah gunung “tinggi” memiliki kondisi batuan yang relatif lebih kompak daripada daerah sesar. Ucapan Terima Kasih Penulis menghaturkan terimakasih kepada PT. Pertamina Geothermal Energy yang telah menyediakan data lapangan pada penelitian ini. 22 Daftar Pustaka Boitnott, G. N., Kirkpatrick, A., 1997. Interpretation of Field Seismic Tomography at The Geysers Geothermal Field, California. in Proc. TwentySecond Workshop on Geothermal Reservoir Engineering., SGP-TR-155, Stanford University, Stanford CA. ELC-Electroconsult.,1983. Feasibility Report Geothermal Area of Mt. Ambang. KomD-5646. Gomberg, J. S., Shedlock, K. M. and Roecker, S. M., 1990. The effect of S-wave arrival times on the accuracy of hypocenter estimation. Bull. Seismol. Soc. Am., 80, 1605–1628. Grandis, H., 2009. Pengantar Pemodelan Inversi Geofisika. Penerbit Himpunan Ahli Geofisika Indonesia., Bandung. Kissling, E., Ellsworth, W. L., Phillips, D. and Kradolfer, Initial reference models earthquake tomography. J. Res., 99, 19635-19646. EberhartU., 1994. in local Geophys. Kissling, E., 1995. Program VELEST USER’S GUIDE – Short Introduction. Institute of Geophysics, ETH Zuerich. 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. rexha.vry@gmail.com 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 -2,98444 0,003117 5 11,05361 6 11,05611 11,01492 -3,0276 0,002157 11,01796 -3,02093 0,001617 7 8 11,01033 11,00713 -3,00043 0,002506 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 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 J. Geofisika Vol. 13 No.1/2012 Pencitraan Struktur 3-D Vp, Vs, Rasio Vp/Vs Menggunakan Tomografi Double Difference di Wilayah Bali Putu Kusuma Yadnya1, Andri Dian Nugraha1, Supriyanto Rohadi2 1) Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan, Institut Teknologi Bandung, Jalan Ganesha 10 Bandung, 40132 2) Balai Besar Meteorologi Klimatologi dan Geofisika, Jakarta, Indonesia Email: kusuma.putu.student@itb.ac.id Abstrak Wilayah Bali memiliki potensi kegempaan yang dipengaruhi oleh zona subduksi di selatan Bali dan back arch thrust fault di utara Bali. Pada studi ini, kami menggunakan katalog data dari jaringan seismograf Onyx yang dioperasikan olah Badan Meteorologi, Klimatologi dan Geofisika (BMKG) menggunakan inversi tomografi double difference (tomoDD) untuk menentukan struktur Vp, Vs dan rasio Vp/Vs di bawah wilayah Bali. Data tersebut meliputi 8 stasiun yang berlokasi di Jawa Timur, Bali dan Lombok yang merekam 1.230 kejadian gempabumi. Kami memilih gempabumi yang minimal direkam oleh 4 stasiun. Pada tahap pertama, lokasi gempa di update menggunakan metoda HypoDD untuk menghilangkan fixed depth gempa dari katalog data dan terpilih 943 kejadian gempa. Pada tahap kedua, kami melakukan inversi tomoDD menggunakan kejadian gempa yang sudah terpilih. Hasil inversi tomoDD menunjukkan adanya anomali kecepatan rendah Vp dan Vs serta nilai Vp/Vs tinggi pada beberapa kedalaman di bawah gunung Agung. Kami menginterpretasikan fitur-fitur ini kemungkinan berhubungan dengan zona pelelehan sebagian dibawah gunung Agung. Slab Indo-Australia yang menunjam di bawah Bali diindikasikan dengan anomali tinggi Vp dan Vs, meskipun anomali ini kurang begitu jelas terlihat. Hiposenter hasil relokasi mengindikasikan keberadaan back arc thrust fault di utara wilayah Bali. Kata Kunci: Bali, gempabumi, tomoDD, struktur kecepatan seismik Abstract Bali regions have seismicity potency influenced by subduction zone in the southern part of Bali and back arc thrust fault in the northern part of Bali. In this study, we used earthquake data catalog from Onyx seismograph network operated by Meteorological Climatological Geophysical Agency of Indonesia (MCGA) for the double difference tomographic (tomoDD) inversion to determine Vp, Vs, and rasio Vp/Vs structure beneath Bali region. The data including 8 stations that were located in East Java, Bali and Lombok areas with1,230 earthquakes event. We selected the event which has recorded at least by 4 stations. The first step, we relocated the events using HypoDD to update the fix depth of the initial location as input for tomoDD and finally we got relocated events of 943 events. The second step, we conducted tomoDD inversion using the relocated events. The results show low velocities anomaly of Vp and Vs are observed beneath Agung and Batur volcanoes and also high Vp/Vs ratio are exhibited at several depths beneath Agung volcano. We interpret this features may be associated with partial melting zone beneath those volcanoes. The subducting Indo-Australia beneath Bali region is indicated by high Vp and high Vs anomalies, however these anomalies are not so clear observed. The relocated events show indication of back arc thrust in northern part of Bali regions. Keywords: Bali, earthquake, tomoDD, seismic velocities structure 32 J. Geofisika Vol. 13 No.1/2012 1. Pendahuluan Pulau Bali merupakan pulau yang terkenal akan keindahan alam dan budaya yang menjadi daya tarik dunia. Dibalik keindahan alam dan budaya tersebut, Bali memiliki tingkat kerawanan gempa yang cukup tinggi. Pulau Bali sendirimemiliki tektonik yang sangat unik. Subduksi lempeng Indo-Australia terhadap lempeng Eurasia dengan kecepatan 7 cm per tahun (Demets dkk., 1994), telah menghasilkan efek berupa struktur geologi sesar aktif di daerah Bali dan sekitarnya. Studi seismisitas lokal daerah Bali hasil pencatatan jaringan seismik lokal yang dilakukan oleh Masturyono (1994) memperoleh hasil analisis bahwa seismisitas gempabumi lokal dan dangkal memberi petunjuk adanya struktur sesar naik belakang busur kepulauan. Data dari jaringan BMKG daerah Bali masih jarang digunakan untuk aplikasi tomografi. Metode tomografi double difference merupakan metode tomografi baru untuk mencitrakan struktur kecepatan dibawah Bali.Studi tomografi di Bali pernah dilakukan oleh Widiyantoro dkk. (2011) menggunakan data regional dan global. Melalui tomografi double difference ini diharapkan dapat menggambarkan dan memahami subduksi dan back arc trust yang ada di Bali. 2. Tektonik Pulau Bali Pulau Bali dan pulau-pulau lain disebelah timur pulau Bali merupakan gugusan kepulauan Sunda kecil yang terbentuk sebagai akibat proses subduksi lempeng IndoAustralia kebawah lempeng Eurasia. Lempeng Indo-Australia diperkirakan bergerak dengan kecepatan 7 cm/tahun dengan arah mendekati normal (Tregoning dkk., 1994).Aktifitas gempa bumi dangkal juga terdapat di daratan Pulau Bali dan cekungan Bali disebelah utara pulau Bali.Cekungan ini terjadi akibat adanya struktur geologi sesar naik belakang busur Silver dkk.(1986). Patahan busur belakang Wetar dan Flores pertama kali dilaporkan oleh Hamilton(1979) berdasarkan beberapa profil refleksi dari Lamont-Doherty. Hamilton (1979) menemukan adanya patahan di utara pulau Alor dan Pantar disisi timur busur belakang zona subduksi Jawa yang biasa dikenal sebagai sesar Sungkup belakang busur Wetar, Flores sampai Sumbawa.Sedangkan Silver dkk.(1986) memperkirakan bahwa patahan tersebut disisi barat berlanjut sampai ke cekungan Bali yang terletak di utara Pulau Bali. 3. Data Gempa Pada penelitian ini digunakan data gempa yang terekam oleh 8 stasiun pada jaringan OnyxBadan Meteorologi Klimatologi dan Geofisika (BMKG) yang berada di Bali, Lombok dan Jawa Timur. Jumlah gempa yang terekam sebanyak 3,457 event yang berasal dari Desember 2005 sampai Februari 2013 (Gambar 1).Oleh karena akurasi hiposenter sangat ditentukan oleh jumlah fasa dan jumlah stasiun yang mencatat maka dilakukan seleksi data pada 3,457 event data Onyx dengan kriteria minimal 4 stasiun yang merekam satu event, sehingga akhirnya diperoleh 1,230 event hasil seleksi. Data 1,230 event ini kemudian dilakukan relokasi gempa menggunakan hypoDD untuk menghilangkan adanya fix depth yang banyak terdapat pada data awal BMKG. Data hasil relokasi hypoDD diperoleh berupa data travel time, dengan jumlah fasa P sebanyak 4,993 dan fasa S sebanyak 3,992. Data diferensialtime yang digunakan dalam penelitian ini setelah diolah dengan program ph2dt menghasilkan data katalog diferensialtime dengan jumlah fasa P sebanyak 6480 dan fasa S sebanyak 5,566. Model kecepatan 1-D yang digunakan dalam penelitian ini mengikuti model kecepatan 1-D Koulakov dkk.(2007) untuk kedalaman< 20 km dan untuk kedalaman> 20 km menggunakan model AK135. Nilai kecepatan di tiap kedalamannya diperlihatkan pada Tabel 1dengan nilai rasio Vp/Vs adalah 1.73 km/s. 33 J. Geofisika Vol. 13 No.1/2012 (1) 34 adalah waktu kejadian (origin Dimana time) dari gempa ke-i, u adalah medan slowness dan ds adalah sebuah elemen dari panjang lintasan. Hubungan antara waktu tiba dan lokasi gempa sangat tidak linier, sehingga penyederhanaan dari deret Taylor biasanya digunakan untuk melinierkan persamaan (1). kelinieran ini menghubungkan misfit antara waktu tiba observasi dan waktu terhadap pertubasi yang sesuai tiba teoritis dengan hiposenter dan parameter struktur kecepatan, dapat dituliskan sebagai: (2) Kedalaman (km) Gambar 1Peta sebaran episenter (bulat berwarna) data katalog Onyx BMKG dari bulan Desember 2005 sampai Februari 2013 sebanyak 3,457 event gempa bumi. Sebaran stasiun(segitiga terbalik warna hitam) dan gunungapi (segitiga merah). Penulisan yang sama juga berlaku untuk gempa ke j yang juga teramati pada stasiun k menjadi: (3) sehingga pengurangan untuk gempa i dan j didapatkan: Tabel 1Model kecepatan 1-D Vp, Vs. Kedalaman (Km) 5 10 15 25 35 45 60 75 90 105 120 165 Vp (km/s) 5 6 6.75 7.11 7.24 7.37 7.6 7.77 7.95 8.01 8.05 8.13 Vs (km/s) 2.89 3.47 3.90 4.11 4.18 4.26 4.39 4.49 4.60 4.63 4.65 4.70 4. Tomografi Double Difference (tomoDD) Persamaan waktu tiba gelombang badan T dari sebuah gempa i ke stasiun k dapat dinyatakan menggunakan teori penjalaran sinar sebagai sebuah integral lintasan (Zhang dan Thurber,2003): 34 (4) Dengan mengasumsikan gempa ini berdekatan satu sama lainnya sehingga lintasan sehingga lintasan dari gempa-gempa ke stasiun tertentu hampir identik dan juga mengasumsikan bahwa struktur kecepatan diketahui, sehingga persamaan (4) dapat disederhanakan menjadi: (5) Dimana adalah selanjutnya disebut double- difference(DD) (Waldhauser dan Ellsworth,2000). adalah perbedaan antara waktu tiba gelombang differensial observasi dan kalkulasi (teoritis) untuk dua buah gempa, dan dapat dituliskan sebagai: J. Geofisika Vol. 13 No.1/2012 )!"# ( )$% (6) Dengan adalah waktu tempuh gelombang gempa i ke stasiun k, dan adalah waktu tempuh gelombang gempa j ke statiun k. Penulisan matriks tomografi mengikuti sebagian notasi Wolfe (2002). Misalkan p = 1,….,P adalah sejumlah gempa dengan Np adalah waktu tiba untuk masing-masing gempa. Untuk tiap-tiap gempa persamaan (2) dapat dituliskan dalam bentuk matriks sebagai berikut: ( &' (' )' * (7) ' Dimana :&' (+' 4)adalah matriks turunan parsial terhadap parameter hiposenter (x,y,z,t0). (' (4 1) adalah vektor pertubasi untuk parameter hiposenter (x,y,z,t0). )' (+' .) adalah matriks turunan model (panjang raypath) terhadap model slowness. *(. 1)adalahVektor pertubasi slowness. Vektor waktu tiba ' (+' 4) adalah residual. Perbedaan utama diantara tomografi konvensional dan DD terutama pada pembentukan matriks yang menggunakan sebuah operator matriks baru QDD . Matriks ini merupakan sebuah operator matriks yang akan memilih event yang akan dihubungkan dalam satu cluster (Zhang dan Thurber,2003). Bentuk matrik dari tomografi doubledifference adalah: /00 &1( /00 )1* /00 1 (8) Penulisan yang lengkap dari persamaan struktur kecepatan dan lokasi gempa dengan menggunakan data absolut dan diferensial / 47 / (9) adalah: 2 3 00 6 3 00 6 1 58 54 dimanaw adalah bobot relatif (relatif weighting) antara waktu tempuh absolut dan relatif, dan I adalah matrik identitas. 5. Hasil dan Interpretasi Hasil inversi tomoDD memiliki dua keluaran yaitu relokasi hiposenter gempa dan struktur kecepatan 3-D (Vp dan Vs) dari daerah penelitian. Pada Gambar 2a-c ditunjukkan posisi hiposenter awal dan setelah relokasi menggunakan tomoDD. Hasil relokasi hiposenter gempa menggunakan tomoDD memperlihatkan bahwa hasil tomoDD dan hasil hypoDD tidak jauh berbeda. Namun dapat dilihat pada posisi awal gempa terdapat fix depth dan setelah relokasi tomoDD nilai fix depth telah hilang. Pada histogram(Gambar 3a-b) ditunjukkan nilai residual time dari metode hypoDD dan tomoDD. Dari histogram ini terlihat bahwa residual time dari tomoDD lebih mendekati nilai nol, hal ini menunjukkan bahwa tomoDD secara statistik mampu meningkatkan akurasi dari posisi hiposenter tersebut.Penampang vertikal melewati daerah penelitian ditunjukkan pada penampang A-B (Gambar 2c) yang melewati Gunung Agung, dimana terlihat gempa mengikuti suatu pola seismisitas. Dari selatan ke utara terlihat adanya suatu pola yang dapat diinterpretasikan sebagai subduksi lempeng Indo-Australia yang menunjam ke lempeng Eurasia. Dari utara ke selatan terlihat adanya suatu pola yang dapat diinterpretasikan sebagai sesar naik belakang busur. 35 J. Geofisika Vol. 13 No.1/2012 Gambar 2a) Posisi gempabumi katalog data Onyx BMKG, b) posisi gempa setelah relokasi dengan hypoDD dan c) posisi gempa setelah relokasi dengan tomoDD. Sebaran stasiun (segitiga hitam terbalik), dan gunungapi (segitiga merah). Gambar 3 HistogramResidualtime (tobs-tcal) a)HypoDD dan b)TomoDD dari inversi data katalog BMKG. Hasil tomogram Vp dan Vs memperlihatkan bahwa adanya anomali negatif yang berada dibawah gunung Agung dan Batur. Sedangkan tomogram rasio Vp/Vs pada daerah tersebut memiliki nilai tinggi pada beberapa kedalaman. Berdasarkan hasil uji checkerboard hanya beberapa kedalaman 36 dibawah gunung Agung dan gunung Batur hasil rekonstruksi checkerboard yang balik ke model awal (Gambar 4). Kekurangan sinar gelombang yang melalui daerah tersebut menyebabkan resolusi pada daerah tersebut kurang baik. J. Geofisika Vol. 13 No.1/2012 Gambar 4 Hasil checkerboard testuntuk Vp (kiri) dan Vs (kanan) pada kedalaman 25 km, 35 km, 45 km, dan 60 km. Warna merah menunjukkan anomali negatif sedangkan warna biru menunjukkan anomali positif, garis putus-putus menunjukkan daerah teresolusi baik. 37 J. Geofisika Vol. 13 No.1/2012 Gambar 5 Penampang horizontal untuk a)Vp,b) Vs,dan c) rasio Vp/Vs pada kedalaman 35 dan 45 km. Warna merah menunjukkan anomali negatif dan warna biru menunjukkan anomali positif. Garis putus-putus hitam menunjukkan daerah dengan resolusi tinggi. a) b) c) d) Gambar 6 a) Sebaran episenter gempa (bulat berwarna), penampang vertikal A - B melalui gunung Agung. Tomogram b) Vp, c) Vs, dan d) rasio Vp/Vs. Warna merah menunjukkan anomali negatif dan warna biru menunjukkan anomali positif untuk Vp dan Vs, sedangkan untuk rasio Vp/Vs warna biru untuk nilai tinggi dan warna merah untuk nilai rendah. Garis putus-putus hitam menunjukkan daerah dengan resolusi tinggi, segitiga merah menunjukkan Gunung Agung, bulatan putih menunjukkan hiposenter gempa. Pada penampang horisontal ditemukan adanya anomali kecepatan negatif pada Vp dan Vs 38 yang berada di bawah gunung Agung dan gunung Batur pada beberapa kedalaman J. Geofisika Vol. 13 No.1/2012 (Gambar 5). Namun pada kedalaman 35 km dan 45 km menunjukkan nilai rasio Vp/Vs tinggi di bawah gunung Agung dan gunung Batur. Hal ini menandakan bahwa pada daerah tersebut kemungkinan terjadi partial melting yang menyebabkan Vp dan Vs turun. Sedangkan pada penampang vertikal yang melewati gunung Agung dan gunung Batur ditemukan adanya anomali Vp dan Vs rendah pada kedalaman 120 km menuju Gunung Agung dan Gunung Batur (Gambar 6). Daerah dengan anomali rendah ini diinterpretasikan sebagai kemungkinan daerah migrasi magma dari slab yang menunjam ke lempeng benua. Material partialmelting lunak dan tidak rigid, sehingga kecepatannya menjadi turun. lempeng Eurasia menerus sampai mantel dalam dan memiliki inklinasi sekitar 70 derajat. Pada penelitian ini sebaran hiposenter yang mengindikasikan adanya subduksi berdasarkan pola event gempa menunjukkan bahwa inklinasi dari tunjaman lempeng sekitar 60 derajat.Sedangkan menurut Daryono(2011) sesar naik belakang busur yang berada di Utara Bali memiliki inklinasi sekitar 30 derajat. Hasil Daryono ini diperkuat oleh penelitian (McCaffrey dan Nabelek, 1987) yang menyatakan kemiringan sesar belakang busur ini sekitar 33 derajat. Hal ini sesuai dengan pola kemenerusan yang berada di Utara Bali yang memiliki inklinasi 30 derajat pada studi ini. Menurut Widiyantoro dkk. (1997) subduksi lempeng Indo-Australia yang menunjam ke Gambar 7 Diagram skematik interpretasi terhadap hasil relokasi dan hasil tomogram Vp di bawah gunung Agung hasil inversi tomoDD pada studi ini. Ilustrasi slab (garis putus-putus biru) diambil dari hasil tomografi global Widiyantoro dkk. (2011). Warna merah menunjukkan anomali negatif dan warna biru menunjukkan anomali positif, sebaran hiposenter gempa (bintang kuning), bulatan biru dan merah menunjukkan aliran magma yang menuju gunung Agung, dan sesar naik belakang busur ditunjukkan oleh garis putus-putus hitam. Pada Gambar 7 ditunjukkan diagram skematik terhadap penampang vertikal yang melewati gunung Agung. Pada gambar tersebut slab yang menunjam ditunjukkan oleh garis putusputus warna biru, gesekan slab kemungkinan menyebabkan adanya partial melting yang terjadipada kedalaman 120 km yang ditunjukkan oleh bulatan biru dan merah sebagai aliran magma yang menuju gunung Agung. Sedangkan sesar naik belakang busur ditunjukkan oleh garis putus-putus hitam yang berada di utara Bali. 39 J. Geofisika Vol. 13 No.1/2012 6. Kesimpulan Dari hasil inversi tomoDD pada studi ini dapat disimpulkan bahwa: 1. Hasil relokasi hiposenter data katalog Onyx BMKG dengan menggunakan metode tomoDD menunjukkan nilai fix depth 15 km, 33 km, 80 km, 160 km, 240 km sudah tidak terlihat lagi. Hasil relokasi hiposenter tersebut menunjukkan adanya pola yang menunjukkan subduksi dari arah selatan Bali dan sesar naik belakang busur yang terjadi di utara Bali. 2. Tomogram Vp dan Vs menunjukkan anomali kecepatan rendah pada daerah dibawah gunung Agung dan gunung Batur, sementara tomogram ratio Vp/Vs menunjukkan anomali tinggi pada beberapa kedalaman. 3. Keberadaan anomali tinggi Vp dan Vs yang menunjukkan adanya subduksi hanya terlihat pada kedalaman 50 sampai 100 km. 4. Anomali kecepatan rendah yang terjadi membentuk jalur migrasi magma antara gunung Agung dan partial melting akibat subduksi dari selatan Bali. Ucapan Terima Kasih Penulis menghaturkan terimakasih kepada BMKG atas katalog data yang diberikan pada studi ini; Dr. H. Zhang untuk program tomoDD;I. Ramadhan atas bantuan dan kerjasamanya pada saat pengerjaan tomoDD; Teknik Geofisika, FTTM, ITB, untuk fasilitas laboratorium yang menunjang penelitian ini. Daftar Pustaka Daryono, 2011, Identifikasi Sesar Naik Belakang Busur (Back Arc Thrust) Daerah Bali Berdasarkan Seismisitas dan Solusi Bidang Sesar. Artikel Kebumian. Badan Meteorologi Klimatologi dan Geofisika. DeMets, C., R. G. Gordon, D.F. Argus and S. Stein, 1994, Effect of Recent to The Geomagnetics Reversal Time Scale on Estimates of Current Plate Motions, Revisions Geophysical Research Letter, 21,2191 - 2194. 40 Hamilton, W., 1979, Tectonic of the Indonesian Region.U.S. Geological Survey Profesional Paper 1078, 345 pp Koulakov, I., M. Bohm, G. Asch, B.-G.Lühr, A. Manzanares, K. S. Brotopuspito, Fauzi, M. A. Purbawinata, N. T. Puspito, A.Ratdomopurbo, H. Kopp, W. Rabbel, dan E. Shevkunova, 2007, P and S velocity structure of the crust and the upper mantle beneath central Java from local tomography inversion, J. Geophys. Res. 112, B08310, doi 10.1029/2006JB004712. Masturyono, 1994, Seismicity of The Bali Region From A Local Seismic Network: Constraints On Bali Back Arc Thrusting. Thesis Master of Science. Rensselaer Polytechnic Institute, New York. McCaffrey, R. dan Nabelek, J., 1987, Earthquakes, Gravity and The Origin of The Bali Basin: An Example of A Naschent Continental Fold and Thrust Belt. J. Geophys. Res., 92, 441-460. Silver,E.A., Breen, N.A. dan H. Prasetyo, 1986, Multibeam Study of the Flores BackArc Thrust Belt, Indonesia. J. Geophys. Res, Vol. 91, No. B3, pp. 3489-3500. Trenggoning, P., F. K. Brunner, Y. Bock, S. S. O. Puntodewo, R. McCaffrey, J. F. Genrich, E. Calais, J. Rais, dan C. First geodetic Subarya, 1994, measurement of convergence across the Geophys, java Trench, Res.lett.,21,2135-2138 Waldhauser, F., dan W. L. Ellsworth, 2000, A double-difference earthquake location algorithm: Method and application to the Northern Hayward Fault, California: Bull. Seism. Soc., 90, 13531368. Widiyantoro, S. dan R. D. Van der Hilst, 1997, Mantle structure beneath Indonesia inferred from highresolution tomographic imaging, J. Geophys. Int.,130,167-182. J. Geofisika Vol. 13 No.1/2012 Widiyantoro, S, Pesicek, J. D., dan Thurber. C. H., 2011, Subducting slab structure below the eastern Sunda arc inferred from non-linear seismic tomographic imaging, Geological Society, London, Special Pub; v.355;p.139-155 doi: 10.1144/SP355.7 Wolfe, C. J., 2002, On the mathematics of using difference operators to relocateEarthquakes, Bull. Seism. Soc. Am., 92, 2879-2892 Zhang, H., dan C. H. Thurber, 2003, Double difference tomography: The method and its application to the Hayward Fault, California: Bull. Seism. Soc. Am., 93,1875-1889. 41 SYARAT DAN FORMAT PENULISAN JURNAL GEOFISIKA Umum Editor menerima artikel berupa hasil penelitian atau hasil studi baik dalam bentuk kajian teroritik maupun eksperimental atau gabungan keduanya dalam bidang seismologi, eksplorasi geofisika, meteorology, oseanografi, geodesi, geologi dan geografi. Naskah harus berisi informasi yang benar, jelas, dan memiliki konstribusi substantif dalam setiap bidang kajian. Penulisan harus singkat dan jelas sesuai dengan format Jurnal Geofisika. Informasi dalam naskah belum pernah dimuat dan tidak sedang dalam proses untuk dimuat di media lain, baik media cetak maupun elektronik. Pengiriman dan Penilaian Naskah Naskah asli yang dikirimkan ke editor harus sesuai dengan format naskah yang ditentukan. Naskah tersebut sebaiknya dikirimkan dalam bentuk softcopy dengan melampirkan biografi singkat penulis pertama, afiliasi, dan alamat lengkap termasuk alamat surat elektronik (email). Setiap makalah yang diterima akan diseleksi oleh Tim Editor untuk kemudian diteruskan kepada Tim Penilai/Penelaah yang memiliki kewenangan untuk mengoreksi, mengembalikan untuk diperbaiki, dan menolak tulisan yang dianggap tidak layak untuk dimuat. Penilaian dilakukan secara objektif dan tertulis oleh minimal dua penilai. Format Penulisan Naskah Format penulisan Jurnal Geofisika dapat dilihat pada halaman berikut di bawah ini. Panduan penulis tersebut sesuai dengan format baku Jurnal Geofisika dan dapat dijadikan sebagai contoh. Bahasa yang digunakan adalah Bahasa Indonesia atau Bahasa Inggris. Bila menggunakan Bahasa Indonesia, menggunakan bahasa yang baku dan benar. Penggunakan bahasa dan istilah asing sebaiknya disertai makna/arti istilah tersebut. J. Geofisika Vol.XX No. XX/20XX Penentuan Judul untuk Jurnal Geofisika dengan Pendekatan Linier dan Eksponensial Penulis Pertama1, Penulis Kedua2, Penulis Ketiga3 1) Afiliasi penulis pertama 2) Afiliasi penulis kedua 3) Afiliasi penulis ketiga Email. penulis.pertama@email.com Abstrak Abstrak berisi latar belakang, tujuan, metodologi, hasil dan kesimpulan secara singkat dan jelas. Jumlah kata dalam abstrak tidak lebih dari 300 kata. Abstrak ditulis dengan huruf Times New Roman dengan ukuran huruf 10. Kata kunci: Tujuan, huruf, jumlah kata Abstract An Abstract consists of background, objective, methodology, results and conclusion. The Abstract should be less then 300 words, in 10 point Italic Times New Roman font. Keywords: Objective, font, words number 1. Struktur Naskah 2.2 Huruf dan Spasi Struktur naskah adalah Judul, Nama Penulis (tanpa gelar), afiliasi tempat kerja, email, abstrak, kata kunci, pendahuluan/latar belakang dan tujuan, data dan metodologi, hasil dan diskusi, kesimpulan, ucapan terima kasih dan daftar pustaka. Jurnal Geofisika menggunakan beberapa ketentuan untuk huruf dan sepasi sebagai berikut: 2. Format Makalah • • 2.1 Tata Letak Naskah diketik dengan format kertas A4 dan setiap halaman diberi nomor dengan panjang naskah antara 8 – 15 halaman. Menggunakan margin kertas sebagai berikut: • • Margin atas 3 cm, bawah 3 cm, kiri ,5 cm dan kanan 2,5 cm. Kecuali pada bagian abstrak margin kiri dan kanan masing-masing 4 dan 3,5 cm. Badan naskah ditulis dengan dua kolom dengan lebar kolom 7,25 dan jarak antar kolom 1 cm. • • Badan naskah diketik 1 spasi dengan huruf Times New Roman dengan ukuran huruf 11 Judul makalah dicetak tebal dengan huruf besar pada awal kata menggunakan Times New Roman dengan ukuran huruf 14 Nama dan afiliasi penulis berturutturut dengan Times New Roman dengan ukuran huruf 10 Abstrak diketik dengan menggunakan Times New Roman dengan ukuran huruf 10, dan untuk bahasa Inggris Abstract diketik dengan huruf miring. 3. Bahasa, Satuan dan Persamaan Bahasa yang digunakan adalah Bahasa Indonesia yang baik dan benar atau juga dapat menggunakan Bahasa Inggris dengan mengikuti aturan tata bahasa (grammar) yang baik dan benar. J. Geofisika Vol.XX No. XX/20XX Penggunaan singkatan dan tanda-tanda diusahkan mengikuti aturan nasional dan/atau internasional dengan mengacu pada system satuan internasional (SI). Setiap persamaan harus diketik dan diberi nomor seperti pada contoh berikut. (1) P( E) =exp dimana E adalah perubahan fungsi objektif atau perubahan misfit akibat perturbasi model tersebut. 4. Tata Letak Gambar Gambar dapat dimasukkan dalam kolom atau meliputi kedua kolom untuk gambar yang lebig besar. Legenda gambar harus terlihat dengan jelas. Keterngan gambar ditulis dengan jelas dibawah gambar dengan ukuran huruf Times New Roman 10 seperti contoh berikut. 2 Vp Vs 1 0 -1 K edalam an (k m ) -2 -3 -4 -5 -6 -7 -8 -9 0 2 4 6 8 Kecepatan Gelombang Seismik (km/s) Gambar 1. Plot model kecepatan untuk gelombang P (Vp) dan gelombang S (Vs) yang digunakan pada studi ini. Ucapan Terima Kasih Ucapan terima kasih ditulis terpisah tanpa nomor 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.