Kontribusi Fisika Indonesia Vol. 12 No. 2, April 2001
Inversi Data Magnetotellurik 1-D Menggunakan Metoda Simulated Annealing Akhmad Syaripudin dan Hendra Grandis Program Studi Geofisika, Departemen Geofisika dan Meteorologi - ITB Jl. Ganesa 10 Bandung – 40132 e-mail :
[email protected] Abstrak Makalah ini membahas metoda inversi non-linier menggunakan teknik simulated annealing yang meniru proses termodinamika pendinginan substansi hingga mencapai keadaan setimbang dengan energi minimum. Metoda tersebut diterapkan pada inversi data magnetotellurik 1-D dimana parameter model adalah resistivitas yang bervariasi terhadap kedalaman. Perturbasi model dilakukan hingga dicapai misfit minimum antara data perhitungan dengan data pengamatan. Inversi data sintetik dan data lapangan menunjukkan hasil yang cukup memuaskan. Dengan memperhitungkan faktor ekivalensi, model sintetik dapat diperoleh kembali dengan cukup baik dengan data misfit sekitar 10%. Hasil inversi data lapangan menunjukkan distribusi resistivitas bawah-permukaan yang konsisten dengan kondisi geologi setempat. Kata kunci : Magnetotellurik, Inversi Abstract This paper describes a non-linear inversion method employing simulated annealing technique that imitates a thermodynamic process in which a substance is cooled down slowly to achieve an equilibrium state with a minimum energy. The method was applied to invert 1-D magnetotelluric data in which the model parameter is resistivity varied with depth. Model perturbations carried out in search for a minimum data misfit. Inversions of synthetic and field data showed satisfactory results. By considering equivalence problems, synthetic models were resolved relatively well with data misfit around 10%. Inversions of field data showed subsurface resistivity distribution that agrees well with the geological condition of the surveyed area. Keywords : Magnetotellurics, Inversion fungsi obyektif pada ruang model tersebut dilakukan secara acak (random) menggunakan metoda simulated annealing. Metoda simulated annealing adalah teknik optimasi global yang meniru proses termodinamika pembentukan kristal hasil pendinginan substansi secara perlahan hingga dicapai kesetimbangan dengan energi minimum. Dalam hal inversi energi dianalogikan sebagai misfit data. Pengujian metoda simulated annealing dilakukan melalui inversi data sintetik yang merupakan respons model sintetik sederhana (3 lapisan). Inversi data lapangan menghasilkan distribusi resistivitas bawah permukaan yang dapat menjelaskan kondisi geologi daerah penyelidikan.
1. Pendahuluan Pemodelan geofisika dilakukan untuk memperkirakan distribusi sifat fisis bawah-permukaan berdasarkan data yang diukur di permukaan bumi. Pada pemodelan inversi dicari model optimum yang berasosiasi dengan minimum suatu fungsi obyektif. Pada umumnya fungsi obyektif merupakan selisih kuadratik respons model dengan data observasi. Pada kasus magnetotellurik (MT) 1-D resistivitas hanya bervariasi terhadap kedalaman sehingga model dapat direpresentasikan oleh lapisan-lapisan horisontal dengan parameter model adalah resistivitas dan ketebalan tiap lapisan. Meskipun pemodelan MT 1-D relatif sederhana, hubungan antara data dengan parameter model sangat tidak linier sehingga inversi dengan pendekatan linier sering dianggap kurang memadai1). Untuk menghindari permasalahan akibat linierisasi maka dilakukan inversi secara non-linier menggunakan teknik optimasi global dalam pencarian minimum fungsi obyektif. Ruang model dimana terdapat solusi (model) yang dicari didefinisikan oleh batas minimum dan maksimum harga setiap parameter model. Pencarian minimum
2. Metodologi 2.1 Forward Modeling Perhitungan forward modeling MT 1-D dilakukan dengan menggunakan hubungan rekursif impedansi 2) atau medan EM di permukaan dua lapisan yang berurutan. Dalam bentuk persamaan matriks medan EM yang saling ortogonal (Ex, Hy)
49
50
KFI Vol. 12 No. 2 2001
di permukaan bumi (E1, H1) yang terdiri dari N lapisan dapat ditulis sebagai berikut 3) :
~ ⎛ E1 ⎞ N −1 ⎛ a j ⎜~ ⎟= ∏ ⎜ ⎜ H ⎟ j =1 ⎜ c j ⎝ ⎝ 1⎠
b j ⎞ ⎛ EN ⎟⎜ a j ⎟⎠ ⎜⎝ H N
⎛ EN = ∏ T j ⎜⎜ j =1 ⎝HN N −1
⎞ ⎟⎟ ⎠
(1)
⎞ ⎟⎟ ⎠
⎛ − ∆E ⎞ ⎟ ⎝ T ⎠
Pa = exp ⎜
dengan a j =1 + exp(−2 k j h j )
b j = Z I , j (1 − exp(−2 k j h j )) c j = Z I−,1j
diperoleh konfigurasi parameter model dengan misfit yang semakin kecil. Dalam penelitian ini digunakan algoritma Metropolis 4, 5). Pada setiap perturbasi model ditentukan apakah hasil perturbasi diterima atau ditolak. Jika misfit menjadi lebih kecil maka model diterima. Jika misfit menjadi lebih besar maka model diterima dengan probabilitas Pa :
(1 − exp(−2 k j h j ))
dimana Z I , j adalah impedansi intrinsik lapisan kej yang merupakan fungsi resistivitas (ρj) dan
~
ketebalan lapisan (hj). Perbedaan E1 dengan E1
~
serta H 1 dengan H 1 tidak perlu diperhitungkan (akan saling menghapuskan) karena hasil akhir yang diperlukan adalah impedansi, yaitu rasio antara medan listrik dan medan magnetik 3). Hasil perkalian N − 1 matriks T j (2×2) adalah matriks
S (2×2) dengan komponen-komponen (s11, s12, s21, s22) sehingga impedansi di permukaan lapisan kesatu (permukaan bumi) adalah : Z1 =
~ E1 E1 = ~ H1 H1
Z1 =
s11 Z I , N + s12 s11 E N + s12 H N = s21 E N + s22 H N s21 Z I , N + s22
(2)
Respons model 1-D dalam bentuk kurva sounding resistivitas dan fasa sebagai fungsi periode dapat dihitung dari impedansi yang merupakan besaran kompleks. 2.2 Metoda Simulated Annealing Metoda simulated annealling adalah prosedur optimasi global yang meniru proses termodinamika pembentukan kristal. Subtansi pada temperatur tinggi dengan konfigurasi acak didinginkan secara perlahan sampai terbentuk konfigurasi kristal dengan energi minimum. Konfigurasi parameter model menggambarkan kondisi sistem dengan tingkat misfit (atau energi) tertentu. Pada temperatur tinggi, perturbasi model yang menghasilkan peningkatan misfit masih dimungkinkan meskipun probabilitasnya relatif kecil dibandingkan dengan perturbasi yang menghasilkan penurunan misfit. Sejalan dengan proses penurunan faktor “temperatur” maka perturbasi yang menghasilkan penurunan misfit akan lebih dominan hingga
(3)
dimana ∆E = E n − E n −1 adalah selisih misfit sesudah dan sebelum perturbasi pada iterasi ke-n. Jika ∆E ≥ 0 maka dibangkitkan bilangan random ξ dengan distribusi uniform pada interval [0, 1]. Jika ξ ≤ Pa maka model diterima, jika sebaliknya maka perturbasi ditolak dan diambil model pada iterasi sebelumnya. Perturbasi yang menghasilkan misfit lebih besar masih memiliki kemungkinan untuk diterima sehingga algoritma ini memiliki kemampuan untuk menghindari minimum lokal. Probabilitas Pa pada persamaan (3) ditentukan oleh misfit dan faktor temperatur T. Pada temperatur tinggi, Pa akan cukup besar meskipun ∆E ≥ 0 sehingga perturbasi hampir selalu diterima. Temperatur diturunkan secara perlahan sehingga probabilitas perturbasi dengan ∆E ≥ 0 akan semakin kecil. Penurunan temperatur tidak boleh terlalu cepat karena akan menyebabkan iterasi terjebak pada minimum lokal. Jika penurunan temperatur terlalu pelan maka diperlukan jumlah iterasi yang sangat besar untuk mencapai misfit minimum. Pada penelitian ini digunakan skema penurunan temperatur yang dinyatakan oleh :
Tn = T0 x n
(4)
dimana T0 dan Tn adalah temperatur awal dan temperatur pada iterasi ke n, dan x < 1 adalah faktor penurunan temperatur. 2.3 Parameterisasi Model Pada model bumi berlapis horisontal dengan parameter model resistivitas dan ketebalan lapisan diperlukan informasi mengenai jumlah lapisan. Untuk menghindari pengaruh pemilihan jumlah lapisan, medium didiskretisasi menjadi sejumlah lapisan yang cukup banyak (N = 20 atau lebih) dengan ketebalan homogen dalam skala logaritmik. Ketebalan lapisan yang meningkat terhadap kedalaman dimaksudkan untuk mengakomodasi berkurangnya resolusi metoda MT terhadap kedalaman. Dalam iterasi ketebalan lapisan dibuat tetap sehingga parameter model adalah resistivitas
KFI Vol. 12 No. 2 2001
51
tiap lapisan yang akan menggambarkan variasi resistivitas terhadap kedalaman. Pada awal iterasi, resistivitas tiap lapisan dibuat sama dengan rata-rata data resistivitas semu (data). Resistivitas lapisan ke i diperturbasi dengan membangkitkan bilangan random uniform dalam interval [0, 1] yang kemudian dipetakan ke dalam interval [ρmin, ρmaks] sehingga diperoleh :
ρ imin
≤ρ ≤ i
ρ imaks
(5)
Resistivitas yang dapat dipilih dalam interval tersebut adalah harga-harga diskret dengan interval uniform dalam skala logaritmik. Harga ρmin = 1 Ohm.m dan ρmaks = 1000 Ohm.m untuk tiap lapisan dipilih secara a priori dan dianggap telah mewakili medium konduktif dan medium resistif. Fungsi misfit dapat dinyatakan oleh :
⎛ ρ cal ND ⎛ a,i ⎜ S e = ∑ ⎜ log ⎜ obs ⎜ i =1 ⎜ ⎝ ρ a,i ⎝
⎞ ⎟ + λ φcal − φobs i i ⎟ ⎠
⎞ ⎟ ⎟⎟ (6) ⎠
dimana ρa adalah resistivitas semu, φ adalah fasa dalam radian dan λ merupakan faktor pembobot untuk fasa. Mengingat ketebalan lapisan yang relatif kecil dan jumlah lapisan cukup banyak maka fungsi misfit pada persamaan (6) kurang sensitif terhadap perubahan resistivitas lapisan. Hal ini disebabkan oleh masalah ekivalensi atau ambiguitas dimana terdapat banyak model dengan respons yang cocok dengan data 6). Oleh karena itu digunakan faktor smoothing atau penghalusan model yang pada dasarnya adalah memberikan probabilitas lebih besar pada resistivitas yang tidak jauh berbeda dengan resistivitas lapisan di atas dan di bawah suatu lapisan. Faktor smoothing untuk lapisan ke i adalah sebagai berikut : ⎛ ρ ⎞ ⎛ ρ ⎞ S m = log ⎜⎜ i ⎟⎟ + log ⎜⎜ i ⎟⎟ ⎝ ρ i −1 ⎠ ⎝ ρ i +1 ⎠
(7)
Dengan demikian misfit total pada iterasi ke n yang digunakan pada perhitungan probabilitas Pa adalah :
En = Se + β S m
(8)
dimana β adalah pembobot untuk faktor smoothing yang dapat dipilih secara coba-coba sampai diperoleh hasil yang memadai. 3. Hasil dan Pembahasan 3.1 Inversi Data Sintetik Pengujian metoda simulated annealing dilakukan melalui inversi data sintetik untuk mengetahui efektivitas metoda tersebut dalam
memperoleh kembali model sintetik. Dua model sintetik yang digunakan mewakili model sederhana yang teridiri dari 3 lapisan yaitu lapisan konduktif pada medium resistif (model-1) dan lapisan resistif pada medium konduktif (model-2). Parameter model sintetik tersebut ditampilkan pada tabel berikut : Tabel model sintetik-1 Lapisan
Kedalaman (m)
1 2 3
100 600 -
Resistivitas (Ohm.m) 250 10 1000
Tabel model sintetik-2 Lapisan
Kedalaman (m)
1 2 3
200 900 -
Resistivitas (Ohm.m) 10 1000 5
Data sintetik adalah resistivitas semu dan fasa sebagai fungsi perioda (10-3 sampai 10+2 detik) yang telah ditambah noise dengan distribusi normal dengan rata-rata 0 dan standar deviasi 10% dari data tanpa noise. Jumlah lapisan yang digunakan adalah 20 lapisan pada interval kedalaman 10 sampai 1000 m. Pada semua inversi digunakan T0 = 5, x = 0,99 dan jumlah iterasi 300 sehingga penurunan temperatur cukup lambat namun waktu eksekusi hingga diperoleh misfit minimum masih dapat dianggap memadai. Gambar 1 menampilkan hasil inversi data sintetik-1 yang diperoleh dengan menggunakan λ = 0,1 dan β = 1,2. Tampak bahwa model inversi dapat merepresentasikan model sintetik dengan cukup baik meskipun terdapat perbedaan yang disebabkan oleh diskretisasi lapisan dan masalah ekivalensi. Secara umum respons model inversi (hanya resistivitas semu yang ditampilkan) cocok dengan data sintetik. Gambar 2 memperlihatkan hasil inversi data sintetik-2 yang diperoleh dengan menggunakan λ = 0,1 dan β = 0,5. Pola yang hampir sama dengan hasil inversi data sintetik-1 diperoleh untuk model inversi, fluktuasi misfit maupun kecocokan data perhitungan dengan data sintetik. Pada inversi data sintetik digunakan faktor pembobot fasa (λ) cukup kecil sehingga misfit untuk fasa menjadi lebih besar. Hal tersebut dimaksudkan untuk mengetahui sejauh mana fasa berpengaruh dalam memperoleh kembali model sintetik. Pada data lapangan umumnya kualitas fasa kurang baik atau bahkan tidak dapat digunakan sama sekali.
52
KFI Vol. 12 No. 2 2001
(a)
(a) resistivitas (Ohm.m) 1 10
10
__
1000
1
100
10
__
1000
....
sintetik
100
1000
inversi sintetik
100
1000
(b)
(b) 1E+3
1E+3
1E+2
1E+1
1E+0
__
1E-1
....
inversi
resistivitas semu (Ohm.m)
resistivitas semu (Ohm.m)
10
inversi
kedalaman (m)
kedalaman (m)
....
resistivitas (Ohm.m)
100
1E+2
1E+1
1E+0
__
1E-1
....
sintetik
1E-2
inversi sintetik
1E-2 1E-3
1E-2
1E-1
1E+0
1E+1
1E+2
1E-3
perioda (detik)
1E-2
1E-1
1E+0
1E+1
1E+2
perioda (detik)
(c)
40
30
30
misfit
misfit
(c) 40
20
10
20
10
0
0 0
100
200
300
iterasi
Gambar 1. (a) Model hasil inversi data sintetik-1 (b) Perbandingan antara data perhitungan dan data sintetik (resistivitas semu). (c) Misfit sebagai fungsi iterasi.
0
100
200
300
iterasi
Gambar 2. (a) Model hasil inversi data sintetik-2 (b) Perbandingan antara data perhitungan dan data sintetik (resistivitas semu). (c) Misfit sebagai fungsi iterasi.
KFI Vol. 12 No. 2 2001
resistivitas (Ohm.m) 1
10
100
1000
kedalaman (m)
100
karbonat
1000
10000
(b) 1E+3
resistivitas semu (Ohm.m)
Data lapangan diperoleh dari suatu daerah eksplorasi hidrokarbon di Pulau Seram, Maluku. Konfigurasi atau struktur geologi daerah tersebut tidak dapat menghasilkan data seismik yang baik. Data log menunjukkan bahwa formasi batuan karbonat yang kemungkinan mengandung reservoir memiliki resistivitas yang cukup tinggi (antara 100 sampai 2000 Ohm.m dengan rata-rata sekitar 700 Ohm.m). Hal ini cukup kontras dengan resistivitas formasi lainnya yang berkisar antara 1 sampai 100 Ohm.m dengan rata-rata sekitar 10 Ohm.m. Pengukuran MT dilakukan pada 3 lokasi (A, B dan C) untuk mencoba mengidentifikasi formasi batuan karbonat tersebut. Koreksi efek statik (static shift) dilakukan dengan menggunakan data time domain EM (TDEM) 7). Pada semua inversi data lapangan digunakan faktor pembobot fasa λ = 0 karena kualitas data fasa umumnya kurang baik. Gambar 3 memperlihatkan model resistivitas dan misfit data hasil inversi untuk lokasi A dengan β = 1,2. Pada kedalaman kurang dari 1000 m terdapat lapisan sangat konduktif (1 – 2 Ohm.m) yang diselingi oleh lapisan dengan resistivitas 10 Ohm.m (antara kedalaman 200 m sampai 600 m). Pada kedalaman sekitar 1500 m terdapat peralihan yang cukup jelas dari lapisan konduktif (kurang dari 10 Ohm.m) ke resisitif (100 – 1000 Ohm.m). Sesuai dengan indikasi data log di dekat lokasi A, lapisan resistif tersebut dapat ditafsirkan sebagai formasi batuan karbonat. Pada Gambar 4 ditampilkan hasil inversi data dari lokasi B dimana β = 2. Tampak bahwa lapisan dekat permukaan (sampai kedalaman 250 m) relatif konduktif yaitu sekitar 10 Ohm.m yang diikuti dengan lapisan sangat konduktif 2 Ohm.m sampai kedalaman sekitar 1200 m. Peralihan ke lapisan resistif terjadi secara gradual dan harga resistivitas yang tidak terlalu tinggi (sekitar 100 Ohm.m). Hal ini mengindikasikan bahwa keberadaan formasi batuan karbonat di lokasi B tidak terlalu signifikan. Hal tersebut sesuai dengan hasil interpretasi geologi yang menunjukkan formasi batuan karbonat di daerah ini mengalami pembajian. Hasil inversi data dari lokasi C (dengan β = 0,8) menunjukkan bahwa lapisan resistif yang dapat diasosiasikan sebagai formasi batuan karbonat terdapat pada kedalaman yang lebih besar dibanding di lokasi A dan B yaitu lebih dari 2500 m (Gambar 5). Kualitas data yang kurang baik (khususnya di lokasi C) dan kedalaman di luar jangkauan data MT menyebabkan ketebalan lapisan formasi batuan karbonat di semua lokasi tidak dapat diperkirakan.
(a)
1E+2
1E+1
1E+0
__
1E-1
....
inversi data
1E-2 1E-2
1E-1
1E+0
1E+1
1E+2
1E+3
perioda (detik)
(c) 40
30
misfit
3.2 Inversi Data Lapangan
53
20
10
0 0
100
200
300
iterasi
Gambar 3. (a) Model hasil inversi data lapangan di lokasi A. (b) Perbandingan antara data perhitungan dan data lapangan (resistivitas semu). (c) Misfit sebagai fungsi iterasi.
54
KFI Vol. 12 No. 2 2001
(a)
(a)
resistivitas (Ohm.m)
resistivitas (Ohm.m) 1
10
1
100
kedalaman (m)
kedalaman (m)
100
1000
100
100
karbonat
1000
10000
1000
karbonat
10000
(b)
(b) 1E+3
1E+3
resistivitas semu (Ohm.m)
resistivitas semu (Ohm.m)
10
1E+2
1E+1
1E+0
__
1E-1
....
inversi data
1E-2
1E+2
1E+1
1E+0
__
1E-1
....
inversi data
1E-2 1E-2
1E-1
1E+0
1E+1
1E+2
1E+3
1E-2
perioda (detik)
1E-1
1E+0
1E+1
1E+2
1E+3
perioda (detik)
(c)
40
30
30
misfit
misfit
(c) 40
20
10
20
10
0
0 0
100
200
300
iterasi
Gambar 4. (a) Model hasil inversi data lapangan di lokasi B. (b) Perbandingan antara data perhitungan dan data lapangan (resistivitas semu). (c) Misfit sebagai fungsi iterasi.
0
100
200
300
iterasi
Gambar 5. (a) Model hasil inversi data lapangan di lokasi C. (b) Perbandingan antara data perhitungan dan data lapangan (resistivitas semu). (c) Misfit sebagai fungsi iterasi.
KFI Vol. 12 No. 2 2001
55
4. Kesimpulan
Daftar Pustaka.
Metoda simulated annealing untuk inversi non-linier dapat mengeksplorasi ruang model secara efektif sehingga menghasilkan solusi yang optimal. Beberapa faktor atau parameter inversi yang harus dipilih dan dapat mempengaruhi hasil inversi adalah : • parameterisasi model dan skema perturbasi model • bentuk fungsi obyektif atau fungsi misfit • penanganan masalah ekivalensi atau ketidakunikan solusi • skema penurunan temperatur
1. Grandis, H., Inversi non-linier menggunakan metoda Monte-Carlo: Kasus magnetotellurik 1dimensi, Kontribusi Fisika Indonesia, 10, 51 – 57, (1999). 2. Kaufmann, A.A., Keller, G.V., The Magnetotelluric Sounding Method, Elsevier. Amsterdam, 595 pp. 1981. 3. Grandis, H., An alternative algorithm for onedimensional magnetotelluric response calculation, Computer and Geosciences, 25, 119 – 125, (1999). 4. Dosso, S.E., Oldenburg, D.W., Magnetotelluric appraisal using simulated annealing, Geophysical Journal International, 106, 379 – 385, (1991). 5. Rothman, D.H., Automatic estimation of large residual static corrections, Geophysics, 51, 332 – 346, (1985). 6. Constable, S.C., Parker, R.L., Constable, C.G., Occam's inversion, A practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289 – 300, (1987). 7. Hendro, A., Grandis, H., Koreksi efek statik pada data magnetotellurik menggunakan data elektromagnetik transien, Prosiding PIT HAGI ke-21, Jakarta, 1996. 8. Sharma, S.P., Kaikkonen, P., Global optimisation of time domain electromagnetic data using very fast simulated annealing, Pure and Applied Geophysics, 155, 149 – 168, (1999).
Faktor-faktor tersebut sangat bergantung pada masalah yang ditinjau. Alternatif yang telah dipilih untuk inversi data MT 1-D dapat memberikan hasil yang cukup baik. Pada metoda yang relatif sederhana ini harga resistivitas a priori memiliki probabilitas uniform dalam interval [ρmin, ρmaks] yang tetap selama iterasi. Hal ini menyebabkan banyak perturbasi model yang ditolak pada saat temperatur cukup rendah. Untuk menghindari hal tersebut maka pada teknik yang dikenal sebagai very fast simulated annealing 8) digunakan fungsi distribusi probabilitas Cauchy di sekitar harga resistivitas lapisan pada setiap iterasi. Lebar interval harga resistivitas menurun sesuai dengan penurunan temperatur mengingat resistivitas model telah mendekati harga optimum pada temperatur yang cukup rendah. Informasi yang diperoleh dari hasil inversi data lapangan cukup signifikan dan sesuai dengan data lain yang tersedia (data log dan hasil interpretasi geologi). Hasil yang diperoleh dapat dijadikan model awal inversi yang menggunakan parameter model dengan jumlah lapisan yang terbatas (tipe blok) agar diperoleh batas lapisan yang lebih tegas.
