Analisis Spasial Kerentanan Kondisi Tanah terhadap Guncangan Gempa Bumi di Jawa Barat Berdasarkan Vs30 Menggunakan universal Kriging

Bab 1. Pendahuluan

1.1 Latar Belakang

Jawa Barat merupakan salah satu wilayah yang perlu mendapat perhatian dalam kajian kegempaan karena memiliki riwayat kejadian gempa bumi yang tersebar pada beberapa lokasi. Informasi gempa bumi yang dipublikasikan BMKG (2026) memuat parameter seperti waktu kejadian, magnitudo, kedalaman, koordinat, dan wilayah terdampak, sehingga data tersebut dapat digunakan untuk melihat sebaran kejadian gempa secara spasial. Dalam kajian bahaya gempa, dampak guncangan tidak hanya berkaitan dengan besarnya magnitudo atau kedalaman pusat gempa, tetapi juga dipengaruhi oleh kondisi tanah di lokasi tertentu. Oleh karena itu, analisis kondisi tanah menjadi penting untuk memahami wilayah mana yang secara karakteristik tapak lebih rentan terhadap penguatan guncangan.

Salah satu indikator yang umum digunakan untuk menggambarkan karakteristik tanah adalah Vs30, yaitu rata-rata kecepatan gelombang geser hingga kedalaman 30 meter. USGS (2025) menjelaskan bahwa Vs30 digunakan sebagai indikator kondisi tapak seismik karena dapat membedakan karakteristik tanah lunak, tanah sedang, tanah keras, hingga batuan. Nilai Vs30 yang lebih rendah umumnya menunjukkan tanah yang lebih lunak, sedangkan nilai Vs30 yang lebih tinggi menunjukkan tanah yang lebih keras. Dalam penelitian ini, kerentanan yang dimaksud bukan kerentanan sosial atau kerusakan bangunan, melainkan kerentanan kondisi tanah terhadap guncangan gempa berdasarkan karakteristik Vs30. Dengan demikian, data gempa bumi digunakan sebagai informasi overlay, sedangkan variabel yang dimodelkan secara spasial adalah nilai Vs30.

Pendekatan spasial diperlukan karena nilai Vs30 tidak berdiri sendiri sebagai angka pada satu titik, tetapi memiliki pola kedekatan antarwilayah. Lokasi yang berdekatan cenderung memiliki karakteristik tanah yang lebih mirip dibandingkan lokasi yang berjauhan, sehingga metode interpolasi spasial dapat digunakan untuk memperkirakan nilai Vs30 pada wilayah yang tidak memiliki titik pengamatan langsung. Esri (2025) menjelaskan bahwa Kriging merupakan metode geostatistik yang menghasilkan permukaan prediksi dari sekumpulan titik pengamatan dengan mempertimbangkan struktur spasial data melalui semivariogram. Penelitian ini menggunakan Universal Kriging karena hasil pemeriksaan awal menunjukkan adanya kecenderungan tren spasial pada nilai Vs30 berdasarkan koordinat lokasi. Dengan pendekatan tersebut, penelitian ini diharapkan dapat menghasilkan peta prediksi Vs30 dan peta ketidakpastian prediksi sebagai dasar untuk memahami variasi spasial kerentanan kondisi tanah terhadap guncangan gempa bumi di Jawa Barat.

1.2 Tujuan Penelitian

Berdasarkan latar belakang, dapat diidentifikasi beberapa permasalahan utama dalam penelitian ini, yaitu:

1. Mendeskripsikan sebaran kejadian gempa bumi dan nilai Vs30 di Jawa Barat.
2. Memetakan nilai Vs30 di Jawa Barat menggunakan Universal Kriging.
3. Mengidentifikasi wilayah dengan karakteristik tanah yang relatif lebih rentan terhadap guncangan gempa berdasarkan nilai Vs30

1.3 Pertanyaan Penelitian

Berdasarkan tujuan penelitian, pertanyaan prediksi pada penelitian ini ialah:

1. Bagaimana sebaran kejadian gempa bumi dan nilai Vs30 di Jawa Barat?
2. Bagaimana hasil pemetaan nilai Vs30 di Jawa Barat menggunakan Universal Kriging?
3. Wilayah mana yang memiliki karakteristik tanah relatif lebih rentan terhadap guncangan gempa berdasarkan nilai Vs30?

Bab 2. Data Dan Metode

2.1 Data Penelitian

