eISSN 2722-9886   pISSN 1411-6669 https://jurnal.unej.ac.id/index.php/MIMS/index
Variasi spasial dan temporal nilai- b pada gempa bumi di wilayah Sulawesi Tengah, Gorontalo, dan sekitarnya menggunakan metode robust fitting
(Spatial and temporal variation of b -value of an earthquake in the Central Sulawesi, Gorontalo, and surrounding areas using robust fitting method)

## Nina Fitriyati*, Madona Yunita Wijaya, M. Alvi Bisyri
Program Studi Matematika, Fakultas Sains dan Teknologi Universitas Islam Negeri Syarif Hidayatullah Jakarta
*korespondensi: nina.fitriyati@uinjkt.ac.id
Received: 01-09-2022, accepted: 27-09-2022

## Abstract
This study discusses variation in seismic and tectonic modeled by a Gutenberg-Richter relationship for earthquakes in the Central Sulawesi, Gorontalo, and surrounding areas using the Robust Fitting Method (RFM) with the weight function of Tukey’s bisquare. The declustering process on earthquake data is carried out using the Reasenberg equation. The  values  for  both  parameters  are  analyzed  spatially  and  temporally.  In  the  spatial analysis, the research area is divided into 43 grids. In the temporal analysis, the research area is divided into zone A and zone B. The data grouping is done using a sliding time window  method,  i.e.,  grouping  50  earthquake  catalogs  with  5  overlapping  events.  The results according to spatial analysis show that the b-values range from 0.38 – 1.19. Areas with low b-values (0.38 – 0.7) occur around the Palu-Koro Fault, i.e., Palu city, Malacca strait, and to Toli-Toli, and also in the northern region of Gorontalo, i.e., the subduction plate of the Sulawesi Sea. Meanwhile, high b-values (0.71 – 1.19) are in the Tomini Bay area which is an area with frequent occurrence of earthquakes but has the small potential to  generate  large-scale  earthquakes.  The  results  of  the  temporal  b-value  estimation  in zones  A  and  B  range  between  values  of  0.38  -  1.25.  The  b-values  appear  to  decrease before the occurrence of major earthquakes in 1996 and 2018 in zone A. The b-values decreased before the occurrence of major earthquakes in 1990, 1991, 2000, and 2008 in zone B. However, the b-values cannot be used as a precursor before the big earthquake in 1997.
Keywords: Tukey’s  bisquare,  Reasenberg  equation,  Gutenberg-Richter  relationship, sliding time window, Robust Fitting Method .
MSC2020: 86A15

