Penyelesaian Ves Menggunakan Metode Monte Carlo.docx

NAMA NIM PAPER

PENYELESAIAN VES MENGGUNAKAN METODE MONTE CARLO

Pada metode VES ada beberapa penyelesaikan atau pendekatan yang dapat diaplikasikan untuk menjadi model salah satunya dengan Monte Carlo. Metode ini terdapat pada pemodelan inversi dimana pemilihan model pada metode pencarian acak (random search) sesuai namanya dilakukan secara acak. Setiap Model dalam ruang model memiliki peluang yang sama untuk dipilih sebagai sample model. Bilangan acak dibangkitkan dengan probabilitas uniform anta 0 dan 1 yang kemudia dipetakan pada interval harga parameter model. Perhitungan pemodelan kedepan dilakukan untuk model terpilih yang jumlahnya tidak terlalu besar jika dibandingkan dengan jumlah keseluruhan model yang mungkin pada ruang model. (Grandis 2009). Model terdiri dari sejumlah lapisan dengan ketebalan tetap dengan resistivitas sebagai parameter model yang dicari. Pilihan harga resistivitas lapisan adalah harga-harga diskret yang dapat mewakili medium konduktif sampai resistif. Algortima MCMC didasarkan pada sampling ruang model secara lebih terarah dengan memanfaatkan probabilitas transisi rantai Markov. Harga resistivitas untuk satu lapisan tertentu memiliki probabilitas terpilih lebih besar jika menghasilkan misfit lebih kecil. Karena adanya ekivalensi, konfigurasi resistivitas lapisan yang bervariasi dapat menghasilkan misfit yang dekat dengan harga optimumnya. Oleh karena itu selain misfit kecil, model disebut optimum jika variasi resistivitasnya minimum (smooth). Penerapan inversi pada data sintetik maupun data lapangan memberikan hasil yang memuaskan. Teknik inversi ini akan

dikembangkan dengan menambahkan kriteria kesinambungan antara satu lapisan dengan lapisan yang sama dari titik sounding yang berdekatan untuk memudahkan korelasi dan interpretasi. Pada metode Monte Carlo jumlah model yang harus digunakan berkorelasi dengan jumlah model yang mungkin tersebut(NM). Harga resistivitas suatu lapisan ρi dipilih dari sejumlah M harga Rj; j = 1, 2, ..., M yang terdistribusi secara homogen dalam skala logaritmik antara Rmin dan Rmax. Probabilitas Rj untuk resistivitas suatu lapisan ρi dinyatakan oleh: P(R j ) = exp(-E(m | ri = R j ))

(1)

E(m|ri = Rj) adalah misfit dari suatu model m dimana ri = Rj, dengan resistivitas lapisan selain lapisan ke-i memiliki harga tertentu yang tetap. Teknik filter Gosh-Koefoed (Koefoed, 1979) digunakan untuk menghitung respons model 1-D yaitu resistivitas-semu sebagai fungsi dari spasi elektroda Schlumberger (AB/2). Resistivitas suatu lapisan dipilih secara acak dari Rj namun dengan probabilitas P(Rj); j = 1, 2, ..., M sebagai bobot. Harga resistivitas yang berasosiasi dengan misfit rendah memiliki probabilitas tinggi untuk terpilih, demikian pula sebaliknya. Prosedur pemilihan acak berbobot (weighted random selection) dilakukan melalui mekanisme roulette wheel sebagaimana diuraikan oleh Grandis (2009). Rantai Markov adalah rangkaian keadaan (state) yang dikontrol oleh probabilitas transisi yang secara umum pada tahap ke-n dinyatakan oleh: Pn (a, b)  P( X n 1  b | X n  a)

(2)

Persamaan (2) adalah probabilitas bersyarat (conditional probability) yang menyatakan probabilitas suatu keadaan b terjadi jika keadaan sebelumnya adalah a. Persamaan (1) pada dasarnya adalah probabilitas transisi suatu rantai Markov. Dalam hal ini state menyatakan