Data penelitian ini terdiri dari data Vs30, data kejadian gempa bumi, dan data batas administrasi wilayah Jawa Barat. Data Vs30 diperoleh dari USGS Global Vs30 Mosaic Map Viewer. USGS menjelaskan bahwa Global Vs30 Mosaic menyediakan informasi Vs30, yaitu rata-rata kecepatan gelombang geser pada kedalaman 30 meter, dan data tersebut dapat diakses melalui aplikasi peta Vs30 USGS. Data gempa bumi diperoleh dari katalog United States Geological Survey (USGS) melalui pustaka ObsPy dengan layanan FDSN Event Web Service. USGS menjelaskan bahwa layanan FDSN Event Web Service dapat digunakan untuk melakukan pencarian informasi gempa berdasarkan parameter tertentu, seperti waktu, lokasi, magnitudo, dan kedalaman kejadian. Data gempa pada penelitian ini diambil untuk wilayah sekitar Jawa Barat dengan batas latitude -8 sampai -5 dan longitude 105 sampai 109 pada periode 2000 sampai 2026. Data batas administrasi diperoleh dari Geoportal Badan Informasi Geospasial (BIG) berupa batas administrasi desa/kelurahan Jawa Barat, yang digunakan sebagai pembatas wilayah analisis dan dasar visualisasi peta.

Unit observasi utama dalam penelitian ini adalah titik lokasi Vs30 di Jawa Barat. Setiap titik memiliki informasi koordinat longitude dan latitude serta nilai vs30_mosaic. Variabel vs30_mosaic digunakan sebagai variabel utama yang diprediksi menggunakan Universal Kriging. Koordinat longitude dan latitude digunakan untuk membentuk objek spasial, kemudian ditransformasikan ke sistem koordinat UTM sehingga proses analisis jarak dapat dilakukan dalam satuan meter. Pada Universal Kriging, koordinat hasil transformasi digunakan sebagai variabel x_utm dan y_utm untuk mewakili tren spasial. Data gempa bumi digunakan sebagai data pendukung dengan variabel longitude, latitude, magnitude, depth_km, tahun, bulan, tanggal, dan kota. Data gempa tidak digunakan sebagai variabel prediksi, tetapi ditampilkan sebagai overlay pada peta hasil prediksi Vs30 untuk melihat posisi kejadian gempa terhadap karakteristik kondisi tanah.

2.2 Metode Universal Kriging

Metode yang digunakan dalam penelitian ini adalah Universal Kriging. Metode ini dipilih karena nilai Vs30 memiliki unsur lokasi dan menunjukkan indikasi tren spasial berdasarkan koordinat. Esri (2025) menjelaskan bahwa Universal Kriging digunakan ketika data memiliki komponen tren yang berubah menurut lokasi. Dalam penelitian ini, variabel yang diprediksi adalah vs30_mosaic, sedangkan koordinat hasil transformasi UTM, yaitu x_utm dan y_utm, digunakan sebagai komponen tren spasial. Model Universal Kriging ialah sebagai berikut:

\[ Z(s) = \beta_0 + \beta_1 x + \beta_2 y + \varepsilon(s) \]

Keterangan:

\(Z(s)\) = nilai Vs30 pada lokasi (s), (x) dan (y) adalah koordinat UTM

\(\varepsilon(s)\) = residual spasial

Pebesma (2025) menjelaskan bahwa dalam package gstat, Universal Kriging dapat dijalankan dengan formula seperti \(z ~ x + y\), sedangkan Ordinary Kriging menggunakan formula \(z ~ 1\). Oleh karena itu, model pada penelitian ini dijalankan dengan formula \(vs30_mosaic ~ x_utm + y_utm\). Struktur spasial residual dimodelkan menggunakan semivariogram. SciKit-GStat (2025) menjelaskan bahwa variogram menggambarkan hubungan antara jarak antar titik dengan tingkat perbedaan nilai observasi. Semivariogram empiris pada penelitian ini dicocokkan dengan beberapa model teoritis, yaitu Spherical, Exponential, dan Gaussian. Model variogram terbaik dipilih berdasarkan nilai error fitting yang paling kecil dan bentuk kurva yang paling mengikuti semivariogram empiris. Dalam Kriging, bobot prediksi tidak ditentukan secara manual, tetapi dihitung dari model variogram. Esri (2025) menjelaskan bahwa Kriging menggunakan bobot dari semivariogram untuk memprediksi nilai pada lokasi yang tidak memiliki pengamatan.