## 1. Pendahuluan
Indonesia  merupakan  salah  satu  negara  dengan  tingkat  kegempaan  yang  cukup  tinggi. Hal ini terjadi karena Indonesia berada di daerah pertemuan tiga lempeng tektonik yaitu lempeng  Eurasia,  Samudra  Pasifik  dan  Indo-Australia.  Hal  ini  menyebabkan  banyak pulau  di  Indonesia  rawan  akan  terjadinya  gempa  bumi,  salah  satunya  adalah  Pulau Sulawesi. Berdasarkan buku peta sumber gempa bumi dan bahayanya di Indonesia yang disusun  oleh  Pusat  Studi  Gempa  bumi  Nasional  (PuSGen)  [ 1 ],  seismisitas  di  Pulau Sulawesi  berkaitan  dengan  sejumlah  sesar  dan  subduksi  di  Sulawesi  Tengah  dan Gorontalo yaitu Sesar Palu Koro, Sesar Gorontalo, Sesar Matano, dan Subduksi Lempeng laut Sulawesi dan Subduksi Sangihe. Peta seismisitas Sulawesi untuk magnitudo gempa ≥ 4,5 SR dan struktur geologi di Pulau Sulawesi dapat dilihat pada Gambar 1. Dari gambar ini dapat diamati bahwa gempa-gempa dangkal di darat secara umum terjadi di sekitar Sesar Palu-Koro dan Sesar Matano serta terjadi juga gempa dengan magnitudo ≥ 6 SR di kedua sesar tersebut. Sementara gempa-gempa yang terjadi di Teluk Tomini dipengaruhi oleh aktivitas subduksi Lempeng Laut Sulawesi dan subduksi Sangihe.
(a) (b)
Gambar 1. (a) Peta seismisitas dari katalog PuSGeN 2017 untuk magnitude gempa ≥ 4,5 SR dan (b) struktur geologi Pulau Sulawesi [ 1 ].
Mengingat tingkat kegempaan yang tinggi di Pulau Sulawesi, diperlukan suatu metode untuk  menggambarkan  kondisi  tektonik  dan  potensi  gempa  bumi  besar  di  wilayah tersebut. Prediksi tingkat kegempaan sangat bermanfaat untuk mengantisipasi kerusakan yang diakibatkan oleh gempa tersebut. Salah satu antisipasi yang dapat dilakukan adalah dengan  mengasuransikan  bangunan  dan  kendaraan  sehingga  kerugian  dapat diminimalisir [ 2 ].
Potensi  terjadinya  gempa  bumi  di  suatu  wilayah  dapat  diketahui  dengan  mempelajari tentang kondisi tektonik dan aktivitas kegempaan melalui pendekatan hubungan antara frekuensi  kejadian  dan  magnitudo  gempa  bumi.  Hubungan  tersebut  masih  didominasi oleh relasi Gutenberg-Richter [ 3 ] . Da Silva [ 3 ] mengeneralisasi relasi Gutenberg-Richter berdasarkan fungsi probabilitas. Relasi Gutenberg-Richter merupakan hubungan antara frekuensi  kejadian  gempa  bumi  dengan  Magnitudo,  yang  dinyatakan  dalam  bentuk berikut [ 4 ]:
log 𝑁(𝑀) = 𝑎 − 𝑏𝑀, (1)
dengan 𝑀 adalah magnitudo gempa bumi, 𝑁(𝑀) adalah jumlah kumulatif gempa bumi dengan magnitudo ≥ 𝑀 , nilai 𝑎 dan 𝑏 adalah parameter. Nilai 𝑎 adalah parameter seismik pada suatu wilayah yang menunjukkan tingkat aktivitas kegempaan di suatu wilayah yang besarnya bergantung pada periode observasi dan luasnya wilayah [ 5 ]. Sedangkan nilai- 𝑏 adalah  parameter  tektonik  yang  menggambarkan  kondisi  tektonik  pada  suatu  wilayah yang  berhubungan  dengan  aktivitas  kegempaan.  Nilai- 𝑏 disebut  juga  dengan  koefisien parameter 𝑀 . Nilai- 𝑏 yang rendah pada suatu wilayah mengindikasikan adanya potensi gempa  besar.  Sedangkan  nilai- 𝑏 yang  tinggi  mengindikasikan  potensi  kecil  terjadinya gempa bumi besar [ 5 ].
Pencarian  nilai- 𝑏 terhadap  ruang  yang  menjadi  parameter  seismotektonik  di  suatu wilayah  disebut  nilai- 𝑏 secara  spasial  dan  dinotasikan  dengan  nilai- 𝑏 (s).  Nilai- 𝑏 (s) menjelaskan  distribusi  frekuensi  magnitudo  yang  terjadi  dalam  data  seismik,  sehingga variasi  nilai- 𝑏 (s)  menjadi  salah  satu  komponen  kunci  dari  peta  bahaya  seismic  [ 6 ]. Nilai- 𝑏 yang  tinggi  berkaitan  intensitas  stres  yang  rendah  atau  heterogenitas  material yang tinggi [ 7 ]. Perubahan nilai- 𝑏 secara temporal baik dalam jangka pendek, menengah dan panjang dapat menjadi suatu prekursor (tanda-tanda awal terjadinya gempa bumi) . Hasil penelitian [ 8 ] menunjukkan bahwa gempa bumi besar didahului oleh peningkatan nilai- 𝑏 dalam  jangka  waktu  menengah  kemudian  diikuti  oleh  penurunan  nilai- 𝑏 dalam beberapa bulan sebelum terjadinya gempa bumi besar. Penelitian mengenai variasi spasial dan temporal nilai- 𝑏 telah dilakukan untuk  meneliti pola kegempaan dan potensi gempa di wilayah Pulau Jawa [ 9 ], [ 10 ], wilayah Papua [ 5 ], wilayah Taiwan [ 11 ].
Beberapa metode estimasi nilai- 𝑏 telah dikembangkan. Salah satunya adalah metode least square . Namun metode ini menghasilkan estimasi nilai- 𝑏 yang bias karena hasil estimasi sangat  rentan  terhadap  data  pencilan  [ 12 ].  Aki  [ 13 ]  mengusulkan  penggunaan  metode Maksimum  Likelihood  Estimation (MLE)  dan  Maslov  dan  Anokhin  [ 14 ]  mengusulkan penggunaan Generalized  Logistic  equation (GLE).  Menurut  Han  dkk.  [ 15 ],  MLE memberi  terlalu  banyak  bobot  untuk  kejadian  gempa  yang  kecil  dan  sensitivitasnya terhadap amplitudo kelengkapan Magnitude of Completeness (Mc) sangat luas. Sehingga nilai- 𝑏 hasil  MLE  dapat  memiliki  penyimpangan  yang  besar  jika  Mc  tidak  ditentukan secara akurat. Oleh karena itu, Han dkk. [ 15 ] mengusulkan metode Robust Fitting Method (RFM) . RFM memberikan nilai- 𝑏 yang lebih stabil dan memiliki sensitivitas yang baik
terhadap  terjadinya  gempa  bumi  besar.  Hasil  penelitian  Han  dkk.  [ 15 ]  menunjukkan bahwa  variasi  nilai- 𝑏 menggunakan  RFM  lebih  efektif  sebagai prekursor gempa  bumi dibandingkan MLE.
RFM  yang  dikembangkan  oleh  Han  dkk.  [ 15 ]  adalah  metode  estimasi  yang  tangguh ( robust )  untuk  parameter 𝑀 yang  disebut Iteratively  Reweight  Least  Square  (IRLS) dengan fungsi pembobot Bisquare Tukey . Metode tersebut telah digunakan di beberapa penelitian  salah  satunya  Bai  dkk.  [ 16 ].  Bai  dkk.  [ 16 ]  merekomendasikan  fungsi  bobot Bisquare Tukey karena kinerja yang baik.
Pada  penelitian  ini  akan  dikaji  estimasi  parameter  tektonik  melalui  nilai- 𝑏 pada  relasi Gutenberg-Richter untuk  Sulawesi  Tengah,  Gorontalo,  dan  sekitarnya  menggunakan Metode RFM dengan bobot Bisquare Tukey . Estimasi parameter dilakukan secara spasial maupun temporal. Variasi nilai- 𝑏 secara spasial digunakan untuk mengetahui nilai- 𝑏 di beberapa  wilayah  Sulawesi  Tengah  dan  Gorontalo  yang  berpotensi  terjadinya  gempa besar.  Sedangkan  variasi  nilai- 𝑏 secara  temporal  digunakan  untuk  mengetahui  naik turunnya  nilai- 𝑏 terhadap  waktu  sebagai prekursor terjadinya  gempa  bumi.  Data  yang digunakan  pada  penelitian  ini  adalah  data  katalog  gempa  bumi  Sulawesi  Tengah  dan Gorontalo tahun 1985 – 2020 dengan batas wilayah 3,8° LS – 1,8° LU dan 119° BT – 123,65° BT yang diambil dari National Earthquake Information Center (NEIC) USGS https://earthquake.usgs.gov/earthquakes [ 17 ]. Data berjumlah 3078 kejadian gempa bumi yang  berupa titik terjadinya gempa bumi ( longitude dan latitude ), waktu kejadian gempa bumi  (tahun,  bulan,  hari,  jam,  menit),  kedalaman  gempa  bumi  dan  magnitudo  gempa bumi.