NAMA NIM

model, dimana model a dan model b masing-masing adalah model sebelum dan sesudah resistivitas suatu lapisan diperbaharui. Dimulai dengan harga resistivitas homogen, pembaharuan (update) harga resistivitas semua lapisan dilakukan secara iteratif. Model yang dihasilkan membentuk rantai Markov yang konvergen baik secara teoritik maupun empirik menuju model optimum (Grandis dkk., 1999; 2002). Fenomena ekivalensi serta jumlah lapisan yang cukup besar menyebabkan model hasil inversi yang berbeda memiliki misfit yang hampir sama. Untuk menstabilkan inversi maka digunakan kendala tambahan dengan meminimumkan variasi resistivitas antar lapisan. Hal ini menyebabkan model hasil inversi bersifat smooth sehingga identik dengan smoothness constrained inversion (Constable dkk., 1987).

NAMA NIM

PENURUNAN RUMUS INTEGRAL FUNGSI BESSEL

Fungsi Bessel merupakan solusi kanonik y(x) dari persamaan diferensial Bessel. umumnya, fungsi bessel juga dikenal sebagai fungsi silinder atau harmonika silindris (cylindrical harmonics), sebab mereka dijumpai dalam penyelesaian persamaan Laplace pada koordinat tabung/silindris. Dengan bentuk umum Persamaan Diferensial Bessel : 𝑥 2 𝑦′′ + 𝑥𝑦′ + (𝑥 2 − 𝑣 2 )𝑦 = 0 ∞

∞

𝑟

𝑦(𝑥) = 𝑥 ∑ 𝑎𝑚 𝑥

𝑚

= ∑ 𝑎𝑚 𝑥 𝑚+𝑟

𝑚=0

𝑚=0

dengan syarat a₀≠0, sehingga : ∞

∞

𝑦′(𝑥) = 𝑥 𝑟 ∑ (𝑚 + 𝑟)𝑎𝑚 𝑥 𝑚+𝑟−1 = 𝑥 𝑟−1 ∑ (𝑚 + 𝑟)𝑎𝑚 𝑥 𝑚 𝑚=0

∞

𝑚=0 ∞

𝑦′′(𝑥) = 𝑥 𝑟 ∑ (𝑚 + 𝑟)(𝑚 + 𝑟 − 1)𝑎𝑚 𝑥 𝑚+𝑟−2 = 𝑥 𝑟−2 ∑ (𝑚 + 𝑟)(𝑚 + 𝑟 − 1)𝑎𝑚 𝑥 𝑚 𝑚=0

𝑚=0

Persamaan Diferensial menjadi : ∞ 2

𝑋 [𝑥

𝑟−2

∞ 𝑚

∑ (𝑚 + 𝑟)(𝑚 + 𝑟 − 1)𝑎𝑚 𝑥 ] + 𝑥 [𝑥 𝑚=0

𝑟−1

∑ (𝑚 + 𝑟)𝑎𝑚 𝑥 𝑚 ] + (𝑥 2 𝑚=0

∞

− 𝑣 2 ) [𝑥 𝑟 ∑ 𝑎𝑚 𝑥 𝑚+𝑟 ] = 0 𝑚=0

Atau, ∞

∞

∑ (𝑚 + 𝑟)(𝑚 + 𝑟 − 1)𝑎𝑚 𝑥 𝑚=0

∞

𝑚+𝑟

+ ∑ (𝑚 + 𝑟)𝑎𝑚 𝑥 𝑚+𝑟 𝑚=0 ∞

+ ∑ 𝑎𝑚 𝑥 𝑚+𝑟+2 − 𝑣 2 ∑ 𝑎𝑚 𝑥 𝑚+𝑟 = 0 𝑚=0

𝑚=0