Prediksi nilai Vs30 dilakukan pada grid wilayah Jawa Barat. Hasil utama dari model adalah nilai prediksi Vs30 dan nilai ketidakpastian prediksi. Esri (2025) menjelaskan bahwa nilai variansi prediksi yang lebih rendah menunjukkan tingkat keyakinan prediksi yang lebih tinggi, sedangkan nilai variansi yang lebih tinggi menunjukkan ketidakpastian yang lebih besar. Evaluasi model dilakukan menggunakan cross-validation dengan ukuran ME, MAE, RMSE, R² CV, dan korelasi antara nilai observasi dan prediksi. Dokumentasi gstat menjelaskan bahwa fungsi krige.cv dapat digunakan untuk validasi silang pada simple, ordinary, maupun universal kriging (Pebesma, 2025). Oada penelitian ini, informasi waktu pada data gempa hanya digunakan sebagai informasi deskriptif, sedangkan model utama difokuskan pada pemetaan nilai Vs30 berdasarkan lokasi.

BAB 3. Hasil Dan Pembahasan

3.1 Deskriptif Data dan Sebaran Awal

Data Vs30 yang digunakan dalam penelitian ini berjumlah 42.727 titik setelah proses pembersihan data. Nilai Vs30 minimum sebesar 180 m/s dan maksimum sebesar 900 m/s, dengan rata-rata sebesar 450,01 m/s dan median sebesar 437 m/s. Nilai standar deviasi sebesar 192,68 menunjukkan bahwa nilai Vs30 di Jawa Barat memiliki variasi yang cukup besar antarwilayah. Berdasarkan rentang nilainya, wilayah kajian tidak hanya didominasi oleh satu jenis karakteristik tanah, tetapi mencakup kondisi tanah yang relatif lebih lunak hingga lebih keras.

Data kejadian gempa bumi yang digunakan berjumlah 234 kejadian. Magnitudo gempa berada pada rentang 0 sampai 6,5, dengan rata-rata magnitudo sebesar 4,40 dan median sebesar 4,4. Kedalaman gempa berkisar antara 1,41 km sampai 329,6 km, dengan rata-rata kedalaman sebesar 83,62 km. Hal ini menunjukkan bahwa data gempa yang digunakan mencakup kejadian dengan kedalaman dangkal hingga dalam. Dalam penelitian ini, data gempa tidak digunakan sebagai variabel yang diprediksi, tetapi digunakan sebagai overlay untuk melihat posisi kejadian gempa terhadap hasil pemetaan Vs30.

Tabel 3.1. Statistik Deskriptif Data Vs30 dan Gempa Bumi

Data	n	Minimum	Rata-rata	Median	Maksimum	Standar Deviasi
Vs30	42.727	180	450,01	437	900	192,68
Magnitudo gempa	234	0	4,40	4,4	6,5	-
Kedalaman gempa	234	1,41	83,62	74,7	329,6	-

Peta sebaran titik Vs30 menunjukkan bahwa titik pengamatan tersebar cukup rapat di wilayah Jawa Barat. Sebaran titik yang rapat ini penting karena proses interpolasi spasial sangat bergantung pada ketersediaan titik pengamatan di sekitar lokasi prediksi. Sementara itu, peta kejadian gempa menunjukkan bahwa titik gempa tidak tersebar merata, tetapi lebih banyak muncul pada beberapa area tertentu, terutama pada bagian barat dan selatan Jawa Barat. Pola ini menunjukkan bahwa analisis spasial diperlukan agar informasi lokasi tidak hilang dalam proses pembahasan data.

Gambar 3.1. Peta Sebaran Titik Vs30 di Jawa Barat

Gambar 3.2. Peta Sebaran Kejadian Gempa Bumi di Jawa Barat

3.2 Pemeriksaan Pola Spasial Data Vs30

Pemeriksaan pola spasial dilakukan untuk memastikan bahwa data Vs30 layak dianalisis menggunakan pendekatan Kriging. Hasil Moran’s I menunjukkan nilai sebesar 0,7808. Nilai tersebut menunjukkan adanya autokorelasi spasial positif yang kuat. Artinya, titik-titik yang saling berdekatan cenderung memiliki nilai Vs30 yang relatif mirip. Kondisi ini mendukung penggunaan metode Kriging karena metode tersebut memanfaatkan hubungan spasial antar titik pengamatan dalam proses prediksi.