## 2. Metodologi
Gempa bumi yang terjadi di Sulawesi tengah dan Gorontalo diakibatkan oleh pergerakan tektonik  yang  berbeda-beda.  Seperti  gempa  besar  di  Gorontalo  pada  tahun  2008 disebabkan  oleh  aktivitas  subduksi  lempeng  laut  Sulawesi.  Kemudian  gempa  besar  di Palu  pada  tahun  2018  disebabkan  oleh  aktivitas  sesar  Palu  Koro.  Untuk  menganalisis nilai- 𝑏 , data perlu dipisahkan menjadi gempa utama dan gempa susulan melalui algoritma declustering .  Algoritma  ini  digunakan  untuk  mengindetifikasi  gempa  susulan  dan mengelompokkannya  sesuai  dengan  zona  interaksi  spasial  dan  temporal.  Salah  satu algoritma declustering adalah melalui persamaan Reasenberg berikut [ 18 ]:
𝜏 =   − ln(1 − 𝑝 1 )𝑡/10 2(∆𝑀−1)/3 , (2)
dengan 𝜏 adalah nilai look-ahead time untuk membangun klaster ketika kejadian gempa bumi pertama tidak berklaster, 𝑝 1 adalah peluang terjadinya gempa utama, 𝑥 𝑚𝑒𝑓𝑓 adalah magnitudo cut-off dari katalog gempa, 𝑡 adalah waktu, dan ∆𝑀 = 𝑀 gempa utama − 𝑥 𝑚𝑒𝑓𝑓 . Berikut  adalah  nilai  parameter  standar  dari  algoritma declustering Reasenberg
berdasarkan  pusat  pengujian  RELM  [ 19 ].  Algoritma declustering versi  terakhir  dapat diunduh  dari website USGS  dengan  nama  CLUSTER2000.  Beberapa  nilai  input parameter pada algoritma tersebut dapat dilihat pada Tabel 1 [ 20 ].
Berikut adalah penjelasan beberapa istilah pada Tabel 1: 𝜏 𝑚𝑖𝑛 dan 𝜏 𝑚𝑎𝑥 , masing-masing, adalah  nilai  minimum  atau  nilai  maksimum  dari look-ahead  time untuk  membangun klaster ketika kejadian gempa bumi pertama tidak berklaster, 𝑝 1 adalah probabilitas untuk mendeteksi kejadian pada klaster berikutnya digunakan untuk menghitung 𝜏 , 𝑥 𝑘 adalah peningkatan  magnitudo cut-off yang  lebih  rendah  selama  pengclusteran  : 𝑥 𝑚𝑒𝑓𝑓 = 𝑥 𝑚𝑒𝑓𝑓 + 𝑥 𝑘 𝑀 , 𝑟 𝑓𝑎𝑐𝑡 adalah jumlah jari-jari retakan yang mengelilingi setiap gempa bumi dalam  peristiwa  baru  yang  dianggap  sebagai  bagian  dari  klaster.  Setelah  proses declustering , langkah selanjutnya adalah menentukan nilai Mc.
Tabel 1. Nilai parameter standar pada algoritma declustering Reasenberg [ 20 ].
Parameter  Nilai Standard  Nilai Minimum  Nilai Maksimum 𝜏 𝑚𝑖𝑛 (hari) 1 0,5 2,5 𝜏 𝑚𝑎𝑥 (hari) 10 3 15 𝑝 1 0,95 0,9 0,99 𝑥 𝑘 0,5 0 1 𝑥 𝑚𝑒𝑓𝑓 1,5 1,6 1,8 𝑟 𝑓𝑎𝑐𝑡 10 5 20
Nilai Mc didefinisikan sebagai magnitudo terendah dari 100% gempa bumi dalam volume ruang-waktu  yang terdeteksi [ 21 ]. Nilai Mc dapat berbeda-beda di setiap wilayah yang diteliti. Penentuan nilai Mc sangat berpengaruh pada nilai- 𝑏 karena nilai Mc berpengaruh pada jumlah sampel yang dapat diambil. Jika nilai Mc diambil terlalu tinggi maka jumlah sampel akan sangat sedikit karena akan banyak data yang dibuang. Namun jika nilai Mc diambil terlalu kecil maka hasil estimasi nilai parameter seismisitas akan menjadi bias. Pada  penelitian  ini,  nilai  Mc  ditentukan  menggunakan  metode maximum  curvature (MAXC). Metode ini merupakan cara cepat dan mudah dalam memperkirakan nilai Mc, yaitu  dengan  menentukan  titik  kelengkungan  maksimum  dengan  menghitung  nilai maksimum  turunan  pertama  dari  kurva  frekuensi  magnitude  [ 22 ].  Meskipun  mudah dalam penerapannya, namun menurut penelitian Mignan dkk. [ 21 ], teknik MAXC valid dalam  pemetaan  Mc  resolusi  tinggi.  Untuk  mencapai  hasil  yang  stabil,  teknik  MAXC memerlukan sampel kejadian yang lebih sedikit dibandingkan teknik lain.
Untuk menganalisis nilai- 𝑏 secara spasial dan temporal, nilai- 𝑏 dihitung menggunakan RFM.  Metode  ini  diusulkan  oleh  Han  dkk.  [ 15 ].  RFM  dapat  meminimalkan  pengaruh pencilan  seperti  yang  terjadi  pada least  square ,  memberikan  nilai- 𝑏 yang  stabil  dan handal tanpa persyaratan ketat untuk Mc. Selain itu RFM juga dapat memberikan variasi nilai- 𝑏 sebagai prekursor gempa besar dengan baik dibanding MLE. Pada penelitian ini,
algoritma  RFM  yang  digunakan  untuk  mengestimasi  nilai- 𝑏 adalah Iteratively Reweighted  Least  Square (IRLS)  dengan  fungsi  pembobot Bisquare  Tukey sebagai berikut [ 16 ]:
𝑤 = { (1 − ( 𝑢 𝑖 𝑘 ) 2 ) 2 |𝑢 𝑖 | < 𝑘 0                       |𝑢 𝑖 | ≥ 𝑘 , (3)
dengan 𝑢 𝑖 adalah  residual  ke- 𝑖 , 𝑘 adalah  konstanta  tuning  yang  digunakan  untuk menentukan  kekuatan  penaksiran  terhadap  pencilan  dan  efisiensi  penaksiran  dalam ketidakadaan pencilan [ 23 ], [ 24 ]. Fungsi pembobot ini memiliki kinerja yang baik dalam kasus outlier lavarage yang rendah dan tinggi [ 16 ].
Berikut  adalah  prosedur  RFM  dengan  estimasi  nilai- 𝑏 menggunakan  IRLS  dan  fungsi pembobot Bisquare Tukey [ 15 ]:
1) Memfit model (1) menggunakan weighted least square (WLS).
2) Menghitung r adj ( adjusted residual ) menggunakan persamaan berikut:
𝑟 𝑎𝑑𝑗 = 𝑟 𝑖 √1 − ℎ 𝑖 ,
dengan 𝑟 𝑖 adalah residual  least  squares ke- 𝑖 dan ℎ 𝑖 adalah  leverages  yang menyesuaikan residual dengan menurunkan bobot poin pada data high-leverage yang memiliki  efek  besar  pada  model  fit  menggunakan  WLS. Adjusted  residual ( 𝑟 𝑎𝑑𝑗 ) terstandarisasi dihitung menggunakan persamaan:
𝑢 = 𝑟 𝑎𝑑𝑗 𝑘.𝑠 ,
dengan 𝑘 adalah  konstanta  tuning dan 𝑠 adalah  variansi  robust  yang  dihitung menggunakan persamaan:
𝑠 = 𝑀𝐴𝐷 0.6745 = 𝑚𝑒𝑑𝑖𝑎𝑛{𝑟 𝑖 − 𝑚𝑒𝑑𝑖𝑎𝑛(𝑟 𝑖 )}
0.6745 ,
dengan  MAD  adalah Median  absolute  deviation .  Konstanta  0,6745  dipilih sedemikain sehingga estimator nilai 𝑠 menjadi unbiased dari 𝜎 jika jumlah sampel besar dan error berdistribusi normal [ 25 ].
3) Menghitung robust weights menggunakan persamaan (3) dengan 𝑘 = 1 .
4) Jika  model  fit  konvergen,  maka  iterasi  selesai,  jika  tidak  maka  lakukan  iterasi berikutnya dengan kembali ke langkah no 1).
Pada  analisis  nilai- 𝑏 secara  spasial,  disebut  nilai- 𝑏 (s),  wilayah  Sulawesi  Tengah, Gorontalo,  dan  sekitarnya  akan  dibagi  menjadi  beberapa  grid  dengan  ukuran  0.49°  × 0.49°. Peta spasial nilai- 𝑏 dilakukan dengan bantuan Software ArcMap . Sedangkan untuk analisis  nilai- 𝑏 secara  temporal,  disebut  nilai- 𝑏 (t),  jangka  waktu  penelitian  dilakukan
pada  tahun  1985  sampai  dengan  2020.  Perhitungan  nilai- 𝑏 (t)  menggunakan  metode sliding  time  window yaitu  nilai- 𝑏 (t)  dihitung  untuk 𝑁 Gempa  bumi  seperti  yang digunakan dalam penelitian [ 5 ], [ 15 ], [ 26 ]. Kemudian time window digeser dengan jumlah gempa bumi tetap dan seterusnya sampai gempa bumi terakhir. Penggeseran window ini disebut overlap . Pada penelitian ini kami memilih 𝑁 = 50 dan overlap 5 gempa bumi.
Untuk analisis nilai- 𝑏 (t), wilayah penelitian akan dibagi menjadi dua zona, yaitu zona A dan zona B (perhatikan Gambar 9). Pembagian zona ini untuk mengetahui  nilai- 𝑏 (t) di zona  A  yang  berada  di  wilayah  sesar  Palu  Koro  dan  sekitarnya.  Zona  ini  berada  pada koordinat dari tiga persegi, yaitu persegi satu dengan titik 3,0° LS – 1,8° LU & 119,22° BT – 120,6° BT, persegi dua dengan titik 1,84° LS – 0,47° LS dan 120,6° BT – 121,18° BT, dan persegi tiga dengan titik 3,0° LS – 1,84° LS dan 120,6° BT – 123,65° BT. Zona B yang berada di wilayah sekitar subduksi laut Sulawesi dan subduksi Sangihe. Zona B terletak pada koordinat dari dua persegi yaitu persegi satu dengan titik 0,47° LS  – 1,8° LU dan 120,6° BT – 123,65° BT, dan persegi 2 dengan titik 1,84° LS  – 0,47° LS dan 121,18° BT – 123,65° BT.

