Programming,

Inversi 1D Occam Magnetotellurik

8:06:00 PM Leo Cahya D 12 Comments

Hello Digital World,

Beberapa waktu lalu saya sempat posting tentang nonlinearitas dari fungsi forward dalam magnetotelluric yang hasil inversinya menunjukkan multi solusi hanya dari sebuah data pengukuran. Well, kali ini saya akan mencoba membuat jenis inversi lain untuk fungsi forward magnetotelluric ini yang namanya Occam[1].



Occam ini salah satu metode inversi Linearized aka Least Square aka kelompok metode inversi ribet yang harus ngitung matriks jacobi mana sering singular nyebelin bikin inversi gagal.

Occam Inversion
Dasarnya sih mirip dengan metode Levenberg-Marquardt aka Damping Gauss-Newton. Hanya saja kita tambahkan parameter delta untuk smoothing (regulasi tikonov orde 1)[3] :


dan update parameternya :


dimana m merupakan parameter yang akan diestimasikan, J matriks jacobi, alpha adalah damping parameter dan d adalah selisih data estimasi dan data pengukuran.

Bagaimana cara membuatnya?
- Pertama kita perlu membuat forward operator-nya yang dapat dilihat di appendix[1]. Ternyata simpel ya?

- Lalu, membuat matriks Jacobi aka Sensitivitas-nya. Apa itu matriks Jacobi? ya matriks dari fungsi forward yang seperti ini :


Ada opsi 3 metode dalam membuat matriks Jacobi [2] :
1. Adjoint equation,
2. Sensitivity equation (Appendix setelah forward MT [1]),
3. Pertubation aka Numerical method.




Saya sih lebih suka menggunakan yang metode numeris (karena malas membuat program analitiknya yang lebih ribet dari forwardnya haha~).
Metode numeris yang saya gunakan adalah forward difference :
 
Rumus dan Skemanya
- Selesaikan rumus update parameter Occam secara iterative. Paling bagus setiap iterasi dibuat multiple damping parameter dan dipilih RMS error yang lebih kecil.


Kira-kira contoh hasilnya seperti ini untuk data sintetik :



Well, kenapa ya harus dibuat smooth ? Setelah program ini saya buat jadi sedikit dapat pencerahan maksud dari Om Constable membuat Occam. Pernah mendengar Occam's Razor ?



Dari pernyataan Occam tersebut, maksud Om Constable membuat metode Inversi Occam ini (dan mungkin asal penamaan metode inversinya juga) adalah mencari solusi paling simpel lewat bentuk model yang smooth (meskipun tidak super fit dengan data pengukuran) untuk menghindari kasus overparameterized seperti pada postingan saya sebelumnya. Maksud saya disini bukan berarti metode lain salah, tetapi metode ini merupakan salah satu opsi untuk menghindari kasus "ke-lebay-an" hasil inversi yang sedang kita tidak harapkan pada suatu pengolahan

So that's it. Lain kali saya update dengan programnya, masih belum rapi nih.

See ya,
-L

Referensi:
[1] Constable, S.C., R.L. Parker, and C.G. Constable, 1987: Occam's Inversion: a practical algorithm for generating smooth models from EM sounding data, Geophysics, 52, pp. 289-300.
[2]  P.  R. McGILLIVRAY and D. W. OLDENBURG. 1990. Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem: a comparative study. Geophysical Prospecting.
[3]  Pei, Donghong, Ph.D. 2007. Modeling and inversion of dispersion curves of surface waves in shallow site investigations. UNIVERSITY OF NEVADA, RENO.

12 comments:

Event

Asisten IESO 2014 (lagi..)

6:51:00 AM Leo Cahya D 1 Comments

Hello digital world,