Selain autokorelasi spasial, pemeriksaan stasioneritas dilakukan secara visual melalui grafik tren Vs30 terhadap koordinat X dan Y UTM. Pada grafik terhadap koordinat X, garis tren tidak sepenuhnya datar dan menunjukkan perubahan nilai Vs30 mengikuti arah barat-timur. Pada grafik terhadap koordinat Y, pola tren terlihat lebih jelas, yaitu nilai Vs30 meningkat pada bagian tertentu lalu menurun kembali pada koordinat Y yang lebih tinggi. Pola tersebut menunjukkan bahwa rata-rata Vs30 tidak sepenuhnya konstan di seluruh wilayah kajian. Oleh karena itu, Universal Kriging lebih sesuai digunakan karena metode ini dapat memasukkan komponen tren spasial melalui koordinat lokasi.

Gambar 3.3. Grafik Tren Vs30 terhadap Koordinat X

Gambar 3.4. Grafik Tren Vs30 terhadap Koordinat Y

Kerapatan titik pengamatan juga diperiksa melalui histogram jarak nearest neighbor. Histogram menunjukkan bahwa sebagian besar titik Vs30 memiliki jarak yang relatif dekat dengan titik terdekatnya, terutama pada kisaran kurang dari 2.000 meter. Hal ini menunjukkan bahwa data pengamatan cukup rapat untuk mendukung proses interpolasi. Namun, masih terdapat beberapa titik dengan jarak yang lebih jauh, sehingga pada wilayah tertentu hasil prediksi tetap perlu dilihat bersama peta ketidakpastian.

3.3 Hasil Pemodelan Variogram Universal Kriging

Pemodelan variogram dilakukan untuk menggambarkan struktur spasial residual pada Universal Kriging. Pada penelitian ini, tiga model variogram teoritis dibandingkan, yaitu Spherical, Exponential, dan Gaussian. Berdasarkan hasil fitting, model Exponential memiliki nilai SSErr paling kecil, yaitu 5.376.743,52. Nilai tersebut lebih kecil dibandingkan model lain yang diuji. Oleh karena itu, model Exponential dipilih sebagai model variogram terbaik untuk proses Universal Kriging.

Tabel 3.2. Perbandingan Model Variogram Universal Kriging

Model	SSErr	Nugget	Partial Sill	Total Sill	Range
Exponential	5.376.743,52	6.557,63	26.344,10	32.901,73	21.164,40
Spherical	11.624.913,28	10.849,90	21.313,42	32.163,32	58.438,94
Gaussian	19.803.820,04	12.542,58	19.228,07	31.770,65	24.661,09

Model Exponential memiliki nugget sebesar 6.557,63, partial sill sebesar 26.344,10, total sill sebesar 32.901,73, dan range sebesar 21.164,40 meter. Nugget menunjukkan variasi pada jarak yang sangat dekat atau variasi yang tidak tertangkap oleh model. Partial sill menunjukkan variasi spasial yang dapat dijelaskan oleh struktur spasial, sedangkan total sill menunjukkan keseluruhan variasi yang terbentuk. Range sebesar sekitar 21,16 km menunjukkan bahwa hubungan spasial antar titik Vs30 masih cukup kuat sampai jarak tersebut. Setelah melewati jarak tersebut, keterkaitan spasial antar titik cenderung semakin melemah.

Gambar 3.5. Perbandingan Model Variogram Spherical, Exponential, dan Gaussian

3.4 Hasil Prediksi Vs30 dan Ketidakpastian Prediksi

Prediksi nilai Vs30 dilakukan pada 1.485 titik grid di wilayah Jawa Barat menggunakan Universal Kriging. Hasil cross-validation menunjukkan nilai ME sebesar 0,965, MAE sebesar 72,19, RMSE sebesar 102,20, R² CV sebesar 0,724, dan korelasi antara nilai observasi dan prediksi sebesar 0,851. Nilai ME yang mendekati nol menunjukkan bahwa model tidak memiliki bias prediksi yang besar. Nilai korelasi sebesar 0,851 menunjukkan hubungan yang kuat antara nilai Vs30 hasil observasi dan hasil prediksi. Dengan demikian, model Universal Kriging dapat dikatakan cukup baik dalam memetakan pola spasial Vs30 di wilayah kajian.

Tabel 3.3. Hasil Cross-Validation Universal Kriging