## 3. Hasil dan Pembahasan
Data yang digunakan dalam penelitian ini adalah data katalog gempa bumi yang diambil dari NEIC-USGS untuk wilayah Sulawesi Tengah, Gorontalo, dan sekitarnya. Soehaimi dkk.  [ 27 ]  menyatakan  bahwa  gempa  bumi  merusak  sering  terjadi  di  sepanjang  lajur patahan  Palu  Koro  seperti  pada  tahun  1927,  1938,  1985,  1998,  2007  dan  2012.  Oleh karena itu, data gempa bumi yang digunakan pada penelitian ini dipilih data pada rentang waktu 1 Januari 1985 sampai dengan Bulan Mei 2020. Batas daerah ditentukan pada 3,8° LS – 1,8° LU dan 119° BT – 123,65° BT. Terdapat 3078 data kejadian gempa bumi. Peta sebaran kejadian gempa bumi tersebut dapat dilihat pada Gambar 2 dengan kedalaman sumber gempa dikelompokkan menjadi 3 kategori z < 184,7 km, z < 369,4 km, dan z < 615,6 km. Tabel 2 menyajikan empat kejadian gempa bumi terbesar pada rentang waktu penelitian.

## Tabel 2. Empat kejadian gempa bumi terakhir
Longitude Latitude Tahun Bulan Hari M Kedalaman (km) Jam Menit 120,2622 -1,6864 2020 5 20 4,3 69,68 20 50 121,0113 0,1937 2020 5 18 4,4 158,47 19 20 119,6699 -1,2324 2020 5 9 4,2 10 23 43 123,5027 -0,0596 2020 5 1 4,1 67,67 21 46

## Gambar 2. Peta Kejadian Gempa Bumi Sulawesi Tengah dan Gorontalo
Untuk  memisahkan  gempa  bumi  utama  dan  susulan,  dilakukan  proses declustering terhadap data katalog gempa bumi di atas menggunakan Metode Reasenberg. Proses ini menghasilkan 2785 kejadian gempa bumi utama yang ditampilkan pada Gambar 3. Hal ini berarti bahwa terdapat sekitar 10% kejadian gempa susulan. Gempa terbanyak berada di  sekitar  wilayah  Gorontalo  yang  merupakan  jalur  subduksi  Lempeng  Sangihe. Sedangkan gempa besar terjadi di wilayah Laut Sulawesi dan Selat Malaka.
Gambar 3. Peta kejadian gempa bumi di Sulawesi Tengah & Gorontalo setelah proses declustering .

## 3.1. Analisis Spasial Nilai- 𝒃
Pada analisis nilai- 𝑏 (s), wilayah di  Sulawesi  Tengah dan Gorontalo  dibagi  menjadi  43 grid  dengan  ukuran  0,49°  ×  0,49°.  Analisis  nilai- 𝑏 (s)  terhadap  data  gempa  hasil declustering dilakukan secara keseluruhan dan secara parsial sesuai dengan banyaknya grid. Penentuan nilai Mc dilakukan menggunakan metode MAXC. Distribusi frekuensi magnitudo pada seluruh data gempa bumi disajikan pada Gambar 4.
Gambar 4. Distribusi frekuensi magnitudo pada keseluruhan data hasil declustering .
Mengikuti  [5],  untuk  memenuhi  linearitas  terhadap  distribusi  frekuensi  magnitudo, estimasi  nilai  Mc  pada  seluruh  wilayah  Sulawesi  Tengah,  Gorontalo  dan  sekitarnya adalah  4,4.  Hal  ini  berarti  katalog  USGS  mampu  merekam  gempa  dengan  magnitudo terkecil 4,4 SR dengan baik. Berdasarkan distribusi frekuensi magnitudo diperoleh nilai- b sebesar 0,921 dan nilai- a sebesar 7,33. Hasil estimasi nilai- a (s) dan nilai- 𝑏 (s) pada relasi Gutenberg-Richter  dapat  dilihat  pada  Tabel  3.  Nilai  Std. Error  adalah nilai  perkiraan deviasi standar dari populasi sampel statistik sedangkan nilai- t adalah nilai yang tertera pada  tabel  distribusi t - student yang  digunakan  untuk  pengujian  signifikansi  dari parameter a dan b .
Tabel 3. Hasil estimasi koefisien relasi Gutenberg-Richter menggunakan RFM untuk
keseluruhan data hasil declustering . Parameter Nilai Std. Error Nilai- t a (s) 7,3305 0,3510 18,6220 b (s) 0,9210 0,0571 -14,3626
Dari hasil estimasi ini maka dapat diketahui relasi Gutenberg-Richter untuk keseluruhan data gempa bumi di wilayah Sulawesi Tengah, Gorontalo dan sekitarnya adalah:
log 𝑁(𝑀) = 7,3305 − 0,9210𝑀.
Scatterplot dari relasi Gutenberg-Richter di atas dapat dilihat pada Gambar 5.
Gambar 5. Scatterplot relasi Gutenberg-Richter menggunakan RFM
Menggunakan  cara  yang  serupa  dengan  di  atas,  akan  diperoleh  estimasi  nilai  Mc  dan estimasi koefisien pada relasi Gutenberg-Richter untuk setiap grid. Tabel 4 menyajikan hasil  estimasi  tersebut  pada  beberapa  grid.  Gambar  6  merupakan  peta  spasial  untuk estimasi nilai- 𝑏 (s) pada 43 grid. Latitude dan longitude dari titik tengah grid-grid tersebut tertera pada Tabel 4.
Tabel 4. Hasil estimasi nilai Mc dan koefisien relasi Gutenberg-Richter menggunakan RFM pada beberapa grid.
No Grid Titik Tengah Grid Banyak Data Mc nilai- 𝒃 (s) nilai- a (s) Relasi Gutenberg-Richter Longitude Latitude 1 122,6 -2,95 31 4,4 0,7104 3,8605 log 𝑁(𝑀) = 3,8605 – 0,7104 𝑀 2 121,7 -2,5 31 4,3 0,6981 4,0906 log 𝑁(𝑀) = 4,0906 – 0,6981 𝑀 3 120,8 -2,275 31 4,2 1,0585 5,5857 log 𝑁(𝑀) = 5,5857 – 1,0585 𝑀 4 120,35 -1,825 33 4,2 0,81 4,4815 log 𝑁(𝑀) = 4,4815 – 0,8100 𝑀 5 119,45 -1,375 35 4,2 0,7243 4,2937 log 𝑁(𝑀) = 4,2937 – 0,7243 𝑀 6 120,35 -1,375 43 4,5 0,6322 3,9084 log 𝑁(𝑀) = 3,9084 – 0,6322 𝑀 7 120,575 -1,375 46 4,6 0,6004 4,052 log 𝑁(𝑀) = 4,0520 – 0,6004 𝑀 8 121,25 -1,15 40 4,1 0,59 3,5948 log 𝑁(𝑀) = 3,5948 – 0,59 00𝑀 9 120,575 -0,925 34 4,2 0,6285 3,7306 log 𝑁(𝑀) = 3,7306 – 0,6285 𝑀 10 122,6 -1,6 31 4,3 0,4936 3,2139 log 𝑁(𝑀) = 3,2139 – 0,4936 𝑀 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 42 121,25 2 48 4,3 0,9762 5,447 log 𝑁(𝑀) = 5,447 – 0,9762 𝑀 43 122,6 1,775 36 4,3 0,825 4,8118 log 𝑁(𝑀) = 4,8118 – 0,825 𝑀
Gambar 6. Peta spasial untuk estimasi nilai- 𝑏 (s) pada 43 grid
Pada Gambar 6 nilai- 𝑏 (s) pada 43 grid berkisar antara 0,38 – 1,19 dan diklasifikasikan menjadi  6 klasifikasi dengan warna  yang berbeda. Nilai- 𝑏 (s)  yang rendah  (0,38  – 0,7) didapatkan di hampir seluruh wilayah Sulawesi tengah mulai dari Palu, Donggala, Teluk Palu,  Selat  Malaka  hingga  Toli-Toli  dan  wilayah  Teluk  Tolo.  Kemudian  di  wilayah Gorontalo nilai- 𝑏 (s) yang rendah terdapat pada laut Sulawesi di wilayah utara Gorontalo. Sedangkan  nilai- 𝑏 (s)  yang  tinggi  (0,71  –  1,19)  berada  di  wilayah  teluk  Tomini  dekat dengan  Kep.  Togian  dan  selatan  Gorontalo  di  sekitar  Kabupaten  Bone  Bolango. Berdasarkan hasil tersebut, diperoleh nilai- 𝑏 (s) yang rendah berada pada wilayah Sesar Palu Koro dan wilayah laut Sulawesi yang merupakan subduksi Lempeng Laut Sulawesi Utara.  Kedua  wilayah  tersebut  memiliki  riwayat  gempa  bumi  besar  seperti  pada  tahun 1996,  2008,  dan  2018.  Hal  ini  sejalan  dengan  hasil  penelitian  Rohadi  [ 5 ]  yang mengindikasikan bahwa nilai- 𝑏 (s) yang rendah memiliki potensi terjadinya gempa besar, sedangkan nilai- 𝑏 (s) yang tinggi berpotensi kecil terjadi gempa besar. Peta nilai- 𝑏 (s) pada Gambar 6 hanya untuk menggambarkan kondisi tektonik dan potensi gempa bumi besar di  wilayah  tersebut  dan  tidak  dapat  digunakan  untuk  memprediksi  terjadinya  gempa. Selain  itu,  peta  ini  tidak  dapat  diprediksi  sampai  berapa  tahun  estimasi  nilai- 𝑏 (s)  ini berlaku.
Selanjutnya  akan  diestimasi  nilai- 𝑏 (s)  untuk  gempa  yang  terjadi  pada  periode  waktu tanggal 1 Juni 2020 – 1 Agustus 2020. Terdapat 9 gempa bumi yang terjadi pada periode waktu tersebut dapat dilihat pada Tabel 5. Bila kita petakan secara spasial, nilai-b(s) untuk periode  waktu  tersebut  dapat  dilihat  pada  Gambar  7.  Titik  gempa  pada  periode  waktu
tersebut disimbolkan dengan titik bulat merah. Dari titik gempa bumi yang terjadi hampir semua terjadi pada semua wilayah dengan nilai- 𝑏 (s) tinggi (0,71 – 1,19). Hanya ada dua gempa bumi yang terletak di wilayah dengan nilai- 𝑏 (s) rendah yaitu pada tanggal 10 Juni 2020 dan 14 Juli 2020.
Tabel 5. Gempa bumi yang terjadi pada 1 Juni 2020 – 1 Agustus 2020.
Longitude Latitude Tahun Bulan Hari M Kedalaman (km) Jam Menit 122,9192 -0,0798 2020 8 1 4,5 149,92 7 56 123,6477 0,5715 2020 7 29 4,2 243,41 0 28 123,0804 -2,0524 2020 7 14 4,6 10 2 59 120,7944 0,2119 2020 7 11 4,8 109,25 7 17 123,5902 0,5232 2020 6 27 4,3 246,42 12 4 121,7377 0,775 2020 6 21 4,2 77,32 7 44 122,2145 -2,7093 2020 6 17 4,9 10 11 10 122,2814 -2,7948 2020 6 15 4,9 10 23 36 122,426 -0,685 2020 6 10 4,9 14,58 20 1
Beberapa gempa bumi untuk periode tanggal 1 Juni 2020 – 1 Agustus 2020 terjadi pada wilayah  dengan  ativitas  kegempaan  tinggi.  Hal  ini  ditunjukkan  dengan  nilai- a (s)  yang tinggi di wilayah tersebut seperti yang ditunjukkan pada Tabel 5. Nilai- a (s) merupakan parameter seismik yang menggambarkan aktivitas kegempaan di wilayah tersebut.
Gambar 7. Peta spasial untuk estimasi nilai- 𝑏 (s) untuk periode tanggal 1 Juni–1 Agustus 2020 .

