Analisis Data Geofisika: Memahami Teori Inversi

Analisis Data Geofisika: Memahami Teori Inversi

Dr. Eng. Supriyanto, M.Sc Edisi I

Departemen Fisika-FMIPA Univeristas Indonesia 2007

Untuk Muflih Syamil, Hasan Azmi dan Farah Raihanah

Mottoku : Tenang, Kalem dan Percaya Diri

Kata Pengantar

Buku ini semula merupakan diktat perkuliahan mata kuliah Analisis Data Geofisika yang diberikan kepada mahasiswa Geofisika pada Departemen Fisika, Fakultas MIPA, Universitas Indonesia. Acuan utama untuk edisi perdana ini adalah buku Geophysical Data Analysis: Understanding

Inverse Problem Theory and Practice yang ditulis oleh Max A. Meju dan diterbitkan oleh Society of Exploration Geophysicists pada tahun 1994. Semoga keberadaan buku ini dapat membantu mahasiswa geofisika untuk memiliki kemampuan memformulasikan masalah, menyusun hipotesis, metode dan solusi sehingga mampu menyelesaikan masalah-masalah geofisika secara mandiri. Terima kasih yang tak terhingga ingin saya sampaikan kepada Dede Djuhana yang telah bersedia berbagi memberikan file format buku dalam LATEX sehingga tampilan buku ini menjadi jauh lebih baik. Depok, 15 Maret 2007 Supriyanto S.

v

Daftar Isi

Lembar Persembahan

i

Kata Pengantar

v

Daftar Isi

vii

Daftar Gambar

ix

Daftar Tabel

xi

1 Pendahuluan

1

1.1 Definisi dan Konsep Dasar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2 Proses geofisika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3 Eksplorasi geofisika dan inversi . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.4 Macam-macam data geofisika . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.5 Deskripsi proses geofisika: Model matematika . . . . . . . . . . . . . . . . . . . .

4

1.6 Diskritisasi dan linearisasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2 Formulasi Masalah Inversi

9

2.1 Klasifikasi masalah inversi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2 Inversi Model Garis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.3 Inversi Model Parabola . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.4 Inversi Model Bidang . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

2.5 Contoh aplikasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

2.5.1 Menghitung gravitasi di planet X . . . . . . . . . . . . . . . . . . . . . . .

24

2.5.2 Analisa data seismik pada reflektor tunggal horizontal . . . . . . . . . . .

28

2.5.3 Analisa data seismik pada reflektor tunggal miring . . . . . . . . . . . . .

29

2.6 Kesimpulan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3 Penyelesaian Masalah Overdetermined

31

3.1 Regresi linear sederhana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

3.2 Metode least square . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

3.3 Aplikasi regresi linear pada analisis data seismik refraksi . . . . . . . . . . . . . .

34

3.4 Regresi linear dengan pendekatan matriks . . . . . . . . . . . . . . . . . . . . . .

36

3.5 Soal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

vii

viii 4 Constrained Linear Least Squares Inversion

41

4.1 Inversi dengan informasi awal . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

4.1.1 Memformulasikan persamaan terkonstrain . . . . . . . . . . . . . . . . . .

43

4.1.2 Contoh aplikasi inversi terkonstrain . . . . . . . . . . . . . . . . . . . . . .

43

4.2 Inversi dengan Smoothness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

4.2.1 Formulasi masalah . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

4.2.2 Solusi masalah . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

Daftar Pustaka

49

Indeks

51

Daftar Gambar

1.1 Alur pemodelan inversi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2 Alur pemodelan forward . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3 Alur eksperimen lapangan dan eksperimen laboratorium . . . . . . . . . . . . . .

4

1.4 Konfigurasi elektroda pada metode Schlumberger . . . . . . . . . . . . . . . . . .

6

2.1 Data observasi perubahan suhu terhadap kedalaman dari permukaan tanah . . .

10

2.2 Hasil inversi atas data observasi perubahan suhu terhadap kedalaman . . . . . .

14

2.3 Data observasi perubahan suhu terhadap kedalaman dari permukaan tanah . . .

15

2.4 Hasil inversi atas data observasi perubahan suhu terhadap kedalaman. Tanda titik merah adalah data observasi sementara kurva biru adalah kurva hasil inversi

18

2.5 Data observasi perubahan nilai terhadap posisi koordinat x dan y . . . . . . . . .

20

2.6 Hasil inversi berbentuk bidang. Saya sebut saja sebagai bidang inversi . . . . . .

23

2.7 Ini adalah hasil inversi yang sama, hanya saja bidang inversi diputar (di-rotasi) pada suatu arah sehingga kita bisa saksikan dan pastikan bahwa sebaran titik data observasi berada di-dekat permukaan bidang inversi . . . . . . . . . . . . .

23

2.8 Grafik data pengukuran gerak batu . . . . . . . . . . . . . . . . . . . . . . . . . .

26

2.9 Grafik hasil inversi parabola . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

2.10 Reflektor mendatar pada kedalaman z. Kecepatan gelombang v dianggap konstan. S adalah sumber gelombang seismik dan R adalah penerima gelombang seismik. Jarak antara S dan R disebut offset (x). Sementara garis merah yang ada panahnya adalah lintasan gelombang seismik. . . . . . . . . . . . . . . . . .

29

2.11 Reflektor miring dengan sudut kemiringan sebesar α. Kecepatan gelombang v dianggap konstan. S adalah sumber gelombang seismik dan R adalah penerima gelombang seismik. Jarak antara S dan R disebut offset (x). Sementara garis merah yang ada panahnya adalah lintasan gelombang seismik. . . . . . . . . . .

30

3.1 Hasil plotting data observasi dalam sumbu-x dan sumbu-y . . . . . . . . . . . . .

32

3.2 Contoh solusi regresi linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

ix

Daftar Tabel

2.1 Data temperatur bawah permukaan tanah terhadap kedalaman . . . . . . . . . .

10

2.2 Data temperatur bawah permukaan tanah terhadap kedalaman . . . . . . . . . .

15

2.3 Distribusi suatu nilai dalam area tertentu . . . . . . . . . . . . . . . . . . . . . . .

20

2.4 Data ketinggian terhadap waktu dari planet X . . . . . . . . . . . . . . . . . . . .

25

2.5 Data variasi offset (x) dan travel time (t) . . . . . . . . . . . . . . . . . . . . . . .

29

2.6 Data variasi offset (x) dan travel time (t) . . . . . . . . . . . . . . . . . . . . . . .

30

3.1 Contoh data observasi yang dapat diolah oleh least squares . . . . . . . . . . . .

32

3.2 Data seismik refraksi: waktu-datang gelombang, ti , dan jarak antara source dan receiver atau jarak offset, xi

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

