Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
DOI: doi.org/10.21009/SPEKTRA.032.05
IMAGING OF 3-D SEISMIC TOMOGRAPHY FOR INTERNAL STRUCTURE UNDER THE MOUNTAIN MERAPI USING THE LOTOS-10 SOFTWARE
Ayu Wita Sari 1, a) , Gede Bayu Suparta 2, b)
1 STIKES Guna Bangsa Yogyakarta, Jl. Padjajaran No.Kel, Ngringin, Condongcatur, Kec. Depok, Kabupaten Sleman, Daerah Istimewa Yogyakarta 55283
2 Jurusan Fisika FMIPA UGM, Bulaksumur Yogyakarta, Senolowo, Sinduadi, Mlati, Kabupaten Sleman, Daerah Istimewa Yogyakarta 55281
Email: a) ayu.0221@gmail.com, b) gbsuparta@gmail.com

## Abstrak
Telah  dilakukan  penelitian  pencitraan  tomografi  seismik  3D  untuk  stuktur  internal  di  bawah  Gunung Merapi  dengan  empat  stasiun  pencatat  gempa  dan  gempa  vulkanik  sebagai  sumber  sinar  gelombang. Penelitian  ini  menggunakan  perangkat  LOTOS-10  (Local  Tomography  Software  10)  untuk  inversi tomografi seismik 3D. Karakteristik medium bawah Gunung Merapi dapat digambarkan oleh parameter fisis seperti kecepatan gelombang primer dan sekunder. Hasil pengolahan data seismograf menunjukkan metoda  tomografi  seismik  dapat  mengungkap  struktur  bawah  permukaan  Gunung  Merapi  melalui distribusi  anomali  deviasi  kecepatan  dan  Vp/Vs  ratio.  Kualitas  citra  yang  dihasilkan  dengan menggunakan  gelombang  primer  lebih  jelas  resolusinya  dan  waktu  yang  digunakan  lebih  efesien, sehingga  dapat  digunakan  sebagai  informasi  mitigasi  bencana  sebelum  gempa  erupsi  terjadi.  Daerah anomali  negatif  yang  diperoleh  terletak  di  bawah  puncak  Gunung  Merapi  pada  kedalaman  3  -  5  km mempunyai  karakter  fisis  yaitu  zona  lemah,  kurang  kompak,  panas  dan  heterogen.  Daerah  anomali tersebut  dapat  diinterpretasikan  sebagai  keberadaan  zona  materi  panas  yang  berasosiasi  dengan  sisa dapur magma dangkal.
Kata-kata kunci : Gunung api Merapi, sifat fisis, tomografi seismic, lotos-10.

## Abstract
3D seismic tomography imaging research conducted for internal structures under Merapi Volcano with four  earthquake  recording  and  volcanic  earthquake  stations  as  a  source  of  wave  rays.  This  study  used LOTOS-10 (Local Tomography Software 10) for 3D seismic tomography  inversion. Characteristics of the  medium  under  Merapi  Volcano  described  by  physical  parameters  such  as  primary  and  secondary wave  velocities.  The  result  of  seismograph  data  processing  shows  seismic  tomography  method  can reveal  the  subsurface  structure  of  the  Merapi  Volcano  through  the  distribution  of  deviation  anomaly speed and Vp / Vs ratio. Image quality generated by using primary wave more clearly the resolution and time  used  more  efficient,  so  that  can  be  used  as  disaster  mitigation  information  before  earthquake eruption happened. The negative anomaly area obtained under the peak of Merapi Volcano at a depth of 3  -  5  km  has  the  physical  characteristics  of  weak,  less  compact,  hot  and  heterogeneous  zones.  The anomalous region can interpret as the existence of a zone of heat material associated with the rest of the shallow magma kitchen.
Keywords : Mount of Merapi, physical character, Seismic of Tomography, lotos-10.
Received: 29 November 2017 Revised: 29 April 2018
Accepted: 18 August 2018 Published: 30 August 2018
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA

## PENDAHULUAN
Aktivitas gunung api merupakan kejadian alami yang sulit diduga. Indikasi aktivitas suatu gunung api  antara  lain  adanya  pertumbuhan  kubah  lava  dan  letusan  paroksimal  yang  disertai  awan  panas. Aktivitas  Gunung  Merapi yang  membahayakan  adalah  letusan  yang  disertai  awan  panas,  semburan gas  beracun,  aliran  lava  pijar  ataupun  bahaya  sekunder  sebagai  efek  letusan.  Letusan  yang  disertai awan panas sukar diperkirakan sehingga banyak yang tidak terselamatkan pada kejadian tersebut.
Berdasarkan  permasalahan  yang  diuraikan  di  atas,  maka  memotivasi  para  peneliti  pendahulu untuk  menemukan  metode-metode  dalam  penyelesaian  permasalahan  tersebut.  Publikasi  tentang metode pencitraan struktur 3-D kecepatan gelombang seismik pertama kali dikenalkan oleh Aki dan Lee  [1].  Metode  ini  selanjutnya  dikenal  dengan  metode  tomografi  seismik.  Tomografi  seismik merupakan metode pemetaan bawah permukaan yang belum lama digunakan dalam eksplorasi migas dan  juga  sudah  lama  dipakai  di  bidang-bidang  lain  di  luar  seismik  seperti,  radio  astronomi,  dunia kedokteran dan lainnya
Metode  tomografi  yang  telah  berhasil  digunakan  untuk  mengungkap  citra  struktur  litosfer  bumi yang menunjam ke dalam lapisan mantel menggunakan gelombang gempa bumi tektonik global juga telah  dicoba  diterapkan  di  daerah  yang  lebih  kecil  yaitu  daerah  gunung  api.  Penerapan  metoda tomografi seismik di daerah gunung api diharapkan dapat memberikan gambaran bawah permukaan secara lebih rinci sehingga dapat melengkapi model bawah permukaan lainnya.
Simulasi tomografi seismik dengan model Gunung api Merapi dengan program LOTOS-10 ( Local Tomography Software 10) digunakan untuk memudahkan proses pengkajian suatu metode tomografi seismik  di  daerah  gunung  api.  Sehingga  dengan  simulasi  tomografi  seismik  yang  dihasilkan diharapkan mampu menggambarkan keaadaan di bawah Gunung Merapi dengan lebih rinci. Program ini  sebelumnya  telah  diuji  coba  pada  gempa  tektonik  secara  global  [4]  untuk  itu  penulis  ingin melakukan percobaan pada gempa Vulkanik dengan objek Gunung api Merapi.