Kemarin tanggal 13-18 Maret 2014, Saya, Theo, dan Putri Portugal menjadi asisten praktek lagi untuk pelatihan olimpiade SMA yang rencananya akan maju di IESO (International Earth Science Olympiad) 2014 di US kalau tidak salah. Seperti tahun sebelumnya, lokasinya di Hotel University, bedanya perwakilan dosen dari geofisika cuma Pak Afif dan Pak Eddy karena Pak Herlan sibuk ke Jerman rumornya.

Jadwal kegiatan untuk sesi Geofisika seperti tahun sebelumnya, 3 hari kelas, 1 hari lapangan, dan 1 hari ujian seleksi. Saya kali ini dapat bagian mengajarkan metode refraksi dan geolistrik. Hal-hal yang diajarkan pun masih simpel, baru asumsi lapisan datar untuk refraksi, dipole-dipole untuk geolistrik. 

Sesi lapangan sesuai ide Pak Eddy, saya ajak anak-anak tersebut profiling 2D menggunakan dipole-dipole di depan Lab Geofisika yang teduh dan kebetulan Pohon Taloknya sedang berbuah hehe. Lintasan yang mereka buat hanya 20 meter dengan a=1 m dan mereka udah bilang capek, hehe yasudahlah namanya cuma hanya untuk pemahaman.

Selain kegiatan belajar, menurut saya makan malam merupakan acara yang menarik juga, pada saat itu saya sempat berkenalan dengan dosen dan mahasiswa Astronomi ITB, ramah-ramah orangnya.




Sepertinya itu saja yang bisa diceritakan. Well, semoga sukses semua deh buat mereka, Aamiin :)

See ya,

L

1 comments:

Programming

Inversi 1D Magnetotelurrik menggunakan metode VFSA

6:24:00 AM Leo Cahya D 7 Comments

Hello Digital World,


First of all, nonlinearity is a b*tch.
Yap, fungsi-fungsi asumsi dalam metode-metode geofisika itu memang nonlinear, saya baru tau maksud dari bercandaan "kalau orang geofisika ditanya 1+1 pasti bakal dijawab, 'Anda maunya berapa?' " setelah saya membuat beberapa program inversi seperti gravitasi dan MT1D ini.  

Jadi ceritanya saya merampungkan project mt1d dengan teman saya theo yang sebenarnya untuk suatu lomba dua tahun silam (tapi karena ada kegiatan lapangan batal deh). Acuan yang saya pakai menggunakan metode forward yang sama dengan buatan Pak Hendra Grandis*, sedangkan untuk inversinya saya menggunakan metode Very Fast Simulated Annealing (VFSA) besutan Ingber (1989).

Kenapa saya menggunakan VFSA?
Ya menurut saya dalam menyelesaikan masalah inversi perlu dilihat dulu berapa jumlah parameter model yang diinversikan dan lama proses perhitungan forwardnya. Jika parameter < 200 jenis (dipost selanjutnya mungkin akan saya jelaskan kenapa saya berpendapat seperti ini melalui experimen dengan problem simpel) dan waktu untuk forward sebuah model kurang dari 30 detik metode random walk semacam monte carlo dan simulated annealing opsi yang baik, dimana anda tak perlu repot-repot menghitung matriks jacobi dan stress apabila terjadi singularitas.

Hasilnya?
Error kecil tapi model beda dari hasil forward awal. Yup, super ultimate mega nonlinear deh. Bisa dilihat hasil 3 kali running untuk suatu data sintetik dengan model yang sama tetapi hasilnya berbeda.
(Jumlah parameter = 21 (11 nilai rho dan 10 ketebalan lapisan) )

Test 1

Test 2

Test 3

Baru ini, belum seismik haha.. Man gotta dream like never seen obstacles~

Berikut program matlabnya jika ingin dipakai :
forward: mt1dfor2.p
inversi : MT_1D_Inversill.m

*Referensi : H. Grandis, M. Menvielle, M. Roussignol. Bayesian inversion with Markov chains—I. The magnetotelluricone‐dimensional case

see ya
- L

7 comments: