0% found this document useful (0 votes)
206 views92 pages

Seismic Unix Program

h

Uploaded by

Nurlia Adu
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
206 views92 pages

Seismic Unix Program

h

Uploaded by

Nurlia Adu
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

Low velocity

High velocity

SEISMIC UNIX ( SU ) AT A GLANCE


VOL.1 INTRODUCTION IN SEISMIC
DATA PROCESSING

RANDY A.K. CONDRONEGORO

DEPARTEMEN TEKNIK GEOFISIKA


INSTITUT TEKNOLOGI BANDUNG
2003
P r e f a c e

You are in a maze of twisty little passages, all alike

First of all i really would like to thank the following peoples who gave me big
influences:

1. DR. Eng T.A Sanny


2. DR. Sony Winardhie
3. Sri Widiyantoro, PhD
4. DR. Bagus Endar B Nurhandoko
5. Gunawan Handayani, MSCE, PhD

I assume that the readers of this book already have prior knowledge about Unix or
linux environment, shell programming and Matlab software. Also a knowledge in
C/C++ programming language would be advantageous because the source code of the
Seismic Unix is originally written in this language, you can always look in the source
code every time you like to understand how the mathematical concept in seismic
processing translated numerically ( believe me, itll be very useful for you to
understand that anything who has infinite and continuous values and parameters in the
real world should be translated in finite and discrete values in digital world, it will
bring you sense of awareness about assumption and limitation in numerical
programming that even will bring you in a deeper perceptive about geophysics ), after
all the purpose of open source is that the one that use it can always look and modify
the source code ( under circumstances described in the GNU General Public License
from Free Software Foundation ).

Well, a book is not a book if it not gives the author any feedback, so feel free to give
me any suggestions and any kind of critics. You can email me at
[email protected] .

Bandung , 2003

Randy Condronegoro
and human knowledge belong to the world
CHAPTER 1
LITTLE ENDIAN ? BIG ENDIAN ?

Sebelum kita memulai dengan data seismik dalam format SEG-Y ada
baiknya untuk mengenal lebih dulu platform komputer yang kita gunakan
dimana perangkat lunak SU di install.
Berdasarkan dari cara meletakkan data dalam komputer, terdapat 2 macam
metode yang telah digunakan yaitu Little Endian (LE) dan Big Endian
(BE). LE biasanya digunakan pada mesin berprosesor Intel dan DEC
sementara BE biasanya digunakan pada mesin berprosesor IBM,
Motorolla, SUN dan SGI.
Metode BE diterapkan dengan cara meletakkan bit dengan nilai eksponen
atau pangkat terbesar pada awal file ( Most Significant Bit ( MSB ) pada
awal file ). Sedangkan pada metode LE, MSB diletakkan pada akhir file.
Untuk lebih jelasnya dapat dilihat pada contoh dibawah :
Misal terdapat suatu nilai desimal 2100 maka apabila ditransformasi
kedalam bentuk heksadesimal :
2100 : 16 = 131 sisa 4
131 : 16 = 8 sisa 3
8 : 16 = 0 sisa 8
sehingga 2100 (desimal) = 834 (heksadesimal)
artinya 8 x 162 = 2048
3 x 161 = 48
4 x 160 = 4 +

2100
Bilangan Heksadesimal ini dapat dinotasikan dalam bentuk bilangan biner
4 bit seperti 1100 , 0001 dan seterusnya ( sebetulnya bilangan desimal
diatas dapat langsung ditransformasikan dalam bentuk biner, penggunaan
heksadesimal dipakai untuk mempermudah pengertian ) .
Sehingga didapat :
8 = 1000 ( 1 x 23 + 0 x 22 + 0 x 21 + 0 x 20 )
3 = 0011 ( 0 x 23 + 0 x 22 + 1 x 21 + 1 x 20 )
4 = 0100 ( 0 x 23 + 1 x 22 + 0 x 21 + 0 x 20 )
Apabila digunakan format 16 bit atau 2 byte ( 1 byte terdiri dari 8 bit )
akan menjadi :

0000 1000 0011 0100

Beginning End of File


of File 0 8 3 4 (EOF)
(BOF)

0 x 163 8 x 162 3 x 161 4 x 160


Most Significant Bit Least Significant Bit

Format diatas dikenal sebagai Big Endian Order dimana MSB atau Most
Significant Bit atau bilangan dengan pangkat terbesar ( dalam hal ini 163 )
terletak pada awal file ( BOF ), sedangkan untuk LE akan tersusun
kebalikannya yaitu :

0011 0100 0000 1000

Beginning End of File


of File 3 4 0 8 (EOF)
(BOF)

4 x 161 4 x 160 0 x 163 8 x 162

1 byte 1 byte
Least Significant Bit Most Significant Bit

dimana terlihat bahwa MSB terletak pada akhir file ( EOF ).