## METODE
Penelitian  yang  dilakukan  yaitu  pengolahan  data  gempa  vulaknik  dari  4  seismograf  yang terpasang  di  sekitar  puncak  Gunung  Merapi  menggunakan  metode  inversi  tomografi  seismik.  Data dikumpulkan  dari  BPPTK  (Balai  Peyelidikan  dan  Pengembangan  Teknologi  Kegunungapian) Yogyakarta  dari  tahun  2007  –  2009  [2].  Pengolahan  data  secara  umum  menggunakan software LOTOS-10  dan  Hypo71  serta  SURFER  9.  Tahapan  penelitian  pencitraan  tomografi  seismik  3D  di bawah Gunung Merapi secara umum diperlihatkan pada (Gambar 1).
Pengolahan  pertama  yaitu  menentukan  waktu  datang  ( arrival  time )  gelombang  P  (primer)  dan gelombang  S  (sekunder)  dari  setiap  stasiun  pada  satu  kejadian  gempa.  Data  gempa  vulkanik  yang diperoleh  yaitu  sebanyak  67  kejadian  dan  terdapat  488  sinar  gelombang.  Selanjutnya memperhitungkan  hiposenter  dari  setiap  kejadian  gempa  menggunakan Software  Hypo71 [6]. Parameter  yang  digunakan  untuk input pada  program  ini  yaitu  daftar  lokasi  stasiun,  daftar  model kecepatan  lapisan  berdasarkan  referensi,  daftar  waktu  tiba  gelombang  P  dan  S  beserta  karakteristik gelombang.
Masalah  tomografi  seismik  sesungguhnya  adalah  masalah  non-linier.  Selama  ini  pemecahan masalah  tomografi  seismik  dilakukan  dengan  pendekatan  linier,  yaitu  dengan  melakukan penyederhanaan  masalah  non-linier  menjadi  masalah  linier.  LOTOS  melakukan  penyelesaian tomografi  dengan  dua  tahap  pekerjaan,  yaitu  optimisasi  kecepatan  1D  dan  inversi  tomografi  dalam 3D [3]. Untuk optimisasi kecepatan 1D, yaitu pencarian model kecepatan 1D terbaik dan preliminer lokasi  sumber  dilakukan  dengan  langkah  awal  perhitungan  daftar  waktu  tempuh  model  1D.  Tahap berikutnya menghitung waktu tempuh yang didasarkan interpolasi bilinier dari nilai yang didapatkan dalam  daftar  referensi  untuk  menghitung  lokasi  sumber  dengan  menggunakan  metode grid  search . Matriks kernel yang menggambarkan variasi kecepatan dalam setiap kedalaman pada waktu tempuh setiap sinar kemudian dikomputasi. Inversi matriks dilakukan untuk anomali kecepatan gelombang P dan  S  dalam  1D,  dan  koreksi  parameter  sumber  (dx,  dy,  dz,  dt),  dengan  melakukan  metode  LSQR [8], dengan damping yang dikontrol oleh blok smoothing khusus.
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
Hasil  yang  telah  didapatkan  dalam  pekerjaan  pertama,  kemudian  digunakan  dalam  inversi tomografi 3D, tahap awal dalam inversi ini, yaitu penentuan lokasi sumber dalam  model kecepatan 3D. Jika lokasi sumber dalam model 1D menggunakan metode grid search , maka dalam 3D [4, 5]. Hasil  yang  telah  didapatkan  dalam  pekerjaan  pertama,  kemudian  digunakan  dalam  inversi tomografi 3D, tahap awal inversi ini, yaitu penentuan lokasi sumber dalam model kecepatan 3D. Jika lokasi  sumber  dalam  model  1D  menggunakan  metode grid  search ,  maka  dalam  3D,  menggunakan metode  gradien  dari pencarian nilai ekstrim  ( goal  function ) fungsi G .  Parameterisasi  grid  dibangun pada  iterasi  pertama  saja,  anomali  kecepatan  dikomputasi  dalam  node  yang  terdistribusi  dalam volume  studi  menurut  kepadatan  sinar,  terhubung  satu  dan  yang  lainnya  melalui  blok smoothing dalam  inversi.  Orientasi  grid  yang  berbeda  digunakan  untuk  menghilangkan  pengaruh  node  yang terdistribusi tidak  merata, dan  di- stack .  Perhitungan derivatif  matriks  awal dilakukan  menggunakan jejak sinar yang dikomputasi setelah lokasi sumber dalam model 3D, yang berisi parmeter kecepatan, koreksi  sumber  ( dx,  dy,  dz,  dt )  dan  koreksi  stasiun.  Amplitudo  dan  kehalusan  dari  penyelesaian dikontrol oleh dua blok tambahan. Kenaikan dari bobot blok ini menyebabkan penurunan amplitudo dan  kehalusan  penyelesaian.  Inversi  untuk  keseluruhan  matriks sparse dilakukan  dengan  metode LSQR ( Least Square Regresion ). Hasil yang didapat dari parameterisasi grid-grid direratakan dalam satu model 3D yang kemudian digunakan untuk iterasi selanjutnya.
GAMBAR 1. Diagram alur penelitian pencitraan tomografi seismik di bawah Gunung Merapi
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA

## HASIL DAN PEMBAHASAN
Hasil citra tomogram dari software LOTOS-10 ini dipresentasikan menjadi dua irisan yaitu irisan horizontal  dan  irisan  vertikal.  Masing-masing  irisan  akan  menampilkan  sinar  gelombang  dengan cakupan yang bervariasi. Cakupan sinar bahwa ray covery dan ray density di bawah kawah Gunung Merapi cukup baik sehingga resolusi tomogram di daerah-daerah ini juga cukup tinggi. Irisan vertikal dilakukan dalam arah Barat-Timur dan Selatan-Utara. Adapun irisan horizontal ditampilkan dengan interval kedalaman 2 km, 4 km, 6 km, 8 km, 10 km di bawah elevasi referensi. Model tomografi yang ditampilkan terdiri dari:

## Model Tomogram Waktu Tunda Gelombang P (Primer)
Hasil  penelitian  yang  didapatkan  adalah  citra  anomali  berdasarkan  waktu  tunda  gelombang  P ditampilkan  dalam  besaran  presentasi  deviasi  kecepatan  gelombang  P  (d V p ),  satuan  besaran  deviasi ini  adalah  km/detik.  Harga  anomali  deviasi  kecepatan  kecepatan  berkisar  -20  %  sampai  +10  % terhadap model kecepatan gelombang P. Dalam pembacaan warna untuk warna merah tua mewakili anomali negatif maksimum dan warna biru tua mewakili anomali positif maksimum.
Keenam  gambar  yang  dihasilkan  adalah  anomali  nagatif  terlihat  jelas  pada  irisan  bidang horizontal z =  0  km  dan z =  -2  km  sedangkan  pada z =  -4  km, z =  -6  km, z =  -8  km, z =  -10  km anomali  mulai  kurang  teresolusi.  Pada  gelombang  P  untuk  irisan  vertikal  arah  Barat-Timur  terlihat jelas  pada  bidang  2A-2B  dan  tampilan  anomali  negatifnya  tampak  semakin  berkurang  pada  sisi kanan,  sedangkan  untuk  irisan  vertikal  arah  Selatan-Utara  bidang  1A-1B  anomali  deviasi kecepatannya  positif  dari  pada  bidang  2A-2B,  semakin  ke  selatan  anomali  negatifnya  semakin berkurang  karena terdapat warna  biru  muda  namun  dipresentasikan  sampai  kedalaman  15  km  yang ditunjukkan pada Gambar 2 .
(a) (b)
GAMBAR 2. Irisan vertikal (a). 1A-1B dan (b). 2A-2B anomali deviasi kecepatan gelombang P.
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
Kedalaman (z =  -10 - 1 km)
Kedalaman (z = 1 - 3 km)
Kedalaman (z = 3 - 5 km)
Kedalaman (z = 5 - 7 km)
Kedalaman (z = 7 - 9 km)
Kedalaman (z = 9 - 11 km)
GAMBAR 3. Irisan horizontal tomogram anomali deviasi kecepatan gelombang P

## Model Tomogram Waktu Tunda Gelombang S (sekunder)
Citra anomali berdasarkan waktu tunda gelombang S juga ditampilkan dalam presentase besaran deviasi  kecepatan  gelombang  S  (dVs),  satuan  besaran  deviasi  ini  adalah  km/detik.  Harga  anomali deviasi  kecepatan  berkisar  antara  -20%  sampai  +18%  terhadap  model  kecepatan  gelombang  S. Warna merah tua mewakili anomali negatif maksimum dan warna biru tua mewakili anomali positif maksimum.
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
Citra  tomogram  yang  dihasilkan,  anomali  negatif  dapat  dilihat  jelas  hanya  pada  irisan  bidang horizontal z = 0 km, sedangkan pada z = -2 km, z = -4 km, z = -6 km, z = -8 km, z = -10 km anomali negatif  tidak  teresolusi  dengan  baik.  Anomali  negatif  pada  z  =  0  km  dominan  mengumpul  dengan luasan  kurang  lebih  30  km2  yang  ditampilkan  pada Gambar  5 .  Pada  gelombang  S  untuk  irisan vertikal arah Barat-Timur terlihat jelas pada bidang 2A-2B dan tampilan anomali negatifnya tampak semakin berkurang pada sisi kanan, sedangkan untuk irisan vertikal arah Selatan-Utara bidang 1A-1B anomali deviasi kecepatannya positif dari pada bidang 2A-2B, semakin ke selatan anomali negatifnya semakin  berkurang  karena  terdapat  warna  biru  muda  namun  dipresentasikan  sampai  kedalaman  15 km yang ditunjukkan pada Gambar 4 .
(a) (b)
Gambar 4 Irisan vertikal (a). 1A-1B dan (b). 2A-2B anomali deviasi kecepatan gelombang S.
Kedalaman (z = -10 – 1 km)
Kedalaman (z = 1 – 3 km)
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
Kedalaman (z = 3 – 5 km)
Kedalaman (z = 5 – 7 km)
Kedalaman (z = 7 – 9 km)
kedalaman (z = 9 – 11 km)
Gambar 5 Irisan horizontal tomogram anomali deviasi kecepatan gelombang S.

