Bab III Metodologi III.1 Perhitungan Waktu Tempuh Menggunakan Metoda Pseudo Bending RayTtracing Ray tracing adalah perunutan lintasan sinar (ray path) antara sumber gempa dengan stasiun penerima. Ray tracing merupakan teknik yang sangat fundamental guna menghitung waktu tempuh dalam memecahkan masalah forward dan inversi model seismologi. Ada beberapa metoda ray tracing antara lain: shooting dengan menggunakan hukum Snellius, dan pseudo bending dengan menggunakan prinsip Fermat. Pemilihan metoda pseudo bending dalam penelitian ini disebabkan oleh waktu komputasi untuk penghitungan waktu tempuh dan lintasan sinar gelombang yang relatif cepat. Metoda pseudo bending tidak secara langsung memecahkan persamaan sinar gelombang, tetapi sebagai penggantinya digunakan minimisasi secara langsung waktu tempuh dengan memberikan gangguan kecil secara bertahap pada lintasan sinar gelombang (Um dan Thurber, 1987). Waktu tempuh sepanjang lintasan sinar antara dua titik, source (i) dan receiver (j) diungkapkan dalam bentuk integral garis sebagai berikut:
j1 ∫ dl i v
(3.1.1)
dimana v adalah kecepatan penjalaran gelombang dan dl segmen lintasan sinar. Tinjau lintasan sinar dalam ruang tiga dimensi, vektor posisi dinyatakan dalam
koordinat
Cartesian (Gambar 3.1). Vektor T dan N masing-masing adalah vektor satuan tangensial dan normal. Vektor satuan normal arahnya selalu ke pusat lengkungan. Vektor satuan T dan N dapat dinyatakan sebagai berikut:
25
sinar
z
N
dr r1
T n
r2 y
x
Gambar 3.1. Lintasan sinar gelombang.
dr dr = =T dr ds d 2r ds 2
(3.1.2)
=N
3.1.3)
Persamaan lintasan sinar dinyatakan dalam persamaan vektor sebagai,
−
d 2r ds 2
=
1 dv dr ∇v − v ds ds
(3.1.4)
Vektor di ruas kiri berlawanan arah dengan vektor satuan N. ∇v merupakan gradien kecepatan dalam ruang. Sedangkan dv/ds juga merupakan besaran gradien kecepatan di dalam arah ds atau lintasan yang dapat dihitung melalui perkalian titik antara ∇v dan T, besarnya sama dengan ∇v cos(α), α adalah sudut antara kedua vektor. Definisikan suatu vektor n yang arahnya juga berlawanan dengan vektor satuan N, vektor ini nantinya berguna dalam menentukan arah gangguan pada saat bending dilakukan pada lintasan sinar.
26
n = −v
d 2r ds
2
dv dr = ∇v − ds ds
(3.1.5)
yang mana,
dv = ∇v.T = ∇v T cos(α) = ∇v 1cos(α) = ∇v cos(α) ds
(3.1.6)
sehingga
n = ∇v − ∇v cos(α )T
(3.1.7)
Langkah pertama penerapan pseudo bending berdasarkan Pers. 3.1.7 di atas adalah menentukan vektor satuan tangensial, kemudian mencari vektor gradien kecepatan dari model kecepatan yang ada, dari kedua vektor ini maka vektor n dapat ditentukan. Bila koordinat sumber (xs,ys,zs) dan stasiun penerima (xr,yr,zr) maka
vektor
satuan T adalah,
T=
(xr − xs , yr − ys , zr − zs )' (xr − xs )2 + (yr − ys )2 + (xr − zs )2
(3.1.8)
Garis lurus dari titik sumber ke penerima ini merupakan lintasan awal sinar. Gangguan pertama dilakukan pada titik tengahnya, di titik tengah ini gradien kecepatan ditentukan dari model kecepatan dengan cara membagi jarak antara sumber dan stasiun penerima menjadi beberapa selang, dari setiap selang gradiennya dihitung kemudian dijumlahkan, hasil penjumlahan ini merupakan gradien rata-rata pada arah garis ini. Gradien rata-rata ini belum tentu merupakan gradien kecepatan di titik tengah di atas, oleh karena itu harus dipilih dalam berbagai arah garis sedemikian hingga diperoleh
27
gradien kecepatan rata-rata maksimum. Gradien maksimum ini merupakan wakil gradien di titik tengah ini. Untuk lebih jelasnya secara grafis dapat dilihat dalam Gambar 3.2.
Gambar 3.2. Penentuan gradien kecepatan maksimum dalam model kecepatan, S=sumber dan R=penerima, garis diputar dalam berbagai arah bidang di dalam ruang 3D.
vk+1Receiver L ∇v cos α vmid
n
R n / n
L vk-1
∇v
vk
Source Gambar 3.3. Geometri gangguan titik tengah sejauh R dalam arah vektor satuan n / n (Um dan Thurber, 1987).
Setelah vektor T, ∇v, dan n dapat ditentukan maka masalah selanjutnya adalah menentukan sejauh mana titik tengah ini dipindahkan dalam arah vektor satuan n/ n. Pergeseran sejauh R ini dapat diturunkan sebagai berikut (Gambar 3.3).
28
Sebagaimana sudah dijelaskan di depan bahwa waktu tempuh gelombang sepanjang lintasan sinar dinyatakan oleh Pers. 3.1.1. Sedangkan di dalam segmen garis seperti gambar di atas menggunakan aturan trapesoidal.
t=
t=
1 2
(L2 + R 2 ) v 1
1 2
(L
2
+
k −1
(
+
k
1 vk + 1
)
C=
(L2 + R 2 ) v1
1 1 1 1 + R 2 + + + vk − 1 vk vk vk + 1
(L2 + R 2 ) 12 v 1
t=
1 1 + vk 2
+
k −1
1 1 + vk + 1 vk
1 1 1 + 2 v k −1 v k +1
(3.1.9)
(3.1.10)
(3.1.11)
(3.1.12)
)2 (C + v-k1). 1
t = L +R 2
2
(3.1.13)
Disini harga vk merupakan fungsi dari (x,y,z), sehingga vk= vk(x,y,z) dapat didekati dengan deret Taylor orde pertama di sekitar (x,y,z)=(xmid,ymid,zmid) menjadi, v (x, y, z) = v (x k
k
mid
,y
(∇v )
k mid
v =v k
v =v k
v =v k
mid
mid mid
+ (∇v )
k mid
mid
,z
mid
) + ...
((x, y, z) − (xmid , ymid , zmid ))
.R
n n =v + (∇v ) . R mid k mid n n
(3.1.14) (3.1.15)
+ (∇vk )mid cos( β )R
(3.1.16)
+v R
(3.1.17)
n
Dimana
29
v = (∇vk )mid cos( β )
(3.1.18)
n
Sudut β adalah sudut antara vektor gradien kecepatan di titik tengah dengan vektor satuan n / n . Untuk mencari harga minimum waktu tempuh T maka,
dT =0 dR
(3.1.19)
(
)
R C + v-k1 + (L + R)
dv-k1 =0 dR
(3.1.20)
Oleh karena
dv-k1 1 dvk =− dR v2 dR
(3.1.21)
k
maka persamaan (3.1.20) menjadi,
( k
)
R.v Cvk + 1 − (L + R)
(Cv )R + (2Cv 3
2 n
)
dvk =0 dR
(3.1.22)
(
)
v R 2 + Cv2mid + vmid R - L2 v = 0 n
mid n
(3.1.23)
Akar-akar persamaan orde tiga di atas dicari dengan cara mereduksi persamaan menjadi orde kedua, sehingga dengan rumus persamaan kuadrat dan dengan hanya memilih satu akar riilnya saja maka R dapat diperoleh.
R
1,2
=
− b ± b2 − 4ac 2a
(3.1.24)
30
yang mana,
a = 2Cv
v ;
mid n
b = (Cv 2
mid
+v
mid
);
c = −L2 v
n
(3.1.25)
Selanjutnya bending antara dua titik terus dilanjutkan dan waktu tempuhnya dihitung, kemudian lintasan satu dengan lintasan lain dibandingkan sampai perbedaan waktu tempuh cukup kecil.
Gambar 3.4. Ray tracing di dalam ruang 3D (kiri atas) di dalam medium dengan model kecepatan yang mengandung anomali positif (+25%). Tampak sinar gelombang berusaha mendekati anomali positif untuk memperoleh waktu tempuh terkecil baik pada irisan horisontal (kanan atas), irisan vertikal barat timur (kiri bawah), maupun pada irisan vertikal selatan utara (kanan bawah).
31
Gambar 3.5. Ray tracing di dalam ruang 3D (kiri atas) di dalam medium dengan model kecepatan yang mengandung anomali negatif (−25%). Tampak sinar gelombang berusaha menjauhi anomali negatif untuk memperoleh waktu tempuh terkecil baik pada irisan horisontal (kanan atas), irisan vertikal barat timur (kiri bawah), maupun pada irisan vertikal selatan utara (kanan bawah). Contoh ray tracing lintasan sinar gelombang dari sumber ke beberapa penerima di Kompleks Gunung Guntur menggunakan metoda pseudo bending di dalam medium yang mengandung anomali positif dan negatif masing-masing dapat dilihat dalam Gambar 3.4 dan Gambar 3.5.
III.2 Hubungan antara Beda Waktu Tiba Fase Gelombang S dan P dengan Jarak Hiposenter Dalam kenyataannya gelombang gempa menjalar di dalam medium yang mempunyai model kecepatan banyak lapis (multi layer velocity) sehingga lintasan gelombang
melengkung
mengikuti
Hukum
Snellius.
Untuk
menyederhanakan
permasalahan dalam menurunkan hubungan antara jarak gempa atau hiposenter (D) dengan beda waktu tiba gelombang S dan P (tsp), maka diambil model kecepatan satu
32
lapis yang bersifat homogen isotropis (dalam segala arah kecepatan sama), sehingga lintasan gelombang dari sumber gempa (F) ke stasiun gempa (S) berupa garis lurus (Gambar 3.6). Perumusan matematikanya dapat diturunkan sebagai berikut.
E
R
D S∗
Vp Vs
Gambar 3.6. Gelombang P dan S menjalar dari S ke R pada t=t0 masing-masing dengan kecepatan Vp dan Vs. Waktu t0 disebut juga origin time atau waktu gempa. Waktu tiba gelombang P di stasiun R adalah tp dan gelombang S adalah ts.
Jarak tempuh baik gelombang P dan S dari pusat gempa S ke stasiun gempa R adalah D,
D = V (t − t )
(3.2.1)
D = Vs (t s − t 0 )
(3.2.2)
p
p
0
Atau
D Vp D Vs
= tp − t0
(3.2.3)
= ts − t0
(3.2.4)
33
Pers. 3.2.4 dikurangi Pers. 3.2.3 maka,
D
−
Vs D(
D(
D(
1 Vs
D
= (t s − t 0 ) − (t p − t 0 )
Vp −
1
) = ts − t0 − tp + t0
Vp
Vp Vp Vs
Vs
−
Vp − V s Vp Vs
Vp Vs
Vp − Vs
Vp
D= (
Vp Vs
(3.2.6)
) = ts − tp
(3.2.7)
) = ts − tp
Vp Vs
D=
(3.2.5)
(3.2.8)
(t s − t p )
(3.2.9)
(t s − t p )
(3.2.10)
− 1)
D = k(t sp )
(3.2.11)
Yang mana,
k=
Vp Vp ( − 1) Vs
(3.2.12)
t sp = t s − t p
(3.2.13)
Jadi ada hubungan linier antara D dan tsp, yaitu D=ktsp. Makin besar harga tsp maka makin jauh sumber gempa tersebut, tetapan k disebut juga sebagai tetapan Omori yang bergantung pada kecepatan Vp, Vs atau Vp/Vs. Untuk lapisan permukaan harga Vp=3 km/det dan Vp/Vs=1.73, maka k=4.1 km/det. Sebagai contoh tsp=3,8 detik maka jarak
34
hiposenternya 15,58 km dari stasiun gempa. Gempa yang mempunyai tsp<3 detik disebut gempa mikro atau gempa vulkanik untuk daerah di gunungapi, 3
4 disebut tektonik jauh. Pusat gempa atau hiposenter dapat ditentukan dari tiga stasiun atau lebih, yaitu mencari perpotongan lingkaran yang berpusat di ketiga stasiun dengan jari-jari masing-masing diambil dari jarak hiposenternya. III.3 Hubungan Waktu Tiba Gelombang P dengan Beda Waktu Tiba Gelombang S dan P Hubungan kedua besaran ini dalam bentuk linier dan berguna untuk mencari ratio Vp/Vs dan waktu gempa t0. Berdasarkan Pers. 3.2.1 dan Pers. 3.2.11 di atas, Vp (t p − t 0 ) = k( t sp )
(3.3.14)
Vp t p − Vp t 0 = k( t sp )
(3.3.15)
Vp t p = k( t sp ) + Vp t 0
(3.3.16)
Vp t p =
Vp ( t sp ) + Vp t 0 Vp −1 Vs
tp =
1 Vp Vs
−1
(3.3.17)
( t sp ) + t 0
(3.3.18)
Bentuk umum liniernya dapat dinyatakan sebagai y=Ax+B, dan x dan y masingmasing adalah tsp dan tp. Masalah selanjutnya adalah mencari tetapan A dan B melalui regresi linier antara besaran tsp dan tp dari setiap stasiun, sehingga Vp/Vs dan t0 dapat dirumuskan sebagai berikut. A=
1 Vp Vs
(3.3.19) −1
atau
35
Vp Vs
=
1 +1 A
(3.3.20)
dan t0 = B
(3.3.21)
Penerapan regresi linier antara tsp dengan tp dapat dilihat dalam salah satu rekaman gempa vulkanik Gunungapi Guntur yang direkam oleh ke tujuh stasiun gempa (Gambar 3.7).
Gambar 3.7. Contoh rekaman gempa vulkanik Gunungapi Guntur oleh tujuh stasiun gempa dengan waktu tiba fase gelombang P dan S masing-masing (kiri). Regresi linier antara beda waktu tiba gelombang P dan S (tsp) dengan waktu tiba gelombang P (tp) menghasilkan rasio Vp/Vs dan waktu terjadinya gempa (t0) (kanan). Berdasarkan data gempa vulkanik Gunung Guntur dari tahun 1995−2007 distribusi harga rasio Vp/Vs dapat dilihat dalam Gambar 3.8.
36
Gambar 3.8. Distribusi harga rasio Vp/Vs data gempa vulkanik Gunung Guntur dari tahun 1995−2007 dominan pada harga 1,8. III.4 Penentuan Hiposenter Menggunakan Metoda Bola Parameter dasar gempa adalah pusat gempa di dalam ruang (x0,y0,z0), waktu terjadinya gempa (t0), dan magnituda gempa. Salah satu cara praktis menentukan pusat gempa menggunakan metoda bola secara grafis analitis dari tiga stasiun gempa. Setiap stasiun dianggap pusat bola dan jari-jari dihitung menggunakan Pers. 3.2.11. Elevasi ketiga pusat bola dianggap terletak pada bidang horisontal yang sama. Perpotongan ketiga bola membentuk tiga bidang tegak dan ketiga bidang tegak akan saling berpotongan satu sama lainnya, ketiga garis potong di bidang horisontal merupakan episenter gempa. Kedalaman gempa dapat ditentukan secara grafis dari lingkaran baru yang dibentuk dari salah satu garis potong sebagai diameternya dan panjang garis dari
episenter sampai ke lingkaran baru dan tegak lurus diameter
merupakan kedalaman gempa (Gambar 3.9).
37
Gambar 3.9. Penentuan hiposenter menggunakan metoda bola. Lokasi stasiun (LGP, PTR, dan CTS) adalah pusat bola dan jarak hiposenter D adalah jari-jari lingkaran. Perpotongan garis potong ketiga lingkaran di E adalah episenter gempa (x0,y0) dan kedalaman gempa adalah panjang garis EF (z0). Dalam hal ini Panjang garis EF 6,48 km dan bidang horisontal di atas terletak 2 km di atas permukaan laut. Oleh karena elevasi referensi z=0 terletak 4 km di atas muka laut maka kedalaman gempa menjadi 8,48 km dari elevasi referensi. III.5 Penentuan Hiposenter Menggunakan Metoda Grid Search Perhitungan hiposenter menggunakan metoda bola tidak begitu teliti karena model kecepatannya hanya satu lapisan kecepatan, sedangkan kondisi medium di lapangan jauh lebih kompleks dan metoda bola hanya melibatkan tiga stasiun gempa. Supaya perhitungan mendekati kondisi medium sebenarnya maka diperlukan metoda yang menggunakan model kecepatan banyak lapis dan melibatkan semua stasiun gempa, metoda itu adalah metoda grid search. Metoda ini merupakan salah satu metoda optimasi dengan menguji secara langsung sekumpulan bakal solusi secara iteratif satu persatu. Pemilihan solusi ditentukan oleh kesalahan terkecil tres antara hasil pengamatan dengan
38
hasil perhitungan. Dalam hal perhitungan hiposenter bakal solusi adalah himpunan koordinat pusat gempa dan waktu gempa, {(x0,y0,z0,t0)i}. Kesalahan terkecilnya adalah perbedaan waktu tiba antara gelombang P hasil pengamatan (observation), (tp)obs, dengan hasil perhitungan (calculation), (tp)cal dari n stasiun gempa.
t res =
1 n
n
(( )obs − (t p )cal )2
∑ tp
(3.5.1)
i =1
Waktu tiba gelombang P hasil perhitungan, (tp)cal, ditentukan berdasarkan waktu tempuh gelombang P dari sumber ke penerima, ttrv, di dalam model kecepatan, v, ditambah waktu terjadinya gempa, t0, yang dihitung menggunakan metoda peseudo bending ray tracing,
(t p )cal = t 0 + t trv n
dl ij
i =1
vi
(t trv ) j = ∑
(3.5.2) (3.5.3)
Besaran vi adalah model kecepatan di dalam model terparameterisasi, indeks j adalah sinar gelombang dari pusat gempa ke stasiun ke-j, dan i adalah elemen volume kei. Model kecepatan untuk perhitungan hiposenter adalah model kecepatan banyak lapis (multi layers) di dalam medium yang telah diparameterisasi ke dalam setiap elemen volume 2x2x2 km3 yang didapat berdasarkan hasil studi tomografi waktu tunda gelombang P Gunungapi Guntur tahun 2002 (Suantika, 2002) dan telah disesuaikan dengan data gempa vulkanik terbaru. Untuk menghemat waktu perhitungan maka ruang uji dibatasi pada ruang volume 4x4x4 km3 dan hiposenter yang di dapat berdasarkan metoda bola dipakai sebagai titik pusat volume, sedangkan waktu gempa diberikan antara t0-dt
39
satu. Hasil perhitungan hiposenter menggunakan metoda grid search dapat dilihat dalam Gambar 3.10.
Gambar 3.10. Hasil hiposenter metoda grid search (bulatan merah) mempunyai tingkat kesalahan (time residual) cukup kecil, kurang dari 0,1 detik. Ruang uji dibatasi pada ruang volume 4x4x4 km3, hasil hiposenter metoda bola (bulatan kuning) dipakai sebagai titik pusat. Gambar ditampilkan masing-masing pada irisan horisontal (kiri atas), irisan vertikal selatan utara (kanan atas), dan irisan vertikal barat timur (kiri bawah). III.6 Perhitungan Model Kecepatan Penjalaran gelombang atau sinar (ray) dari sumber S ke penerima R di dalam medium banyak lapisan dengan model kecepatan yang bervariasi terhadap kedalaman, v=v0+kz, mengalami refraksi atau pembiasan arah penjalaran (Gambar 3.11). Besarnya sudut pembiasan atau sudut pergi i menurut Hukum Snellius tergantung pada sudut awal
40
atau sudut datang i0, nilai kecepatan di lapisan pertama adalah v0, dan kecepatan di lapisan kedua adalah v. Secara matematis dapat dinyatakan sebagai berikut,
Gambar 3.11. Penjalaran gelombang dari sumber S ke penerima R di dalam medium banyak lapis dengan model kecepatan yang bervariasi terhadap kedalaman, v=v0+kz.
sin(i) sin(i 0 ) = =p v v0
(3.6.1)
Harga p disebut juga parameter sinar (ray parameter ) yang harganya tetap untuk setiap sinar. Dari Pers. 3.6.1 dapat diturunkan hubungan sinus, kosinus, dan tangen sudut pembiasan i terhadap parameter sinar p dan kecepatan gelombang v sebagai berikut (Telford et al., 1976),
sin(i) = pv
(3.6.2)
cos(i) = 1 − (pv) 2
(3.6.3)
tan(i) =
pv
(3.6.4)
1 − (pv) 2
41
Untuk menentukan model kecepatan di setiap lapisan telebih dahulu dicari total waktu tempuh gelombang t dari S ke R dan jarak episenter x sebagai berikut. Tinjau elemen sinar ds dengan sudut pembiasan i,
cos(i) =
dz ds
(3.6.5)
cos(i) =
dz vdt
(3.6.6)
dz dt
(3.6.7)
dt =
dz v cos(i)
(3.6.8)
t=∫
dz v cos(i)
(3.6.9)
v cos(i) =
Dan jarak episenter adalah,
tan(i) =
dx dz
(3.6.10)
x = ∫ tan(i) dz
(3.6.11)
Masukkan Pers. 3.6.3 dan Pers. 3.6.4 masing-masing ke dalam Pers. 3.6.9 dan Pers. 3.6.11 sehingga,
t=∫ x=∫
dz
(3.6.12)
v 1 − (pv) 2 pv dz
(3.6.13)
1 − (pv) 2
Masukkan model kecepatan v=v0+kz ke dalam Pers. 3.6.12 dan Pers. 3.6.13 maka,
42
dz
t=∫
(3.6.14)
(v 0 + kz) 1 − (pv 0 + pkz) 2
x=∫
(pv 0 + pkz)dz 1 − (pv 0 + pkz)
(3.6.15)
2
Untuk menyelesaikan integral diatas, maka umpamakan u=pv0+pkz dan turunannya adalah du=pk dz. Maka Pers. 3.6.14 dan Pers. 3.6.15 menjadi,
1 du pk
t=∫
1 u 1− u2 p u.
x=∫
1 du pk
1− u
2
=
=
1 du ∫ k u 1− u2
(3.6.16)
1 udu ∫ pk 1 − u 2
(3.6.17)
Dalam hal ini u=pv0+pkz sama dengan u=p(v0+kz)=pv dan pv=sin(i), sehingga u=sin(i). Terapkan teknik integral menggunakan substitusi trigonometri u=sin(i) dan turunannya, du=cos(i)di, serta persamaan trigonometri cos2(i)+sin2(i)=1, maka Pers. 3.6.16 dan Pers. 3.6.17 menjadi,
t=
1 cos(i)di ∫ k sin(i) 1 − sin 2 (i)
(3.6.18)
x=
1 sin(i)cos(i)di ∫ pk 1 − sin 2 (i)
(3.6.19)
Penyelesaian lebih lanjut integral waktu tempuh dan jarak episenter dilakukan masing-masing. Untuk waktu tempuh adalah,
43
t=
1 cos(i)di 1 cos(i)di 1 di 1 = ∫ = ∫ = ∫ csc(i)di ∫ k sin(i) cos 2 (i) k sin(i)cos(i) k sin(i) k
(3.6.20)
t=
1 ln (csc(i) − cot(i) ) + C k
(3.6.21)
t=
1 1 cos(i) 1 1 − cos(i) + C = ln +C ln − k sin(i) sin(i) k sin(i)
(3.6.22)
Untuk menyederhanakan Pers. 3.6.22 tinjau persamaan trigonometri,
1 − cos 2 (i) = sin 2 (i)
(3.6.23)
(1 + cos(i))(1 − cos(i)) = sin(i).sin(i) (1 − cos(i)) sin(i) = (1 + cos(i)) sin(i)
(3.6.24) (3.6.25)
Dan tinjau juga persamaan trigonometri lain, cos(2i) = cos 2 (i) − sin 2 (i)
(3.6.26)
cos(2i) = 1 − 2 sin 2 (i)
(3.6.27)
cos(2i) = 2 cos 2 (i) − 1
(2.6.28)
Atau berdasarkan Pers. 3.6.27 dan Pers. 3.6.28 didapat,
sin 2 (i) =
1 − cos(2i) 2
(3.6.29)
cos 2 (i) =
1 + cos(2i) 2
(3.6.30)
Dan kuadrat tangen sudut i menjadi,
44
tan 2 (i) =
1 − cos(2i) 1 + cos(2i)
(3.6.31)
tan 2 (i) =
1 − cos(2i) sin(2i) 1 − cos(2i) sin(2i) x = x 1 + cos(2i) sin(2i) sin(2i) 1 + cos(2i)
(3.6.32)
Atau
Substitusikan Pers. 3.6.25 ke dalam Pers. 3.6.32 maka,
1 − cos(2i) 1 − cos(2i) 1 − cos(2i) tan (i) = x = sin(2i) sin(2i) sin(2i)
2
2
1 − cos(2i) tan(i) = sin(2i)
(3.6.33)
(3.6.34)
Bila sudut i diganti dengan i/2 maka,
1 − cos(i) i tan( ) = 2 sin(i)
(3.6.35)
Jadi Pers. 3.6.22 menjadi,
t=
1 i ln tan( ) + C k 2
(3.6.36)
Harga sudut i berkisar antara i0 dan i, maka Pers. 3.6.36 menjadi,
i
t=
i 1 i 1 i ln tan( ) = ln tan( ) − ln tan( 0 ) k 2 i k 2 2 0
45
(3.6.37)
i tan( ) 1 2 t = ln i0 k tan( ) 2
(3.6.38)
Sedangkan penyelesaian jarak episenter adalah,
x=
1 sin(i)cos(i)di 1 sin(i)cos(i)di = ∫ ∫ pk pk 1 − sin 2 (i) cos 2 (i)
(3.6.39)
x=
1 sin(i)cos(i)di 1 1 (− cos(i)) = ∫ ∫ sin(i)di = pk cos(i) pk pk
(3.6.40)
Oleh karena sudut i berkisar dari i0 sampai i maka,
x=
1 [(− cos(i))]ii0 = 1 [(− cos(i)) − (− cos(i 0 ) )] pk pk
(3.6.41)
x=
1 [cos(i 0 ) − cos(i)] pk
(3.6.42)
Penerapan Pers. 3.6.38 dan Pers. 3.6.42 untuk memperoleh model kecepatan satu dimensi (1−D) dilakukan dengan cara memparameterisasi terlebih dahulu tebal lapisan. Dalam penelitian ini digunakan 10 lapisan kecepatan dan setiap lapisan tebalnya 2 km, dua lapisan pertama terletak di atas permukaan laut dan delapan lapisan terakhir terletak di bawah permukaan laut. Nilai awal perhitungan terdiri dari kecepatan awal v0 di lapisan pertama, gradient kecepatan k, dan sudut datang awal i0 di lapisan pertama. Harga p diperoleh dari Persamaan Snellius, p=sin(i0)/v0. Pada lapisan berikutnya dihitung kecepatan v menggunakan persamaan v=v0+kz, yang mana z=2 km. Sudut pergi i di lapisan kedua dihitung menggunakan rumus i=sin-1(pv). Kemudian waktu tempuh dan jarak tempuh sinar dihitung menggunakan Pers. 3.6.38 dan Pers. 3.6.42, begitu seterusnya sampai
46
lapisan kesepuluh. Nilai awal kecepatan gelombang P diberikan secara try and error, sedangkan kecepatan nilai awal gelombang S diperoleh dari nilai awal kecepatan gelombang P dan rasio Vp/Vs=1,8 dari Gambar 3.8. Selanjutnya membandingkan kurva jarak waktu tempuh berdasarkan hasil perhitungan di atas dengan kurva hasil pengamatan yang diperoleh melalui data hiposenter. Lakukan penyesuaian harga v0, k, dan i0 sedemikian rupa sehingga kedua kurva sesuai (Gambar 3.12) sehingga model kecepatan 1−D dapat diperoleh (Gambar 3.13).
Gambar 3.12. Kurva jarak-waktu tempuh hasil pengamatan (bulatan) dan perhitungan (garis) untuk gelombang P (kiri) dan gelombang S (kanan).
47
Gambar 3.13. Model kecepatan 1−D yaitu model kecepatan hanya bervariasi terhadap kedalaman. Model kecepatan gelombang P (biru) dan gelombang S (merah) diperoleh setelah kedua kurva jarak-waktu tempuh sesuai. III.7 Tomografi Waktu Tunda (Delay Time) Gelombang P dan S Waktu tempuh gelombang tij sepanjang lintasan sinar antara dua titik, source (i) dan receiver (j) di dalam medan kecepatan v diungkapkan dalam bentuk integral garis Pers. 3.1.1. Tomografi delay time (waktu tunda) gelombang P atau S adalah selisih waktu tempuh gelombang dari sumber ke−i ke stasiun penerima ke−j yang menjalar dalam medium dengan struktur kecepatan sebenarnya (v0) dengan waktu tempuh gelombang yang menjalar pada struktur kecepatan model referensi (v). Waktu tempuh dalam struktur kecepatan v0 adalah t0ij disebut juga waktu tempuh yang diamati atau observed. Waktu tempuh dalam struktur kecepatan v adalah t disebut juga sebagai waktu tempuh yang dihitung atau calculated. Dalam hal ini t0ij dianggap terdeviasi ke t dan v0 ke v masing-masing sebesar δtij dan δv (Indira dan Gupta, 1998). Dengan mengabaikan suku orde kedua maka didapat hubungan, 48
(t 0 ) = t + δt ij
(3.7.1)
ij
v = v + δv
(3.7.2)
0
j
(t ) = ∫ 0 ij ( L )
ij obs
t=
j
s 0 (r) dl 0
cal
(3.7.3)
obs
j 1 dl = ∫ s(r) dl v(r) ( Lij )
∫ ( Lij )
j 1 dl 0 = ∫ v 0 (r) ( Lij )
(3.7.4)
cal
j
δt = ( t ) − t = ∫ ij 0 ij (L )
s 0 (r) dl 0 −
ij obs
j
∫ s(r) dl ( Lij )
(3.7.5)
cal
s(r) menyatakan perlambatan atau slowness sebagai fungsi posisi di titik (x,y,z) atau vektor r dan dl adalah panjang segmen sinar seismik atau ray-segment length. Bila t0ij di dalam Pers. 3.7.5 diganti dengan Pers. 3.7.1 dengan mempertimbangkan Pers. 3.7.2 maka
δtij menjadi,
δt = ij
j
∆s(r) dl + dt ij
∫ ( Lij )
(3.7.7)
cal
Suku dtij merupakan perubahan waktu tempuh yang disebabkan oleh kesalahan hiposenter, kesalahan data dan suku-suku orde dua atau lebih dari pendekatan yang digunakan. Menurut prinsif Fermat lintasan sinar gelombang di dalam medan kecepatan dan titik ujung tetap maka perubahan waktu tempuh dtij akibat perubahan lintasan sinar dapat diabaikan, dan lintasan (Lij)obs dapat didekati oleh (Lij)cal (Widiyantoro, 2000), Sehingga δtij menjadi,
δt = ij
j
∫ ∆s(r) dl ( Lij )
(3.7.8)
cal
49
Persamaan waktu tunda di atas dalam bentuk model parameterisasi atau model n elemen blok volume dapat dinyatakan dalam, n
∆t j = ∑ ∆s i dl ij
(3.7.9)
i =1
Bila diuraikan untuk semua sinar gelombang ke−j maka Pers. 3.7.9 menjadi persamaan linier simultan,
∆t 1 = ∆s1dl11 + ∆s 2 dl 21 + ... + ∆s n dl n1
(3.7.10)
∆t 2 = ∆s1dl12 + ∆s 2 dl 22 + ... + ∆s n dl n2
(3.7.11)
... ∆t j = ∆s1dl1j + ∆s 2 dl 2j + ... + ∆s n dl nj
(3.7.12)
Dalam bentuk notasi matriks [b]=[A].[x], ∆t 1 dl11 ∆t dl 2 = 12 ... ... ∆t j dl1j
dl 21 ... dl n1 ∆s1 dl 22 ... dl n2 ∆s 2 . ... ... ... ... dl 2j ... dl nj ∆s n
(2.7.13)
Dalam hal ini indeks i adalah blok ke-i, j adalah sinar seismik ke−j, dan dl adalah panjang sinar di dalam elemen blok. Matriks sebelah kiri Pers. 2.7.13 merupakan matriks data, [b] dan matriks pertama sebelah kanan persamaan merupakan Matriks Kernel, [A], dan yang kedua adalah matriks parameter yang dicari, [x]. Matriks parameter dicari menggunakan metoda inversi. Deviasi perlambatan (slowness deviation) ∆s hasil inversi tomografi terhadap model kecepatan referensi (v) dapat dikonversikan ke dalam deviasi kecepatan (velocity perturbation) ∆v. Untuk mencari deviasi kecepatan dapat diturunkan sebagai berikut,
50
s=
1 v
(3.7.14)
∆s = s ∆s =
obs
−s
cal
=
1 −1 v + ∆v v
(3.7.15)
v − v − ∆v − ∆v = 2 2 v + ∆v v v + ∆v v
(3.7.16)
∆s −1 = ∆v v2 + ∆v v
(3.7.17)
∆s − 1 = ∆v v2
(3.7.18)
Oleh karena ∆v mendekati nol atau kecil sekali, maka suku ∆v v di dalam penyebut dapat diabaikan, maka ∆v menjadi,
∆v = −∆s v 2
(3.7.19)
Rumus tomografi waktu tunda gelombang P dalam bentuk matriks dapat diturunkan dari Pers. 3.7.13 dan konversi deviasi perlambatan ke deviasi kecepatan dapat diturunkan melalui Pers. 3.7.19. ∆t p1 dl p11 ∆t dl p2 = p12 ... ... ∆t pj dl p1j
dl p21 ... dl pn1 ∆s p1 dl p22 ... dl pn2 ∆s p2 . ... ... ... ... dl p2j ... dl pnj ∆s pn
∆V = −∆s V 2 p
p
(3.7.20)
(3.7.21)
p
Untuk gelombang S adalah,
51
∆t s1 dl s11 ∆t dl s2 = s12 ... ... ∆t sj dl s1j
dl s21 dl s22 ... dl s2j
... dl sn1 ∆s s1 ... dl sn2 ∆s s2 ... ... ... ... dl snj ∆s sn
∆V = −∆s V 2 s
(3.7.22)
(3.7.23)
s s
Harga Vp dan Vs masing-masing merupakan model kecepatan gelombang P dan S. III.8 Tomografi Q Model Spectral Fitting Studi tomografi menggunakan metoda atenuasi (kebalikan daripada harga Q) cukup banyak dilakukan oleh para ahli. Untuk memahami atenuasi maka perlu ditinjau karakteristik medium. Medium dapat dipandang sebagai sistem osilasi yang mana amplituda getaran akan meluruh terhadap waktu menuju nol. Sistem ini disebut sistem osilasi teredam. Peredam mengubah energi kinetik menjadi energi panas karena ada gaya gesekan. Gelombang gempa menjalar di dalam bumi kehilangan amplituda dan teredam atau teratenuasi, dan atenuasi di dalam seismogram ditunjukkan oleh amplituda getaran terlihat meluruh. Ada dua jenis atenuasi yaitu, atenuasi intrising dan atenuasi hamburan (scattering). Atenuasi intrising adalah hilangnya energi di dalam medium karena ada sifat kentalan (viscousity), porositas (porosity), dan permeabelitas (permeability), makin tinggi kekentalan, porositas, dan permeabelitas makin tinggi atenuasinya. Sedangkan atenuasi hamburan adalah terdistribusinya energi getaran di dalam medium akibat adanya fluktuasi acak kecepatan, efek geometri (pemantulan, pembelokan, dan penghamburan gelombang), dan keragaman (heterogeneity) batuan. Makin beragam batuan maka makin tinggi atenuasinya (Latchman et al., 1996 dan Eberhart et al., 2002). Keberadaan harga Q dalam medium menyebabkan waktu tempuh gelombang terbobotkan dan juga menyebabkan dispersi gelombang yaitu kecepatan gelombang dalam penjalarannya bergantung pada frekuensi (Stein dan Wysession, 2003). Di sini akan dibahas cara mencari harga Q menggunakan analisa spektral dari data domain waktu gempa. Gelombang gempa vulkanik Kompleks Guntur direkam
52
menggunakan seismograf kecepatan short periode (frekuensi natural 1−2 Hz). Data domain waktu (time domain ) kecepatan dan dibentuk oleh empat domain waktu yang berasal dari empat faktor, yaitu: 1. Domain waktu dari sumber getaran (source time function) 2. Domain waktu dari pengaruh medium (propagation effect) 3. Domain waktu dari kondisi batuan di stasiun gempa (site response) 4. Domain waktu dari tanggapan alat seismograf (instrumental response) Keempat faktor ini berkonvolusi membentuk rekaman seismogram yang dapat diamati di permukaan. Dalam domain frekuensi (frequency domain) keempat faktor bergabung dalam bentuk perkalian (Scherbaum, 1996). S(f) = A(f) I(f) R(f) B(f)
(3.8.1)
Yang mana,
S(f)
= Spektrum gempa yang diamati
A(f) = Spektrum sumber I(f)
= Spektrum intrument
R(f) = Spektrum stasiun B(f) = Spektrum medium.
Keempat spektrum masing-masing dapat diterangkan sebagai berikut.
II.8.1 Spektrum Sumber Gempa Model source time function berbentuk segi empat dan transfomasi Fouriernya berbentuk fungsi sinc. Dalam sumbu linier amplituda spektral fungsi sinc akan meluruh pada frekuensi lebih besar daripada fc, A(f)~f-n, n=1,2,3,...,k. Dalam skala log-log mempunyai gradien −n. Pengembangan model lebih lanjut adalah membuat amplitudanya tetap atau datar pada frekuensi lebih kecil daripada fc, (Stein dan Wysession, 2003, Aki
53
dan Richards, 1980, dan Bath, 1974). Amplituda pada f
D(f) =
M o R(θ , ϕ )
f cn n n 4π ρ s V 3 f c + f
(3.8.2)
Di sini,
Mo
= Momen seismik
R
= Pola radiasi gelombang merupakan fungsi daripada azimuth (θ) dan incident angle (φ)
ρ
= Densitas medium di sumber
s
= Jarak hiposenter
V
= Kecepatan gelombang (P atau S)
fc
= Frekuensi sudut sumber (corner frequency)
f
= Frekuensi
n
= Pangkat bilangan bulat pada frekuensi penyebab amplituda spektral meluruh pada frekuensi lebih besar daripada fc (n=2 atau n=3).
Model yang lebih sesuai dengan data pengamatan dalam penelitian ini adalah mengambil harga n=2, sehingga model spektrum menjadi,
D(f) = Ω o
Ωo =
f c2
(3.8.3)
f c2 + f 2
M o R(θ , ϕ )
(3.8.4)
4π ρ s V 3
54
Domain waktu perpindahan d(t) ditransformasikan ke dalam domain frekuensi menjadi D(ω), dan domain waktu kecepatan atau turunan d(t) terhadap waktu ditransformasikan ke domain frekuensi menjadi spektrum kecepatan, A(ω)=|iωD(ω)|. Yang mana, ω=2πf, f adalah frekuensi dan i adalah bilangan imaginer. Sehingga spektrum kecepatan menjadi, A(f) = 2π f D(f)
A(f) = 2π f Ω o
(3.8.5)
f c2
(3.8.6)
f c2 + f 2
III.8.2 Spektrum Atenuasi Medium Pengaruh medium pada penjalaran gelombang dari sumber ke penerima adalah adanya penyerapan energi gelombang oleh medium karena medium tidak elastis sempurna. Fraksi energi yang terserap per siklus osilasi disebut quality factor atau Qfactor (Lay et al., 1995). Kebalikan dari harga Q disebut atenuasi. Adanya nilai Q medium menyebabkan waktu tempuh gelombang dari sumber ke penerima, t, terbobotkan menjadi, t*=t/Q. Model spektrum atenuasi diberikan sebagai berikut.
B(f) = exp(
−π f s 1 s 1 ) = exp(−π f ) = exp(−π f t) VQ QV Q
B(f) = exp(−π f
t ) = exp(−π f t * ) Q
(3.8.7) (3.8.8)
III.8.3 Spektrum Pengaruh Stasiun Spektrum pengaruh stasiun diperoleh dari rekaman gempa tektonik jauh di setiap stasiun pengamatan. Pengaruh medium pada penjalaran gelombang ke jaringan stasiun gempa dianggap sama karena jarak hiposenter jauh lebih besar daripada jarak antar stasiun di dalam jaringan. Jadi perbedaan spektrum di setiap stasiun dianggap berasal dari pengaruh kondisi batuan di bawah stasiun tersebut.
55
Gambar 3.14. Spektrum pengaruh stasiun (kanan atas) diperoleh dari spektrum gempa tektonik jauh yang terekam di stasiun target (Stasiun CSP, kiri atas) dibagi dengan spektrum gempa yang sama direkam di stasiun referensi (Stasiun K74) . Spektrum referensi dipilih berdasarkan kestabilan nilai amplituda pada daerah frekuensi rendah (tengah atas). Spektrum pengaruh stasiun (tengah bawah) digunakan untuk membagi spektrum gempa vulkanik (kiri bawah) dan hasilnya adalah spektrum gempa vulkanik yang telah terkoreksi (kanan bawah). Pilih salah satu spektrum di salah satu stasiun sebagai spektrum referensi, yaitu spektrum yang tidak dipengaruhi kondisi stasiun yang ditandai oleh spektrum dengan amplituda paling stabil pada daerah frekuensi rendah. Spektrum pengaruh stasiun didapat dengan cara membagi spektrum setiap stasiun dengan spektrum referensi. Selanjutnya spektrum hasil bagi ini ini digunakan untuk membagi spektrum gempa vulkanik yang terekam di masing-masing stasiun (Gambar 3.14) (Triastuty, 2001).
III.8.4 Spektrum Pengaruh Instrumen. Dalam penelitian di Kompleks Gunungapi Guntur digunakan dua jenis seismometer kecepatan yaitu seismometer yang memiliki frekuensi natural 1 Hz dan 2
56
Hz. Masing-masing mempunyai tetapan redaman kritis, h=0,7. Spektrum tanggapan instrumen tampak berbeda pada frekuensi lebih kecil daripada 5 Hz. Tanggapan seismometer 1 Hz lebih baik sekitar 4,7 kali daripada seismometer 2 Hz (Gambar 3.15).
Gambar 3.15. Spektrum tanggapan seismometer 1 Hz dan 2 Hz. Amplituda spektrum maksimum dinormalisasi menjadi satu. III.8.5 Spektrum Model. Spektrum model kecepatan yang dipakai dalam penelitian ini adalah perkalian antara spektrum sumber dengan spektrum atenuasi medium.
A(f).B(f) = 2πfΩ 0 M(f) = 2πfΩ 0
f c2 f c2
f c2 f c2 + f 2
+f
2
exp(−πft * )
exp(−πft * )
(3.8.9)
(3.8.10)
57
Spektrum model akan disesuaikan (best fit) dengan spektrum data pengamatan setelah dikoreksi atau dibagi dengan spektrum pengaruh stasiun dan spektrum tanggapan instrumen.
Gambar 3.16. Rekaman gempa vulkanik dalam time domain direkam di Stasiun CSP mempunyai waktu tiba gelombang P dan S cukup jelas (atas). Kurva spektrum best fit gempa vulkanik baik gelombang P dan S (bawah). Time window diambil 1,28 detik atau 128 data dengan sampling rate 0,01 detik. Bila beda waktu tiba gelombang P dan S (SP) lebih kecil daripada 1,28 detik maka mulai data ke 128– (SP/0.01) di beri nilai nol. Kualitas data dapat dilihat dari jarak vertikal antara spektrum noise dan sinyal cukup jauh. Ada tiga parameter yang mempengaruhi bentuk spektrum model (Pers. 3.8.10) diantaranya adalah Ω0, fc, dan t*. Harga Ω0 menyebabkan posisi kurva bergerak dalam sumbu vertikal dan secara fisis ditentukan oleh momen seismik sumber dan jarak hiposenter, harga fc menentukan kurva mulai meluruh pada f>fc, secara fisik merupakan merupakan karakter sumber, dan harga t* merupakan waktu tempuh yang terbobot oleh harga Q. Dalam prakteknya penentuan ketiga harga di atas ditentukan secara trial and
58
error melalui metoda grid search sedemikian hingga kurva spektrum model atau hasil perhitungan cocok dengan kurva spektrum hasil pengamatan (mempunyai faktor kesalahan kecil). Dalam penelitian tomografi ini digunakan analisa spektral baik untuk fase gelombang P maupun fase gelombang S. Pemilihan data gempa didasarkan pada perbedaan jarak yang cukup besar di sumbu vertikal antara spektrum sinyal dengan spektrum noise (Gambar 3.16). Harga fc dari setiap stasiun belum tentu sama harganya walaupun berasal dari satu sumber gempa dan oleh karena fc mewakili efek sumber maka harga fc dari setiap stasiun dirata-ratakan dan perhitungan diulang menggunakan fc rata-rata. III.8.6 Tomografi Q Metoda Spectral Fitting Tomografi Q model ini dapat diturunkan berdasarkan waktu tempuh gelombang dari sumber ke penerima melalui medium kontinyu yang mempunyai kecepatan gelombang V dan faktor kualitas Q dirumuskan sebagai berikut (Bath, 1974). t = t1
dl t = t 0 VQ
t* = ∫
(3.8.11)
Penerapan di dalam medium yang diskrit atau medium terparametrisasi ke dalam n elemen blok volume di daerah penelitian menggunakan rumus yang diturunkan berdasarkan Pers. 3.8.11, sebagai berikut.
n
t *j = ∑
dl ij
(3.8.12)
i =1 Vi Q i
dimana,
tj*
= Waktu tempuh terbobot sinar gelombang ke−j
59
dlij
= Panjang lintasan sinar gelombang ke−j di elemen volume ke−i
Vi
= Kecepatan gelombang di elemen volume ke−i
Qi
= Faktor kualitas di elemen volume ke−i
n
= Jumlah elemen blok volume di daerah penelitian.
Rumus tomografi Q spectral fitting untuk gelombang P dapat diturunkan berdasarkan Pers. 3.8.12 sebagai berikut.
n
t *pj = ∑
dl pij
(3.8.13)
i =1 Vpi Q pi
atau
t *p1 = t *p2 =
dl p11 Vp1Q p1 dl p12 Vp1Q p1
+ +
dl p21 Vp2 Q p2 dl p22 Vp2 Q p2
+ ... + + ... +
dl pn1 Vpn Q pn dl pn2 Vpn Q pn
(3.8.14)
... t *pj =
dl p1j Vp1Q p1
+
dl p2j Vp2 Q p2
+ ... +
dl pnj Vpn Q pn
Dan untuk gelombang S adalah,
n
t *sj = ∑
dl sij
(2.8.15)
i =1 Vsi Q si
atau
60
t *s1 = t *s2 =
dl s11 Vs1Q s1
+
dl s12
+
Vs1Q s1
dl s21 Vs2 Q s2 dl s22 Vs2 Q s2
+ ... + + ... +
dl sn1 Vsn Q sn dl sn2 Vsn Q sn
3.8.16)
... t *sj =
dl s1j Vs1Q s1
+
dl s2j Vs2 Q s2
+ ... +
dl snj Vsn Q sn
Dalam bentuk matriks [b]=[A][x], yang mana [b] matriks data, [A] matriks Kernel, dan [x] matriks parameter yang dicari, persamaan tomografi Q metoda spectral fitting gelombang P dan S masing-masing adalah: dl p11 t *p1 Vp1 * dl p12 t p2 ... = Vp1 * ... t pj dl p1j Vp1
dl s11 t *s1 Vs1 * dl s12 t s2 = ... Vs1 * ... t sj dl s1j V s1
dl p21 Vp2 dl p22 Vp2 ... dl p2j Vp2
dl s21 Vs2 dl s22 Vs2 ... dl s2j Vs2
dl pn1 1 Vpn Q p1 dl pn2 1 ... Vpn . Q p2 ... ... ... dl pnj 1 ... Vpn Q pn ...
dl sn1 1 Vsn Q s1 dl sn2 1 ... Vsn . Q s2 ... ... ... dl snj 1 ... Vsn Q sn
3.8.17)
...
(3.8.18)
Jadi harga Qp dan Qs di setiap elemen blok volume dicari melalui proses inversi masingmasing Pers. 3.8.17 dan Pers. 3.8.18.
61
III.9 Tomografi Q Metoda Spectral Ratio Menurut Pers. 3.8.1 spektrum gelombang P dan S yang terekam dalam bentuk domain waktu di dalam seismogram dapat ditampilkan masing-masing sebagai berikut.
S p (f) = A p (f) I p (f) R p (f) B p (f)
(3.9.1)
Ss (f) = A s (f) I s (f) R s (f) B s (f)
(3.9.2)
Metoda spektral ratio adalah membagi atau membandingkan antara amplituda spektrum gelombang S (Pers. 3.9.1) dengan amplituda spektrum gelombang P (Pers. 3.9.2). Rasio spektral ini dapat menghilangkan efek sumber, efek stasiun, dan efek intrumen. Yang tersisa adalah efek atenuasi dan efek pola radiasi R(θ,φ) yang tercakup di dalam efek sumber (Pers. 3.8.2) (Bath, 1974). Pola radiasi tidak bergantung pada frekuensi f, sehingga persamaan rasio spektral menjadi,
Ss (f) R s (θ , ϕ ) B s (f) = S p (f) R p (θ , ϕ ) B p (f) Ss (f) S p (f)
=
(3.9.3)
* R s (θ , ϕ ) exp(−πft s ) R p (θ , ϕ ) exp(−πft * ) p
(3.9.4)
Bila kedua sisi dari Pers. 20 diubah dalam bentuk logaritma natural, * R s (θ , ϕ ) exp(−πft s ) Ss (f) = ln ln S p (f) R p (θ , ϕ ) exp(−πft * ) p
(3.9.5)
* exp(−πft s ) Ss (f) R (θ , ϕ ) = ln ln + ln s S p (f) R p (θ , ϕ ) * exp( − ft ) π p
(3.9.6)
S (f) R (θ , ϕ ) ln s = ln exp(−πft * ) − ln exp(−πft * ) + ln s S p (f) R p (θ , ϕ ) s p
(3.9.7)
62
S (f) R (θ , ϕ ) ln s = − πft * − − πft * + ln s S p (f) R p (θ , ϕ ) s p
(3.9.8)
S (f) ln s = −π t * − t * f + C S p (f) s p
(3.9.9)
Jadi berdasarkan Pers. 3.9.9 ada hubungan linier antara logaritma rasio spektral dengan beda waktu tempuh terbobot antara gelombang S dengan P. Untuk rekaman gempa vulkanik Gunung Guntur hubungan linier di atas masih bergantung pada frekuensi, artinya hubungan linier hanya berlaku pada lebar pita frekuensi tertentu sehingga perlu dilakukan pemilihan lebar pita. Pemilihan lebar pita frekuensi dilakukan dengan cara mendata frekuensi yang paling sering muncul dari 2686 rekaman gempa kemudian dibuat grafik distribusi frekuensi (Gambar 3.17). Dari grafik ini diputuskan untuk mengambil lebar pita frekuensi antara 10−20 Hz dan jumlah rekaman gempa yang dipakai dalam studi ini menjadi 2520 rekaman. Contoh rekaman gempa vulkanik Gunung Guntur beserta spectral ratio masing-masing dapat dilihat dalam Gambar 3.18 dan Gambar 3.19.
63
Gambar 3.17. Distribusi frekuensi linieritas spectral ratio dari 2520 rekaman gempa vulkanik Gunung Guntur terletak pada pita frekuensi 10−20 Hz.
Gambar 3.18. Rekaman gempa vulknik Gunung Guntur. Analisis spektral menggunakan 128 data (1,28 detik) pada gelombang P dibatasi oleh dua garis merah, Gelombang S dibatasi oleh dua garis biru, dan noise dibatasi oleh dua garis hitam. Bila beda waktu tiba antara S dan P lebih kecil dari pada 1,28 detik maka sisa data diberi nilai nol.
64
Gambar 3.19. Spektrum gelombang P (kiri atas) dan spektrum gelombang S (kanan atas) dari gempa vulkanik yang direkam di stasiun PSC. Spektral noise adalah warna hijau dan letaknya cukup jauh dari spektral gelombang P dan S. Rasio spektral gelombang antara S dan P (kiri bawah) mempunyai kecenderungan linier pada interval frekuensi 5−35 Hz (kiri bawah). Berdasarkan distribusi frekuensi semua rekaman gempa maka dipilih pada interval 10−20 Hz dan harga (ts*−tp*) dapat dihitung dari gradien persamaan garis hasil regresi linier (kanan bawah). Persamaan tomografi rasio spektral untuk gelombang P dan S dapat diturunkan melalui pengurangan Pers. 3.8.15 oleh Pers. 3.8.13 dan menerapkan model hubungan linier antara harga Qp dan Qs rata-rata. Harga Qp dan Qs rata-rata di Kompleks Gunung Guntur dapat diperoleh dari hasil regresi waktu tempuh terbobot menggunakan metoda best fitting dengan data waktu tempuh (t) yang didapat dari hasil perhitungan hiposenter. Model hubungan adalah linier, t=Qt*. Gradien garis regresi merupakan harga Q, sehingga Qp=84,5214, Qs=147,5667, dan Qs/Qp=1,7459 (Gambar 2.18). Jadi Qs=1,7459Qp atau Qp=0.5728Qs.
65
Gambar 3.20. Harga rata-rata Qp dan Qs kompleks Gunung Guntur diperoleh melalui hasil regresi waktu tempuh terbobot (t*) dengan waktu tempuh (t) (atas). Harga Qs=1.7459Qp atau Qp=0.5728Qs (bawah). Jadi persamaan tomografi rasio spektral untuk n elemen blok volume adalah,
n
t *sj − t *pj = ∑
dl sij
i =1 Vsi Q si
n
−∑
dl pij
(3.9.10)
i =1 Vpi Q pi
Untuk rasio spektral gelombang P adalah harga Qs diganti dengan 1,7459Qp,
66
n
t *sj − t *pj = ∑
dl sij
i =1 Vsi 1,7459 Q pi
n
−∑
dl pij
(3.9.11)
i =1 Vpi Q pi
dl s1j dl s2j dl snj t *sj − t *pj = + + ... + Vs1 1,7459 Q p1 Vs2 1,7459 Q p2 Vsn 1,7459 Q pn
dl p1j dl p2j dl pnj + + ... + Vp1Q p1 Vp2 Q p2 Vpn Q pn
−
(3.9.12)
dl s1j dl p1j dl s2j dl p2j + + ... t *sj − t *pj = − − Vs1 1,7459 Q p1 Vp1Q p1 Vs2 1,7459 Q p2 Vp2 Q p2
dl snj dl pnj + − Vsn 1,7459 Q pn Vpn Q pn
(3.9.13)
dl s2j dl s1j dl p1j 1 dl p2j 1 t *sj − t *pj = − + − + ... 1,7459 Vs1 Vp1 Q p1 1,7459 Vs2 Vp2 Q P2 dl pnj 1 dl snj − 1,7459 Vsn Vpn Q pn
+
(3.9.14)
Bila Pers. 3.9.14 diuraikan terhadap semua sinar gelombang (j=1,2,3...) maka bentuk lengkap persamaan adalah,
dl dl dl p11 1 dl p21 1 s11 s21 t *s1 − t *p1 = − + − + ... 1,7459 Vs1 Vp1 Q p1 1,7459 Vs2 Vp2 Q p2 dl pn1 1 dl sn1 + − 1,7459 Vsn Vpn Q pn dl dl dl p12 1 dl p22 1 s12 s22 t *s2 − t *p2 = − + − + ... 1,7459 Vs1 Vp1 Q p1 1,7459 Vs2 Vp2 Q p2
67
(3.9.15)
dl pn2 1 dl sn2 + − 1,7459 Vsn Vpn Q pn
(3.9.16)
dl s1j dl s2j dl p1j 1 dl p2j 1 t *sj − t *pj = − + − + ... 1,7459 Vs1 Vp1 Q p1 1,7459 Vs2 Vp2 Q p2 dl pnj 1 dl snj + − 1,7459 Vsn Vpn Q pn
2.9.17)
Dalam bentuk matriks adalah,
t *s1 * t s2
1 dl dl dl dl p11 pn1 s11 sn1 ... Q p1 − − 1,7459 Vs1 Vp1 1,7459 Vsn Vpn 1 ... ... ... . Q p2 (3.9.18) = dl dl pnj ... dl s1j dl snj p1j 1 ... − − 1,7459 V V 1,7459 V V sn pn s1 p1 Q pn
− t *p1 − t *p2
... * t sj − t *pj
Untuk gelombang S maka harga Qp diganti dengan 0.5728Qs, dan berdasarkan Pers. 3.9.18 rumus tomografinya rasio spektralnya menjadi,
t *s1 * t s2
1 dl dl p11 dl pn1 dl sn1 s11 Q s1 − ... − Vs1 0,5728 Vp1 Vsn 0,5728 Vpn 1 . (3.9.19) ... ... ... = Q s2 dl dl dl dl s1j p1j pnj snj ... V − 0,5728 V ... V − 0,5728 V 1 pn s1 p1 sn Q sn
− t *p1 − t *p2
... * t sj − t *pj
Selanjutnya harga Qp dan Qs di setiap elemen blok volume dicari melalui proses inversi masing-masing Pers. 3.9.18 dan Pers. 3.9.19.
68
III.10 Inversi Tomografi Masalah inversi tomografi dengan menggunakan pendekatan linier adalah memperoleh solusi persamaan linier yang dinyatakan dalam bentuk matriks,
[A][x] = [δt ]
(3.10.1)
yang mana [A] adalah matriks Kernel (m x n) dan elemen matriksnya berupa panjang sinar di setiap blok, [x] adalah vektor solusi (n x 1) mempunyai elemen berupa parameter slowness deviation atau atenuasi yang akan dicari, dan [δ δt] merupakan vektor data (m x 1) berlemen waktu tunda atau waktu tempuh terbobot. Dalam inversi tomografi biasanya jumlah sinar seismik atau persamaannya (m−persamaan) lebih banyak daripada jumlah blok atau bilangan yang tidak diketahui (n−parameter), sehingga persamaan di atas menjadi over determined dan matriks [A] tidak bujur sangkar, hanya matriks bujursangkar yang mempunyai invers matriks, oleh karena itu maka matriks [A] dibuat bujur sangkar dengan cara mengalikan dengan matriks transposnya.
A T A [x] = A T δt
(3.10.2)
Sehingga persamaan di atas menjadi persamaan normal dari persoalan least square linier. Sudah umum dalam tomografi seismik bahwa sinar seismik tidak selalu dapat melintasi semua blok, sehingga matriks [A] banyak mempunyai elemen nol dan matriks [A] sangat renggang atau sparse dan persamaannya menjadi tidak bebas atau tidak independent, kondisi matriks ini disebut under determined. Sebagai akibatnya adalah tidak semua elemen vektor [x] dapat diperoleh jawabannya. Disamping itu data selalu mengandung kesalahan sehingga persamaan menjadi tidak konsisten dan solusinya menjadi tidak tunggal atau non-unique. Untuk itu perlu dilakukan minimisasi
69
norm[x−(ATA)-1(ATδt)] atau minimisasi x−(ATA)-1(ATδt)ξ, yang mana ξ adalah Euclidean norm, untuk problem least square ξ=2. Dalam LSQR pemakaian [ATA] yang banyak memakan waktu secara jelas dapat dihindarkan (Paige dan Saunder dalam Widiyantoro, 2000). Metoda LSQR memecahkan vektor [x] dengan jalan meminimisasikan Ax−δ δtmelalui minimisasi [xTx]. Oleh karena matriks [ATA] mendekati singular maka perkiraan solusi dalam LSQR menggunakan dua jenis damping yaitu explicit damping dan implicit damping. Explicit damping ada dua yaitu minimum norm dan gradient damping. Minimum norm membiaskan solusi hasil inversi ke model awal. Dengan adanya damping maka algoritma LSQR secara iteratif mendekati minimum norm dari solusi least square persamaan [A][ x] =[δ δt] dengan cara menyelesaikan [x] melalui,
[
min ( Ax − δt )T ( Ax − δt ) + λ 2xT x
]
(3.10.3)
yang mana λ2 adalah parameter damping. Gradient damping membiaskan solusi hasil inversi ke suatu model yang smooth. Cara kerja damping ini adalah dengan menggunakan roughening matrix. Roughening matrix mengambil selisih harga deviasi slowness dari suatu blok dengan blok-blok di sekitarnya. Disini berlaku syarat sebagai berikut, N M i 1 2 min ∑ ∑ (x − x j ) i = 1 2Ni j = 1 i
(3.10.4)
Ni adalah jumlah blok di sekitar blok ke−i dan M adalah jumlah semua blok. Damping ini menambahkan sebanyak M persamaan ke sistem persamaan linier, N
i
xi − ∑ x j = 0
(3.10.5)
j=1
70
Yang mana i=1,2,3,...,M. Tujuan penerapan gradient damping adalah untuk mendapatkan tomogram yang relatif smooth sehingga diharapkan dapat memudahkan dalam interpretasi. Dalam proses inversi dengan pendekatan iteratif, maka jumlah iterasi berlaku sebagai implicit damping (Widiyantoro, 1997). Solusi tomografi yang didapat dalam penelitian ini dilakukan dalam satu langkah inversi atau one step inversion.
III.11 Hubungan antara Konstanta Elastisitas Medium Batuan dengan Kecepatan Gelombang P dan S Batuan sebagai medium elastik dengan densitas, ρ, mempunyai parameter elastisitas seperti shear modulus (µ), Young modulus (E), bulk modulus (K), Lame constant (µ dan λ), Poisson’s ratio (σ). Shear modulus didefinisikan sebagai rasio antara tekanan geser (shear stress) dengan perubahan panjang (strain),
Young modulus
didefinisikan sebagai rasio antara tekanan normal (normal stress) dengan perubahan panjang, bulk modulus didefinisikan sebagai rasio antara tekanan normal dengan perubahan volume, Konstanta Lame terdiri dari dua kontanta yaitu kontanta elastisitas geser (µ) yang disebut juga sebagai rigiditas dan kontanta elastisitas volume (λ), dan konstanta elastisitas Poisson’s ratio didefinisikan sebagai rasio antara perubahan panjang (strain) arah transversal dengan perubahan panjang arah longitudinal. Shear modulus, Young modulus, bulk modulus, dan Lame constant mempunyai satuan sama dengan tekanan. Sedangkan Poisson’s ratio tidak mempunyai satuan. Konstanta Lame mempunyai hubungan dengan konstanta elastisitas lainnya sebagai berikut,
λ=
σE (1 + σ )(1 − 2σ )
(3.11.1)
µ=
E 2(1 + σ )
(3.11.2)
71
Begitu pula hubungan antara kecepatan gelombang P (Vp) dan kecepatan gelombang S (Vs) dengan beberapa parameter elastisitas di atas dan dengan densitas (ρ) adalah sebagai berikut (Howell, 1959),
λ + 2µ V = = p ρ V = s
4 K+ µ 3 ρ
(3.11.3)
µ ρ
(3.11.4)
Gelombang P adalah gelombang longtudinal yaitu arah getaran partikel medium searah dengan arah penjalaran gelombang dan gelombang S adalah gelombang transversal yaitu arah getaran partikel tegak lurus dengan arah penjalaran gelombang. Dari persamaan di atas baik gelombang P maupun S tergantung pada shear modulus yang mencerminkan komponen kecepatan gelombang geser (Pers. 3.11.2). Gelombang P juga tergantung pada bulk modulus, sedangkan gelombang S tidak bergantung pada bulk
modulus dan hanya bergantung pada shear modulus sehingga gelombang S disebut juga sebagai shear wave velocity. Kecepatan gelombang yang tidak bergantung pada shear
modulus disebut juga bulk sound velocity dan dapat diturunkan dari Pers. 3.11.1 dengan memberikan shear modulus atau µ=0,
ϕ=
K ρ
(3.11.5)
Besaran φ adalah bulk sound velocity yang mempunyai dimensi sama dengan kecepatan dan hanya bergantung pada bulk modulus dan densitas. Sifat fisik medium dapat dipelajari melalui konstanta elastisitas di atas dan konstanta elastisitas dapat diketahui melalui pengukuran di laboratorium atau melalui pengukuran kecepatan gelombang P dan S di lapangan. Jadi besaran seperti Poisson’s
ratio
dan bulk sound velocity dapat diukur menggunakan gelombang P dan S
72
berdasarkan perumusan yang diturunkan dari Pers. 3.11.1 sampai Pers. 3.11.4 sebagai berikut, λ + 2µ σE λ 2σ ρ (1 + σ )(1 − 2σ ) p = = +2= +2= +2 µ E µ (1 − 2σ ) V2 s 2(1 + σ ) ρ V2
V2 p
V2
−2=
2σ (1 − 2σ )
3.11.6)
(3.11.7)
s
V2 V2 V2 p p p − 2 (1 − 2σ ) = − 2 − 2σ − 2 = 2σ V2 V2 V2 s s s
(3.11.8)
V2 V2 V2 V2 p p p p − 2 = 2σ − 2 + 2σ = 2 − 2 + 2 σ = 2 − 2 σ V2 V2 V2 V2 s s s s
(3.11.9)
V2 V2 p p − 2 = 2 − 1σ V2 V2 s s
(3.11.10)
Dan Poisson’s ratio (σ) menjadi,
73
V 2 p − 2 Vs σ= V 2 p 2 − 1 Vs
(3.11.11)
Poisson’s ratio merupakan fungsi dari rasio Vp/Vs. Poisson’s ratio juga merupakan fungsi kontanta elastisitas geser (µ) dan kontanta elastisitas volume (λ), dapat diturunkan berdasarkan Pers. 3.11.6,
λ 2σ +2= +2 µ (1 − 2σ )
(3.11.12)
λ 2σ = µ (1 − 2σ )
(3.11.13)
λ (1 − 2σ ) = 2σ µ
(3.11.14)
λ λ − 2σ = 2σ µ µ
(3.11.15)
λ λ λ = 2σ + 2σ = 2σ 1 + µ µ µ
(3.11.16)
λ µ
= 2σ
(3.11.17)
λ = 2σ µ+λ
(3.11.18)
λ 1+ µ
σ=
λ 2(µ + λ )
(3.11.19)
74
Untuk fluida harga µ=0, maka Poisson’s ratio mencapai harga maksimum, σ =0,5. Sedangkan besaran bulk sound velocity hanya bergantung pada bulk modulus dan densitas dapat diturunkan dari Pers. 3.11.3, Pers. 3.11.4, dan Pers. 3.11.5 sebagai berikut,
4 K+ µ 3 = K + 4 µ =ϕ2 + 4 µ V2 = p ρ ρ 3ρ 3ρ V2 = ϕ 2 + p
4µ 3ρ
(3.11.20) (3.11.21)
dan, V2 = s
µ ρ
(3.11.22)
Substitusikan Pers. 3.11.22 ke Pers. 3.11.21,
V2 = ϕ 2 + p
4 2 V 3 s
(3.11.23)
4 3
s
(3.11.24)
4 3
s
ϕ 2 = V2 − V2 p
ϕ = V2 − V2 p
(3.11.25)
III.12 Magnituda Gempa Vulkanik Gunung Guntur Secara umum besar kecilnya energi suatu gempa sangat bergantung pada amplituda maksimum (peak to peak) dan durasi gempa. Makin besar dan lama gempa itu maka makin besar energinya. Cara menentukan lama gempa adalah mencari perbedaan antara waktu tiba gelombang P dan waktu akhir gempa dari coda wave yang berakhir pada latar belakang noise yang sama dengan sebelum gelombang P tiba. Dalam hal ini energi gempa vulkanik Gunung Guntur dihitung menggunakan besaran di atas. Gempa vulkanik mempunyai SP antara 0−1.6 detik atau jarak dari stasiun gunungapi ke pusat gempa vulkanik sekitar 0−5,4 km (ambil kecepatan gelombang P di
75
lapisan permukaan Vp=2,5 km/detik, Vp/Vs=1,73, dan konstanta Omori k=3,424 km/detik). Magnituda dan energi gempa vulkanik yang mempunyai jarak hiposenter antara 0-5 km dirumuskan oleh Richter dan Guthenberg (Richter et al., 1958) sebagai berikut,
M = log A + 1,4
(3.12.1)
r
Ar =
A pp 2
x
2800 I
(3.12.2)
Harga M adalah magnituda dalam skala Richter, App adalah hasil pengukuran amplituda peak to peak di seismogram dalam mm, Ar merupakan amplituda getaran tanah dalam mikron, dan I merupakan pembesaran seismograf merupakan fungsi dari frekuensi. Angka 2800 merupakan pembesaran seismograf mekanik standard Richter. Untuk seismograf sistem telemetri PS−2 Kinemetrics yang dipakai oleh stasiun permanen CTS di Gunung Guntur pada frekuensi puncak 12,5 Hz mempunyai pembesaran sekitar 110.000 kali. Stasiun CTS digunakan sebagai stasiun standard di dalam menghitung magnituda gempa karena rekamannya stabil dan bisingnya kecil. Kemudian energi (E) gempa vulkanik dalam erg dihitung menggunakan rumus sebagai berikut, log E = 11,8 + 1,5 M
(3.12.3)
Menurut Lee (Lee et al., 1972) menunjukkan ada hubungan antara magnituda (M) dengan logaritma lama gempa (t) secara linier, dan dirumuskan seperti di bawah ini, M = C + C log(t) 1
(3.12.4)
2
yang mana, C1 dan C2 merupakan tetapan.
76
Gambar 3.21. Regresi linier logaritma lama gempa dengan magnituda skala Richter gempa vulkanik di kompleks Gunung Guntur. Besarnya tetapan di atas untuk gempa vulkanik Gunung Guntur dicari melalui regresi linier antara logaritma lama gempa dengan magnituda skala Richter yang besarnya adalah, C1=−0,870 dan C2=1,049 (Gambar 3.21). Sedangkan magnituda gempa skala Richter dihitung menggunakan Pers. 3.12.1 dan Pers. 3.12.2. Selanjutnya perhitungan magnituda gempa vulkanik Gunung Guntur menggunakan Pers. 3.12.4.
77