Jika x tidak selalu nol, maka yang pasti = 0 adalah koefisien-koefisien dari 𝑥 𝑟+𝑠 dengan m=s dalam jumlah pertama, kedua dan keempat, dan dengan m=s-2 pada jumlah yang ketiga. Sehingga dapat dijabarkan sebagai berikut : Koefisien 𝑥 𝑟 : (𝑟 − 1)𝑟𝑎0 + 𝑟𝑎0 − 𝑣 2 𝑎0 = 0 (𝑟 2 − 𝑟 + 𝑟 − 𝑣 2 )𝑎0 = 0 (𝑟 2 − 𝑣 2 )𝑎0 = 0; 𝑎0 ≠ 0 Persamaan indical : 𝑟 2 − 𝑣 2 = 0; 𝑟12 = ±𝑣 Koefisien 𝑥 𝑟+1 : 𝑟(𝑟 + 1)𝑎1 + (𝑟 + 1)𝑎1 − 𝑣 2 𝑎1 = 0 (𝑟 2 + 𝑟 + 𝑟 − 𝑣 2 )𝑎1 = 0 (2𝑟 + 1 + 𝑟 2 − 𝑣 2 )𝑎1 = 0

NAMA NIM

(2𝑟 + 1)𝑎1 = 0; (2𝑟 + 1)𝑡𝑖𝑑𝑎𝑘𝑠𝑒𝑙𝑎𝑙𝑢0 𝑎1 = 0 𝑟+1 (𝑠 Koefisien 𝑥 : + 𝑟 − 1)(𝑠 + 𝑟)𝑎𝑠 + (𝑠 + 𝑟)𝑎𝑠−2 − 𝑣 2 𝑎𝑠 = 0 [(𝑠 + 𝑟)(𝑠 + 𝑟 − 1 + 1) − 𝑣 2 ]𝑎𝑠 = −𝑎𝑠−2 [(𝑠 + 𝑟)2 − 𝑣 2 ]𝑎𝑠 = −𝑎𝑠−2 𝑎𝑠−2 𝑎𝑠 = (𝑠 + 𝑟)2 − 𝑣 2 Untuk r = v : 𝑎𝑠−2 𝑎𝑠−2 𝑎𝑠−2 𝑎𝑠 = = −𝑎𝑠 = 2 =− 2 2 2 2 (𝑠 + 𝑟) − 𝑣 𝑠 + 2𝑠𝑣 + 𝑣 − 𝑣 𝑠(𝑠 + 2𝑣) Karena a₁ = 0 ; v ≥ 0, maka untuk s ganjil 𝑎𝑠 = 0 dan s genap = 2m ; m = 1, 2, 3, ... 1 1 𝑎2𝑚 = − 𝑎2𝑚−2 = − 2 𝑎 2𝑚(2𝑣 + 2𝑚) 2 𝑚(𝑣 + 𝑚) 2𝑚−2 1

Karena a₀ sebarang dan a₀≠0, maka bisa dipilih 𝑎0 = 2𝑣 Γ(𝑣+1) Dimana Γ(𝑣 + 1) adalah fungsi gamma. Untuk keperluan di sini cukup ketahui bahwa Γ(𝛼) didefinisikan oleh integral ∞

Γ(𝛼) = ∫0 𝑒 −𝑡 𝑡 𝛼−1 𝑑𝑡

(α ≥ 0)

Dengan integrasi parsial diperoleh Γ(𝛼 + 1) = 𝛼Γ(𝛼) Karena ∞

Γ(1) = ∫ 𝑒 −𝑡 𝑑𝑡 = 1 0