## Model Vp/Vs Ratio
Berdasarkan  model  Vp/Vs  yang  dihasilkan  untuk  citra  tomogram  ditampilkan  tanpa  satuan, besaran  ini  bergantung  pada  kecepatan  gelombang  Vp  dan  Vs.  Besaran  ini  mempunyai  harga  dari 1,60-1,88.  Warna  merah  tua  mewakili  Vp/Vs  ratio  maksimum  dan  warna  biru  tua  mewakili  Vp/Vs ratio minimum. Untuk selanjutnya akan dibahas Vp/Vs ratio minimum  yang  merupakan anomali di sekitar Vp/Vs ratio  yang  lebih besar  karena  anomali ini lebih  mencerminkan  daerah lemah (kurang padat) di bawah permukaan yang berkaitan langsung dengan kegiatan vulkanik Gunung Merapi. Pada irisan bidang horizontal z = 0 km dapat terlihat jelas dan paling luas area Vp/Vs ratio minimumnya. Untuk  Vp/Vs  ratio  minimum  pada  irisan  bidang  horizontal  selanjutnya  mulai  berkurang  dan  kecil, bahkan tak terlihat untuk Vp/Vs ratio minimum pada irisan bidang horizontal z = 8 km dapat dilihat pada Gambar 7 .
Pada irisan vertikal terdiri dari dua bagian yaitu Selatan-Utara (1A-1B) dan Barat-Timur (2A-2B) terlihat  berbeda  luasan  daerah  Vp/Vs  ratio  minimumnya.  Hasil  citra  yang  diperoleh  untuk  bidang vertikal  arah  Selatan-Utara  Vp/Vs  ratio  minimum  terlihat  kurang  jelas  (berwarna  biru  muda),  pada arah  ini  tampak  Vp/Vs  ratio  minimum  sampai  di  kedalaman  15  km  seperti  yang  ditampilkan  pada Gambar 6 .
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
(a) (b)
Gambar 6. Irisan vertikal bidang (a). 1A-1B dan (b). 2A-2B tomogram Vp/Vs ratio .
Kedalaman (z = -10 - 1 km)
Kedalaman (z = 1 - 3 km)
Kedalaman (z = 3 - 5 km)
Kedalaman (z = 5 - 7 km)
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
Kedalaman (z = 7 - 9 km) Kedalaman (z = 9 - 11 km)
Gambar 7 Irisan horizontal tomogram Vp/Vs ratio .

## Interpretasi Tomogram Bawah Permukaan Gunung Merapi
Model  tomografi  gelombang  P  dan  S  sama-sama  menunjukkan  deviasi  negatif  terhadap  model kecepatan,  deviasi  negatif  maksimum  mencapai  -20%.  Deviasi  negatif  atau  anomali  negatif kecepatan  menurut  studi  terdahulu  sangat  berkaitan  dengan  sumber  panas  di  bawah  permukaan  [9] atau adanya daerah lemah (sesar atau struktur bawah permukaan hasil kegiatan tektonik) [1, 10]. Dari hasil penelitian ini terlihat jelas bahwa tomografi gelombang P lebih jelas atau lebih teresolusi bila dibandingkan dengan anomali negatif model tomografi gelombang S. Ini berarti anomali gelombang P lebih negatif daripada anomali gelombang S.
Pada  model  tomografi  Vp/Vs  ratio  terlihat  juga  anomali  negatif  sangat  jelas  dan  terletak  pada posisi yang sama dengan model gelombang P dan S, harga Vp/Vs ratio terendah di daerah anomali di bawah kawah Merapi pada kedalaman z = 0 adalah Vp/Vs = 1,60. Berdasarkan hasil studi di Gunung Nevado  del  Ruiz  (Columbia)  harga  Vp  dan  Vp/Vs  yang  rendah  di  bawah  kawah  yang  aktif berasosiasi dengan sistem geotermal yang didominasi oleh uap [7].
Model bawah permukaan Gunung Merapi dapat diinterpretasikan berdasarkan anomali tomogram dan distribusi pusat gempa vulkanik. Karakter fisis anomali adalah daerah lemah, heterogen, kurang kompak, panas, mengandung uap air tetapi tidak jenuh, dan kemungkinan di pusat anomali sebagian kecil ada yang bersifat partial melting.
Sumber gempa vulkanik di bawah Kawah Merapi cukup dangkal yaitu 0–2 km di bawah elevasi referensi, di bawah kedalaman ini merupakan zona aseismik yang diinterpretasikan sebagai kantong magma  terbesar  dan  dangkal.  Sementara  pada  kedalaman  2-15  km  merupakan  saluran  magma  dari kantong  magma  yang  paling  dalam.  Pada Gambar  8 merupakan  interpretasi  dari  bawah  Gunung Merapi  yang  dilakukan  penelitian  tomografi  jawa  tengah  dan  menemukan  sebuah  kantong  magma yang besar di bawah merapi bagian utara yang disebut juga MLA (Merapi-Lawu Anomaly).
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
Gambar 8. Intrerpretasi kecepatan di bawah Jawa Tengah pada referensi [3]

