Aplikasi Metode Inversi Simulated Annealing pada Penentuan Hiposenter Gempa Mikro dan Tomografi Waktu Tunda 3-D Struktur Kecepatan Seismik untuk Studi Kasus Lapangan Panas Bumi “RR” Tugas Akhir Diajukan sebagai salah satu syarat untuk menempuh ujian strata satu Program Studi Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan, Institut Teknologi Bandung
Disusun oleh : Rexha Verdhora Ry 123 09 019
Program Studi Teknik Geofisika Fakultas Teknik Pertambangan dan Perminyakan Institut Teknologi Bandung 2013
Aplikasi Metode Inversi Simulated Annealing pada Penentuan Hiposenter Gempa Mikro dan Tomografi Waktu Tunda 3-D Struktur Kecepatan Seismik untuk Studi Kasus Lapangan Panas Bumi “RR” Tugas Akhir Diajukan sebagai salah satu syarat untuk menempuh ujian strata satu Program Studi Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan, Institut Teknologi Bandung
Disusun oleh : Rexha Verdhora Ry 123 09 019
disetujui dan disahkan: Bandung, September 2013
Pembimbing
Dr. Andri Dian Nugraha NIP. 19780908 200912 1 001
ABSTRAK
Pengamatan aktivitas gempa mikro pada area geotermal digunakan untuk mendeteksi permeabilitas struktur. Untuk itu diperlukan penentuan lokasi hiposenter, dimana proses tersebut melibatkan pencarian solusi hiposenter yang memiki error minimum antara waktu tempuh observasi dan kalkulasi. Sebelumnya, lokasi hiposenter pada lapangan panas bumi “RR” telah ditentukan menggunakan metode Geiger, dimana teknik optimisasi lokal seperti ini bergantung pada model awal. Bagaimanapun, pada penelitian ini, metode simulated annealing diterapkan pada data dan model kecepatan 1-D yang sama untuk merelokasi hiposenter dan meminimalisasi fungsi error. Waktu tempuh dihitung menggunakan ray tracing metode shooting. Hasil penelitian ini menunjukkan bahwa lokasi hiposenter memiliki error RMS yang lebih kecil dibandingkan dengan penelitian sebelumnya, yang secara statistik berasosiasi dengan solusi yang lebih baik. Selanjutnya, data lokasi hiposenter yang baru akan digunakan sebagai input untuk membuat struktur kecepatan seismik 3-D dengan tomografi seismik untuk Vp, Vs, dan rasio Vp/Vs. Waktu tempuh pada model kecepatan 3-D dihitung menggunakan ray tracing metode pseudo-bending. Metode inversi yang digunakan pada pemodelan tomografi adalah iterative damped least square dan simulated annealing. Hasilnya menunjukkan bahwa metode iterative damped least square menghasilkan solusi model kecepatan yang lebih baik, ditinjau secara visualisasi dan statistik. Hasil pemodelan tomografi menunjukkan keberadaan Vp/Vs rendah yakni sekitar 1.45 – 1.65 pada area reservoir di elevasi -1 sampai -3 km (MSL = 0 km). Hal ini diinterpretasikan sebagai steam-saturated rock pada area reservoir lapangan panas bumi “RR”. Keberadaan area reservoir ini didukung oleh data well-trajectory, dimana zona Vp/Vs tinggi berada di sekitar sumur injeksi dan zona Vp/Vs rendah berada di sekitar sumur produksi. Jika dibandingkan dengan penelitian lain pada lapangan ini, ada kemungkinan bahwa sistem fasa reservoir tersebut mengalami perubahan dari water-saturated ke steam-saturated. Kata kunci : gempa mikro, seismik tomografi, inversi geofisika, ray tracing. i
ABSTRACT
Observation of micro-earthquake activity in the area of geothermal is used to detect the permeability of the structure. It is necessary for determining the location of hypocenter which the process involves finding a hypocenter location that has minimum error between the observed and the calculated travel times. Previously, hypocenter location at “RR” Geothermal Field has been determined by Geiger’s method, where this local optimization technique depends on the initial model. However, in this study, simulated annealing method was applied on same data and 1D velocity model to relocate hypocenter and minimize error function. The travel times were calculated using ray tracing shooting method. Our results show hypocenter location has smaller RMS error compared to the previous study that can be statistically associated with better solution. Futhermore, the new hypocenter location data will be used as input to produce 3-D seismic velocity structure for Vp, Vs, and the ratio Vp/Vs using seismic tomography. The travel times through 3-D velocity model were calculated using ray tracing pseudo-bending method. Inversion methods which used for tomographic modeling are iterative damped least square and simulated annealing. Our results show that iterative damped least square produced better velocity model, reviewed by visualization and statistic. Tomographic modeling results indicate the presence of low Vp/Vs at around 1.45 – 1.65 in the reservoir area at elevation of -1 to -3 km (MSL = 0 km). This is interpreted as steam-saturated rock in the reservoir area of “RR” geothermal field. The existences of the reservoir area are supported by the data of well-trajectory, where the zones of high Vp/Vs are around the injection wells and the zones of low Vp/Vs are around the production wells. When compared with other studies in this field, it is possible that the reservoir’s phase system has changed from watersaturated to steam-saturated. Keywords: micro-earthquakes, seismic tomography, geophysical inversion, ray tracing.
ii
KATA PENGANTAR
Puji syukur penulis panjatkan kehadirat Allah SWT, karena atas Rahmat-Nya penyusunan tugas akhir ini dapat diselsesaikan pada waktu yang tepat. Tugas akhir yang berjudul Aplikasi Metode Inversi Least Square dan Simulated Annealing pada Penentuan Lokasi Hiposenter dan Tomografi 3-D Struktur Kecepatan Seismik untuk Studi Kasus Lapangan Panas Bumi “RR” ini disusun sebagai salah satu syarat kelulusan untuk jenjang sarjana Program Studi Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan, Institut Teknologi Bandung. Tak lupa penulis juga mengucapkan terima kasih kepada pihak-pihak yang telah membantu penulis dalam menyelesaikan tugas akhir ini. Ucapan terima kasih penulis diberikan kepada: 1. Dr. Andri Dian Nugraha sebagai dosen dan pembimbing penulis yang membantu dalam memberikan ide dan hal-hal baru tentang tugas akhir ini. 2. Keluarga penulis, Ayah, Bunda, dan kedua Abang dari penulis yang memberikan dukungan serta doa sepenuhnya sejak kecil hingga sekarang. 3. Seluruh Bapak dan Ibu dosen Program Studi Teknik Geofisika ITB atas ilmu yang telah diberikan selama penulis menempuh pendidikan. 4. Dr. Hendra Grandis yang telah mengenalkan dan mengajarkan penulis kepada pemograman dari permasalahan geofisika yang sederhana hingga permasalahan yang kompleks, serta berdiskusi perihal tugas akhir ini. 5. Bu Lilik, Bu Ning, Pak Pur, Pak Agus, Pak Dedi dan seluruh karyawan Program Studi Teknik Geofisika. 6. Tania, yakni sahabat, inspirasi, teman diskusi dan tempat berbagi ilmu yang tiada hentinya memberikan masukan dan semangat hingga akhir penyusunan tugas akhir ini. 7. Teman seperjuangan laboratorium Geofisika Dekat Permukaan, Naufal, Kusuma, Wieke, Lukman, Nando, Fajar, Ruth, Yuni, Rani, Qadry, Kak Vega, Mas Rendy, dan Mas Billy yang saling memberi masukan dan semangat. 8. Kakak-kakak kelas yang membantu penyusunan tugas akhir ini, Ame, Saladin, Waskito, dan Kak Mia. iii
9. Teman-teman seperjuangan Oktober, Joshua, Vanno, Anno, Yoyo, Gadih, Harya, Firdaus, Dodi, Fandy, Aza, Fikka, Tina, dan lain-lain yang memacu dalam pengerjakan tugas akhir ini. 10. Sahabat-sahabat Beskem, Dhiar, Rafi, dan Aldo, yang selalu setia mendukung pengerjaan tugas akhir ini. 11. Teman-teman satu angkatan Teknik Geofisika ITB 2009, keluarga terindah penulis di Bandung. 12. Teman-teman HIMA-TG “TERRA” ITB. Semoga menjadi himpunan yang selalu bersinar dan jaya. 13. Serta semua pihak yang tidak dapat disebutkan namanya satu-persatu. Penulis menyadari bahwa tugas akhir ini masih jauh dari kesempurnaan, oleh karenanya penulis mengharapkan saran ataupin kritik yang dapat membangun penulis menjadi lebih baik. Semoga tugas akhir ini dapat memberikan manfaat dalam pengembangan ilmu pengetahuan di tanah air tercinta ini, Indonesia. Terima kasih.
Bandung, Septermber 2013
Rexha Verdhora Ry
iv
DAFTAR ISI
HALAMAN JUDUL HALAMAN PENGESAHAN ABSTRAK …………………………………………………………………………... i ABSTRACT ………………………………………………………………………... ii KATA PENGANTAR ……………………………………………………………... iii DAFTAR ISI ………………………………………………………………………... v DAFTAR GAMBAR …………………………………………………………….... vii DAFTAR TABEL ………………………………………………………………….. xi BAB I PENDAHULUAN ………………………………………………………….. 1 1.1
Latar Belakang ……………………………………………………………… 1
1.2
Tujuan Penelitian …………………………………………………………… 2
1.3
Ruang Lingkup Masalah ……………………………………………………. 3
1.4
Metodologi Penelitian ………………………………………………………. 3
1.5
Sistematika Penyajian ………………………………………………………. 4
BAB II TEORI DASAR ……………………………………………………………. 5 2.1
Sistem Geotermal dan Hidrotermal ………………………………………… 5
2.2
Penentuan Lokasi Hiposenter ………………………………………………. 8
2.3
2.2.1
Ray Tracing Metode Shooting ……………………………………… 9
2.2.2
Inversi Simulated Annealing ………………………………………. 10
Seismik Tomografi ………………………………………………………… 13 2.3.1
Parameterisasi Model ……………………………………………… 14
2.3.2
Ray Tracing Metode Pseudo-Bending …………………………….. 14
2.3.3
Tomografi Waktu Tunda …………………………………………... 16
2.3.4
Inversi Iterative Damped Least Square …………………………… 17
2.3.5
Tes Resolusi ……………………………………………………….. 20
v
BAB III PENENTUAN LOKASI HIPOSENTER ………………………………... 21 3.1
Penelitian Terdahulu ………………………………………………………. 21
3.2
Persiapan Data …………………………………………………………….. 24
3.3
Perhitungan Waktu Tempuh Gelombang ………………………………….. 26
3.4
Algoritma Inversi ………………………………………………………….. 27
3.5
Pengolahan Data Sintetis ………………………………………………….. 30
3.6
Analisis dan Pembahasan ………………………………………………….. 32
BAB IV SEISMIK TOMOGRAFI 3-D …………………………………………… 37 4.1
Persiapan Data …………………………………………………………….. 37
4.2
Parameterisasi Model ……………………………………………………… 39
4.3
Perhitungan Waktu Tempuh Gelombang ………………………………….. 40
4.4
Algoritma Inversi ………………………………………………………….. 42
4.5
4.6
4.7
4.4.1
Algoritma Iterative Damped Least Square ………………………... 42
4.4.2
Algoritma Simulated Annealing …………………………………… 42
Tes Resolusi ……………………………………………………………….. 43 4.5.1
CRT Iterative Damped Least Square ……………………………… 44
4.5.2
CRT Simulated Annealing ………………………………………… 48
Pemodelan Inversi Tomografi ……………………………………………... 52 4.6.1
Aplikasi Iterative Damped Least Square ………………………….. 52
4.6.2
Aplikasi Simulated Annealing ……………………………………... 59
Analisis dan Pembahasan ………………………………………………….. 66 4.7.1
Perbandingan Solusi Model Kecepatan Kedua Metode Inversi …... 66
4.7.2
Analisis Struktur Kecepatan Vp, Vs, dan rasio Vp/Vs ……………. 72
BAB V PENUTUP ………………………………………………………………... 85 5.1
Kesimpulan ………………………………………………………………... 85
5.2
Saran ………………………………………………………………………. 86
DAFTAR PUSTAKA ……………………………………………………………... 87
vi
DAFTAR GAMBAR
Gambar 2.1
Skema umum dari suatu sistem hidrotermal ………………………... 6
Gambar 2.2
Ilustrasi sudut datang dan transmisi ray pada batas lapisan ………… 9
Gambar 2.3
Ilustrasi iterasi sudut dari sumber ke stasiun ………………………. 10
Gambar 2.4
Diagram alir metode simulated annealing secara umum ..……….... 13
Gambar 2.5
Diagram alir perhitungan waktu tempuh gelombang menggunakan ray tracing pseudo-bending pada penelitian ini .…................................. 15
Gambar 2.6
Diagram alir pemodelan inversi tomografi dengan metode iteratve damped least square pada penelitian ini .…...................................... 19
Gambar 2.7
Diagram alir tes resolusi checkerboard pada penelitian ini .…......... 20
Gambar 3.1
Peta persebaran event yang dihasilkan oleh penelitian sebelumnya dan stasiun ………………………............................................................ 22
Gambar 3.2
Peta penampang U-S terhadap kedalaman dari persebaran event dan stasiun ………………………............................................................ 23
Gambar 3.3
Peta penampang B-T terhadap kedalaman dari persebaran event dan stasiun ………………………............................................................ 23
Gambar 3.4
Plot model kecepatan untuk gelombang P (Vp) dan gelombang S (Vs) ………………………........................................................................ 25
Gambar 3.5
Diagram alir perhitungan waktu tempuh gelombang (metode shooting) pada penelitian ini .…........................................................................ 26
Gambar 3.6
Skema perhitungan waktu tempuh gelombang total dari sumber ke stasiun penerima ………………………............................................ 27
Gambar 3.7
Diagram alir metode inversi simulated annealing pada penentuan lokasi hiposenter penelitian ini .…...……….……………………… 29
Gambar 3.8
Peta persebaran hiposenter hasil inversi, hiposenter sebenarnya, dan stasiun pengamat ………………………........................................... 31
Gambar 3.9
Histogram error RMS dari 20 kali perhitungan …………………… 31
Gambar 3.10 Peta persebaran event yang dihasilkan oleh penelitian ini menggunakan simulated annealing dan stasiun …………………… 33 Gambar 3.11 Peta penampang U-S terhadap kedalaman dari persebaran event dan stasiun ………………………............................................................ 34
vii
Gambar 3.12 Peta penampang B-T terhadap kedalaman dari persebaran event dan stasiun ………………………............................................................ 34 Gambar 3.13 Peta persebaran event yang dihasilkan oleh penelitian sebelumnya dan penelitian ini menggunakan simulated annealing (bulat merah) beserta stasiun ……………………………………………………………… 35 Gambar 3.14 Histogram distribusi waktu residual dan error RMS solusi hiposenter gempa mikro yang dihasilkan dari penelitian sebelumnya dan simulated annealing pada penelitian ini …….................................... 36 Gambar 4.1
Distribusi event dan area data hiposenter yang digunakan ………... 37
Gambar 4.2
Area penelitian dan parameterisasi model blok tiga dimensi pada penelitian ini ………………………………………………..……… 39
Gambar 4.3
Peta distribusi cakupan ray path, hiposenter, dan stasiun hasil ray tracing metode pseudo-bending melalui model kecepatan awal gelombang P ……………………….................................................. 40
Gambar 4.4
Peta distribusi cakupan ray path, hiposenter, dan stasiun hasil ray tracing metode pseudo-bending melalui model kecepatan awal gelombang S ……………………….................................................. 41
Gambar 4.5
Diagram alir metode inversi simulated annealing pada pemodelan tomografi penelitian ini ……..…....................................................... 43
Gambar 4.6
Model awal Checkerboard Resolution Test ……………...…........... 44
Gambar 4.7
Tomogram slice horizontal dari tes resolusi checkerboard (Iterative Damped Least Square) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km … 45
Gambar 4.8
Tomogram slice vertikal Barat - Timur dari tes resolusi checkerboard (Iterative Damped Least Square) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km beserta daerah yang teresolusi baik (garis hitam putusputus) ………………………............................................................. 46
Gambar 4.9
Tomogram slice vertikal Utara - Selatan dari tes resolusi checkerboard (Iterative Damped Least Square) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km beserta daerah yang teresolusi baik (garis hitam putus-putus) ………………………................................................... 47
Gambar 4.10 Histogram distribusi delay time dari CRT iterative damped least square ………………………............................................................ 48 Gambar 4.11 Histogram distribusi delay time dari CRT simulated annealing …... 48 Gambar 4.12 Tomogram slice horizontal dari tes resolusi checkerboard (Simulated Annealing) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km ……………... 49 Gambar 4.13 Tomogram slice vertikal Barat - Timur dari tes resolusi checkerboard (Simulated Annealing) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km ………………………........................................................................ 50 viii
Gambar 4.14 Tomogram slice vertikal Utara - Selatan dari tes resolusi checkerboard (Simulated Annealing) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km ………………………............................................................. 51 Gambar 4.15 Tomogram slice horizontal dari struktur kecepatan gelombang P (Iterative Damped Least Square) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km ……………………….................................................................. 53 Gambar 4.16 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang P (Iterative Damped Least Square) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km ……………………….............................. 54 Gambar 4.17 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang P (Iterative Damped Least Square) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km ……………………….......................................... 55 Gambar 4.18 Tomogram slice horizontal dari struktur kecepatan gelombang S (Iterative Damped Least Square) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km ……………………….................................................................. 56 Gambar 4.19 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang S (Iterative Damped Least Square) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km ……………………….............................. 57 Gambar 4.20 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang S (Iterative Damped Least Square) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km ……………………….......................................... 58 Gambar 4.21 Tomogram slice horizontal dari struktur kecepatan gelombang P (Simulated Annealing) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km … 60 Gambar 4.22 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang P (Simulated Annealing) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km ………………………............................................. 61 Gambar 4.23 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang P (Simulated Annealing) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km ………………………............................................................. 62 Gambar 4.24 Tomogram slice horizontal dari struktur kecepatan gelombang S (Simulated Annealing) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km … 63 Gambar 4.25 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang S (Simulated Annealing) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km ………………………............................................. 64 Gambar 4.26 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang S (Simulated Annealing) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km ………………………............................................................. 65 Gambar 4.27 Cuplikan tomogram slice vertikal Barat – Timur ………………….. 67
ix
Gambar 4.28 Cuplikan tomogram inversi iterative damped least square dan simulated annealing ………………………...................................... 68 Gambar 4.29 Perbandingan tomogram horizontal untuk iterative damped least square, simulated annealing, dan selisih absolut dari kedua solusi .. 69 Gambar 4.30 Perbandingan tomogram vertikal B-T untuk iterative damped least square, simulated annealing, dan selisih absolut dari kedua solusi .. 70 Gambar 4.31 Histrogram distribusi dt (delay time) untuk solusi model kecepatan Vp dari inversi menggunakan iterative damped least square dan simulated annealing ………………………....................................................... 71 Gambar 4.32 Histrogram distribusi dt (delay time) untuk solusi model kecepatan Vs dari inversi menggunakan iterative damped least square dan simulated annealing ………………………....................................................... 72 Gambar 4.33 Tomogram slice horizontal dari struktur Vp, struktur Vs, rasio Vp/Vs pada z = 1, 0, -1, dan -2 km ………………………........................... 73 Gambar 4.34 Tomogram slice horizontal dari struktur Vp, struktur Vs, rasio Vp/Vs pada z = -3, -4, -5, dan -6 km ………………………........................ 74 Gambar 4.35 Tomogram slice vertikal Barat - Timur dari struktur Vp, struktur Vs, rasio Vp/Vs pada B-T 8, 9.5, 11, dan 12.5 km ……………….......... 75 Gambar 4.36 Tomogram slice vertikal Barat - Timur dari struktur Vp, struktur Vs, rasio Vp/Vs pada B-T 14, 15.5, 17, dan 18.5 km …………….......... 76 Gambar 4.37 Tomogram slice vertikal Utara - Selatan dari struktur Vp, struktur Vs, rasio Vp/Vs pada U-S 7, 8, 9, dan 10 km …………………….......... 77 Gambar 4.38 Tomogram slice vertikal Utara - Selatan dari struktur Vp, struktur Vs, rasio Vp/Vs pada U-S 11, 12, 13, dan 14 km ……………………… 78 Gambar 4.39 Peta persebaran event, stasiun, well- trajectory, dan batas reservoir 80 Gambar 4.40 Tomogram slice vertikal Utara - Selatan dari struktur Vp, struktur Vs, rasio Vp/Vs pada U-S 9, 10, 11, dan 12 km beserta well-trajectory ………………………........................................................................ 81 Gambar 4.41 Tomogram slice vertikal dari rasio Vp/Vs pada U-S 10 km, beserta interpertasi zona kondensasi uap, cap rock, zona reservoir ……….. 82 Gambar 4.42 Model volumetri zona Vp/Vs rendah dan tinggi beserta well-trajectory dari sumur produksi dan sumur injeksi ………………………………………... 83
Gambar 4.43 Model volumetri zona Vp/Vs rendah beserta sumur produksi …….……... 84 Gambar 4.44 Model volumetri zona Vp/Vs tinggi beserta sumur injeksi ….....……….... 84
x
DAFTAR TABEL
Tabel 3.1 Contoh katalog data fasa gelombang P ……………….......................... 24 Tabel 3.2 Katalog model kecepatan 1-D Vp dan Vs ……………………….......... 25 Tabel 3.3 Katalog data sintetis waktu tempuh gelombang ……………………… 30 Tabel 3.4 Katalog solusi hasil inversi pada data sintetis ……………………….... 30 Tabel 4.1 Contoh format katalog data fasa gelombang P ……….......................... 38 Tabel 4.2 Contoh format katalog data fasa gelombang S ……….......................... 38
xi
BAB I PENDAHULUAN
1.1
Latar Belakang
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, lokasi hiposenter gempa 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. 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 satu dimensi yang sama untuk menentukan lokasi hiposenter dan meminimumkan fungsi misfit. Kedua solusi akan dibandingkan dan lokasi hiposenter yang memiliki root mean square (RMS) error lebih kecil akan berasosiasi dengan solusi yang lebih baik secara statistik.
1
Analisis gempa mikro kemudian dilakukan dengan pemodelan tomografi. Pemodelan tomografi tersebut mengaplikasikan metode inversi iterative damped least square dan simulated annealing. Metode iterative damped least square sendiri telah banyak digunakan untuk menyelesaikan pemodelan tomografi dan telah berhasil mencitrakan struktur kecepatan seismik bawah permukaan dengan baik, sedangkan metode simulated annealing masih jarang digunakan untuk menyelesaikan permasalahan inversi non-linier untuk banyak parameter seperti tomografi ini. Pada penelitian ini, metode simulated annealing akan diimplementasikan pada pemodelan tomografi. Kemudian, solusinya akan dibandingkan dengan solusi yang diperoleh dari metode inversi yang lebih umum digunakan yaitu metode iterative damped least square. Hasil inversi tomografi yang diperoleh merupakan model kecepatan tiga dimensi dari bawah permukaan untuk gelombang P (Vp) dan gelombang S (Vs). Keberadaan fasa fluida atau leburan pada daerah penelitian dapat diindikasi dari anomali kecepatan dan pengamatan rasio Vp/Vs.
1.2
Tujuan Penelitian
Tujuan dari penelitian tugas akhir ini adalah sebagai berikut: 1. Menerapkan metode inversi simulated annealing pada penentuan lokasi hiposenter. 2. Membandingkan solusi lokasi hiposenter yang diperoleh dari metode simulated annealing dengan penelitian sebelumnya (metode Geiger). 3. Menerapkan metode inversi iterative damped least square dan simulated annealing untuk mencitrakan struktur tiga dimensi bawah permukaan kecepatan gelombang seismik (Vp, Vs, dan rasio Vp/Vs) pada lapangan geotermal “RR”. 4. Membandingkan solusi model struktur kecepatan seismik yang diperoleh dari metode iterative damped least square dengan metode simulated annealing. 5. Menginterpretasikan
struktur
kecepatan
gelombang
seismik
untuk
mengetahui kondisi dan karakteristik reservoir lapangan geotermal “RR”.
2
1.3
Ruang Lingkup Masalah
Ruang lingkup masalah dalam penelitian tugas akhir ini adalah sebagai berikut: 1. Data gempa mikro yang digunakan merupakan data lapangan geotermal “RR” yang direkam selama satu tahun. 2. Asumsi bahwa hasil picking waktu tiba gelombang dari penelitian sebelumnya sudah akurat sehingga penyelesaian masalah difokuskan pada metode inversi yang digunakan. 3. Penerapan inversi tidak menggunakan pembobotan nilai. 4. Perbandingan solusi inversi dilihat dari error RMS. 5. Penentuan waktu tempuh gelombang dari sumber ke penerima pada proses penentuan hiposenter menggunakan metode ray tracing shooting, sedangkan pada tomografi model kecepatan gelombang seismik menggunakan metode ray tracing pseudo-bending. 6. Luas area penelitian untuk pemodelan tomografi adalah 15 x 10 x 8 km3. 7. Pada penelitian ini, nama lapangan geotermal, koordinat lokasi hiposenter, dan koordinat lokasi stasiun tidak ditampilkan atau disamarkan.
1.4
Metode Penelitian
Metode penelitian yang digunakan selama penyusunan tugas akhir ini secara sistematis ialah: 1. Studi Pustaka Berapa literatur terkait yang dipelajari dalam penelitian ini antara lain: (a) teori dan konsep metode inversi least square dan simulated annealing, (b) teori dan algoritma ray tracing untuk metode shooting dan pseudo-bending, (c) konsep inversi penentuan lokasi hiposenter dan pemodelan tomografi 3 dimensi, dan (d) konsep sistem panas bumi dan interpretasi. 2. Persiapan Data Pengumpulan data gempa mikro di lapangan geotermal “RR” selama satu tahun berupa waktu tiba gelombang P dan S, model kecepatan 1-D untuk fasa gelombang P dan S, serta solusi hiposenter hasil penelitian sebelumnya. 3
3. Pemrograman Pemrograman dilakukan dengan menggunakan perangkat lunak MatLab. Pemrograman tersebut terdiri dari perhitungan waktu tempuh gelombang menggunakan ray tracing 1-D shooting method, proses inversi penentuan lokasi hiposenter menggunakan metode simulated annealing, perhitungan waktu tempuh gelombang menggunakan ray tracing 3-D metode pseudobending, proses inversi tomografi menggunakan metode damped least square, dan proses inversi tomografi menggunakan metode simulated annealing.
1.5
Sistematika Penyajian
Penulisan laporan tugas akhir ini dibagi menjadi 5 bab, yaitu: 1. Bab I Pendahuluan Bab ini membahas latar belakang, tujuan penelitian, ruang lingkup masalah, metode penelitian, serta sistematika penyajian. 2. Bab II Teori Dasar Bab ini membahas landasan teori yang digunakan, terdiri dari sistem geotermal dan hidrotermal, penentuan lokasi hiposenter, seismik tomografi, perhitungan waktu tempuh gelombang, inversi simulated annealing dan inversi iterative damped least square, serta tes resolusi 3. Bab III Penentuan Lokasi Hiposenter Bab ini membahas pengolahan data dan hasil dalam penentuan lokasi hiposenter, terdiri dari tahapan penelitian terdahulu, persiapan data, perhitungan waktu tempuh gelombang, algoritma inversi, pengolahan data sintetis, dan analisis pembahasan. 4. Bab IV Seismik Tomografi 3-D Bab ini membahas pengolahan data dan hasil dalam pemodelan tomografi 3D, terdiri dari tahapan persiapan data, parameterisasi model, perhitungan waktu tempuh gelombang, algoritma inversi, tes resolusi, pemodelan inversi tomografi, dan analisis pembahasan. 5. Bab V Kesimpulan dan Saran Bab ini membahas kesimpulan yang diperoleh dari aplikasi 2 metode inversi untuk penentuan lokasi hiposenter dan pemodelan tomografi 3-D, serta saran untuk penelitian selanjutnya.
4
BAB II TEORI DASAR
2.1
Sistem Geotermal dan Hidrotermal
Sistem geotermal merupakan suatu istilah umum yang menggambarkan proses perpindahan panas secara alami dalam suatu suatu volume tertutup di kerak bumi, dimana panas berpindah dari suatu sumber panas menuju daerah pelepasan panas seperti permukaan bebas (Hochstein dan Browne, 2000). Terdapat lima tipe utama dari sistem geotermal berdasarkan kriteria geologi, geofisika, hidrologi, dan teknik (Goff dan Janik, 2000). Tiga tipe yang pertama memiliki reservoir berisi air panas secara alami dan biasa disebut sistem hidrotermal. Dua tipe lainnya membutuhkan pemompaan fluda terutama air ke bawah permukaan (injeksi fluida) kemudian baru diproduksi untuk mengekstraksi air panas. Kelima tipe tersebut adalah: 1. Sistem young igneous, berasosiasi dengan vulkanisme kuarter dan intrusi magma. Sekitar 95% dari aktivitas vulkanik terjadi di batas lempeng dan zona hot spot. Lingkungan seperti itu berasosiasi dengan kebanyakan aktivitas tektonik dan seismisitas. Sistem geotermal ini secara umum merupakan yang terpanas dengan temperatur ≤370oC dan kedalaman reservoir secara umum ≤1.5 km. 2. Sistem tektonik, berasosiasi dengan peningkatan seismisitas karena sesar dan peningkatan aliran panas oleh kerak tipis, namun tidak memiliki batuan beku. Sistem ini terjadi pada lingkungan backarc, daerah crustal extension, zona kolisi, dan sepanjang zona sesar. Sistem ini secara umum memiliki temperatur reservoir ≤250oC dan terdapat di kedalaman ≥1.5 km. 3. Sistem geopressure, ditemukan di cekungan sedimen dimana penurunan dan penguburan dalam dari lapisan berisi fluida menjadi panas karena reservoir overpressured. Aliran panas dan seismistas secara umum rendah sampai normal. Kebanyakan sistem ini memiliki karakteristik menyerupai lapangan minyak dan gas. Sistem ini terdapat di kedalaman yang cukup dalam yaitu 1.5 sampai 3 km dan memiliki temperatur 50 sampai 190oC. 5
4. Sistem hot dry rock, memanfaatkan panas yang tersimpan dalam batuan berporositas rendah dan batuan tidak permeabel pada variasi kedalaman dan temperatur. Sistem ini berasosiasi dengan gradien termal di bawah permukaan. Sistem ini memiliki temperatur berkisar antara 120 sampai 225 oC dengan kedalaman 2 sampai 4 km. 5. Sistem magma tap, memanfaatkan panas yang keluar dari tubuh magma dangkal. Sama seperti pada sistem hot dry rock, sistem ini membutuhkan instalasi sistem injeksi fluida untuk memproduksi fluida panas. Pada sistem ini, magma merupakan bentuk paling murni panas alami yang memiliki temperatur ≤1200oC. Sistem hidrotermal merupakan salah satu bentuk lain dari sistem geotermal dimana transfer panas dari sumber panas menuju permukaan merupakan proses konveksi bebas melibatkan fluida meteorik dengan atau tanpa fluida magmatik (Hochtein dan Browne, 2000). Fluida terlepas (discharge) di dekat permukaan lewat suatu rekahan, kemudian diisi kembali oleh air meteorik yang berasal dari luar (recharge). Suatu sistem hidrotermal terdiri dari sumber panas (heat source), reservoir dengan fluida panas, batuan penutup (cap rock), daerah pengisian kembali (recharge), dan daerah pelepasan (discharge) dengan manifestasi di permukaan.
Gambar 2.1 Skema umum dari suatu sistem hidrotermal (dalam Dickson, 2004).
6
Berhubungan dengan sistem geotermal, perubahan kecepatan gelombang seismik akan dipengaruhi oleh sifat fisis dari batuan penyusun sistem tersebut. Peningkatan temperatur akan menurunkan kecepatan gelombang seismik baik kecepatan gelombang P maupun kecepatan gelombang S (Wang dkk., 1990). Penurunan kecepatan ini disebabkan oleh pelenturan dan melelehnya batuan pada temperatur tinggi serta perbedaan ekspansi termal pada mineral-mineral penyusun batuan. Selain itu, perubahan kecepatan gelombang tidak hanya dipengaruhi oleh temperatur, tetapi juga oleh tekanan dan karakteristik pori batuan. Tekanan memberikan efek yang equivalent terhadap perubahan kecepatan gelombang, dimana peningkatan tekanan akan memberikan efek pada peningkatan kecepatan gelombang seismik akibat terjadinya kompaksi batuan (Trampert, 2001). Sementara itu, keberagaman karakteristik batuan memberikan pengaruh yang bervariasi terhadap kecepatan gelombang seismik. Porositas memberikan efek yang besar terhadap perubahan kecepatan gelombang seismik, dan perubahan kecepatan akan bergantung pada variasi bentuk pori dan fluida pengisi batuan (Takei, 2002). Kecepatan gelombang S dipengaruhi fluida pengisi pori dimana modulus shear fluida adalah nol. Peningkatan densitas juga akan mempengaruhi perubahan kecepatan. Sedangkan pada gelombang P, modulus bulk memberikan efek relatif lebih besar pada perubahan fluida dan ruang pori dibandingkan perubahan densitas. Kecepatan gelombang P akan meningkat seiring peningkatan modulus bulk dari fluida pengisi pori (Tatham dan McCormack, 1991). Secara umum, fluida yang mengisi pori akan memberikan pengaruh pada kecepatan gelombang P dan kecepatan gelombang S sebagai berikut. Pada batuan steamsaturated, baik kecepatan gelombang P maupun S cenderung menurun. Pada batuan water-saturated, saturasi air meningkatkan kecepatan gelombang P relatif dibandingkan pada steam-saturated dan kecepatan gelombang S menurun. Peningkatan gelombang P terjadi karena air cenderung lebih incompressible dibandingkan udara sehingga keberadaan air dalam pori akan meningkatkan modulus bulk batuan tersebut. Sedangkan penurunan kecepatan gelombang S terjadi karena peningkatan densitas saat batuan tersaturasi oleh air (Wang dkk., 1990).
7
2.2
Penentuan Lokasi Hiposenter
Penentuan lokasi hiposenter melibatkan suatu proses inversi untuk mencari suatu lokasi hiposenter yang memiliki error minimum antara waktu tempuh observasi dengan waktu tempuh kalkulasi. Secara garis besar, proses inversi ini merupakan upaya untuk mendapatkan informasi sifat fisis bawah permukaan bumi berdasarkan respons dan informasi rekaman gempa di permukaan. Adapun pada proses tersebut, informasi atau parameter yang telah diketahui ialah waktu tiba gelombang dan posisi stasiun perekam gempa, sedangkan informasi atau parameter yang belum diketahui ialah waktu terjadi gempa, posisi sumber gempa (hiposenter), ray path, dan model kecepatan bawah permukaan. Maka karena itu, permasalahan inversi yang dihadapi ialah kasus inversi non-linier. Untuk mempermudah perhitungan, model kecepatan gelombang seismik yang digunakan merupakan model kecepatan 1-D dimana kecepatan homogen secara lateral dan hanya bervariasi tiap kedalaman. Kemudian, waktu tempuh adalah waktu yang ditempuh gelombang dari awal terjadinya gempa sampai terekam pada seismogram, ditunjukkan oleh persamaan: (2.2.1) dan waktu tempuh antara sumber i dan stasiun pengamat j dapat diformulasikan oleh persamaan: √
(2.2.2) Keterangan:
tj
= waktu tiba (arrival time) di stasiun pengamat j,
τi
= waktu terjadinya gempa (origin time) dari sumber i,
Tij
= waktu tempuh antara sumber i dan stasiun pengamat j,
xi , yi , zi = koordinat posisi dari sumber i, xj , yj , zj = koordinat posisi dari stasiun pengamat j, V
= kecepatan gelombang seismik.
8
2.2.1
Ray Tracing Metode Shooting Untuk menentukan waktu tempuh gelombang P dan S yang merambat dari sumber gempa ke stasiun penerima dilakukan perhitungan menggunakan metode ray tracing. Pada penentuan lokasi hiposenter ini, penulis mengunakan metode shooting. Metode ini menerapkan hukum Snell untuk menentukan sudut datang dan transmisi pada setiap batas lapisan kecepatan gelombang seismik, ditunjukkan oleh Gambar 2.2. Cara kerja metode ini adalah dengan cara melakukan tebakan sudut berangkat yang bisa tepat mengenai stasiun.
Gambar 2.2 Ilustrasi sudut datang dan transmisi ray pada batas lapisan.
Skema di atas mengikuti persamaan berikut:
(2.2.3) Untuk penerapan di model kecepatan gelombang seismik 1-D, perhitungan lebih praktis karena hanya melakukan iterasi terhadap sudut datang, sedangkan sudut transmisi pada setiap batas lapisan ditentukan secara langsung dengan menggunakan hukum Snell dari hasil penentuan sudut berangkat ray. Ilustrasi metode ini diperlihatkan pada Gambar 2.3.
9
Gambar 2.3 Ilustrasi iterasi sudut dari sumber ke stasiun.
2.2.2
Inversi Simulated Annealing Tujuan utama metode inversi umumnya bukan untuk mengetahui secara keseluruhan bentuk permukaan obyektif pada ruang model melainkan mencari nilai minimum. Nilai minimum yang dicari adalah minimum global dari fungsi obyektif yang berasosiasi dengan model optimum. Pada metode pencarian global, pola fungsi obyektif sebenarnya dapat diperkirakan berdasarkan harga fungsi obyektif pada beberapa sampel model yang dipilih secara acak dari ruang model. Pemilihan model pada metode ini dilakukan secara acak. Setiap model dalam ruang model memiliki peluang yang sama untuk dipilih sebagai sampel model. Bilangan acak dibangkitkan dengan probabilitas uniform antar 0 dan 1 yang kemudian dipetakan pada interval harga parameter model. Perhitungan pemodelan kedepan dilakukan untuk model yang terpilih yang jumlahnya tidak terlalu besar jika dibandingkan dengan jumlah keseluruhan model yang mungkin pada ruang model. Metode ini sering disebut sebagai metode Monte-Carlo karena mengambil analogi dengan perjudian yang umumnya bersifat acak. 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 10
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 energi sistem yang minimum. 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 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 masing-masing. 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 )
(2.2.4)
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: 11
P(∆E) = exp (
(2.2.5)
)
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 ke model 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. Hal ini bertujuan untuk menghindari konvergensi solusi ke minimum lokal. 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
(2.2.6)
dimana 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. Diagram alir dari metode simulated annealing ditunjukkan oleh Gambar 2.4. Suatu perumusan atau algoritma lain dapat dikembangkan untuk mencari solusi dari suatu permasalahan inversi yang kompleks ataupun unik. 12
Gambar 2.4 Diagram alir metode simulated annealing secara umum.
2.3
Seismik Tomografi
Tomografi didefinisikan sebagai suatu rekonstruksi sebuah model dari observasi besaran fisis yang merepresentasikan efek dari penjalaran suatu bentuk radiasi melalui benda yang diamati. Penerapan dalam penelitian ini adalah berupa pemodelan tomografi pada bidang seismik untuk mencitrakan struktur 3-D kecepatan gelombang seismik dengan menggunakan data gempa. Konsep dasar dari tomografi seismik adalah memperhitungkan data waktu tempuh gelombang (travel time). Waktu tempuh ditunjukkan oleh persamaan: (2.3.1) dan waktu tempuh antara sumber i dan stasiun pengamat j dapat diformulasikan (dalam Stewart, 1991): ∫
(2.3.2) 13
di mana Perlambatan (slowness) sebagai fungsi dari posisi mempunyai persamaan: (2.3.3) Keterangan:
tj
= waktu tiba (arrival time) di stasiun pengamat j,
τi
= waktu terjadinya gempa (origin time) dari sumber i,
Tij
= waktu tempuh antara sumber i dan stasiun pengamat j,
s(r)
= perlambatan (slowness) sebagai fungsi dari posisi,
dl
= panjang segmen sinar seismik dari lintasan integrasi Lij 3-D bergantung pada lokasi hiposenter dan penerima serta struktur 3-D bumi,
V 2.3.1
= kecepatan gelombang seismik.
Parameterisasi Model Parameterisasi yang digunakan dalam inversi tomografi umumnya adalah blok atau kotak. Penentuan besarnya blok tergantung pada banyak data yang tersedia dan kebutuhan resolusi dalam investigasi bawah permukaan. Besar blok dapat seragam dan tidak seragam. Berikut ini beberapa contoh parameterisasi model dengan skala seismik tomografi dalam bidang ilmu kebumian (Havskov dan Ottemӧller, 2010): a. Skala global dengan ukuran blok 20 – 50. b. Skala regional dengan ukuran blok 10 – 20. c. Skala lokal dengan ukuran blok 20 km. d. Skala gunungapi dengan ukuran blok 2 km. e. Skala eksplorasi dan geoteknik dengan ukuran 1-5 m.
2.3.2
Ray Tracing Metode Pseudo-Bending Untuk menentukan waktu tempuh gelombang P dan S yang merambat dari sumber gempa ke stasiun penerima dilakukan perhitungan menggunakan metode ray tracing. Pada metode ini, diperlukan suatu model awal kecepatan baik berupa model tiga dimensi atau satu dimensi berupa blok-blok baik 14
homogen maupun heterogen sebesar area penelitian. Pada penelitian ini digunakan metode ray tracing pseudo-bending (Um dan Thurber, 1987) untuk menghitung waktu tempuh kalkukasi dari sumber ke penerima dalam proses inversi tomografi. Metode pseudo-bending menggunakan prinsip Fermat dimana gelombang merambat pada lintasan dengan waktu tempuh tercepat. Waktu tempuh (T) sepanjang lintasan gelombang diekspresikan dalam sebuah persamaan integral di antara dua titik (Um dan Thurber, 1987). Dalam perhitungan waktu tempuh gelombang secara penjumlahan numerik sepanjang ray segment, persamaan waktu tempuh gelombang dapat dituliskan kembali dengan menggunakan cara trapezoidal (Um dan Thurber, 1987):
∑| ⃑
⃑
(
)
| (2.3.4)
dimana merupakan nomor titik yang mendefinisikan ray dan adalah vektor posisi titik ke-k, sedangkan merupakan kecepatan gelombang pada titik ke-k. Dalam perhitungan waktu tempuh, perlu dilakukan pembangunan matriks Kernell. Untuk itu, panjang ray di tiap blok perlu dihitung terlebih dahulu sehingga diperoleh ray trace dan waktu tempuh untuk setiap source dan receiver. Diagram alir sederhana dari pemograman ini ditunjukkan oleh Gambar 2.5.
Gambar 2.5 Diagram alir perhitungan waktu tempuh gelombang menggunakan ray tracing pseudo-bending pada penelitian ini.
15
2.3.3
Tomografi Waktu Tunda Dalam meyelesaikan suatu pemodelan inversi tomografi, kita memerlukan tomografi delay time atau waktu tunda (dt).
Pada dasarnya delay time
merupakan perbedaan antara waktu tempuh (travel time) observasi dengan waktu tempuh (travel time) kalkulasi. (2.3.5) di mana tobs merupakan waktu tempuh (travel time) observasi yang didapat dari hasil observasi di lapangan, sedangkan tcal merupakan waktu tempuh (travel time) pada suatu medium dari model kecepatan yang telah diketahui. Persamaan matriks tomografi didapat dari persamaan delay time seperti di bawah ini (Van der Hilst dan Engdahl, 1991): ∫ (2.3.6) Persamaan waktu tunda (delay time) di atas dalam bentuk model parameterisasi atau jika dalam model n elemen blok dapat dinyatakan sebagai berikut: ∑
(2.3.7)
Dalam persamaan di atas index i merupakan blok ke-i, j adalah sinar seismik ke-j, dan dl adalah panjang gelombang di dalam elemen blok. Selanjutnya persamaan (x) dapat diuraikan menjadi persamaan linier simultan untuk semua sinar gelombang ke-j sebagai berikut:
(2.3.8)
16
Persamaan 2.3.8 dapat dibuat dalam bentuk matriks sebagai berikut: [d] = [A] [x]
[
]
[
][
]
(2.3.9) di mana [d] merupakan matriks yang berisi delay time untuk setiap sinar gelombang ke-j yang direkam oleh stasiun, [A] merupakan matriks kernel yaitu matriks yang berisi panjang ray path untuk setiap blok yang dilewati gelombang dan matriks [A] mempunyai ukuran j x n, sedangkan [x] merupakan matriks yang berisi parameter deviasi slowness yang dimiliki setiap elemen blok dan matriks [x] merupakan parameter yang dicari. 2.3.4
Inversi Iterative Damped Least Square Matriks parameter [x] pada persamaan 2.3.9 dapat dicari dengan menggunakan metode inversi. Penelitian ini menggunakan inversi least square untuk mendapatkan parameter deviasi slowness. Least square adalah metode penentuan parameter fitting suatu fungsi untuk disesuaikan dengan sebaran data pengukuran sehingga memiliki eror simpangan yang terkecil dan diharapkan membentuk suatu model yang dapat menggambarkan sebaran data dengan baik (Grandis, 2009). Parameter deviasi slowness didapat dengan menyelesaikan persamaan inversi seperti berikut (Paige dan Saunders, 1982) : (2.3.10) Proses inversi tomografi menggunakan damping pada matriks tomografi agar hasil inversinya dapat diterima. Terdapat dua tipe damping eksplisit (Widiyantoro dkk., 2000) yang digunakan pada penelitian ini, yaitu: a. Norm Damping Damping ini memberikan solusi untuk blok yang tidak dilewati sinar seismik menjadi bias terhadap model awal. Matriks [x] pada 17
persamaan 2.3.11 merupakan bentuk matriks ketika diberikan norm damping, dengan A merupakan matriks Kernell, I merupakan matriks identitas, dt merupakan delay time, dan α adalah parameter norm damping. (
)
(
)
(2.3.11)
b. Gradient Damping Damping ini memberikan solusi untuk blok yang tidak dilewati sinar seismik menjadi bias terhadap model yang smooth. Pada persamaan 2.3.12, matriks [x] merupakan bentuk matriks gradient damping dimana γG adalah parameter gradient damping. (
)
(
) (2.3.12)
Kedua damping tersebut dapat digabungkan pada solusi inversi tomografi sehingga persamaannya menjadi:
(
)
( ) (2.3.13)
Setelah didapat matriks [x] yang berisi parameter deviasi slowness (ΔS), maka parameter ini nantinya akan dikonversi untuk mendapatkan parameter deviasi kecepatan (ΔV). Persamaan konversi untuk mendapatkan deviasi kecepatan adalah sebagai berikut:
S = S – So S= S= S= S=
18
V= V= V
=
V= (2.3.14) Sehingga dari konversi turunan deviasi slowness (Δs) didapatkan persamaan kecepatan V pada iterasi ke-i hasil inversi seperti persamaan di bawah ini: Vi = Vo +
= Vo (2.3.15)
Inversi dilakukan secara iteratif sehingga liniearisasi terjadi pada setiap iterasi. Diagram alir dari proses inversi yang terjadi hingga mendapat solusi model kecepatan dengan selisih waktu tempuh kalkulasi dan observasi minimum, ditunjukkan oleh Gambar 2.6.
Gambar 2.6 Diagram alir pemodelan inversi tomografi dengan metode iteratve damped least square pada penelitian ini.
19
2.3.5
Tes Resolusi Tes resolusi atau yang dikenal Checkerboard Resolution Test (CRT) merupakan metode forward modelling yang bertujuan untuk menguji kehandalan teknik inversi yang digunakan dalam inversi tomografi dan resolusi di seluruh ruang model. Tes resolusi ini dilakukan dengan cara mengalikan anomali positif dan negatif secara selang seling arah vertikal dan horisontal dengan model kecepatan awal yang akan digunakan dalam inversi tomografi. Besarnya anomali yang tergantung terhadap prediksi model pertubasi yang akan dihasilkan saat inversi tomografi. Ukuran anomali sebaiknya sama dengan ukuran dari parameterisasi model dalam inversi tomografi. Hasil forward modeling dari perkalian antara model CRT dengan model kecepatan awal tersebut kemudian menjadi data input observasi untuk inversi tomografi. Daerah yang memiliki hasil inversi baik (memiliki anomali positif dan negatif) merupakan daerah yang teresolusi dengan baik. Distribusi pola anomali tersebut akan sesuai dengan daerah yang dilewati banyak sinar gelombang.
Gambar 2.7 Diagram alir tes resolusi checkerboard pada penelitian ini.
20
BAB III PENENTUAN LOKASI HIPOSENTER
3.1
Penelitian Terdahulu
Penelitian sebelumnya telah dilakukan untuk mengidentifikasi event dan menentukan lokasi hiposenter pada lapangan tersebut menggunakan metode Geiger. Penelitian tersebut menghasilkan 591 event gempa mikro (Chevron Salak Geothermal). Data tersebut direkam oleh 19 stasiun dan terjadi selama satu tahun di lapangan geotermal “RR” dari periode Januari 2007 sampai Desember 2007. Gambar 3.1 menunjukkan persebaran hiposenter dan stasiun hasil penelitian tersebut. Menurut penelitian sebelumnya, gempa mikro pada lapangan geotermal “RR” terjadi pada rentang antara elevasi 0.75 sampai -8.43 km (MSL = 0 km). Apabila dirincikan, persebaran event dapat dibagi menjadi 4 zona sebagai berikut (Gambar 3.2 dan 3.3) : 1. Zona dangkal, terdapat 67 event gempa mikro dari total 591 event pada rentang antara elevasi 1 sampai 0 km (kotak abu-abu pada Gambar 3.2 dan 3.3). Ada kemungkinan bahwa gempa mikro tersebut terjadi terlalu dangkal. 2. Zona menengah, pada rentang antara elevasi 0 sampai -6.5 km (kotak jingga). 3. Zona dalam, pada rentang antara elevasi -6.5 sampai -9 km (kotak biru tua). 4. Zona jauh, terdapat 9 event gempa mikro dari total 591 event terjadi jauh dari persebaran event (lingkaran ungu). Setiap metode inversi memiliki kelebihan dan kekurangannya masing-masing. Solusi hiposenter tersebut bisa saja terjebak pada minimum lokal karena pada penelitian sebelumnya, metode inversi yang digunakan ialah metode least square dimana metode itu sangat bergantung pada model awal yang baik. Oleh karena itu, penulis akan menerapkan metode inversi lain yaitu metode simulated annealing, dimana solusi dari metode tersebut tidak bergantung pada model awal. Dengan asumsi bahwa picking waktu tiba gelombang pada penelitian sebelumnya sudah benar dan akurat, penelitian ini difokuskan pada metode inversi yang digunakan. 21
U
MSL
U
Gambar 3.1 Peta persebaran event yang dihasilkan oleh penelitian sebelumnya (bulat hijau) dan stasiun (segitiga biru terbalik).
22
zona dangkal
MSL
zona jauh zona menengah
zona dalam
Gambar 3.2 Peta penampang U-S terhadap kedalaman dari persebaran event (bulat hijau) dan stasiun (segitiga biru terbalik).
MSL
zona dangkal
zona menengah zona jauh
zona dalam
Gambar 3.3 Peta penampang B-T terhadap kedalaman dari persebaran event (bulat hijau) dan stasiun (segitiga biru terbalik).
23
3.2
Persiapan Data
Data yang digunakan merupakan data rekaman gempa mikro yang direkam oleh 19 stasiun pengamat dan terjadi selama satu tahun di lapangan geotermal “RR” dari periode Januari 2007 sampai Desember 2007. Secara rinci, data tersusun dari 591 event dengan 3184 fasa gelombang P dan 2760 fasa gelombang S. Keseluruhan data dari penelitian sebelumnya mencakup data waktu terjadi gempa (origin time), waktu tempuh gelombang P dan S (travel time), dan lokasi stasiun pengamat. Contoh cuplikan dari keseluruhan data tersebut untuk fasa gelombang P ditunjukkan oleh Tabel 3.1. Tabel 3.1 Contoh katalog data fasa gelombang P No. event 3 3 3 5 5 5 5 6 6 6 6 6 6 6 6 6
Origin Time hh.mm
ss
10:25 10:25 10:25 17:26 17:26 17:26 17:26 23:29 23:29 23:29 23:29 23:29 23:29 23:29 23:29 23:29
4.997 4.997 4.997 1.396 1.396 1.396 1.396 7.067 7.067 7.067 7.067 7.067 7.067 7.067 7.067 7.067
No. Stasiun ST06 ST07 ST15 ST15 ST06 ST13 ST07 ST02 ST05 ST13 ST12 ST03 ST06 ST10 ST07 ST15
hh.mm
ss
Travel Time ss
10:25 10:25 10:25 17:26 17:26 17:26 17:26 23:29 23:29 23:29 23:29 23:29 23:29 23:29 23:29 23:29
6.307 6.407 6.45 3.039 4.101 4.418 4.601 7.617 7.808 7.799 7.961 8.05 8.038 8.249 8.41 8.684
1.31 1.41 1.453 1.643 2.705 3.022 3.205 0.55 0.741 0.732 0.894 0.983 0.971 1.182 1.343 1.617
Arrival Time
Pada penelitian ini, penulis menggunakan data fasa gelombang P tersebut serta model kecepatan 1-D fasa gelombang P (Vp) dari penelitian sebelumnya (Chevron Salak Geothermal) sebagai input data. Model kecepatan 1-D dari penelitian sebelumnya ditunjukkan oleh Tabel 3.2 dan Gambar 3.4.
24
Tabel 3.2 Katalog model kecepatan 1-D Vp dan Vs (Chevron Salak Geothermal) Lapisan
Ketebalan Lapisan (km)
Elevasi Top Lapisan (km)
1 2 3 4 5 6 7 8
2 0.4 0.3 0.4 0.5 0.9 2.5 20
2 0 -0.4 -0.7 -1.1 -1.6 -2.5 -5
Kecepatan Gelombang P (km/s) 3.35 3.99 4.37 4.62 5.25 5.83 6.6 8
Kecepatan Gelombang P (km/s) 1.93 2.30 2.52 2.67 3.03 3.37 3.81 4.62
2 Vp Vs
1 0 -1
Kedalaman (km)
-2 -3 -4 -5 -6 -7 -8 -9
0
2 4 6 8 Kecepatan Gelombang Seismik (km/s)
Gambar 3.4 Plot model kecepatan untuk gelombang P (Vp) dan gelombang S (Vs).
(Chevron Salak Geothermal)
25
3.3
Perhitungan Waktu Tempuh Gelombang
Pada penelitian ini, waktu tempuh observasi (tobs) merupakan selisih antara waktu tiba (arrival time) dengan waktu terjadi gempa (origin time) dari data observasi. Sementara itu, waktu tempuh kalkulasi (tcal) merupakan hasil ray tracing dengan metode shooting. Ray tracing tersebut menggunakan model kecepatan 1-D untuk menentukan jejak sinar yang dilewati gelombang dari sumber gempa ke stasiun penerima, kemudian menghitung waktu tempuh gelombang. 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. Gambar 2.3 menunjukkan algoritma perhitungan waktu tempuh gelombang dan Gambar 2.4. menunjukkan skema perhitungan waktu tempuh total dari suatu sumber ke suatu stasiun.
Gambar 3.5 Diagram alir perhitungan waktu tempuh gelombang (metode shooting) pada penelitian ini. 26
Gambar 3.6 Skema perhitungan waktu tempuh gelombang total dari sumber (bintang) ke stasiun penerima (segitiga hitam terbalik). [V1, V2, …., V7] adalah kecepatan lapisan ke-i, [dl1, dl2, …., dl7] adalah panjang ray di lapisan ke-i, dan [t1, t2, …., t7] adalah waktu tempuh gelombang di lapisan ke-i. (i = 1, 2, …., 7).
3.4
Algoritma Inversi
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. Dalam proses inversi ini, 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 27
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 bisa disebut dengan perturbasi model. Berdasarkan Vakil-Baghmishes dan Navarbaf (2008), parameter model akan diperbarui sesuai dengan persamaan: (
)
]
(3.4.1) sehingga pada simulasi ke-k akan dihasilkan dk yang bernilai [-1,1] dan besar pembaruan model di simulasi ke-k adalah: (3.4.2) Pada persamaan 3.4.2 , 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 Ak menggambarkan 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 28
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 U akan 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 Umax x Lmax kali. Berdasarkan algoritma di atas, parameter annealing yang digunakan pada penelitian ini adalah: Amin = 10, Lmax = 15, dan Umax = 20. Gambar 2.5 menggambarkan algoritma keseluruhan prosedur inversi yang dilakukan sampai memperoleh solusi.
Gambar 3.7 Diagram alir metode inversi simulated annealing pada penentuan lokasi hiposenter penelitian ini. 29
3.5
Pengolahan 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 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 3.3 Tabel 3.3 Katalog data sintetis waktu tempuh gelombang Posisi Hiposenter X (km) Y (km) Z (km)
11.0
11.0
-3.0
No. 1 2 3 4 5
Posisi Stasiun X (km) Y (km) 16.61666 9.3853 16.51366 11.5613 12.29366 11.4463 14.31366 9.6153 14.00566 8.2193
Z (km) 1.335 1.42 1.05 1.22 1.097
Waktu Tempuh (s) 1.561148715 1.540462621 0.943305874 1.211803692 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 3.4 dan Gambar 3.8. Tabel 3.4 Katalog solusi hasil inversi pada data sintetis Percobaan ke1 2 3 4 5 6 7 8 9 10
Posisi Hiposenter X (km) Y (km) Z (km) 11.03191 10.99055 -2.99346 11.02536 11.04477 -3.01915 11.06267 11.06747 -3.0462 11.01087 11.03581 -2.98444 11.05361 11.01492 -3.0276 11.05611 11.01796 -3.02093 11.01033 11.00713 -3.00043 11.02544 10.99621 -2.99265 11.01741 11.00869 -3.02176 11.03251 11.00522 -3.02504
error RMS (s) 0.002495 0.003215 0.002582 0.003117 0.002157 0.001617 0.002506 0.002915 0.002446 0.001925 30
11 12 13 14 15 16 17 18 19 20
11.0046 11.00836 11.00018 11.04784 11.00541 11.0129 11.0188 10.97239 11.0188 11.04076
10.98935 10.95518 10.98141 10.97016 11.00144 10.9679 10.97034 10.98823 11.02949 10.98443
-2.98496 -3.0037 -2.98632 -3.06595 -2.97761 -3.00386 -3.01281 -2.9635 -2.99187 -3.04058
0.00305 0.002366 0.00247 0.003328 0.002881 0.002379 0.002822 0.002069 0.002202 0.002941
Gambar 3.8 Peta persebaran hiposenter hasil inversi (lingkaran merah), hiposenter sebenarnya (bulat hijau), dan stasiun pengamat (segitiga biru terbalik).
x 10-3
Gambar 3.9 Histogram error RMS dari 20 kali perhitungan. 31
Dapat dilihat algoritma inversi simulated annealing 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.
3.6
Analisis dan Pembahasan
Berikut ialah 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 (Gambar 3.9 dan 3.10) : 1. Zona dangkal, berada pada rentang antara elevasi 1 sampai 0 km (kotak abuabu). Terdapat 18 event gempa mikro. 2. Zona menengah, berada pada rentang antara elevasi 0 sampai -6.5 km (kotak jingga). Terdapat 556 event gempa mikro. 3. Zona dalam, berada pada rentang antara elevasi -6.5 sampai -9 km (kotak biru tua). Terdapat 9 event gempa mikro. 4. Zona jauh, berada pada rentang yang jauh dari persebaran event lain (lingkaran ungu). Terdapat 9 event gempa mikro. 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. Gempa-gempa 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. 32
U
MSL
U
Gambar 3.10 Peta persebaran event yang dihasilkan oleh penelitian ini menggunakan simulated annealing (bulat merah) dan stasiun (segitiga biru terbalik).
33
zona dangkal
MSL
zona menengah
zona jauh
zona dalam
Gambar 3.11 Peta penampang U-S terhadap kedalaman dari persebaran event (bulat merah) dan stasiun (segitiga biru terbalik).
MSL
zona dangkal
zona menengah
zona jauh zona dalam
Gambar 3.12 Peta penampang B-T terhadap kedalaman dari persebaran event (bulat merah) dan stasiun (segitiga biru terbalik).
34
U
MSL
U
Gambar 3.13 Peta persebaran event yang dihasilkan oleh penelitian sebelumnya (bulat hijau) dan penelitian ini menggunakan simulated annealing (bulat merah) beserta stasiun (segitiga biru terbalik). Garis biru muda menunjukkan perubahan posisi hiposenter pada event yang sama.
35
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 3.14. Dapat disimpulkan bahwa metode simulated annealing menghasilkan solusi yang lebih baik secara statistik.
a)
b)
Gambar 3.14 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).
Persebaran gempa mikro pada lapangan panas bumi “RR” akan berhubungan dengan struktur 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 pemodelan tomografi seismik 3-D untuk mencitrakan struktur bawah permukaan pada lapangan panas bumi “RR”. 36
BAB IV SEISMIK TOMOGRAFI 3-D
4.1
Persiapan Data
Data yang digunakan merupakan data rekaman gempa mikro yang direkam oleh 19 stasiun pengamat dan terjadi selama satu tahun di lapangan geotermal “RR” dari periode Januari 2007 sampai Desember 2007. Data hiposenter yang digunakan ialah hasil inversi dengan metode simulated annealing pada penelitian ini karena solusi yang dihasilkan lebih baik secara statistik. Secara rinci, data tersusun dari 591 event dengan 3184 fasa gelombang P dan 2760 fasa gelombang S. Keseluruhan data dari penelitian sebelumnya mencakup data lokasi hiposenter, waktu terjadi gempa (origin time), waktu tempuh gelombang P dan S (travel time), dan lokasi stasiun pengamat. Persebaran event dapat dilihat pada Gambar 3.10. Data hiposenter tersebut kemudian disaring untuk meningkatkan akurasi perhitungan. Hiposenter yang digunakan ialah hiposenter yang berlokasi pada rentang area 4 – 23 km Barat-Timur, 5 – 17 km Utara-Selatan, dan kedalaman maksimum adalah -6.5 km dari MSL. Penyaringan data tersebut ditunjukkan pada Gambar 4.1. Contoh cuplikan data yang digunakan sebagai input ditunjukkan oleh Tabel 4.1 dan Tabel 4.2. U MSL
Gambar 4.1 Distribusi event (bulat merah) dan area data hiposenter yang digunakan (kotak hijau) pada pemodelan tomografi. 37
Tabel 4.1 Contoh format katalog data fasa gelombang P No. event 10 10 10 10 11 11 11 11 11 12 12 12 12 12 13 13 13
Posisi Hiposenter X (km) Y (km) Z (km) 11.8927 11.6264 -2.28777 11.8927 11.6264 -2.28777 11.8927 11.6264 -2.28777 11.8927 11.6264 -2.28777 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.80351 6.3916 0.072705 15.80351 6.3916 0.072705 15.80351 6.3916 0.072705
Posisi Stasiun X (km) Y (km) 12.44266 11.4643 14.49666 9.5183 9.28366 9.2213 11.32566 8.0413 12.44266 11.4643 11.32566 8.0413 14.49666 9.5183 9.28366 9.2213 16.52066 11.5663 14.84966 12.2513 16.52066 11.5663 12.44266 11.4643 14.49666 9.5183 11.32566 8.0413 14.49666 9.5183 16.52066 11.5663 9.28366 9.2213
Z (km) 1.087 1.298 1.04 0.985 1.087 0.985 1.298 1.04 1.473 1.285 1.473 1.087 1.298 0.985 1.298 1.473 1.04
Waktu Tempuh (s) 0.79 1.199 1.201 1.014 0.507 0.714 0.918 0.94 1.429 0.811 1.051 1.185 1.356 1.727 1.047 1.632 2.155
Tabel 4.2 Contoh format katalog data fasa gelombang S No. event 10 10 10 10 11 11 11 11 11 12 12 12 12 12 13 13
Posisi Hiposenter X (km) Y (km) Z (km) 11.8927 11.6264 -2.28777 11.8927 11.6264 -2.28777 11.8927 11.6264 -2.28777 11.8927 11.6264 -2.28777 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 11.94034 10.331 -0.08764 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.40857 14.0998 -1.239 15.80351 6.3916 0.072705 15.80351 6.3916 0.072705
Posisi Stasiun X (km) Y (km) 12.44266 11.4643 14.49666 9.5183 9.28366 9.2213 11.32566 8.0413 12.44266 11.4643 11.32566 8.0413 14.49666 9.5183 9.28366 9.2213 16.52066 11.5663 14.84966 12.2513 16.52066 11.5663 12.44266 11.4643 14.49666 9.5183 11.32566 8.0413 14.49666 9.5183 16.52066 11.5663
Z (km) 1.087 1.298 1.04 0.985 1.087 0.985 1.298 1.04 1.473 1.285 1.473 1.087 1.298 0.985 1.298 1.473
Waktu Tempuh (s) 1.313 2.04 2.128 1.914 1.205 1.481 1.741 1.784 2.617 1.306 1.777 2.038 2.186 3.024 2.018 2.942
38
Setelah disaring, jumlah event gempa mikro yang diperoleh ialah 573 event dan menghasilkan 3099 raypath fasa gelombang P dan 2499 raypath fasa gelombang S. Kemudian inversi membutuhkan model kecepatan awal pada pengolahan data. Model kecepatan awal yang digunakan adalah model kecepatan 1-D fasa gelombang P (Vp) dan fasa gelombang S (Vs) dari penelitian sebelumnya. Model kecepatan 1-D yang digunakan ditunjukkan oleh Tabel 3.2 dan Gambar 3.4.
4.2
Parameterisasi Model
Pada penelitian ini dilakukan parameterisasi model blok tiga dimensi. Penentuan jumlah dan besarnya tiap blok model ini bergantung pada luas area dan kedalaman daerah penelitian, serta distribusi data. Daerah penelitian memiliki luas area 15 x 10 x 8 km3 pada rentang area 5 – 20 km Barat-Timur, 5 – 15 km Utara-Selatan, dan elevasi 1.5 sampai -6.5 km (MSL = 0 km). Parameterisasi model dibangun dengan dimensi berjumlah 10 x 10 x 16 blok (1600 blok) dan ukuran dimensi blok 1.5 x 1 x 0.5 km3.
U
Gambar 4.2 Area penelitian dan parameterisasi model blok tiga dimensi pada penelitian ini.
39
4.3
Perhitungan Waktu Tempuh Gelombang
Perhitungan waktu tempuh gelombang dilakukan terhadap data distribusi event gempa mikro dan stasiun penerima. Waktu tempuh observasi (tobs) merupakan selisih antara waktu tiba (arrival time) dengan waktu terjadi gempa (origin time) dari data observasi. Sementara itu, waktu tempuh kalkulasi (tcal) merupakan hasil ray tracing dengan metode pseudo-bending (Um dan Thurber, 1987). Ray tracing tersebut menggunakan model kecepatan 3-D untuk menentukan jejak sinar yang dilewati gelombang dari sumber gempa ke stasiun penerima, kemudian menghitung waktu tempuh gelombang dari jejak sinar tersebut.
U
U
MSL
Gambar 4.3
Peta distribusi cakupan ray path (garis biru), hiposenter (bulat merah), dan stasiun (segitiga biru terbalik) hasil ray tracing metode pseudo-bending melalui model kecepatan 3-D gelombang P.
40
Ray tracing pada penelitian ini menggunakan kode program MATLAB yang telah dibuat dan dikembangkan sebelumnya oleh Syahputra (2011),
kemudian
dimodifikasi dan disesuaikan dengan penelitian ini. Ray tracing ini menghasilkan panjang ray tiap segmen dan panjang ray keseluruhan untuk setiap perambatan gelombang P dan S dari hiposenter ke stasiun. Penjalaran gelombang P dan S dari hiposenter menuju stasiun ditunjukkan oleh Gambar 4.2 dan Gambar 4.3.
U
U
MSL
Gambar 4.4
Peta distribusi cakupan ray path (garis hijau), hiposenter (bulat merah), dan stasiun (segitiga biru terbalik) hasil ray tracing metode pseudo-bending melalui model kecepatan 3-D gelombang S.
41
4.4
Algoritma Inversi
4.4.1
Algoritma Iterative Damped Least Square Berdasarkan Widiyantoro dkk. (2000), dengan menggunakan data waktu tunda (selisih antara waktu tempuh observasi dan waktu tempuh kalkukasi) dan panjang ray path tiap segmen model tiga-dimensi, dapat dibangun matriks tomografi seperti pada persamaan (2.3.8) – (2.3.9). Untuk menghindari nilai determinan matriks sama dengan nol dan memberikan solusi untuk blok yang tidak dilewati sinar seismik, digunakan norm damping (α) dan gradient damping (γ). Matriks yang terbentuk ditunjukkan oleh persamaan 2.3.13. Pada penelitian ini, nilai norm damping (α) yang digunakan adalah 10, sedangkan nilai gradient damping (γ) yang digunakan adalah 1. Selanjutnya, setelah matriks tomografi dibangun, proses inversi least square dapat dilakukan. Inversi dilakukan untuk matriks [A] terhadap matriks [dt] sehingga akan didapatkan matriks [x] yang merupakan nilai perubahan dari parameter slowness (ΔS). Proses ini dilakukan secara iteratif. Nilai perubahan kecepatan pada suatu iterasi ditunjukkan oleh persamaan (2.3.15). Algoritma keseluruhan proses inversi ditunjukkan oleh Gambar 2.6. Kemudian, pada penelitian ini, iterasi maksimum yang digunakan adalah 12. Selain itu, iterasi akan berhenti jika nilai selisih antara RMS pada iterasi ke i dan RMS pada iterasi sebelumnya lebih kecil daripada 0.001.
4.4.2
Algoritma Simulated Annealing Secara garis besar, algoritma inversi simulated annealing yang digunakan sama seperti pada penentuan lokasi hiposenter (bab 3.4) berdasarkan White (1984) dan Weber (2000). Hal yang berbeda ialah parameter annealing yang digunakan karena parameter model yang dicari jauh lebih banyak.. Pada
42
pemodelan tomografi ini, parameter annealing yang digunakan adalah: Amin = 20, Lmax = 40, dan Umax = 100. Untuk menghemat proses inversi, ray tracing hanya dilakukan pada saat penurunan temperatur. Kemudian, blok model yang tidak dilewati sinar seismik dibuat bias terhadap model awal dan blok model yang hanya dilewati sedikit sinar seismik diberikan bobot perturbasi kecil. Gambar 4.5 menunjukkan prosedur inversi yang dilakukan sampai memperoleh solusi model kecepatan.
Gambar 4.5 Diagram alir metode inversi simulated annealing pada pemodelan tomografi.
4.5
Tes Resolusi
Pada penelitian ini, tes resolusi atau biasa dikenal Checkerboard Resolution Test (CRT) akan dilakukan pada konfigurasi data hiposenter dan stasiun yang sama. Model checkerboard didapat dengan mengalikan anomali +10% dan -10% terhadap
43
model kecepatan awal (model kecepatan 1-D). Model ini memiliki dimensi jumlah blok 10 x 10 x 16.
Gambar 4.6 Model awal Checkerboard Resolution Test.
Kemudian dilakukan proses inversi sesuai dengan algoritma inversi yang telah dibuat. Metode inversi yang digunakan yaitu metode iterative damped least square dan simulated annealing. 4.5.1
CRT Iterative Damped Least Square Hasil tes resolusi menggunakan metode iterative damped least square ditunjukkan oleh Gambar 4.7 sampai 4.9. Kemudian, daerah yang teresolusi baik diinterpretasikan berada pada daerah yang memiliki anomali positif dan negatif dengan jelas, ditandai dengan garis hitam putus-putus. Pada proses inversi, iterasi berhenti dengan error RMS sebesar 0.0104 detik. Histogram persebaran delay time ditunjukkan oleh Gambar 4.10.
44
Gambar 4.7 Tomogram slice horizontal hasil dari tes resolusi checkerboard (Iterative Damped Least Square) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali positif ditunjukkan oleh warna biru dan anomali negatif ditunjukkan oleh warna merah. 45
U
Gambar 4.8 Tomogram slice vertikal Barat - Timur hasil dari tes resolusi checkerboard (Iterative Damped Least Square) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali positif ditunjukkan oleh warna biru dan anomali negatif ditunjukkan oleh warna merah. 46
U
Gambar 4.9 Tomogram slice vertikal Utara - Selatan hasil dari tes resolusi checkerboard (Iterative Damped Least Square) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali positif ditunjukkan oleh warna biru dan anomali negatif ditunjukkan oleh warna merah.
47
RMS error = 0.0104 s
Gambar 4.10 Histogram distribusi delay time dari CRT iterative damped least square.
4.5.2
CRT Simulated Annealing Pada proses inversi, iterasi berhenti dengan error RMS sebesar 0.0104 detik. Histogram persebaran delay time ditunjukkan oleh Gambar 4.11.
RMS error = 0.0153 s
Gambar 4.11 Histogram distribusi delay time dari CRT simulated annealing.
Hasil tes resolusi menggunakan metode simulated annealing ditunjukkan oleh Gambar 4.12 sampai 4.14. Kemudian, daerah yang teresolusi baik diinterpretasikan berada pada daerah yang memiliki anomali positif dan negatif dengan jelas, ditandai dengan garis hitam putus-putus.
48
Gambar 4.12 Tomogram slice horizontal hasil dari tes resolusi checkerboard (Simulated Annealing) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali positif ditunjukkan oleh warna biru dan anomali negatif ditunjukkan oleh warna merah. 49
U
Gambar 4.13 Tomogram slice vertikal Barat - Timur hasil dari tes resolusi checkerboard (Simulated Annealing) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali positif ditunjukkan oleh warna biru dan anomali negatif ditunjukkan oleh warna merah. 50
U
Gambar 4.14 Tomogram slice vertikal Utara - Selatan dari hasil tes resolusi checkerboard (Simulated Annealing) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali positif ditunjukkan oleh warna biru dan anomali negatif ditunjukkan oleh warna merah.
51
4.6
Pemodelan Inversi Tomografi
Pemodelan inversi tomografi diterapkan pada data rekaman gempa mikro di lapangan geotermal “RR” menggunakan algoritma inversi yang telah dibuat. Metode inversi yang digunakan ada 2 yaitu metode iterative damped least square dan simulated annealing. 4.6.1
Aplikasi Iterative Damped Least Square Pemodelan inversi tomografi pada data waktu tempuh gelombang P menghasilkan struktur kecepatan gelombang P (Vp) bawah permukaan lapangan geotermal “RR” yang ditunjukkan oleh tomogram struktur Vp pada Gambar 4.15 sampai 4.17. Gambar 4.15 merupakan tomogram penampang horizontal, sedangkan Gambar 4.16 dan 4.17 merupakan tomogram penampang vertikal. Pada proses inversi, iterasi berhenti pada iterasi ke-4 dengan nilai RMS error sebesar 0.0405 detik. Kemudian pemodelan tomografi juga dilakukan untuk gelombang S. Pemodelan inversi tomografi pada data waktu tempuh gelombang S menghasilkan struktur kecepatan gelombang S (Vs) bawah permukaan lapangan geotermal “RR” yang ditunjukkan oleh oleh tomogram struktur Vs pada Gambar 4.18 sampai 4.20. Gambar 4.18 merupakan tomogram penampang horizontal, sedangkan Gambar 4.19 dan 4.20 merupakan tomogram penampang vertikal. Pada proses inversi, iterasi berhenti pada iterasi ke-6 dengan nilai RMS error sebesar 0.0987 detik.
52
Gambar 4.15 Tomogram slice horizontal dari struktur kecepatan gelombang P (Iterative Damped Least Square) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 53
U
Gambar 4.16 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang P (Iterative Damped Least Square) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 54
U
Gambar 4.17 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang P (Iterative Damped Least Square) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah.
55
Gambar 4.18 Tomogram slice horizontal dari struktur kecepatan gelombang S (Iterative Damped Least Square) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 56
U
Gambar 4.19 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang S (Iterative Damped Least Square) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 57
U
Gambar 4.20 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang S (Iterative Damped Least Square) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah.
58
4.6.2
Aplikasi Simulated Annealing Proses
inversi
dengan
menggunakan
metode
simulated
annealing
membutuhkan waktu yang relatif jauh lebih lama dibandingkan dengan menggunakan metode iterative least square. Hal ini disebabkan banyaknya iterasi yang dilakukan untuk menguji setiap parameter acak yang dihasilkan. Pemodelan inversi tomografi pada data waktu tempuh gelombang P menghasilkan struktur kecepatan gelombang P (Vp) bawah permukaan lapangan geotermal “RR” yang ditunjukkan oleh tomogram struktur Vp pada Gambar 4.21 sampai 4.23. Gambar 4.21 merupakan tomogram penampang horizontal, sedangkan Gambar 4.22 dan 4.23 merupakan tomogram penampang vertikal. Pada proses inversi tomografi ini, iterasi berhenti dengan nilai RMS error sebesar 0.0425 detik. Kemudian pemodelan tomografi juga dilakukan untuk gelombang S. Pemodelan inversi tomografi pada data waktu tempuh gelombang S menghasilkan struktur kecepatan gelombang S (Vs) bawah permukaan lapangan geotermal “RR” yang ditunjukkan oleh oleh tomogram struktur Vs pada Gambar 4.24 sampai 4.26. Gambar 4.24 merupakan tomogram penampang horizontal, sedangkan Gambar 4.25 dan 4.26 merupakan tomogram penampang vertikal. Pada proses inversi ini, iterasi berhenti dengan nilai RMS error sebesar 0.1081 detik.
59
Gambar 4.21 Tomogram slice horizontal dari struktur kecepatan gelombang P (Simulated Annealing) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 60
U
Gambar 4.22 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang P (Simulated Annealing) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 61
U
Gambar 4.23 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang P (Simulated Annealing) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah.
62
Gambar 4.24 Tomogram slice horizontal dari struktur kecepatan gelombang S (Simulated Annealing) pada z = 0, 1, -1. -2, -3, -4, -5, dan -6 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 63
U
Gambar 4.25 Tomogram slice vertikal Barat - Timur dari struktur kecepatan gelombang S (Simulated Annealing) pada B-T 8, 9.5, 11, 12.5, 14, 15.5, dan 17 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 64
U
Gambar 4.26 Tomogram slice vertikal Utara - Selatan dari struktur kecepatan gelombang S (Simulated Annealing) pada U-S 7, 8, 9, 10, 11, 12, dan 13 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah.
65
4.7
Analisis dan Pembahasan
4.7.1
Perbandingan Solusi Model Kecepatan Kedua Metode Inversi Proses
inversi
dengan
menggunakan
metode
simulated
annealing
membutuhkan waktu yang relatif jauh lebih lama dibandingkan dengan menggunakan metode iterative least square. Pada komputer dengan processor AMD Athlon II X2 250, perbandingannya adalah sebagai berikut. Proses inversi menggunakan metode iterative least square untuk gelombang P membutuhkan waktu ± 20 menit, sedangkan untuk gelombang S membutuhkan waktu ± 30 menit. Proses inversi menggunakan metode simulated annealing untuk gelombang P membutuhkan waktu ± 18 jam, sedangkan untuk gelombang S membutuhkan waktu ± 19 jam. Pemodelan tomografi dengan menggunakan metode inversi iterative damped least square dan simulated annealing menghasilkan tomogram dengan pola keseluruhan yang cukup mirip. Anomali kecepatan rendah muncul pada elevasi 1 sampai -1 km, sedangkan anomali kecepatan tinggi lebih mendominasi pada elevasi -1 sampai -3 km. Contoh pola ini dapat dilihat pada beberapa tomogram slice vertikal Barat – Timur yang ditunjukkan oleh Gambar 4.27. Namun secara rinci, pada tomogram simulated annealing, terdapat beberapa anomali kecepatan yang muncul secara acak atau tidak berpola. Beberapa blok model muncul dengan nilai anomali kecepatan yang tidak menerus terhadap sekitarnya. Jika dibandingkan dengan tomogram iterative damped least square, solusi model kecepatan memiliki pola menerus terhadap sekitarnya. Hal ini dapat dilihat pada Gambar 4.28.
66
Vp
(a) Iterative Damped Least Square
Vs
(b) Simulated Annealing (c) Iterative Damped Least Square
(d) Simulated Annealing
Gambar 4.27 Cuplikan tomogram slice vertikal Barat – Timur dari: (a) struktur Vp untuk iterative damped least square, (b) struktur Vp untuk simualated annealing, (c) struktur Vs untuk iterative damped least square, (d) struktur Vs untuk simulated annealing. Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 67
Iterative Damped Least Square (a)
Simulated Annealing (b)
Gambar 4.28 Cuplikan tomogram inversi iterative damped least square (a) dan simulated annealing (b). Lingkaran ungu menunjukkan blok model yang tidak berpola. Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah.
Solusi model kecepatan dari inversi simulated annealing bersifat lebih acak (tidak berpola) dibandingkan solusi model kecepatan dari inversi iterative damped least square. Hal ini bisa disebabkan karena inversi simulated annealing melibatkan proses pencarian acak untuk mencari solusi parameter model. Parameter model yang dicari secara acak tersebut pada inversi tomografi terbilang sangat banyak, yaitu 1600 parameter model kecepatan. Ada kemungkinan bahwa solusi tersebut belum berhasil mencapai minimum global dari fungsi misfit-nya. Perbandingan dan perbedaan dari kedua solusi dapat dilihat melalui beberapa cuplikan tomogram horizontal dan vertikal dari struktur Vp dan Vs seperti yang ditunjukkan oleh Gambar 4.29 dan 4.30. 68
Struktur kecepatan Gelombang P (%Vp)
Struktur kecepatan Gelombang S (%Vs)
Gambar 4.29 Perbandingan tomogram horizontal untuk iterative damped least square (a), simulated annealing (b), dan selisih absolut dari kedua solusi (c).
69
Struktur kecepatan Gelombang P (%Vp)
Struktur kecepatan Gelombang S (%Vs)
Gambar 4.30 Perbandingan tomogram vertikal B-T untuk iterative damped least square (a), simulated annealing (b), dan selisih absolut dari kedua solusi (c).
70
Penulis menyimpulkan bahwa pada pemodelan tomografi ini, solusi model kecepatan yang dihasilkan oleh metode iterative damped least square lebih baik daripada solusi yang dihasilkan oleh metode simulated annealing. Hal ini juga didukung oleh sebaran nilai misfit kedua solusi, dimana solusi model kecepatan yang dihasilkan oleh metode iterative damped least square memiliki sebaran nilai delay time dan RMS error yang lebih kecil sehingga solusi tersebut terbilang lebih baik secara statistik.
RMS error = 0.0405 s
RMS error = 0.0425 s
Gambar 4.31 Histrogram distribusi dt (delay time) untuk solusi model kecepatan Vp dari inversi menggunakan iterative damped least square (kiri) dan simulated annealing (kanan).
RMS error = 0.0987 s
RMS error = 0.1081 s
Gambar 4.32 Histrogram distribusi dt (delay time) untuk solusi model kecepatan Vs dari inversi menggunakan iterative damped least square (kiri) dan simulated annealing (kanan). 71
4.7.2
Analisis Struktur Kecepatan Vp, Vs, dan rasio Vp/Vs Pada struktur kecepatan gelombang P, anomali kecepatan rendah muncul pada elevasi 1 sampai -1 km. Zona tersebut juga memiliki anomali kecepatan rendah pada struktur kecepatan gelombang S. Anomali kecepatan rendah ini berkisar pada -5% sampai -10%. Sedangkan pada rasio Vp/Vs, zona tersebut memiliki nilai yang relatif tinggi, yaitu berkisar pada 1.8 sampai 2. Pada struktur kecepatan gelombang P dan gelombang S, anomali kecepatan tinggi lebih mendominasi pada elevasi -1 sampai -3 km. Anomali kecepatan tinggi ini berkisar pada 5% sampai 10%. Sedangkan pada rasio vp/vs, zona tersebut memiliki nilai yang relatif rendah, yaitu berkisar pada 1.45 sampai 1.65.
72
Vp
Vs
Vp Vs
Gambar 4.33 Tomogram slice horizontal dari struktur Vp, struktur Vs, rasio Vp/Vs pada z = 1, 0, -1, dan -2 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 73
Vp
Vs
Vp Vs
Gambar 4.34 Tomogram slice horizontal dari struktur Vp, struktur Vs, rasio Vp/Vs pada z = -3, -4, -5, dan -6 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 74
U
Vp
Vs
Vp Vs
Gambar 4.35 Tomogram slice vertikal Barat - Timur dari struktur Vp, struktur Vs, rasio Vp/Vs pada B-T 8, 9.5, 11, dan 12.5 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 75
U
Vp
Vs
Vp Vs
Gambar 4.36 Tomogram slice vertikal Barat - Timur dari struktur Vp, struktur Vs, rasio Vp/Vs pada B-T 14, 15.5, 17, dan 18.5 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 76
U
Vp
Vs
Vp Vs
Gambar 4.37 Tomogram slice vertikal Utara - Selatan dari struktur Vp, struktur Vs, rasio Vp/Vs pada U-S 7, 8, 9, dan 10 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 77
U
Vp
Vs
Vp Vs
Gambar 4.38 Tomogram slice vertikal Utara - Selatan dari struktur Vp, struktur Vs, rasio Vp/Vs pada U-S 11, 12, 13, dan 14 km beserta daerah yang teresolusi baik (garis hitam putus-putus). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 78
Pada sistem panas bumi, keberadaan zona dengan temperatur tinggi pada lapisan bawah permukaan memberikan pengaruh yang bervariasi pada nilai Vp dan Vs. Pada keadaan water-saturated rock, nilai kecepatan gelombang P akan cenderung lebih tinggi dibandingkan nilai kecepatan gelombang S sehingga nilai rasio Vp/Vs relatif lebih tinggi. Sedangkan pada keadaan steam-saturated rock, nilai kecepatan gelombang P akan cenderung lebih rendah dibandingkan nilai kecepatan gelombang S sehingga nilai rasio Vp/Vs relatif lebih rendah. (Wang dkk., 1990 dan Takei, 2002). Berdasarkan penjelasan tersebut, maka pada elevasi 1 sampai -1 km, zona dengan nilai rasio Vp/Vs relatif tinggi (berkisar 1.8 – 2) dapat diinterpretasikan sebagai zona yang berasosiasi dengan water-saturated rock. Zona ini berada pada kedalaman yang dangkal, kemungkinan merupakan zona kondensasi dari uap reservoir pada ground water. Selain itu, bisa juga merupakan zona dari batuan penutup (cap rock), dimana umumnya lempung teralterasi mengandung banyak air. Pada elevasi -1 sampai -3 km, zona dengan nilai rasio Vp/Vs relatif rendah (berkisar 1.45 – 1.65) dapat diinterpretasikan sebagai zona yang berasosiasi dengan steam-saturated rock. Zona ini dapat diidentifikasi sebagai reservoir lapangan geotermal “RR”. Lapisan atau zona reservoir ini berada pada 12 – 18 km arah Barat-Timur, 8 – 12 km arah Utara-Selatan, dan kedalaman 1 – 3 km dari MSL. Studi sebelumnya menyatakan bahwa reservoir pada lapangan geotermal “RR” merupakan reservoir dominasi air (Stimac dkk., 2008). Hal ini menjadi menarik ketika pencitraan tomografi mengidentifikasi suatu zona reservoir berupa steam-saturated rock. Gunasekara dkk. (2003) menyebutkan bahwa eksploitasi dan produksi pada lapangan geotermal dapat menyebabkan perubahan sistem fasa pada reservoir. Penipisan pada pori fluida menyebabkan air digantikan oleh uap. Selain itu, penurunan tekanan pada reservoir menyebabkan penurunan titik didih, sehingga terjadi boiling dan
79
fasa uap terbentuk. Namun, apakah perubahan sistem fasa terjadi pada reservoir lapangan geotermal “RR”, hal ini belum dapat disimpulkan. Area reservoir berserta well-trajectory ditunjukkan oleh Gambar 4.39. Batas reservoir diperoleh dari Stimac dkk. (2008). Kemudian dilakukan slice pada model tomografi di sekitar area reservoir untuk melihat lebih jelas keadaan bawah permukaan. Keberadaan zona reservoir yang telah diidentifikasi sebelumnya, didukung oleh keberadaan well-trajectory pada zona anomali Vp/Vs rendah.
Gambar 4.39 Peta persebaran event (bulat merah), stasiun (segitiga biru terbalik), welltrajectory (garis hitam), dan batas reservoir (garis hitam putus-putus).
Tomogram penampang vertikal C-C’, D-D’, E-E’, dan F-F di sekitar area reservoir beserta well-trajectory ditunjukan oleh Gambar 4.40.
80
Vp
Vs
Vp/Vs
(a)
(b)
(c)
Gambar 4.40 Tomogram slice vertikal Utara - Selatan dari: (a) struktur Vp, (b) struktur Vs, (c) rasio Vp/Vs pada U-S 9, 10, 11, dan 12 km beserta well-trajectory (garis hitam). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah.
81
Zona dengan nilai Vp/Vs tinggi dan rendah terlihat cukup jelas pada tomogram penampang vertikal U-S 10 km. Kemudian dari tomogram ini, zona kondensasi uap, cap rock, dan reservoir dapat diinterpretasi seperti yang ditunjukkan pada Gambar 4.41. zona kondensasi D
cap rock Injection Production Well Well
D’
zona reservoir
Gambar 4.41 Tomogram slice vertikal dari rasio Vp/Vs pada U-S 10 km, beserta interpertasi zona kondensasi uap (garis putus-putus abu-abu), cap rock (garis putus-putus biru tua), zona reservoir (garis putus-putus merah), dan welltrajectory (garis hitam). Anomali tinggi ditunjukkan oleh warna biru dan anomali rendah ditunjukkan oleh warna merah. 82
Melalui Gambar 4.41, dapat dilihat bahwa pada penampang U-S km, zona dengan nilai Vp/Vs tinggi berada di sekitar sumur injeksi, sedangkan zona dengan nilai Vp/Vs rendah berada di sekitar sumur produksi. Bentuk volume dari zona dengan nilai Vp/Vs rendah dan tinggi di sekitar area reservoir ditunjukkan oleh Gambar 4.42 sampai 4.44. Secara umum, dapat dilihat bahwa sumur produksi pada lapangan panas bumi “RR” berada di sekitar zona Vp/Vs rendah.
U
Gambar 4.42 Model volumetri zona Vp/Vs rendah (volume merah) dan tinggi (volume biru) beserta well-trajectory dari sumur produksi (garis hitam) dan sumur injeksi (garis magenta).
83
U
Gambar 4.43 Model volumetri zona Vp/Vs rendah (volume merah) beserta well-trajectory dari sumur produksi (garis hitam).
U
Gambar 4.44 Model volumetri zona Vp/Vs tinggi (volume rendah) beserta well-trajectory dari sumur injeksi (garis magenta).
84
BAB V PENUTUP
5.1
Kesimpulan
Pada penelitian ini, melalui penerapan metode inversi least square dan simulated annealing untuk penentuan lokasi hiposenter dan pemodelan seismik tomografi 3-D pada lapangan geotermal “RR”, dapat disimpulkan bahwa: 1. Metode inversi simulated annealing dapat diterapkan pada kasus penentuan lokasi hiposenter, dimana metode ini tidak bergantung pada model awal. 2. 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 terjebak pada minimum lokal karena model awal yang kurang baik. 3. Metode inversi iterative damped least square merupakan metode inversi stabil dan relatif lebih cepat (membutuhkan waktu ± 20 - 30 menit) untuk diterapkan pada pemodelan tomografi dalam mencitrakan struktur tiga dimensi bawah permukaan, sedangkan metode simulated annealing belum maksimal dan relatif lebih lama (membutuhkan waktu ± 18 - 19 jam) untuk diterapkan pada pemodelan tomografi, sebagai perbandingan pada processor AMD Athlon II X2 250 untuk data yang sama. 4. Solusi model kecepatan yang dihasilkan oleh metode iterative damped least square lebih baik daripada solusi yang dihasilkan oleh metode simulated annealing. Solusi model kecepatan dari inversi simulated annealing bersifat lebih acak (tidak berpola). Selain itu, solusi model 85
kecepatan yang dihasilkan oleh metode iterative damped least square memiliki sebaran nilai delay time dan RMS error yang lebih kecil sehingga solusi tersebut terbilang lebih baik secara statistik. 5. Terdapat zona dengan rasio Vp/Vs relatif rendah (berkisar 1.45 – 1.65) yang kemungkinan diinterpretasikan sebagai steam-saturated rock pada area reservoir lapangan geotermal “RR”. Zona reservoir ini berada pada 12 – 18 km arah Barat-Timur, 8 – 12 km arah Utara-Selatan, dan kedalaman 1 – 3 km dari MSL. Keberadaan zona reservoir didukung oleh data well-trajectory, dimana zona Vp/Vs tinggi berada di sekitar sumur injeksi dan zona Vp/Vs rendah berada di sekitar sumur produksi. 6. Ada kemungkinan bahwa sistem fasa pada reservoir lapangan geotermal “RR” mengalami perubahan dari water-saturated ke steam-saturated.
5.2
Saran
Saran untuk penelitian di masa yang akan datang adalah: 1. Algoritma inversi simulated annealing perlu dikembangkan lagi untuk menyelesaikan kasus non-linier dengan banyak parameter seperti tomografi. 2. Ada baiknya proses inversi yang dilakukan menggunakan pembobotan untuk setiap data first break gelombang sehingga solusi inversi yang diperoleh semakin baik. 3. Pada penelitian ini hanya dilakukan inversi tomografi terhadap parameter kecepatan gelombang seismik P dan S. Untuk penelitian selanjutnya, dapat dilakukan penelitian untuk mencari parameter seismik lain seperti atenuasi ataupun konstanta anisotropi sehingga dapat diperoleh hasil interpretasi yang lebih baik dan detail mengenai lapangan geotermal yang diteliti. 4. Agar dapat mengidentifikasi perubahan sistem fasa dari water-saturated ke steam-saturated pada reservoir lapangan geotermal “RR”, baiknya dapat dilakukan pemodelan 4-D tomografi secara time-lapse untuk melihat karakteristik struktur kecepatan bawah permukaan dari awal produksi sampai sekarang.
86
DAFTAR PUSTAKA Afnimar, 2009. Seismologi, Bandung: ITB. Chevron Salak Geothermal. Dickson, M.H dan Fanelli, M., 2004. What is Geothermal Energy?, International Geothermal Association. UNESCO publication. Goff, F. dan Janik, C.J., 2000. Geothermal systems, In Encyclopedia of Volcanoes, H. Sigurdsson, B.F. Houghton, S.R. McNutt, H. Rymer dan J. Stix (eds.), Academic Press. Grandis, H., 2009. Pengantar pemodelan inversi geofisika, Penerbit Himpunan Ahli Geofisika Indonesia, Jakarta. Gunasekara, R.C., Foulger, G.R., dan Julian, B.R., 2003. Reservoir depletion at The Geysers geothermal area, California, shown by four-dimensional seismic tomography, J. Geophys. Res., Vol. 109, No. B3, 2134. Havskov, J. dan Ottemӧller, L., 2010. Routine data processing in earthquake seismology, Springer. Hochstein, M.P. dan Browne, P.R.L., 2000. Surface manifestation of geothermal systems with volcanic sources, In Encyclopedia of Volcanoes, H. Sigurdsson, B.F. Houghton, S.R. McNutt, H. Rymer dan J. Stix (eds.), Academic Press. Kissling, E., Ellsworth, W.L., Eberhart-Phillips, D., dan Kradolfer, U., 1994. Initial reference models in local earthquake tomography, J. Geophys. Res., Vol. 99. Menke, W., 1984. Geophysical data analysis: Discrete inverse theory, Academic Press. Nugraha, A. D., Suantika, G., dan Widiyantoro, S., 2006. Relocation of the volcano earthquakes using 3-D seismic velocity models for Guntur volcano in Garut, Indonesia, Indonesian J. of Geophys., No. 2, hal 20-26.
87
Paige, C.C. dan Saunders, M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Soft., 8(43-71), 195-209. Stewart, R.R., 1991. Exploration Seismic Tomography: Fundamentals, Society of Exploration Geophysicist, Alberta. Stimac, J., Nordquist, G., Suminar, A., dan Azwar, L.S., 2008. An Overview of The Awibengkok Geothermal System, Indonesia, Geothermics 37, Elsevier. Syahputra, A., 2011. Pengembangan perangkat lunak tomografi seismik 2-D dan 3D: aplikasi pada tomografi lubang bor dan gunungapi, Tugas Akhir Sarjana, Program Studi Teknik Geofisika, ITB, Bandung. Takei, Y., 2002. Effect of pore geometry on VP/VS: From equilibrium geometry to crack, J. Geophys. Res. Vol. 107, No. B2, 2043. Tatham, R. H. dan McCormack, M. D., 1991. Multicomponent seismology in petroleum exploration, Edited by E. B. Neitzel and D. F. Winterstein, Soc. Expl. Geophys., Tulsa, hal 248. Trampert, J., Vacher1, P., Vlaar, N., 2001. Sensitivities of seismic velocities to temperature, pressure and composition in the lower mantle, Phys. Earth Planet. Int. 124, hal 255–267 Um, J. dan Thurber, C., 1987. A fast algorithm for two-point seismic ray tracing, Bull. siesm. Soc. Am., 77, hal 972-986. Vakil-Baghmishes, M. dan Navarbaf, A., 2008. Modified very fast simulated annealing algorithm. 2008 International Symposium on Telecommunications. Van der Hilst, R.D dan Sparkman, W., 1989. Importance of the reference model in the linearized tomography and images of subduction below the Caribbean plate. Geophys. Res. Lett., 16, 1093-1096. Van der Hilst, R.D. dan Engdahl, E.R., 1991. On ISC PP dan pP data and their use in delay time tomography of the Caribbean region. Geophys. J. Int., 106, 169188. 88
Wang, Z., M. L. Batze, dan A. M. Nur., 1990. Effect of different pore fluids on seismic velocities in rock, Can. J. Explor. Geophys., Vol. 26 NOS. 1 & 2, hal 104-112. Weber, Z., 2000. Seismic travel time tomography: a simulated annealing approach, Elsevier. Physics of Earth and Planetary Interiors, 199, hal 149-159. White, S.R., 1984. Concepts of scale in simulated annealing, Proc. IEEE Int. Conference on Computer Design, Port Chester, hal 646–651. Widiyantoro, S., Gorbatoc, A., Kennett, B.L.N., dan Fukao, Y., 2000. Improving Global Shear Wave Traveltime Tomography Using Three Dimensional Ray Tracing and Iterative Inversion. Geophys. J. Int., 141, hal 747-758.
89
Penentuan Hiposenter Gempa Mikro Menggunakan Metode Inversi Simulated Annealing untuk Studi Kasus Lapangan Panas Bumi “RR” Micro-earthquake Hypocenter Determination Using Simulated Annealing Inversion Method for “RR” Geothermal Field Case Study Rexha Verdhora RY1, Andri Dian NUGRAHA2 1
Program Studi Teknik Geofisika, Fakultas Teknik Pertambangan dan Perminyakan, ITB Kelompok Keahlian Geofisika Global, Fakultas Teknik Pertambangan dan Perminyakan, ITB
2
Sari Penentuan lokasi hiposenter melibatkan proses pencarian solusi hiposenter yang memiki error minimum antara waktu tempuh observasi dan kalkulasi. 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 panas bumi “RR” telah ditentukan menggunakan metode Geiger. Bagaimanapun, pada makalah ini, metode simulated annealing diterapkan pada data dan model kecepatan 1-D yang sama untuk merelokasi hiposenter dan meminimalisasi fungsi error. Waktu tempuh dihitung menggunakan ray tracing metode shooting. Hasil penelitian ini menunjukkan bahwa lokasi hiposenter memiliki error RMS yang lebih kecil dibandingkan dengan penelitian sebelumnya, yang secara statistik berasosiasi dengan solusi yang lebih baik. Kata-kata kunci: penentuan lokasi hiposenter, inversi, metode Geiger, simulated annealing, metode shooting. Abstract Hypocenter determination involves finding a hypocenter location that has minimum error 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 take 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 of the initial model. Previously, hypocenter location at “RR” Geothermal Field has been determined by 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 shooting method. Our results show hypocenter location has smaller root means square (RMS) error compared to the previous study that can be statistically associated with better solution. Keywords:hypocenter determination, inversion, Geiger’s method, simulated annealing, shooting method. *Jl. Ganesa 10 Bandung 40132, Tel: +62-22-2534137, Faks: +62-22-2534137, E-mail:
[email protected]
I. 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, lokasi hiposenter gempa 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. Menurut penelitian tersebut, gempa mikro pada lapangan geotermal “RR” terjadi pada rentang antara elevasi 0.75 sampai -8.43 km (MSL = 0 km). Apabila dirincikan, persebaran event dapat dibagi menjadi 4 zona yaitu zona dangkal
1
pada rentang antara elevasi 1 sampai 0 km, zona menengah pada rentang antara elevasi 0 sampai -6.5 km, zona dalam pada rentang antara elevasi -6.5 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
U
MSL
Gambar 1. Peta persebaran event yang dihasilkan oleh penelitian sebelumnya (bulat hijau) dan stasiun (segitiga biru terbalik).
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 1-D (Gambar 2) yang sama untuk menentukan lokasi hiposenter dan meminimumkan fungsi misfit. Kedua solusi akan dibandingkan dan lokasi hiposenter yang memiliki RMS error lebih kecil akan berasosiasi dengan solusi yang lebih baik secara statistik. 2 Vp Vs
1 0 -1
Kedalaman (km)
-2 -3 -4 -5 -6 -7 -8 -9
0
2 4 6 8 Kecepatan Gelombang Seismik (km/s)
Gambar 2. Plot model kecepatan untuk gelombang P (Vp) dan gelombang S (Vs)
2
II. METODOLOGI 2.1 Ray Tracing Metode Shooting Untuk menentukan waktu tempuh gelombang P dan S yang merambat dari sumber gempa ke stasiun penerima 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 Aanealing 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 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 (
)
(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 ke model 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. III. 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 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.
3
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 Vakil-Baghmishes dan Navarbaf (2008), model akan diperbarui sesuai dengan persamaan: (
)
]
(3)
sehingga pada simulasi ke-k dihasilkan dk yang bernilai [-1,1] dan besar pembaruan model di simulasi ke-k adalah: (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 Ak menggambarkan 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 U akan 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 Umax x Lmax kali. 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.
Gambar 3. Diagram alir metode inversi simulated annealing pada penentuan lokasi hiposenter.
4
IV. HASIL DAN ANALISIS 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 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. Tabel 1. Katalog data sintetis waktu tempuh gelombang Posisi Hiposenter Posisi Stasiun X (km)
Y (km)
11
11
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
4
11.01087
11.03581
-2.98444
0.003117
5
11.05361
11.01492
-3.0276
0.002157
6
11.05611
11.01796
-3.02093
0.001617
7
11.01033
11.00713
-3.00043
0.002506
8
11.02544
10.99621
-2.99265
0.002915
9
11.01741
11.00869
-3.02176
0.002446
10
11.03251
11.00522
-3.02504
0.001925
11
11.0046
10.98935
-2.98496
0.00305
12
11.00836
10.95518
-3.0037
0.002366
13
11.00018
10.98141
-2.98632
0.00247
14
11.04784
10.97016
-3.06595
0.003328
15
11.00541
11.00144
-2.97761
0.002881
16
11.0129
10.9679
-3.00386
0.002379
17
11.0188
10.97034
-3.01281
0.002822
18
10.97239
10.98823
-2.9635
0.002069
19
11.0188
11.02949
-2.99187
0.002202
20
11.04076
10.98443
-3.04058
0.002941
Dapat dilihat algoritma inversi simulated annealing 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.
5
4.2 Data Riil Data 591 event gempa mikro dan model kecepatan 1-D dari penelitian sebelumnya akan digunakan pada penelitian ini. Berikut ialah 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 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. Gempa-gempa 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.
U
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. Persebaran gempa mikro pada lapangan panas bumi “RR” akan berhubungan dengan struktur 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”.
6
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).
V. 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 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. DAFTAR PUSTAKA 1. Grandis, H., 2009. Pengantar Permodelan Inversi Geofisika. Penerbit Himpunan Ahli Geofisika Indonesia, pp. 93-100. 2. Vakil-Baghmishes, M. dan Navarbaf, A., 2008. Modified Very Fast Simulated Annealing Algorithm. 2008 International Symposium on Telecommunications. 3. Weber, Z., 2000. Seismic travel time tomography: a simulated annealing approach. Elsevier. Physics of Earth and Planetary Interiors, 199, pp 149-159. 4. White, S.R., 1984. Concepts of scale in simulated annealing. Proc. IEEE Int. Conference on Computer Design, Port Chester, pp. 646–651.
7
Applications of Least Square and Simulated Annealing Inversion Method on 3-D Seismik Velocity Tomography for “RR” Geothermal Field Case Study R. V. Ry1 and A. D. Nugraha2 1
Study Program of Geophysical Engineering, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, Jalan Ganesa No.10, Bandung 40132, Indonesia Email:
[email protected] 2 Global Geophysical Group, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, Jalan Ganesa No.10, Bandung 40132, Indonesia Email:
[email protected] ABSTRACT Observations of micro-earthquake activity in the area of geothermal can be used to detect the permeability of the structure. Micro-earthquake analyses were performed with seismic tomography for imaging the structure of 3-D seismic velocity (Vp, Vs, and the Vp/Vs ratio). The travel times through 3-D velocity model were calculated using the ray tracing pseudobending method. Inversion methods which used in tomographic modeling were iterative damped least square and simulated annealing method. The results show that the iterative damped least square method produces better solutions of velocity models, visualization and statistics are reviewed. Tomographic modeling results indicate the presence of low Vp/Vs at around 1.45 – 1.65 in the reservoir area at elevation of -1 to -3 km (MSL = 0 km). This is interpreted as steam-saturated rock in the reservoir area of “RR” geothermal field. The existences of the reservoir area are supported by the data of well-trajectory, where the zones of high Vp/Vs are around the injection wells and the zones of low Vp/Vs are around the production wells. When compared with other studies in this field, it is possible that the reservoir’s phase system has changed from water-saturated to steam-saturated.
Keywords: Seismic tomography, least square, simulated annealing, pseudo-bending. Mathematics Subject Classification: Computing Classification System:
1. INTRODUCTION Observations of micro-earthquake activity in the geothermal area can be used to detect the permeability of the structure. Changes in the permeability of the structure can be caused by pore pressure and temperature changes due to interactions between the circulating fluid reservoir with hot rock, or by the stimulation of the geothermal reservoir associated with fluid injection. Micro-earthquake analyses were performed with seismic tomography for imaging the structure of 3-D seismic velocity (Vp, Vs, and the Vp/Vs ratio). The tomographic modelings apply iterative damped least square and simulated annealing inversion method. Iterative damped least square method itself has been widely used to solve tomographic modeling and has been successfully imaged seismic velocity structure beneath the surface, while the simulated annealing method is seldom used to solve non-linear inversion problem for many parameters like tomography. In this study, simulated annealing method will be implemented in tomographic modeling. Then, the solutions will be compared with the solutions obtained from the common used inversion methods such as iterative damped least square method. Tomographic inversion results are the three-dimensional velocity model of the subsurface for the P wave (Vp) and S
1
wave (Vs). The presence of a fluid phase or partial melting in the study area can be indicated by the observation of velocity anomalous and Vp/Vs ratio. The data used are the micro-earthquakes data which recorded by 19 monitoring stations and occurred during one year in the geothermal field "RR" from the period January 2007 to December 2007. In detail, the data are composed of 573 events with 3099 phase micro seismic P wave and S wave phase 2499. The data are included the hypocenter location, the origin time of the microearthquake, the travel time of P and S waves, and the location of monitoring stations. Then the initial velocity model is required for the inversion processing. Initial velocity model used is the 1-D velocity model of the P wave phase (Vp) and S wave phase (Vs). In this study, a three-dimensional block model parameterization was used. Determination of the number and size of each block of the model depends on the area and depth of the study area, as well 3 as data distribution. Study area’s dimension was 15 x 10 x 8 km in the range area of 5 - 20 km WestEast, 5 - 15 km North-South, and elevation of 1.5 to -6.5 km (MSL = 0 km). Parameterization model 3 block was used which had dimension 1.5 x 1 x 0.5 km , so we had 10 x 10 x 16 blocks (1600 blocks). 2
(a)
(b)
Vp Vs
1 0 -1
Kedalaman (km)
-2 -3 Depth (km)
Depth (km)
-4 -5 -6 -7 -8 -9
North - South (km)
West - East (km)
0
2 4 6 8 Kecepatan Gelombang Seismik (km/s)
Seismic Velocity (km/s)
Figure 1. (a) The map of the distributions of the micro-earthquake events (red dots) and stations (blue reverse triangles). (b) The 1D velocity model of P waves phase and S waves phase.
2. METHODOLOGY 2.1 Ray Tracing Pseudo-bending The pseudo-bending method was used in this seismic tomography study as ray tracing method to determine possible ray path in 3D velocity model and calculate synthetic travel time from hypocenter to receiver. Pseudo-bending (Um & Thurber, 1987) is an approach in minimization of travel time based on Fermat’s Principle by giving small perturbations gradually on ray paths. We modified MATLAB code that has been developed by Syahputra (2011).
2.2 Iterative Damped Least Square For resolve tomography inversion, iterative damped least square is implemented. We also added norm and gradient damping to constrain blocks without ray and to produce smooth solution model, respectively (Widiyantoro et.al, 2000). The equation of matrix tomography can be defined: (
)X=( )
2
where the rectangular matrix of A contains ray-segment lengths in each of the block defined by the model parameterization, αI is the norm damping matrix, γG is the gradient damping matrix, the vector model X contains slowness perturbation from the initial model used of ray tracing and the data vector dt of delay time. The vector model X was then solved by using LSQR method (Paige & Saunders, 1972). This process was done iteratively. In this study, the values of norm damping (α) and gradient damping (γ) used were 10 and 1. The iteration will stop at maximum number of iterations which used is 12, or in addition, if the difference between the RMS value at i-th iteration and RMS on the previous iteration is smaller than 0.001.
2.3 Simulated Annealing The algorithms were developed according to White (1984) and Weber (2000), and then modified to solve this case. The initial value of T0 should be as small as possible but in such a way that virtually all simulation steps are accepted. To estimate the optimal value of T0 we used a procedure based on trial error. Then, we got the values of the cost function for a predefined number of system configurations randomly chosen around the initial configuration. T0 was chosen to be greater than or equal to the standard deviation of the previously calculated cost values. The decrement should be chosen such that a few number of simulation steps suffice to reestablish at each temperature Tk. We used simple decrement rule such as Tk+1 = 0.99Tk. Updating parameter formulated as follow (Vakil-Baghmishes dan Navarbaf, 2008): (
)
]
So the produced dk belongs to [-1,1] and the rate of change in k-th iteration will be:
The interval [xmax, xmin] represents permissive range of change in the k-th iteration. Xmax and Xmin values will be depended on possible range from first iteration. In our implementation Lk is determined such that for each values of Tk a minimum number of simulation steps, Amin is accepted, where Amin is a predefined fix number. However, as Tk approaches 0, acceptance probability decreases and eventually the value of Lk reach infinity. Consequently, Lk is not allowed to exceed some constant Lmax to avoid extremely large number of simulation steps at low values of Tk . Let U denote the number of the consecutive temperatures for which the number of the accepted simulation steps is zero. If U exceeds a predefined number, Umax, the annealing procedure stops. In other words, the procedure stops if the system configuration is not altered for Umax number of consecutive temperatures, or if the configuration remains unchanged for Umax x Lmax number of consecutive simulation steps. According to the above discussion, the implementation of the simulated annealing method used annealing parameters which took the following values: Amin = 20, Lmax = 40 , and Umax= 100.
3. RESULTS AND DISCUSSION 3.1 Checkerboard Resolution Test Checkerboard Resolution Test (CRT) is a forward modeling method that aims to test the reliability of the inversion technique used in tomographic inversion and to see the resolution throughout the model space. In this study, the resolution test will be performed on the same data configuration of the
3
hypocenter and station. Model checkerboard anomalies obtained by multiplying the +10% and -10% of the initial velocity model (model 1-D velocity). This model has dimensions of 10 x 10 x 16 blocks. Then the inversion process is carried out in accordance with the inversion algorithms that have been made. Inversion methods used are iterative damped least square and simulated annealing. Resolution test results shown in Figure 2(a) for iterative damped least squares and Figure 2(b) for simulated annealing. Then, area with good resolution was interpreted on the region which has positive and negative anomalies clearly, marked with dashed black lines. (a)
(b)
Figure 2. Tomogram of horizontal cross-section of Checkerboard Resolution Test from (a) iterative damped least square dan (b) simulated annealing at z = 0, 1, -1. -2, -3, -4, -5, dan -6 km, and area with good resolution (dashed black lines). Positive anomalies are shown by blue colour and negative anomalies are shown by red colour. 3.2 Tomographic Modeling Tomographic modeling using iterative damped least square method produced the subsurface velocity structures of P wave (Vp) and S wave (Vs) at “RR” geothermal field by the tomogram shown in Figure 3(a) and 3(c). In the process of inversion for P wave, the iteration stops at the 4th iteration with the RMS error of 0.0405 seconds. While in the process of inversion for S waves, the iteration stops at the 6th iteration with the RMS error of 0.0987 seconds.
4
(a)
(b)
(c)
(d)
Figure 3. Tomogram of horizontal cross-section for (a) Vp structures from iterative damped least square and (b) simulated annealing, (c) Vs structures from iterative damped least square and (d) simulated annealing, at z = 0, 1, -1. -2, -3, -4, -5, dan -6 km, and area with good resolution (dashed black lines). Positive anomalies are shown by blue colour and negative anomalies are shown by red colour.
5
Tomographic modeling using simulated annealing method produced the subsurface velocity structures of P wave (Vp) and S wave (Vs) at “RR” geothermal field by the tomogram shown in Figure 3(b) and 3(d). In the process of inversion for P wave, the iteration stops with the RMS error of 0.0425 seconds. While in the process of inversion for the S wave, the iteration stops with the RMS error of 0.1081 seconds. Inversion process using simulated annealing method takes relatively much longer than using iterative least square method. On a computer with a processor AMD Athlon II X2 250, the comparisons are as follows. Inversion process using iterative damped least square method for P wave takes ± 20 minutes, while for the S wave takes ± 30 minutes. Inversion process using simulated annealing method for P wave takes ± 18 hours, while for the S wave takes ± 19 hours. Tomographic modeling using iterative damped least square and simulated annealing produced the tomogram with quite similar pattern for overall. The low velocity anomalies appear at elevation of 1 to -1 km, while the high velocity anomaly is more dominant at elevation of -1 to -3 km. However in detail, in the tomogram of simulated annealing, there are some anomalies that appear randomly pace or not patterned. The solutions of velocity model from simulated annealing inversion are more random (not patterned) compared to the solutions from iterative damped least square inversion. Simulated annealing inversion involved random search process to find solutions of the model parameter. There is a possibility that the solutions have not yet reached the global minimum of the misfit function. We concluded that in this tomographic modeling, the solutions of velocity model produced by iterative damped least square method are better than the solutions produced by simulated annealing method. It is also supported by the distribution of misfit value from both solutions, where the solutions of velocity model generated by iterative damped least square method have smaller value of delay time distribution and RMS error, so the solutions are better statistically. (a) RMS error = 0.0405 s
Occurences
RMS error = 0.0425 s
Occurences
dt (s)
dt (s)
(b) RMS error = 0.0987 s
Occurences
RMS error = 0.1081 s
Occurences
dt (s)
dt (s)
Figure 4. Histrogram of delay time distribution for the solutions of (a) Vp and (b) Vs from inversion results using iterative damped least square (left) dan simulated annealing (right). 3.3 Analysis of Vp, Vs, dan Vp/Vs Ratio In the geothermal system, the existence of zones with high temperatures in the subsurface varying influence on the values of Vp and Vs. On the state of water-saturated rock, the values of P wave velocity will tend to be higher than the values of the S wave velocity so that the values of Vp/Vs ratio
6
are relatively higher. While in the state of steam-saturated rock, the values of P wave velocity will tend to be lower than the values of the S wave velocity so that the values of Vp/Vs ratio are relatively lower. (Wang et al., 1990 and Takei, 2002). Previous studies conducted on the “RR” geothermal field has described reservoir area boundary along with well-trajectory of production and injection wells (Stimac et al., 2008). In order to view the velocity structures around the reservoir area, the vertical slices (cross-section C-C ', D-D', E-E ', and F-F') are performed on tomographic models around the reservoir area, shown in Figure 5. Tomographic models used are the velocity models produced by iterative damped least square inversion. (a) Vp
(b) Vs
(c) Vp/Vs
Vp and Vs
Vp/Vs
Figure 5. Tomogram of vertical cross-section (North-South) from: (a) Vp, (b) Vs, and (c) Vp/Vs ratio at U-S 9, 10, 11, dan 12 km , and well-trajectory (black lines). Positive anomalies are shown by blue colour and negative anomalies are shown by red colour.
7
At elevation of 1 to -1 km, the zones of high Vp/Vs ratio (around 1.8 - 2) can be interpreted to be associated with water-saturated rock. These zones are located at a shallow depth. It could be the zones of condensation of reservoir’s steam in the ground water or the zones of cap rock which generally, altered clays contain a lot of water. At elevation of -1 to -3 km, the zones of low Vp/Vs ratio (arround 1.45 - 1.65) can be interpreted to be associated with steam-saturated rock. These zones can be identified as the reservoir of “RR” geothermal field. The reservoir zones are located at 12 - 18 km West-East, 8 - 12 km North-South, and depth of 1 - 3 km from the MSL. Previous studies stated that the reservoir in the "RR" geothermal field is the dominance of the water reservoir (Stimac et al., 2008). It becomes interesting when tomographic imaging identify the zones in reservoir area to be associated with steam-saturated rock. Gunasekara et al. (2003) mention that the exploitation and production of the geothermal field can cause changes in the reservoir phase system. The progressive depletion of pore fluid causes the replacement of pore liquid with vapor. In addition, the pressure drop in the reservoir causes a decrease in the boiling point, resulting in boiling and vapor phase is formed. However, whether changes occur in the system phase of “RR” geothermal field reservoir, it could not be concluded yet. Volume forms of the zones of Vp/Vs ratio are shown in Figure 6. In general, it can be seen that the production wells are around the zone of low Vp/Vs and the injection wells are arround the zone of high Vp/Vs ratio. Low Vp/Vs
High Vp/Vs
Depth (km)
Depth (km)
N
N
West - East (km)
West - East (km)
Figure 6. (a)Volume forms of low Vp/Vs zones and well-trajectory of production wells (black lines). (b)Volume forms of high Vp/Vs zones and well-trajectory of injection wells (magenta lines).
4. CONCLUSION The iterative damped least square method is a stable and relatively faster inversion method (± takes 20-30 minutes) to be applied in tomographic modeling for imaging three-dimensional velocity structure, while the simulated annealing method is relatively slower (± takes 18-19 hours) to be applied in tomographic modeling, as a comparison in the processor AMD Athlon II X2 250 for the same data. The solutions of velocity model produced by iterative damped least square method are better than the solutions produced by simulated annealing method. The solution of velocity model from simulated annealing inversion is more random (not patterned). In addition, the solutions produced by iterative damped least square method have smaller values of delay time distribution and RMS error, so that the solutions are better statistically. Tomographic modeling results indicate the presence of low Vp/Vs at around 1.45 – 1.65 in the reservoir area at elevation of -1 to -3 km (MSL = 0 km). This is interpreted as steam-saturated rock in the reservoir area of “RR” geothermal field. The existences of the reservoir area are supported by the
8
data of well-trajectory, where the zones of high Vp/Vs are around the injection wells and the zones of low Vp/Vs are around the production wells. When compared with other studies in this field, it is possible that the reservoir’s phase system has changed from water-saturated to steam-saturated.
5. REFERENCES Gunasekara, R.C., Foulger, G.R., dan Julian, B.R., 2003. Reservoir depletion at The Geysers geothermal area, California, shown by four-dimensional seismic tomography, J. Geophys. Res., Vol. 109, No. B3, 2134. Nugraha, A. D., Suantika, G., dan Widiyantoro, S., 2006. Relocation of the volcano earthquakes using 3-D seismic velocity models for Guntur volcano in Garut, Indonesia, Indonesian J. of Geophys., No. 2, hal 20-26. Stimac, J., Nordquist, G., Suminar, A., dan Azwar, L.S., 2008. An Overview of The Awibengkok Geothermal System, Indonesia, Geothermics 37, Elsevier. Syahputra, A., 2011. Pengembangan perangkat lunak tomografi seismik 2-D dan 3-D: aplikasi pada tomografi lubang bor dan gunungapi, Tugas Akhir Sarjana, Program Studi Teknik Geofisika, ITB, Bandung. Takei, Y., 2002. Effect of pore geometry on VP/VS: From equilibrium geometry to crack, J. Geophys. Res. Vol. 107, No. B2, 2043. Um, J. dan Thurber, C., 1987. A fast algorithm for two-point seismic ray tracing, Bull. siesm. Soc. Am., 77, hal 972-986. Vakil-Baghmishes, M. dan Navarbaf, A., 2008. Modified very fast simulated annealing algorithm. 2008 International Symposium on Telecommunications. Van der Hilst, R.D. dan Engdahl, E.R., 1991. On ISC PP dan pP data and their use in delay time tomography of the Caribbean region. Geophys. J. Int., 106, 169-188. Wang, Z., M. L. Batze, dan A. M. Nur., 1990. Effect of different pore fluids on seismic velocities in rock, Can. J. Explor. Geophys., Vol. 26 NOS. 1 & 2, hal 104-112. Weber, Z., 2000. Seismic travel time tomography: a simulated annealing approach, Elsevier. Physics of Earth and Planetary Interiors, 199, hal 149-159. White, S.R., 1984. Concepts of scale in simulated annealing, Proc. IEEE Int. Conference on Computer Design, Port Chester, hal 646–651.
9