## 3.2. Analisis Temporal Nilai- b
Pada  analisis  temporal,  nilai- 𝑏 (t)  dihitung  menggunakan  metode sliding  time  window yaitu  nilai- 𝑏 (t)  dihitung  berdasarkan  data  gempa  dari  katalog  gempa  yang  berisi  50 kejadian gempa bumi dan digeser setiap 5 kejadian. Data katalog gempa yang digunakan adalah  3078  kejadian  gempa  bumi.  Pada  Tabel  6  disajikan  hasil  estimasi  nilai- 𝑏 (t) menggunakan RFM pada time window pertama. Dari estimasi parameter ini dapat dibuat persamaan log 𝑁(𝑀) = 4,6058 − 0,7148𝑀 .  Dengan  menggunakan  metode  yang  sama, estimasi  nilai- 𝑏 (t)  untuk  seluruh  data  katalog  gempa  tahun  1985  –  2020  di  Sulawesi Tengah  dan  Gorontalo  dapat  dilihat  pada  Gambar  6.  Gambar  ini  memperlihatkan perubahan  nilai- 𝑏 (t)  berdasarkan  waktu  di  mana  nilai- 𝑏 (t)  berkisar  antara  0,3  sampai dengan  1,6.  Kita  akan  membandingkan  estimasi  nilai- 𝑏 (t)  pada  Gambar  6  dengan beberapa gempa besar yang terjadi pada kurun waktu tahun 1985 – 2020 yang disajikan pada Tabel 7.
Tabel 6. Hasil estimasi nilai-b(t) menggunakan RFM pada time window pertama.
Parameter Nilai Std. Error Nilai- t a (t) 4,6058 0,9194 5,0095 b (t) -0,7148 0,1667 -4,2888
Gambar 8. Variasi temporal nilai- 𝑏 (t) tahun 1985-2020 di Sulawesi Tengah dan Gorontalo.
Tabel 7 menunjukkan semua gempa besar (dengan magnitudo ≥ 7 SR) yang terjadi di wilayah Sulawesi tengah dan Gorontalo. Beberapa gempa besar yang terjadi di wilayah tersebut didahului oleh penurunan nilai- 𝑏 (t) seperti yang ditunjukkan pada Gambar 8. Hal ini sejalan dengan pendapat [ 8 ] bahwa perubahan nilai-b secara temporal dapat menjadi prekursor gempa bumi. Dari Gambar 8 dan Tabel 6 dapat dilihat hanya gempa besar pada tahun 1997 yang kurang jelas prekursornya.