N	ME	MAE	RMSE	R² CV	Korelasi Observasi-Prediksi
2.000	0,965	72,19	102,20	0,724	0,851

Peta prediksi Vs30 menunjukkan bahwa bagian utara Jawa Barat cenderung memiliki nilai Vs30 yang lebih rendah dibandingkan beberapa wilayah tengah dan selatan. Nilai Vs30 yang lebih rendah menunjukkan karakteristik tanah yang relatif lebih lunak, sedangkan nilai Vs30 yang lebih tinggi menunjukkan tanah yang relatif lebih keras. Berdasarkan hasil klasifikasi grid prediksi, sebagian besar wilayah berada pada kelas SC atau tanah keras, yaitu sebanyak 958 grid atau 64,51%. Kelas SD atau tanah sedang/lunak mencakup 470 grid atau 31,65%. Sementara itu, kelas SB atau batuan mencakup 53 grid atau 3,57%, dan kelas SE atau tanah sangat lunak hanya mencakup 4 grid atau 0,27%.

Gambar 3.6. Peta Prediksi Vs30 di Jawa Barat

Peta ketidakpastian prediksi menunjukkan bahwa sebagian besar wilayah memiliki nilai simpangan baku Kriging yang relatif rendah. Ketidakpastian yang lebih tinggi tampak pada beberapa area pinggir wilayah, terutama pada bagian yang lebih jauh dari sebaran titik pengamatan. Hal ini wajar karena prediksi Kriging akan lebih stabil pada wilayah yang memiliki banyak titik pengamatan di sekitarnya. Oleh karena itu, peta prediksi Vs30 perlu dibaca bersama peta ketidakpastian agar interpretasi tidak hanya didasarkan pada nilai prediksi, tetapi juga tingkat keyakinan model.

Gambar 3.7. Peta Ketidakpastian Prediksi Vs30

3.5 Overlay Kejadian Gempa terhadap Peta Prediksi Vs30

Overlay kejadian gempa dilakukan untuk melihat posisi titik gempa terhadap peta prediksi Vs30. Hasil overlay menunjukkan bahwa sebagian besar titik gempa berada pada wilayah dengan kelas SC atau tanah keras, yaitu sebanyak 185 kejadian atau 79,06%. Sebanyak 39 kejadian atau 16,67% berada pada kelas SD atau tanah sedang/lunak. Sementara itu, 9 kejadian atau 3,85% berada pada kelas SB atau batuan, dan 1 kejadian atau 0,43% berada pada kelas SE atau tanah sangat lunak.

Hasil ini menunjukkan bahwa lokasi kejadian gempa pada data penelitian tidak hanya berada pada wilayah dengan Vs30 rendah, tetapi tersebar pada beberapa kelas kondisi tanah. Namun, dari sisi kerentanan kondisi tanah, wilayah dengan Vs30 yang lebih rendah tetap perlu diperhatikan karena tanah yang lebih lunak secara konsep lebih berpotensi mengalami penguatan guncangan dibandingkan tanah yang lebih keras. Dengan demikian, overlay ini tidak digunakan untuk menyimpulkan bahwa kelas tanah tertentu menyebabkan gempa, tetapi untuk melihat posisi kejadian gempa terhadap karakteristik tanah hasil interpolasi.

Gambar 3.8. Overlay Kejadian Gempa Bumi terhadap Peta Prediksi Vs30

3.6 Interpretasi Substantif dan Keterbatasan

Secara substantif, hasil penelitian menunjukkan bahwa kerentanan kondisi tanah terhadap guncangan gempa bumi di Jawa Barat tidak tersebar secara merata. Wilayah dengan nilai Vs30 lebih rendah menunjukkan karakteristik tanah yang relatif lebih lunak dan dapat dianggap lebih rentan terhadap penguatan guncangan. Sebaliknya, wilayah dengan nilai Vs30 lebih tinggi menunjukkan kondisi tanah yang relatif lebih keras. Hasil klasifikasi menunjukkan bahwa sebagian besar grid prediksi berada pada kelas tanah keras, tetapi masih terdapat wilayah tanah sedang/lunak yang cukup luas, terutama pada bagian utara dan beberapa area pinggir wilayah kajian.

Penggunaan Universal Kriging sesuai dengan karakteristik data karena nilai Vs30 menunjukkan autokorelasi spasial dan indikasi tren berdasarkan koordinat. Model Exponential dipilih sebagai model variogram terbaik karena memiliki nilai error fitting paling kecil. Hasil validasi silang juga menunjukkan performa model yang cukup baik, ditunjukkan oleh korelasi observasi-prediksi yang kuat dan nilai R² CV sebesar 0,724. Dengan demikian, hasil pemetaan dapat digunakan sebagai gambaran awal untuk melihat variasi spasial karakteristik tanah di Jawa Barat.

Penelitian ini tetap memiliki beberapa keterbatasan. Pertama, variabel yang dimodelkan hanya Vs30, sehingga hasil penelitian hanya menggambarkan kerentanan dari sisi kondisi tanah. Kedua, data gempa hanya digunakan sebagai overlay dan tidak dimasukkan sebagai variabel prediksi. Oleh karena itu, penelitian ini tidak memprediksi kejadian gempa, magnitudo, PGA, maupun kerusakan. Ketiga, hasil prediksi dipengaruhi oleh kerapatan dan sebaran titik Vs30, sehingga wilayah dengan titik pengamatan lebih jarang memiliki ketidakpastian yang lebih besar. Keempat, penelitian ini menggunakan pendekatan spasial saja dan belum memasukkan struktur waktu, sehingga belum dapat menggambarkan perubahan pola kegempaan secara temporal.

BAB 4. Kesimpulan

Berdasarkan hasil analisis, data Vs30 di Jawa Barat menunjukkan variasi spasial yang cukup jelas, dengan nilai minimum 180 m/s, maksimum 900 m/s, dan rata-rata 450,01 m/s. Hasil Moran’s I sebesar 0,7808 menunjukkan adanya autokorelasi spasial positif yang kuat, sehingga pendekatan spasial layak digunakan. Grafik tren terhadap koordinat X dan Y juga menunjukkan adanya kecenderungan perubahan nilai Vs30 mengikuti arah lokasi, sehingga Universal Kriging sesuai digunakan dalam penelitian ini.

Model variogram terbaik pada Universal Kriging adalah model Exponential karena memiliki nilai SSErr paling kecil dibandingkan model lainnya. Hasil cross-validation menunjukkan performa model yang cukup baik, dengan RMSE sebesar 102,20, R² CV sebesar 0,724, dan korelasi observasi-prediksi sebesar 0,851. Peta prediksi menunjukkan bahwa sebagian besar wilayah Jawa Barat berada pada kelas SC atau tanah keras, tetapi masih terdapat wilayah dengan kelas SD atau tanah sedang/lunak yang perlu diperhatikan karena secara karakteristik lebih rentan terhadap penguatan guncangan gempa. Dengan demikian, hasil penelitian ini dapat digunakan sebagai gambaran awal kerentanan kondisi tanah terhadap guncangan gempa bumi berdasarkan nilai Vs30, bukan untuk memprediksi kejadian gempa, magnitudo, PGA, atau tingkat kerusakan.

Daftar Pustaka

Badan Meteorologi, Klimatologi, dan Geofisika. (2026). Gempa bumi dirasakan. https://www.bmkg.go.id/gempabumi/gempabumi-dirasakan

Esri. (2025a). How kriging works. ArcGIS Pro Documentation. https://doc.esri.com/en/arcgis-pro/latest/tool-reference/spatial-analyst/how-kriging-works.html

Esri. (2025b). Kriging (Spatial Analyst tools). ArcGIS Pro Documentation. https://doc.esri.com/en/arcgis-pro/latest/tool-reference/spatial-analyst/kriging.html?tabs=dialog

Esri. (2025c). Kriging in Geostatistical Analyst. ArcGIS Pro Documentation. https://doc.esri.com/en/arcgis-pro/latest/help/analysis/geostatistical-analyst/kriging-in-geostatistical-analyst.html

Esri. (2025d). Understanding universal kriging. ArcGIS Pro Documentation. https://doc.esri.com/en/arcgis-pro/latest/help/analysis/geostatistical-analyst/understanding-universal-kriging.html

Pebesma, E. (2025). gstat: Spatial and spatio-temporal geostatistical modelling, prediction and simulation. CRAN. https://cran.r-project.org/web/packages/gstat/gstat.pdf

SciKit-GStat. (2025). Variography. SciKit-GStat Documentation. https://scikit-gstat.readthedocs.io/en/latest/userguide/variogram.html

U.S. Geological Survey. (2025). Vs30 models and data. USGS Earthquake Hazards Program. https://earthquake.usgs.gov/data/vs30/