Dengan Γ(𝑣 + 1) = 𝑣Γ(𝑣) = 𝑣 Untuk v = 0, 1, 2, 3, ... sehingga : 𝑎0 1 1 == − 2 2(2𝑣 + 2) 2 (𝑣 + 1) 2𝑣 Γ(𝑣 + 1) −1 −1 1 = 𝑣+2 = 𝑣+2 = 2 (𝑣 + 1)Γ(𝑣 + 1) 2 Γ(𝑣 + 2) 1/2𝑣+2 Γ(𝑣 + 2) 𝑎2 1 −1 𝑚 = 2 → 𝑎4 = − 2 =− 2 2+𝑣 2 2(2𝑣 + 2) 22 (𝑣 + 1) 2 Γ(𝑣 + 1) −1 1 = 𝑣+4 = 𝑣+2 2 2Γ(𝑣 + 3) 2! 2 Γ(𝑣 + 3) 𝑎4 1 −1 𝑚 = 3 → 𝑎6 = − 2 =− 2 4+𝑣 2 3(𝑣 + 2) 32 (𝑣 + 2) 2! 2 Γ(𝑣 + 1) 1 𝑚 = 𝑚 → 𝑎2𝑚 = (−1)𝑚 𝑣+2𝑚 𝑚! 2 Γ(𝑣 + 𝑚 + 1) 𝑚 = 1 → 𝑎2 = −

∞

𝑦 = ∑ (−1)𝑚 𝑚=0

1 𝑚! 2𝑣+2𝑚 Γ(𝑣

+ 𝑚 + 1)

𝑥 𝑣+2𝑚

fungsi y yang merupakan penyelesaian persamaan diferensial berbentuk deret tak hingga ini disebut Fungsi Bessel Jenis Pertama orde v dan dinotasikan dengan 𝐽𝑣 (𝑥) 1

𝑚 𝑣+2𝑚 Jadi 𝐽𝑣 (𝑥) = ∑∞ 𝑚=0(−1) 𝑚!2𝑣+2𝑚 Γ(𝑣+𝑚+1) 𝑥

NAMA NIM 1

𝑚 2𝑚 𝐽𝑣 (𝑥) = 𝑥 𝑣 ∑∞ 𝑚=0(−1) 𝑚!2𝑣+2𝑚 Γ(𝑣+𝑚+1) 𝑥

(1)

Nilai bilangan bulat v biasanya dinyatakan dengan n. Jadi untuk n≥0 (−1)𝑚 𝑥 2𝑚

𝐽𝑛 (𝑥) = 𝑥 𝑛 ∑∞ 𝑚=0 22𝑚+𝑛 m!(𝑛+𝑚)! Deret di ruas kanan pada persamaan (1) konvergen mutlak untuk setiap x (uji dengan tes hasil bagi). Fungsi ini merupakan solusi persamaan diferensial (1) untuk v bukan bilangan bulat negatif. Khususnya untuk v = 0, dari (1) diperolah: 𝑥2 𝑥4 𝑥6 𝐽0 (𝑥) = 1 − 2 + 2 2 − 2 2 2 +. . . ., 2 2 4 2 4 6 Yaitu fungsi Bessel orde nol. Unruk memberikan gambaran mengenai aplikasi pemodelan inversi non-linier yang agak lebih kompleks maka pemodelan yang dilakukan 1D dimana pada pemdelan ini model bumi dianggap berlapis horisontal sehingga resistivitas hanya bervariasi terhadap kedalaman dan pendekatan ini dianggap memadai untuk kondisi lengkungan geologi sedimen yang tidak terlalu dalam.Fungsi pemodelan kedepan pada metoda geolistrik dengan model 1-D diformulasikan sebagai persamaan integral Hankel yang menyatakan resistivitas-semu ρa sebagai fungsi dari resistivitas dan ketebalan (ρ, hk) tiap lapisan, k = 1, … , n dan n adalah jumlah lapisan: ∞

𝜌𝑎 = 𝑠 2 ∫ 𝑇(𝜆)𝐽1 (𝜆𝑠)𝜆 𝑑𝜆 0

NAMA NIM

SCRIPT SVD Script function [ t, r, rms ] = VES1dinv ( v2, rhoa, t, r ) ab2 = [ 1 2 5 10 20 50 100 200 500 1000]; rhoa = [ 999 998 979 869 517 130 103 110 120 130]; r = [500 500 500 500]; t = [5 10 15]; v2 = reshape(ab2,length(ab2),1); rhoa = reshape(rhoa,length(rhoa),1); r = reshape(r,1,length(r)); t = reshape(t,1,length(t)); mat = [v2, rhoa]; m = [r t]; rinitial = r; tinitial = t; lr = length(r); lt = length(t); kr = 10e-10; iteration = 1; maxiteration = 500; dfit = 1; clf

while iteration<maxiteration r = m(1:lr); t = m(1+lr:lr+lt); for i = 1:length(ab2) s = ab2(i); [g] = VES1dmod (r,t,s); rhoa1(i,:) = g; end

e1 = [log(rhoa)-log(rhoa1)]; dd = e1; misfit1 = e1'*e1; if misfit1
Keterangan  Fungsi untuk VES inversi 1D  ab2 merupakan nilai vector dari AB/2  rhoa adalah nilai resistivitas semu yang didapat  t merupakan kedalaman tiap layernya  r adalah vector yang dari resistivitas sebenarnya tiap layer  v2 pengulangan baris ab2  rhoa = pengulangan baris rhoa  mat = menjadikan kolom menjadi baris  m = parameter model dimana kasus ini menyatakan r itu adalah 3 interface dan t mennyatakan 4 lapisan  lr adalah banyaknya data r sepanjang bentangan  lt adalah banyaknya data t sepanjang bentangan  kradalah koefisien r (10e-10)  iterasi dimulai dari 1  iterasi maksimalnya 500   r=lr, t=1+lr:lr+lt  Untuk setiap length ab2 dinyatakan s  g adalah rumus ves1dmod  Ketika hasil iterasi kurang dari iterasimaksimum maka hasil setiap rhoa1 = g, dan hasilnya akan masuk ke iterasi 4 karena hasilnya tidak memenuhi iterasi ke 2 dan 3, dan akan dilanjutkan ke iterasi 5  Dan ketika hasilnya melebihi iterasi maksimum maka akan dilanjutkan ke iterasi 2 sampai iterasi 5  e1 menyatakan error pada iterasi yang pertama dan dinyatakan dengan rumus [log(rhoa)log(rhoa1)]

NAMA NIM

% % % %

loglog(ab2,rhoa,'k.',ab2,rhoa1,'k'); axis([1 1000 1 1000]) xlabel('AB/2 (m)'); ylabel('Apparent Resistivity (Ohm-m)'); break end

[A] = jacobian(mat,ab2,r,t,lr,lt,rhoa,rhoa1); [U S V] = svd(A,0); ss = length(S); say = 1; k = 0; while say<ss diagS = diag(S); beta = S(say)*(dfit^(1/say)); if beta<10e-5 beta = 0.001*say; end for i4 = 1:ss SS(i4,i4) = S(i4,i4)/(S(i4,i4)^2+beta); end dmg = V*SS*U'*dd; mg = exp(log(m)+dmg'); r = mg(1:lr); t = mg(1+lr:lr+lt); for i5 = 1:length(ab2) s = mat(i5); [g] = VES1dmod (r,t,s); rhoa4(i5,:) = g; end