## Tabel 7. Gempa besar dengan magnitudo ≥ 7 SR yang terjadi di Sulawesi Tengah dan Gorontalo
Longitude Latitude Tahun Bulan Hari M Kedalaman (km) Jam  Menit 122,884 1,188 1990 4 18 7,8 25,72 13 39 122,962 1,166 1991 5 19 7 33,06 0 58 122,788 1,196 1991 6 20 7,5 31,35 5 18 119,931 0,729 1996 1 1 7,9 24,04 8 5 120,450 1,001 1996 7 22 7 33 14 19 122,536 1,241 1997 11 25 7 24,02 12 14 123,571 -1,108 2000 5 4 7,6 26,11 4 21 122,090 1,271 2008 11 16 7,4 30,07 17 2 119,848 -0,284 2018 9 28 7,5 19,03 10 2
Untuk mengetahui perubahan nilai- b (t) di wilayah Sulawesi Tengah dan Gorontalo secara detail,  kami  membagi  wilayah  penelitian  menjadi  yaitu  Zona  A  dan  Zona  B  seperti tampak pada Gambar 9. Pembagian wilayah ini bertujuan untuk mengetahui perubahan nilai- 𝑏 (t) di wilayah Sulawesi Tengah yang berada di sekitar Sesar Palu Koro dan Sesar Matano di Zona A, serta wilayah Gorontalo yang berada di sekitar subduksi lempeng laut Sulawesi dan subduksi Sangihe pada Zona B. Pada zona A terdapat 3 gempa bumi besar dan mayoritas gempa buminya diakibatkan oleh aktivitas sesar Palu Koro. Pada zona B terdapat 6 gempa bumi besar dan mayoritas gempa buminya diakibatkan oleh aktivitas subduksi lempeng laut Sulawesi dan subduksi Sangihe.
Gambar 9. Pembagian zona A dan B
Sebelum mengestimasi nilai- 𝑏 (t) untuk zona A dan B, terlebih dahulu ditentukan nilai Mc pada kedua zona. Gambar 10 menyajikan grafik perubahan nilai Mc berdasarkan waktu untuk  zona  A.  Dapat  dilihat  bahwa  nilai  Mc  berubah-ubah  sesuai  dengan  besarnya magnitudo gempa pada waktu tersebut. Gambar 11 menyajikan grafik perubahan nilai- 𝑏 (t) untuk zona A. Nilai- 𝑏 (t) pada Gambar 11 dihitung menggunakan RFM di setiap time window . Berdasarkan Gambar 11, variasi temporal nilai- 𝑏 (t) pada zona A berkisar antara 0,38 - 1,25. Nilai- 𝑏 (t) terlihat menurun sebelum terjadinya gempa bumi besar pada tahun 1996 dan 2018.
Gambar 10. Grafik perubahan nilai Mc berdasarkan waktu untuk zona A
Gambar 11. Grafik variasi nilai- 𝑏 ( t ) berdasarkan waktu untuk zona A
Gambar 12 menyajikan grafik perubahan nilai Mc berdasarkan waktu untuk zona B dan Gambar  13  menyajikan  grafik  perubahan  nilai- b ( t )  untuk  zona  B.  Nilai- 𝑏 ( t )  dihitung menggunakan  RFM  di  setiap time  window .  Berdasarkan  Gambar  13,  variasi  temporal nilai- 𝑏 ( t ) pada zona B berkisar antara 0,36 - 1,36. Selain itu, berdasarkan gambar tersebut nilai- 𝑏 ( t ) menurun sebelum terjadinya gempa bumi besar pada tahun 1990, 1991, 2000, dan 2008. Gempa bumi besar  yang kurang jelas prekursornya terjadi pada tahun 1997, nilai- 𝑏 ( t ) terlihat mengalami peningkatan sebelum terjadinya gempa bumi besar.
Gambar 12. Grafik perubahan nilai Mc berdasarkan waktu untuk zona B
Gambar 13. Grafik variasi nilai- 𝑏 (t) berdasarkan waktu untuk zona B
Dari variasi temporal nilai- 𝑏 (t) dapat diketahui bahwa hampir semua gempa bumi besar yang terjadi pada kedua zona didahului oleh peningkatan nilai- 𝑏 (t) pada jangka menengah dan mengalami penurunan dalam waktu singkat sesuai dengan hasil penelitian Sammonds dkk. [ 8 ]. Hal ini menunjukkan bahwa variasi nilai- 𝑏 (t) di beberapa wilayah gempa bumi mampu menjadi prekursor sebelum terjadinya gempa bumi besar. Oleh karena itu, dapat dilihat juga bahwa RFM menghasilkan sensitivitas yang baik terhadap terjadinya gempa bumi besar di Zona A dan B sejalan dengan hasil penelitian Han dkk. [ 15 ]. Terbukti dari penurunan nilai- 𝑏 (t) yang cepat dalam waktu singkat sebelum gempa bumi besar terjadi dan peningkatan  yang cepat juga setelah gempa bumi besar terjadi seperti yang terjadi pada tahun 2008 di zona B dan tahun 2018 di zona A.

## 4. Kesimpulan
Hasil pemetaan secara spasial, nilai- 𝑏 (s) pada 43 grid berkisar di antara nilai 0,38 – 1,19. Wilayah dengan nilai- 𝑏 (s) rendah (0,38 - 0,7) terdapat di sekitar Sesar Palu Koro, yaitu Kota Palu, Selat Malaka, hingga Toli-Toli dan juga berada di wilayah utara Gorontalo yaitu  laut  Sulawesi  yang  merupakan  subduksi  lempeng  laut  Sulawesi.  Hal  ini mengindikasikan bahwa beberapa wilayah tersebut berpotensi terjadi gempa bumi besar. Sedangkan nilai- 𝑏 (s)  yang cukup tinggi  (0,71  – 1,19) berada di  wilayah Teluk  Tomini yang merupakan wilayah yang sering terjadi gempa bumi namun berpotensi kecil terjadi gempa bumi besar.
Hasil  pemetaan  secara  temporal,  nilai- 𝑏 (t)  pada  zona  A  dan  B  menunjukkan  hampir semua  kejadian  gempa  bumi  besar  didahului  oleh  penurunan  nilai- 𝑏 (t).  Pada  zona  A, nilai- 𝑏 (t) terlihat  menurun sebelum terjadinya  gempa bumi besar pada tahun 1996 dan 2018. Pada zona B, nilai- 𝑏 (t) menurun sebelum terjadinya gempa bumi besar pada tahun 1990,  1991,  2000,  dan  2008.  Namun  pada  tahun  1997,  nilai- 𝑏 (t)  tidak  dapat  dijadikan prekursor sebelum gempa besar terjadi.