Pengertian mengenai BE dan LE ini menjadi penting apabila pertukaran
data antar platform dilakukan, misal dari platform IBM ke Intel, biasanya
data seismik ( SEG-Y) ditulis oleh mesin berprosesor IBM atau SUN dan
SGI ( karena kebanyakan menggunakan environment Unix ) sehingga data
akan ditulis dengan format BE ,apabila kita ingin membaca data tersebut
dibaca pada mesin berprosesor Intel ( seperti kebanyakan PC di Indonesia
dan dengan menggunakan program selain SU , karena SU pada saat
penginstalan sudah terdapat opsi untuk mesin berformat BE atau LE ) dan
sistem operasinya DOS/Windows maka data yang ditulis dalam format BE
akan diperlakukan sebagai data berformat LE, hal ini tentunya akan
menghasilkan error pada saat membaca data.
Penulis mempunyai sedikit pengalaman mengenai hal ini, pada saat
penulis akan membaca data yang didownload dari Centre for Wave
Phenomena ( data Ozdogan Yilmaz yang dapat didownload gratis dari
www.cwp.mines.edu ) dengan menggunakan Matlab dengan perintah
fopen dan opsi penulis gunakan adalah n atau native ( local machine
format ) artinya LE karena penulis menggunakan PC berprosesor Intel )
pada saat pembacaan jumlah sampel per trace nilai yang diberikan berbeda
dengan yang terdapat pada header trace, setelah penulis menggunakan opsi
b atau menggunakan format BE maka nilai jumlah sampel per trace
terbaca dengan benar dan data dapat dibaca oleh Matlab dengan lancar.
(untuk lebih jelas mengenai Pembacaan data SEG-Y dengan Matlab dapat
mendownload rutin program SeismicLab dari website M.D.Sacchi Dept.
of Physics, University of Alberta,Canada ( www.ualberta.ca ), atau
SegyMAT ( http://segymat.sourceforge.net ) oleh Thomas Mejer Hansen,
Niels Bohr Institute for Astronomy, Physics and Geophysics (NBIfAFG),
University of Copenhagen, Denmark.

SU Gunakan suswapbytes atau swapbytes untuk menswap atau membalik


tips data per bytenya ( dari LE ke BE dan vice versa )

Historical Note : Mengapa sih ada BE dan LE, wellBE digunakan karena dari sononya
umat manusia terbiasa dengan format ini, coba saja bilangan tahun, misal tahun 2003
terlihat kan bahwa 2003 terdiri dari 2 x 103 + 0 x 102 + 0 x 101 + 3 x 100 jadi pangkat
tertinggi ada disebelah kiri atau diawal suatu bilangan, tetapi kemudian pada praktek
penulisan data terdapat kesulitan apabila data bit yang lebih rendah akan
ditransformasikan menjadi bit bit yang lebih tinggi, misal: 8 bit menjadi 16, atau 16 bit
menjadi 32 bit., maka format BE akan menghasilkan data yang berbeda padahal
semestinya nilainya sama, sehingga kemudian timbullah format LE ( misal pada contoh
diatas nilai 2100=08 34 heksa 16 bit apabila dijadikan 32 bit menjadi 08 34 00 00, nilai
ini tentunya tidak sama lagi dengan 2100 karena sekarang 8 heksa berkaitan dengan 166,
nilai 3 berkaitan dengan 165 dst, sementara dengan format LE 2100=34 08 menjadi 34 08
00 00 akan menghasilkan nilai yang sama yaitu 2100 karena disini nilai 8 heksa tetap
berkaitan dengan 162 dst , lihat juga http://www.princeton.edu/~kazad/index.htm atau
http://www.noveltheory.com/TechPapers/Little Endian vs_ Big Endian.htm ) .
CHAPTER 2
FLOATING POINT REPRESENTATION, IN THE SHADOW OF THE
INTERNATIONAL BUSSINESS MACHINE (IBM)

Dan kalian berpikir bahwa masalah sudah selesai, dengan sangat menyesal
penulis katakan bahwa permasalahan LE/BE diatas hanya merupakan
sebagian dari keseluruhan masalah yang ada. Permasalahan lain yang
harus dimengerti sebelum kita dapat membaca data seismik (SEG-Y)
dengan benar adalah mengambil nilai bilangan floating point atau bilangan
percahan yang benar.
Kita mulai dengan sedikit pelajaran sejarah..
Pada saat zaman pencerahan terjadi (yaitu zaman komputer mulai
mendominasi dan era digital mulai menggusur sistem analog) antara tahun
1970 dan 1980, jagad komputer dikuasai sebagian besar oleh perusahaan
IBM dimana perusahaan ini memproduksi mainframe-mainframe yang
digunakan dimana-mana (belum lagi booming PC yang dimulai oleh
perusahaan yang sama). IBM juga akhirnya banyak menghasilkan standar-
standar yang secara de facto digunakan oleh sebagian vendor-vendor
komputer yang lain dan alhasil juga dijadikan standar data oleh suatu
komunitas tertentu yang menggunakan mainframe IBM dalam
penyimpanan datanya. Salah satu dari hasil kerja dari IBM adalah
standarisasi representasi bilangan dalam komputer, yang salah satunya
kemudian dipakai oleh komunitas geofisika adalah bilangan mengambang
atau pecahan atau floating point yaitu 32 bits IBM Floating Point.
Tetapi selain IBM, kemudian banyak perusahaan komputer yang
mengeluarkan standar mereka sendiri-sendiri yang mengakibatkan
kompatibilitas data semakin berkurang sebagai hasil perbedaan format data
ini.Xerox Data Systems (XDS), Modcomp,GE,CDC adalah sebagian
contoh format floating point dari vendor masing-masing (walaupun untuk
dunia geofisika format data lain tidak terlalu berpengaruh karena relatif
tidak digunakan tetapi untuk mahasiswa Teknik Informatika pasti menjadi
hal yang cukup menjengkelkan).
Dikemudian hari karena masalah kompatibilitas maka perlu dibentuk suatu
format data standar untuk merepresentasikan suatu tipe bilangan yang
sama, karenanya kemudian oleh IEEE (Institute of Electronic and
Electrical Engineering) dikeluarkan standar-standar baru untuk setiap tipe
dan bagi bilangan mengambang yang digunakan oleh kalangan geofisika
adalah yang kemudian dikenal sebagai 32 bit IEEE Floating Point (dalam
IEEE 754 Standard Formats) dimana Intel juga mengadopsi standar ini.
Karena kalangan geofisika khususnya di industri perminyakan telah
menggunakan IBM Floating Point cukup lama sehingga sebagian besar
(mungkin hampir seluruhnya) masih dalam format IBM floating point, dan
juga penggunaan mesin-mesin dengan format dan platform yang selaras
dengan IBM masih tetap digunakan hingga sekarang (karena di Industri
perminyakan kebanyakan masih menggunakan IBM dan SGI sementara
SUN dan IBM kebanyakan digunakan oleh kalangan akademik luar negeri,
kalangan akademik kita dan negara berkembang lainnya didominasi oleh
PC -karena harga yang relatif jauh lebih murah dan kemampuan PC yang
semakin berkembang pesat- ) maka kedua representasi bilangan floating
point tadi digunakan keduanya sebagai standar data SEG-Y (standar yang
baru di revisi tahun 2002, sebelum 2002 yang digunakan adalah 32 bit
IBM floating point).
Cukup sekian pelajaran sejarahnya..
Representasi Bilangan Mengambang atau Pecahan
Apabila kita mempunyai suatu bilangan pecahan, misal: 375,414 maka
terdapat banyak cara untuk merepresentasikan bilangan tersebut diatas,
salah satunya adalah dalam bentuk Notasi Exponensial. Secara umum
notasi ini akan menyimpan informasi-informasi sebagai berikut
1. Tanda apakah bilangan positif atau negatif
2. Bilangan yang disimpan atau disebut Mantissa
3. Pangkat yang akan disimpan atau Eksponen
4. Tanda positif atau negatif dari Eksponen
5. Basis dari pangkat (apakah 10,8 atau 2 dsb)
lokasi nilai tanda eksponen
desimal (koma)
tanda mantissa eksponen

+375,414 x 10+0

mantissa basis pangkat

Dalam bentuk inipun terdapat banyak cara untuk menyatakan bilangan


diatas yaitu :
1. 375,414 x 100
2. 3,75414 x 102
3. 375414 x 10-3
4. 0,375414 x 103
5. 375414000 x 10-6
6. 0,000375414 x 106

atau kalau disertai dengan pembulatan


7. 37541 x 10-2
8. 0,003754 x 105
dan lain-lainnya. Tetapi semisal kita dibatasi hanya bisa menyimpan
sejumlah data bilangan saja maka dari pilihan diatas mana yang akan kita
pilih sebagai format kita? tentunya pilihan akan jatuh pada yang paling
sedikit menyita tempat untuk bilangan-bilangan tersebut yaitu pada nomor
1,2 dan 3. Tetapi kemudian apabila kita menggunakan pilihan 1 dan 2
maka komputer perlu diberitahu caranya menyimpan tanda koma
(informasi dimana pecahan dimulai) maka pilihan yang paling tepat
tentunya kemudian jatuh pada nomor 3, betul? tidak juga, terdapat satu
cara lagi untuk menyimpan data tersebut diatas dengan tingkat akurasi
yang lebih baik. Misal data yang akan kita simpan adalah nilai 375,4145
dan bukannya 375,414 dan seperti diatas kita dibatasi hanya mempunyai
tempat 6 digit untuk menyimpan bilangan kita ditambah beberapa digit
untuk eksponennya maka pilihan 1,2 dan 3 akan mengurangi tingkat
akurasi data kita karena digit desimal terakhir tidak dapat dimasukkan
dalam penyimpanan, hal ini bisa dihindari apabila kita bisa menyimpan
salah satu nilai diantara 375,4145 secara implisit atau secara tak langsung,
yaitu misal kita bisa menyimpan nilai 300 secara implisit maka data yang
harus disimpan hanya 75,4145 yang sesuai untuk 6 digit alokasi
penyimpanan. Atau yang disimpan secara implisit adalah nilai 0,0005
maka yang perlu disimpan adalah 375,414 yang juga sesuai dengan alokasi
kita. Seperti dijabarkan di Bab 1 nilai 300 adalah tempat Most Significant
Bit sementara nilai 0.0005 adalah nilai Least Significant Bit. 32 bit Single
Precision IEEE Floating Point menerapkan metoda yang pertama untuk
meningkatkan tingkat akurasinya yaitu most significant bit disimpan secara
implisit, dalam IEEE yang digunakan adalah perpangkatan dengan basis 2
dan bukannya basis 10 seperti contoh diatas.
Untuk menyimpan bilangan eksponen, maka diterapkan metode titik
tengah atau biasa disebut bias, misalkan bilangan eksponen mempunyai
tempat sebanyak 2 digit yang berarti nilai terendah 0 dan tertinggi adalah
99 maka nilai eksponen negatif diletakkan di 50 bilangan pertama dan nilai
eksponen positif diletakkan di 50 bilangan berikutnya, misal :
10+2 mempunyai bilangan 50+2 = 52
10-3 mempunyai bilangan 50-3 = 47
Untuk contoh bilangan 375,414 representasi dalam format IEEE adalah :
untuk 375 ->
1 x 28 + 0 x 27 + 1 x 26 + 1 x 25 + 1 x 24 + 0 x 23 + 1 x 22 + 1 x 21 + 1x 20
= 256 + 0 + 64 + 32 + 16 + 0 + 4 + 2 + 1
dalam bentuk biner :
101110111
untuk 0,414 ->
= 0x2-1 + 1x2-2 + 1x2-3 + 0x2-4 + 1x2-5 + 0x2-6 + 0x2-7 + 1x2-8 + 1x2-9 +
1x2-10 + 1x2-11 + 1x2-12 + 1x2-13 0x2-14 + 1x2-15
= 0.25 + 0.125 + 0.03125 + 0.00390625 + 0.001953125+ 0.0009765625
+0.00048828125+0.000244140625+0.0001220703125+0.0000305175781
25 (dst tetapi karena hanya disediakan 23 bit maka hanya diambil sampai
titik ini, sisa menjadi round off error/truncation error)
dalam bentuk biner :
011010011111101
total dalam bentuk biner :
101110111,011010011111101
apabila titik desimal kita geser (dinormalisasi) sehingga menjadi :
1,01110111011010011111101 x 28
dan karena Most Significant Bits disimpan secara implisit, maka nilai biner
dari mantissa diatas yang disimpan adalah :
01110111011010011111101
untuk IEEE eksponen disediakan sebanyak 8 bit sehingga mencakup nilai
0 hingga 256 dan titik bias di 127 sehingga nilai eksponen diatas apabila
dimasukkan bersama biasnya :
127 + 8 = 135
dalam bentuk biner 135 adalah :
1 x 27 + 0 x 26 + 0 x 25 + 0 x 24 + 0 x 23 + 1 x 22 + 1 x 21 + 1x 20
= 128+0+0+0+0+4+2+1
= 10000111
tanda atau sign bernilai 0 (karena positif,nilai 1 untuk negatif)
sehingga nilai yang dimasukkan secara total adalah (format Big Endian):
0 10000111 01110111011010011111101

mantissa

Eksponen+bias(127)
sign
berikut adalah format IEEE dan IBM untuk kategori 32 bit,

Format 32 bits IEEE Single Precision Floating Point:


Sign (1 bit) Exponent (8 bits) Mantissa (23 bits)
S EEEEEEEE MMMMMMMMMMMMMMMMMMMMMMM

Most Significant Bits = Implicit/Hidden


Exponent base = 2

Format 32 bits IBM Single Precision Floating Point:


Sign (1 bit) Exponent (7 bits) Mantissa (24 bits)
S EEEEEEE MMMMMMMMMMMMMMMMMMMMMMMM

Most Significant Bits =Not Implicit/Hidden


Exponent base = 16
CHAPTER 3
SEG-Y AND SU DATA STRUCTURE

Format data yang digunakan dalam perekaman data seismik biasanya


digunakan format data SEG-Y. Format data ini membagi data keseluruhan
menjadi 3 bagian yaitu :
1. EBCDIC ( Extended Binary-Coded-Decimal Interchange Code ) card
sebesar 3200 byte yang terdiri dari 40 baris dengan 80 karakter tiap
barisnya ( saat ini juga digunakan standar dari ASCII ( American
Standard for Information Interchange ) untuk menulis informasi text
sebesar 3200 byte ini sehingga bagian awal dari data SEG-Y ini dapat
dibaca oleh text editor biasa atau NC (DOS)/MC (Linux) ).
2. Binary header sebesar 400 byte yang berisi informasi tentang data
SEG-Y
3. Bagian data dari SEG-Y dimana terbagi menjadi 2 bagian yaitu :
a. Trace header sebesar 240 byte pada masing-masing trace.
Trace header ini berisi macam-macam informasi mengenai
data trace yang didahuluinya seperti : tracl, tracr, cdp, offset,
sx,gx, ns ( number of sample atau jumlah sampel tiap trace,
data ini sangat penting bagi SU atau program lainnya dalam
membaca dan menampilkan data, jadi pastikan data ns ini
ada!! ), dll . Besarnya data pada masing-masing variabel
diatas berbeda, ada yang bertipe short ( 16 bit = 2 byte) dan
ada yang bertipe int ( 32 bit = 4 byte ) dimana keduanya
bertipe integer / bilangan bulat. (lihat Appendix 1)
b. Data SEG-Y itu sendiri yang mempunyai format data 32 bit
floating point. Pada awalnya format yang dipakai adalah IBM
32 bit floating point tetapi untuk saat ini format dari IEEE
(Institute of Electronic and Electrical Engineering) yaitu
IEEE 32 bit floating point juga digunakan dan merupakan
standar dari SEG ( Society of Exploration Geophysics ).
Apabila digambarkan secara skematis dapat dilihat sebagai berikut :

3200 byte 400 byte 1st 1st Nth Nth


EBCDIC/ binary 240 byte Data 240 byte Data
ASCII file trace trace trace trace
textual header header header
file
header

Menurut SEG-Y rev 1 Data Exchange Format, SEG Technical Standards


Committee Release 1.0 May 2002, terdapat beberapa perubahan dari
format data diatas ( standard 1975 ), ada yang sifatnya opsional seperti
terlihat pada skema berikut ( untuk lebih jelasnya lihat pada www.seg.org
subsection standard ) :

3200 byte 400 byte 1st Nth 1st 1st Nth Nth
EBCDIC/ binary 3200 3200 240 Data 240 Data
byte byte byte trace byte trace
ASCII file trace trace
extend extend
textual header ed ed head head
file textual textual er er
header file file
header header

Data SU adalah merupakan data SEG-Y tanpa 3200 byte EBCDIC/ASCII


textual file header dan 400 byte binary file header, sehingga secara
skematis terlihat seperti :

1st 1st Nth Nth


240 byte Data 240 byte Data
trace trace trace trace
header header
CHAPTER 4
SYNTHETIC DATA, THE MAKING OF

Terdapat berbagai macam cara dalam SU untuk membuat data sintetik


diantaranya adalah suplane, suspike, sufdmod2, susyncz, susynlv,
susynvxz, susynlvcw, susynvxzcs dan lainnya. Yang akan disinggung
disini adalah pembuatan data sintetik dengan menggunakan fungsi susynlv
atau sintetik seismogram yang dibuat berdasarkan asumsi bahwa fungsi
kecepatan didalam bumi merupakan suatu fungsi linear ( fungsi kecepatan
merupakan orde pertama dengan jarak/kedalaman sebagai variabel
bebasnya ).
Sebelum membuat sebuah data sintetik perlu diketahui terlebih dahulu
koordinat-koordinat yang dapat digunakan dalam sintetik seismogram.
Terdapat dua macam koordinat yang dapat digunakan yaitu koordinat
berdasarkan posisi sumber dan penerima ( source dan receiver ) atau sx
dan gx serta yang kedua adalah koordinat berdasarkan posisi midpoint dan
offset atau mx dan ox. Ingat disini semuanya diasumsikan mempunyai arah
yang sama atau pada aksis yang sama yaitu x, penggambaran dibawah
hanya menggambarkan nilai posisinya dan bukan posisi secara nyata.

sx
ox

gx

mx
Beberapa masukan atau input yang dibutuhkan dalam fungsi susynlv
diantaranya adalah nxs ( jumlah shot yang dilakukan ), fxs ( first shot in x
direction atau posisi awal shot dalam km arah x ), dxs ( selang jarak untuk
setiap tembakan atau shot ), nxo ( jumlah offset yang ada antara source dan
receiver ), fxo ( nilai jarak offset terkecil atau pertama ), dxo ( selang
perbedaan antara satu offset dengan offset selanjutnya ), nxm ( jumlah
midpoint yang ada ), fxm ( posisi midpoint pertama ), dxm ( jarak antara
satu midpoint dengan midpoint selanjutnya ). Juga diperlukan input berupa
posisi reflektor dan amplitudonya serta nilai fungsi linear kecepatan
berupa perubahan kecepatan terhadap kedalaman ( v / z ) dan perubahan
kecepatan terhadap jarak horisontal ( v / x ). Untuk perekaman sintetik
seismogram maka juga diperlukan data input berupa nt ( jumlah sampling
waktu atau jumlah titik tiap tracenya ), ft ( waktu awal perekaman ) dan dt
( selang tiap titik perekaman ). Beberapa kemungkinan susunan sumber
dan penerima dapat dilihat pada skema dibawah :

Common shot

Common
receiver

geophone
source

Common
midpoint

Common
offset
Contoh dibawah menuliskan shell script untuk membuat susunan common
shot ( untuk single shot atau satu kali penembakan ) :
#! /bin/sh

# Examples of producing single shot records with susynlv


# We can add some noise to make data a bit more
# realistic, in particular, velocity analyses (with suvelan)
# seem very "blocky" without a little noise.
# routine are adopted from su/src/demos

WIDTH1=300
WIDTH2=600
HEIGHT=400
WIDTHOFF1=0 1
WIDTHOFF2=315
HEIGHTOFF1=100
HEIGHTOFF2=550

####Work in km (synlv puts 10^3 in header fields)####


# ref="1:1,2;4,2"
# reflector(s):"amplitude:x1,z1;x2,z2;x3,z3;..."
REF1="1.0:0.0,1.0;5.0,1.0"

n1=300 ft=0.0 dt=0.008 # time information: number,


# first, spacing
nxs=1 fxs=0.0 dxs=0.025 # shot information: number,
# first, spacing 2
nxo=30 fxo=0.05 dxo=0.05 # offset information: number,
# first, spacing

dvdz=0.0 # velocity gradient (dv/dz)


dvdx=0.0 # velocity gradient (dv/dx)
v00=2.0 # P-wave velocity at surface in km/s
ob=1 # include obliquity (both sides of
# interface are reflecting)
er=1 # use exploding reflector amps
smooth=0 # connect picks on reflectors linearly

susynlv nt=$n1 ft=$ft dt=$dt \


nxo=$nxo fxo=$fxo dxo=$dxo \
nxs=$nxs fxs=$fxs dxs=$dxs \ 3
v00=$v00 dvdx=$dvdx dvdz=$dvdz \
ob=$ob er=$er smooth=$smooth \
ref=$REF1 > dataSingleShot
suxwigb < dataSingleShot perc=99 \
title="Single Shot Records" windowtitle="Single Shot" \
label1="Time (sec)" label2="Receiver (km)" \ 4
f2=1 d2=0.02 \
-geometry ${WIDTH2}x${HEIGHT}+${WIDTHOFF2}+${HEIGHTOFF1} &

suchart < dataSingleShot | xgraph n=250 linecolor=1 \ 5


linewidth=0 marksize=3 mark=8 &

exit 0

Keterangan :
1. Input parameter untuk ukuran window dimana tampilan data akan
digambarkan
2. Input parameter untuk fungsi susynlv
3. pemanggilan fungsi susynlv dengan data keluaran berupa file dengan
nama datasingleshot
4. penggambaran data dengan fungsi suxwigb dengan masukan file
datasingleshot dengan percent clip sebesar 99 % dan input parameter
lain untuk legenda window
5. pemanggilan fungsi suchart untuk mendapatkan informasi posisi sx
dan gx dan xgraph untuk menggambarkan dalam koordinat sx,gx.
Hasil dari shell script diatas dapat dilihat pada gambar dibawah :
Common
shot

Untuk menggambarkan pola Common Receiver, Common Offset dan


Common Midpoint maka input parameter yang perlu diubah adalah pada
input nxo,nxm,nxs,fxo,fxm,fxs,dxo,dxm,dxs. Contoh shell script untuk
masing-masing pola diatas dapat ditulis sebagai berikut ( nama file
keluaran dapat berbeda, penulisan hanya pada bagian yang perlu diubah
saja, sementara yang lainnya sama dengan shell script Common Shot
diatas) :

Common Receiver
n1=300 ft=0.0 dt=0.008 # time information: number,
# first, spacing
nxs=30 fxs=0.0 dxs=0.05 # shot information: number,
# first, spacing
nxo=30 fxo=0.05 dxo=0.05 # offset information: number,
# first, spacing

susynlv nt=$n1 ft=$ft dt=$dt \


nxo=$nxo fxo=$fxo dxo=$dxo \
nxs=$nxs fxs=$fxs dxs=$dxs \
v00=$v00 dvdx=$dvdx dvdz=$dvdz \
ob=$ob er=$er smooth=$smooth \
ref=$REF1 | suwind key=gx min=750 \
max=750 > dataSingleReceiver

Common Offset

nt=300 ft=0.0 dt=0.008 # time information: number,


# first, spacing
nxm=30 fxm=0.05 dxm=0.025 # midpoint information: number,
# first, spacing
nxo=1 fxo=0.0 dxo=0.05 # offset information: number,
# first, spacing
Common Midpoint

nt=300 ft=0.0 dt=0.008 # time information: number,


# first, spacing
nxm=1 fxm=0.05 dxm=0.025 # midpoint information: number,
# first, spacing
nxo=30 fxo=0.0 dxo=0.05 # offset information:
# number, first, spacing

Perlu diperhatikan disini bahwa pada script untuk common receiver


terdapat fungsi SU lain yang digunakan yaitu suwind, karena untuk
menghasilkan data common receiver diperlukan pembuatan sintetik data
secara multiple shot untuk kemudian di pilah ( windowing ) menjadi
common receiver ( untuk lebih jelasnya lihat Bag. 4. Windowing dan
Sorting ), Hasil Script diatas dapat dilihat pada gambar dibawah :
Common Common
receiver midpoint

Common
offset

Sedangkan untuk menghasilkan data multiple shot baik menggunakan


konfigurasi berdasarkan susunan source dan receiver, maupun
menggunakan konfigurasi berdasarkan susunan midpoint dan offset dapat
dilihat pada perubahan shell script seperti dibawah :

Susunan Source dan Receiver

REF2="0.5:-1.0,2.0;5.0,2.0"

n1=300 ft=0.0 dt=0.008 # time information: number,


# first, spacing
nxs=10 fxs=0.0 dxs=0.025 # shot information: number,
# first, spacing
nxo=30 fxo=0.05 dxo=0.05 # offset information: number,
# first, spacing

susynlv nt=$n1 ft=$ft dt=$dt \


nxo=$nxo fxo=$fxo dxo=$dxo \
nxs=$nxs fxs=$fxs dxs=$dxs \
v00=$v00 dvdx=$dvdx dvdz=$dvdz \
ob=$ob er=$er smooth=$smooth \
ref=$REF1 ref=$REF2 > datamultiple

Susunan Midpoint dan Offset

REF2="0.5:-1.0,2.0;5.0,2.0"

nt=300 ft=0.0 dt=0.008 # time information: number,


# first, spacing
nxm=10 fxm=0.05 dxm=0.025 # midpoint information: number,
# first, spacing
nxo=30 fxo=0.0 dxo=0.05 # offset information: number,
# first, spacing

susynlv nt=$n1 ft=$ft dt=$dt \


nxo=$nxo fxo=$fxo dxo=$dxo \
nxs=$nxs fxs=$fxs dxs=$dxs \
v00=$v00 dvdx=$dvdx dvdz=$dvdz \
ob=$ob er=$er smooth=$smooth \
ref=$REF1 ref=$REF2 > dataCMPCOPmultiple
( dengan menggunakan Matlab dapat juga dibuat data sintetik, untuk
lebih jelasnya lihat pada bagian peredaman noise dengan dekonvolusi,
pada chapter VI )

Untuk membuat data sintetik mengenai groundroll atau even-even


SU
tips yang bersifat linear dapat dilakukan di SU dengan menggunakan
rutin suaddevent .
CHAPTER 5
WINDOWING AND SORTING ( TRACE HEADER MANIPULATION IN

REVIEW )

Dalam pengolahan data seismik biasanya diperlukan pemilahan trace baik


itu dalam time atau dalam urutan trace. Dalam SU terdapat beberapa
routine yang dapat melakukan pemilahan trace data seismik, diantaranya
adalah suwind, susort, sushw, suchw, suaddhead, suedit dll. Yang akan
dibahas disini adalah suwind, susort, sushw dan suchw dan beberapa
routine pembatu untuk melihat rentang isi trace header. Sushw adalah
suatu routine untuk menetapkan isi suatu variabel dalam trace dengan nilai
yang baru sedangkan suchw menetapkan isi suatu nilai dari variabel
dengan mengacu pada nilai yang sudah ada baik itu dari variabel itu
sendiri atau dari variabel lain. Sedangkan suwind adalah suatu routine
untuk memilah trace berdasarkan nilai variabel dari trace tersebut dan
susort adalah routine untuk mengurutkan trace berdasarkan nilai suatu
variabel dalam trace ( baik kecil ke besar dan vice versa ).
Formulasi perhitungan yang terdapat pada sushw adalah :
i = itr + d
val (key) = a + b * ( i % j ) + c * ( i / j )
dimana % adalah operasi modulo dan / adalah operasi pembagian.
Formulasi perhitungan yang terdapat pada suchw adalah :
Val (key1) = ( a + b * val (key2) + c * val (key3) ) / d
Kembali pada contoh untuk menghasilkan data multiple diatas dimana
shell scriptnya adalah sbb :
#! /bin/sh

WIDTH1=300
WIDTH2=600
HEIGHT=400
WIDTHOFF1=0
WIDTHOFF2=315
HEIGHTOFF1=100
HEIGHTOFF2=550

######Work in km (synlv puts 10^3 in header fields)######

REF1="1.0:0.0,1.0;5.0,1.0"
REF2="0.5:-1.0,2.0;5.0,2.0"

n1=300 ft=0.0 dt=0.008 # time information: number, first,


spacing
nxs=50 fxs=0.0 dxs=0.05 # shot information: number,
first, spacing
nxo=10 fxo=0.05 dxo=0.05 # offset information: number,
first, spacing

dvdz=0.0 # velocity gradient (dv/dz)


dvdx=0.0 # velocity gradient (dv/dx)
v00=2.0 # P-wave velocity at surface in km/s
ob=1 # include obliquity (both sides of interface are
# reflecting)
er=1 # use exploding reflector amps
smooth=0 # connect picks on reflectors linearly

susynlv nt=$n1 ft=$ft dt=$dt \


nxo=$nxo fxo=$fxo dxo=$dxo \
nxs=$nxs fxs=$fxs dxs=$dxs \
v00=$v00 dvdx=$dvdx dvdz=$dvdz \
ob=$ob er=$er smooth=$smooth \
ref=$REF1 ref=$REF2 > datamultiple
suximage < datamultiple perc=99 title="Multiple Shot \
windowtitle="Multiple Shot"\
label1="Time (sec)" label2="Shot Gather" \
f1=0 d1=0.08 f2=0 d2=0.05 \
-geometry ${WIDTH2}x${HEIGHT}+${WIDTHOFF2}+${HEIGHTOFF1}&

suchart < datamultiple | xgraph n=500 linecolor=1 \


linewidth=0 marksize=3 mark=8 &

dari shell script diatas terlihat bahwa data sintetik dibuat dengan
pengaturan adanya 50 shot dengan masing-masing shot terdapat 10
receiver dengan penambahan jarak shot dan offset receiver sebesar 0.05
km atau 50 m. hasil dapat dilihat pada gambar sbb :
misal saat ini ingin dapat melihat data dimana tidak ada data yang
overlaping atau model susunan dimana tiap shot mempunyai rentang jarak
sebesar susunan receiver ( lihat gambar )

maka dapat digunakan routine suwind dengan menggunakan nilai


parameter :
j = melewatkan trace j sebesar ( rentang jarak )
suwind < datamultiple > datawindowmul0 key=sx s=0 j=500
apabila dilihat data seismiknya maka menjadi :
Sebelum kita melakukan proses selanjutnya, ada baiknya kita melihat dulu
susunan nilai parameter pada trace header, hal ini dilakukan dengan
mengunakan routine surange pada datamultiple dan datamultiple0 hasil
windowing diatas :
surange < datamultiple

#original
500 traces:
tracl=(1,500) tracr=(1,500) fldr=(1,50) tracf=(1,10)
trid=1 offset=(50,500) sx=(0,2450) gx=(50,2950) counit=1
ns=300 dt=8000 d2=0.050000 f2=0.050000

surange < datawindowmul0

#window multiple
50 traces:
tracl=(1,410) tracr=(1,410) fldr=(1,41) tracf=(1,10)
trid=1 offset=(50,500) sx=(0,2000) gx=(50,2500) counit=1
ns=300 dt=8000 d2=0.050000 f2=0.050000

pada datamultiple atau data ariginal terlihat bahwa terdapat 500 trace
dalam data dan nilai tracl ( trace in line ) dan tracr ( trace in reel )
mempunyai nilai 1,2,3,4,5.500, untuk melihat nilai variabel ini dapat
kita gunakan sugethw seperti :
sugethw < datamultiple key=tracl,tracr,tracf | more
dimana hasilnya adalah :
tracl=1 tracr=1 tracf=1

tracl=2 tracr=2 tracf=2

tracl=3 tracr=3 tracf=3

tracl=4 tracr=4 tracf=4

tracl=5 tracr=5 tracf=5

tracl=6 tracr=6 tracf=6

tracl=7 tracr=7 tracf=7

tracl=8 tracr=8 tracf=8

tracl=9 tracr=9 tracf=9

tracl=10 tracr=10 tracf=10

tracl=11 tracr=11 tracf=1

tracl=12 tracr=12 tracf=2

tracl=13 tracr=13 tracf=3

tracl=14 tracr=14 tracf=4

tracl=15 tracr=15 tracf=5

tracl=16 tracr=16 tracf=6

tracl=17 tracr=17 tracf=7

tracl=18 tracr=18 tracf=8


. . .
. . .
tracl=500 tracr=500 tracf=10

dari tampilan diatas terlihat bahwa nilai tracl dan tracr berurutan dari 1
sampai 500 dan nilai tracf ( trace in field ) terdiri dari 1 sampai 10 ( ini
melambangkan dalam satu field atau grup atau satu shot terdapat 10
receiver ). Sekarang akan kita lihat nilai yang didapat dari datamultiple0 :
sugethw < datawindowmul0 key=tracl,tracr,tracf | more
dimana hasilnya adalah :
tracl=1 tracr=1 tracf=1

tracl=2 tracr=2 tracf=2

tracl=3 tracr=3 tracf=3

tracl=4 tracr=4 tracf=4


tracl=5 tracr=5 tracf=5

tracl=6 tracr=6 tracf=6

tracl=7 tracr=7 tracf=7

tracl=8 tracr=8 tracf=8

tracl=9 tracr=9 tracf=9

tracl=10 tracr=10 tracf=10

tracl=101 tracr=101 tracf=1

tracl=102 tracr=102 tracf=2

tracl=103 tracr=103 tracf=3

tracl=104 tracr=104 tracf=4

tracl=105 tracr=105 tracf=5

tracl=106 tracr=106 tracf=6

tracl=107 tracr=107 tracf=7

tracl=108 tracr=108 tracf=8


. . .
. . .
tracl=410 tracr=410 tracf=10

dari kenampakan diatas terlihat bahwa tracl dan tracr walaupun berjumlah
50 buah ( pada datamultiple0 ) tetapi mempunyai nilai yang tidak
berurutan, nilai ini merupakan nilai dari datamultiple dan karena kita
memotong sebagian datanya maka nilai tracr dan tracl juga terpotong.
Untuk suatu alasan tertentu mungkin kita ingin agar nilai tracl dan tracr
kembali menjadi nomor yang berurutan sehingga apabila datamultiple0
diatas memiliki 50 trace maka nilai tracr dan tracl akan berurutan dari 1
sampai 50. hal ini dapat dilakukan dengan menggunakan sushw yaitu :
sushw < datawindowmul0 > datawindowmul key=tracl a=1 b=1 \
c=10 j=10 d=0
seperti telah dituliskan sebelumnya bahwa sushw ( SU Set Header Word )
mempunyai formulasi perhitungan :
i = itr + d
val (key) = a + b * ( i % j ) + c * ( i / j )
dimana % adalah modulo dan / adalah pembagian, yang perlu diperhatikan
disini adalah bahwa masing-masing variabel diatas adalah bilangan integer
atau bilangan bulat sehingga tidak ada bilangan pecahan dan setiap
pembagian akan menghasilkan bilangan bulat, contoh :
0/5=0
1/5=0
3/5=0
5/5=1
7/5=1
8/6=1
23 / 7 = 3
sementara modulo atau % adalah perhitungan dimana yang dicari adalah
nilai sisa dari suatu pembagian :
0 % 5 = 0 / 5 = 0 sisa 0
1 % 5 = 1 / 5 = 0 sisa 1
3 % 5 = 3 / 5 = 0 sisa 3
5 % 5 = 5 / 5 = 1 sisa 0
7 % 5 = 7 / 5 = 1 sisa 2
8 % 6 = 8 / 6 = 1 sisa 2
23 % 7 = 23 / 7 = 3 sisa 2
kembali ke formulasi diatas, misal kita mempunyai 50 trace yang akan
dihitung maka akan terdapat 50 iterasi atau itr yang dimulai dari 0 sampai
49. penjabaran formulasi diatas dapat dilihat sebagai berikut :
i = itr + d
val (key) = a + b*(i%j) + c*(i/j)

Nilai awal Pertambahan Pertambahan


nilai data nilai grup
dalam grup

j adalah nilai jumlah data dalam satu grup ( atau trace dalam satu grup )
dalam kasus datamultiple0 diatas j bernilai 10, kita juga ingin pertambahan
nilai dalam grup sebesar 10 jadi c bernilai 10, pertambahan nilai data
dalam suatu grup sebesar 1 sehingga b bernilai 10. Nilai awal adalah 1
sehingga a bernilai 1. val ( key ) adalah nilai yang akan kita hitung dan
ditaruh pada variabel yang kita tetapkan ( dalam hal ini tracr dan tracl )
Apabila dilihat maka akan menjadi

iterasi itr i Nilai dalam Nilai Grup Nilai Nilai total = Nilai
Grup, b=1, awal dalam grup + Nilai
j=10 grup + Nilai awal
1 0 0 0%10=0 0/10=0 1 1
1*(0)=0 10*(0)=0
2 1 1 1%10=1 1/10=0 1 2
1*(1)=1 10*(0)=0
3 2 2 2%10=2 2/10=0 1 3
1*(2)=2 10*(0)=0
4 3 3 3%10=3 3/10=0 1 4
1*(3)=3 10*(0)=0
5 4 4 4%10=4 4/10=0 1 5
1*(4)=4 10*(0)=0
6 5 5 5%10=5 5/10=0 1 6
1*(5)=5 10*(0)=0
7 6 6 6%10=6 6/10=0 1 7
1*(6)=6 10*(0)=0
8 7 7 7%10=7 7/10=0 1 8
1*(7)=7 10*(0)=0
9 8 8 8%10=8 8/10=0 1 9
1*(8)=8 10*(0)=0
10 9 9 9%10=9 9/10=0 1 10
1*(9)=9 10*(0)=0
11 10 10 10%10=0 10/10=1 1 11
1*(0)=0 10*(1)=10
12 11 11 11%10=1 11/10=1 1 12
1*(1)=1 10*(1)=10
13 12 12 12%10=2 12/10=1 1 13
1*(2)=2 10*(1)=10
14 13 13 13%10=3 13/10=1 1 14
1*(3)=3 10*(1)=10
.. .. .. .. .. ..

50 49 49 49%10=9 49/10=4 1 50
1*(9)=9 10*(4)=40
Apabila kita lihat dengan surange kembali maka didapat :
Surange < datawindowmul

50 traces:
tracl=(1,50) tracr=(1,410) fldr=(1,41) tracf=(1,10)
trid=1 offset=(50,500) sx=(0,2000) gx=(50,2500) counit=1
ns=300 dt=8000 d2=0.050000 f2=0.050000

terlihat bahwa nilai tracl sudah berubah menjadi nilai 1 sampai 50, untuk
merubah nilai tracr dapat digunakan cara yang sama.
Dari hasil surange baik untuk datamultiple maupun untuk data multiple0
diatas terlihat bahwa trace header mempunyai informasi nilai posisi sx
(posisi sumber atau source) dan posisi gx (posisi penerima atau receiver),
dari kedua informasi ini kita dapat menghitung posisi CDP atau common
depth point ( atau CMP ( common mid point )). Hal ini dilakukan dengan
menggunakan routine suchw ( SU Compute Header Word ), formulasinya
adalah :
Val (key1) = ( a + b * val (key2) + c * val (key3) ) / d
Dimana kita mempunyai dua variabel yang diketahui nilainya yaitu sx dan
gx sebagai val2 dan val3 sehingga val1 sebagai CDP dapat dihitung. Nilai
cdp merupakan transformasi linear dari sx dan gx yaitu :
CDP = ( sx + gx ) / 2
Sehingga kita dapat meletakkan nilai a=0, b=1, c=1, key2=sx, key3=gx
dan d=2 serta key1=cdp :
suchw < datawindowmul > datawindowmul2 key1=cdp key2=gx \
key3=sx b=1 c=1 d=2
dengan menggunakan data sebelum di window :
suchw < datamultiple > datamul2 key1=cdp key2=gx key3=sx \
b=1 c=1 d=2
dan dengan surange dapat dilihat nilai CDP yang telah dihitung:
surange < datawindowmul2
tracl=(1,50) tracr=(1,410) fldr=(1,41) tracf=(1,10)
cdp=(25,2250) trid=1 offset=(50,500) sx=(0,2000)
gx=(50,2500) counit=1 ns=300 dt=8000 d2=0.050000
f2=0.050000

surange < datamul2


tracl=(1,500) tracr=(1,500) fldr=(1,50) tracf=(1,10)
cdp=(25,2250) trid=1 offset=(50,500) sx=(0,2000)
gx=(50,2500) counit=1 ns=300 dt=8000 d2=0.050000
f2=0.050000
Setelah perhitungan CDP dilakukan maka data seismik diatas dapat dilihat
dalam sudut pandang yang berbeda, apabila sebelumnya data seismik
dilihat dalam susunan shot gather ( atau mungkin receiver gather ) maka
saat ini kita dapat melihatnya dalam bentuk CDP gather ( dan juga offset
gather ). Pertama-tama yang harus dilakukan adalah melakukan sorting
pada data, yaitu mengurutkan nilai CDP dari yang terkecil sampai yang
terbesar dan hal ini dilakukan dengan menggunakan routine susort :

susort < datawindowmul2 > datasortmul cdp offset


suchart < datasortmul | xgraph n=50 linecolor=1 linewidth=0\
marksize=3 mark=8 &
suximage < datasortmul perc=99 title=" Sorting on windowed\
multiple shot" windowtitle="Window + Sort" label1="Time \
(sec)" label2="cdp gather"f1=0.0 d1=0.08 f2=0.025 d2=0.05 &

( script diatas akan melakukan sorting pada data yang telah di window dan
kemudian dilihat susunan posisi sumber dan penerima dalam bentuk CDP
gather )
susort < datamul2 > datasortmul3 cdp offset
suchart < datasortmul3 | xgraph n=500 linecolor=1 \
linewidth=0 marksize=3 mark=8 &
suximage < datasortmul3 perc=99 title=" Sorting on non \
windowed multiple shot" windowtitle="cdp Sort" label1="Time\
(sec)" label2="cdp gather" f1=0.0 d1=0.08 f2=0.025 d2=0.05 &

( script diatas akan melakukan sorting pada data yang belum di window
(data aslinya) dan kemudian dilihat susunan posisi sumber dan penerima
dalam bentuk CDP gather )
hasilnya dapat dilihat pada gambar berikut :
cdp gather

apabila dilihat dari hasil penggambaran sususan sumber dan penerima


pada grafik diatas terlihat bahwa tidak ada perbedaan antara cdp gather
dengan shot gather baik untuk yang telah di window maupun yang belum
(data asli). Untuk bagian yang telah di window hal ini terjadi karena
apabila dilihat dari susunan sumber dan penerima maka pada data ini
untuk satu CDP hanya terdapat satu dan hanya satu saja offset-nya ( full
fold hanya satu ) sehingga kenampakan susunan shot gather dan CDP
gather akan nampak sama ( tidak berubah ). Untuk data asli yang tersusun
sebagai CDP gather untuk dapat melihat dalam bentuk CDP dan offset
maka harus dilakukan windowing kembali untuk hanya mengambil sx dan
gx yang termasuk dalam bagian full fold CDP, dengan kata lain CDP yang
tidak mempunyai bagian dalam full fold tidak akan dipakai. Hal ini dapat
dilakukan dengan menggunakan script sebagai berikut :

suwind < datamul2 > datamulcdpfold key=cdp min=225 max=2475


susort < datamulcdpfold > datamulcdpfold2 cdp offset
suchart < datamulcdpfold2 | xgraph n=500 title="cdp gather \
full fold" label1="offset" label2="cdp" linecolor=1 \
linewidth=0 marksize=3 mark=8 &
suximage < datamulcdpfold2 perc=99 title="CDP full fold \
sorting on non windowed multiple shot" label1="Time (sec)" \
label2="CDP gather"
kenampakan grafik mempunyai sumbu aksis gx dan sx, tampilan diatas
dapat dirubah menjadi dalam bentuk aksis CDP ( CMP ) dan offset, untuk
lebih lengkapnya dapat dilihat pada script dibawah :
#! /bin/sh

WIDTH1=300
WIDTH2=600
HEIGHT=400
WIDTHOFF1=0
WIDTHOFF2=315
HEIGHTOFF1=100
HEIGHTOFF2=550

######Work in km (synlv puts 10^3 in header fields)######


# ref="1:1,2;4,2" reflector(s):"amplitude:x1,z1;x2,z2;..."
REF1="1.0:0.0,1.0;5.0,1.0"
REF2="0.5:-1.0,2.0;5.0,2.0"

n1=300 ft=0.0 dt=0.008 # time information: number,


# first, spacing
nxs=50 fxs=0.0 dxs=0.05 # shot information: number,
# first, spacing
nxo=10 fxo=0.05 dxo=0.05 # offset information: number,
# first, spacing

dvdz=0.0 # velocity gradient (dv/dz)


dvdx=0.0 # velocity gradient (dv/dx)
v00=2.0 # P-wave velocity at surface in km/s
ob=1 # include obliquity (both sides of interface are
# reflecting)
er=1 # use exploding reflector amps
smooth=0 # connect picks on reflectors linearly

susynlv nt=$n1 ft=$ft dt=$dt \


nxo=$nxo fxo=$fxo dxo=$dxo \
nxs=$nxs fxs=$fxs dxs=$dxs \
v00=$v00 dvdx=$dvdx dvdz=$dvdz \
ob=$ob er=$er smooth=$smooth \
ref=$REF1 ref=$REF2 > datamultiple
suximage <datamultiple perc=99 title="Multiple ShotRecords"\
windowtitle="Multiple Shot" label1="Time(sec)" \
label2="Shot Gather" f1=0 d1=0.08 f2=0 d2=0.05 \
-geometry ${WIDTH2}x${HEIGHT}+${WIDTHOFF2}+${HEIGHTOFF1} &
suchart < datamultiple key1=sx key2=gx | xgraph \
title="shot gather on raw data" label1="sx" label2="gx" \
n=500 linecolor=1 linewidth=0 marksize=3 mark=8 &

#-----------------------------------------------------------

suwind < datamultiple > datawindowmul0 key=sx s=0 j=500


suximage < datawindowmul0 perc=99 \
title="Windowing on Multiple Shots" \
windowtitle="Windowing" label1="Time (sec)"\
label2="Shot Gather" f1=0 d1=0.08 f2=0 d2=0.05 &
sushw < datawindowmul0 > datawindowmul key=tracl a=1 b=1 \
c=10 j=10 d=0 suchart < datawindowmul key1=sx key2=gx |
xgraph title=" shot gather on windowed data" \
label1="sx" label2="gx" n=50 linecolor=1 linewidth=0 \
marksize=3 mark=8 &

suchw < datawindowmul > datawindowmul2 key1=cdp key2=gx \


key3=sx b=1 c=1 d=2
susort < datawindowmul2 > datasortmul cdp offset
suchart < datasortmul key1=offset key2=cdp | xgraph \
title="offset gather windowed data" label1="offset" \
label2="cdp" n=50 linecolor=1 linewidth=0 marksize=3 \
mark=8 &
suximage < datasortmul perc=99 \
title=" CDP Sorting on windowed multiple shot" \
windowtitle="Window + Sort" label1="Time (sec)" label2="cdp\
gather" f1=0.0 d1=0.08 f2=0.025 d2=0.05 &

susort < datawindowmul2 > datasortmul2 offset cdp


suchart < datasortmul2 key1=offset key2=cdp | xgraph \
title="offset gather windowed data" label1="offset" \
label2="cdp" n=50 linecolor=1 linewidth=0 marksize=3 \
mark=8 &
suximage < datasortmul2 perc=99 \
title=" offset sorting on windowed multiple shot" \
windowtitle=" offset sorting" label1="time" label2="offset \
gather" f1=0 d1=0.08 f2=0.05 d2=0.05 &

#-----------------------------------------------------------

suchw < datamultiple > datamul2 key1=cdp key2=gx key3=sx \


b=1 c=1 d=2
susort < datamul2 > datasortmul3 cdp offset
suchart < datasortmul3 key1=offset key2=cdp | xgraph n=500 \
title="offset gather on raw data " label1="offset" \
label2="cdp" linecolor=1 linewidth=0 marksize=3 mark=8 &
suximage < datasortmul3 perc=99 \
title=" CDP Sorting on non windowed multiple shot" \
windowtitle="cdp Sort" label1="Time (sec)" \
label2="cdp gather" f1=0.0 d1=0.08 f2=0.025 d2=0.05 &

susort < datamul2 > datasortmul4 offset cdp


suchart < datasortmul4 key1=offset key2=cdp | xgraph n=500 \
title="offset gather on raw data" \
label1="offset" label2="cdp" linecolor=1 linewidth=0 \
marksize=3 mark=8 &
suximage < datasortmul4 perc=99 \
title="Offset sorting on non windowed multiple shot" \
windowtitle=" offset sorting" label1="time" \
label2="offset gather" f1=0 d1=0.08 f2=0.05 d2=0.05 &

#-----------------------------------------------------------

suwind < datamul2 > datamulcdpfold key=cdp min=225 max=2475


susort < datamulcdpfold > datamulcdpfold2 cdp offset
suchart < datamulcdpfold2 key1=offset key2=cdp | xgraph \
n=500 title="cdp/offset gather full fold" label1="offset" \
label2="cdp" linecolor=1 linewidth=0 marksize=3 mark=8 &
suximage < datamulcdpfold2 perc=99 \
title="CDP full fold sorting on non windowed multiple shot"\
label1="Time (sec)" label2="CDP gather"

exit 0
hasil dari suchart dari script diatas adalah sebagai berikut :
offset
CHAPTER 6
SEG-Y,SU,GPR ,OYO DATA , READ..AND SO THEY SAY

Setelah kita mengetahui susunan struktur data dalam SEG-Y dan SU maka
akan dengan mudah kita dapat membaca data tersebut dengan berbagai
program yang ada baik itu kita susun sendiri atau dengan program
komersial. Data lain juga dapat disusun dalam format SU asalkan format
data lain tersebut juga kita ketahui, sehingga dapat ditransformasikan
dalam format SU untuk selanjutnya di olah dengan SU.
Contoh untuk membaca data lain dan diolah dengan SU adalah data GPR,
misal terdapat data GPR format RAMAC dimana dalam format ini data
terdiri dari 2 file, yaitu file berekstensi *.rad dan *.rd3, file berekstensi
*.rad adalah file header yang berisi keterangan-keterangan mengenai data
GPR yang diambil, misal berapa Hz frekwensi yang digunakan, nilai
sample per trace ( ini penting !! ) atau ns dan lainnya, file ini ditulis dalam
format teks atau kode ASCII sehingga dapat dibaca dengan teks editor
biasa ( notepad, wordpad, nc di windows/dos. mc, kwrite di linux ).
Sedangkan file berekstensi *.rd3 adalah file data GPR itu sendiri dimana
data GPR ditulis dalam format 16 bit integer.
Untuk bisa diolah dan dibaca oleh SU tentunya yang perlu dirubah adalah
format penulisan data GPR yaitu ditransformasikan dari 16 bit integer ke
32 bit floating point dan kemudian ditambahkan header untuk tiap tracenya
dimana dalam tiap header ini ditambahkan keterangan jumlah sample per
tracenya ( ns ) yang diketahui dari file *.rad tadi.
File header Data GPR
Ekstensi *.rad Ekstensi *.rd3
Format 16 bit integer

16 bit integer 32 bit


floating point

Data GPR
Format 32 bit floating point

Informasi dari header file


+ header tiap trace
( number of sample / ns ( 240 byte )
, dll )

Data GPR
Format SU

Apabila berada dalam lingkungan Linux dan menggunakan SU maka


prosedur diatas dapat dengan mudah dilakukan dengan cara :
1. mengubah data GPR menjadi format 32 bit floating point
recast < fileinput.rd3 in=short out=float > fileout32bit
2. menambahkan header pada tiap trace
( nilai ns harus diketahui, dalam GPR format RAMAC nilai ns ini
tersimpan dalam file teks berekstensi *.rad )
suaddhead < fileout32bit ns=520 > fileout.su
Dalam lingkungan DOS/Windows dengan menggunakan Matlab hal diatas
dapat dilakukan dengan membuat suatu fungsi ( *.m file ) untuk
menjalankan prosedur diatas. Dibawah ini terdapat contoh m file untuk
menjalankan prosedur diatas dimana file input adalah data GPR baik
berformat little endian atau big endian ( biasanya little endian ) ,
sedangkan file output adalah file yang sudah berformat SU yang juga bisa
berformat little endian atau big endian ( biasanya big endian )

function gpr2su(ns,filenamein,filenameout,endian)
% fungsi untuk transformasi data GPR RAMAC menjadi SU format
% created by randy on may 16th 2003
%
% endian=1, input litle output litle
% endian=2, input litle output big
% endian=3, input big output litle
% endian=4, input big output big
% data type = IEEE 32 bits floating point

% open file
if ~exist(filenamein),
error('File tidak ada atau file harus berextension .rd3');
return;
end

if nargin == 4
if endian==1 | endian==2 % input file format litlle endian
FID = fopen(filenamein,'r','l');
elseif endian==3 | endian==4 % input file format big endian
FID = fopen(filenamein,'r','b');
end
elseif nargin ==3 % secara default dalam litlle endian
FID = fopen(filenamein,'r','l');
else
error('jumlah input tidak mencukupi ( tak ada file out ) atau terlalu
banyak input');
return;
end
% read the file to the end
% A adalah matriks berdimensi 1 x ( ntrace x ns )
% dimana ( ntrace x ns ) adalah jumlah kolom matriks A
[A] = fread(FID,inf,'short');
fclose(FID)
% karena ns sudah diketahui maka untuk mengubah matriks baris A
% harus diketahui jumlah jumlah trace ( ntrace ) dengan cara membagi
% jumlah kolom matriks A dengan ns
ntrace=size(A,1)./ns;
% memformulasikan header trace
% 240 bytes header trace dianggap format 2 bytes short
% ( 16 bit integer )
% karena hanya akan digunakan untuk menset ns ( number of
% samples )
% jadi 240 bytes / 2 bytes = 120 elemen
header=zeros(1,120);
% ns ditaruh pada bytes ke (114./2)+1=58
header(58)=ns;

% membuat file output


if nargin==4
if endian==1 % format litlle endian
FID = fopen(filenameout,'w','l');
elseif endian==2 % format big endian
FID = fopen(filenameout,'w','b');
elseif endian==3 % format litle endian
FID = fopen(filenameout,'w','l');
elseif endian==4 % format big endian
FID = fopen(filenameout,'w','b');
end
elseif nargin==3
FID = fopen(filenameout,'w','b'); % secara default format big endian
end

for i=1:ntrace
awal=((i-1).*ns)+1;
akhir=i.*ns;
data=A(awal:akhir);
fwrite(FID,header,'int16');
fwrite(FID,data,'float');
end

fclose(FID);

function OYO2segy
% fungsi untuk membaca file OYO dan mentransformasikan
% dalam format SEGY
% created by randy on 2003

% masukan nama file OYO yang akan dibaca


prompt = {'Masukan nama file oyo:','Nama file output:'};
lines= 1;
%def = {'*.org','*.sgy'};
def = {'1l1ps24.org','tes1.sgy'};
answer = inputdlg(prompt,'Input Nama File',lines,def);
% jika cancel ditekan
if size(answer,1)==0 & size(answer,2)==0
return
else
answer1=char(answer(1,1));

% opening file
fid=fopen(answer1,'r','l');
[message,errnum] = ferror(fid);

% if file can be opened


if errnum==0
% calculate the bytes of data
SIZE_MAKS=9999999999999999;
fseek(fid,0,'bof');
BytesData=0;
for i=1:SIZE_MAKS
if feof(fid)~=1
fread(fid,1,'int8');
BytesData=BytesData+1;
else
break
end
end
% karena perhitungan mulai dari indeks 0
BytesData=BytesData-1;

% we know that format data OYO for mcseis OYO-3 is 16 bit


% integers
% get the information of nt, ns and dt
fseek(fid,2875,'bof');
tempdt=fread(fid,4,'uchar');
tempdt=char(tempdt');
dt=eval(tempdt);

% read the power of tenth


fseek(fid,2882,'bof');
pow=0;
for i=1:10
zeroval=fread(fid,1,'uint8');
if zeroval==46
break
else
pow=pow+1;
end
end

dt=dt*(10.^pow);

% read the nt value from file


fseek(fid,2867,'bof');
temp=fread(fid,2,'uchar');
temp=char(temp');
nt=eval(temp);

% read the ns value from file


fseek(fid,2870,'bof');
temp=fread(fid,4,'uchar');
temp=char(temp');
ns=eval(temp);

traces=zeros(nt,ns);

% set the pointer at the end of file


fseek(fid,0,'eof');

% read each trace from the end backwards to the beginning


for i=1:nt
% read each sample in each trace
% set the pointer backwards
fseek(fid,-(ns.*2.*i),'eof');
% allocate buffer
buff=zeros(1,ns);
% read the file and set to the buffer
buff=fread(fid,ns,'int16');
% put to trace from last to beginning
trace((nt-(i-1)),:)=buff';
end

figure;wigb(trace');
% Apakah kenampakan sudah betul atau ingin di flip
button = questdlg('Apakah trace ingin di Flip?','Question of The
Day..','Yes','No','Help','No');
if strcmp(button,'Yes')
figure;wigb((flipud(trace))');
button2 = questdlg('Is it Like this?','Question of The
Month..','Yes','No','Help','No');
if strcmp(button2,'Yes')
trace=flipud(trace);
elseif strcmp(button2,'No')
trace=trace;
elseif strcmp(button2,'Help')
msgbox('Did i told you to learn?? , i believe i did..','..just
joking guys');
trace=trace;
end

elseif strcmp(button,'No')
trace=trace;
elseif strcmp(button,'Help')
msgbox('Sorry, You must learn by your self','Help','error');
trace=trace;
end

% writing header trace


% 240 bytes , 2 bytes each element = 120 element
header=zeros(1,120);
% the add 1 operation is included because of
% matlab is start with index 1
% ns ditaruh pada bytes ke (114./2)+1=58
header(58)=ns;
header(59)=dt;

% open file for write


answer2=char(answer(2,1));
button3 = questdlg('Little Endian file format?','Question of The
Year..','Yes','No','Help','Yes');
if strcmp(button3,'Yes')
FID = fopen(answer2,'w','l');
elseif strcmp(button3,'No')
FID = fopen(answer2,'w','b');
elseif strcmp(button3,'Help')
msgbox('well ask me instead??..','Whats on earth is this');
FID = fopen(answer2,'w','b');
end

% write trace header


textHeader=zeros(1,3200);
fseek(FID,0,'bof');
fwrite(FID,textHeader,'uint8');
% write binary header
binHeader=zeros(1,400);
%fseek(FID,0,'eof');
fwrite(FID,binHeader,'uint8');

% read output type

answer3 = inputdlg('Tipe data keluaran? 1. Float 4 bytes 2. Fixed


4 bytes 3. Fixed 2 bytes','Output File Format');

switch eval(char(answer3))
case 1
type=1;
fseek(FID,0,'bof');
fseek(FID,3224,'bof');
fwrite(FID,type,'uint8');
% write trace to file
fseek(FID,0,'eof');
for i=1:nt
fwrite(FID,header,'uint16');
data=trace(i,:);
fwrite(FID,data,'float32');
end
case 2
type=2;
fseek(FID,0,'bof');
fseek(FID,3224,'bof');
fwrite(FID,type,'uint8');
% write trace to file
fseek(FID,0,'eof');
for i=1:nt
fwrite(FID,header,'uint16');
data=trace(i,:);
fwrite(FID,data,'int32');
end
case 3
type=3;
fseek(FID,0,'bof');
fseek(FID,3224,'bof');
fwrite(FID,type,'uint8');
% write trace to file
fseek(FID,0,'eof');
for i=1:nt
fwrite(FID,header,'uint16');
data=trace(i,:);
fwrite(FID,data,'int16');
end
end

% write the information of data


% you do need the add 1 operation here
fseek(FID,0,'bof');
fseek(FID,3220,'bof');
fwrite(FID,ns,'uint16');
fseek(FID,0,'bof');
fseek(FID,3216,'bof');
fwrite(FID,dt,'uint16');

fclose(fid);
fclose(FID);

% jika file tak dapat dibuka


else
errordlg('File Tidak ada atau tidak dapat dibuka');
end
end
CHAPTER 7
NOISE , THE PERFECT UNWANTED

Dalam SU terdapat berbagai macam cara untuk berperang dengan Noise.


Dengan menggunakan sufilter maka berbagai macam filtering yang
berkaitan dengan frekwensi dapat dilakukan, baik itu low pass, high pass
maupun band pass untuk menghilangkan ground roll atau mungkin guided
waves. Dengan menggunakan sufxdecon atau supef maka kita dapat
menekan efek dari multiple yang ada pada data dan meningkatkan resolusi
dari sinyal atau data.
Apabila data diambil dalam bentuk gather dimana satu shot diterima oleh
banyak receiver ( seperti layaknya data seismik ) maka terdapat satu
prosedur dimana tingkat rasio sinyal dapat ditingkatkan. Prosedur ini
adalah stacking dan dalam SU dapat dilakukan dengan menggunakan
sustack ( tetapi rutin ini biasanya dilakukan setelah NMO dan analisa
kecepatan, untuk GPR terutama saat metoda pengambilan dengan profiling
atau zero offset biasanya proses stacking juga dilakukan saat pengambilan
data dilakukan ). Dengan menerapkan stacking maka signal to noise ratio
dapat ditingkatkan dan random noise dapat ditekan.
Pada contoh dibawah akan ditampilkan proses dekonvolusi yang
diterapkan pada data 1 dari kumpulan data Ozdogan Yilmaz ( data dapat di
download pada website cwp.mines.edu ). Proses dekonvolusi yang
diterapkan dibagi menjadi 2 bagian, yaitu bagian pertama adalah proses
dekonvolusi yang digunakan untuk menekan adanya multiple atau
reverberasi, sedangkan proses yang kedua adalah untuk meningkatkan
resolusi dengan menerapkan spiking deconvolution (lihat Appendix 3).
Proses ini kemudian dilanjutkan dengan menerapkan band pass filter
untuk menghilangkan efek ground roll yang biasanya berfrekwensi
rendah.
Shell script untuk melakukan routine diatas adalah sebagai berikut :
#! /bin/sh

# fungsi dekonvolusi pada data ozdata.1


# first look at the data
suxwigb < ozdata.1 perc=90 title="Data Ozdoglan Yilmas #1" \
label1="time (sec)" label2="km" &
# then look the autocorrelation of the data
suacor < ozdata.1 | suxwigb \
title="Autocorrelation of the data" &
# we can see that wavelet packet around 0.18 sec
# so first we attack the multiple
supef < ozdata.1 maxlag=0.18 | suxwigb perc=90 \
title="Multiple supression" &
# then we spike it up
supef < ozdata.1 maxlag=0.18 | supef minlag=0.01 \
maxlag=0.04 | suxwigb perc=90 title="multiple & spiking" &
# and band pass filtering
supef < ozdata.1 maxlag=0.18 | supef minlag=0.01 \
maxlag=0.04| sufilter f=5,20,30,35 | suxwigb perc=90 \
title="Multiple+spiking+band pass" &

hasil dari shell script diatas adalah :


Dengan menggunakan cara yang sama, data GPR juga dapat di olah seperti
halnya data seismik. Pada contoh dibawah adalah hasil proses dekonvolusi
dari data sample yang disertakan dalam program Gradix ( tentunya
sebelum dilakukan pemrosesan data perlu ditranformasikan dulu format
data GPR ke dalam format data SU, lihat bagian V untuk file gpr2su.m
atau dalam lingkungan SU dengan recast dan suaddhead ).

Data GPR
Reverberasi sudah
berkurang dan resolusi
meningkat

Data GPR setelah dekonvolusi

Dalam contoh berikut akan dibuat data sintetik dengan menggunakan


Matlab dimana fungsi traces.m ini diadopsi dari fungsi traces.c dari
su\src\demos\deconvolution\wiener-levinson dan dalam fungsi traces.m ini
ditambahkan proses dekonvolusi ( lihat appendix 3 ).
function traces

% fungsi yang diadopsi dari traces.c


% $CWPROOT\src\demos\deconvolution\wiener-levinson

x=[]; % trace value


ns=512; % number of samples
ntr=8; % number of traces
dtime=16; % one way time in samples through "ocean" layer
rbot=0.8; % reflection strength of "ocean" bottom
h=100; % location in samples of two way reverb train
amp=0.2; % strength of reflector
loc=170; % location of reflector on trace 1 in time samples
dip=12; % dip of reflector in time samples
x=zeros(ntr,ns);
for j=1:ntr
for i=1:ns
% mencari modulo dari (i-h) % dtime
sisa=((i-h)./dtime)-((floor((i-h)./dtime)));
if i >= h & sisa == 0
k = (i-h)/dtime;
if floor(k./2)==k./2
sgn=1;
else
sgn=-1;
end
x(j,i) = sgn .* (k+1) .* (rbot.^k);
end
% reflektor miring
if i == (loc + (j.*dip))
x(j,i) = x(j,i) + amp;
end
end
end
% jika ingin menambah noise maka tanda komentar dibawah
% dapat dihapus
% add noise
% noise=randn(ntr,ns).*0.2;
% x=x+noise;
x=x';
% doing convolution with source wavelet
% ingat minimum phase wavelet yang digunakan dalam konvolusi
% hal ini penting dalam asumsi dekonvolusi ( baca Yilmaz, seismic data
% proccesing , lihat juga appendix 3 )
w=ricker(40,0.004);w=w(19:37);
x=x';
xorigin=x;
for i=1:ntr
temp=x(i,:);
xkonv=konv(w,temp);
xnew(i,:)=xkonv;
end
x=xnew';
figure;wigb(x)

% doing autocorrelation

for i=1:ntr
temp=autocor(x(:,i));
acor(:,i)=temp';
end
figure;wigb(acor)

% terlihat adanya multiple wavelet setelah 35 ms


% now go after reverberation
% kita hanya mengambil sampai n ke 35
acor=acor';
for i=1:ntr
a_windowed=acor(i,1:35);
% crosscorelation of input dengan desired output ( spike )
spike=zeros(1,size(a_windowed,2));
spike(1)=1;
c=kroscor(a_windowed,spike);
% membentuk matriks Toeplitz
r=toeplitz(a_windowed);
filt=lsqr(r,c',0.0000001,100);
data=konv(filt,xorigin(i,:));
newx(i,:)=data;
end
newx=newx';
figure;wigb(newx);

Hasil dari fungsi diatas adalah sebagai berikut :


Data sintetik dimana terdapat Hasil autokorelasi data
banyak reverberasi atau
multiple yang mengubur
reflektor sebenarnya

Hasil proses dekonvolusi dari data sintetik dimana


terlihat reflektor miring dan reflektor horisontal
yang tadinya terkubur oleh reverberasi dan mungkin
noise
CHAPTER 8
VELOCITY ANALYSIS,NORMAL MOVEOUT AND VICE VERSA

Setelah melakukan sorting, windowing, filtering dan deconvolution


langkah selanjutnya adalah melakukan analisa kecepatan dari data seismik
( data lainnya ) yang ada. Analisa kecepatan ini dilakukan berulang-ulang
hingga tercapai bentuk model kecepatan yang ( menurut sang pemroses
data ) paling ideal atau mendekati keadaan sebenarnya. Terdapat beberapa
macam bentuk analisa kecepatan diantaranya adalah :
1. Velocity analysis
2. Velocity Scan ( NMO Scan )
3. Constant Velocity Stack
Pada gambar dibawah terdapat contoh sederhana dimana terdapat sebuah
lapisan horisontal, sebuah sumber S dan sebuah penerima P, dimana pada
jarak offset x dan sebuah titik tengah atau midpoint M akan menghasilkan
waktu tempuh menurut persamaan :
t2(x) = t2(o) + x2 / v2
t(o) = 2 MD / v
persaman diatas adalah persamaan waktu tempuh yang tergantung pada
offset dan kecepatan media yang dilalui oleh gelombang dimana hasilnya
akan membentuk suatu hiperbola.
x

S M P

D
Geometri Normal Moveout
t
t2(x) = t2(o) + x2 / v2

t = t2(x) ( t2(o) + x2 / v2 )

t2(o)

Travel Time Hyperbola

Dalam SU velocity analysis dapat dilakukan dengan menggunakan routine


suvelan, ingat bahwa terdapat hal penting yang harus diingat dalam
melakukan suvelan yaitu : data merupakan data CMP gather , apabila
tidak ada keterangan mengenai CMP atau nilai CMP tidak sama pada trace
header maka suvelan akan menghasilkan data kecepatan yang sukar
diinterpretasi ( it might still gonna work, but with bizarre result ! ).
Dibawah ini terdapat contoh data sintetik yang dibuat untuk kemudian di
sorting dan dipilih bagian CMP yang full fold dan kemudian dilakukan
analisa kecepatan.
#! /bin/sh

WIDTH1=300
WIDTH2=600
HEIGHT=400
WIDTHOFF1=0
WIDTHOFF2=315
HEIGHTOFF1=100
HEIGHTOFF2=550

# reflector(s):"amplitude:x1,z1;x2,z2;x3,z3;..."
REF1="1.0:0.0,1.0;5.0,1.0"
REF2="0.5:-1.0,2.0;5.0,2.0"

n1=300 ft=0.0 dt=0.008


nxs=50 fxs=0.0 dxs=0.05
nxo=10 fxo=0.05 dxo=0.05

dvdz=0.5
dvdx=0.0
v00=2.0
ob=1
er=1
smooth=0

susynlv nt=$n1 ft=$ft dt=$dt \


nxo=$nxo fxo=$fxo dxo=$dxo \
nxs=$nxs fxs=$fxs dxs=$dxs \
v00=$v00 dvdx=$dvdx dvdz=$dvdz \
ob=$ob er=$er smooth=$smooth \
ref=$REF1 ref=$REF2 > datamultiple
suximage < datamultiple perc=99 title="Multiple Shot
Records" windowtitle="Multiple Shot"\
label1="Time (sec)" label2="Shot Gather" \
f1=0 d1=0.08 f2=0 d2=0.05 \
-geometry ${WIDTH2}x${HEIGHT}+${WIDTHOFF2}+${HEIGHTOFF1} &

suchart < datamultiple | xgraph n=500 linecolor=1


linewidth=0 marksize=3 mark=8 &
#-----------------------------------------------------------

#-----------------------------------------------------------

suchw < datamultiple > datamul2 key1=cdp key2=gx \


key3=sx b=1 c=1 d=2
susort < datamul2 > datasortmul3 cdp offset
suchart < datasortmul3 | xgraph n=500 \
title="cdp gather" linecolor=1 \
linewidth=0 marksize=3 mark=8 &
suximage < datasortmul3 perc=99 \
title=" CDP Sorting on non windowed multiple shot"\
windowtitle="cdp Sort"\
label1="Time (sec)" label2="cdp gather" \
f1=0.0 d1=0.08 f2=0.025 d2=0.05 &

#-----------------------------------------------------------
suwind < datamul2 > datamulcdpfold key=cdp min=225 max=2475
susort < datamulcdpfold > datamulcdpfold2 cdp offset
suchart < datamulcdpfold2 key1=offset key2=cdp | xgraph \
n=500 title="cdp gather full fold" \
linecolor=1 linewidth=0 marksize=3 mark=8 &
suximage < datamulcdpfold2 perc=99 title="CDP full fold
sorting multiple shot" \
f2=1 d2=1 label1="Time (sec)" label2="CDP gather" &

exit 0

kemudian analisa kecepatan dilakukan dengan menggunakan script


sebagai berikut :
#! /bin/sh

WIDTH=450
HEIGHT=600
WIDTHOFF=40
HEIGHTOFF=50

fv=1000 # first velocity were going to use


nv=100 # jumlah data kecepatan
dv=25 # beda antara kec satu dengan lainnya
# 1000,1025,1050,1075 dst

fold=10
nout=501
dxout=0.004

# velan on cdp 225


cdp=225
suwind < datamulcdpfold2 key=cdp min=$cdp max=$cdp |
suvelan nv=$nv dv=$dv fv=$fv >veldata.$cdp
suximage < veldata.$cdp label1="time (sec)" \
label2="Velocity (m/sec)" \
title="Velocity analysis for CMP $cdp" cmap=hsv5 \
windowtitle="V-Scan= $cdp" legend=1 units=semblance \
wbox=$WIDTH hbox=$HEIGHT xbox=$WIDTHOFF \
ybox=$HEIGHTOFF f2=1000 d2=25 &

# velan on cdp 1200


cdp=1200
suwind < datamulcdpfold2 key=cdp min=$cdp max=$cdp |
suvelan nv=$nv dv=$dv fv=$fv >veldata.$cdp
suximage < veldata.$cdp label1="time (sec)" \
label2="Velocity (m/sec)" \
title="Velocity analysis for CMP $cdp" cmap=hsv5 \
windowtitle="V-Scan= $cdp" legend=1 units=semblance \
wbox=$WIDTH hbox=$HEIGHT xbox=$WIDTHOFF \
ybox=$HEIGHTOFF f2=1000 d2=25 &

# velan on cdp 2400


cdp=2400
suwind < datamulcdpfold2 key=cdp min=$cdp max=$cdp |
suvelan nv=$nv dv=$dv fv=$fv >veldata.$cdp
suximage < veldata.$cdp label1="time (sec)" \
label2="Velocity (m/sec)" \
title="Velocity analysis for CMP $cdp" cmap=hsv5 \
windowtitle="V-Scan= $cdp" legend=1 units=semblance \
wbox=$WIDTH hbox=$HEIGHT xbox=$WIDTHOFF \
ybox=$HEIGHTOFF f2=1000 d2=25 &

exit 0

hasil dari shell script diatas adalah :


Sedangkan apabila kita menggunakan data ozdogan Yilmaz ( contoh
dalam hal ini ozdata.10 ) maka perlu dilakukan perubahan nilai CMP
karena pada data tersebut bukan merupakan data CMP gather ( mungkin
data shot gather dll ). Hal ini dapat dilakukan dalam SU dengan perintah
suchw, pada contoh dibawah selain data CMP dihitung ulang, data
offset,sx dan gx juga dihitung kembali ( untuk mengetahui perubahan yang
terjadi maka dapa dilihat dengan perintah surange ) :

suchw < ozdata.10 > oz10.3 key1=cdp key2=sx b=1 d=1


suchw < oz10.3 > oz10.4 key1=offset key2=tracl b=50 d=1
surange < oz10.4
suchw < oz10.4 > oz10.3 key1=sx key2=cdp key3=offset \
b=1 c=0.5 d=1
surange < oz10.3
suchw < oz10.3 > oz10.4 key1=gx key2=cdp key3=offset \
b=1 c=-0.5 d=1

dengan menggunakan proses analisa kecepatan yang sama seperti pada


data sintetik diatas maka akan didapat seperti gambar berikut:
setelah melakukan analisa tentunya kemudian dapat ditentukan nilai
kecepatan pada rentang waktu tertentu pada data ( seismik ) tersebut, hal
ini dilakukan dengan memilih nilai maksimum pada data velan diatas
contoh :

Picking Velocity and Intercept Time

Pada data sintetik juga dapat kita pick nilai kecepatan dan nilai intercept
timenya sebagai berikut :
V RMS
V Hasil (Interpolasi)
Pick

Picking Velocity and Intercept Time On


Synthetic Data

sehingga didapat data masukan untuk NMO adalah :


cdp=225,1200,2400
tnmo=0,0.9,1.6
vnmo=2200,2200,2600
tnmo=0,0.9,1.6
vnmo=2200,2200,2600
tnmo=0,0.9,1.6
vnmo=2200,2200,2600
Format data diatas adalah pada CDP mana saja yang akan dijadikan
referensi, dalam hal ini terdapat tiga CDP yang dijadikan referensi yang
juga berarti terdapat tiga pasang intercept time (tnmo) dan velocity (vnmo)
dan juga jumlah data pada masing-masing tnmo dan vnmo haruslah sama.
Dengan berbekal ketiga referensi diatas maka sisa intercept time dan CDP
yang lain akan mempunyai nilai vnmo yang merupakan hasil interpolasi
sehingga menghasilkan RMS Velocity ( Root Mean Squared Velocity ).
Data diatas disimpan dalam suatu file yang kemudian akan diakses oleh
routine sunmo untuk melakukan koreksi normal moveout pada file input.
Dalam contoh diatas file input data kecepatan tersebut dinamakan stkvel.p,
untuk melakukan NMO dapat dilakukan dengan membuat shell script
seperti dibawah ini :
#! /bin/sh
# do and display nmo and stack

WIDTH=300
WIDTHOFF3=630
WIDTHOFF4=945
HEIGHT=400
HEIGHTOFF=50

# Do NMO
sunmo < datamulcdpfold2 par=stkvel.p > nmodata

# Display NMO
suximage < nmodata label1="Time" label2="CMP" title="NMO \
Data" f2=1 d2=1 &
# Sort to CDP's and stack
susort < nmodata |
sustack normpow=1.0 > stackdata

# Display Stack (cdps from 350-4600, spacing is 50, full


fold from 1450-3500)
sugain < stackdata tpow=2 gpow=0.5 |
suximage label1="Time" label2="CMP" title="Stack Data" \
windowtitle="Stack" f2=1 d2=1 legend=1 units="amplitude" \
wbox=$WIDTH hbox=$HEIGHT xbox=$WIDTHOFF4 ybox=$HEIGHTOFF &

exit 0
Pada shell script diatas selain dilakukan NMO kemudian berlanjut dengan
proses stacking dengan menggunakan routine sustack (perhatikan bahwa
setelah stacking maka jumlah trace akan berkurang tetapi S/N Ratio
meningkat) sehingga didapat hasil seperti berikut :

NMO corected Data


Stacked Data
Appendix 1 ( Setunggal )
Byte order, a review

Pengetahuan mengenai alamat byte suatu data dapat membantu dalam mengambil
informasi yang terkandung dalam alamat tersebut, hal ini juga akan membantu kita
untuk tidak tergantung pada satu perangkat lunak tertentu untuk membacanya (apalagi
apabila perangkat lunak itu harus beli !..), apabila struktur alamat suatu data telah kita
ketahui maka kita bisa saja menggunakan bermacam perangkat lunak yang telah ada
tergantung kesenangan atau membuat suatu rutin sendiri dari bahasa pemrograman
yang kita kuasai.
Informasi yang terkandung dalam suatu data SEG-Y terkandung dalam binary header
nya dan juga terdapat dalam masing-masing trace header, disini yang akan
disinggung adalah trace headernya karena berdasarkan informasi dari trace header
inilah sebagian besar program SU akan berjalan.
Dibawah ini terdapat alamat/indeks byte dari tiap-tiap variabel dalam trace header :

Alamat/indeks
Byte

.tracl = 'int'; % 32 bits = 4 bytes 0 (alamat byte yang terpakai adalah


0,1,2,3)
.tracr= 'int'; % 32 bits = 4 bytes 4 (alamat byte = 4 5 6 7 byte ke-8
adalah permulaan untuk .fldr)
.fldr= 'int'; % 32 bits = 4 bytes 8
.tracf= 'int'; % 32 bits = 4 bytes 12
.ep= 'int'; % 32 bits = 4 bytes 16
.cdp= 'int'; % 32 bits = 4 bytes 20
.cdpt= 'int'; % 32 bits = 4 bytes 24
.trid= 'short'; % 16 bits = 2 bytes 28
.nvs= 'short'; % 16 bits = 2 bytes 30
.nhs= 'short'; % 16 bits = 2 bytes 32
.duse= 'short'; % 16 bits = 2 bytes 34
.offset= 'int'; % 32 bits = 4 bytes 38
.gelev= 'int'; % 32 bits = 4 bytes 42
.selev= 'int'; % 32 bits = 4 bytes 46
.sdepth= 'int'; % 32 bits = 4 bytes 50
.gdel= 'int'; % 32 bits = 4 bytes 54
.sdel= 'int'; % 32 bits = 4 bytes 58
.swdep= 'int'; % 32 bits = 4 bytes 62
.gwdep= 'int'; % 32 bits = 4 bytes 66
.scalel= 'short'; % 16 bits = 2 bytes 68
.scalco= 'short'; % 16 bits = 2 bytes 70
.sx= 'int'; % 32 bits = 4 bytes 74
.sy= 'int'; % 32 bits = 4 bytes 78
.gx= 'int'; % 32 bits = 4 bytes 82
.gy= 'int'; % 32 bits = 4 bytes 86
.counit= 'short'; % 16 bits = 2 bytes 88
.wevel= 'short'; % 16 bits = 2 bytes 90
.swevel= 'short'; % 16 bits = 2 bytes 92
.sut= 'short'; % 16 bits = 2 bytes 94
.gut= 'short'; % 16 bits = 2 bytes 96
.sstat= 'short'; % 16 bits = 2 bytes 98
.gstat= 'short'; % 16 bits = 2 bytes 100
.tstat= 'short'; % 16 bits = 2 bytes 102
.laga= 'short'; % 16 bits = 2 bytes 104
.lagb= 'short'; % 16 bits = 2 bytes 106
.delrt= 'short'; % 16 bits = 2 bytes 108
.muts= 'short'; % 16 bits = 2 bytes 110
.mute= 'short'; % 16 bits = 2 bytes 112
.ns= 'unsigned short'; % 16 bits =2 bytes 114
.dt= 'unsigned short'; % 16 bits =2 bytes 116
.gain= 'short'; % 16 bits = 2 bytes 118
.igc= 'short'; % 16 bits = 2 bytes 120
.igi= 'short'; % 16 bits = 2 bytes 122
.corr= 'short'; % 16 bits = 2 bytes 124
.sfs= 'short'; % 16 bits = 2 bytes 126
.sfe= 'short'; % 16 bits = 2 bytes 128
.slen= 'short'; % 16 bits = 2 bytes 130
.styp= 'short'; % 16 bits = 2 bytes 132
.stas= 'short'; % 16 bits = 2 bytes 134
.stae= 'short'; % 16 bits = 2 bytes 136
.tatyp= 'short'; % 16 bits = 2 bytes 138
.afilf= 'short'; % 16 bits = 2 bytes 140
.afils= 'short'; % 16 bits = 2 bytes 142
.nofilf= 'short'; % 16 bits = 2 bytes 144
.nofils= 'short'; % 16 bits = 2 bytes 146
.lcf= 'short'; % 16 bits = 2 bytes 148
.hcf= 'short'; % 16 bits = 2 bytes 150
.lcs= 'short'; % 16 bits = 2 bytes 152
.hcs= 'short'; % 16 bits = 2 bytes 154
.year= 'short'; % 16 bits = 2 bytes 156
.day= 'short'; % 16 bits = 2 bytes 158
.hour= 'short'; % 16 bits = 2 bytes 160
.minute= 'short'; % 16 bits = 2 bytes 162
.sec= 'short'; % 16 bits = 2 bytes 164
.timbas= 'short'; % 16 bits = 2 bytes 166
.trwf= 'short'; % 16 bits = 2 bytes 168
.grnors= 'short'; % 16 bits = 2 bytes 170
.grnofr= 'short'; % 16 bits = 2 bytes 172
.grnlof= 'short'; % 16 bits = 2 bytes 174
.gaps= 'short'; % 16 bits = 2 bytes 176
.otrav= 'short'; % 16 bits = 2 bytes 178
.d1= 'float'; % 32 bits = 4 bytes 182
.f1= 'float'; % 32 bits = 4 bytes 186
.d2= 'float'; % 32 bits = 4 bytes 190
.f2= 'float'; % 32 bits = 4 bytes 194
.ungpow= 'float'; % 32 bits = 4 bytes 198
.unscale= 'float'; % 32 bits = 4 bytes 202
.ntr= 'int'; % 32 bits = 4 bytes 206
.mark= 'short'; % 16 bits = 2 bytes 208
.unass= 'short'; % 16 bits = 2 bytes 210 ( tidak teralokasikan / bebas )
.trace= 'float'; % 32 bits = 4 bytes 240

Pada struktur alamat byte diatas misalnya terlihat bahwa ns atau jumlah sampel per
trace berada pada byte ke-114 berarti untuk membaca informasi ns kita dapat melihat
pada byte ini. Untuk melihat isi byte ke-114 ini dapat digunakan bermacam cara
seperti :
1. Dengan Matlab dapat digunakan perintah
FID = fopen(ozdata.25,'r','b');
% membuka data ozdata.25 dalam mode read only dan data dibaca secara Big Endian
status = fseek(FID,114,'bof');
% meletakkan pointer pada byte ke-114 dari awal file (bof / beginning of file)
ns = fread(FID,1,ushort);
% membaca data sebanyak 1 elemen yang bertipe unsigned short
% sepanjang 2 byte dan dibaca sebagai integer / bilangan bulat.
( ozdata.25 dapat didownload di www.cwp.mines.edu, data ini sudah dalam format
SU sehingga 3200 byte teks data dan 400 byte binary header sudah tidak ada lagi,
apabila data masih dalam format SEG-Y maka alamat offset bytenya tinggal
ditambahkan, misal untuk alamat ns adalah 3200 + 400 + 114 = 3614 )
2. Dengan menggunakan KhexEdit di Linux kita dapat melihat langsung isi data
dalam bentuk biner ( binary code ), pada contoh gambar dibawah terlihat data
biner dari Ozdata.25 ,
Data dalam Bentuk
Ala bentuk biner ASCII
mat
byte

KhexEdit terbagi menjadi 3 bagian, bagian kiri adalah alamat byte ( byte ke- ) dari
data, bagian tengah adalah data dari file yang kita buka dalam bentuk biner, sedang
bagian kanan adalah terjemahan dalam bentuk ASCII nya. Pada gambar diatas terlihat
pada byte ke-114 sampai 115 ( ns mempunyai panjang data 2 byte,lihat kotak merah
pada gambar diatas ) terbaca 00001000 00110100 sehingga apabila diterjemahkan
menjadi 08 34 ( heksadesimal ) atau 2100 ( desimal ) sehingga ns = 2100 (terdapat
2100 data dalam tiap tracenya). Dapat dilihat juga bahwa sistem penulisan data diatas
menganut format Big Endian ( seperti kebanyakan data seismik ).
3. Dengan menggunakan NC (DOS/Windows) atau MC (Linux)
Bentuk ASCII
dari data

Data dalam
bentuk
heksadesimal

Alamat byte dalam


bentuk heksadesimal
Appendix 2 ( Kalih )
Susynlv, some of the mathematical concept
Keuntungan dari program SU dan program-program open source lainnya adalah
kebebasan yang diberikan oleh sang pembuat algoritma pemrograman untuk melihat
isi program tersebut sehingga kita dapat mengetahui cara program tersebut berjalan
atau dalam hal SU kita dapat mengetahui dasar teori matematika yang digunakan di
program tersebut. Fungsi susynlv sendiri dapat dilihat pada susynlv.c dan didalam
susynlv terdapat beberapa subroutine yang memanggil fungsi yang terdapat pada
modelling.c.
Pada prinsipnya program susynlv akan menghitung waktu tempuh antara sumber ke
reflektor dan reflektor ke penerima. Dalam fungsi susynlv kita diharuskan
memasukkan input berupa koordinat reflektor, dan didalam fungsi susynlv ini tiap
titik dalam koordinat ini akan dihubungkan dalam segmen-segmen kecil dengan
interpolasi dengan cubic spline method ( misal REF1=amp:a;b;c;d;e;f;g;h;i dengan
masing-masing titik terdiri dari koordinat (x,z), i.e a=(x1,y1) dan seterusnya, maka
dengan interpolasi ini akan dibentuk segmen-segmen kecil yang merupakan
interpolasi antara titik-titik a,b,c,d tadi )

d e Interpolasi
segmen

f
c

g
b h
a i

Setelah reflektor diinterpolasi menjadi segmen-segmen maka untuk membuat satu


trace dihitung waktu tempuh dari tiap-tiap segmen untuk kemudian dikonvolusikan
dengan wavelet sumber ( misal wavelet ricker ).
Dalam susynlv perhitungan waktu dilakukan dengan menggunakan persamaan :
V
t = ln
Vo
persamaan diatas didapat dari :
x = v t
x
=t
v
atau
x
t=
v
x
t =
v(x)
v
dimana dijabarkan bahwa v( x) = ( x ) = a ( x) , yaitu perubahan kecepatan ter-
x
hadap jarak , sehingga persamaan diatas akan menjadi :
x
t =
v( x) + a ( x)
1
t = v( x) + a( x)x
t + c = ln(V + a)
c = kons tan ta

pada saat t = 0, maka a = 0 dan V = Vo, sehingga


c = ln(Vo)
t + ln(Vo) = ln(V + a )
t = ln(V + a ) ln(Vo)
V + a
t = ln
Vo
apabila a sangat kecil dan bisa diabaikan akan didapat :
V
t = ln
Vo
persamaan diatas digunakan pada saat asumsi ray path berupa garis lurus, apabila ray
path tidak berupa garis lurus maka digunakan persamaan lain ( Ill leave it to you to
discover what kind of equation is that and how did it happen that way .. )
Appendix 3 ( Tigo )
Deconvolution, the Resolution Booster
Dalam lampiran ini akan dijabarkan secara sederhana penggunaan dekonvolusi untuk
meningkatkan resolusi dari data. Dasar teori yang lebih lengkap dapat dibaca pada
buku karangan Ozdogan Yilmaz, Seismic Data Processing, SEG.
Sebuah trace dalam seismogram dapat dijabarkan dari hasil konvolusi :
x(t) = w(t) * e(t) 3.1
dimana x(t) adalah trace data seismogram, w(t) adalah wavelet dan e(t) adalah
impulse response. Apabila ditambah dengan adanya noise maka menjadi :
x(t) = w(t) * e(t) + n(t) 3.2
Pada proses dekonvolusi terdapat beberapa asumsi yang digunakan sehingga proses
dekonvolusi dapat dilakukan yaitu :
1.a. Bumi terdiri dari lapisan-lapisan horisontal yang mempunyai kecepatan konstan.
b. Sumber getaran menghasilkan plane waves yang mengenai batas perlapisan pada
posisi tegak lurus terhadap perlapisan sehingga tidak terjadi gelombang geser
(shear waves).
2. Bentuk gelombang tidak berubah saat gelombang tersebut menjalar di dalam bumi
(prinsip Stationary).
3. Tidak terdapat Noise ( bising ), n(t) = 0 .
4. Bentuk gelombang sumber ( wavelet )diketahui .
5. Reflektifitas adalah proses yang acak sehingga seismogram mempunyai
karakteristik seperti wavelet pada spektrum amplitudo dan autokorelasinya.
Dengan berdasarkan pada asumsi diatas maka persamaan 3.2 dapat ditulis :
x(t) = w(t) * e(t) 3.3
sementara impulse resmponse dari bumi dapat disusun berdasarkan persamaan :
e(t) = a(t) * x(t) 3.4
a(t) adalah suatu operator inversi, sehingga dengan menggabungkan persamaan 3.3
dan 3.4 didapat :
x(t) = w(t) * a(t) * x(t) 3.5
dengan menghilangkan x(t) diatas maka :
(t) = w(t) * a(t) 3.6
(t) = 1 untuk t = 0 dan
(t) = 0 untuk t 0
dengan menyelesaikan persamaan 3.6 maka didapat filter operator untuk persamaan
3.4 yaitu :
a(t) = (t) * w(t) 3.7
dimana w(t) adalah invers dari seismic wavelet. Sehingga dengan mengetahui invers
wavelet maka didapatkan filter operatornya.
Karena wavelet dari suatu trace tidak diketahui, maka dapat disubtitusi dengan trace
itu sendiri, hal ini berdasar asumsi 5 dimana trace seismik mempunyai karakteristik
seperti wavelet sumbernya pada proses spektrum dan autokorelasinya ( lebih
lengkapnya baca Ozdogan Yilmaz, Seismic Data Processing ). Autokorelasi dari trace
seismik kemudian akan dibentuk menjadi suatu matriks , yaitu bentuk matriks
Toeplitz ( matriks yang mempunyai sifat simetri pada arah diagonal ) dan apabila
yang diinginkan adalah meningkatkan resolusi dimana wavelet menjadi spike maka
dapat dijabarkan menjadi bentuk matriks seperti :
r0 r1 r 2 ... rn 2 rn 1 a 0 g 0
r1 r0 r1 r 2 ... rn 2 a1 0

r2 r1 r 0 r1 r2 ... a 2 0
= 3.8
... r 2 r 1 r 0 r1 r 2 ... ...
rn 2 ... r 2 r1 r0 r1 an 2 ...

rn 1 rn 2 ... r 2 r1 r 0 an 1 0
dimana r adalah autokorelasi dari data seismik, a adalah operator filter yang dicari, g
adalah hasil crosscorelation antara input yang diinginkan (spike ,dll) dengan input
data seismik. Persamaan diatas adalah bentuk umum dari matriks persamaan A.x = b
yang dapat diselesaikan dengan metode least square ( bisa juga dengan metode
lainnya seperti SVD, conjugate gradient, atau rekursi Levinson untuk menghemat
memori, dll ).
Panjang filter operator a tergantung dari banyaknya data yang diambil dari
autokorelasi dari data seismik, semakin mirip karakteristik autokorelasi data seismik
dengan wavelet maka akan semakin baik filter operator yang didapat dan resolusi
yang dihasilkan semakin baik pula..
Pada contoh dibawah akan dijabarkan proses dekonvolusi dengan menggunakan
Matlab , disini digunakan wavelet ricker yang bersifat minimum phase ( sifat
minimum phase ini sangat penting! , baca ozdogan Yilmaz, Seismic Data Processing )
, fungsi ricker.m dapat didownload pada website M.D.Sacchi Dept. of Physics,
University of Alberta,Canada ( www.ualberta.ca ).
function decontest
% e = impulse response (dt = 0.004)
e=[10 0 0 0 0 0 0 10 0 10 0 0 10 0 0 0 0 0 0 5 10 0 0 0 0 0 0 0 0 0 0 0 0];
% w = wavelet ricker f dominan=40Hz dan dt=0.004
% diambil yang minimum phase
w=ricker(40,0.004);w=w(19:37);
% x adalah seismogram hasil konvolusi w dan e
x=konv(w,e);figure;plot(x)
% first kita cari autokorelasi dari seismogram
a=autocor(x);figure;plot(a);
for i=1:5
if i>1
pembagi=2*i;
else
pembagi=1;
end
% kita hanya mengambil sampai n ke
a_windowed=a(1:floor((size(a,2))./pembagi));
% crosscorelation of input dengan desired output ( spike )
spike=zeros(1,size(a_windowed,2));
spike(1)=1;
c=kroscor(a_windowed,spike);
% membentuk matriks Toeplitz
r=toeplitz(a_windowed);
filt=lsqr(r,c',0.0000001,100)
figure
e_new=konv(filt,x);
plot(e_new)
hold
plot(e,'r')
end

function r=toeplitz(in)
% membentuk matriks Toeplitz
kol=size(in,2);
r=zeros(kol,kol);
for i=1:kol
r(i:-1:1,i)=in(1:i)';
r(i:kol,i)=in(1:kol-i+1)';
end

function c=konv(a,b)
% Fngsi konvolusi antara a dan b
% matriks a yang dibalik dan bergerak dari kiri ke kanan
if size(a,1)>1
a=a';
end
if size(b,1)>1
b=b';
end
[b1,k1]=size(a);[b2,k2]=size(b);
b_new=zeros(k2,k2);
for i=1:k2
for j=1:k2
if i==j
b_new(i,j)=b(i);
end
end
end
a_new=zeros(2.*k1,k2);
for i=1:k2
a_new(i:i+(k1-1),i)=a(:);
end
c=sum(a_new*b_new,2)';

function c=autocor(a)
% Fungsi autokorelasi antara a dan a
b=a;
if size(a,1)>1
a=a';
end
if size(b,1)>1
b=b';
end
[b1,k1]=size(a);[b2,k2]=size(b);
b_new=zeros(k2,k2);
for i=1:k2
for j=1:k2
if i==j
b_new(i,j)=b(i);
end
end
end
a_new=zeros(k1,k1);
for i=1:k1
a_new(i:-1:1,i)=a(1:i)';
end
c=sum(a_new*b_new,2)';

function c=kroscor(a,b)
% Fungsi korelasi antara a dan b
% matriks a yang bergerak dari kiri ke kanan
if size(a,1)>1
a=a';
end
if size(b,1)>1
b=b';
end
[b1,k1]=size(a);[b2,k2]=size(b);
b_new=zeros(k2,k2);
for i=1:k2
for j=1:k2
if i==j
b_new(i,j)=b(i);
end
end
end
a_new=zeros(k2,k2);
for i=1:k2
a_new(i:-1:1,i)=a(1:i)';
end
c=sum(a_new*b_new,2)';

function [w] = ricker(f,dt)


%RICKER: Ricker wavelet
% This function designs a Ricker wavelet with
% peak frequence f.
% [w] = ricker(f,dt);
% IN f: central freq. in Hz (f <<1/(2dt))
% dt: sampling ineterval in sec
% OUT w: the ricker wavelet
% Example:
% w = ricker(40,0.004); plot(w);
% SeismicLab
% Version 1
% written by M.D.Sacchi, last modified December 10, 1998.
% [email protected]
% Copyright (C) 1998 Seismic Processing and Imaging Group
% Department of Physics
% The University of Alberta
nw=6./f/dt;
nw=2*floor(nw/2)+1;
nc=floor(nw/2);
i=1:nw;
alpha=(nc-i+1).*f*dt*pi;
beta=alpha.^2;
w=(1.-beta.*2).*exp(-beta);

Hasil dari fungsi diatas dapat dilihat pada gambar berkut :


Impulse Response

wavelet ricker(minimum phase)

seismogram ( konvolusi wavelet dan impulse response, perhatikan bahwa impulse


response pada titik 9 dan 10 terkubur menjadi satu bagian )
autocorelation of seismogram

Hasil Dekonvolusi dimana panjang operator = panjang autokorelasi


Merah=impulse response,biru=hasil dekonvolusi

Hasil Dekonvolusi dimana panjang operator = panjang autokorelasi / 2


Hasil Dekonvolusi dimana panjang operator = panjang autokorelasi / 4

Hasil Dekonvolusi dimana panjang operator = panjang autokorelasi / 6

Hasil Dekonvolusi dimana panjang operator = panjang autokorelasi / 8

Apabila wavelet ricker digunakan dalam proses autokorelasi, maka hasilnya akan
kembali seperti impulse response semula :
Penggunaan Wavelet Ricker dalam autokorelasi
Menghasilkan data seperti bentuk awal
Appendix 4 ( Sekawan )
Normal MoveOut, the coherence and semblance
Pada proses normal moveout diperlukan analisa kecepatan agar setiap hiperbola yang
ada dapat menjadi satu single event yang horisontal. Untuk melihat hal tersebut maka
diperlukan evaluasi coherence sepanjang kurva moveout. Analisa ini dilakukan pada
sekelompok nilai kecepatan tertentu sehingga kemudian suatu nilai kecepatan yang
berkaitan dengan puncak atau peak dari coherence inilah yang akan dapat
menghasilkan event yang flat / horisontal.
Jika C(to,v) melambangkan coherence pada sebuah window yang berpusat pada t = to ,
maka perkiraan kecepatan yang akan membuat kuva hiperbola menjadi flat ada pada
nilai peak dari C(to,v). perhitungan standar untuk pengukuran nilai coherence adalah
semblance, untuk menghasilkan sebuah semblance pertama-tama yang perlu
dilakukan adalah menentukan data window sepanjang 2M + 1 yang berpusat pada
suatu intercept time ti dan data tersebut telah di moveout dengan kecepatan vj sehingga
didapat ( Sacchi,1999 ) :

X ti Mt ,1 X ti Mt , 2 ... X ti Mt , N
M (ti , v j ) = ... ... ... ...
X ti + Mt ,1 X ti + Mt , 2 ... X ti + Mt , N
N adalah jumlah trace dan t adalah sampling waktu dari seismogram.
Semblance dapat didapat dengan menggunakan persamaan :

t (k X tk )
2

S (t i , v j ) =
X k t
2
tk

dimana t adalah indeks dari window data ( sejumlah 2M +1 ) dan k adalah indeks dari
jumlah trace ( sejumlah N ). Perhitungan diatas akan dihitung sebanyak jumlah
intercept time yang dipilih ( ti ) dan banyaknya kecepatan yang dipakai ( vj ).
Berikut adalah contoh data sintetik yang akan di analisa kecepatan dengan
membentuk semblance dengan 3 pilihan jumlah M, yaitu 0,5,10,15,.,50.
function sintetik
%
% fungsi untuk membuat data sintetik
% created by randy on 2003

% kedalaman = 1000 m
depth=1000;
% offset
offset=0:50:500;
% kecepatan konstan
velocity=1000;
% membuat array 2D dari data time
% dengan maksimal time=2 det
% jumlah data per trace = 1000
% sampling = dt = 2/1000 = 0.002 detik
ns=1000;
time_maks=2;
sampling=time_maks./ns;
trace=zeros(ns,size(offset,2));
% membuat wavelet ricker
w=ricker(20,sampling);

% menghitung jarak dari reflektor ke receiver


% kemudian menghitung waktu tempuh apabila kecepatan
% adalah konstan = 1500 m/det
for i=1:size(offset,2)
x=offset(i);
z=depth;
jarak=sqrt((x.^2)+(z.^2));
time=jarak./velocity;
indeks=floor(time./sampling);
temp=trace(:,i)';
temp(indeks)=1;
temp=konv(w,temp);
temp=temp(1:ns);
tracenew(:,i)=temp';
end
trace=tracenew;
imagesc(trace);colormap(gray)

Kenampakan data sintetik

function tesSemblance(trace,offset,sampling)
% fungsi untuk melihat kenampakan semblance
% untuk berbagai jumlah window (M)
% created by randy on 2003

% jika fungsi ini dijalankan setelah


% fungsi sintetik.m maka tidak perlu input
if nargin==0
load trace
offset=0:50:500;
ns=1000;
time_maks=2;
sampling=time_maks./ns;
end
% Looping sampai 10 macam jumlah M
for i=1:10
% menentukan jumlah M (window)
M=(i-1).*5;
% menentukan velocity moveout
v=1000;
% melakukan NMO
[trace_out] = cvnmo(trace,offset,sampling,v,1);
sem(:,i)=semblance(trace_out,M);
end
figure;
imagesc(sem)

function [D_out] = cvnmo(D_in,offset,dt,v,conj);


% CVNMO: constant velocity NMO and NMO^{-1} correction
%
% Apply constant velocity NMO and inverse NMO
% correction to CMP gatherts
%
% [D_out] = cvnmo(D_in,offset,dt,v,conj)
%
% IN D_in: data to NMO or un-NMO
% v: velocity in m/sec
% offset: source-receicer distance in meters
% dt: sampling in secs
% conj: flag to NMO conj=1
% conj: flag to un-NMO conj=-1
%
% OUT D_out: output data after NMO or un-NMO
%
%
%
% SeismicLab
% Version 1
%
% written by M.D.Sacchi, last modified December 10, 1998.
% [email protected]
%
% Copyright (C) 1998 Seismic Processing and Imaging Group
% Department of Physics
% The University of Alberta
[nt,nh]=size(D_in);
D_out=zeros(nt,nh);
n=1:1:nt;
for k=1:nh
d=sqrt( ((n-1)*dt).^2+(offset(k)/v)^2);
nn=floor(d/dt)+1;
if conj == 1;
for kk=1:nt
if nn(kk)<nt; D_out(n(kk),k) = D_in(nn(kk),k);end;
end
end
if conj ==-1;
for kk=1:nt
if nn(kk)<nt; D_out(nn(kk),k) = D_in(n(kk),k);end;
end
end
end
return
Kenampakan Semblance dari M=0
hingga M=50

Dari gambar diatas terlihat bahwa semblance dengan jumlah M=50 terlihat lebih
smooth karena jumlah perata-rataan window yang lebih banyak, perlu diingat bahwa
data (seismik) yang real mempunyai banyak noise sehingga semakin banyak M tentu
akan semakin baik dalam memperhalus kenampakan semblance dan peak atau puncak
dimana akan ditetapkan kecepatan terbaik tentu akan semakin baik (walaupun perlu
diingat bahwa semakin besar nilai M maka waktu komputasi akan semakin
meninggkat) .
Dibawah kembali akan diperlihatkan apabila penggunaan semblance untuk melakukan
analisa kecepatan, disini ditetapkan nilai M=20.
function tesVelan(trace,offset,sampling)
% fungsi untuk melihat kenampakan semblance
% untuk berbagai macam kecepatan
% created by randy on 2003

% jika fungsi ini dijalankan setelah


% fungsi sintetik.m maka tidak perlu input
if nargin==0
load trace
offset=0:50:500;
ns=1000;
time_maks=2;
sampling=time_maks./ns;
end
% menentukan nilai M
M=20;
% looping dimana variabel kecepatan dirubah
% untuk setiap perulangan
for i=1:10
% menentukan velocity moveout
% mulai dari 200 hingga 2000
v=200.*i;
% melakukan NMO
[trace_out] = cvnmo(trace,offset,sampling,v,1);
% melakukan penggabungan data setelah moveout
if i==1
nmotrace=trace_out;
else
nmotrace=[nmotrace trace_out];
end
sem(:,i)=semblance(trace_out,M);
end
figure;
imagesc(sem)
figure;
imagesc(nmotrace);
colormap(gray);

Dari fungsi diatas dapat terlihat bahwa pada nilai kecepatan v=1000 lah terlihat nilai
koherensi yang paling tinggi sehingga dapat disimpulkan bahwa dengan kecepatan
inilah moveout dari data dapat menjadi event yang paling flat atau datar.

Nilai V=1000

Analisa Kecepatan dari data


sintetik
Nilai V=1000

Kenampakan trace setelah


dimoveout
References

Scales, Jhon.A., Theory of Seismic Imaging, Samizdat Press, Center for Wave
Phenomena and New England Research, 1997.

Barryman, Jmes .G ., Non Linear Inversion and Tomography: 1. Borehole Seismic


Tomography, A Lecture Notes, Massachusetts Institute of Technology, 1991.

Borchers, Brian ., Geophysical Inverse Theory Notes ., http://www.nmt.edu , 2000.

Vermeer, Gijs .J.O ., Seismic Wavefield Sampling, A Wavenumber Approach to


Acquisition Fundamentals ., Society of Exploration Geophysicists ., Shell Research
B.V ., 1990.

Yilmaz , Ozdogan ., Seismic Data Processing ., Society of Exploration Geophysicists


., 1987.

Li, Qiang ., High Resolution Hyperbolic Radon Transform Multiple Removal .,


Master of Science Thesis ., Dept. of Physics .,University of Alberta ., Canada ., 2001.

Sacchi, Mauricio .D ., Statistical and Transform Methods for Seismic Signal


Processing ., Dept. of Physics .,University of Alberta ., Canada., 1999.

Sacchi, Mauricio .D ., Refraction and Reflection Seismology ; Lecture Notes .,


Dept. of Physics .,University of Alberta ., Canada ., 1999.

Edwards Jr, C.H ., Penney , David.E ., Elementary Linear Algebra ., Prentice-Hall


,Inc ., New Jersey ., 1988.

Boas, Mary.L., Mathematical Methods in The Physical Sciences, Jhon Wiley &
Sons, Inc, Canada, 1983.

Englander, Irv., The Architecture of Computer Hardware and Systems Software,


Jhon Wiley & Sons, Inc, USA, 2003

Curran, James.M, Little Endian Vs. Big Endian, http://noveltheory.com , 1997.

SEG Technical Standards Committee ., SEG Y Rev 1 Data Exchange Format .,


Society of Exploration Geophysicists ., 2002.

Stockwell Jr, Jhon .W, Cohen, Jack K, The New SU User Manual, Society of
Exploration Geophysicists, The Gas Research Institute and Center for Wave
Phenomena, 1998

Welsh, Matt, Slackware Linux, Installation & Getting Started, Walnut Creek
CDROM Books, 1994.

Zhang, Tony, Teach Yourself C in 24 Hours, Sams Publishing, Indianapolis,


Indiana, 1997
Husain, Kamran, Teach Yourself Unix Shell Programming in 14 days, Sams
Publishing, Indianapolis, Indiana, 1994.

You might also like