e2 = [log(rhoa)-log(rhoa4)]; misfit2 = e2'*e2; if misfit2>misfit1 % ('Beta control') say = say+1; k = k+1; if k == ss-1

 dd mendefinisikan e1  misfit adalah selisih model dengan parameter data yang dirumuskan dengan e1'*e1  jika misfit pada iterasi 1 kurang dari kr maka fungsi langsung berhenti  [A] mendefinisikan iterasi jacobi  [U S V] mendefinisikan svd pada A  ss = banyak data s  say = 1  k=0  diagS mendefinisikan diag(s)  beta mendefinisikan rumus S(say)*(dfit^(1/say))  jika beta kurang dari 10e-5 hasilnya = 0.001*say;  untuk iterasi ke 4 berlaku pada semua data s, yang dirumuskan S(i4,i4)/(S(i4,i4)^2+beta)  dmg mendefinisikan perkalian matriks kernel dengan g yang dirumuskan V*SS*U'*dd  mg merupakan perkalian antara parameter model dengan g yang dirumuskan dengan exp(log(m)+dmg')  r kemudian dinyatakan dengan mg(1:lr)  t kemudian dinyatakan dengan mg(1+lr:lr+lt)  kemudian hasil dari i4 masuk ke i5 untuk mendapatkan error yang lebih sehingga semakin mendekati model yang sebenarnya  i5 berlaku untuk semua data ab2  s akan sama dengan transpose dari hasil matriks i5 dan nilai iterasi ke 4 dari resistivitas semu =g  perhitungan error 2 ini dilakukan untuk mengoreksi hasil rhoa

NAMA NIM

iteration = maxiteration; say = ss+1; end else say = ss+1; m = mg; dfit = (misfit1-misfit2)/misfit1; iteration = iteration+1; a = iteration; if dfit
 Mengeplot grafik hubungan antara AB/2 terhadap resistivitas semu

function [g] = VES1dmod(r,t,s) q = 13; f = 10;

 bagian ini menytakan tiap parameter yang digunakan untuk

 Dan grafik hubungan resistivitas sebenarnya terhadap kedalaman

NAMA NIM

m = 4.438; x = 0; e = exp(0.5*log(10)/m); h = 2*q-2; u = s*exp(-f*log(10)/m-x); l = length(r); n = 1; for i = 1:n+h w = l; v = r(l); while w>1 w = w-1; aa = tanh(t(w)/u); v = (v+r(w)*aa)/(1+v*aa/r(w)); end a(i) = v; u = u*e; end i = 1; g = 105*a(i)-262*a(i+2)+416*a(i+4)746*a(i+6)+1605*a(i+8); g = g-4390*a(i+10)+13396*a(i+12)27841*a(i+14); g= g+16448*a(i+16)+8183*a(i+18)+2525*a(i+20); g = (g+336*a(i+22)+225*a(i+24))/10000; return function [A] = jacobian(mat,ab2,r,t,lr,lt,rhoa,rhoa1) par = 0.1; r2 = r; for i2 = 1:lr r2(i2) = (r(i2)*par)+r(i2); for ii = 1:length(ab2) s = mat(ii); [g] = VES1dmod (r2,t,s); rhoa2(ii,:) = g; end A1(:,i2) = [(rhoa2rhoa1)/(r(i2)*par)]*r(i2)./rhoa; r2 = r; end t2 = t; for i3 = 1:lt t2(i3) = (t(i3)*par)+t(i3); for ii = 1:length(ab2) s = mat(ii); [g] = VES1dmod (r,t2,s); rhoa3(ii,:) = g;

perhitungan model (g) pada tiap lapisan

 di script ini hanya melakukan perhitungan A untuk iterasi ke 2 dan ke 3  ditunjukkan bahwa pada perhitungan iterasi ke 2 dan 3 menggunakan fungsi yang sama, sementara pada iterasi ke 4 dan 5 menggunakan fungsi yang berbeda

NAMA NIM

end A2(:,i3) = [(rhoa3rhoa1)/(t(i3)*par)]*t(i3)./rhoa; t2 = t; end A = [A1 A2]; return

FLOWCHART Membuat lapisan-lapisan atau model perkiraan dengan nilai resistivitas semunya

Memprediksi model dengan matriks kedalaman dan menentukan resistivitas pada tiap layernya.

Membuat matriks s berdasarkan data jarak antar geophone

Menghitung nilai resistivitas semu dengan teori koefisien filter dan menghitung nilai error

Membuat matriks Jacobi

Proses SVD pada matriks Jacobi

Menghitung nilai eigen value

Mengekstraksi matriks eigen value yang sudah didapat

Mendapatkan nilai ∆m