## SIMPULAN
Berdasarkan  hasil  penelitian  pencitraan  tomografi  di  bawah  Gunung  api  Merapi  dapat  ditarik kesimpulan:
• Citra  3-D  yang  dihasilkan dari  LOTOS  sudah  memadai  yang  menghasilkan  tiga model  citra anomali,  yaitu  model  citra  anomali  gelombang  P,  citra  anomali  gelombang  S,  dan  model Vp/Vs ratio .  Posisi anomali negatif  dari  ketiga  model  tersebut  berada  pada area  yang  sama, yaitu  anomali  lebih  luas  berada  pada  kedalaman z =  0  km  di  bawah  elevasi  referensi  dan terdalam yaitu pada kedalaman z = 15 km di bawah elevasi referensi.
• Hasil  model  citra  anomali  yang  didapatkan  menjelaskan  bahwa  Gelombang  P  memberikan gambaran  tomogram  yang  lebih  jelas  daripada  gelombang  S,  karena signal  to  noise  ratio gelombang S lebih rendah dan waktu tiba gelombang S tidak sejelas waktu tiba gelombang P. Sifat  fisis  material  pada  bagian  anomali  negatif  di  bawah  Gunung  Merapi  adalah  kurang kompak, panas, heterogen, pada pusat anomali kemungkinan sebagian bersifat partial melting dan  bersifat  fluida.  Posisi  anomali  negatif  tepat  berada  di  bawah  pusat-pusat  erupsi sebelumnya  yang  mencerminkan  sistem  kantong  magma  Gunung  Merapi.  Anomali  negatif sangat berkaitan dengan daerah lemah yang dihasilkan oleh kegiatan vulkanik Gunung Merapi yang  selalu  aktif  di  area  kantong  magma.  Karena  itu  perubahan  sifat  fisis  anomali  secara berkala  berdasarkan  studi  tomografi  ini  kemungkinan  dapat  digunakan  untuk  memantau tingkat kegiatan Gunung Merapi.

## REFERENSI
[1] K. Aki and W. H. K. Lee, 1976, “Determination Of Three-Dimensional Velocity Anomalies Under  a  Seismic  Array  Using  First  P  Arrival  Times  from  Local  Earthquake,” Journal  of Geophysical Research, No. 23, Vol. 81, 4381-4399.
[2] BPPTK, 2000, Leaflet: Pesona Merapi , BPPTK, Yogyakarta.
[3] I. Koulakov and S. V. Sobolev, 2006, “A Tomographic Image of Indian Lithosphere Break- off  Beneath  the  Pamir  Hindukush  Region”, Geophysics  Journal  International, Vol.  164,  p. 425-440.
[4] I.  Koulakov  et  al.,  2007,  “Joint  Inversion  of  Active  and  Passive  Seismic  Data  in  Central Java,” Geophysics Journal International, Vol. 170, p. 923-932.
[5] I.  Koulakov,  2009,  “LOTOS  Code  For  Local  Earthquake    Tomographic  Inversion  : Benchmarks  for  Testing  Tomographic  Algorithms,” Bulletin  of  Seismological  Society  of America (BSSA), No. 1, Vol.99, p. 194-214.
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA
[6] W.  H.  K.  Lee  and  J.  C.  Lahr,  1972,  “Hypo71:  A  computer  Program  for  Determining Hypocenter, Magnitude, and First Motion Pattern of Local Earthquakes,” Open-File Report U. S. Geological Survey , p. 100. California.
[7] J.  M.  Londono,  2002,  “A  Seismic  Model  for  the  Volcanic  Activity  of  Nevado  del  Ruiz Vulcano Colombia,” Thesis, Kyoto University, Kyoto.
[8] C. C. Paige and M. A. Saunders, 1982. “LSQR: An algorithm for sparse linear equations and sparse least squares,” ACM trans. Math. Soft ., 8, 43–7
[9] L.  K.  Steck,  et  al.,  1998.  “Crust  and  Upper  Mantle  P  Wave  Velocity  Structure  Beneath Valles Caldera, New Mexico: Result from the James Teleseismic Tomography Experiment,” Journal Geophys. Res., 103, 24.301-24.320.
[10] J.  Um  and  C.  Thurber,  1987,  “A  Fast  Algorithm  for  Two-Point  Seismic  Ray  Tracing,” Bulletin of The Seismological Society of America , Vol. 77, No. 3, p. 972-966.
Spektra: Jurnal Fisika dan Aplikasinya http://doi.org/10.21009/SPEKTRA