## Daftar Pustaka
[1]   M.  Irsyam,  D.  Hilman  dan  I.  Meilano,  Peta  Sumber  Gempa  dan  Bahaya  Gempa Bumi  Indonesia  Tahun  2016,  Jakarta:  Pusat  Studi  Gempa  Nasional  (PuSGeN), 2017.
[2]   R. Cahyandari, A. S. Awalluddin, D. S. Mariani, Sukono dan P. L. Ghazali, “Table of  integration  model  for  motor  vehicle  sharia  insurance,” InPrime:  Indonesian Journal  of  Pure  and  Applied  Mathematics, vol.  2,  no.  1,  pp.  59-64,  2020,  doi: https://doi.org/10.15408/inprime.v2i1.14811
[3]   S.L.E.F. Da Silva, “κ-generalised gutenberg–richter law and the self-similarity of earthquakes,” Chaos  Solitons  &  Fractals, p.  143(C):110622,  2021,  doi: https://doi.org/10.1016/j.chaos.2020.110622
[4]   B. Gutenberg dan C. Richter, “Frequency of earthquakes in California,” Bulletin of the  Seismological  Society  of  America, vol.  34,  no.  4,  pp.  185-188,  1944,  doi: https://doi.org/10.1785/BSSA0340040185
[5]   S.  Rohadi,  “Distribusi  spasial  dan  temporal  parameter  seismotektonik  sebagai indikasi tingkat aktivitas kegempaan di wilayah Papua,” Jurnal  Meteorologi  dan Geofisika, vol.
16, no. 3, pp. 189-198, 2015, doi:
http://dx.doi.org/10.31172/jmg.v16i3.289
[6]   J. Wesseloo, “Evaluation of the spatial variation of b -value,” The Journal of The Southern African Institute of Mining and Metallurg, vol. 114, pp. 823-828, 2014.
[7]   C.H. Scholz, “The frequency-magnitude relation of microfacturing in rock and its relation to earthquakes,” Bulletin of the Seismological Society of America, vol. 58, no. 1, pp. 399-415, 1968, doi: https://doi.org/10.1785/BSSA0580010399
[8]   P.R. Sammonds, P.G. Meredith dan I.G. Main, “Role of pore fluids in the generation of seismic precursors to shear fracture,” Letters to Nature, vol. 359, pp. 228-230, 1992.
[9]   S.  Rohadi,  H.  Grandis  dan  M.A.  Ratag,  “Studi  variasi  spasial  seismisitas  zona subduksi Jawa,” Jurnal Meteorologi dan Geofisika, vol. 8, no. 1, pp. 42-47, 2007, doi: http://dx.doi.org/10.31172/jmg.v8i1.4
[10]  Linda,  N.  Ihsan  dan  P.  Palloan,  “Analisis  distribusi  spasial  dan  temporal seismotektonik berdasarkan nilai b -value dengan menggunakan metode likelihood di Pulau Jawa,” Jurnal Sains dan Pendidikan Fisika (JSPF), vol. 15, no. 1, pp. 16- 31, 2019, doi: https://doi.org/10.35580/jspf.v15i1.9403
[11]  C.-H. Lin, “Temporal b -value variations throughout a seismic faulting process: the 2008  taoyuan  earthquake  in  Taiwan,” Terrestrial  Atmospheric  and  Oceanic Sciences, vol. 21, no. 2, 2010. doi: https://doi.org/10.3319/TAO.2009.02.09.01(T)
[12]  W. Marzocchi dan L. Sandri, “A review and new insights on estimation of the b - value and its uncertainty,” ANNALS  OF  GEOPHYSICS, vol.  46,  pp.  1271-1282, 2003, doi: https://doi.org/10.4401/ag-3472
[13]  K.  Aki,  “Maximum  likelihood  estimate  of b in  the  formula  logN=a-bM  and  its confidence limits,” Bulletin of Earthquake Research Institute of the University of Tokyo, vol. 43, pp. 237-239, 1965.
[14]  L.A. Maslov dan V.M. Anokhin, “Derivation of the Gutenberg-Richter empirical formula from the solution of the generalized logistic equation,” Natural  Science, vol. 4, pp. 648-651, 2012, doi: https://doi.org/10.4236/ns.2012.428085
[15]  Q.  Han,  L.  Wang,  J.  Xu,  A.  Carpinteri,  dan  G.  Lacidogna,  “A  robust  method  to estimate  the b -value  of  the  magnitude-frequency  distribution  of  earthquake,” Chaos,  Solitions  and  Fractals, vol.  81,  pp.  103-110,  2015,  doi: https://doi.org/10.1016/j.chaos.2015.09.004
[16]  X.  Bai,  W.  Yao  dan  J.E.  Boyer,  “Robust  fitting  of  mixture  regression  models,” Computational Statistics and Data Analysis, vol. 56, no. 7, pp. 2347-2359, 2012, doi: https://doi.org/10.1016/j.csda.2012.01.016
[17]  NEIC,  “USGS,  https://earthquake.usgs.gov/earthquakes.,”  1  Agustus  2022. [Online].
[18]  P. Reasenberg, “Second-order moment of Central California Seismicity,” Journal of  Geophysical  Research, vol.  90,  no.  B7,  pp.  5479-5495,  1985,  doi: https://doi.org/10.1029/JB090iB07p05479
[19]  D. Schorlemmer dan M.C. Gerstenberger,  “RELM  testing  center,” Seismological Research Letters, vol. 78, no. 1, pp. 30-36, 2007, doi: https://doi.org/10.1785/gssrl.78.1.30
[20]  T. van Stiphout, J. Zhuang dan D. Marsan, “Seismicity declustering,” Community Online Resource for Statistical Seismicity Analysis, 2012. doi: https://doi.org/10.5078/corssa-52382934
[21]  A.  Mignan  dan  J.  Woessner,  “Bayesian  estimation  of  the  spatially  varying completeness  magnitude  of  earthquake  catalogs,” Bulletin  of  the  Seismological Society  of  America, vol.  101,  no.  3,  pp.  1371-1385,  2011,  doi: https://doi.org/10.1785/0120100223
[22]  S.  Wiemar  dan  M.  Wyss,  “Minimum  magnitude  of  completeness  in  earthquake catalogs: examples from Alaska, the Western United States, and Japan,” Bulletin of the  Seismological  Society  of  America  vol.  90,  no.  4,  pp.  859-869,  2000,  doi: https://doi.org/10.1785/0119990114
[23]  E.D.  Pradewi  dan  Sudarno,  “Kajian  estimasi-M  IRLS  menggunakan  fungsi pembobot Huber dan Bisquare Tukey pada data ketahanan pangan di Jawa Tengah,” Media Statistika, vol.
5, no. 1, pp. 1-10, 2012, doi: https://doi.org/10.14710/medstat.5.1.1-10
[24]  I.G.  Main,  P.G.  Meredith  dan  P.R.  Sammonds,  “Temporal  variations  in  seismic event  rate  and b-values from stress  corrosion  constitutive laws,” Tectonophysics, vol. 211, pp. 233-246, 1992, doi: https://doi.org/10.1016/0040-1951(92)90061-A
[25]  D.C. Montgomery, E.A. Peck dan G.G. Vining, Introduction to Linear Regression Analysis, New York: John Wiley & Sons, 2012.
[26]  P. Nuannin, O. Kulhanek dan L. Persson, “Spatial and temporal b value anomalies preceding the devastating off coast of NW Sumatra Earthquake of December 26, 2004,” Geophysical  Research  Letters, vol.  32,  no.  11,  p.  11307,  2005,  doi: https://doi.org/10.1029/2005GL022679
[27]  A.  Soehaimi,  Marjiyono  dan  H.  Wahyono,  Peta  seismotektonik  daerah  Palu  dan sekitarnya,  Sulawesi  Tengah,  skala  1  :  150.000, Bandung:  Pusat  Survei  Geologi, 2013.