3.3 Data seismik refraksi: waktu-datang gelombang (ti ) pada empat posisi geophone (xi ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.1 Data seismik refraksi: waktu-datang gelombang (ti ) pada empat posisi geophone (xi ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

45

Bab 1

Pendahuluan

1.1 Definisi dan Konsep Dasar Dalam geofisika, kegiatan pengukuran lapangan selalu dilakukan berdasarkan prosedur yang sudah ditentukan. Kemudian, hasil pengukuran dicatat dan disajikan dalam bentuk tabel angkaangka pengukuran. Hasil pengukuran tersebut sudah barang tentu sangat tergantung pada kondisi dan sifat fisis batuan bawah permukaan. Tabel angka-angka itu selanjutnya disebut data observasi atau juga biasa disebut data lapangan . Kita berharap data eksperimen dapat memberi informasi sebanyak-banyaknya, tidak sekedar mengenai sifat fisis batuan saja, melainkan juga kondisi geometri batuan bawah permukaan dan posisi kedalaman batuan tersebut. Informasi itu hanya bisa kita dapat bila kita mengetahui hubungan antara sifat fisis batuan tersebut dan data observasinya. Penghubung dari keduanya hampir selalu berupa persamaan matematika atau kita menyebutnya sebagai model matematika. Maka dengan berdasarkan model matematika itulah, kita bisa mengekstrak parameter fisis batuan dari data observasi. Proses ini disebut proses inversi atau istilah asingnya disebut inverse modelling, lihat Gambar 2.9. Sementara proses kebalikannya dimana kita ingin memperoleh data prediksi hasil pengukuran berdasarkan parameter fisis yang sudah diketahui, maka proses ini disebut proses forward atau forward modelling, lihat Gambar 1.2. Proses inversi adalah suatu proses pengolahan data lapangan yang melibatkan teknik penyelesaian matematika dan statistik untuk mendapatkan informasi yang berguna mengenai distribusi sifat fisis bawah permukaan. Di dalam proses inversi, kita melakukan analisis terhadap data lapangan dengan cara melakukan curve fitting (pencocokan kurva) antara model matematika dan data lapangan. Tujuan dari proses inversi adalah untuk mengestimasi parameter fisis batuan yang tidak diketahui sebelumnya (unknown parameter). Proses inversi terbagi dalam level-level tertentu mulai dari yang paling sederhana seperti fitting garis untuk data seismik refraksi sampai kepada level yang rumit seperti tomografi akustik dan matching (pencocokan) kurva resistivity yang multidimensi. Contoh problem inversi dalam bidang geofisika adalah 1. Penentuan struktur bawah tanah 2. Estimasi parameter-parameter bahan tambang 1

BAB 1. PENDAHULUAN

2

Gambar 1.1: Alur pemodelan inversi

Gambar 1.2: Alur pemodelan forward

1.2. PROSES GEOFISIKA

3

3. Estimasi parameter-parameter akumulasi sumber energi 4. Penentuan lokasi gempa bumi berdasarkan waktu gelombang datang 5. Pemodelan respon lithospere untuk mengamati proses sedimentasi 6. Analisis sumur bor pada hidrogeologi

1.2 Proses geofisika Perambatan gelombang seismik, perambatan gelombang elektromagnetik di bawah tanah dan juga aliran muatan (arus listrik) ataupun arus fluida pada batuan berpori adalah contoh-contoh proses geofisika. Data lapangan tak lain merupakan refleksi dari kompleksitas sistem geofisika yang sedang diamati, yang dikontrol oleh distribusi parameter fisis batuan berikut struktur geologinya.

1.3 Eksplorasi geofisika dan inversi Tujuan utama dari kegiatan eksplorasi geofisika adalah untuk membuat model bawah permukaan bumi dengan mengandalkan data lapangan yang diukur bisa pada permukaan bumi atau di bawah permukaan bumi atau bisa juga di atas permukaan bumi dari ketinggian tertentu. Untuk mencapai tujuan ini, idealnya kegiatan survey atau pengukuran harus dilakukan secara terus menerus, berkelanjutan dan terintegrasi menggunakan sejumlah ragam metode geofisika. Seringkali –bahkan hampir pasti– terjadi beberapa kendala akan muncul dan tak bisa dihindari, seperti kehadiran noise pada data yang diukur. Ada juga kendala ketidaklengkapan data atau malah kurang alias tidak cukup. Namun demikian, dengan analisis data yang paling mungkin, kita berupaya memperoleh informasi yang relatif valid berdasarkan keterbatasan data yang kita miliki. Dalam melakukan analisis, sejumlah informasi mengenai kegiatan akuisisi data juga diperlukan, antara lain: berapakah nilai sampling rate yang optimal? Berapa jumlah data yang diperlukan? Berapa tingkat akurasi yang diinginkan? Selanjutnya –masih bagian dari proses analisis– model matematika yang cocok mesti ditentukan yang mana akan berperan ketika menghubungkan antara data lapangan dan distribusi parameter fisis yang hendak dicari. Setelah proses analisis dilalui, langkah berikutnya adalah membuat model bawah permukaan yang nantinya akan menjadi modal dasar interpretasi. Ujung dari rangkaian proses ini adalah penentuan lokasi pemboran untuk mengangkat sumber daya alam bahan tambang/mineral dan

oil-gas ke permukaan. Kesalahan penentuan lokasi berdampak langsung pada kerugian meteril yang besar dan waktu yang terbuang percuma. Dari sini terlihat betapa pentingnya proses analisis apalagi bila segala keputusan diambil berdasarkan data eksperimen.

1.4 Macam-macam data geofisika Data geofisika bisa diperoleh dari pengukuran di lapangan atau bisa juga dari pengukuran di laboratorium. Gambar 1.3 memperlihatkan alur pengambilan data dari masing-masing pen-

BAB 1. PENDAHULUAN

4

Gambar 1.3: Alur eksperimen lapangan dan eksperimen laboratorium gukuran. Pada pengukuran lapangan, data geofisika yang terukur antara lain bisa berupa densitas, kecepatan gelombang seismik, modulus bulk, hambatan jenis batuan, permeabilitas batuan, suseptibilitas magnet dan lain sebagainya yang termasuk dalam besaran fisis sebagai karakteristik bawah permukaan bumi. Pada pengukuran di laboratorium, model lapisan bumi ataupun keberadaan anomali dalam skala kecil dapat dibuat dan diukur respon-nya sebagai data geofisika. Diharapkan hasil uji laboratorium tersebut bisa mewakili kondisi lapangan yang sesungguhnya yang dimensinya jauh lebih besar. Jika suatu pengukuran diulang berkali-kali, entah itu di lapangan maupun di laboratorium, seringkali kita temukan hasil pengukuran yang berubah-ubah, walaupun dengan variasi yang bisa ditolerir. Variasi ini umumnya disebabkan oleh kesalahan instrumen pengukuran (instru-

mental error) atau bisa juga dikarenakan kesalahan manusia (human error). Seluruh variasi ini bila di-plot kedalam histogram akan membentuk distribusi probabilistik.

1.5 Deskripsi proses geofisika: Model matematika Seluruh proses geofisika dapat dideskripsikan secara matematika. Sebagaimana yang telah disebutkan diawal, suatu formulasi yang bisa menjelaskan sistem geofisika disebut model. Namun perlu ditekankan juga bahwa istilah model memiliki ragam konotasi berbeda di kalangan geosaintis. Misalnya, orang geologi kerapkali menggunakan istilah model konseptual, atau istilah model fisik yang digunakan untuk menyebutkan hasil laboratorium, atau dalam catatan ini kita menggunakan istilah model matematika yang merupakan istilah umum dikalangan para ahli

1.6. DISKRITISASI DAN LINEARISASI

5

geofisika. Kebanyakan proses geofisika dapat dideskripsikan oleh persamaan integral berbentuk di =

z

Z

Ki (z)p(z)dz

(1.1)

0

dimana di adalah respon atau data yang terukur, p(z) adalah suatu fungsi yang berkaitan dengan parameter fisis yang hendak dicari (misalnya: hambatan jenis, densitas, kecepatan, dan lain-lain) yang selanjutnya disebut parameter model, dan Ki disebut data kernel. Data kernel menjelaskan hubungan antara data dan parameter model p(z). Parameter model (misalnya kecepatan, resistivitas dan densitas) bisa jadi merupakan fungsi yang kontinyu terhadap jarak atau posisi. Sebagai contoh, waktu tempuh t antara sumber gelombang seismik dengan penerimanya sepanjang lintasan L dalam medium, yang distribusi kecepatan gelombangnya kontinyu v(x, z), ditentukan oleh

Z

t=

L

1 dl v(x, z)

(1.2)

Deskripsi matematika terhadap sistem geofisika seperti contoh di atas disebut forward modelling. forward modelling digunakan untuk memprediksi data simulasi berdasarkan hipotesa kondisi bawah permukaan. Data simulasi tersebut biasanya dinamakan data teoritik atau data sintetik atau data prediksi atau data kalkulasi. Cara seperti ini disebut pendekatan forward atau lebih dikenal sebagai pemodelan forward (Gambar 1.2).

1.6 Diskritisasi dan linearisasi Dalam banyak kasus, model bumi selalu fungsi kontinyu terhadap jarak dan kedalaman. Mari kita ambil kasus massa dan momen inersia bumi. Keduanya terkait dengan densitas bawah permukaan sesuai rumus-rumus berikut M assa = 4π

Z

R

r2 ρ(r)dr

(1.3)

0

M oment inersia =

8π 3

Z

R

r4 ρ(r)dr

(1.4)

0

dimana R adalah jejari bumi dan ρ(r) merupakan fungsi densitas terhadap jarak r. ρ(r) juga berhubungan dengan p(z) pada persamaan (1.1). Persamaan (1.4 dan 1.4) dapat dinyatakan dalam formulasi yang lebih umum yaitu di =

Z

R

Ki (r)p(r)dr

(1.5)

0

sama persis dengan persamaan (1.1). Integral ini relatif mudah dievaluasi secara komputasi dengan matematika diskrit. Pendekatan komputasi memungkinkan kita untuk menyederhanakan ρ(r)dr menjadi m, sementara Ki menjadi Gi sehingga persamaan (1.5) dapat dinyatakan sebagai di =

X

Gij mj

(1.6)

BAB 1. PENDAHULUAN

6

Gambar 1.4: Konfigurasi elektroda pada metode Schlumberger Ini adalah bentuk diskritisasi. Secara umum, memang pada kenyataannya ketika melakukan eksperimen di lapangan, data pengukuran maupun paremeter model selalu dibatasi pada interval tertentu. Kita sering berasumsi bahwa bawah permukaan bumi terdiri dari lapisan-lapisan yang masing-masing memiliki sifat fisis atau parameter fisis p(z) yang seragam. Misalnya lapisan tertentu memiliki densitas sekian dan ketebalan sekian. Langkah praktis ini yang terkesan menyederhanakan obyek lapangan disebut langkah parameterisasi. Dalam kuliah ini, kita akan selalu memandang model yang diskrit dan juga parameter yang diskrit daripada model dan paremeter yang kontinyu. Sehingga proses inversi yang akan kita lakukan disebut sebagai teori inversi diskrit dan bukan teori inversi kontinyu. Dalam bentuk diskrit, persamaan (1.2) bisa dinyatakan sebagai ti =

p X Lij j=1

(1.7)

vj

Perlu dicatat disini bahwa waktu tempuh t tidak berbanding lurus dengan parameter model v, melainkan berbanding terbalik. Hubungan ini dinamakan non-linear terhadap v. Namun demikian, jika kita mendefinisikan parameter model c = 1/v, dimana c adalah slowness gelombang seismik, maka problem ini bisa dinyatakan sebagai p X

ti =

Lij cj

(1.8)

j=1

Hubungan ini disebut linear. Persamaan memenuhi bentuk d = Gm. Operasi transformasi seperti itu dinamakan linearisasi parameter. Dan proses menuju kesana dinamakan linearisasi. Sekarang mari kita lihat problem dari pengukuran resistivitas semu dengan metode Schlumberger untuk mengamati lapisan bawah permukaan yang diasumsikan terdiri dari dua lapisan. Formula model yang diturunkan oleh Parasnis, 1986 adalah

2

ρa (L) = ρ1 1 + 2L

Z

0

∞

K(λ)J1 (λL)λdλ

(1.9)

1.6. DISKRITISASI DAN LINEARISASI

7

dimana L = AB/2 adalah jarak masing-masing elektroda terhadap titik tengah, J1 adalah fungsi Bessel orde 1 dan K(λ) adalah fungsi parameter (resistivitas masing-masing lapisan yaitu ρ1 dan ρ2 serta ketebalan lapisan paling atas t) dari sistem yang kita asumsikan. K(λ) dinyatakan sebagai (−2λt)

K(λ) =

−k1,2

(−2λt)

1 + −k1,2

dimana k1,2 =

ρ1 − ρ2 ρ1 + ρ2

Kita bisa lihat bahwa persamaan (1.9) tidak bisa didekati dengan d = Gm sebagaimana yang dilakukan pada persamaan (1.2). Oleh karena itu persamaan resistivitas semu di atas disebut highly non-linear.

Bab 2

Formulasi Masalah Inversi

2.1 Klasifikasi masalah inversi Dalam masalah inversi, kita selalu berhubungan dengan parameter model (M ) dan data (N ); yang mana jumlah dari masing-masing akan menentukan klasifikasi permasalahan inversi dan cara penyelesaiannya. Bila jumlah model parameter lebih sedikit dibandingkan data observasi (M < N ), maka permasalahan inversi ini disebut overdetermined. Umumnya masalah ini diselesaikan menggunakan pencocokan (best fit ) terhadap data observasi. Dalam kondisi yang lain dimana jumlah parameter yang ingin dicari (M ) lebih banyak dari pada jumlah datanya (N ), maka masalah inversi ini disebut underdetermined. Dalam kasus ini terdapat sekian banyak model yang dapat sesuai kondisi datanya. Inilah yang disebut dengan masalah non-uniqness. Bagaimana cara untuk mendapatkan model yang paling mendekati kondisi bawah permukaan? Menurut Meju, 1994 persoalan ini bisa diselesaikan dengan model yang parameternya berbentuk fungsi kontinyu terhadap posisi. Kasus yang terakhir adalah ketika jumlah data sama atau hampir sama dengan jumlah parameter. Ini disebut evendetermined. Pada kasus ini model yang paling sederhana dapat diperoleh menggunakan metode inversi langsung. Pada bab ini, saya mencoba menyajikan dasar teknik inversi yang diaplikasikan pada model garis, model parabola dan model bidang. Uraian aplikasi tersebut diawali dari ketersediaan data observasi, lalu sejumlah parameter model (unknown parameter) mesti dicari dengan teknik inversi. Mari kita mulai dari model garis.

2.2 Inversi Model Garis Secara teori, variasi temperatur bawah permukaan akan semakin meningkat ketika temperatur tersebut diukur semakin kedalam permukaan bumi. Misalnya telah dilakukan sebanyak sepuluh kali (N = 10) pengukuran temperatur (Ti ) pada kedalaman yang berbeda beda (zi ) sebagaimana ditunjukan datanya pada Tabel 2.6.

9

BAB 2. FORMULASI MASALAH INVERSI

10

Tabel 2.1: Data temperatur bawah permukaan tanah terhadap kedalaman Pengukuran ke-i Kedalaman (m) Temperatur (O C) 1 z1 = 5 T1 = 35.4 2 z2 = 16 T2 = 50.1 3 z3 = 25 T3 = 77.3 4 z4 = 40 T4 = 92.3 5 z5 = 50 T5 = 137.6 6 z6 = 60 T6 = 147.0 7 z7 = 70 T7 = 180.8 8 z8 = 80 T8 = 182.7 9 z9 = 90 T9 = 188.5 10 z10 = 100 T10 = 223.2

Variasi Suhu vs Kedalaman 250

Suhu (Celcius)

200

150

100

50

0

0

20

40 60 Kedalaman (m)

80

100

Gambar 2.1: Data observasi perubahan suhu terhadap kedalaman dari permukaan tanah Source code untuk menggambar grafik tersebut dalam Matlab adalah 1 2 3

clc clear close

4 5 6 7

% Data observasi z = [5 16 25 40 50 60 70 80 90 100]; T = [35.4 50.1 77.3 92.3 137.6 147.0 180.8 182.7 188.5 223.2];

8 9 10 11 12 13 14

% Plot data observasi plot(z,T,’*r’); grid; xlabel(’Kedalaman (m)’); ylabel(’Suhu (Celcius)’); title(’\fontsize{14} Variasi Suhu vs Kedalaman’);

2.2. INVERSI MODEL GARIS

11

Lalu kita berasumsi bahwa variasi temperatur terhadap kedalaman ditentukan oleh rumus berikut ini: m1 + m2 zi = Ti

(2.1)

dimana m1 dan m2 adalah konstanta-konstanta yang akan dicari. Rumus di atas disebut model matematika. Sedangkan m1 dan m2 disebut parameter model atau biasa juga disebut unknown parameter. Pada model matematika di atas terdapat dua buah parameter model, (M = 2). Sementara jumlah data observasi ada empat, (N = 10), yaitu nilai-nilai kedalaman, zi , dan temperatur, Ti . Berdasarkan model tersebut, kita bisa menyatakan temperatur dan kedalaman masing-masing sebagai berikut: m1 + m2 z1 = T1 m1 + m2 z2 = T2 m1 + m2 z3 = T3 m1 + m2 z4 = T4 m1 + m2 z5 = T5 m1 + m2 z6 = T6 m1 + m2 z7 = T7 m1 + m2 z8 = T8 m1 + m2 z9 = T9 m1 + m2 z10 = T10 Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini:                     





T1

1

z1

1

   T2 z2     T  z3   3    z4  " #  T4  T  z5  m1  5 =   T6 m z6  2    T  z7   7    T8  z8    T  z9   9 T10 z10

1 1 1 1 1 1 1 1

                    

(2.2)

Lalu ditulis secara singkat Gm = d

(2.3)

dimana d adalah data yang dinyatakan dalam vektor kolom, m adalah model parameter, juga dinyatakan dalam vektor kolom, dan G disebut matrik kernel. Lantas bagaimana cara menda-


12

patkan nilai m1 dan m2 pada vektor kolom m? Manipulasi1 berikut ini bisa menjawabnya GT Gm = GT d

(2.4)

dimana T disini maksudnya adalah tanda transpos matrik. Selanjutnya, untuk mendapatkan elemen-elemen m, diperlukan langkah-langkah perhitungan berikut ini: GT Gm = GT d [GT G]−1 GT Gm = [GT G]−1 GT d m = [GT G]−1 GT d

(2.5)

1. Tentukan transpos dari matrik kernel, yaitu GT 

          G=         



1

z1

1

 z2   z3    z4   z5    z6   z7    z8   z9   z10

1 1 1 1 1 1 1 1

T

⇒

G =

"

1

1

1

1

1

1

1

1

1

1

z1 z2 z3 z4 z5 z6 z7 z8 z9 z10

#

2. Lakukan perkalian matriks GT G 

GT G =

"

1

1

1

1

1

1

1

1

1

z1 z2 z3 z4 z5 z6 z7 z8 z9

       #  1   z10         



1

z1

1

 z2   z3    z4   " P # N zi z5   = P P 2 zi zi z6    z7   z8   z9   z10

1 1 1 1 1 1 1 1

dimana N = 10 atau sesuai dengan jumlah data observasi; sementara i = 1, 2, 3, ..., 10.

1

Matrik G biasanya tidak berbentuk bujursangkar. Akibatnya tidak bisa dihitung nilai invers-nya. Dengan mengalikan matrik G dan transpose matrik G, maka akan diperoleh matrik bujursangkar

2.2. INVERSI MODEL GARIS

13

3. Kemudian tentukan pula GT d 

GT d =

"

1

1

1

1

1

1

1

1

1

z1 z2 z3 z4 z5 z6 z7 z8 z9

       #  1   z10         

T1



 T2   T3    T4   " P # T5  Ti  = P T6  z i Ti   T7   T8   T9   T10

4. Dengan menggunakan hasil dari langkah 2 dan langkah 3, maka persamaan (2.4) dapat dinyatakan sebagai GT Gm = GT d "

# " P # P #" m1 Ti N zi = P P P 2 m2 z i Ti zi zi

(2.6)

Berdasarkan data observasi pada tabel di atas, diperoleh m = [GT G]−1 GT d

# P #−1 " P Ti N zi = P 2 P P m2 zi z i Ti zi # " #−1 " # " 10 536 1314 m1 = 536 38006 88855 m2 "

m1

#

"

Operasi matriks di atas akan menghasilkan nilai m1 = 24.9 dan m2 = 2.0. Perintah di Matlab untuk menghitung elemen-elemen m, yaitu m=inv(G’*G)*G’*d Secara lebih lengkap, source code Matlab untuk melakukan inversi data observasi adalah 1 2 3

clc clear close

4 5 6 7 8

% Data observasi z = [5 16 25 40 50 60 70 80 90 100]; T = [35.4 50.1 77.3 92.3 137.6 147.0 180.8 182.7 188.5 223.2];

(2.7)


14 9 10 11 12 13 14

% Plot data observasi plot(z,T,’*r’); grid; xlabel(’Kedalaman (m)’); ylabel(’Suhu (Celcius)’); title(’\fontsize{14} Variasi Suhu vs Kedalaman’);

15 16 17 18 19 20 21 22

% Membentuk matrik kernel G dan vektor d n = length(z); for k = 1:n G(k,1) = 1; G(k,2) = z(k); end d = T’;

23 24 25

% Perhitungan inversi dengan general least-squares m = inv(G’*G)*G’*d;

26

28 29 30 31

% Plot hasil inversi (berupa garis least-squares) hold on; zz = 0:0.5:z(n); TT = m(1) + m(2)*zz; plot(zz,TT);


200

Suhu (Celcius)

27

150

100

50

0

0

20

40 60 Kedalaman (m)

80

100

Gambar 2.2: Hasil inversi atas data observasi perubahan suhu terhadap kedalaman

Demikianlah contoh aplikasi teknik inversi untuk menyelesaikan persoalan model garis. Anda bisa mengaplikasikan pada kasus lain, dengan syarat kasus yang anda tangani memiliki bentuk model yang sama dengan yang telah dikerjakan pada catatan ini, yaitu model garis: y = m1 + m2 x. Selanjutnya mari kita pelajari inversi model parabola.

2.3. INVERSI MODEL PARABOLA

15

2.3 Inversi Model Parabola Kembali kita ambil contoh variasi temperatur terhadap kedalaman dengan sedikit modifikasi data. Misalnya telah dilakukan sebanyak delapan kali (N = 8) pengukuran temperatur (Ti ) pada kedalaman yang berbeda beda (zi ). Tabel pengukuran yang diperoleh adalah: Data observasi Tabel 2.2: Data temperatur bawah permukaan tanah terhadap kedalaman Pengukuran ke-i Kedalaman (m) Temperatur (O C) 1 z1 = 5 T1 = 20, 8 2 z2 = 8 T2 = 22, 6 3 z3 = 14 T3 = 25, 3 4 z4 = 21 T4 = 32, 7 5 z5 = 30 T5 = 41, 5 6 z6 = 36 T6 = 48, 2 7 z7 = 45 T7 = 63, 7 8 z8 = 60 T8 = 74, 6 tersebut selanjutnya di-plot ke dalam grafik variasi suhu terhadap kedalaman.


70

Suhu (Celcius)

60

50

40

30

20

5

10

15

20

25 30 Kedalaman (m)

35

40

45

50

Gambar 2.3: Data observasi perubahan suhu terhadap kedalaman dari permukaan tanah Lalu kita berasumsi bahwa variasi temperatur terhadap kedalaman memenuhi model matematika berikut ini: m1 + m2 zi + m3 zi2 = Ti

(2.8)

dimana m1 , m2 dan m3 adalah unknown parameter. Jadi pada model di atas terdapat tiga buah model parameter, (M = 3). Adapun yang berlaku sebagai data adalah nilai-nilai temperatur T1 , T2 ,..., dan T8 . Berdasarkan model tersebut, kita bisa menyatakan temperatur dan kedalaman sebagai sistem persamaan simultan yang terdiri atas 8 persamaan (sesuai dengan jumlah data


16 observasi):

m1 + m2 z1 + m3 z12 = T1 m1 + m2 z2 + m3 z22 = T2 m1 + m2 z3 + m3 z32 = T3 m1 + m2 z4 + m3 z42 = T4 m1 + m2 z5 + m3 z52 = T5 m1 + m2 z6 + m3 z62 = T6 m1 + m2 z7 + m3 z72 = T7 m1 + m2 z8 + m3 z82 = T8 Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini:                

1 z1 z12





T1



    T2  1 z2 z22        1 z3 z32   T3      m1 T4  1 z4 z42        m2  =   1 z5 z52   T5      m 3  T6  1 z6 z62      T  1 z7 z72   7   T8 1 z8 z82

(2.9)


(2.10)

dimana d adalah data yang dinyatakan dalam vektor kolom, m adalah model parameter, juga dinyatakan dalam vektor kolom, dan G disebut matrik kernel. Lantas bagaimana cara mendapatkan nilai m1 , m2 dan m3 pada vektor kolom m? Manipulasi berikut ini bisa menjawabnya GT Gm = GT d

(2.11)

dimana t disini maksudnya adalah tanda transpos matrik. Selanjutnya, untuk mendapatkan elemen-elemen m, diperlukan langkah-langkah perhitungan berikut ini:

2.3. INVERSI MODEL PARABOLA

17

1. Tentukan transpos dari matrik kernel, yaitu GT 

       G=       

1 z1 z12



 1 z2 z22   1 z3 z32    2 1 z4 z4   1 z5 z52    2 1 z6 z6   1 z7 z72   1 z8 z82

⇒





1

1

1

1

1

1

1

1

 GT =  z1

z2

z3

z4

z5

z6

z7

 z8  z82

z12 z22 z32 z42 z52 z62 z72

2. Tentukan GT G 



1 1 1 1 1 1 1  G G =  z1 z2 z3 z4 z5 z6 z7 z12 z22 z32 z42 z52 z62 z72 T

     1   z8     z82     

1 z1 z12



 1 z2 z22   1 z3 z32  P P 2    N zi zi  2 1 z4 z4  P P P   2 = zi zi zi3  1 z5 z52  P P P   zi2 zi3 zi4 2 1 z6 z6   1 z7 z72   2 1 z8 z8

dimana N = 8 dan i = 1, 2, 3, ..., 8.

3. Kemudian tentukan pula GT d 



1 1 1 1 1 1 1  G d =  z1 z2 z3 z4 z5 z6 z7 z12 z22 z32 z42 z52 z62 z72 T

     1   z8     z82     

T1



 T2    T3    P Ti  T4  P  = z i Ti    T5  P 2  z i Ti T6   T7   T8

4. Sekarang persamaan (2.16) dapat dinyatakan sebagai    P P P 2  N zi zi m1 Ti P 2 P 3   P    P zi zi zi   m2  =  z i Ti   P 2 P 3 P 4 P 2 zi zi zi m3 z i Ti 

(2.12)


18

Berdasarkan data observasi pada tabel di atas, diperoleh 

8

219

8547



m1





349, 89



     8547 393423   m2  =  12894, 81   219 8547 393423 19787859 m3 594915, 33

(2.13)

Matlab telah menyediakan sebuah baris perintah untuk menghitung elemen-elemen m, yaitu m=inv(G’*G)*G’*d Sehingga operasi matriks di atas akan menghasilkan nilai m1 = 21, m2 = 0, 05 dan m3 = 0, 02. Gabungan antara data observasi dan kurva hasil inversi diperlihatkan oleh Gambar 2.4 Script


70

Suhu (Celcius)

60

50

40

30

20

10

0

10

20 30 Kedalaman (m)

40

50

Gambar 2.4: Hasil inversi atas data observasi perubahan suhu terhadap kedalaman. Tanda titik merah adalah data observasi sementara kurva biru adalah kurva hasil inversi Matlab yang lengkap untuk tujuan ini adalah 1 2 3

clc clear close

4 5 6 7

% Data observasi z = [5 8 14 21 30 36 45 50]; T = [20.8 22.6 25.3 32.7 41.5 48.2 63.7 74.6];

8 9 10

% Plot data observasi plot(z,T,’*r’);

2.3. INVERSI MODEL PARABOLA 11 12 13 14

19

grid; xlabel(’Kedalaman (m)’); ylabel(’Suhu (Celcius)’); title(’\fontsize{14} Variasi Suhu vs Kedalaman’);

15 16 17 18 19 20 21 22 23

% Membentuk matrik kernel G dan vektor d n = length(z); for k = 1:n G(k,1) = 1; G(k,2) = z(k); G(k,3) = z(k).^2; end d = T’;

24 25 26


27 28 29 30 31 32

% Plot hasil inversi (berupa garis least-squares) hold on; zz = 0:0.5:z(n); TT = m(1) + m(2)*zz + m(3)*zz.^2; plot(zz,TT);

Demikianlah contoh aplikasi teknik inversi untuk menyelesaikan persoalan model parabola. Anda bisa mengaplikasikan pada kasus lain, dengan syarat kasus yang anda tangani memiliki bentuk model yang sama dengan yang telah dikerjakan pada catatan ini, yaitu model garis: y = m1 + m2 x + +m3 x2 . Selanjutnya mari kita pelajari inversi model bidang atau model 2dimensi (2-D).


20

2.4 Inversi Model Bidang Misalnya telah dilakukan sebanyak sepuluh kali (N = 10) pengukuran suatu nilai parameter pada suatu area sebagaimana ditunjukan datanya pada Tabel 2.3. Tabel 2.3: Distribusi suatu nilai dalam area tertentu Pengukuran ke-i X (m) Y (m) Nilai 1 2 3 10,6 2 5 6 23,5 3 7 2 27,3 4 4 7 20,8 5 1 8 11,1 6 3 9 18,9 7 6 4 25,4 8 9 1 33,5 9 8 5 33,2 10 4 5 24,1

Sebaran nilai terhadap X dan Y

35 30

Nilai

25 20 15 10 10 10 8

5

6 4

Y (m)

0

2 0

X (m)

Gambar 2.5: Data observasi perubahan nilai terhadap posisi koordinat x dan y Source code untuk menggambar grafik tersebut dalam Matlab adalah 1 2 3

clc clear close

4 5 6 7 8

% Data observasi x = [2 5 7 4 1 3 6 9 8 4]; y = [3 6 2 7 8 9 4 1 5 5]; nilai = [10.6 23.5 27.3 20.8 11.1 18.9 25.4 33.5 33.2 24.1];

2.4. INVERSI MODEL BIDANG

21

9 10 11 12 13 14 15 16

% Plot data observasi plot3(x,y,nilai,’*r’); grid; xlabel(’X (m)’); ylabel(’Y (m)’); zlabel(’Nilai’); title(’\fontsize{14} Sebaran nilai terhadap X dan Y’);

Model matematika untuk 2-dimensi berikut ini digunakan untuk analisa data tersebut: m1 + m2 xi + m3 yi = di

(2.14)

dimana m1 , m2 dan m3 merupakan unknown parameter yang akan dicari. Adapun yang berlaku sebagai data adalah d1 , d2 , d3 , ..., dN . Berdasarkan model matematika tersebut, kita bisa nyatakan m1 + m2 x1 + m3 y1 = d1 m1 + m2 x2 + m3 y2 = d2 m1 + m2 x3 + m3 y3 = d3 .. .. .. .. .. . . . . . m1 + m2 xN + m3 yN = dN Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini:         







1

x1

y1

1

x2

1 .. .

x3 .. .

     d2  y2  m 1        y3    m2  =  d3     ..   ...  .  m3   yN dN

1 xN

d1


(2.15)

dimana d adalah data yang dinyatakan dalam vektor kolom, m adalah unknown parameter, juga dinyatakan dalam vektor kolom, dan G disebut matrik kernel. Lantas bagaimana cara mendapatkan nilai m1 , m2 dan m3 pada vektor kolom m? Sama seperti sebelumnya, kita harus membuat persamaan matriks berikut ini GT Gm = GT d

(2.16)

dimana T disini maksudnya adalah tanda transpos matrik. Selanjutnya, untuk mendapatkan elemen-elemen m, diperlukan langkah-langkah perhitungan berikut ini:


22 1. Tentukan transpos dari matrik kernel, yaitu GT 

    G=   



1

x1

y1

1

x2

1 .. .

x3 .. .

 y2   y3   ..  .   yN

1 xN



1

1

1

···

1



  Gt =  x1 x2 x3 · · · xN  y1 y2 y3 · · · yN

⇒

2. Tentukan GT G



1

1

1

···



       

1

 GT G =  x1 x2 x3 · · · xN y1

y2

y 3 · · · yN



1

x1

y1

1

x2

1 .. .

x3 .. .

   P P y2  N xi yi  P 2 P  P  y3  = xi xi xi yi   P P P ..  yi2 xi yi yi .  yN

1 xN

dimana N = jumlah data. dan i = 1, 2, 3, ..., N . 3. Kemudian tentukan pula GT d



1

1

1

···

1

 Gt d =  x1 x2 x3 · · · xN y1

y2

y3

· · · yN



       

d1



  P  d2  d i    P d3  = xi di   P ..  yi d i .  dN

4. Sekarang, persamaan (2.16) dapat dinyatakan sebagai 

   P  P P N xi yi m1 di P 2 P  P    P  xi xi xi di  xi yi   m2  =   P P P P 2 m3 yi d i yi xi yi yi

(2.17)

5. Sampai disini, jika tersedia data observasi, maka anda tinggal memasukan data tersebut ke dalam persamaan di atas, sehingga nilai elemen-elemen m dapat dihitung dengan perintah matlab m=inv(G’*G)*G’*d Langkah-langkah selanjutnya akan sama persis dengan catatan sebelumnya (model linear dan model parabola) Gabungan antara data observasi dan kurva hasil inversi diperlihatkan oleh Gambar 2.6

2.4. INVERSI MODEL BIDANG

23

Sebaran nilai terhadap X dan Y

40

Nilai

30

20

10

0 10 10 8

5

6 4

Y (m)

0

2 0

X (m)

Gambar 2.6: Hasil inversi berbentuk bidang. Saya sebut saja sebagai bidang inversi

Gambar 2.7: Ini adalah hasil inversi yang sama, hanya saja bidang inversi diputar (di-rotasi) pada suatu arah sehingga kita bisa saksikan dan pastikan bahwa sebaran titik data observasi berada di-dekat permukaan bidang inversi Script Matlab yang lengkap untuk tujuan ini adalah 1

clc


24 2 3

clear close

4 5 6 7 8

% Data observasi x = [2 5 7 4 1 3 6 9 8 4]; y = [3 6 2 7 8 9 4 1 5 5]; nilai = [10.6 23.5 27.3 20.8 11.1 18.9 25.4 33.5 33.2 24.1];

9 10 11 12 13 14 15 16

% Plot data observasi plot3(x,y,nilai,’*r’); grid; xlabel(’X (m)’); ylabel(’Y (m)’); zlabel(’Nilai’); title(’\fontsize{14} Sebaran nilai terhadap X dan Y’);ndata = length(nilai);

17 18 19 20 21 22 23 24 25

% Membentuk matrik kernel G dan vektor d ndata = length(x); for k = 1:ndata G(k,1) = 1; G(k,2) = x(k); G(k,3) = y(k); end d = nilai’;

26 27 28


29 30 31 32 33 34

% Plot hasil inversi (berupa garis least-squares) hold on; [X,Y] = meshgrid(min(x):max(x),min(y):max(y)); Z = m(1) + X.*m(2) + Y.*m(3); surf(X,Y,Z);

Anda bisa mengaplikasikan data pengukuran yang anda miliki, dengan syarat kasus yang anda tangani memiliki bentuk model yang sama dengan yang telah dikerjakan pada catatan ini, yaitu memiliki tiga buah model parameter yang tidak diketahui dalam bentuk persamaan bidang (atau 2-dimensi): d = m1 + m2 x + m3 y.

2.5 Contoh aplikasi 2.5.1

Menghitung gravitasi di planet X

Seorang astronot tiba di suatu planet yang tidak dikenal. Setibanya disana, ia segera mengeluarkan kamera otomatis, lalu melakukan ekperimen kinematika yaitu dengan melempar batu vertikal ke atas. Hasil foto-foto yang terekam dalam kamera otomatis adalah sebagai berikut Plot data pengukuran waktu vs ketinggian diperlihatkan pada Gambar 2.8. Anda diminta untuk membantu proses pengolahan data sehingga diperoleh nilai konstanta gravitasi di planet tersebut dan kecepatan awal batu. Jelas, ini adalah persoalan inversi, yaitu mencari unkown parameter (konstanta gravitasi dan kecepatan awal batu) dari data observasi (hasil foto gerak sebuah batu). Langkah awal untuk memecahkan persoalan ini adalah dengan mengajukan asumsi model

2.5. CONTOH APLIKASI

25

Tabel 2.4: Data ketinggian terhadap waktu dari planet X Waktu (dt) Ketinggian (m) Waktu (dt) Ketinggian (m) 0,00 5,00 2,75 7,62 0,25 5,75 3,00 7,25 0,50 6,40 3,25 6,77 0,75 6,94 3,50 6,20 1,00 7,38 3,75 5,52 1,25 7,72 4,00 4,73 1,50 7,96 4,25 3,85 1,75 8,10 4,50 2,86 2,00 8,13 4,75 1,77 2,25 8,07 5,00 0,58 2,50 7,90

matematika, yang digali dari konsep-konsep fisika, yang kira-kira paling cocok dengan situasi pengambilan data observasi. Salah satu konsep dari fisika yang bisa diajukan adalah konsep tentang Gerak-Lurus-Berubah-Beraturan (GLBB), yang formulasinya seperti ini 1 ho + vo t − gt2 = h 2 Berdasarkan tabel data observasi, ketinggian pada saat t = 0 adalah 5 m. Itu artinya ho = 5 m. Sehingga model matematika (formulasi GLBB) dapat dimodifikasi sedikit menjadi 1 vo t − gt2 = h − ho 2

(2.18)

Selanjut, didefinisikan m1 dan m2 sebagai berikut m1 = vo

1 m2 = − g 2

(2.19)

sehingga persamaan model GLBB menjadi m1 ti + m2 t2i = hi − 5

(2.20)

dimana i menunjukkan data ke-i. Langkah berikutnya adalah menentukan nilai tiap-tiap elemen matrik kernel, yaitu dengan memasukan data observasi kedalam model matematika (persamaan (2.20)) m1 t1 + m2 t21 = h1 − 5 m1 t2 + m2 t22 = h2 − 5 m1 t3 + m2 t23 = h3 − 5 .. .. . . . = .. m1 t20 + m2 t220 = h20 − 5


26 9

8

7

Tinggi (meter)

6

5

4

3

2

1

0

0

0.5

1

1.5

2

2.5 Waktu (detik)

3

3.5

4

4.5

5

Gambar 2.8: Grafik data pengukuran gerak batu

Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini:              

t1 t2 t3 t4 .. . t19 t20

t21



 h1 − 5      h2 − 5 " #    h3 − 5  m1   =  ..  m  2 . ..   .     h19 − 5 t219   h20 − 5 t220 t22 t23 t24

          

Operasi matrik di atas memenuhi persamaan matrik Gm = d Seperti yang sudah dipelajari pada bab ini, penyelesaian masalah inversi dimulai dari proses manipulasi persamaan matrik sehingga perkalian antara Gt dan G menghasilkan matriks bujursangkar Gt Gm = Gt d Selanjutnya, untuk mendapatkan m1 dan m2 , prosedur inversi dilakukan satu-per-satu

(2.21)


27

1. Menentukan transpos matrik kernel, yaitu Gt 

      G=      



t1

t21

t2

 t22   t23    2 t4   ..  .   t219   2 t20

t3 t4 .. . t19 t20

⇒

Gt =

"

t1 t2 t3 t4 . . . t19 t20 t21 t22 t23 t24 . . . t219 t220

#

2. Menentukan Gt G 

t

GG=

"

   #         

t1 t2 t3 t4 . . . t19 t20 t21 t22 t23 t24 . . . t219 t220



t1

t21

t2

 t22   t23   " P 2 P 3 #  ti ti 2 t4   = P t3 P t4 ..  i i .   t219   t220

t3 t4 .. . t19 t20

dimana N = 20 dan i = 1, 2, ..., N . 3. Kemudian menentukan hasil perkalian Gt d 

Gt d =

"

t1 t2 t3 t4 . . . t19 t20 t21 t22 t23 t24 . . . t219 t220

   #         

h1



 h2   # h3   " P  t h i i h4   = P t2 h ..  i i .   h19   h20

4. Sekarang persamaan (2.21) dapat dinyatakan sebagai " P

P

t2i t3i

P

P

t3i t4i

#"

m1 m2

#

=

" P

#

"

P

Berdasarkan data observasi, diperoleh "

179, 4

689, 1

689, 1 2822, 9

#"

m1 m2

=

ti hi t2i hi

273, 7 796, 3

#

(2.22)

#

Hasil operasi matriks ini dapat diselesaikan dengan satu baris statemen di matlab yaitu


28 9

8

7

Ketinggian (m)

6

5

4

3

2

1

0

0

0.5

1

1.5

2

2.5 Waktu (dt)

3

3.5

4

4.5

5

Gambar 2.9: Grafik hasil inversi parabola m=inv(G’*G)*G’*d Hasil inversinya adalah nilai kecepatan awal yaitu saat batu dilempar ke atas adalah sebesar m1 = vo = 3,2009 m/dt. Adapun percepatan gravitasi diperoleh dari m2 dimana m2 = − 21 g = -0,8169; maka disimpulkan nilai g adalah sebesar 1,6338 m/dt2 . Gambar 2.9 memperlihatkan kurva hasil inversi berserta sebaran titik data observasi. Garis berwarna biru merupakan garis kurva fitting hasil inversi parabola. Sedangkan bulatan berwarna merah adalah data pengukuran ketinggian (m) terhadap waktu (dt). Jelas terlihat bahwa garis kurva berwarna biru benar-benar cocok melewati semua titik data pengukuran. Ini menunjukkan tingkat akurasi yang sangat tinggi. Sehingga nilai kecepatan awal dan gravitasi hasil inversi cukup valid untuk menjelaskan gerak batu di planet X. 2.5.2

Analisa data seismik pada reflektor tunggal horizontal

Suatu survei seismik dilakukan untuk mengetahui kedalaman sebuah reflektor mendatar sebagaimana tampak pada gambar Waktu tempuh gelombang (t), yang bergerak sesuai dengan lintasan warna merah, memenuhi model matematika berikut ini

4z 2 x2 + 2 = t2 v2 v

(2.23)

Data observasi yang berhasil dihimpun dari survei tersebut adalah Berdasarkan data tersebut, tuliskan langkah demi langkah proses inversi selengkap mungkin untuk menentukan: • kecepatan gelombang seismik (v) pada lapisan • kedalaman reflektor mendatar (z) terhadap permukaan (surface)


S

29

R1

R2 R3 R4 R5 R6 R7 R8

Surface

z v Reflektor

Gambar 2.10: Reflektor mendatar pada kedalaman z. Kecepatan gelombang v dianggap konstan. S adalah sumber gelombang seismik dan R adalah penerima gelombang seismik. Jarak antara S dan R disebut offset (x). Sementara garis merah yang ada panahnya adalah lintasan gelombang seismik. Tabel 2.5: Data variasi offset (x) dan travel time (t) Receiver R ke-i Offset (x), meter Travel time (t), detik 1 60 0, 5147 2 80 0, 5151 3 100 0, 5155 4 120 0, 5161 5 140 0, 5167 6 160 0, 5175 7 180 0, 5183 8 200 0, 5192 Petunjuk: sederhanakan model matematika di atas menjadi m1 + m2 x2 = t2

(2.24)

Hasil yang benar adalah: kecepatan = 2797 m/dt ; kedalaman = 719 meter. Jika anda tidak mendapatkan kedua angka tersebut, maka proses inversi anda mungkin masih keliru! 2.5.3

Analisa data seismik pada reflektor tunggal miring

Suatu survei seismik dilakukan untuk mengetahui keberadaan sebuah reflektor miring sebagaimana tampak pada gambar Waktu tempuh gelombang (t), yang bergerak sesuai dengan lintasan warna merah, memenuhi model matematika berikut ini t2 =

4z 2 4xz sin α x2 + + 2 v2 v2 v

(2.25)

Data observasi yang berhasil dihimpun dari survei tersebut adalah Berdasarkan data tersebut, tentukan: • kecepatan gelombang seismik (v) pada lapisan • kedalaman reflektor miring (z) terhadap permukaan (surface) — jarak terdekat ke sumber gelombang seismik • sudut kemiringan reflektor (α)


30

x S

Surface

R2

R1

R3

R4

R5

z v a

Reflektor

O Gambar 2.11: Reflektor miring dengan sudut kemiringan sebesar α. Kecepatan gelombang v dianggap konstan. S adalah sumber gelombang seismik dan R adalah penerima gelombang seismik. Jarak antara S dan R disebut offset (x). Sementara garis merah yang ada panahnya adalah lintasan gelombang seismik. Tabel 2.6: Data variasi offset (x) dan travel time (t) Receiver R ke-i Offset (x), meter Travel time (t), detik 1 60 0, 4877 2 80 0, 4900 3 100 0, 4924 4 120 0, 4949 5 140 0, 4974 Petunjuk: sederhanakan model matematika di atas menjadi m1 + m2 x + m3 x2 = t2

(2.26)

2.6 Kesimpulan Dari sejumlah contoh pada bab ini, terlihat bahwa matrik kernel kerap kali berubah-ubah, sesuai dengan model matematika. Jadi, model matematika secara otomatis akan mempengaruhi bentuk rupa matrik kernelnya.

Bab 3

Penyelesaian Masalah Overdetermined

3.1 Regresi linear sederhana Jika suatu masalah inversi dapat direpresentasikan kedalam persamaan d = Gm, maka ia disebut linear. Kita dapat menjalankan prosedur yang sederhana untuk memperoleh nilai m dari data observasi. Dalam kenyataannya, tidak semua data observasi berhimpit dengan satu garis lurus. Jika kita mencoba melakukan fitting terhadap semua titik data observasi kepada satu garis, maka garis yang didapat disebut garis regresi. Misalnya, ada satu set data observasi yang ditulis sebagai (x1 , y1 ),(x2 , y2 ),...,(xn , yn ), garis regresi dinyatakan sebagai y = a0 + a1 x

(3.1)

yi = a0 + a1 xi + ei

(3.2)

dan setiap data memenuhi relasi berikut

dimana ei disebut error, residual, atau sering juga disebut misfit atau kesalahan prediksi (pre-

diction error ). Garis regresi tidak akan berhimpit dengan setiap data observasi dan biasanya untuk kasus inversi seperti ini selalu overdetermined. Secara umum, tipe masalah inversi seperti ini diselesaikan dengan metode least squares. Dengan metode least squares, kita mencoba meminimalkan error, ei , dengan cara menentukan nilai a0 dan a1 sedemikian rupa sehingga diperoleh jumlah-kuadrat-error, (S), yang minimal.

3.2 Metode least square Diketahui data eksperimen tersaji pada Tabel 3.1 Lalu data tersebut di-plot dalam sumbu x dan y. Sekilas, kita bisa melihat bahwa data yang telah di-plot tersebut dapat didekati dengan sebuah persamaan garis, yaitu a1 xi + a0 . Artinya, kita melakukan pendekatan secara linear, dimana fungsi pendekatan-nya adalah P (xi ) = a1 xi + a0 31

(3.3)

BAB 3. PENYELESAIAN MASALAH OVERDETERMINED

32

Tabel 3.1: Contoh data observasi yang dapat diolah oleh least squares xi yi xi yi 1 1,3 6 8,8 2 3,5 7 10,1 3 4,2 8 12,5 4 5,0 9 13,0 5 7,0 10 15,6 16 14 12 10

Y 8 6 4 2 0 1

2

3

4

5

6

7

8

9

10

X

Gambar 3.1: Hasil plotting data observasi dalam sumbu-x dan sumbu-y

Problemnya adalah berapakah nilai konstanta a1 dan a0 yang sedemikian rupa, sehingga posisi garis tersebut paling mendekati atau bahkan melalui titik-titik data yang telah di-plot di atas? Dengan kata lain, sebisa mungkin yi (pada persamaan (3.2)) sama dengan P (xi ) (pada persamaan (3.3)) atau dapat diformulasikan sebagai m X

yi − P (xi ) = 0

(3.4)

yi − (a1 xi + a0 ) = 0

(3.5)

i=1

m X i=1

dimana jumlah data, m = 10. Suku yang berada disebelah kiri dinamakan fungsi error (error function), yaitu E(a0 , a1 ) =

m X

yi − (a1 xi + a0 )

(3.6)

i=1

Semua data yang diperoleh melalui eksperimen, fungsi error-nya tidak pernah bernilai nol. Jadi, tidak pernah didapatkan garis yang berhimpit dengan semua titik data ekperimen. Namun demikian, kita masih bisa berharap agar fungsi error menghasilkan suatu nilai, dimana nilai tersebut adalah nilai yang paling minimum atau paling mendekati nol. Harapan tersebut diwujudkan oleh metode least square dengan sedikit modifikasi pada fungsi error-nya sehingga

3.2. METODE LEAST SQUARE

33

menjadi E(a0 , a1 ) =

m X

[yi − (a1 xi + a0 )]2

(3.7)

i=1

Agar fungsi error bisa mencapai nilai minimum, maka syarat yang harus dipenuhi adalah: ∂E(a0 , a1 ) =0 ∂ai

(3.8)

dimana i = 0 dan 1, karena dalam kasus ini memang cuma ada a0 dan a1 . Maka mesti ada dua buah turunan yaitu: m ∂ X ∂E(a0 , a1 ) [yi − (a1 xi + a0 )]2 = 0 = ∂a0 ∂a0 i=1

m X (yi − a1 xi − a0 )(−1) = 0 2 i=1

a0 .m + a1

m X

xi =

m X

yi

(3.9)

xi yi

(3.10)

i=1

i=1

dan m ∂E(a0 , a1 ) ∂ X = [yi − (a1 xi + a0 )]2 = 0 ∂a1 ∂a1 i=1

2

m X

(yi − a1 xi − a0 )(−xi ) = 0

i=1

a0

m X

xi + a1

i=1

m X

x2i

=

i=1

m X i=1

Akhirnya persamaan (3.9) dan (3.10) dapat dicari solusinya berikut ini: a0 =

Pm

dan a1 =

2 i=1 xi

m

m

Pm

Pm Pm i=1 yi − i=1 xi yi i=1 xi Pm 2 Pm 2 i=1 xi − ( i=1 xi )

Pm

Pm yi i=1 xi Pm i=12 − ( i=1 xi )

Pm

i=1 xi yi − Pm 2 m i=1 xi

(3.11)

(3.12)

Berdasarkan data ekperimen yang ditampilkan pada tabel diawal catatan ini, maka didapat: a0 =

385(81) − 55(572, 4) = −0, 360 10(385) − (55)2

dan a1 =

10(572, 4) − 55(81) = 1, 538 10(385) − (55)2


34

Jadi, fungsi pendekatan-nya, P (xi ), adalah P (xi ) = 1, 538xi − 0, 360 Nilai a0 dan a1 disebut koefisien regresi. Lebih jauh lagi a0 disebut intercept (titik perpotongan) terhadap sumbu y sedangkan a1 adalah gradient atau slope (kemiringan garis). Gambar di bawah ini menampilkan solusi regresi linear tersebut berikut semua titik datanya 16

P(x) = 1.538*x − 0.36

14 12 10 8 6 4 2 0 −2 0

2

4

6

8

10

Gambar 3.2: Contoh solusi regresi linear Teknik diatas diterapkan secara rutin dalam analisis data geofisika, khususnya ketika kita mencoba meng-esktrak satu atau dua parameter model dari data observasi. Teknik ini disebut analisis regresi linear (linear regression analysis ) atau classical least squares fitting. Teknik ini pertama kali dipakai oleh Gauss pada tahun 1809. Teknik ini pada mulanya digunakan untuk mencari solusi dari masalah overdetermined namun pada perkembangannya teknik ini diterapkan juga pada underdetermined problem setelah dimodifikasi. Ketika kita ingin mendapatkan lebih dari 2 parameter maka teknik ini disebut analisis regresi multipel (multiple regression

analysis ).

3.3 Aplikasi regresi linear pada analisis data seismik refraksi Tabel 3.2 menampilkan data survei seismik refraksi yang dilakukan dengan jarak offset xi dan waktu tempuh gelombang dari source ke receiver dicatat dalam ti . Anggap saja persamaan travel time adalah ti =

xi + Th v

(3.13)

Parameter xi dan ti berperan sebagai known parameter (parameter yang diketahui). Sementara kecepatan gelombang pada medium, v, dan waktu tempuh gelombang secara vertikal, Th bertindak sebagai unknowm parameter (parameter yang tidak diketahui). Metode regresi linear, atau yang juga akrab disebut least square, dilakukan untuk memecahkan unknown parameter tersebut, artinya metode regresi linear berupaya mencari nilai v dan Th . Pertama-tama proses

3.3. APLIKASI REGRESI LINEAR PADA ANALISIS DATA SEISMIK REFRAKSI

35

Tabel 3.2: Data seismik refraksi: waktu-datang gelombang, ti , dan jarak antara source dan receiver atau jarak offset, xi xi (m) 2 4 6 8

Trace 1 2 3 4

ti (ms) 5,1 9,2 11,9 14,9

linearisasi dilakukan terhadap persamaan (3.13) sehingga menjadi ti = a0 + a1 xi

a0 = Th

dimana

dan a1 =

1 v

(3.14)

Kesalahan (error ) diasumsikan hanya berasal dari cuplikan waktu gelombang datang. Penerapan metode regresi linear yang berusaha meminimalkan jumlah kuadrat dari error, ei = ti − (a0 + a1 xi ), dapat dinyatakan sesuai persamaan (3.11) dan (3.12). Dengan jumlah data, m = 4, maka a0 =

P4

dan a1 =

P P − 4i=1 xi ti 4i=1 xi = 2, 25 2 P P 4 4 2 4 i=1 xi i=1 xi −

2 i=1 xi

4

P4

P4

i=1 ti

i=1 xi ti

4

P

−

4 2 i=1 xi

P4

i=1 xi

−

P4

i=1 ti 2 i=1 xi

P 4

= 1, 605

Dengan demikian kecepatan gelombang seismik pada medium tersebut1 adalah v = 1/a1 = 623, 053 m/dt. Adapun standard error χ2a1 dan χ2a0 ditentukan oleh rumus berikut χ2 D P

χ2a1

= m

χ2a0

= χ2

(3.15)

x2 D

(3.16)

dimana D = m

m X

x2i

i=1

χ2 =

m

!

−

1 X 2 Ei m−2

m X i=1

xi

!2

i=1

Sebagai catatan tambahan, χ2 adalah nilai deviasi rms (root mean square ) dari data ti terhadap garis regresi hasil analisis (a0 + a1 xi ) dengan faktor (m − 2) karena dalam masalah ini hanya dicari 2 parameter model (a0 dan a1 ). 1

Nilai kecepatan yang diperoleh adalah nilai hasil analisis regresi linear. Namun itu bukan nilai yang mutlak, karena bisa jadi, pengukuran kecepatan gelombang menggunakan (misalnya) sonic-log memperoleh hasil yang berbeda, walaupun perbedaan tersebut tidak boleh terlalu jauh atau masih bisa diterima


36

3.4 Regresi linear dengan pendekatan matriks Metode regresi linear dapat didekati dengan operasi matriks. Caranya sama persis dengan yang telah dibahas pada Bab 2. Pendekatan matrik ini dimulai dari upaya membuat persamaan matrik d = Gm dari persamaan model (persamaan (3.13)). Unknown parameter, yang hendak dipecahkan, terkandung pada elemen-elemen vektor m. Jika data yang kita miliki sangat ideal dalam arti tidak ada error sama sekali, maka m bisa diperoleh sebagai berikut m = G−1 d Akan tetapi, pada kenyataannya semua data pengukuran pasti memiliki error yang besarnya relatif bervariasi. Karenanya, data observasi tak akan pernah fit secara sempurna dengan model. Itu artinya keberadaan misfit tidak pernah bisa dihindari. Konsekuensinya, misfit (baca: error) tersebut mesti disertakan pada d = Gm d = Gm + ei

(3.17)

dan selanjutnya, solusi regresi linear diupayakan dengan cara meminimalkan jumlah kuadrat dari error, ei . Cara ini tak lain berupaya untuk memperoleh misfit terkecil yaitu jarak perbedaan terkecil antara data survei dan model. Dalam formulasi matematika, kuadrat error tersebut dinyatakan dengan q = eT e = (d − Gm)T (d − Gm)

(3.18)

Dimana simbol T maksudnya adalah operasi transpose. Agar kuadrat error di atas menghasilkan nilai minimal, maka persamaan (3.18), diturunkan terhadap m dan hasilnya harus sama dengan nol

∂ dT d − dT Gm − mT GT d + mT GT Gm ∂q = =0 ∂mj ∂mj

sehingga −dT G − GT d + GT Gm + mT GT G = 0 akhirnya diperoleh 2GT Gm = 2GT d GT Gm = GT d

(3.19)

Persamaan (3.19) disebut persamaan normal2 . Dengan persamaan normal, estimasi unknown parameter yang terkandung pada vektor m ditentukan oleh −1 T m = GT G G d

(3.20)

Persamaan (3.20) disebut unconstrained least squares terhadap masalah inversi d = Gm. −1 T Bagian GT G G dinamakan Generalized Inverse yang mengolah data d untuk memperoleh 2

Sebetulnya kita sudah menuliskannya di Bab 2 sebagai persamaan (2.4), namun tanpa penurunan rumus

3.4. REGRESI LINEAR DENGAN PENDEKATAN MATRIKS

37

parameter model m. Untuk menyelesaikan persamaan (3.20) dengan operasi matriks secara numerik atau komputasi bisa menggunakan beberapa metode, diantaranya metode Eliminasi Gauss, LU-Decomposition, Iterasi Gauss-Seidel, dan Singular Value Decomposition. 3.4.0.1

Aplikasi pada data survei seismik refraksi

Kembali kita tampilkan Tabel 3.3 merupakan data survei seismik refraksi yang dilakukan dengan jarak offset xi , dan persamaan travel time adalah ti =

xi + Th v

Tabel 3.3: Data seismik refraksi: waktu-datang gelombang (ti ) pada empat posisi geophone (xi ) xi (m) 2 4 6 8

Trace 1 2 3 4

ti (ms) 5,1 9,2 11,9 14,9

Pertama-tama proses linearisasi dilakukan terhadap persamaan travel time sehingga menjadi ti = a0 + a1 xi

dimana

a0 = Th

dan a1 =

1 v

Berdasarkan linearisasi tersebut, kita bisa menyusun sistem persamaan linear sebagai berikut: t1 = a0 + a1 x1 t2 = a0 + a1 x2 t3 = a0 + a1 x3 t4 = a0 + a1 x4 Semua persamaan tersebut dapat dinyatakan dalam operasi matrik (d = Gm) berikut ini: 

t1





1 x1



#    "  t2   1 x2  a0      t = 1 x  a 3  1  3   t4 1 x4 dimana 

t1



   t2    d=   t3  t4



1 x1

  1 x2 G=  1 x 3  1 x4

     

dan

m=

"

a0 a1

#


38

Selanjutnya, unknown parameter, a0 dan a1 , dipecahkan dengan menggunakan persamaan (3.20). Untuk menuju kesana, langkah-langkah penyelesaian berikutnya adalah 1. Menentukan transpos dari matrik G, yaitu GT 

1 x1



   1 x2    G=  1 x 3   1 x4

⇒

GT =

"

1

1

1

1

x1 x2 x3 x4

#

2. Melakukan perkalian matriks GT G

T

G G=

"

1

1



1 x1

#  1 x2  x4   1 x3 1 x4

1

1

x1 x2 x3



# P  "  N x i = P P 2  xi xi 

dimana N = 4 dan i = 1, 2, 3, 4. 3. Kemudian menentukan perkalian GT d

GT d =

"

1

1

1

x1 x2 x3



t1



# #  " P  t2  t  = P i  x4  t xi ti  3  t4 1

4. Sekarang persamaan (3.20) dapat dinyatakan sebagai "

a0 a1

#

=

"

#−1 " P # P N xi ti P P 2 P xi xi xi ti

(3.21)

Berdasarkan data survei pada Tabel 3.3 di atas, diperoleh "

a0 a1

#

=

"

4

20

20 120

#−1 "

41, 1 237, 6

#

(3.22)

Operasi matriks di atas akan menghasilkan nilai a0 = 2, 25 dan a1 = 1, 605. Dengan demikian v = 1/a1 = 623, 053 m/dt

3.5. SOAL

39

3.5 Soal Diketahui data temperatur borehole sebagai berikut. Tentukan slope dan intercept pada z-axis dan kemudian perkirakan temperatur pada kedalaman 390m. Depth, z(m) 30 70 180 250 300

Temp, t(o C) 25,0 26,2 29,7 34,3 35,5

Bab 4

Constrained Linear Least Squares Inversion

Geofisika merupakan ilmu yang bertujuan untuk memahami kondisi bawah permukaan bumi dengan menggunakan prinsip-prinsip ilmu fisika. Karena letaknya yang berada dibawah permukaan bumi, maka obyek geofisika tidak bisa diamati dengan mata telanjang. Sebagai gantinya, mau tidak mau obyek tersebut dipelajari dengan mengamati respon bawah permukaan bumi ketika suatu gangguan fisis diberikan padanya. Selanjutnya, respon tersebut dianalisis hingga diperoleh solusi yang nantinya siap menjadi acuan interpretasi. Masalahnya, pada kebanyakan pengamatan geofisika, sangat mungkin mendapatkan solusi yang berbeda-beda, yang semuanya itu bisa saja dipakai untuk menjelaskan data observasi. Inilah yang disebut dengan istilah non-uniqueness problem atau non-unik. Namun demikian pada akhirnya, kita harus memilih satu buah solusi yang paling baik menurut kita. Karena kita sadar bahwa obyek yang sedang kita incar dibawah sana, pastilah hanya berada dalam satu kondisi tertentu yang unik. Untuk mendapatkan solusi yang unik, kita harus menambahkan sejumlah informasi yang sebelumnya tidak ada pada persamaan least square d = Gm. Informasi tambahan ini disebut a priori information, yang selanjutnya akan digunakan untuk meng-constrain solusi sehingga diperoleh solusi yang dianggap paling tepat untuk menggambarkan kondisi bawah permukaan. A priori information, atau yang saya indonesiakan menjadi informasi awal, ini didapat dari data geofisika yang lainnya, atau bisa juga dari data borehole, atau juga bersumber dari data geologi.

4.1 Inversi dengan informasi awal Kita dapat menambahkan informasi awal (h) kepada parameter model dalam suatu proses inversi. Secara umum, informasi awal (h) tersebut diharapkan membantu pemodelan inversi sehingga diperoleh hasil yang unik dari sejumlah kemungkinan solusi. Sekali lagi, proses ini disebut meng-constrain, sementara solusi yang diperoleh disebut solusi ter-constrain. Constrain 41

BAB 4. CONSTRAINED LINEAR LEAST SQUARES INVERSION

42

terhadap suatu data dirumuskan sebagai berikut Dm = h

(4.1)

Dimana D adalah matrik -dengan seluruh elemen selain diagonal bernilai nol- yang beroperasi pada parameter model m sedemikian rupa sehingga parameter m "dipaksa" untuk sama persis dengan informasi awal (h). Kalau kita menyelesaikan persamaan (4.1), berarti kita telah melakukan apa yang disebut dengan linear equality constraints. Formulasi matematika-nya adalah sebagai berikut φ = (d − Gm)T (d − Gm) + β 2 (Dm − h)T (Dm − h)

(4.2)

Untuk mendapatkan error minimum maka turunan φ terhadap parameter model m dibuat sama dengan nol 2GT Gm − 2GT d + 2β 2 DT Dm − 2β 2 DT h = 0 sehingga diperoleh persamaan normal (GT G + β 2 DT D)m = GT d + β 2 DT h

(4.3)

Ketika D adalah matrik identitas, maka persamaan normal berubah menjadi (GT G + β 2 I)m = (GT d + β 2 h)

(4.4)

Dari sini solusi ter-constrain didapat sebagai berikut mˆc = (GT G + β 2 I)−1 (GT d + β 2 h)

(4.5)

Formula ini dinamakan inversi linear terkonstrain atau disebut juga the biased linear estimation technique. Keuntungannya adalah formula ini dapat membantu mendapatkan satu solusi yang unik dari sejumlah solusi yang mungkin pada masalah overdetermined, dimana didalamnya terdapat ketidakpastian sebagai akibat dari kesalahan pengukuran (observational errors) dan ketidakpastian (uncertainties). Parameter β ditentukan secara coba-coba (trial and error), namun biasanya ia memiliki nilai antara nol dan satu. Konstanta β disebut faktor pengali undetermined atau disebut juga faktor pengali Lagrange. Sehingga metode ini disebut Lagrange multiplier method (metode pengali Lagrange).

4.1. INVERSI DENGAN INFORMASI AWAL 4.1.1

43

Memformulasikan persamaan terkonstrain

Persamaan Dm = h secara umum memiliki bentuk 

1

   .  .  .

1 .. .





m1 m2 .. .

   ..   .  1



    =    

mp

h1 h2 .. . hp

     

(4.6)

Namun demikian, persamaan matrik di atas dapat dimodifikasi sesuai kebutuhan. Misalnya jika informasi awal yang diketahui hanya ada satu, modifikasinya menjadi

h



i  1 0 ... 0   

m1 m2 .. . mp



 h i   = hknown  

(4.7)

Jika pada kasus lain, kita punya informasi awal yaitu parameter pertama dan parameter ke empat, maka persamaan matrik terkonstrain menjadi       4.1.2

1 0



m1





h1



      m2   0   =   m   0  0  3    1 m4 h4

(4.8)

Contoh aplikasi inversi terkonstrain

Contoh 1: Least squares garis terkonstrain Sekarang kita ingin menerapkan inversi terkonstrain pada pengolahan data first arrivals dari seismik refraksi. Persamaan least square garis adalah ti = m1 + m2 xi

(4.9)

dengan parameter model yang terdiri atas m1 dan m2 . Sementara data lapangan merupakan pasangan dari jarak offset xi dan waktu datang gelombang (first arrival time) ti . Jika ada 4 data lapangan, sistem persamaan linear yang bisa disusun adalah t1 = m1 + m2 x1 t2 = m1 + m2 x2 t3 = m1 + m2 x3 t4 = m1 + m2 x4 Sekarang kita berasumsi memiliki informasi konstrain dari kegiatan explorasi sebelumnya bahwa solusi least square harus melewati titik koordinat (xc , tc ). Jadi kita harus meng-konstrain


44

solusi least square untuk mengakomodasi informasi awal tersebut. h

1 xc

i

"

#

m1 m2

=

h

tc

i

(4.10)

Persamaan matrik di atas harus diintegrasikan dengan d = Gm sehingga solusi akhir merupakan solusi terkonstrain yang kita harapkan lebih akurat dibandingkan jika tidak terkonstrain.

Tentu anda masih ingat pada least square tidak terkonstrain, dimana dikenal dua komponen berikut GT G dan GT d. GT G =

"

# P n xi P 2 P xi xi

dan GT d =

" P P

ti

xi ti

#

Agar menjadi terkonstrain kedua komponen itu mesti dimodifikasi sedemikian rupa sehingga menjadi 

 P n xi 1 P 2  P  (GT G + β 2 I) =  xi xi xc  1 xc 0

(4.11)

dan  P

ti



  P (GT d + β 2 h) =  xi ti  tc

(4.12)

Sebelum dilanjut, mari kita tengok lagi persamaan normal (4.4) secara utuh, yaitu (GT G + β 2 I)m = (GT d + β 2 h) Sekarang mari kita masukkan kedua komponen kedalam persamaan (4.4)    P  P ti m1 n xi 1 P 2    P   P xi ti  xi xi xc   m2  =   tc β 1 xc 0 

(4.13)

dimana nilai β mesti kita tentukan secara coba-coba namun dalam batas antara nol dan satu.

4.1. INVERSI DENGAN INFORMASI AWAL

45

Tabel 4.1: Data seismik refraksi: waktu-datang gelombang (ti ) pada empat posisi geophone (xi ) Trace 1 2 3 4

xi (m) 2 4 6 8

ti (ms) 5,1 9,2 11,9 14,9

Akhirnya, solusi terkonstrain terhadap inversi garis yang melewati (xc , tc ) adalah  −1  P P n xi 1 ti P 2    P    P mˆc =  m2  =  xi ti  xi xi xc   β 1 xc 0 tc 

m1





(4.14)

Sekarang, kita melangkah pada kasus nyata. Tabel 4.1 menunjukkan data pengukuran seismik refraksi. Tentukan parameter model untuk persamaan garis yang melewati titik (xc = 8, yc = 14, 9)! Langkah pertama kita hitung komponen matrik 

   P n xi 1 4 20 1 P 2  P    (GT G + β 2 I) =  xi xi xc  =  20 120 8  1 xc 0 1 8 0

(4.15)

Kemudian kita hitung komponen matrik yang lain  P

ti





41, 1



 P    (GT d + β 2 h) =  xi ti  =  237, 6  tc 14, 9

(4.16)

Problem ini dapat diselesaikan secara numerik dengan metode Eliminasi Gauss. Solusi yang diperoleh adalah m1 = 2, 3857, m2 = 1, 5643 dan β = 0, 2714. Berikut ini adalah script lengkapnya dalam Matlab clear all clc; a(1,1)=4; a(1,2)=20; a(1,3)=1; a(1,4)=41.1; a(2,1)=20; a(2,2)=120; a(2,3)=8; a(2,4)=237.6; a(3,1)=1; a(3,2)=8;


46 a(3,3)=0; a(3,4)=14.9; a n=3; % n = jumlah baris

% berikut ini proses triangularisasi ----------for j=1:(n-1) kk=j+1; for k=kk:n m(k,j)=a(k,j)/a(j,j) for i=1:(n+1) a(k,i)=a(k,i)-m(k,j)*a(j,i); end end end % akhir dari proses triangularisasi ------------a x(n,1)=a(n,n+1)/a(n,n); % berikut ini proses substitusi mundur ---------for k=1:n-1 i=n-1-k+1; j=i+1; sum=0; for k=j:n sum = sum + a(i,k) * x(k,1); end x(i,1)=(a(i,n+1)-sum)/a(i,i); end % akhir dari proses substitusi mundur ------------x

4.2 Inversi dengan Smoothness Cara yang paling efektif untuk menginversi data yang terbatas adalah dengan menentukan konstrain sehingan solusi yang diinginkan menjadi smooth (halus). Tingkatan smooth dapat diukur berdasarkan kondisi fisis atau kondisi geologi.

4.2.1

Formulasi masalah

Mari kita pelajari bagaimana suatu masalah dapat diformulasikan menuju solusi yang smooth. Jika diinginkan parameter model bervariasi dengan jarak selisih yang kecil, maka lakukanlah proses minimalisasi perbedaan paramter yang berdekatan (m1 −m2 ), (m2 −m3 ),..., (mp−1 −mp ).

4.2. INVERSI DENGAN SMOOTHNESS

47

Perbedaan ini dinyatakan sebagai persamaan konstrain Dm = h 



1 −1

    

1

−1 1 −1

    

m1 m2 .. .





0



     0  =    0     0 mp

(4.17)

dimana D adalah operator selisih yang bertindak sebagai matrik smoothness dan Dm disebut penghalus (flatness) solusi vektor parameter model m. Jika parameter model tidak bervariasi secara smooth terhadap posisi, maka gunakanlah persamaan konstrain berbentuk 

1

   .  .  .

1 .. .



   ..   .  1

m1 m2 .. .





0



     0  =    0     0 mp

(4.18)

Dalam kasus ini, D adalah matrik identitas dengan dimensi p × p dan dimensi h adalah p × 1. Operasi ini akan mendorong proses inversi menuju kondisi stabil. Untuk mendapatkan solusi yang smooth, kita gunakan ukuran selisih seperti persamaan (4.2), dinyatakan sebagai q2 (m) = (Dm − h)T (Dm − h) = mT DT Dm = mT Hm

(4.19)

dimana H = DT D. Kita nyatakan mengenai masalah terkonstrain adalah: Dimulai dari data lapangan yang tidak komplit, tidak lengkap, tidak cukup, kita mencari seluruh kemungkinan solusi dengan residual q1 = |d − Gm|2 dan solusi yang paling smoothness dengan judgement dari ukuran q2 (m). Secara matematik, pernyataan di atas memiliki maksud: meminimalkan q2 = mT Hm dibawah kondisi |d − Gm|2 = q1 atau secara umum |d − Gm|2 ≤ qT , dimana qT adalah nilai toleransi maksimum dari residual atau misfit. Masalah konstrain membutuhkan minimalisasi kd − Gmk2 dan q2 (m) secara bersamaan, φ = (d − Gm)T (d − Gm) + β 2 (mT DT Dm) 4.2.2

(4.20)

Solusi masalah

Untuk mendapatkan solusi parameter model, perlu dilakukan minimalisasi terhadap persamaan (4.20),

∂ dT d − mT GT d − dT Gm + mT GT Gm + β 2 mT Hm =0 ∂mj

sehingga GT G + β 2 H m = GT d

48


Ini adalah persamaan normal yang baru. Dan akhirnya solusi smoothness diturunkan sebagai berikut ms = GT G + β 2 H Dan bila D = I, ms = GT G + β 2 I

−1

GT d

(4.21)

−1

GT d

(4.22)

Persamaan (4.22) lebih populer disebut Damped Least Squares solution atau solusi Least Square Teredam. Nama lainya yang juga cukup terkenal adalah inversi Marquardt.

Daftar Pustaka

[1] Meju, A Max., Geophysical Data Analysis: Understanding Inverse Problem Theory and Practice, (1994), Society of Exploration Geophysicists (SEG)

49

Indeks

akurasi, 3

model konseptual, 4

arus fluida, 3

model matematika, 1, 4

bahan tambang, 1 bessel, 7 data lapangan, 1 data observasi, 1

modulus bulk, 4 momen inersia, 5 noise, 3 non-linear, 6

densitas, 4

observasi, 31

eksplorasi, 3

parameterisasi, 6

elektroda, 7

permeabilitas, 4

fitting, 1

resistivity, 1

forward, 1 sampling rate, 3 gelombang elektromagnetik, 3

Schlumberger, 6

gelombang seismik, 3

sedimentasi, 3

gempa bumi, 3

seismik refraksi, 1

geologi, 4

slowness, 6

hambatan jenis, 4 hidrogeologi, 3

sumur bor, 3 unknown parameter, 1

highly non-linear, 7 human error, 4 instrumen, 4 instrumental error, 4 inversi, 1 inversi diskrit, 6 jejari bumi, 5 kernel, 5 komputasi, 5 laboratorium, 3 linearisasi, 6 lithospere, 3 matematika diskrit, 5 model fisik, 4 51

Analisis Data Geofisika: Memahami Teori Inversi

Recommend Documents