Programming,
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.
Inversi 1D Occam Magnetotellurik
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.
- 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 :
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.
maaf mas mau tanya
ReplyDeletecara ngeliat jacobian benar atau ga gimana ya?
sya bikin yg forward difference dimatlab tpi pas digunain di inversi, hasil nya ga bagus
makasih mas
Untuk melihat benar atau tidaknya bisa dilihat isi matriks jacobi itu nol semua atau tidak (paling simpel menggunakan perintah imagesc() di matlab, umumnya tampilan gambar dari matriks jacobi ada garis2 diagonal ).
ReplyDeleteKesalahan program inversi bisa karena :
- jacobinya nol semua, solusinya tinggikan sedikit nilai h (pertubasi) pada forward difference. Umumnya paper-paper inversi yg menggunakan FD nilai h sekitar 0.5-5% dari nilai parameternya. Tapi ingat, makin tinggi nilai h, makin rendah akurasinya terhadap nilai analitik, jadi cari seminim mungkin. (sumber : "Göktürkler, 2010. Traveltime tomography of crosshole radar data without ray tracing")
- kalau menggunakan damped gauss newton atau metode yg ada parameter alpha(dampingnya), inversinya jalan, tapi nilai dampingnya kebesaran (silahkan cek buku Paremeter Estimationnya R.Aster Excersice 10.1untuk lebih jelasnya). Saran saya buat sederet nilai damping yg logaritmik dan coba yg paling cocok yg mana.
Maaf mas, mau tanya, maksud isi matriks jacobi itu nol bagaimana ya?, apakah determinannya yg nol?, atau diagonalnya atau bagaimana?, terima kasih
Deleteisi matriks jacobinya mbak, kadang kalau salah menghitung atau pada saat iterasi tertentu inversi tidak stabil matriksnya nol semua.
Deleteyang benar itu terdapat nilai yang susunannya diagonal (jika di imagesc pada matlab).
contoh :
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Trima kasih mas penjelasannya
ReplyDeleteSy gunain h 5%, nilai jacobiannya ga ada yg 0, cuma nilai diagonalnya semakin kebawah semakin kecil (salah ga?), terus faktor redaman memang dibuat banyak mas sampai prosesnya konvergen,
Kalo model dibuat homogen, kurva resistivitaa semu dan fasa (45) sudah benar, kalo model nya tidak homogen kurva resistivitas nya sudah hmpir fit, tapi fasa belum, dan model hasil inversinya salah.
Kalau pemilihan ketebalan waktu inversi bagaimana mas? (ketebalan yg digunakan waktu forward dipilih beda2).
Well, bisa dampingnya kebesaran, memang nonlinearitas problem MT, atau kesalahan pemrograman, infonya masih minim, saya belum bisa menyimpulkan solusi enaknya bagaimana. Begini saja mas coba membuat data sintetik dan inversi dari contoh animasi saya di atas, hasilnya bagaimana nanti bisa sharing. Pakai metode apa ya?
Deleteyup, sya udah kerjain pake occam juga
Deletealamat emailnya mas? nanti sya kirim sama coding juga, mungkin bisa dikoreksi, hehe
makasih mas..
flowchart saudara saja di share linknya bagaimana? tiap orang jalan pemikiran beda2, kalau dikoreksi dari kode gak selesai2 nanti hehe..
Deleteudh bisa mas, cuma ada beberapa nilai yang minus
Deletetrima kasih mas
Wah menarik nih, bisa share coding ny mas? :D
ReplyDeletewah belum sempet saya rapiin mbak, mending membaca bukunya Sen & Stoffa atau Richard Aster saja sudah banyak contoh2 aplikasinya disana :)
DeleteDear Mas Leo,
ReplyDeleteThank you on all inverse and forward model that Mas Leo was posting and almost 90 % is running well ( have tried).
If Mas Leo don't mind , mohon dikirimkan "Very Fast Simulated Annealing (VFSA) dan metode inversinya dari si Last & Kubik untuk saya coba lagi.
Salam Mas Leo
Dari Ari