Inversi Tak-Linier Magnetotelurik Dua-Dimensi Menggunakan Algoritma Monte Carlo Rantai Markov Nugroho D. Hananto dan Djedi S. Widarto Pusat Penelitian Geoteknologi – Lembaga Ilmu Pengetahuan Indonesia Komplek LIPI Jl. Sangkuriang Bandung - 40135 e-mail :
[email protected]
ABSTRAK Dalam formulasi masalah inversi menggunakan inferensi Bayes, solusi inversi adalah fungsi peluang model posterior yang merupakan aktualisasi peluang model a priori. Pemecahan masalah inversi tak-linier dilakukan dengan algoritma Monte Carlo Markov Chain (MCMC) dimana peluang transisi dari rantai Markov homogen digunakan untuk mengeksplorasi ruang model secara efektif. Penyelesaian inversi taklinier ini memiliki kemampuan untuk menghindari jebakan minimum lokal yang sering terjadi dalam inversi dengan pendekatan linier. Penerapan pemodelan inversi MCMC pada data sintetik magnetotelurik 2-D menunjukkan bahwa inversi data mode TE atau TM secara mandiri menghasilkan citra yang kurang fokus karena sensitivitas masing-masing mode. Namun demikian, resolusi model inversi dari data mode TE dan TM secara terpadu menghasilkan citra yang memuaskan. ABSTRACT In formulation of inverse problem using Bayesian inference, solution to the inverse problem is expressed as posterior probability of model obtained as actualisation of given prior model probabilty. The non-linear inverse problem is resolved by using Monte Carlo Markov Chain (MCMC) algorithm where the transition probability of a homogenous Markov chain is employed to explore effectively the model space. This nonlinear approach has the ability to overcome local minima entraptment, which occasionally occurs in the linearized inversion. Application of the algorithm to invert separate TE and TM mode of MT synthetic data results in non-focused images due to different sensitivities of TE and TM data. However, simultaneous inversion of TE and TM data gives satisfactorily results.
1. PENDAHULUAN Inversi secara teoritis adalah suatu teori matematika yang mencoba menjelaskan bagaimana informasi tentang suatu sistem fisis terparameterisasi dapat diturunkan dari sekumpulan data hasil pengamatan (Mosegaard & Tarantola, 1995). Dalam geofisika, teori inversi telah secara luas dikembangkan dengan tujuan untuk mengambil kesimpulan tentang interior bumi dari sekumpulan data fisika yang diamati di permukaan. Pemecahan inversi biasanya dilakukan dengan mencari suatu model optimum dimana respon yang dihasilkan oleh model itu mendekati data hasil pengamatan. Hal ini dilakukan dengan meminimumkan suatu fungsi obyektif tertentu yang menggambarkan seberapa dekat data pengamatan dengan respon hasil penghitungan suatu model (Menke, 1993). Fungsi obyektif yang diminimumkan dalam inversi geofisika biasanya adalah suatu fungsi yang tak-linier. Sifat tak-linier ini timbul dari kompleksnya fenomena fisis yang diamati. Pendekatan linier dapat dilakukan untuk
2
memecahkan fungsi tak-linier ini dengan cara mendekati gradien fungsi obyektif dengan menerapkan ekspansi deret Taylor orde pertama. Informasi itu selanjutnya digunakan untuk melakukan modifikasi pada parameter model secara iteratif hingga diperoleh suatu model optimal. Pendekatan linier pada fungsi obyektif yang tak-linier seringkali memberikan solusi yang tidak konvergen dan terjebak pada minimum lokal. Kelemahan pendekatan linier ini dapat diatasi dengan melakukan pendekatan global (tanpa linierisasi) pada fungsi obyektif. Metoda stokastik berbasis Monte Carlo seperti UMC (Uniform Monte Carlo), SA (Simulated Annealing) dan GA (Genetic Algorithm) sangat menarik diterapkan untuk melakukan optimasi fungsi misfit dalam ruang parameter multi-dimensi (Sambridge, 1999). Dalam makalah ini dibahas pemecahan masalah inversi tak-linier pada data magnetotelurik dua-dimensi (MT 2-D) menggunakan skema Monte Carlo dalam formulasi statistika Bayes.
JURNAL GEOFISIKA 2004/2
Simulasi rantai Markov ergodik dilakukan untuk menyatakan distribusi (probabilty density function atau PDF) a posteriori marginal dari model parameter. Selanjutnya algoritma pemrograman diuji coba pada data respon model sintetik MT untuk mengetahui resolusi hasil inversi.
dengan mi = ρ k (dan harga elemen model lain tetap) dinyatakan oleh : p(mi = ρk d ) =
g (d mi = ρk )h( mi = ρk ) N
∑ g (d m = ρ ) h ( m = ρ ) i
j
i
(3)
j
j =1
2. TEORI Masalah inversi yang akan diselesaikan adalah menentukan suatu himpunan model parameter yang tidak diketahui m = [mi ] (i = 1,2,..., M ) dari
himpunan
data
hasil
d = [d i ]
pengamatan
(i = 1,2,..., ND). Setiap parameter model memiliki nilai tertentu yang diberikan oleh suatu himpunan nilai a priori yang mungkin ρ = [ρi ]( j = 1,2,...N ) . Dalam kasus ini, perumusan Bayes dapat dituliskan sebagai Π (m d ) =
g (d m ) h (m )
∑ g (d m) h(m)
(1)
m∈ X
dimana
g ( d m)
mengacu
pada
peluang
kondisional dari data yang dihasilkan oleh model (forward modelling), sedangkan h(m) adalah rapat peluang a priori model. Parameter X berisi himpunan NM model yang mungkin (semua kombinasi N nilai a priori yang mungkin untuk M parameter model) (Grandis dkk., 1998). Masalah yang hendak ditinjau secara khusus adalah peluang marginal model dengan elemen tertentu dari parameter model mi = ρ k tanpa mempedulikan elemen model yang lain yang dapat dituliskan sebagai :
∑ g (d m = ρ )h( m = ρ ) i
Π( mi = ρk d ) =
m∈ X i ,k
k
i
k
∑ g (d m)h(m)
m∈ X
(2) Rapat peluang posterior fungsi Π pada persamaan (1) dan (2) tidak dapat dihitung secara langsung. Namun demikian dapat dilakukan simulasi untuk memilih sampel secara efisien dari bagian ruang model yang berkontribusi terbanyak pada denominator persamaan (2). Dimulai dari elemen manapun dalam ruang model, ditentukan peluang relatif dari satu parameter, misalkan mi = ρk untuk k = 1,2,..., N dengan seluruh harga elemen lain tetap. Peluang relatif ini pada kenyataannya adalah kasus khusus dari persamaan (2). Peluang kondisional model
dimana (mi = ρk ) menyatakan secara implisit bahwa parameter model m selain mi adalah tetap. Denominator persamaan (3) dapat dihitung untuk sejumlah N nilai a priori untuk parameter model. Peluang relatif ini dapat digunakan sebagai bobot bagi nilai a priori parameter model yang dipilih untuk memperbaharui nilai mi . Peluang kondisional data yang diberikan oleh model g (d mi = ρk ) dalam persamaan (3) adalah fungsi likelihood yang berkaitan dengan distribusi galat dari data. Dengan mengasumsikan galat yang bersifat Gaussian dan peluang yang seragam untuk setiap nilai a priori dapat dituliskan bentuk yang lebih eksplisit dari persamaan (3) sebagai berikut : p( mi = ρ k ) =
exp( − E ( mi = ρ k )) N
∑ exp(− E (m = ρ )) i
(4)
j
j =1
dimana notasi vektor d telah dihilangkan untuk penyederhanaan dan E ( mi = ρk ) adalah fungsi misfit. Fungsi misfit ini didefinisikan sebagai perbedaan antara data hasil pengamatan dan data hasil penghitungan. Evaluasi persamaan (4) memerlukan pemecahan masalah pemodelan kedepan untuk vektor model m dimana pada suatu iterasi elemen ke-i adalah sama dengan ρ k untuk k = 1,2,..., N . Nilai parameter model mi diperbaharui dengan melakukan seleksi acak dari nilai-nilai a priori yang telah diberi bobot dengan peluang relatifnya. Pada tahapan ini diperlukan penyelesaian pemodelan kedepan sejumlah nilai a priori yang diberikan untuk setiap parameter model m. Parameterisasi model dan proses seleksi nilai a priori untuk m1 secara skematis diperlihatkan pada Gambar 1. Peluang relatif pada persamaan (4) merefleksikan kecenderungan pemilihan sampel pada ruang model. Harga parameter model yang menghasilkan misfit kecil akan lebih mungkin terpilih daripada parameter model dengan misfit yang besar. Dengan demikian pemilihan nilai ρ k tertentu untuk harga parameter mi dengan nilai misfit yang besar tetap dapat terjadi hanya saja peluangnya lebih kecil. Sifat ini memungkinkan
3
JURNAL GEOFISIKA 2004/2
MODEL SINTETIK
algoritma inversi non-linier ini untuk lepas dari jebakan minimum lokal (Rothman, 1986).
0
3.0
Kedalaman (m)
-150 2.5
-300 -450
2.0
-600 1.5
-750 -900
0
300
600
900
1200
1.0
Jarak (m)
Gambar 2. Model sintetik dengan anomali konduktif (10 Ωm) dan anomali resistif (1000 Ωm) pada medium 100 Ωm (skala warna adalah dalam log resistivitas).
(a) Parameterisasi model dan (b) Pemilihan harga “a priori” untuk satu elemen model.
3. HASIL DAN PEMBAHASAN Algoritma inversi diuji dengan melakukan inversi data MT sintetik 2-D. Skema pemodelan kedepan yang digunakan adalah algoritma elemen hingga yang dikembangkan oleh Uchida (1993). Model yang ditinjau berupa model blok sederhana dengan dua anomali. Benda anomali pertama bersifat konduktif (10 Ωm) sedangkan benda anomali kedua bersifat resistif (1000 Ωm) yang berada dalam lingkungan setengah ruang (halfspace) homogen dengan tahanan-jenis 100 Ωm. Tahanan-jenis semu dan fasa respon model dihitung pada 24 stasiun dengan jarak antar stasiun serba sama 50 m. Penghitungan respon sintetik dilakukan dalam mode TM (transverse magnetic), mode TE (transverse electric). Pada mode TM arah medan magnet sejajar arah struktur (strike). Frekuensi yang digunakan adalah sebanyak 11 frekuensi antara 2 Hz hingga 2048 Hz. Pada seluruh data sintetik diberikan derau Gaussian sebesar 0.3 %. Tipper pada mode TE tidak digunakan. Pada Gambar 2 ditampilkan model sintetik yang digunakan dalam penelitian ini.
Kurva misfit terhadap iterasi selanjutnya digambarkan untuk melihat konvergensi inversi. Pada selang dimana konvergensi telah tercapai, dilakukan penghitungan harga parameter model rata-rata (berdasarkan estimasi Bayes) sebagai model inversi. Penghitungan variansi model tidak dilakukan karena jumlah iterasi yang sedikit sehingga hasilnya dianggap tidak dapat merepresentasikan ketidakpastian model inversi. Citra model hasil inversi untuk mode TM ditampilkan pada Gambar 3 sedangkan kurva misfit sebagai fungsi iterasi diperlihatkan pada Gambar 4.
MODEL INVERSI TM 0
3.0
-150 Kedalaman (m)
Gambar 1.
Inversi dilakukan pada data sintetik yang dihasilkan dari pemodelan kedepan respon model pada mode TE, mode TM, dan mode paduan TE dan TM. Inversi ini bertujuan untuk mendapatkan kembali citra model dari data sintetik. Iterasi dibatasi hanya sampai 10 iterasi dengan pertimbangan waktu komputasi yang panjang dan kecenderungan pola misfit yang diharapkan telah mengalami kejenuhan.
2.5
-300 -450
2.0
-600 1.5
-750 -900
1.0
0
300
600 Jarak (m)
900
1200 Log Tahanan-jenis (Ohm.m)
Gambar 3. Citra model hasil inversi mode TM.
4
JURNAL GEOFISIKA 2004/2
12 9 6 3 0 1
4
7
10
iterasi
Gambar 4. Kurva misfit sebagai fungsi iterasi pada inversi mode TM.
Dari citra diatas terlihat bahwa algoritma inversi telah dapat meresolusi daerah anomali resisitif sedangkan daerah anomali konduktif tidak teresolusi dengan baik. Hal ini disebabkan sifat mode TM yang lebih sensitif pada anomali anomali resistif (Uchida, 1993). Kurva misfit menunjukkan penurunan nilai misfit yang signifikan hingga iterasi ke-5. Pada iterasi berikutnya kurva misfit sedikit berosilasi hingga akhirnya konvergen pada iterasi ke-10. Osilasi ini menunjukkan kemampuan algoritma inversi untuk mengenali adanya minimum lokal yang ditemui pada iterasi tertentu dan kemudian melakukan upaya pembaharuan parameter model sebagai usaha untuk lepas dari jebakan minimum lokal tersebut. Hasil inversi pada mode TE ditampilkan pada Gambar 5 sedangkan kurva misfit sebagai fungsi iterasi inversi mode TE ditampilkan pada Gambar 6. Dari citra hasil inversi dapat diamati bahwa daerah anomali konduktif relatif dapat terresolusi sedangkan daerah anomali resistif tidak dapat terresolusi dengan baik. Hal ini disebabkan sifat inheren dari mode TE yang relatif tidak sensitif terhadap anomali resistif. Kurva misfit yang diperoleh menunjukkan adanya penurunan yang tajam pada iterasi kedua yang kemudian berosilasi hingga konvergen pada iterasi ke-10. Pola ini menunjukkan bahwa nilai misfit pada iterasi kedua merupakan salah satu minimum lokal yang dapat dihindari algoritma ini dengan jalan melakukan pembaharuan nilai parameter model. Citra hasil inversi data sintetik pada mode paduan TE dan TM ditampilkan pada Gambar 7 sedangkan kurva misfit sebagai fungsi iterasi diperlihatkan pada Gambar 8. Dari citra mode
Kurva misfit sebagai fungsi iterasi juga menunjukkan adanya kemampuan algoritma ini menghindari minimum lokal pada iterasi ke-5 dan kemudian mencoba merekonstruksi citra kembali melalui pembaharuan nilai parameter model hingga konvergen pada iterasi ke-10.
MODEL INVERSI TE 0
3.0
-150 Kedalaman (m)
15
2.5
-300 -450
2.0
-600 1.5
-750 -900
0
300
600
900
1200
1.0
Log Tahanan-jenis (Ohm.m)
Jarak (m)
Gambar 5. Citra model hasil inversi mode TE.
MODE TE DATA SINTETIK 18 misfit(x1.0e+3)
misfit (x1000)
18
paduan TE dan TM tersebut terlihat bahwa proses inversi telah menghasilkan citra yang mendekati model sintetik. Hal ini disebabkan karena informasi yang digunakan pada mode paduan ini lebih lengkap dan saling melengkapi karena data mode TE dan mode TM digunakan secara bersama-sama untuk melakukan evaluasi model parameter a priori. Penggunaan kedua mode secara bersamaan akan menambah informasi yang dapat disarikan dari kedua mode itu. Namun demikian, waktu komputasi yang diperlukan juga menjadi dua kali lipat lebih lama daripada inversi pada satu mode saja.
15 12 9 6 3 0 1
4
7
10
iterasi
Gambar 6. Kurva misfit sebagai fungsi iterasi pada inversi mode TE.
5
JURNAL GEOFISIKA 2004/2
MODEL INVERSI TE & TM 0
3.0
Kedalaman (m)
-150 2.5
-300 -450
2.0
-600 1.5
-750 -900
0
300
600
900
1200
1.0
Dalam penelitian ini juga telah diperlihatkan bahwa algoritma inversi tak-linier ini memiliki kemampuan untuk lepas dari jebakan minimum lokal. Disamping itu algoritma inversi Monte Carlo Rantai Markov bersifat umum (generik) sehingga dapat diterapkan pada penyelesaian masalah estimasi parameter berdasarkan data geofisika lainnya. Untuk itu perlu diperhitungkan penyesuaian yang menyangkut pemodelan kedepan dan adanya ambiguitas inheren pada setiap metoda atau data geofisika yang ditinjau.
Log Tahanan-jenis (Ohm.m)
Jarak (m)
UCAPAN TERIMA KASIH Gambar 7. Citra model hasil inversi mode TE dan TM.
MODE TE&TM DATA SINTETIK
Penulis menyampaikan ucapan terima kasih kepada Dr. Hendra Grandis dan rekan-rekan dari Lab. Geofisika Terapan, Program Studi Geofisika Dept. GM-ITB yang telah memberikan masukan dan diskusi untuk penulisan makalah ini.
misfit(x1.0e+3)
18 15
DAFTAR PUSTAKA
12
Grandis, H., Mogi, T., Widarto, D.S., 1998, 3-D magnetotelluric inversion using Markov chain algorithm, Proceeding of 99th SEG Japan Conference, Tokyo, Japan.
9 6 3 0 1
4
7
10
iterasi
Menke, W., 1984, Introduction to geophysical data analysis: discrete inverse theory, Academic Press, Orlando.
Gambar 8. Kurva misfit sebagai fungsi iterasi pada inversi mode TE dan TM.
Mosegaard, K. and Tarantola, A. , 1995, Monte Carlo sampling of solutions to inverse problems, J. Geophys. Res., 100, 12.431 – 12.447.
4. KESIMPULAN
Rothman,D.H., 1986, Automatic estimation of large residual statics corrections, Geophysics, 51, 332 – 346.
Penerapan algoritma Monte Carlo Rantai Markov pada inversi data sintetik menghasilkan citra tahanan-jenis yang baik untuk inversi paduan TE dan TM. Inversi mode TE menghasilkan citra tahanan-jenis daerah anomali tahanan-jenis rendah yang baik sedangkan inversi mode TM menghasilkan citra daerah anomali tahanan-jenis tinggi yang baik. Temuan ini berkaitan dengan sensitivitas data mode TE dan TM dimana mode TE lebih sensitif terhadap anomali tahanan-jenis rendah sedangkan mode TM lebih sensitif terhadap anomali tahanan-jenis tinggi.
6
Sambridge, M., 1999, Geophysical inversion with a neighbourhood algorithm – I. Searching a parameter space, Geophys. J. Int, 138, 479494. Uchida, T., 1993, Smooth 2-D inversion of magnetotelluric data based on statistical criterion ABIC, J. Geomag. Geoelectr., 45, 841 – 858.
