Ensiklopedi Seismik Online

Advanced SU: Post-Stack Depth Migration

2012-01-13T14:40:00.000-08:00

by Befriko Murdianto

Seismic Unix menyediakan beberapa fungsi untuk melakukan Post-Stack Depth Migration (PoSDM), antara lain sumigsplit yang menggunakan metode split-step Fourier migration (Stoffa et al, 1990) dan sumigpspi yang menggunakan metode phase-shift plus interpolation – PSPI (Gazdag and Sguazzero, 1984).

Untuk mendemonstrasikan PosDM pada SU, data sintetik SEG/EAGE salt model akan digunakan, data tersebut yang bisa diunduh di sini. Data ini berupa penampang stack dalam two-way-time yang sebelum migrasi. Velocity model dalam depth untuk data sintetik ini bisa diunduh di sini. Baik data sintetik maupun velocity model disimpan dalam format C-style binary format, yaitu SU format tanpa trace headers. Karena tidak ada informasi trace headers, maka berikut ini diberikan informasi yang diperlukan untuk menampilkan data tersebut:

Grid spacing velocity model arah lateral (dx): 40 ft
Grid spacing velocity model arah vertical (dz): 40 ft
Sample interval: 8 msec
Jumlah time sample data sintetik: 626
Jumlah depth sample velocity model: 300
Jumlah CDP data sintetik dan velocity model: 1290

Untuk menampilkan velocity model dari bentuk salt model ini bisa kita gunakan perintah ximage.

ximage <velocities.le.bin n1=300 d1=40 d2=40 label1="(ft)"="(ft)" title="SEG/EAGE Salt Model" legend=1 units="ft/sec" cmap=hsv6 &

Tampilan dari perintah di atas adalah sebagai berikut:

Agar data sintetik dapat ditampilkan dengan perintah SU, kita harus menambahkan trace header dulu supaya binary file tersebut menjadi SU format.

suaddhead ns=626 <data.le.bin | sushw key=dt a=8000 >data.le.su

Lalu dapat kita tampilkan data sintetik tersebut menggunakan perintah suximage.

suximage <data.le.su label1="Time (sec)" label2="CDP" title="Stack Section of SEG/EAGE Salt Model" legend=1 units="Amplitude" bclip=0.02 wclip=-0.02 &

Tampilannya adalah sebagai berikut:

Untuk melakukan PoSDM di SU, velocity file yang akan dipakai untuk PoSDM harus memiliki struktur vfile[iz][ix], yang berarti sumbu x merupakan “fast direction” (vertical) sedangkan sumbu z merupakan “slow direction” (lateral). Hal ini terbalik dengan velocity model yang kita miliki sekarang, yaitu sumbu z merupakan “fast direction” dan sumbu x merupakan “slow direction”. Untuk mengatasinya, kita dapat menggunakan perintah transp di SU untuk mentranspose matriks dari velocity model tersebut

transp <velocities.le.bin n1=300 >velocities.le.transp

Apabila kita tampilkan hasilnya menggunakan ximage

ximage <velocities.le.transp n1=1290 d1=40 d2=40 label1="Distance (ft)" label2="Depth (ft)" title="SEG/EAGE Salt Model (Transposed)" legend=1 units="ft/sec" cmap=hsv6 &

Tampilannya adalah seperti di bawah ini:

Sekarang kita sudah siap untuk melakukan PoSDM. Pertama-tama kita akan menggunakan metode split-step Fourier migration.

sumigsplit <data.le.su vfile=velocities.le.transp nz=300 dz=40 dx=40 >migsplit.su

Tampilkan hasilnya.

suximage <migsplit.su label1="Depth (ft)" label2="CDP" title="Split-step Fourier" legend=1 units="Amplitude" bclip=0.03 wclip=-0.03 &

Selanjutnya kita coba PoSDM menggunakan metode PSPI.

sumigpspi <data.le.su vfile=velocities.le.transp nz=300 dz=40 dx=40 >migpspi.su

Tamplikan hasilnya.

suximage <migpspi.su label1="Depth (ft)" label2="CDP" title="PSPI" legend=1 units="Amplitude" bclip=0.03 wclip=-0.03 &

Terlihat bahwa PSPI memberikan hasil yang sedikit lebih baik untuk subsalt imaging, meskipun steep dip imagingnya tidak sebaik split-step Fourier.

Berikut ini adalah perbandingannya dalam animated gif.

Referensi:
Stoffa, P. L., Fokkema, J. T., Freire, R. M. and Kessinger, W. P., 1990, Split-step Fourier migration, Geophysics, 55, 410-421.
Gazdag, J. and Sguazzero, P., 1984, Migration of seismic data by phase-shift plus interpolation, Geophysics, 49, 124-131.

Advanced SU: FX Decon

2012-01-10T03:45:00.000-08:00

Walaupun secara teknis FX Decon tidak memberikan hasil sebaik Curvelet, tetapi pada kondisi tertentu FX Decon dapat dimanfaatkan untuk meredam noise yang bersifat random (random noise), linear noise bahkan sisa groundroll dan multiple yang berfrekuensi tinggi.

Secara teknis, masing-masing trace dalam domain T-X (waktu-offset) ditransform kedalam frekuensi domain, untuk masing-masing frekuensi, Winner Filter digunakan untuk memprediksi sampel berikutnya, lalu pada akhir proses data ditransformasi balik ke dalam domain T-X.

Pada script berikut saya menerapkan sufxdecon (FX Decon Seismic Unix) untuk sebuah CMP dari data marine setelah koreksi NMO.

Parameter penting yang berpengaruh terhadap hasil FX Decon adalah fmin, fmax, ntrf, ntrw dan taper. Anda silakan melakukan test dengan nilai parameter-parameter yang disesuaikan dengan kondisi data yang anda miliki dengan syarat jangan sampai menghilangkan reflektor dan difraksi.

#!/bin/sh

sunmo < cdp800.su tnmo=0,1,2,4 vnmo=1500,1600,2500,3000 | suwind tmin=3 tmax=4.5 > junk1.su
sunmo < cdp800.su tnmo=0,1,2,4 vnmo=1500,1600,2500,3000 | sufxdecon fmin=3 fmax=60 ntrf=77 ntrw=78 taper=0.01 | suwind tmin=3 tmax=4.5 > junk2.su

suximage < junk1.su perc=95 title='before' &
suximage < junk2.su perc=95 title='after' &
suop2 junk1.su junk2.su op=diff | suximage perc=95 title='diff' &

Sebelum FX Decon:

Setelah FX Decon:

Noise yang dihilangkan:

Perbandingannya:

Advanced SU: Median Filter

2012-01-09T06:14:00.000-08:00

Seismic Unix menyediakan fungsi sumedian yang berguna untuk mengeliminasi event seismik dengan moveout (baca: kemiringan) tertentu.

Fungsi ini sangat berguna baik dalam pengolahan data seismic refleksi maupun Vertical Seismic Profiling (VSP).

Untuk mendemonstrasikan kegunaan sumedian pada data VSP silakan download data VSP di sini.

Link di atas direlease oleh Pemerintah Negara Bagian Victoria, dimana data tersebut hanya boleh digunakan untuk keperluan Studi dan Penelitian.

Selain data VSP terdapat juga link yang berisi data seismik (link1-link2), well-log, dll.

Untuk data VSP, saya akan menggunakan data VSI_007_A_gac_wavefield_z.sgy

Convert sgy ke su:

segyread tape=VSI_007_A_gac_wavefield_z.sgy verbose=1 endian=0 > VSI_007_A_gac_wavefield_z.su

Gunakan shell berikut untuk melakukan windowing, menampilkan dan median filter.

#!/bin/sh
suwind < VSI_007_A_gac_wavefield_z.su key=tracl | suwind key=tracl min=20 max=120 tmin=0.1 tmax=0.6 | suximage perc=90 d1=0.002 title='ORIGINAL'&
suwind < VSI_007_A_gac_wavefield_z.su key=tracl | sumedian tshift=1.25,0.6 xshift=1.11,149.436 median=1 | sumedian tshift=0.396025,0.972893 xshift=3.0093, 47.5488 median=1 | suwind key=tracl min=20 max=120 tmin=0.1 tmax=0.6 | suximage perc=90 d1=0.002 title='DOWNGOING WAVES'&

suwind < VSI_007_A_gac_wavefield_z.su key=tracl | sumedian tshift=0.11,0.76 xshift=1.65,215.53 median=1 | sumedian tshift=0.396025,0.972893 xshift=3.0093, 47.5488 median=1 | suwind key=tracl min=20 max=120 tmin=0.1 tmax=0.6 | suximage perc=90 d1=0.002 title='UPGOING WAVES'&

Parameter penting sumedian adalah tshift (nilai sumbu Y) dan xshift (nilai sumbu X). tshift-xshift merupakan pasangan beberapa point untuk event-event dengan kemiringan tertentu. tshift dapat berupa waktu, depth, dll. dan xshift dapat berupa offset, tracl, dll.

Di bawah ini adalah rekaman VSP yang masih memiliki upgoing waves dan downgoing waves.

Downgoing waves:

Upgoing waves:
Perbandingannya:

Well-Log Data

2012-01-06T02:19:00.000-08:00

Pembahasan mengenai beberapa konsep well logging telah dijelaskan pada artikel Gamma Ray, Resistivity, Neutron Porosity dan Density Logging.

Untuk keperluan studi dan penelitian, well log data dapat anda peroleh di sini, yang dipublikasikan oleh Department of the Interior U.S. Geological Survey untuk 14 Wildcat Wells di Alaska.

Ke-14 data well tersebut memiliki data log dengan format LAS yang anda dapat buka dengan text editor biasa seperti notepad, nedit atau gedit.

Berikut adalah contoh well-log untuk Tulageak-1 (klik untuk memperbesar):

Data well-log dengan format LAS dapat ditampilkan dengan software gratis dengan platform Windows seperti LASReader atau dengan software gratis lainnya seperti octave atau gnuplot.

Berikut adalah contoh cara menampilkan dan memanipulasi Tulageak-1 well log data dengan octave atau matlab.

Dengan Text editor bukalah TL1.LAS sehingga anda memperoleh penampilan seperti di bawah ini:

Baris 1 s/d 35 merupakan LAS header yang berisikan informasi mengenai data ini i.e. kolom pertama merupakan measured depth dengan satuan ft, kolom kedua merupakan data Spontaneous Potential (MV), dst.

Perhatikan pada tabel tersebut terdapat angka -999.00000, angka ini merupakan null value yang artinya tidak ada informasi pengukuran.

Pada terminal Linux, gunakan perintah sed untuk mengganti nilai -999.0000 dengan NaN (Not a Numeric).

sed -e 's/-999.00000/NaN/g' TL1.LAS > TL1_edit.LAS

Lalu gunakan kode berikut untuk membaca dan menampilkan. Download function hdrload.m terlebih dahulu.

clear; clc
[h, d] = hdrload('TL1_edit.LAS');
%plot gammaray
subplot(1,4,1);
plot(d(:,3),d(:,1),'r');
set (gca (), 'ydir', 'reverse','XGrid','on','YGrid','on');
axis([min(d(:,3)) max(d(:,3)) min(d(:,1)) max(d(:,1))]);
xlabel('GR(API)');
ylabel('Measured Depth(ft)');

%plot LLD
subplot(1,4,2);
semilogx(d(:,6),d(:,1),'m');
set (gca (),'ydir', 'reverse','XGrid','on','YGrid','on');
axis([min(d(:,6)) max(d(:,6)) min(d(:,1)) max(d(:,1))]);
xlabel('LLD(ohmM)');

%plot density
subplot(1,4,3);
plot(d(:,8),d(:,1),'b');
set (gca (), 'ydir', 'reverse','XGrid','on','YGrid','on');
axis([min(d(:,8)) max(d(:,8)) min(d(:,1)) max(d(:,1))]);
xlabel('Density(g/cc)');

%plot sonic
subplot(1,4,4);
plot(d(:,10),d(:,1),'g');
set (gca (), 'ydir', 'reverse','XGrid','on','YGrid','on');
axis([min(d(:,10)) max(d(:,10)) min(d(:,1)) max(d(:,1))]);
xlabel('Sonic(us/ft)');

Lalu kita akan melakukan smoothing atau despike dengan metoda Whittaker Method, dimana kodenya bisa didownload di sini.

Perhatikan manipulasi NaN values sebelum diterapkan whittf.

figure;
sonic=d(:,10);
depth=d(:,1);
loc=find(isnan(sonic) == 0);
sonic=sonic(loc);
depth=depth(loc);
plot(sonic,depth,'g',whittf(sonic,10^3),depth,'r',whittf(sonic,10^9),depth,'b');
legend('original','smooth-10^3','smooth-10^9')
set (gca (), 'ydir', 'reverse', 'xdir', 'reverse','XGrid','on','YGrid','on');
axis([min(sonic) max(sonic) min(depth) max(depth)]);
xlabel('Sonic(us/ft)');
ylabel('Measured Depth(ft)');

Di bawah ini saya melakukan despike untuk data sonic dan density, konversi data sonic menjadi kecepatan serta kedalaman dari feet ke meter dan lalu menyimpannya dalam format ascii.

%sonic
sonic=d(:,10);
depth=d(:,1);
loc_sonic=find(isnan(sonic) == 0);
sonic_edited=whittf((1000000./(3.28084.*sonic(loc_sonic))),10^4);
depth_sonic=depth(loc_sonic)./3.28084;
subplot(1,2,1)
plot(sonic_edited,depth_sonic);set (gca (), 'ydir', 'reverse','XGrid','on','YGrid','on');
xlabel('Velocity(m/s)');
ylabel('Measured Depth(m)');

%%%density
density=d(:,8);
depth=d(:,1);
loc_density=find(isnan(density) == 0);
density_edited=whittf(density(loc_density),10^4);
depth_density=depth(loc_density)./3.28084;
subplot(1,2,2)
plot(density_edited,depth_density);set (gca (), 'ydir', 'reverse','XGrid','on','YGrid','on');
xlabel('density(g/cc)');
ylabel('Measured Depth(m)');

dvfile=[depth_sonic,sonic_edited];
drfile=[depth_density,density_edited];

save dvfile.asc dvfile -ascii
save drfile.asc drfile -ascii

Konversi ascii ke binary, dan membuat seismogram sintetik dengan Seismic Unix.

a2b n1=2 < dvfile.asc > dvfile.bin
a2b n1=2 < drfile.asc > drfile.bin

suwellrf dvfile=dvfile.bin drfile=drfile.bin ntr=100 nval=7591 | suxwigb &

Check-Shot Survey

2012-01-03T06:33:00.000-08:00

Check-Shot survey adalah survey pengukuran waktu tempuh gelombang seismik, dimana posisi sumber gelombang diletakkan di permukaan-dekat lubang bor sementara perekam berada di dalam lubang bor.

Gambar di bawah ini adalah ilustrasi geometri Check-Shot Survey.

Courtesy Halliburton

Sumber gelombang biasanya berupa air-gun (500 cu in). Jika survey-nya di darat, air gun disimpan di dalam sebuah bak air atau kolam. Air-gun digunakan karena lebih repeatable dibandingkan dinamit. Sedangkan alat perekam sendiri biasanya berupa geophone accelerator 3 komponen (X,Y,Z) dan berada pada 4 lokasi (4 shuttles) untuk satu rangkaian alat.

Perekaman dilakukan pada beberapa titik kedalaman lubang bor baik sebelum ataupun setelah dipasang casing. Interval kedalaman biasanya sekitar 100m atau sesuai dengan kebutuhan.

Gambar di bawah ini adalah ilustrasi gelombang langsung (direct wave) yang berupa minimum phase dengan first break (garis hitam miring) berupa waktu tempuh gelombang seismik .

Courtesy Halliburton

QC utama yang harus dilakukan pada rekaman Check-Shot adalah kejelasan first break yang bersih dari noise baik untuk komponen X,Y ataupun Z dan peningkatan waktu tempuh sejalan dengan penambahan kedalaman (jika posisi lubang bor vertikal atau miring). Hal ini akan berbeda jika pada horizontal well. Kualitas data dipengaruhi oleh kesehatan alat, coupling antara alat dan lubang bor, kehadiran gerowong yang berada dibalik casing, dll.

Produk utama dari Check-Shot adalah kurva hubungan waktu tempuh dengan kedalaman (TD Curve) yang sangat berguna untuk konversi waktu ke kedalaman, mengkoreksi sonic-sonic corrected check-shot *) untuk keperluan pembuatan seismogram sintetik, memperbaiki kecepatan seismic (velocity scaling), estimasi Q, dll.

Gambar di bawah ini menunjukkan kurva waktu tempuh dan kedalaman yang di-overlay dengan first break Check-Shot (kiri), sonic corrected checkshot (tengah), sintetik seismogram dan data seismik di sekitar lubang bor.

Courtesy Halliburton

*) Sonic corrected checkshot (data sonic dikoreksi oleh check shot): hal ini dilakukan karena diantaranya terjadi perbedaan orientasi ray path (jejak sinar). Jejak sinar seismik pada survey checkshot relatif tegak lurus terhadap perlapisan batuan sedangkan sonic cenderung sejajar. Perbedaan ini menghasilkan perbedaan waktu tempuh yang kita kenal dengan drift curve. Kalkulasi drift curve sudah saya jelaskan di sini.

Seismic Refraction Ray Tracing

2011-12-31T19:27:00.000-08:00

Metoda Seismik Refraksi merupakan metoda yang cukup populer digunakan baik untuk keperluan eksplorasi yang bersifat bisnis ataupun untuk keperluan riset.

Seperti yang saya jelaskan di sini, seismik refraksi digunakan untuk pendeteksian zona dekat permukaan (near surface zone) ataupun Refraction Tomo yang digunakan untuk koreksi statik yang mendukung keberhasilan pengolahan data seismik refleksi.

Seismik Refraksi kerap digunakan didalam global seismology untuk pendeteksian lapisan bumi dalam seperti zona subduksi, moho, mantel, dll.

Teknik dan kualitas Ray Tracing merupakan merupakan kunci keberhasilan pencitraan bawah permukaan. Terdapat berbagai macam metoda ray tracing, diantaranya Fast Marching Method, Finite Difference, Shooting Method, Bending Method, dll.

Untuk mempelajari Ray Tracing silakan anda download kode refract_raytracer.m serta data pendukungnya ak135.txt dan attenuationdata.txt.

Kode tersebut akan membaca lokasi sumber dan penerima serta kedalaman sumber dalam latitude dan longitude dengan model kecepatan ak135. Pada kode tersebut saya menggunakan Shooting Method yang melakukan iterasi beberapa kali sampai selisih error dibawah threshold tertentu.

Gambar di bawah ini merupakan penampang Jejak Sinar dalam 2D yang di-overlay dengan model kecepatan ak135.

Penampang Jejak Sinar dalam 3D (click untuk memperbesar)

Untuk mempelajari lebih jauh mengenai penerapan metoda ini, silakan download paper saya di sini.

Advanced SU: Seismic Attributes

2011-12-27T07:15:00.000-08:00

Seismic Unix menyediakan fungsi suattributes yang berguna untuk menghitung Instantaneous Attributes seperti Amplitude Envelope, Instantaneous Phase, Instantaneous Frequency, Instantaneous Q Factor, dll.

Untuk mendemonstrasikan penggunaaan suattributes saya akan menggunakan porsi data marine yang telah saya bahas di sini. Proses yang telah saya terapkan terhadap data marine tersebut adalah radon demultiple, curvelet transform , PreStack Time Migration dan Outer Mute.

Berikut adalah contoh-contoh penerapan suattributes.

Data original
suximage < stackdata.su perc=95 wbox=1000 hbox=600 title='Original' cmap=rgb1 legend=1 &

Instantaneous Phase
suattributes < stackdata.su mode=phase | suximage perc=95 wbox=1000 hbox=600 title='Instantaneous Phase' cmap=hsv2 legend=1&

Instantaneous Frequency
suattributes < stackdata.su mode=freq | suximage perc=95 wbox=1000 hbox=600 title='Instantaneous Frequency' cmap=hsv1 legend=1&

Amplitude Envelope
suattributes < stackdata.su mode=amp | suximage perc=95 wbox=1000 hbox=600 title='Amplitude Envelope' cmap=rgb1 legend=1&

The 2004 BP Velocity-Analysis Benchmark

2011-12-21T06:02:00.001-08:00

Model ini merupakan courtesy British Petroleum (BP) yang pada awalnya dibuat oleh Frederic Billette dan Sverre Brandsberg-Dahl pada tahun 2004.

Bagian kiri dari model 2D ini berdasarkan model geologi bagian barat Teluk Meksiko. Bagian tengahnya merupakan penyederhanaan dari model bagian Timur atau Tengah Teluk Meksiko dan lepas pantai Angola. Sedangkan bagian kanannya merepresentasikan model di Laut Caspia, North Sea atau Trinidad.

Model kecepatan dengan format SEGY dapat didownload di sini. Selain model kecepatan terdapat juga data seismik shot gathers serta data density.

Model kompleks ini sangat berguna bagi anda para mahasiswa Tugas Akhir atau Seismic Geek untuk bereksperimentasi Teknologi Seismik Lanjut.

Gambar di bawah ini merupakan hasil Reverse Time Migration, courtesy GX Technology.

Advanced SU: Mute Design

2011-12-20T04:14:00.000-08:00

Pada bagian ini saya akan mendemonstrasikan bagaimana caranya melakukan desain dan mengimplementasikan mute dengan Seismic Unix.
Dalam contoh ini, mute desain berubah dari satu cdp ke cdp yang lainnya.

Pastikan anda memiliki data dalam domain cdp-offset yang menghindari stretch mute pada saat nmo yang mana dapat dilakukan dengan meningkatkan nilai smute pada parameter sunmo.

1. Lakukan pemecahan data su dengan perintah suplit dengan kunci cdp.

susplit < data_cdp_nmo.su > cdp_split key=cdp

2. Pindahkan data yang telah dipecah ke sebuah folder terpisah, lalu lakukan rename sehingga sequential:

#!/bin/sh
a=1
for i in *.su;
do
new=$(printf "%0d.su" ${a})
mv ${i} ${new}
let a=a+1
done

3. Ekstrak informasi header dari masing-masing cdp lalu meng-outputkannya ke sebuah file.

#! /bin/sh
datamin=300
datamax=500
k=$datamin
while [ $k -le $datamax ]
do
surange < $k.su > surange_$k
k=$((k+1))
done

4. Lakukan display mute dan pick

#!/bin/sh
datamin=300
datamax=350
interval=10
k=$datamin
while [ $k -le $datamax ]
do
offset=`sed -n "11p" surange_$k |awk '{print $2} '`
suximage < $k.su d1=0.002 f1=0 d2=100 f2=$offset perc=95 title=$k mpicks=mute$k.txt wbox=500 hbox=800 &
read response
k=$((k+$interval))
done

Pada kode di atas data minimum adalah 300 dan maksimum 350 dengan interval pick 10.

Pastikan display data seismik dalam sumbu y adalah TWT dalam (second) dan sumbu x dalam offset (meter)-lihat gambar di bawah Hal ini bisa dikontrol dengan d1=0.002 (sampling interval 2ms), f1=0 waktu pertama, d2=100 penambahan offset (m) dan f2 adalah offset trace pertama dari cdp yang bersangkutan yang diekstrak oleh perintah sed dan awk dengan input surange_$k. Proses mute biasanya dilakukan untuk keperluan stack, sehingga eliminasi noise seperti groundroll harus dilakukan terlebih dahulu. Pada gambar tersebut saya belum menerapkan groundroll removal.

Saya menyarankan untuk melakukan 3 buah picks pada masing-masing display (lihat contoh di atas) dengan mengarahkan kursor pada pasangan offset-waktu tertentu dengan menekan huruf 's' jika selesai tekan huruf 'q' lalu enter, maka display cdp berikutnya akan muncul. Lakukan hal serupa!

Jika proses tersebut selesai, maka data pick akan disimpan pada file mute$k.txt.

5. Untuk memperoleh desain mute pada seluruh CDP anda bisa melakukan interpolasi dari data yang anda pick atau anda bisa menggunakan nilai pick yang sama untuk cdp yang berdekatan. Misal pick pada cdp300 akan digunakan oleh data 300 s/d 305. Sedangkan pick 310 akan digunakan pada cdp 306 s/d 315, dan seterusnya. Hal ini bisa dilakukan dengan perintah berikut:

Buat list perintah copy:

#!/bin/sh
min=300
max=350
n=$min
while [ $n -le $max ]
do
echo cp mutenew$n.txt >>junk1
n=$((n+1))
done

Kumpulkan mute$k.txt menjadi sebuah file, gandakan 10 kali dan membuang 4 data pertama dan 5 data terakhir.
ls *.txt* | awk '{ for (i = 1; i <=10 ; i++) print }' | awk '(NR>4 && NR<56) {print $0} '> junk2

Lalu gabungkan:

paste junk1 junk2 | awk '{print $1,$3,$2}' > junk3

Copy-paste text yang terbentuk pada junk3 pada terminal.
more junk3
cp mute300.txt mutenew300.txt
cp mute300.txt mutenew301.txt
cp mute300.txt mutenew302.txt
cp mute300.txt mutenew303.txt
cp mute300.txt mutenew304.txt
cp mute300.txt mutenew305.txt
cp mute310.txt mutenew306.txt
cp mute310.txt mutenew307.txt
cp mute310.txt mutenew308.txt
cp mute310.txt mutenew309.txt
cp mute310.txt mutenew310.txt
cp mute310.txt mutenew311.txt
cp mute310.txt mutenew312.txt
cp mute310.txt mutenew313.txt
.
.
.

6. Untuk menghindari edge effect pada desain mute, gantilah nilai offset pertama dan terakhir dengan nilai offset min-max yang terbaca oleh surange:

#!/bin/sh
datamin=300
datamax=350
interval=1
k=$datamin
while [ $k -le $datamax ]
do
negoffset=`sed -n "11p" surange_$k | awk '{print $2} '`
posoffset=`sed -n "11p" surange_$k | awk '{print $3} '`
awk -v line=1 -v col=2 -v var=$negoffset '{ if(NR==line) { $col=var}; print $0}' mutenew$k.txt | awk -v line=3 -v col=2 -v var=$posoffset '{ if(NR==line) { $col=var}; print $0}' > muteok$k.txt
k=$((k+$interval))
done

7. Menerapkan mute untuk seluruh cdp dengan parameter mute masing-masing.

#!/bin/sh
datamin=300
datamax=350
interval=1
k=$datamin
while [ $k -le $datamax ]
do
mkparfile < muteok$k.txt > muteok$k.par string1=tmute string2=xmute
sumute mode=0 < $k.su par=muteok$k.par > cdp_after_mute$k.su
k=$((k+$interval))
done

cdp yang telah dimute

8. Untuk menggabungkan hasilnya lakukan perintah cat seperti yang telah dijelaskan disini.

e-book ensiklopediseismik

2011-12-14T03:53:00.000-08:00

e-book ensiklopediseismik adalah kompilasi artikel yang telah saya publish di blog ini.

Format e-book tersebut adalah chm yang mana anda bisa membaca artikel-artikel secara offline dari sebuah file tunggal.

chm dapat dibuka secara otomatis pada platform windows, untuk ubuntu anda bisa menginstall xchm terlebih dahulu dengan perintah sudo apt-get install xchm

Untuk display yang lebih baik, saya sarankan untuk membukanya dengan platform Windows.

e-book tersebut (37.2MB) dapat di download di sini.

Berikut adalah screenshot e-book tersebut.

Animated gif

2011-12-08T08:44:00.000-08:00

Pada bagian ini saya akan membahas bagaimana caranya menghasilkan animated gif seperti contoh di bawah ini pada Ubuntu-Linux (click untuk memperbesar).

Waveform Modeling untuk Model Marmousi

Capture suximage dari Model Marmousi dengan screen-capture tool favorit anda seperti Take Screenshot (save marmousi.png), lalu lakukan waveform modeling seperti yang telah saya bahas di sini.

Jika suxmovie telah aktif, maka lakukan langkah-langkah sbb:

1. sudo apt-get install recordmydesktop gtk-recordmydesktop
2. Applications > Sound & Video > Desktop recorder (sebuah file ogv akan terbentuk setelah anda menekan stop-kotak merah di atas monitor)

3. Convert ogv ke avi
mencoder -idx out.ogv -ovc lavc -oac mp3lame -o out.avi

4. Lakukan pemilihan rentang waktu yang akan dipilih (contoh mulai 0:02 detik sampai 8 detik berikutnya)
mplayer -ao null -loop 0 -ss 0:02 -endpos +8 out.avi

5. Jika anda setuju dengan tahap 4 lakukan destilasi file jpg
mplayer -ao null -ss 0:02 -endpos +8 out.avi -vo jpeg:outdir=moviepics

6. Ubah format jpg ke gif
gedit junk1.sh
#!/bin/sh
echo "#!/bin/sh" >> junk2.sh
for file in *.jpg
do echo "convert $file $(basename $file .jpg).gif" >> junk2.sh
done

Jalankan dengan sh junk2.sh

7. rename data sehingga filename menjadi sequential

a=1
for i in *.gif; do
new=$(printf "%0d.gif" ${a})
mv ${i} ${new}
let a=a+1
done

8. Lakukan merging antara gambar yang telah didestilasi dengan model marmousi. (args=50 adalah derajat transparansi).

gedit junk3.sh
#! /bin/sh
datamin=1
datamax=61
n=$datamin
while [ $n -le $datamax ]
do
convert marmousi.png $n.gif -alpha on -compose blend -define compose:args=50 -gravity South -composite merge$n.gif
n=$((n+1))
done

Jalankan dengan sh junk3.sh

9. Buat perintah animasi

gedit junk4.sh
#!/bin/sh
echo "#!/bin/sh" >> junk5.sh
echo "gifsicle --delay=10 --loop=1000 --colors 256 \\" >> junk5.sh
min=1
max=61
n=$min
while [ $n -le $max ]
do
echo "merge$n.gif \\" >>junk5.sh
n=$((n+1))
done
echo "> anim.gif" >> junk5.sh

10. sudo apt-get install gifsicle

11. Jalankan sh junk5.sh

Curvelet Transform

2011-12-05T02:43:00.000-08:00

Reduksi noise (denoising) merupakan salah satu pekerjaan penting didalam pengolahan data seismik. Dengan mereduksi noise diharapkan kita memperoleh hasil dengan S/N (Signal to Noise Ratio) yang tinggi. Pekerjaan ini cukup menantang terutama jika noise dan sinyal (dalam berbagai domain) tidak terpisahkan dengan baik.

Curvelet Transform merupakan teknologi yang diklaim mampu memisahkan coherent-random noise dari sinyal yang kita kehendaki, karena Transformasi Curvelet akan memisahkan noise dari sinyal dalam dimensi frekuensi, dip, azimuth dan lokasi. Karena kelebihannya ini Curvelet Transform dikenal dengan Transformasi multi dimensi yang tidak dimiliki oleh teknik konvensional seperti Median Filtering dan FX-Decon.

Curvelet Transform diadopsi oleh geoscientist dari Image Processing dan Scientific Computing, dimana pada kedua bidang tersebut telah lama diterapkan untuk berbagai keperluan diantaranya denoising.

Gambar di bawah ini merupakan contoh denoising image Monas dengan kadar noise tertentu. Anda dapat melakukan ekperimentasi ini dengan mendownload kode matlab atau scilab Curvelet Transform di sini dan di sini atau penerapannya pada data seismik yang bisa didownload di sini.

Aplikasi Curvelet pada Image Processing sebelum (kiri) dan setelah denoising (kanan)

Prinsip denoising dengan Curvelet Transform adalah sbb: melakukan FFT 2D pada data seismic (Stack atau CMP gather yang telah di NMO), melakukan window-ing dari hasil pertama dengan Polar Grid (gambar di bawah) dalam domain domain frekuensi dengan demikian diperoleh Curvelet Coefficient, melakukan scalling (denoising) dari koefisien curvelet tersebut, melakukan inversi untuk memperoleh kembali data seismik yang telah di-denoise.

Polar grid Curvelet dalam domain frekuensi (a) dan domain spatial (b) (Courtesy Douma and de Hoop)

Gambar di bawah ini adalah Curvelet dalam domain spatial (kolom pertama) dan amplitude spectra pada polar grid (kolom ketiga) untuk masing-masing Curvelet. Perhatikan bahwa curvelet untuk masing-masing lokasi polar grid memiliki ukuran(size), arah(azimuth), lokasi dan kemiringan tertentu.

Curvelet dan Amplitude Spectra pada Polar Grid (Courtesy Douma and de Hoop)

Gambar di bawah ini adalah aplikasi Curvelet Transform pada data seismik, cube AI yang diterapkan pada noisy data (kiri) dan AI yang diterapkan pada data setelah denoise dengan Curvelet Tranform (kanan).

Courtesy Neelami et. al.

SEGY Header Dumper Version 1.1 and 2.1

2011-11-24T04:36:00.000-08:00

Versi ini merupakan modifikasi kosmetik dari versi sebelumnya sehingga kita memperoleh tampilan seluruh informasi lokasi byte dalam satu screen.
Silakan download Version 1.1 dan Version 2.1

Berikut tampilannya:

Version 1.1

Version 2.1

SEGY Header Dumper Version 2

2011-11-22T05:34:00.000-08:00

Versi ini memiliki kapabilitas untuk melakukan scanning seluruh trace header yang dimiliki oleh data SEGY, lalu menampilkan nilai minimum dan maksimum untuk setiap lokasi byte.

Berikut Contohnya:

./headerdumper_v2 gathertest.segy
------------------------------------------------
Sample interval :6000
No of samples per trace :2002
Number of Trace :13008
---------------------------------------------
--byte#: min-max----------------------------
1- 4: 1-13008
5- 8: 0-0
9- 12: 0-0
13- 16: 0-0
17- 20: 0-0
21- 24: 9010-9820
25- 28: 1-48
29- 32: 65536-65536
33- 36: 0-0
37- 40: 1050-24179
73- 76: 1261066-1593250
77- 80: 0-0
81- 84: 1261066-1593250
85- 88: 0-0
89- 92: 131072-131072
185-188: 0-0
189-192: 1175666186-1175998370
193-196: 0-0
197-200: 1175666186-1175998370
201-204: 0-0

Silakan download kode PERL di sini, lalu ikuti petunjuk yang telah saya jelaskan di sini.

SEGY Header Dumper Version 1

2011-11-22T03:28:00.000-08:00

Pada artikel SEGY, saya telah membahas penggunaan Kode PERL untuk mengekstrak informasi yang dimiliki oleh binary header dan trace header yang pertama dari sebuah data SEGY.

Pada bagian ini anda dapat menggunakan kode PERL yang telah saya kembangkan untuk melakukan ekstrak trace header untuk setiap trace.

Berikut adalah contoh hasil dump sebuah file SEGY.

./segyheaddump1 gathertest.segy
This data is 32-bit IBM floating point
Sample interval :6000
No of samples per trace :2002
Number of Trace :13008
-----INSTRUCTION---------------
-----Press ENTER for the Next Trace--------
-----Press 9 then ENTER to EXIT------------
-----Trace Header Trace#:1 out of 13008
byte# 1- 4:1
byte# 5- 8:0
byte# 9- 12:0
byte# 13- 16:0
byte# 17- 20:0
byte# 21- 24:9010
byte# 25- 28:1
byte# 29- 32:65536
byte# 33- 36:0
byte# 37- 40:1050
byte# 41- 44:0
byte# 45- 48:0
byte# 49- 52:0
byte# 53- 56:0
byte# 57- 60:0
byte# 61- 64:0
byte# 65- 68:0
byte# 69- 72:0
byte# 73- 76:1261066
byte# 77- 80:0
byte# 81- 84:1261066
byte# 85- 88:0
byte# 89- 92:131072
byte# 93- 96:0
byte# 97-100:0
byte#101-104:0
byte#105-108:0
byte#109-112:0
byte#113-116:2002
byte#117-120:393216000
byte#181-184:0
byte#185-188:0
byte#189-192:1175666186
byte#193-196:0
byte#197-200:1175666186
byte#201-204:0

Download kode PERL di sini, lalu simpan dalam directory tempat data seismik anda berada.

Coba jalankan dengan perintah:
perl headerdumper_v1.pl file.segy

Anda dapat meng-compile kode perl dengan perintah pp. Jika anda belum memilikinya install (UBUNTU 10.10) dengan perintah:
sudo apt-get install libpar-packer-perl

Setelah itu compile dengan perintah:
pp -o segyheaderdumper headerdumper_v1.pl

Lalu run:
./segyheaderdumper file.sgy

SEGY

2011-11-19T04:16:00.001-08:00

SEGY merupakan format data seismik yang sangat populer digunakan didalam industri migas.

Sebuah file SEGY terdiri atas susunan sebagai berikut:
Blok #1: EBCDID Header dengan ukuran 3200 byte
Blok #2: Binary header dengan ukuran 400 byte
Blok #3: Trace identification header trace pertama dengan ukuran 240 byte
Blok #4: Data trace pertama (ukuran bervariasi)
Blok #3: Trace identification header trace kedua dengan ukuran 240 byte
Blok #4: Data trace kedua (ukuran bervariasi)
Dst...

Berikut adalah skema file SEGY:

Courtesy www.geo.uib.no

Blok #1 berisikan berbagai informasi seperti Nama Perusahaan, Tahun Pengambilan, Jenis Instrumen, Sampling Rate, No Trace, dll.
Blok #2 berisikan informasi seperti line id, sample int, no sample, Fold, format data, dll.
Blok #3 berisikan informasi seperti Line No, CDP, Koordinat Sumber/Penerima, dll.
Blok #4 berisikan data seismik itu sendiri.

Untuk lebih jelasnya silakan buka link berikut.

Software-software komersial menyediakan fasilitas yang dapat menampilkan semua informasi di atas.

Pada UNIX anda dapat mengekstrak informasi tersebut dengan perintah dd dan sebuah kode program.

Untuk menampilkan EBCDIC header anda dapat menggunakan perintas sbb:

dd if=data.segy conv=ascii ibs=3200 count=1

Atau

dd if=data.segy conv=ascii ibs=3200 count=1 | awk 'BEGIN{RS="C[0-9 ][0-9]"}{printf "C%2d%s\n",NR,$0}'

Berikut adalah contoh EBCDIC header dari perintah di atas (berisikan 40 baris dan 80 kolom i.e. 40*80=3200):

Untuk menampilkan Binary header dan Trace Header yang pertama, anda dapat menggunakan kode PERL di bawah ini.

Bukalah gedit lalu copy kode tersebut lalu berinama headerdumper.pl
Simpan file headerdumper.pl pada directory yang sama dengan data seismik anda, lalu pada terminal Linux ketik perintah:

perl headerdumper.pl gathertest.segy

Sehingga diperoleh hasil berikut:

#!/usr/bin/perl

die "perl headerdumper.pl file.sgy\n " if $#ARGV !=0;
$fname=$ARGV [0 ];
sysopen(IN,$fname,O_RDONLY) or die "Can't Open $fname: $! ";
sysread (IN,$binHeader,3840);
$size= -s IN;
$Samplingrate= unpack( "n ",substr($binHeader,3216,2));
$Noofsamples = unpack( "n ",substr($binHeader,3220,2));
$formatcode= unpack( "n ",substr($binHeader,3224,2));

if($formatcode==1)
{
$bytewidth=4;
$tracecount=($size-3600)/(240+$bytewidth*$Noofsamples);
} elsif($formatcode==2)
{
$bytewidth=4;
$tracecount=($size-3600)/(240+$bytewidth*$Noofsamples);
} elsif($formatcode==3)
{
$bytewidth=2;
$tracecount=($size-3600)/(240+$bytewidth*$Noofsamples);
} elsif($formatcode==4)
{
$bytewidth=4;
$tracecount=($size-3600)/(240+$bytewidth*$Noofsamples);
} elsif($formatcode==5)
{
$bytewidth=4;
$tracecount=($size-3600)/(240+$bytewidth*$Noofsamples);
} elsif($formatcode==6)
{
$bytewidth=4;
$tracecount=($size-3600)/(240+$bytewidth*$Noofsamples);
}
else {print " I dont know this type of format code\n "; }

$head_16bit_001 = unpack( "n ",substr($binHeader,3600,2));
$head_16bit_002 = unpack( "n ",substr($binHeader,3602,2));
$head_16bit_003 = unpack( "n ",substr($binHeader,3604,2));
$head_16bit_004 = unpack( "n ",substr($binHeader,3606,2));
$head_16bit_005 = unpack( "n ",substr($binHeader,3608,2));
$head_16bit_006 = unpack( "n ",substr($binHeader,3610,2));
$head_16bit_007 = unpack( "n ",substr($binHeader,3612,2));
$head_16bit_008 = unpack( "n ",substr($binHeader,3614,2));
$head_16bit_009 = unpack( "n ",substr($binHeader,3616,2));
$head_16bit_010 = unpack( "n ",substr($binHeader,3618,2));
$head_16bit_011 = unpack( "n ",substr($binHeader,3620,2));
$head_16bit_012 = unpack( "n ",substr($binHeader,3622,2));
$head_16bit_013 = unpack( "n ",substr($binHeader,3624,2));
$head_16bit_014 = unpack( "n ",substr($binHeader,3626,2));
$head_16bit_015 = unpack( "n ",substr($binHeader,3628,2));
$head_16bit_016 = unpack( "n ",substr($binHeader,3630,2));
$head_16bit_017 = unpack( "n ",substr($binHeader,3632,2));
$head_16bit_018 = unpack( "n ",substr($binHeader,3634,2));
$head_16bit_019 = unpack( "n ",substr($binHeader,3636,2));
$head_16bit_020 = unpack( "n ",substr($binHeader,3638,2));
$head_16bit_021 = unpack( "n ",substr($binHeader,3640,2));
$head_16bit_022 = unpack( "n ",substr($binHeader,3642,2));
$head_16bit_023 = unpack( "n ",substr($binHeader,3644,2));
$head_16bit_024 = unpack( "n ",substr($binHeader,3646,2));
$head_16bit_025 = unpack( "n ",substr($binHeader,3648,2));
$head_16bit_026 = unpack( "n ",substr($binHeader,3650,2));
$head_16bit_027 = unpack( "n ",substr($binHeader,3652,2));
$head_16bit_028 = unpack( "n ",substr($binHeader,3654,2));
$head_16bit_029 = unpack( "n ",substr($binHeader,3656,2));
$head_16bit_030 = unpack( "n ",substr($binHeader,3658,2));
$head_16bit_031 = unpack( "n ",substr($binHeader,3660,2));
$head_16bit_032 = unpack( "n ",substr($binHeader,3662,2));
$head_16bit_033 = unpack( "n ",substr($binHeader,3664,2));
$head_16bit_034 = unpack( "n ",substr($binHeader,3666,2));
$head_16bit_035 = unpack( "n ",substr($binHeader,3668,2));
$head_16bit_036 = unpack( "n ",substr($binHeader,3670,2));
$head_16bit_037 = unpack( "n ",substr($binHeader,3672,2));
$head_16bit_038 = unpack( "n ",substr($binHeader,3674,2));
$head_16bit_039 = unpack( "n ",substr($binHeader,3676,2));
$head_16bit_040 = unpack( "n ",substr($binHeader,3678,2));
$head_16bit_041 = unpack( "n ",substr($binHeader,3680,2));
$head_16bit_042 = unpack( "n ",substr($binHeader,3682,2));
$head_16bit_043 = unpack( "n ",substr($binHeader,3684,2));
$head_16bit_044 = unpack( "n ",substr($binHeader,3686,2));
$head_16bit_045 = unpack( "n ",substr($binHeader,3688,2));
$head_16bit_046 = unpack( "n ",substr($binHeader,3690,2));
$head_16bit_047 = unpack( "n ",substr($binHeader,3692,2));
$head_16bit_048 = unpack( "n ",substr($binHeader,3694,2));
$head_16bit_049 = unpack( "n ",substr($binHeader,3696,2));
$head_16bit_050 = unpack( "n ",substr($binHeader,3698,2));
$head_16bit_051 = unpack( "n ",substr($binHeader,3700,2));
$head_16bit_052 = unpack( "n ",substr($binHeader,3702,2));
$head_16bit_053 = unpack( "n ",substr($binHeader,3704,2));
$head_16bit_054 = unpack( "n ",substr($binHeader,3706,2));
$head_16bit_055 = unpack( "n ",substr($binHeader,3708,2));
$head_16bit_056 = unpack( "n ",substr($binHeader,3710,2));
$head_16bit_057 = unpack( "n ",substr($binHeader,3712,2));
$head_16bit_058 = unpack( "n ",substr($binHeader,3714,2));
$head_16bit_059 = unpack( "n ",substr($binHeader,3716,2));
$head_16bit_060 = unpack( "n ",substr($binHeader,3718,2));
$head_16bit_061 = unpack( "n ",substr($binHeader,3720,2));
$head_16bit_062 = unpack( "n ",substr($binHeader,3722,2));
$head_16bit_063 = unpack( "n ",substr($binHeader,3724,2));
$head_16bit_064 = unpack( "n ",substr($binHeader,3726,2));
$head_16bit_065 = unpack( "n ",substr($binHeader,3728,2));
$head_16bit_066 = unpack( "n ",substr($binHeader,3730,2));
$head_16bit_067 = unpack( "n ",substr($binHeader,3732,2));
$head_16bit_068 = unpack( "n ",substr($binHeader,3734,2));
$head_16bit_069 = unpack( "n ",substr($binHeader,3736,2));
$head_16bit_070 = unpack( "n ",substr($binHeader,3738,2));
$head_16bit_071 = unpack( "n ",substr($binHeader,3740,2));
$head_16bit_072 = unpack( "n ",substr($binHeader,3742,2));
$head_16bit_073 = unpack( "n ",substr($binHeader,3744,2));
$head_16bit_074 = unpack( "n ",substr($binHeader,3746,2));
$head_16bit_075 = unpack( "n ",substr($binHeader,3748,2));
$head_16bit_076 = unpack( "n ",substr($binHeader,3750,2));
$head_16bit_077 = unpack( "n ",substr($binHeader,3752,2));
$head_16bit_078 = unpack( "n ",substr($binHeader,3754,2));
$head_16bit_079 = unpack( "n ",substr($binHeader,3756,2));
$head_16bit_080 = unpack( "n ",substr($binHeader,3758,2));
$head_16bit_081 = unpack( "n ",substr($binHeader,3760,2));
$head_16bit_082 = unpack( "n ",substr($binHeader,3762,2));
$head_16bit_083 = unpack( "n ",substr($binHeader,3764,2));
$head_16bit_084 = unpack( "n ",substr($binHeader,3766,2));
$head_16bit_085 = unpack( "n ",substr($binHeader,3768,2));
$head_16bit_086 = unpack( "n ",substr($binHeader,3770,2));
$head_16bit_087 = unpack( "n ",substr($binHeader,3772,2));
$head_16bit_088 = unpack( "n ",substr($binHeader,3774,2));
$head_16bit_089 = unpack( "n ",substr($binHeader,3776,2));
$head_16bit_090 = unpack( "n ",substr($binHeader,3778,2));
$head_16bit_091 = unpack( "n ",substr($binHeader,3780,2));
$head_16bit_092 = unpack( "n ",substr($binHeader,3782,2));
$head_16bit_093 = unpack( "n ",substr($binHeader,3784,2));
$head_16bit_094 = unpack( "n ",substr($binHeader,3786,2));
$head_16bit_095 = unpack( "n ",substr($binHeader,3788,2));
$head_16bit_096 = unpack( "n ",substr($binHeader,3790,2));
$head_16bit_097 = unpack( "n ",substr($binHeader,3792,2));
$head_16bit_098 = unpack( "n ",substr($binHeader,3794,2));
$head_16bit_099 = unpack( "n ",substr($binHeader,3796,2));
$head_16bit_100 = unpack( "n ",substr($binHeader,3798,2));
$head_16bit_101 = unpack( "n ",substr($binHeader,3800,2));
$head_16bit_102 = unpack( "n ",substr($binHeader,3802,2));
$head_16bit_103 = unpack( "n ",substr($binHeader,3804,2));
$head_16bit_104 = unpack( "n ",substr($binHeader,3806,2));
$head_16bit_105 = unpack( "n ",substr($binHeader,3808,2));
$head_16bit_106 = unpack( "n ",substr($binHeader,3810,2));
$head_16bit_107 = unpack( "n ",substr($binHeader,3812,2));
$head_16bit_108 = unpack( "n ",substr($binHeader,3814,2));
$head_16bit_109 = unpack( "n ",substr($binHeader,3816,2));
$head_16bit_110 = unpack( "n ",substr($binHeader,3818,2));
$head_16bit_111 = unpack( "n ",substr($binHeader,3820,2));
$head_16bit_112 = unpack( "n ",substr($binHeader,3822,2));
$head_16bit_113 = unpack( "n ",substr($binHeader,3824,2));
$head_16bit_114 = unpack( "n ",substr($binHeader,3826,2));
$head_16bit_115 = unpack( "n ",substr($binHeader,3828,2));
$head_16bit_116 = unpack( "n ",substr($binHeader,3830,2));
$head_16bit_117 = unpack( "n ",substr($binHeader,3832,2));
$head_16bit_118 = unpack( "n ",substr($binHeader,3834,2));
$head_16bit_119 = unpack( "n ",substr($binHeader,3836,2));
$head_16bit_120 = unpack( "n ",substr($binHeader,3838,2));

$head_32bit_001 = unpack( "N ",substr($binHeader,3600,4));
$head_32bit_003 = unpack( "N ",substr($binHeader,3604,4));
$head_32bit_005 = unpack( "N ",substr($binHeader,3608,4));
$head_32bit_007 = unpack( "N ",substr($binHeader,3612,4));
$head_32bit_009 = unpack( "N ",substr($binHeader,3616,4));
$head_32bit_011 = unpack( "N ",substr($binHeader,3620,4));
$head_32bit_013 = unpack( "N ",substr($binHeader,3624,4));
$head_32bit_015 = unpack( "N ",substr($binHeader,3628,4));
$head_32bit_017 = unpack( "N ",substr($binHeader,3632,4));
$head_32bit_019 = unpack( "N ",substr($binHeader,3636,4));
$head_32bit_021 = unpack( "N ",substr($binHeader,3640,4));
$head_32bit_023 = unpack( "N ",substr($binHeader,3644,4));
$head_32bit_025 = unpack( "N ",substr($binHeader,3648,4));
$head_32bit_027 = unpack( "N ",substr($binHeader,3652,4));
$head_32bit_029 = unpack( "N ",substr($binHeader,3656,4));
$head_32bit_031 = unpack( "N ",substr($binHeader,3660,4));
$head_32bit_033 = unpack( "N ",substr($binHeader,3664,4));
$head_32bit_035 = unpack( "N ",substr($binHeader,3668,4));
$head_32bit_037 = unpack( "N ",substr($binHeader,3672,4));
$head_32bit_039 = unpack( "N ",substr($binHeader,3676,4));
$head_32bit_041 = unpack( "N ",substr($binHeader,3680,4));
$head_32bit_043 = unpack( "N ",substr($binHeader,3684,4));
$head_32bit_045 = unpack( "N ",substr($binHeader,3688,4));
$head_32bit_047 = unpack( "N ",substr($binHeader,3692,4));
$head_32bit_049 = unpack( "N ",substr($binHeader,3696,4));
$head_32bit_051 = unpack( "N ",substr($binHeader,3700,4));
$head_32bit_053 = unpack( "N ",substr($binHeader,3704,4));
$head_32bit_055 = unpack( "N ",substr($binHeader,3708,4));
$head_32bit_057 = unpack( "N ",substr($binHeader,3712,4));
$head_32bit_059 = unpack( "N ",substr($binHeader,3716,4));
$head_32bit_061 = unpack( "N ",substr($binHeader,3720,4));
$head_32bit_063 = unpack( "N ",substr($binHeader,3724,4));
$head_32bit_065 = unpack( "N ",substr($binHeader,3728,4));
$head_32bit_067 = unpack( "N ",substr($binHeader,3732,4));
$head_32bit_069 = unpack( "N ",substr($binHeader,3736,4));
$head_32bit_071 = unpack( "N ",substr($binHeader,3740,4));
$head_32bit_073 = unpack( "N ",substr($binHeader,3744,4));
$head_32bit_075 = unpack( "N ",substr($binHeader,3748,4));
$head_32bit_077 = unpack( "N ",substr($binHeader,3752,4));
$head_32bit_079 = unpack( "N ",substr($binHeader,3756,4));
$head_32bit_081 = unpack( "N ",substr($binHeader,3760,4));
$head_32bit_083 = unpack( "N ",substr($binHeader,3764,4));
$head_32bit_085 = unpack( "N ",substr($binHeader,3768,4));
$head_32bit_087 = unpack( "N ",substr($binHeader,3772,4));
$head_32bit_089 = unpack( "N ",substr($binHeader,3776,4));
$head_32bit_091 = unpack( "N ",substr($binHeader,3780,4));
$head_32bit_093 = unpack( "N ",substr($binHeader,3784,4));
$head_32bit_095 = unpack( "N ",substr($binHeader,3788,4));
$head_32bit_097 = unpack( "N ",substr($binHeader,3792,4));
$head_32bit_099 = unpack( "N ",substr($binHeader,3796,4));
$head_32bit_101 = unpack( "N ",substr($binHeader,3800,4));
$head_32bit_103 = unpack( "N ",substr($binHeader,3804,4));
$head_32bit_105 = unpack( "N ",substr($binHeader,3808,4));
$head_32bit_107 = unpack( "N ",substr($binHeader,3812,4));
$head_32bit_109 = unpack( "N ",substr($binHeader,3816,4));
$head_32bit_111 = unpack( "N ",substr($binHeader,3820,4));
$head_32bit_113 = unpack( "N ",substr($binHeader,3824,4));
$head_32bit_115 = unpack( "N ",substr($binHeader,3828,4));
$head_32bit_117 = unpack( "N ",substr($binHeader,3832,4));
$head_32bit_119 = unpack( "N ",substr($binHeader,3836,4));

print "--Segy Header Dumper by Agus Abdullah, PhD----\n ";
print "-----modified from Peter Zahuczki-------------\n ";
print "--------------reel header---------------------\n ";
print "Sampling Rate : ",$Samplingrate, "\n ";
print "No of Samples per Trace : ",$Noofsamples, "\n ";
print "No of Trace : ",$tracecount, "\n ";
print "Format Code : ",$formatcode, "\n ";
print "Size : ",$size, "\n ";
print "--------------First Trace Header 16 BIT-------\n ";
print " byte#1-2: ",$head_16bit_001, "\n ";
print " byte#3-4: ",$head_16bit_002, "\n ";
print " byte#5-6: ",$head_16bit_003, "\n ";
print " byte#7-8: ",$head_16bit_004, "\n ";
print " byte#9-10: ",$head_16bit_005, "\n ";
print " byte#11-12: ",$head_16bit_006, "\n ";
print " byte#13-14: ",$head_16bit_007, "\n ";
print " byte#15-16: ",$head_16bit_008, "\n ";
print " byte#17-18: ",$head_16bit_009, "\n ";
print " byte#19-20: ",$head_16bit_010, "\n ";
print " byte#21-22: ",$head_16bit_011, "\n ";
print " byte#23-24: ",$head_16bit_012, "\n ";
print " byte#25-26: ",$head_16bit_013, "\n ";
print " byte#27-28: ",$head_16bit_014, "\n ";
print " byte#29-30: ",$head_16bit_015, "\n ";
print " byte#31-32: ",$head_16bit_016, "\n ";
print " byte#33-34: ",$head_16bit_017, "\n ";
print " byte#35-36: ",$head_16bit_018, "\n ";
print " byte#37-38: ",$head_16bit_019, "\n ";
print " byte#39-40: ",$head_16bit_020, "\n ";
print " byte#41-42: ",$head_16bit_021, "\n ";
print " byte#43-44: ",$head_16bit_022, "\n ";
print " byte#45-46: ",$head_16bit_023, "\n ";
print " byte#47-48: ",$head_16bit_024, "\n ";
print " byte#49-50: ",$head_16bit_025, "\n ";
print " byte#51-52: ",$head_16bit_026, "\n ";
print " byte#53-54: ",$head_16bit_027, "\n ";
print " byte#55-56: ",$head_16bit_028, "\n ";
print " byte#57-58: ",$head_16bit_029, "\n ";
print " byte#59-60: ",$head_16bit_030, "\n ";
print " byte#61-62: ",$head_16bit_031, "\n ";
print " byte#63-64: ",$head_16bit_032, "\n ";
print " byte#65-66: ",$head_16bit_033, "\n ";
print " byte#67-68: ",$head_16bit_034, "\n ";
print " byte#69-70: ",$head_16bit_035, "\n ";
print " byte#71-72: ",$head_16bit_036, "\n ";
print " byte#73-74: ",$head_16bit_037, "\n ";
print " byte#75-76: ",$head_16bit_038, "\n ";
print " byte#77-78: ",$head_16bit_039, "\n ";
print " byte#79-80: ",$head_16bit_040, "\n ";
print " byte#81-82: ",$head_16bit_041, "\n ";
print " byte#83-84: ",$head_16bit_042, "\n ";
print " byte#85-86: ",$head_16bit_043, "\n ";
print " byte#87-88: ",$head_16bit_044, "\n ";
print " byte#89-90: ",$head_16bit_045, "\n ";
print " byte#91-92: ",$head_16bit_046, "\n ";
print " byte#93-94: ",$head_16bit_047, "\n ";
print " byte#95-96: ",$head_16bit_048, "\n ";
print " byte#97-98: ",$head_16bit_049, "\n ";
print " byte#99-100: ",$head_16bit_050, "\n ";
print " byte#101-102: ",$head_16bit_051, "\n ";
print " byte#103-104: ",$head_16bit_052, "\n ";
print " byte#105-106: ",$head_16bit_053, "\n ";
print " byte#107-108: ",$head_16bit_054, "\n ";
print " byte#109-110: ",$head_16bit_055, "\n ";
print " byte#111-112: ",$head_16bit_056, "\n ";
print " byte#113-114: ",$head_16bit_057, "\n ";
print " byte#115-116: ",$head_16bit_058, "\n ";
print " byte#117-118: ",$head_16bit_059, "\n ";
print " byte#119-120: ",$head_16bit_060, "\n ";
print " byte#121-122: ",$head_16bit_061, "\n ";
print " byte#123-124: ",$head_16bit_062, "\n ";
print " byte#125-126: ",$head_16bit_063, "\n ";
print " byte#127-128: ",$head_16bit_064, "\n ";
print " byte#129-130: ",$head_16bit_065, "\n ";
print " byte#131-132: ",$head_16bit_066, "\n ";
print " byte#133-134: ",$head_16bit_067, "\n ";
print " byte#135-136: ",$head_16bit_068, "\n ";
print " byte#137-138: ",$head_16bit_069, "\n ";
print " byte#139-140: ",$head_16bit_070, "\n ";
print " byte#141-142: ",$head_16bit_071, "\n ";
print " byte#143-144: ",$head_16bit_072, "\n ";
print " byte#145-146: ",$head_16bit_073, "\n ";
print " byte#147-148: ",$head_16bit_074, "\n ";
print " byte#149-150: ",$head_16bit_075, "\n ";
print " byte#151-152: ",$head_16bit_076, "\n ";
print " byte#153-154: ",$head_16bit_077, "\n ";
print " byte#155-156: ",$head_16bit_078, "\n ";
print " byte#157-158: ",$head_16bit_079, "\n ";
print " byte#159-160: ",$head_16bit_080, "\n ";
print " byte#161-162: ",$head_16bit_081, "\n ";
print " byte#163-164: ",$head_16bit_082, "\n ";
print " byte#165-166: ",$head_16bit_083, "\n ";
print " byte#167-168: ",$head_16bit_084, "\n ";
print " byte#169-170: ",$head_16bit_085, "\n ";
print " byte#171-172: ",$head_16bit_086, "\n ";
print " byte#173-174: ",$head_16bit_087, "\n ";
print " byte#175-176: ",$head_16bit_088, "\n ";
print " byte#177-178: ",$head_16bit_089, "\n ";
print " byte#179-180: ",$head_16bit_090, "\n ";
print " byte#181-182: ",$head_16bit_091, "\n ";
print " byte#183-184: ",$head_16bit_092, "\n ";
print " byte#185-186: ",$head_16bit_093, "\n ";
print " byte#187-188: ",$head_16bit_094, "\n ";
print " byte#189-190: ",$head_16bit_095, "\n ";
print " byte#191-192: ",$head_16bit_096, "\n ";
print " byte#193-194: ",$head_16bit_097, "\n ";
print " byte#195-196: ",$head_16bit_098, "\n ";
print " byte#197-198: ",$head_16bit_099, "\n ";
print " byte#199-200: ",$head_16bit_100, "\n ";
print " byte#201-202: ",$head_16bit_101, "\n ";
print " byte#203-204: ",$head_16bit_102, "\n ";
print " byte#205-206: ",$head_16bit_103, "\n ";
print " byte#207-208: ",$head_16bit_104, "\n ";
print " byte#209-210: ",$head_16bit_105, "\n ";
print " byte#211-212: ",$head_16bit_106, "\n ";
print " byte#213-214: ",$head_16bit_107, "\n ";
print " byte#215-216: ",$head_16bit_108, "\n ";
print " byte#217-218: ",$head_16bit_109, "\n ";
print " byte#219-220: ",$head_16bit_110, "\n ";
print " byte#221-222: ",$head_16bit_111, "\n ";
print " byte#223-224: ",$head_16bit_112, "\n ";
print " byte#225-226: ",$head_16bit_113, "\n ";
print " byte#227-228: ",$head_16bit_114, "\n ";
print " byte#229-230: ",$head_16bit_115, "\n ";
print " byte#231-232: ",$head_16bit_116, "\n ";
print " byte#233-234: ",$head_16bit_117, "\n ";
print " byte#235-236: ",$head_16bit_118, "\n ";
print " byte#237-238: ",$head_16bit_119, "\n ";
print " byte#239-240: ",$head_16bit_120, "\n ";
print "--------------First Trace Header 32 BIT-----------\n ";
print " byte# 1- 4: ",$head_32bit_001, "\n ";
print " byte# 5- 8: ",$head_32bit_003, "\n ";
print " byte# 9- 12: ",$head_32bit_005, "\n ";
print " byte# 13- 16: ",$head_32bit_007, "\n ";
print " byte# 17- 20: ",$head_32bit_009, "\n ";
print " byte# 21- 24: ",$head_32bit_011, "\n ";
print " byte# 25- 28: ",$head_32bit_013, "\n ";
print " byte# 29- 32: ",$head_32bit_015, "\n ";
print " byte# 33- 36: ",$head_32bit_017, "\n ";
print " byte# 37- 40: ",$head_32bit_019, "\n ";
print " byte# 41- 44: ",$head_32bit_021, "\n ";
print " byte# 45- 48: ",$head_32bit_023, "\n ";
print " byte# 49- 52: ",$head_32bit_025, "\n ";
print " byte# 53- 56: ",$head_32bit_027, "\n ";
print " byte# 57- 60: ",$head_32bit_029, "\n ";
print " byte# 61- 64: ",$head_32bit_031, "\n ";
print " byte# 65- 68: ",$head_32bit_033, "\n ";
print " byte# 69- 72: ",$head_32bit_035, "\n ";
print " byte# 73- 76: ",$head_32bit_037, "\n ";
print " byte# 77- 80: ",$head_32bit_039, "\n ";
print " byte# 81- 84: ",$head_32bit_041, "\n ";
print " byte# 85- 88: ",$head_32bit_043, "\n ";
print " byte# 89- 92: ",$head_32bit_045, "\n ";
print " byte# 93- 96: ",$head_32bit_047, "\n ";
print " byte# 97-100: ",$head_32bit_049, "\n ";
print " byte#101-104: ",$head_32bit_051, "\n ";
print " byte#105-108: ",$head_32bit_053, "\n ";
print " byte#109-112: ",$head_32bit_055, "\n ";
print " byte#113-116: ",$head_32bit_057, "\n ";
print " byte#117-120: ",$head_32bit_059, "\n ";
print " byte#121-124: ",$head_32bit_061, "\n ";
print " byte#125-128: ",$head_32bit_063, "\n ";
print " byte#129-132: ",$head_32bit_065, "\n ";
print " byte#133-136: ",$head_32bit_067, "\n ";
print " byte#137-140: ",$head_32bit_069, "\n ";
print " byte#141-144: ",$head_32bit_071, "\n ";
print " byte#145-148: ",$head_32bit_073, "\n ";
print " byte#149-152: ",$head_32bit_075, "\n ";
print " byte#153-156: ",$head_32bit_077, "\n ";
print " byte#157-160: ",$head_32bit_079, "\n ";
print " byte#161-164: ",$head_32bit_081, "\n ";
print " byte#165-168: ",$head_32bit_083, "\n ";
print " byte#169-172: ",$head_32bit_085, "\n ";
print " byte#173-176: ",$head_32bit_087, "\n ";
print " byte#177-180: ",$head_32bit_089, "\n ";
print " byte#181-184: ",$head_32bit_091, "\n ";
print " byte#185-188: ",$head_32bit_093, "\n ";
print " byte#189-192: ",$head_32bit_095, "\n ";
print " byte#193-196: ",$head_32bit_097, "\n ";
print " byte#197-200: ",$head_32bit_099, "\n ";
print " byte#201-204: ",$head_32bit_101, "\n ";
print " byte#205-208: ",$head_32bit_103, "\n ";
print " byte#209-212: ",$head_32bit_105, "\n ";
print " byte#213-216: ",$head_32bit_107, "\n ";
print " byte#217-220: ",$head_32bit_109, "\n ";
print " byte#221-224: ",$head_32bit_111, "\n ";
print " byte#225-228: ",$head_32bit_113, "\n ";
print " byte#229-232: ",$head_32bit_115, "\n ";
print " byte#233-236: ",$head_32bit_117, "\n ";
print " byte#237-240: ",$head_32bit_119, "\n ";

SEGY2ASCII

2011-11-12T17:42:00.001-08:00

Konversi SEGY ke ASCII diperlukan karena format ASCII merupakan bahasa manusia (bukan mesin) serta memiliki fleksibilitas untuk melakukan manipulasi matematika lanjut dengan software-software yang ada (seperti octave, matlab, shell script, dll.).

Di luar sana, banyak sekali beredar aplikasi atau tool gratis untuk melakukan konversi tersebut. Diantaranya SegyMAT yang berbasis matlab dan octave yang bisa didownload disini.

Anda akan memperoleh file SegyMAT.tgz, lalu ekstraklah (pada Linux) dengan perintah tar -xvvf SegyMAT.tgz sehingga diperoleh beberapa file dengan ektensi .m.

Untuk mendemonstrasikannya pada data segy, downloadlah data density.segy model Marmousi yang telah saya jelaskan disini.

Lakukan konversi segy ke su dengan perintah:
segyread tape=density.segy verbose=1 endian=0 > density.su

surange < density.su
751 traces:
tracl 752 1502 (752 - 1502)
tracr 752 1502 (752 - 1502)
fldr 1
tracf 752 1502 (752 - 1502)
trid 1
nvs 1
scalco -100
ns 2301
dt 2304

Tampilkan dengan perintah:
suflip flip=0 < density.su | suximage legend=1 perc=90 &

Selanjutnya saya akan menggunakan fungsi ReadSegy.m yang telah ekstrak dari SegyMAT. Simpanlah semua file .m serta data segy pada direktori atau path yang sama, lalu dari direktori tersebut panggillah octave. Jika belum memiliki octave installah dengan sudo get -apt octave

Setelah memanggil octave dengan mengetikkan octave pada terminal linux, anda akan memperoleh octave>

Ekstrak segy dengan perintah:

[density,TraceHeaderInfo,SegyTraceHeaders,SegyHeader]=ReadSegy('density.segy');

Lalu mendapatkan pesan:
Number of Samples Per Trace=2301 (banyaknya sampel per trace)
Number of Traces=751 (banyaknya trace)
Reading trace 200 of 751
Reading trace 400 of 751
Reading trace 600 of 751
header=200 of 751
header=400 of 751
header=600 of 751

Jika anda mengetikkan density pada octave, anda akan memperoleh hal seperti di bawah ini:

density =

Columns 1 through 5:

1000.0 1000.0 1000.0 1000.0 1000.0
1000.0 1000.0 1000.0 1000.0 1000.0
1000.0 1000.0 1000.0 1000.0 1000.0
1000.0 1000.0 1000.0 1000.0 1000.0
1000.0 1000.0 1000.0 1000.0 1000.0
1000.0 1000.0 1000.0 1000.0 1000.0
1000.0 1000.0 1000.0 1000.0 1000.0

Sekarang saya akan men-save density menjadi density.txt saveascii

saveascii(density,'density.txt',1);

Buka terminal linux yang baru, lalu masuklah ke direktori ketika anda memanggil octave.

Pada direktori tersebut anda memiliki density.txt dimana jumlah kolom menunjukkan banyaknya trace dan jumlah baris menunjukkan sampel.

Buatlah matrix kedalaman dengan perintah:
seq 0 2.304 5300 > depth.txt
Sisipkan depth.txt ke-kolom pertama density, dengan perintah paste:

paste depth.txt density.txt > depthdensity.ascii

Coba plot dengan gnuplot:
set xrange [1800:2500]
set yrange [0:5300] reverse
set xlabel 'Densitas (kg/m^3)'
set ylabel 'Kedalaman (ft)'
plot "depthdensity.ascii" using 300:1 title 'density trace#300' with lines

Kembali ke terminal octave, berdasarkan data density saya akan menghitung overburden (psi), lalu outputnya disimpan sebagai segy dengan perintah WriteSegy.

depth=repmat(2.304,2301,751);
depthdens=depth.*(density./1000);
overburden= 0.433.*cumsum(depthdens);
saveascii(overburden,'overburden.ascii',1);
WriteSegy('overburden.segy',overburden);

Keluar dari octave (dengan mengetikkan exit), konversi dari segy ke su:
segyread tape=overburden.segy verbose=1 endian=0 > overburden.su

surange < overburden.su
751 traces:
ns 2301
dt 2000
year 2011
hour 20
minute 34 35 (34 - 35)
sec 0 60 (46 - 2)

Perbaiki sampling rate menjadi 2304
sushw key=dt a=2304 < overburden.su > overburden_new.su

Tampilkan:
suximage < overburden_new.su perc=90 legend=1 title='overburden(psi)' &
Tekan 'h' untuk mengubah skala warna.

gnuplot

2011-11-08T04:22:00.000-08:00

Plotting merupakan pekerjaan yang tidak bisa dipisahkan dari seorang geoscientist.
Sebagai profesional anda dapat melakukan plotting pada software komersial seperti Petrel, Geoframe, Omega2-Attribute Display atau MS-Excel.
Untuk keperluan quick-qc, Unix menyediakan berbagai macam perangkat plotting gratis diantaranya xgraph, ygraph, gnuplot, dll.

Pada bagian ini, saya akan menjelaskan penggunaan gnuplot untuk memplot data dengan format XYZ (data well, peta struktur, gravity, dll).

Anda dapat men-check apakah anda memiliki gnuplot atau tidak dengan which gnuplot atau man gnuplot
Jika anda belum memiliki, anda bisa mengistallnya dengan perintah sudo apt-get install gnuplot-x11 (untuk Ubuntu 10.10).

Berikut adalah contoh plotting data.txt

head -5 data.txt
0 0.0000000000 1.000000000
1 0.0174524064 0.999847695
2 0.0348994967 0.999390827
3 0.0523359562 0.998629534
4 0.0697564737 0.997564050

Panggil gnuplot dengan perintah gnuplot, sehingga diperoleh gnuplot>

set xlabel 'XLABEL'
set ylabel 'YLABEL'
set xrange [0:360]
set yrange [-1:1]
plot "data.txt" using 1:2 title 'SIN' pointtype 3 pointsize 1 , "data.txt" using 1:3 title 'COS' pointtype 6 pointsize 1

Jika saya copy code di atas pada terminal gnuplot (gnuplot>), maka saya kan memperoleh plot dibawah ini.

set xlabel 'XLABEL' membuat label X, set ylabel 'YLABEL' membuat label Y, set xrange [0:360] membatasi range X dari 0 s/d 360, set yrange [-1:1] membatasi rentang Y. using 1:2 (menggunakan kolom 1 untuk X dan kolom 2 untuk Y), pointtype dan pointsize menunjukkan tipe point dan ukurannya. Anda bisa mencoba berbagai tipe 1 s/d 7 atau size dari 0.1 s/d 6 (misalnya).

Untuk memperoleh plot garis, saya menggunakan with lines.
plot "data.txt" using 1:2 title 'SIN' with lines, "data.txt" using 1:3 title 'COS' with lines

gnuplot dapat dipanggil melalui shell (tanpa harus berada di terminal gnuplot).
Berikut contoh shell untuk plotting data di atas.

#!/bin/sh
gnuplot -persist <<PLOT
set xlabel 'XLABEL'
set ylabel 'YLABEL'
set xrange [0:360]
set yrange [-1:1]
plot 'data.txt' using 1:2 title 'SIN' with lines, \
'data.txt' using 1:3 title 'COS' with lines
quit
PLOT

Berikutnya, saya akan menunjukkan memplot peta atau peta 3D dengan perintah splot.

Untuk data grav.xyz anda dapat mengestraknya di sini.

Sehingga saya memperoleh:
head -8 grav.xyz
100.0083 -0.0082 36.20
100.0250 -0.0082 30.80
100.0417 -0.0082 26.40
100.0583 -0.0082 23.90
100.0750 -0.0082 23.80
100.0917 -0.0082 25.50
100.1083 -0.0082 28.40
100.1250 -0.0082 31.20

Panggil gnuplot, lalu copas (copy-paste) kode berikut pada terminal gnuplot

set key off
set view 0,0
set xrange [100:102]
set yrange [-2:0]
set xlabel 'X-Coordinate'
set title "Gravity Map"
set palette model RGB rgbformulae 7,5,15
splot 'grav.xyz' using 1:2:3 with points palette pointsize 1 pointtype 7

Anda bisa merotasinya dengan meng-klik icon dua-panah yang menyerupai yin-yang, sehingga diperoleh gambar 3-D.

gnuplot memiliki fleksibilitas yang tak kalah dari matlab. Berikut contohnya:

Berikut adalah kode untuk gambar di atas:
unset border
set dummy u,v
unset key
set palette model HSV rgbformulae 3,2,2
set parametric
set view 60, 30, 1.1, 1.33
set isosamples 50, 20
set noxtics
set noytics
set noztics
set title "Interlocking Tori - PM3D surface with no depth sorting"
set urange [ -3.14159 : 3.14159 ] noreverse nowriteback
set vrange [ -3.14159 : 3.14159 ] noreverse nowriteback
set xrange [ * : * ] noreverse nowriteback
set yrange [ * : * ] noreverse nowriteback
set zrange [ * : * ] noreverse nowriteback
set cbrange [ * : * ] noreverse nowriteback
set pm3d scansbackward
set colorbox vertical origin screen 0.9, 0.2, 0 size screen 0.05, 0.6, 0 front bdefault
f(x,y) = sin(-sqrt((x+5)**2+(y-7)**2)*0.5)
splot cos(u)+.5*cos(u)*cos(v),sin(u)+.5*sin(u)*cos(v),.5*sin(v) with pm3d, 1+cos(u)+.5*cos(u)*cos(v),.5*sin(v),sin(u)+.5*sin(u)*cos(v) with pm3d
#courtesy: gnuplot.sourceforge.net

Unix untuk Geoscientists

2011-11-07T06:51:00.000-08:00

Unix memiliki ratusan fungsi yang sangat berguna untuk mendukung pekerjaan seorang geoscientist. Tetapi pada kenyataannya, hanya sebagian kecil yang termanfaatkan.

Hal ini mungkin karena terbatasnya pengetahuan, keenganan atau beranggapan bahwa Unix sebagai suatu yang rumit.

Sebagian perintah Unix telah saya jelaskan di sini. Untuk memudahkan memahami perintah Unix yang lainnya, saya jelaskan berdasarkan contoh-contoh berikut:

man
manual dari perintah sebuah fungsi, untuk keluar dari man tekan q

man tee
NAME
tee - read from standard input and write to standard output and files

SYNOPSIS
tee [OPTION]... [FILE]...

Buatlah dua file yang berisi angka-angka dengan perintah seq. Pada terminal Linux ketik seq 1 1 5 > file1 lalu enter dan seq 2 2 10 > file2 lalu enter.

Cek kedua file
more file1
1
2
3
4
5

more file2
2
4
6
8
10

paste
menggabungkan dua file
paste file1 file2 > file3
more file3
1 2
2 4
3 6
4 8
5 10

cut
mengekstrak bagian dari data i.e. kolom, karakter, dll.
cut -f 1,2 file3
1 2
2 4
3 6
4 8
5 10

cut -f 2 file3
2
4
6
8
10

sort
men-sorting baris-baris dari file
sort -n file1
1
2
3
4
5

sort -r file1
5
4
3
2
1

Buat file4 sbb:
more file4
1 5
2 4
3 3
4 2
5 1

sort -k 2 file4
5 1
4 2
3 3
2 4
1 5

wc
menghitung jumlah baris (-l), jumlah kata (-w) dan jumlah karakter (-c)
wc -l file4
5 file4 .....ada 5 baris dalam file4

cat
concatenate file
cat file1 file2 > file5

more file5
1
2
3
4
5
2
4
6
8
10

uniq
membuang atau memprint nilai yang sama

more file5
1
2
2
3
4
4
5
6
8
10

uniq -u file5
1
3
5
6
8
10

uniq -d file5
2
4

grep
print-out berdasarkan karakter atau kata tertentu
more file5
1
2
2
3
4
4
5
6
8
10

grep '2' file5
2
2

sed -e 's/2/20/g' file5
1
20
20
3
4
4
5
6
8
10

sed '5,10 w junk1' file5 (perintah print baris ke 5 s/d 10 dari file5).
more junk1
4
4
5
6
8
10

tac
mengurutkan secara terbalik
tac file5
10
8
6
5
4
4
3
2
2
1

head
menampilkan bagian awal sebuah file
head -3 file5
1
2
2

tail
menampilkan bagian akhir dari file
tail -3 file5
6
8
10

fold
membatasi jumlah karakter perbaris

contoh berikut membatasi hanya 12 karakter dari file yang bernama test
fold -w 12 test
Fold ' is a
Unix command
used for ma
king a file
with long li

tee
printout di screen dan file secara sekaligus

which
adalah perintah untuk mengetahui lokasi perintah yang bisa dieksekusi

which surange
/home/akmal/seismic_unix/bin/surange

history
printout perintah yang telah dipakai di terminal
history > junk
more junk
1 gv modelfile.eps
2 sh psmerge1a.sh
3 ls
4 gv psmerge1a.eps
5 gv modelfile.eps
6 sh raytracing.sh

Kombinasi berbagai perintah
Contoh: hanya meng-output-kan satu kali dari duplikasi yang terjadi pada kolom2

more input.txt
1 3 2
3 1 7
5 3 4
2 1 1
3 1 4
5 3 6
2 3 8
2 2 8
2 2 5

more output.txt
2 1 1
2 2 5
1 3 2

sort -k 2n < input.txt| awk -F\ '{ print NR, $0, count[$2]++ }' | awk '$5 == 0 '| awk '{print $2,$3,$4}' > output.txt

sort -k 2n < input.txt
2 1 1
3 1 4
3 1 7
2 2 5
2 2 8
1 3 2
2 3 8
5 3 4
5 3 6

sort -k 2n < input.txt| awk -F\ '{ print NR, $0, count[$2]++ }'
1 2 1 1 0
2 3 1 4 1
3 3 1 7 2
4 2 2 5 0
5 2 2 8 1
6 1 3 2 0
7 2 3 8 1
8 5 3 4 2
9 5 3 6 3

sort -k 2n < input.txt| awk -F\ '{ print NR, $0, count[$2]++ }' | awk '$5 == 0 '
1 2 1 1 0
4 2 2 5 0
6 1 3 2 0

sort -k 2n < input.txt| awk -F\ '{ print NR, $0, count[$2]++ }' | awk '$5 == 0 '| awk '{print $2,$3,$4}' >output.txt

more output
2 1 1
2 2 5
1 3 2

Advanced SU Part-9: Synthetic Gathers

2011-10-31T07:31:00.001-07:00

Pada Seismic Unix, gathers seismik sintetik dapat diperoleh dengan menggunakan perintah triseis.

Input utama untuk triseis adalah model kecepatan. Model kecepatan yang akan saya gunakan adalah model karbonat, dimana proses pembuatannya telah saya jelaskan di sini.

Untuk membuat jumlah reflektor yang lebih banyak, saya memodifikasi jumlah lapisan sehingga pada model ini saya memiliki reservoir karbonat dengan beberapa lapisan lain di atasnya.

Berikut adalah Kode Bourne Shell untuk model di atas:

#!/bin/sh
trimodel xmin=0 zmin=0 xmax=10.0 zmax=4.0 \
1 xedge=0,10 \
zedge=0,0 \
sedge=0,0 \
2 xedge=0,10 \
zedge=0.5,0.3 \
sedge=0,0 \
3 xedge=0,10 \
zedge=1,0.4 \
sedge=0,0 \
4 xedge=0,10 \
zedge=1.5,0.5 \
sedge=0,0 \
5 xedge=0,1,2,3,4,5,6,7,8,9,10 \
zedge=2.5,2.5,2.4,2.3,1.8,1.2,1.1,1.0,1.2,2.5,2.5 \
sedge=0,0,0,0,0,0,0,0,0,0,0 \
6 xedge=0,10 \
zedge=3,3 \
sedge=0,0 \
7 xedge=0,10 \
zedge=4,4 \
sedge=0,0 \
sfill=0.1,0.3,0,0,0.44,0,0 \
sfill=0.1,0.7,0,0,0.36,0,0 \
sfill=0.1,1.3,0,0,0.32,0,0 \
sfill=0.1,2.0,0,0,0.37,0,0 \
sfill=0.1,2.7,0,0,0.06,0,0 \
sfill=0.1,3.5,0,0,0.11,0,0 \
kedge=1,2,3,4,5,6,7 \
>model.bin

spsplot <model.bin >model.ps \
title="Carbonate Model" \
labelz="Depth [km]" labelx="Distance [km]" \
labelsize=10 \
gedge=1.0 gtri=2.0 \
gmin=0.2 \
gmax=0.8 \
wbox=6.0 hbox=2.0 &

gv model.ps &

Untuk melakukan shooting pada model di atas, saya memilih jumlah shot gathers (konfigurasi split spread):170, nangle=101, fangle=-65, langle=65, nt=1200 (jumlah sampel) dan dt=0.004 (sampling rate). Shot interval: 50m, receiver interval:50m, Offset maksimum:1475m dan minimum -1475m. Shot pertama berada pada lokasi 1.5 dan shot terakhir pada lokasi 9.95km. Karena offfset minimum adalah 1.475, maka lokasi receiver pertama adalah 1.5-1.475=0.025 dan receiver terakhir (seharusnya) 9.95+1.475=11.425. Namun saya membatasi sampai 10km, sehingga anda akan melihat shot gather yang ke 170 tidak akan utuh.

Berikut ini adalah kode shooting data seismik untuk model tersebut.

#!/bin/sh
##modified from mines.edu
/bin/rm -f tmp*
nangle=101
fangle=-65
langle=65
nt=1200 dt=0.004
datafile=model.bin
seismic=carbonate.su

i=0
while [ "$i" -ne "170" ]
do

fs=`echo "$i * 0.05" | bc -l`
sx=`echo "$i * 50" | bc -l`
fldr=`echo "$i + 1" | bc -l`

j=0
while [ "$j" -ne "60" ]
do

fg=`echo "$i * 0.05 + $j *0.05" | bc -l`
gx=`echo "$i * 50 + $j * 50 -1475" | bc -l`
offset=`echo "$j * 50 -1475" | bc -l`
tracl=`echo "$i * 60 + $j + 1" | bc -l`
tracf=`echo "$j + 1" | bc -l`

echo "sx=$sx gx=$gx trace_number=$tracl"

k=2
while [ "$k" -ne "7" ]
do

triseis <$datafile xs=1.5,9.95 zs=0,0 \
xg=0.025,10 zg=0,0 \
nangle=$nangle fangle=$fangle langle=$langle \
kreflect=$k krecord=1 fpeak=12 lscale=0.5 \
ns=1 fs=$fs ng=1 fg=$fg nt=$nt dt=$dt |
suaddhead nt=$nt |
sushw key=dt,tracl,tracr,fldr,tracf,trid,offset,sx,gx \
a=4000,$tracl,$tracl,$fldr,$tracf,1,$offset,$sx,$gx >> tmp$k
k=`expr $k + 1`

done
j=`expr $j + 1`

done
i=`expr $i + 1`

done

susum tmp2 tmp3 >tmpa
susum tmpa tmp4 >tmpb
rm -f tmpa
susum tmpb tmp5 >tmpa
rm -f tmpb
susum tmpa tmp6 >$seismic
rm -f tmpa
rm tmp*

exit

Jika perintah tersebut dijalankan, maka anda akan memperoleh sebuah file seismic gathers carbonate.su

surange < carbonate.su
10200 traces:
tracl 1 10200 (1 - 10200)
tracr 1 10200 (1 - 10200)
fldr 1 170 (1 - 170)
tracf 1 60 (1 - 60)
trid 1
offset -1475 1475 (-1475 - 1475)
sx 0 8450 (0 - 8450)
gx -1475 9925 (-1475 - 9925)
ns 1200
dt 4000

Berikut adalah ke 120
suwind < carbonate.su key=fldr min=120 max=120 | suxwigb perc=95 &

Berikut adalah ke 170 (tidak utuh).
suwind < carbonate.su key=fldr min=170 max=170 | suxwigb perc=95 &

Gunakan perintah berikut untuk quick QC (menampilkan near trace ) untuk masing-masing gather.
suwind key=tracf min=1 max=1 < carbonate.su | suximage perc=95 &

Untuk memperoleh penampang seismik yang utuh, maka anda harus melakukan analisa kecepatan, NMO dan stack. Saya telah menjelaskan semua caranya di sini atau di sini.

Anda mungkin bertanya mengapa saya tidak memperoleh citra flank carbonate sebelah kanan? Jawabannya adalah karena geometri pengambilan data seismic yang saya lakukan tidak memiliki aperture yang memadai. Hal ini bisa di-QC dengan melakukan raytracing pada ujung lokasi sebelah kanan, terlihat bahwa raypath yang terefleksikan oleh flank carbonate tersebut tidak ditangkap near offset, namun jauh terpantulkan ke sebelah kanan.

Jika saya memperlebar aperture sampai 20 km ke kanan, maka dari modeling ray tracing di bawah ini, gelombang yang terpantulkan oleh flank sebelah kanan akan terekam dengan baik.

Berikut adalah model karbonat di atas dengan aperture yang lebih lebar.

Berikut adalah hasil perekaman, jika saya memperlebar sampai 20km.

Dari hasil ini kita melihat bagaimana modeling memiliki peranan yang sangat penting di dalam akuisisi data seismik, jangan sampai jutaan dollar terbuang sia-sia karena kita salah memilih parameter akuisisi.

Advanced SU Part-8: Ray Tracing

2011-10-29T02:51:00.001-07:00

Setelah kita melakukan pemodelan kecepatan, pada bagian ini saya akan mendemonstrasikan cara melakukan ray tracing (triray) pada model tersebut.

Pada desain survey seismik, ray tracing penting dilakukan untuk mempelajari perilaku jejak gelombang seismik, derajat iluminasi, serta penentuan parameter lapangan seperti panjang kabel, geometri, dll.

Gambar di bawah ini menunjukkan jejak gelombang dari 3 shots yang direfleksikan oleh horison 2 (merah), horison 3 (hijau) dan horizon 4 (biru).

Berikut adalah kode Bourne Shell yang digunakan untuk menghasilkan gambar di atas.

#! /bin/sh
set -x

triray < modelfile.bin > rayend1.bin rayfile=ray1.bin \
nangle=20 fangle=-75 langle=75 \
xs=8 zs=0.0 nxz=600 \
refseq=2,1

triray < modelfile.bin > rayend2.bin rayfile=ray2.bin \
nangle=22 fangle=-25 langle=25 \
xs=21 zs=0.0 nxz=600 \
refseq=2,0 refseq=3,1

triray < modelfile.bin > rayend3.bin rayfile=ray3.bin \
nangle=20 fangle=-15 langle=15 \
xs=27 zs=0.0 nxz=600 \
refseq=2,0 refseq=3,0 refseq=4,1

psgraph < ray1.bin > ray1.eps \
nplot=`cat outpar` n=600 hbox=4.0 wbox=10.0 \
x1beg=0.0 x1end=4.0 x2beg=0 x2end=30 \
d1num=0.5 d2num=1.0 style=seismic linecolor=red \

psgraph < ray2.bin > ray2.eps \
nplot=`cat outpar` n=600 hbox=4.0 wbox=10.0 \
x1beg=0.0 x1end=4.0 x2beg=0 x2end=30 \
d1num=0.5 d2num=1.0 style=seismic linecolor=green \

psgraph < ray3.bin > ray3.eps \
nplot=`cat outpar` n=600 hbox=4.0 wbox=10.0 \
x1beg=0.0 x1end=4.0 x2beg=0 x2end=30 \
d1num=0.5 d2num=1.0 style=seismic linecolor=blue \

spsplot < modelfile.bin> modelfile.eps \
gedge=0.5 gtri=2.0 gmin=0 gmax=1 \
labelz="Depth (km)" labelx="Distance (km)" \
dxnum=1.0 dznum=0.5 wbox=10 hbox=4

psmerge in=modelfile.eps in=ray1.eps in=ray2.eps in=ray3.eps > output.eps

gv output.eps
exit

Kode di atas memiliki 3 paket perintah triray, dimana triray pertama digunakan untuk menghasilkan jejak gelombang pertama (merah). modelfile.bin adalah input model kecepatan yang telah saya jelaskan disini. rayend1.bin rayfile=ray1.bin adalah output binary. nangle=20 adalah jumlah sudut tembak (20), fangle adalah sudut tembak pertama -75 derajat dan langle adalah sudut tembak terakhir 75 derajat. xs dan zs adalah lokasi shot (jarak dan kedalaman), nxz=600 adalah jumlah pasangan (x,z) pada rayfile. refseq=2,1 menunjukkan bahwa gelombang harus terefleksikan oleh horison ke 2. Pelajari parameter yang saya pakai untuk masing-masing shot. psgraph adalah merubah format bin ke eps. psmerge adalah menggabungkan semua file eps.

Advanced SU Part-7: Velocity Model

2011-10-28T19:40:00.000-07:00

Seismic Unix memiliki kapabilitas untuk membuat model kecepatan. Pekerjaan ini dilakukan dengan menggunakan perintah trimodel.

Gambar di bawah ini adalah contoh model kecepatan untuk sebuah reservoir karbonat dimana windward di sebelah kanan (ditandai dengan kemiringan yang curam) serta leeward di sebelah kiri (kemiringan landai). Di atasnya terdapat lapisan klastik dan di bawahnya terdapat base karbonat serta basement dan/atau source rocks.

Berikut adalah kode Bourne Shell yang digunakan untuk menghasilkan gambar di atas:

1 #! /bin/sh
2 set -x
3
4 trimodel xmin=0 xmax=30 zmin=0 zmax=4 \
5 1 xedge=0,30 \
6 zedge=0,0 \
7 sedge=0,0 \
8 2 xedge=0,12,24,30 \
9 zedge=1,0.6,0.3,0 \
10 sedge=0,0,0,0 \
11 3 xedge=0,3,6,9,12,15,18,21,24,27,30 \
12 zedge=2.5,2.3,2.1,1.5,1.5,1.3,1.2,1.1,2.1,2.3,2.5 \
13 sedge=0,0,0,0,0,0,0,0,0,0,0 \
14 4 xedge=0,12,24,30 \
15 zedge=3,3,3,3 \
16 sedge=0,0,0,0 \
17 5 xedge=0,30 \
18 zedge=4,4 \
19 sedge=0,0 \
20 kedge=1,2,3,4,5 \
21 sfill=0,0.5,0,0,0.44,0,0 \
22 sfill=0,2.0,0,0,0.16,0,0 \
23 sfill=0,2.7,0,0,0.06,0,0 \
24 sfill=0,3.5,0,0,0.11,0,0 > modelfile.bin
25 spsplot < modelfile.bin > modelfile.eps \
26 gedge=0.5 gtri=2.0 gmin=0 gmax=1 \
27 labelz="Depth (km)" labelx="Distance (km)" \
28 dxnum=1.0 dznum=0.5 wbox=10 hbox=4
29
30 gv modelfile.eps
31 exit

Baris ke-4 memberikan batasan model dengan sumbu X antara 0 s/d 30km dan kedalaman 0 s/d 4km.
Baris 5-7: membuat permukaan bumi yang datar antara 0 s/d 30 dengan kedalaman 0 dan gradient sloth(1/kecepatan²)
Baris 8-10: membuat horizon kedua, dst. sampai baris ke 17-19 untuk horizon ke-5.
Baris 20: menunjukkan jumlah interface (batas).
Baris 21-24: mempopulasi masing-masing lapisan dengan nilai sloth.
Dimana masing-masing angka mencerminkan sfill=x,z,x0,z0,sloth,ds/dx,ds/dz
Baris 23 menunjukkan x=0 dan z=2.7 yang mempopulasi nilai karbonat dengan sloth=0.06 (4000m/s).
Baris 25-28: memplot model binary menjadi eps, baris 30 menampilkan file eps dengan perintah gv.
Jika anda belum memiliki gv, maka install dengan perintah sudo apt-get install gv

Untuk menampilkan model binary, gunakan perintah berikut:
sxplot < modelfile.bin tricolor=none cmap=hsv1 label1=depth edgecolor=red

AWK Part-3: Advanced & Miscellaneous

2011-10-06T07:40:00.000-07:00

1. Mengulangi setiap baris beberapa kali (4 kali)
more input.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3

more output.txt
2.3 3.2 3.6
2.3 3.2 3.6
2.3 3.2 3.6
2.3 3.2 3.6
4.2 1.3 6.5
4.2 1.3 6.5
4.2 1.3 6.5
4.2 1.3 6.5
5.1 4.2 9.3
5.1 4.2 9.3
5.1 4.2 9.3
5.1 4.2 9.3

awk '{ for (i = 1; i <=4 ; i++) print }' input.txt > output.txt

Dengan Bourne Shell

#!/bin/sh
max=4
while read line
do
n=1
while [ $n -le $max ]
do
echo "$line"
n=$((n+1))
done
done < input.txt > output.txt

2. Mengulangi setiap baris sejumlah bilangan yang berada pada kolom pertama

more input.txt
3 2.3 3.2 3.6
2 4.2 1.3 6.5
5 5.1 4.2 9.3

Untuk contoh ini, baris pertama akan diulangi 3 kali, baris kedua 2 kali dan baris ketiga 5 kali.

more output.txt
3 2.3 3.2 3.6
3 2.3 3.2 3.6
3 2.3 3.2 3.6
2 4.2 1.3 6.5
2 4.2 1.3 6.5
5 5.1 4.2 9.3
5 5.1 4.2 9.3
5 5.1 4.2 9.3
5 5.1 4.2 9.3
5 5.1 4.2 9.3

Copy kode berikut lalu berinama apasaja.sh

#!/bin/sh
awk 'BEGIN{
}
{
for(i=0;i<$1;i++) print
}
END {
}' $1

Ubah permission dengan: chmod +x apasaja.sh (atau chmod 777 filename)
dengan ls -l :
sebelum chmod +x
-rw-r--r-- 1 agus agus 64 2011-10-06 05:23 apasaja.sh
setelah chmod +x
-rwxr-xr-x 1 agus agus 64 2011-10-06 05:23 apasaja.sh

Untuk menjalankan kode di atas:
sh apasaja.sh input.txt > output.txt

3. Membandingkan dua file dengan kolom pertama dan kolom kedua sebagai acuan. Baris dari file kedua di-outputkan jika kolom pertama dan kedua pada file1 sama dengan kolom pertama dan kedua pada file2.

more file1
-4.2 1.0 18.12 13.78 0.36 0.56 0.001
-3.8 9.0 24.02 18.91 0.41 0.48 0.001
-0.5 3.3 19.06 16.84 0.29 0.12 0.003

more file2
-3.8 9.0 -1.78172 11.3782 0.01 0.01 0.01 31704
-0.5 3.3 -1.83964 11.2841 0.01 0.01 0.01 31701
-2.3 3.6 -1.82152 11.3152 0.01 0.01 0.01 31702
-7.2 1.0 -1.80205 11.3466 0.01 0.01 0.01 31703
-2.0 3.6 -1.76099 11.4199 0.01 0.01 0.01 31705
29.9 3.2 -1.74029 11.4417 0.01 0.01 0.01 31706

more output.txt
-3.8 9.0 -1.78172 11.3782 0.01 0.01 0.01 31704
-0.5 3.3 -1.83964 11.2841 0.01 0.01 0.01 31701

awk 'FNR==NR{A[$1,$2]=1;next} A[$1,$2]' file1 file2 > output.txt

4. Memilih elemen dan menjadikannya sebagai variabel.

more input.txt
3 2.3 3.2 3.6
2 4.2 1.3 6.5
5 5.1 4.2 9.3

Copy kode berikut lalu berinama kode.sh

#!/bin/sh
angka=`awk 'FNR == 2 {print $3}' < input.txt`
echo elemen tersebut adalah: $angka
kali=`echo "$angka * 2.0" | bc -l`
echo elemen tersebut dikali 2 = $kali

sh kode.sh
elemen tersebut adalah: 1.3
elemen tersebut dikali 2 = 2.60

Perhatikan bahwa simbol acute (`) mengapit semua perintah yang mengoutputkan variabel.
Perintah bc -l di atas digunakan untuk memanggil basic calculator (bc) dengan menggunakan long precision (-l).
Anda dapat mencoba mengetikkan bc -l pada terminal linux untuk menjadikannya sebagai kalkulator.

5. Menggabungkan beberapa proses, contoh berikut adalah memilih baris ganjil lalu mengkalikannya dengan 0.001

awk 'NR%2 ' input.txt | awk '{for(i=1;i <=NF;++i)$i*=0.001;print}' > output.txt

more input.txt
3 2.3 3.2 3.6
2 4.2 1.3 6.5
5 5.1 4.2 9.3

more output.txt
0.003 0.0023 0.0032 0.0036
0.005 0.0051 0.0042 0.0093

6. Berbagai perintah Unix yang berguna untuk Geoscientist
more input.txt
3 2.3 3.2 3.6
2 4.2 1.3 6.5
5 5.1 4.2 9.3

menampilkan 2 baris pertama dari sebuah file: head -2 input.txt
3 2.3 3.2 3.6
2 4.2 1.3 6.5

menampilkan 2 baris terakhir dari sebuah file: tail -2 input.txt
2 4.2 1.3 6.5
5 5.1 4.2 9.3

7. ps -ef | grep agus menampilkan semua proses lalu menseleksinya yang mengandung nama agus

agus 1959 1 0 04:39 ? 00:00:17 mono /us
agus 1969 1874 2 04:41 pts/1 00:03:15 gedit program1.sh
agus 3999 1874 0 07:11 pts/1 00:00:00 ps -ef

Lalu saya ingin meng-kill program gedit: kill -9 1969

8. Untuk melihat aktifitas cpu: top
Perintah top cukup penting untuk memonitor proses berat seperti migrasi data seismik

agus@agus-Satellite-M305:~/programming$ top

Swap: 6416380k total, 0k used, 6416380k free, 634416k cached
PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND
1862 agus 20 0 246m 56m 22m S 16 1.5 18:44.44 plugin-containe
994 root 20 0 83804 38m 15m S 13 1.0 12:14.39 Xorg
1625 agus 20 0 91760 14m 10m S 9 0.4 0:43.72 gnome-terminal
1402 agus 20 0 157m 5680 4472 S 8 0.1 5:26.83 pulseaudio

plugin-containe memakan 16% dari CPU

9. mengurutkan file dari yang paling kecil sampai besar: ls -alSr -h

10. Mengecek kapasitas hard disk df -k -h .

agus@agus-Satellite-M305:~/programming$ df -k -h .
Filesystem Size Used Avail Use% Mounted on
/dev/sda6 56G 13G 40G 24

11. Menghitung jumlah baris dalam sebuah file: wc -l filename
12. Melihat isi file: more filename

AWK-Part 2: Matrix Manipulation

2011-10-01T03:56:00.000-07:00

1. Transpose

more input.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3

more output.txt
2.3 4.2 5.1
3.2 1.3 4.2
3.6 6.5 9.3

awk '{
for (f = 1; f <= NF; f++)
a[NR, f] = $f
}
NF > nf { nf = NF }
END {
for (f = 1; f <= nf; f++)
for (r = 1; r <= NR; r++)
printf a[r, f] (r==NR ? RS : FS)
}' input.txt > output.txt

Jika dengan Bourne Shell, copy kode berikut lalu berinama transpose,
Jalankan dengan perintah sh transpose input.txt > output.txt

#! /bin/sh
exec awk '
NR == 1 {
n = NF
for (i = 1; i <= NF; i++)
row[i] = $i
next
}
{
if (NF > n)
n = NF
for (i = 1; i <= NF; i++)
row[i] = row[i] " " $i
}
END {
for (i = 1; i <= n; i++)
print row[i]
}' ${1+"$@"}

2. Mengurangi Matrix dengan bilangan (0.4)
awk '{for(i=1;i <=NF;++i)$i-=0.4;print}' input.txt > output.txt

more input.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3

more output.txt
1.9 2.8 3.2
3.8 0.9 6.1
4.7 3.8 8.9

3. Menambahkan Matrix dengan bilangan (0.8)
awk '{for(i=1;i <=NF;++i)$i+=0.8;print}' input.txt > output.txt

more output.txt
3.1 4 4.4
5 2.1 7.3
5.9 5 10.1

4. Mengkalikan matrix dengan bilangan (10)
awk '{for(i=1;i <=NF;++i)$i*=10;print}' input.txt > output.txt

Untuk membagi dengan angka 10:
awk '{for(i=1;i <=NF;++i)$i/=10;print}' input.txt > output.txt

5. Pengurangan antara dua matrix
output.txt=matrix2-matrix1

more matrix1
3 3 3
3 3 3
3 3 3

more matrix2
1 1 1
1 1 1
1 1 1

awk '{n1=split($0,m1);getline <"matrix1";for(i=1;i <n1;++i)printf "%d ",m1[i]-$i;printf "%d\n",m1[n1]-$(n1)}' matrix2 > output.txt

more output.txt
-2 -2 -2
-2 -2 -2
-2 -2 -2

Modified from www.tek-tips.com

6. Menjumlahkan dua matrix
awk '{n1=split($0,m1);getline <"matrix1";for(i=1;i <n1;++i)printf "%d ",m1[i]+$i;printf "%d\n",m1[n1]+$(n1)}' matrix2

7. Sorting Matrix: ascending, dengan kolom ke 3 sebagai acuan

sort -k 3n < input.txt > output.txt

more input.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3
4.8 6.4 1.2
6.1 2.4 5.3
8.8 8.4 6.1
8.1 0.6 8.7

more output.txt
4.8 6.4 1.2
2.3 3.2 3.6
6.1 2.4 5.3
8.8 8.4 6.1
4.2 1.3 6.5
8.1 0.6 8.7
5.1 4.2 9.3

8. Sorting Matrix: descending, dengan kolom ke 3 sebagai acuan
sort -k 3nr < input.txt > output.txt

more output.txt
5.1 4.2 9.3
8.1 0.6 8.7
4.2 1.3 6.5
8.8 8.4 6.1
6.1 2.4 5.3
2.3 3.2 3.6
4.8 6.4 1.2

9. Membuat matrix dari 0 s/d 100 dengan penambahan 10

awk 'BEGIN {a=0;b=10; while(a <=100){print a;a=a+b}}' > numbers.txt

Dengan perintah seq
seq 0 10 100 > numbers.txt

more numbers.txt
0
10
20
30
40
50
60
70
80
90
100

10. Print input beberapa kali (contoh 3 kali):

more input.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3

more output.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3

#!/bin/sh
n=1
while [ $n -le 3 ]
do
awk '{print $0}' input.txt > > output.txt
n=$((n+1))
done

AWK Part-1: Manipulasi sebuah file ASCII

2011-09-30T19:13:00.000-07:00

Sebagai seorang Exploration Geoscientist, kita sering berhadapan dengan handling data yang cukup menantang.
Hal ini disebabkan karena ukuran data yang dihadapi sangat besar dengan jumlah baris dapat mencapai ribuan dan serta jumlah kolom mencapai jutaan.
Tentu saja data ini tidak bisa diatasi dengan software komersial seperti Microsoft Excell.

Pada bagian ini saya akan menjelaskan bagaimana memanipulasi sebuah file ASCII dengan AWK.

input.txt
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4

1. output yang diharapkan tersususun atas kolom 2 kolom 2 dan kolom 3
awk '{print $2,$2,$3} ' input.txt > output.txt

more output.txt
2 2 3
2 2 3
2 2 3
2 2 3

2. output yang diharapkan harus memiliki ketelitian. Kolom 1 adalah integer, kolom 2 dan 3 ketelitian 2 di belakang koma.
awk '{printf "%.1d %.2f %.2f\n ",$2,$2,$3} ' input.txt > output.txt

more output.txt
2 2.00 3.00
2 2.00 3.00
2 2.00 3.00
2 2.00 3.00

3. Operasi matematika antara kolom 2 dengan sebuah angka, dalam hal ini adalah perkalian (*)
awk '{printf "%.1d %.2f %.2f\n ",$2,$2*100,$3} ' input.txt > output.txt

more output.txt
2 200.00 3.00
2 200.00 3.00
2 200.00 3.00
2 200.00 3.00

4. Operasi matematika antar kolom. Dalam hal ini, setiap elemen di kolom 2 dibagi elemen kolom 3
awk '{printf "%.1d %.2f %.2f\n ",$2,$2/$3,$3} ' input.txt > output.txt

more output.txt
2 0.67 3.00
2 0.67 3.00
2 0.67 3.00
2 0.67 3.00

5. Output tersusun dari baris 3 sampai baris 5
more input.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3
4.8 6.4 1.2
6.1 2.4 5.3
8.8 8.4 6.1
8.1 0.6 8.7

awk '(NR>2 && NR<6) {print $0} ' input.txt > output.txt

more output.txt
5.1 4.2 9.3
4.8 6.4 1.2
6.1 2.4 5.3

6. Memilih baris ganjil
awk 'NR%2 ' input.txt > output.txt

more output.txt
2.3 3.2 3.6
5.1 4.2 9.3
6.1 2.4 5.3
8.1 0.6 8.7

7. Memilih baris genap
awk '(NR+1)%2 ' input.txt > output.txt

more output.txt
4.2 1.3 6.5
4.8 6.4 1.2
8.8 8.4 6.1

8. Membalikkan susunan baris:
awk '{a[i++]=$0} END {for (j=i-1; j>=0;) print a[j--] } ' input.txt > output.txt

more input.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3
4.8 6.4 1.2
6.1 2.4 5.3
8.8 8.4 6.1
8.1 0.6 8.7

more output.txt
8.1 0.6 8.7
8.8 8.4 6.1
6.1 2.4 5.3
4.8 6.4 1.2
5.1 4.2 9.3
4.2 1.3 6.5
2.3 3.2 3.6

9. Output diprint jika kolom 3 pada input lebih besar dari kolom 2
awk '$3 > $2 ' input.txt > output.txt

more output.txt
2.3 3.2 3.6
4.2 1.3 6.5
5.1 4.2 9.3
6.1 2.4 5.3
8.1 0.6 8.7

10. Menjumlahkan nilai-nilai pada setiap kolom
awk '{s=0; for (i=1; i<=NF; i++) s=s+$i; print s} ' input.txt > output.txt

11. Menjumlahkan semua nilai yang ada pada file
awk '{for (i=1; i<=NF; i++) s=s+$i}; END{print s} ' input.txt

12. Print sampai baris ke 2
awk 'NR < 3 ' input.txt

simpan hasil output.txt

awk 'NR < 3 ' input.txt > output.txt

13. Output baris ke 5
awk 'NR==5 ' input.txt
awk 'NR==5 {print;exit} ' input.txt # lebih efisien untuk file yang besar

14. Berikut adalah contoh meng-outputkan data yang memiliki header serta tidak memiliki TAB pemisah antar kolom. Untuk membaca data seperti ini, kita dapat menggunakan FIELDWIDTHS. Untuk melewati header, dapat digunakan NR.

more input.txt
1Ini adalah header
2Ini adalah header
3Ini adalah header
4.2x1.3Wx1.8wddwd
5.3x4.24x4.1hywkw
6.7x1.2Rx1.6wdwd
12345678901234567

more output.txt
4.2 1.3 1.8
5.3 4.2 4.1
6.7 1.2 1.6

awk 'BEGIN{FIELDWIDTHS="3 1 3 2 3"} NR>3 && NR<7 {print $1,$3,$5}' input.txt > output.txt

FIELDWIDTHS="3 1 3 2 3" menunjukkan:
kolom 1: 3 karakter pertama (karakter 1-3)
kolom 2: 1 karakter berikutnya(karakter ke 4)
kolom 3: 3 karakter berikutnya(karakter 5-7)
...dst
NR>3 && NR<7 menunjukkan hanya meng-outputkan baris ke 4 sampai 6
{print $1,$3,$5} menunjukkan hanya print kolom 1, 3 dan 5

Data Gathers Marine 2D: Radon Demultiple

2011-08-22T06:05:00.001-07:00

Sebagaimana yang telah saya jelaskan di sini, data ini memiliki multiple yang cukup dominan.

Salah satunya adalah water bottom multiple yang terlihat pada semblance yang rendah (~1500m/s), perulangan event dan frekuensi yang tinggi pada data gathers maupun brute stack. Salah satu cara untuk mengeliminasi multiple adalah Radon Demultiple.

Berikut adalah cara melakukan Radon Demultiple pada data marine_cmp_nmo.su (Ingat Radon Demultiple diterapkan setelah koreksi NMO)

1. Dengan menggunakan text editor (vi, gedit, kedit, dll.), copy kode Bourne Shell berikut lalu berinama radon. Jalankan dengan perintah sh radon

#!/bin/sh
##set -x
cdpmin=171
cdpmax=2142
filein=marine_cmp_nmo.su
n=$cdpmin
while [ $n -le $cdpmax ]
do
suwind <$filein key=cdp min=$n max=$n | suradon choose=3 igopt=4 pmin=-600 pmax=150 nmax=400 offref=-6000 > marine_radon_cmp$n.su
n=$((n+1))
done

Kode tersebut memecah data marine ke masing-masing cdp lalu menjalankan perintah radon. Saya menggunakan choose=3 untuk menghitung forward dan reverse tranform sekaligus dan igopt=4 (abs linear), pmin dan pmax mengontrol kemiringan event yang akan dieliminasi. pmin mengontrol event yang miring ke atas dan pmax event yang miring ke bawah.

Karena multiple adalah event yang miring ke bawah maka saya melakukan pembatasan yang cukup agresif yakni 150. nmax=400 adalah jumlah trace maksimum per cdp dan offref=-6000 adalah offset reference yang jauh lebih besar dari offset data tersebut. Untuk mempelajari lebih jauh anda bisa melakukan test untuk merubah parameter tersebut dan mengetikkan suradon pada terminal linux.

Gambar di bawah ini adalah data input dan hasil transformasi radon dengan option di atas. Perhatikan bahwa setiap sample data dengan amplitudo nol akan menghasilkan artifact pada Radon. Untuk menghilangkan artifact tersebut saya harus melakukan mute.

Data Input

Hasil Radon Demultiple

2. Ekstrak mute dari NMO data yang memiliki stretch mute.

Kode berikut memecah data menjadi masing-masing cdp.

#!/bin/sh
##set -x
cdpmin=171
cdpmax=2142
filein=marine_cmp_nmo.su
n=$cdpmin
while [ $n -le $cdpmax ]
do
suwind <$filein key=cdp min=$n max=$n > marine_nmo_cmp$n.su
n=$((n+1))
done

3. Ekstrak mute

Perintah berikut menghasilkan dummy gathers dimana sample data dengan nilai amplitudo nol akan nol dan yang lainnya akan bernilai 1.

#!/bin/sh
##set -x
cdpmin=171
cdpmax=2142
n=$cdpmin
while [ $n -le $cdpmax ]
do
suop2 marine_nmo_cmp$n.su marine_nmo_cmp$n.su op=quo > marinecmp_nmo_mute_cmp$n.su
n=$((n+1))
done

Contoh dummy gather dari perintah di atas.

4. Menerapkan mute:
Intinya meng-kalikan hasil Radon dengan dummy gather.

#!/bin/sh
##set -x
cdpmin=171
cdpmax=2142
filein=angolacmp_nmo.su
n=$cdpmin
while [ $n -le $cdpmax ]
do
suop2 marine_radon_cmp$n.su marinecmp_nmo_mute_cmp$n.su op=prod > marine_radon_muted_cmp$n.su
n=$((n+1))
done

Hasil Radon yang telah dimute.

Perbandingannya:

5. Menggabungkan semua cdp radon

Copy kode berikut lalu beri nama gabung dan jalankan dengan sh gabung

#!/bin/sh

cdpmin=171
cdpmax=2142

echo "#!/bin/sh" >> gabungradon
echo "##set -x" >> gabungradon
echo "cat \\" >> gabungradon
n=$cdpmin
while [ $n -le $((cdpmax)) ]
do
echo "marine_radon_muted_cmp$n.su \\" >>gabungradon
n=$(( n+1 ))
done
echo " > marine_radon_all_muted.su" >>gabungradon

6. Jalankan sh gabungradon dan stack.
sustack < marine_radon_all_muted.su > stack_radon.su
suximage < stack_radon.su perc=90 d1=0.004 f1=0

Data Gathers Marine 2D

2011-08-13T21:56:00.000-07:00

Data gathers marine 2D berikut ini telah dipresentasikan pada Geophysics Meeting dan sudah dibuka untuk publik. Saya memperoleh copy data tersebut dari sebuah training dari Third Party. Tentu saja penggunaannya hanyalah untuk kepentingan pendidikan bukan untuk kepentingan komersial.

Berikut adalah cara memperoleh data tersebut serta beberapa teknik processing dengan Seismic Unix dan software gratis lainnya:

1. Download data (format su) melalui link-link berikut:
marine1
marine2
marine3
marine4
marine5
marine6
marine7
marine8
marine9

2. Gabungkan semua data dengan menjalankan kode Bourne Shell berikut:

#!/bin/sh
##set -x
cat marine1.su marine2.su marine3.su marine4.su marine5.su marine6.su marine7.su marine8.su marine9.su > marine.su

Untuk mengkonversi su ke sgy:
sushw < marine.su key=d2 a=0.000 > marine_rest.su
segyhdrs ns=??? dt=??? < marine_rest.su
segywrite tape=marine.sgy endian=0 verbose=1 < marine_rest.su
Untuk mengetahui ns dan dt, lakukan surange < marine.su

3. Lakukan surange untuk melihat informasi header
surange < marine.su

Jika anda mendapatkan pesan error:

surange: fgettr.c: on trace #2 number of samples in header (0) differs from number for first trace (55393)

Maka anda harus melakukan perintah berikut terlebih dahulu:

suswapbytes < marine.su format=0 ns=1250 > marine_ok.su

mv marine_ok.su marine.su

Error di atas terjadi karena kompatibilitas dari satu komputer ke komputer lainnya.

4. Sorting menjadi cmp domain:
susort < marine.su > marine_cmp.su cdp offset

5. Lakukan pemilihan untuk setiap 200 cmp (persiapan analisa kecepatan)
suwind < marine_cmp.su key=cdp min=400 max=1800 j=200 > marine_every_200cmp.su

6. Lakukan analisa kecepatan dengan Velocity Picker

File > Open Trace File atau dengan ctrl+o lalu pilihlah file marine_every_200cmp.su

Berikut adalah tampilan Velocity Picker untuk analisa kecepatan data ini:

Setelah selesai melakukan picking,

Save Velocity Field (Su Format) ctrl+s maka sebuah file .par akan terbentuk.

7. Lakukan koreksi NMO

sunmo < marine_cmp.su > marine_cmp_nmo.su par=marine_every_200cmp.par

8. Lakukan Stack dan tampilkan

sustack < marine_cmp_nmo.su > marine_stack.su key=cdp
suwind < marine_stack.su key=tracl min=1 max=1972 | suximage perc=80 &

Berikut adalah stack data tersebut:

Dari gambar di atas terlihat bahwa data ini masih didominasi oleh multiple serta difraksi. Anda dapat mengatasinya dengan menerapkan Radon Demultiple dan Migrasi baik PSTM, PreSTM dan bahkan PreSDM. Perlu juga diterapkan proses-proses lain seperti source designature, rotasi wavelet dari minimum phase ke zero phase,median filtering, dsb.

Jika anda telah berhasil dan selesai melakukan processing (khususnya dengan Seismic Unix) anda dapat men-share hasil, gambar serta pengalamannya pada milist ensiklopediseismik.

Advanced Seismic Unix Processing Part-6

2011-01-02T22:47:00.000-08:00

Super CMP Gather

Super CMP Gather yang prinsipnya merupakan penjumlahan dari CMP yang berdekatan kerap sekali digunakan dalam analisa kecepatan. Dibandingkan dengan Single CMP gather, Super CMP gather memiliki S/N (Signal to Noise Ratio) yang lebih tinggi sehingga memberikan resolusi semblance yang lebih baik. Hal ini mudah dipahami karena sifat signal yang koheren dibandingkan dengan noise yang berubah secara waktu dan tempat. Penjumlahan akan mengguatkan yang koheren dan meng-cancel out yang tidak koheren.

Untuk mendemonstrasikan pembuatan Super CMP Gather dengan Seismic Unix, saya akan mengunakan CMP gather Marmousi yang terkontaminasi oleh noise.

Berikut adalah beberapa CMP gathers model Marmousi sebelum terkontaminasi noise:

suwind < data_cmp.su key=cdp min=100 max=103 | suxwigb perc=90 &
Lakukan suaddnoise untuk mengkontaminasi model Marmuosi dengan noise:
suaddnoise < data_cmp.su > data_cmp_noise.su sn=10
suwind < data_cmp_noise.su key=cdp min=100 max=103 | suxwigb perc=90 &

Jalankan Kode Bourne Shell berikut untuk menjumlahkan 2 CMP yang berdekatan:

#!/bin/sh
##set -x
cdpmin=1
cdpmax=574
filein=data_cmp_noise.su

# Memecah data menjadi beberapa CMP gathers:
n=$cdpmin
while [ $n -le $cdpmax ]
do
suwind < $filein key=cdp min=$n max=$n > cedepe_1$n.su
n=$(( n+1 ))
done
# Menjumlahkan 2 CMP yang berdekatan dengan suop2 op=sum:
n=$cdpmin
while [ $n -le $((cdpmax-1)) ]
do
suop2 cedepe_1$n.su cedepe_1$((n+1)).su op=sum > cedepe2_$n.su
n=$(( n+1 ))
done

# Untuk menjumlahkan 4 CMP yang berdekatan, lakukan proses
# di atas untuk Super CMP gather yang baru terbentuk.
# Banyaknya CMP yang dijumlahkan tentunya harus mempertimbangkan
# karakter lapisan dan geologi. Untuk lapisan yang datar menjumlahkan 6 CMP
# masih bisa dilakukan. Sedangkan untuk lapisan miring serta perubahan geologi
# yang tiba-tiba menjumlahkan 2 atau 4 CMP mungkin sudah cukup.
# Saya sarankan untuk selalu membandingkannya dengan Single Gathernya.

rm *cedepe_1*
fileout=data_cmp_noise_supergather.su

# Membuat Bourne Shell yang bernama gabung untuk menggabungkan Super gather:
echo "#!/bin/sh" >> gabung
echo "##set -x" >> gabung
echo "cat \\" >> gabung
n=$cdpmin
while [ $n -le $((cdpmax-2)) ]
do
echo "cedepe2_$n.su \\" >> gabung
n=$(( n+1 ))
done
echo "| sushw key=tracl a=1 > $fileout" >> gabung
echo "rm *cedepe*" >> gabung

Setelah kode di atas dijalankan, run shell gabung, dengan sh gabung.

Tampilkan beberapa Single CMP gather:
suwind < data_cmp_noise.su key=cdp min=100 max=103 | suxwigb perc=90 title="Single Gather" &

Tampilkan beberapa Super CMP gather:
suwind < data_cmp_noise_supergather.su key=cdp min=100 max=103 | suxwigb perc=90 title="Super Gather (sum 2 adjacent CMPs)" &

Perbandingannya:

Berikut adalah semblance dari Single Gather:
suwind < data_cmp_noise.su key=cdp min=100 max=102 | suvelan | suximage cmap=hsv1 legend=1 title="Single Gather" &

Semblance Super Gather:
suwind < data_cmp_noise_supergather.su key=cdp min=100 max=102 | suvelan | suximage cmap=hsv1 legend=1 title="Super Gather (sum 2 adjacent CMPs)"&
Perbandingannya:

Advanced Seismic Unix Processing Part-5

2011-01-01T05:30:00.000-08:00

Continuous Wavelet Transform

Continuous Wavelet Transform (CWT) sangat bermanfaat untuk menganalisa lapisan tipis, analis anomali 'sub-seismic', ataupun untuk mempelajari karakteristik kandungan frekuensi dari seismik.

Pada Seismic Unix, CWT dapat dilakukan dengan menggunakan fungsi sugabor. Input untuk sugabor adalah trace-trace seismik pada domain waktu, sedangkan outputnya adalah data seismik pada domain waktu-frekuensi yang dihasilkan melalui transformasi Gabor.

Pada bagian ini, saya akan menunjukkan CWT pada data seismik yang bisa diperoleh di sini (Line A).

Lakukan konversi dari format sgy menjadi su:
segyread tape=lineA.sgy conv=0 endian=0 > lineA.su

Tampilkan:
suximage < lineA.su
Lakukan surange untuk mempelajari headers dari line A:
surange < lineA.su

Lakukan suwind, untuk mengetest sugabor:
suwind < lineA.su key=cdp min=100 max=100 > lineA_cdp100.su

Lakukan sugabor:
sugabor < lineA_cdp100.su > lineA_cdp100_cwt.su
suximage < lineA_cdp100_cwt.su
Run kode Bourne Shell berikut untuk melakukakn CWT pada keseluruhan line A:

#!bin/sh
n3=390 #no of seismic traces
n1=1001 #no of time samples
d1=0.004 #time sampling interval
d3=1 #trace increment
n2=20 #no of frequency traces
ntr=7800 #no of seismic traces * freq traces
d2=6.25 #frequency sample interval

filein=lineA.su

sugabor < $filein | suxmovie n1=$n1 n2=$n2 ntr=$ntr n3=$n3 d1=$d1 d2=$d2 d3=$d3 title="CWT" label1="Time (sec)"
label2="Freq (hz)" loop=1 bclip=4 wclip=0 title="trace no %g"
loop=1 sleep=20000 -geometry 300x500 exit 0 &

Berikut adalah hasil CWT untuk beberapa trace dari line A:

Dari hasil di atas, saya melakukan time-variant bandpass filter dengan sutvband:
sutvband < lineA.su tf=.5,1.5,2.5,3.5 f=15,17.5,50,55 f=15,17.5,40,45 f=15,17.5,35,40 f=15,17.5,30,35 > lineA_tvf.su

Gambar di bawah ini menunjukkan penampang sebelum sutvband:
suwind < lineA.su tmax=1.5 | suximage cmap=rgb1 perc=90 title="Before TVF" &

Sesudah sutvband:
suwind < lineA_tvf.su tmax=1.5 | suximage cmap=rgb1 perc=90 title="After TVF" &

Perbandingannya:

Advanced Seismic Unix Processing Part-4

2010-12-26T10:31:00.000-08:00

Waveform Modelling

Seismic Unix memiliki fungsi sufdmod2 yang berguna untuk melakukan Waveform Modelling.

sufdmod2 memerlukan model kecepatan dalam domain kedalaman dengan format "C-style binary floating point numbers".

Berikut ini adalah script untuk melakukan Waveform Modelling untuk model kecepatan Marmousi.

Oleh karena sufdmod2 memerlukan memory yang cukup excessive, maka langkah pertama yang saya lakukan adalah mereduksi jumlah sampel (suresamp) velocity.su menjadi 25% (rf=0.25). Anda dapat menggunakan surange untuk mengevaluasi jumlah sampel sebelum dan setelah suresamp:

suresamp < velocity.su > velocity1.su rf=0.25

Selanjutnya melakukan suwind dengan j=4 untuk mereduksi jumlah trace menjadi 25%:
suwind < velocity1.su key=tracl j=4 > velocity2.su
mv velocity2.su velocity.su

Mengubah format dari su ke binary:
sustrip < velocity.su > velocity.h@ ftn=0

Konversi dari time ke depth domain:
suttoz < velocity.su > velocitydepth.su vfile=velocity.h@

Membatasi depth dengan maksimum depth=2400m:
suwind < velocitydepth.su > velocitydepth1.su tmax=2400
mv velocitydepth1.su velocitydepth.su

Mengubah format dari su ke binary untuk model kecepatan dalam domain kedalaman:
sustrip < velocitydepth.su > velocitydepth.h@ ftn=0

Ketika anda melakukan sustrip anda akan memperoleh informasi parameter berikut:
n1=348 n2=575 nt=348 ntr=575 ns=348

Gunakan text editor untuk membuat kode Waveform Modelling:

#! /bin/sh
##gunakan surange velocitydepth.su untuk memperoleh informasi d1 dan d2
n1=348
n2=575
d1=6.912000
nt=348
ntr=575
d2=9.216
ns=348
velfile=velocitydepth.h@

ximage < $velfile wbox=400 hbox=200 n1=$n1 n2=$n2 d1=$d1 d2=$d2 \
wbox=800 hbox=400 legend=1 title="Marmousi model" label1="Depth (m)" cmap=hsv2 \
label2="distance m " units="m/s" &

# input sufdmod2
xs=2000 zs=10 hsz=10 vsx=1000 verbose=2
vsfile="vsfile.su" ssfile="ssfile.su" hsfile="hsfile.su"
tmax=2.5 mt=18
label1="Depth m"
label2="Distance m"

sufdmod2 < $velfile nz=$n1 dz=$d1 nx=$n2 dx=$d2 verbose=1 \
xs=$xs zs=$zs hsz=$hsz vsx=$vsx dt=0.002 hsfile=$hsfile fmax=70 mono=0 \
vsfile=$vsfile ssfile=$ssfile verbose=$verbose \
tmax=$tmax abs=1,1,1,1 mt=$mt |
suxmovie clip=1.0 fframe=70 dframe=10 title=" %g ms" \
label1=$label1 label2=$label2 \
n1=$n1 d1=$d1 f1=$f1 n2=$n2 d2=$d2 f2=$f2 \
loop=1 sleep=20000 -geometry 800x400 &
exit 0

Gambar di bawah ini menunjukkan model Marmousi yang telah di-downsampling dengan parameter di atas:

Gambar di bawah ini menunjukkan snapshot untuk wavefield pada Model Marmousi dengan posisi source pada X=1000m dan depth=10m.

Gambar di bawah ini menunjukkan animasi wavefield pada Model Marmousi:

Advanced Seismic Unix Processing Part-3

2010-12-18T22:30:00.000-08:00

PreStack Time Migration

Pada bagian ini saya akan menunjukkan bagaimana caranya melakukan PreStack Time Migration (PreSTM) untuk data Marmousi dengan Seismic Unix. Seismic Unix menyediakan fungsi suktmig2d untuk melakukan PreSTM dengan metoda Kirchhoff. Data input untuk PreSTM adalah RMS Velocity atau Average velocity atau Stacking Velocity serta data Seismic pada domain Common Offset.

Marilah kita gunakan data interval velocity dan data seismic yang sebelumnya saya gunakan di sini.

Oleh karena data kecepatan yang kita miliki adalah interval velocity, maka langkah pertama adalah melakukan konversi dari interval velocity menjadi Average Velocity.

Gunakan perintah berikut untuk mengkonversi format su menjadi ascii:
suascii < velocity.su bare=1 > velocity.ascii

Selanjutnya buatlah kode octave dengan gedit dan berilah nama rmsvelgen.m:

clear; clc
nocmp=574; %%no cmp veocity harus sama dengan no cmp seismic
nosamp=751; %%no samples
nocmporig=2301; %%no cmp origial pada model kecepatan (kita harus melakuka desimasi)
data1=load('velocity.ascii');
data2=reshape(data1,nosamp,[]);
ind=[4:4:nocmporig-4]; %% desimasi cmp
data3=data2(:,ind);

for i=1:nocmp
trace{i}=data3(:,i);
end

for j=1:nocmp
trace{j}=data(:,j);
trace{j}(1)=trace{j}(1);
for i=1:nosamp-1
trace{j}(i+1)=trace{j}(i)+trace{j}(i+1);
end
end

for j=1:nocmp
for i=1:nosamp-1
traceb{j}(i)=trace{j}(i)/i;
end
end

vel2=cell2mat(traceb');
vel3=vel2';
save -ascii velrms.txt vel3

Run kode octave di atas dengan mengetik rmsvelgen pada terminal octave.

Pada terminal linux, ketiklah perintah berikut untuk mengkonversi dari ascii ke binary:
a2b < velrms.txt > velrms.bin

Convert dari binary ke su, flip dan tampilkan lalu QC:
suaddhead ns=574 < velrms.bin | sushw key=dt a=4000 > velrms.su
suflip flip=0 < velrms.su > velrms1.su
mv velrms1.su velrms.su
suximage < velrms.su legend=1 &

Ubah menjadi binary kembali, lalu smoothing, batasi hingga t=2.75s (menghindari noise) dan QC dengan menampilkannya:

sustrip < velrms.su > velrms.h@ ftn=0
smooth2 < velrms.h@ n1=750 n2=574 r1=20 r2=20 > smoothvelrms
suaddhead ns=750 < smoothvelrms | sushw key=dt a=4000 > smoothvelrms.su
suwind < smoothvelrms.su > smoothvelrms1.su tmax=2.75
mv smoothvelrms1.su smoothvelrms.su
suximage < smoothvelrms.su legend=1 &

Setelah window ximage aktif, tekan huruf 'h' pada keyboard untuk mengubah skala warna.

Ubah kembali ke binary:
sustrip < smoothvelrms.su > smoothvelrms.h@ ftn=0

Untuk data seismik, lakukan sorting untuk membuat domain Common Offset dan membatasinya hingga t=2.75s

susort offset gx < data.su > data_co.su
suwind < data_co.su > data_co1.su tmax=2.75
mv data_co1.su data_co.su

QC Offset -800
suwind < data_co.su key=offset min=-800 max=-800 | suximage perc=90 title="offset -800m"

PreSTM dilakukan offset by offset, sehingga buatlah kode Bourne Shell untuk menjalankan PreSTM. Lakukan surange untuk mempelajari parameter-parameternya.

#!/bin/sh
##set -x
offsetmin=-2575
offsetmax=-200
offdist=25
filein=data_co.su
vfile=smoothvelrms.h@
#################################################
n=$offsetmin
while [ $n -le $offsetmax ]
do
suwind < $filein key=offset min=$n max=$n | suktmig2d vfile=$vfile dx=25 > data_pstm_co$n.su
n=$((n+$offdist))
done

Setelah kode Bourne Shell tersebut dijalankan, maka anda akan memperoleh hasil PreSTM untuk masing-masing offset.

Buatlah Kode Bourne Shell Berikut untuk melakukan penggabungan file Common Offset:

#!/bin/sh
##set -x
offsetmin=-2575
offsetmax=-200
offdist=25
fileout=data_pstm_co.su
echo "#!/bin/sh" >> gabung
echo "##set -x" >> gabung
echo "cat \\" >> gabung
n=$offsetmin
while [ $n -le $((offsetmax)) ]
do
echo "data_pstm_co$n.su \\" >> gabung
n=$((n+$offdist))
done
echo "| sushw key=tracl a=1 > $fileout" >> gabung

Setelah selesai jalankan sh gabung, sehingga diperoleh file data_pstm_co.su

Lakukan sorting data sehingga diperoleh cdp gather:
susort cdp offset < data_pstm_co.su > data_pstm_cdp.su

Lakukan QC sehingga diperoleh gambar berikut:
suwind < data_pstm_cdp.su key=cdp min=100 max=103 | suxwigb perc=90 &

Lakukan Mute:
sumute < data_pstm_cdp.su > data_pstm_cdp_mute.su xmute=-2575,-200 tmute=1.1,0.1 key=offset

Lakukan QC:
suwind < data_pstm_cdp_mute.su key=cdp min=100 max=103 | suxwigb perc=90 &

Lakukan Stacking untuk memperoleh hasil Stack PreSTM:
sustack < data_pstm_cdp_mute.su | suximage perc=90 &

Analisa kecepatan biasanya dilakukan pada interval CMP tertentu serta sampling waktu yang tidak beraturan. Untuk mengatasi hal ini, saya membuat kode Octave untuk melakukan interpolasi CMP serta regularisasi time sampling, sehingga hasilnya dapat digunakan untuk melakukan PreSTM.

clear; clc
nocmp=22; %%% katakanlah pada model Marmousi kita hanya melakukan analisa kecepatan untuk 22 cmp
cmpint=20; %%% dengan interval cmp=20
ts=0.004; %%%sampling waktu
tmin=0; %%%waktu t0
tmax=3-ts; %%%waktu maksimum pada data Marmousi awal
cmp1=100; %%%CMP pertama dari analisa kecepatan
cmp2=520; %%%CMP terakhir dari analisa kecepatan
cmp1x=1; %%%CMP awal pada seismik
cmp2x=574; %%%CMP akhir pada seismik
ns=750; %%%jumlah sample pada data seismik

%File dari hasil analisa kecepatan atau Velocity Picker,
%hapus 2 baris pertama, tnmo=, vnmo=, dan \ (spasi \).
%dengan demikian filenya bisa dibaca oleh octave

fid=fopen('vpick.txt','rt');
for i=1:(nocmp*2)
line=fgets(fid);
tv{i}=str2num(line);
end

for i=1:nocmp
time{i}=[tmin,tv{i*2-1},tmax];
end

for i=1:nocmp
a=size(tv{i*2})(1,1); b=size(tv{i*2})(1,2);
vel{i}=[tv{i*2}(a),tv{i*2},tv{i*2}(b)];
end

%interpolasi dan regularisasi
tnew=[tmin:ts:tmax];

for i=1:nocmp
vnew{i}=interp1(time{i},vel{i},tnew);
end

vall=cell2mat(vnew');
vel1=vall';

velearly=repmat(vel1(:,cmp1x),1,(cmp1-cmp1x));
vellast=repmat(vel1(:,nocmp),1,(cmp2x-cmp2));

for i=1:ns
valhornew{i}=interp1([cmp1:cmpint:cmp2]', vall(:,i),[cmp1:cmp2]');
end

velmid=cell2mat(valhornew);
vel2=velmid';

velocity=[velearly,vel2,vellast];
save -ascii velocity.txt velocity

Lakukan proses di atas untuk melakukan smoothing dan konversi ke binary.
Gambar berikut adalah RMS Velocity dengan interpolasi dari hasil analisa kecepatan.

Interactive Velocity Analysis SU

2010-12-12T19:59:00.000-08:00

Analisa kecepatan lebih nyaman dilakukan secara interaktif. Jika anda bekerja dengan Seismic Unix, anda dapat menggunakan Velocity Picker 1.1 yang bisa didownload secara gratis.

Berikut adalah cara untuk mendownload, menginstall dan mempergunakan Velocity Picker 1.1:

Download Velocity Picker 1.1 di sini , sehingga anda memperoleh file yang bernama: VelocityPicker V1.1.zip
Oleh karena file tersebut mengandung spasi, rename file tersebut menjadi VelocityPicker (contoh)
Pada terminal linux, ketik unzip VelocityPicker . Jika belum memiliki fungsi unzip, ketik sudo apt-get install unzip.
Setelah diunzip, semua file akan tersimpan dalam directory yang bernama dist
Masuklah ke dalam directory dist, lalu ketik java -jar "VelocityPicker_V2.jar"
Sebuah windows Velocity Picker akan muncul
File > Open Trace File (CMP Gathers dengan format SU)
Nomor CMP yang akan dianalisa, dapat dipilih pada scroll box
Pick pada puncak semblance dengan MB1
Setelah selesai picking, Save Velocity Field (SU Format)
File yang telah disave tersebut, selanjutnya anda copy dan paste di Bourne Shell nmo.sh
Untuk mempelajari tentang fitur aplikasi ini, lihatlah Help > Instructions

Berikut adalah snapshot dari aplikasi Velocity Picker 1.1. yang saya gunakan untuk analisa CMP gathers Marmousi:

Advanced Seismic Unix Processing Part-2

2010-12-11T06:09:00.000-08:00

PreStack Depth Migration

Seismic Unix (SU) memiliki kapabilitas untuk melakukan PreStack Depth Migration (PreSDM). Fungsi dalam SU untuk melakukan PreSDM adalah sukdmig2d, dimana metoda yang digunakannya adalah metoda Kirchhoff.

Untuk mendemonstrasikan bagaimana caranya melakukan PreSDM, saya akan menggunakan shot gathers serta model kecepatan Marmousi.

Data Marmousi bisa didownload di Gzipped tar file (46 MB)
Ekstrak data yang telah didownload dengan perintah tar -xvvf ieee.tar.gz
Semua file yang telah diekstrak akan berada di folder ieee, file-file tersebut adalah: data.segy (seismic shot gathers), velocity.segy dan density.segy
Lakukan konversi dari format segy ke su dengan perintah:
segyread tape=data.segy conv=0 endian=0 > data.su
segyread tape=velocity.segy conv=0 endian=0 > velocity.su

Coba tampilkan model kecepatan dengan perintah suximage < velocity.su maka anda akan memperoleh model kecepatan yang terotasi 90 derajat.
Lakukan rotasi -90 derajat untuk model kecepatan tersebut dengan perintah:
suflip flip=0 < velocity.su > velocity1.su
mv velocity1.su velocity.su
Sekarang coba tampilkan kembali, lalu tekan huruf 'r' untuk mengubah skala warna. sehingga anda akan memperoleh gambar sbb (Model Marmousi dibuat berdasarkan tipikal model geologi di West Africa-Angola):

Sekarang lihatlah data seismik, lakukan sorting dari shot gather ke CMP gather dengan teknik yang saya jelaskan pada bagian ini, lakukan analisa kecepatan, koreksi NMO dan stack. Sehingga anda memperoleh gambar sbb:

Gambar di atas merupakan stack dari Model Marmousi dengan menggunakan stacking velocity. Terlihat dominasi difraksi akibat zona yang kompleks pada Model Marmousi, demikian juga dengan pencitraan yang kurang baik pada bagian yang lebih dalam. Pada tahapan berikutnya kita akan berupaya untuk memperbaiki citra tersebut dengan menggunakan PreSDM.
PreSDM memiliki tiga tahapan utama: smoothing (smooth2) model kecepatan, raytracing (rayt2d) dan depth migration (sukdmig2d). Dengan data input: shot gathers (domain time) dengan format su, velocity dengan format "C-style binary floating point numbers". Data output adalah: Stack dalam domain depth.
Ubahlah format velocity dari su ke "C-style binary floating point numbers" dengan perintah: sustrip < velocity.su > velocity.h@ ftn=0
Mengingat data seismik memiliki noise di atas 2.9s, maka kita harus menghilangkannya dengan perintah:

suwind < data.su > data1.su tmax=2.9
mv data1.su data.su

Dengan gedit, copy-lah kode Bourne Shell berikut lalu save dengan nama presdm.sh:

#!/bin/sh
# Parameter Model Kecepatan
nz=751
nx=2301
dz=4
dx=4
# Parameter Seismik
nt=750
ntr=96
dt=4000
ifile=data.su #data input seismik shot gathers
ofile=datamig.su # data output stack depth migrated
tfile=tfile # traveltimes
vfile=velocity.h@ # input kecepatan dengan format C floating point

# smoothing model kecepatan
smooth2 < $vfile n1=$nz n2=$nx r1=20 r2=20 >smoothvel

#raytracing untuk meghasilkan tfile
rayt2d <smoothvel dt=0.004 nt=751 dz=$dz nz=$nz dx=$dx nx=$nx \
fxo=0 dxo=25 nxo=369 fxs=0 dxs=100 nxs=93 >$tfile

#migrasi
sukdmig2d infile=$ifile datain=$ifile outfile=$ofile dataout=$ofile \
ttfile=$tfile fzt=0 dzt=4 nzt=751 fxt=0 dxt=25 nxt=369 \
fs=0 ns=93 ds=100 nzo=751 dzo=4 dxm=25 mtr=1

Pada terminal linux ketiklah sh presdm.sh
Setelah selesai (Intel Atom, 1GB MEM: kurang lebih 60 menit), tampilkan hasilnya dengan suximage < datamig.su, sehingga diperoleh hasil sbb:

Bandingkan hasilnya dengan stack di atas, menakjubkan?

Advanced Seismic Unix Processing Part-1

2010-12-06T16:38:00.000-08:00

Radon Demultiple

Pada bagian ini, saya akan menunjukkan bagaimana caranya melakukan eliminasi multiple melalui Transformasi Radon.

Eliminasi multiple dengan teknik Radon atau dikenal juga dengan Radon Demultiple memiliki 4 tahapan utama:

Koreksi NMO (CMP gather)
Transformasi ke domain τ-slope index
Desain Mute
Inverse τ-slope index

Setelah koreksi NMO, karakteristik multiple akan berbeda dengan reflektor. Jika nilai kecepatan NMO yang dipilih tepat maka reflektor akan menjadi datar (flat) sedangkan multiple masih melengkung ke arah bawah (karena nilai kecepatannya lebih rendah). Jika ditransformasi ke dalam domain τ-slope index (baca: sumbu y adalah waktu dan x adalah kemiringan event), maka reflektor yang datar akan memiliki slope index bernilai sekitar nol dan multiple akan bernilai positif. Sedangkan untuk event yang melengkung ke atas (over NMO-ed, difraksi, side swipe) akan memiliki nilai slope index negatif.

Gambar di atas (kiri) menunjukkan CMP sintetik dengan reflektor pada waktu 0.2s dan reflektor yang bersamaan dengan multiple pada 0.4s, sedangkan gambar sebelah kanan adalah CMP dalam domain τ-slope index. Perhatikan reflektor yang flat pada waktu 0.2s dan 0.4s memiliki nilai slope index di sekitar nol, sedangkan multiple pada 0.4s memiliki nilai slope index sekitar 100.

Berikut adalah tahapan untuk menghasilkan gambar di atas:

Run kode Bourne Shell berikut:

#!/bin/sh
npl=3 # number of planes
nt=200 # number of time samples
ntr=100 # number of traces
taper=0 # no end-of-plane taper = 1 taper planes to zero at the end
offset=400 # offset
dt=0.004 # time sample interval in seconds
suplane npl=$npl nt=$nt ntr=$ntr taper=$taper offset=$offset \
dip1=0 len1=128 ct1=100 cx1=50 \
dip2=1 len2=128 ct2=115 cx2=50 \
dip3=0 len3=128 ct3=50 cx3=50 > planemul.su

2. Tampilkan:

sushw < planemul.su key=offset a=0 b=10 | sufilter f=0,5,50,60 > planemul_fil.su
suximage < planemul_fil.su title="Sebelum Radon Demultiple" key=offset label1="time (s)" label2="offset (meters)" wbox=500 hbox=500 &

3. Lakukan Tranformasi Radon dan tampilkan:

suradon < planemul_fil.su igopt=3 interp=4 choose=0 depthref=1000 interoff=0 offref=1190 pmin=-300 pmax=300 > radon.su
suximage < radon.su perc=99 label1="tau (s)" label2="slope index" title="suplane data Radon transformed" wbox=500 hbox=500 &

Gambar di atas adalah CMP dalam domain τ-slope index, sebelum (kiri) dan setelah mute.

Untuk melakukan Mute:

sumute < radon.su key=f2 \
xmute=-300,30,31,300 \
tmute=0.8,0.79,0.01,0 \
mode=1 ntaper=10 > radon_muted.su

suximage < radon_muted.su label1="tau (s)" label2="slope index" title="suplane data Radon transformed-muted" wbox=500 hbox=500 perc=99 &

Berikut adalah perintah untuk melakukan inversi τ-slope index setelah dimute serta menampilkan hasilnya:

suradon < radon_muted.su igopt=3 interp=4 choose=4 depthref=1000 interoff=0 offref=1190 pmin=-300 pmax=300 | suximage title="Setelah Radon Demultiple" key=offset label1="time (s)" label2="offset (meters)" wbox=500 hbox=500 &

Dari gambar di atas terlihat bahwa, multiple pada 0.4s telah dieliminasi dengan sempurna.

Gambar di bawah ini menunjukkan muting untuk data real pada domain τ-slope index . Kiri (sebelum mute) kanan (setelah mute).

Pemilihan batas pengambilan mute, harus dilakukan dengan hati-hati, jangan sampai melukai sinyal reflektor, mengingat berbagai hal seperti error pemilihan nilai kecepatan NMO, anisotropic, dll. Jadi, sebaiknya kita melakukan test dengan mengembalikannya ke domain waktu-offset, lalu lihatlah perilaku reflektornya. Gambar di bawah ini menunjukkan perbandingan CMP gather sebelum (kiri) dan setelah Radon Demultiple(kanan).

Seismic Processing with Seismic Unix - Part9

2010-12-01T20:37:00.000-08:00

Tahap sebelumnya

Tahap 18:

Pada tahapan ini, saya akan menunjukkan bagaimana melakukan PoSTM (Post Stack Time Migration) dengan menggunakan Seismic Unix.

Seismic Unix menyediakan fungsi migrasi dengan menggunakan beberapa metodologi diantaranya Stolt Migration, Gazdag atau Phase-Shift Migration, Claerbout's Migration, dll. Teknik-teknik migrasi tersebut memiliki kelebihan dan kekurangannya masing-masing.

Saya akan menerapkan PoSTM pada data telah distack dengan menggunakan velocity analysis dan juga telah mengalami koreksi statik. Saya memilih Metoda Stolt, karena metoda ini sangat cepat dan cukup robust.

Sebelumnya, saya akan menerapkan mute dan tapering (walaupun pada sustolt sendiri diterapkan tapering) sehingga diperoleh ujung lintasan kiri, kanan dan atas yang lebih gradual. Hal ini penting dilakukan karena perbedaan amplitudo yang tiba-tiba akan menghasilkan migration artifacts.

Setelah itu migrasi diterapkan dengan menggunakan kecepatan RMS (yang saya ambil dari CMP ke 1000) dari hasil velocity analysis sebelumnya.

Berikut ini adalah kode Bourne Shell untuk melakukan mute, menampilkan stack sebelum migrasi, melakukan Stolt Migration dan sekaligus menampilkan hasilnya:

#!/bin/sh

sumute < Line001_stack_vel2_elev_rstat.su key=tracl xmute=1,150,1132,1132 tmute=3.0,0.2,0.2,3.0 ntaper=50 > Line001_stack_vel2_elev_rstat_mute.su

suximage < Line001_stack_vel2_elev_rstat_mute.su key=cdp perc=90 title='Before Stolt Migration' &

time=0.0187891,0.494781,0.914405,1.37787,1.94781,2.90605 vels=1992.35,2211.92,2488.77,2765.61,2975.64,3319.31

sustolt < Line001_stack_vel2_elev_rstat_mute.su \
cdpmin=1 cdpmax=1282 dxcdp=50 \
tmig=$time vmig=$vels \
smig=0.6 vscale=1 lstaper=50 lbtaper=50 | suximage title="After Stolt Migration" key=cdp min=100 max=1100 perc=90 verbose=0 &

exit

Copy-lah code di atas dengan menggunakan text editor (vi, pico, gedit, atau nedit) lalu save dengan nama tertentu katakanlah stoltmig. Lalu ketik sh stoltmig

Berikut adalah stack sebelum migrasi (untuk mengubah skala warna, tekan hurup r pada ximage yang aktif, untuk mengeksplorasi lebih jauh ketik ximage pada terminal linux).

Setelah Migrasi:

Animasi (klik untuk memperbesar)

Seismic Processing with Seismic Unix - Part8

2010-11-28T04:34:00.000-08:00

Lihat tahap sebelumnya

Pada bagian ini, saya akan menunjukkan bagaimana caranya melakukan elevation statics dan residual statics untuk data seismik yang kita miliki.
Elevation statics umumnya dilakukan sebelum koreksi NMO pada tahap 13. Akan tetapi untuk melihat efek elevation statics terhadap citra seismik, saya lakukan setelah memperoleh citra yang terbaik. Hal ini sah-sah saja untuk dilakukan, mengingat kita masih memiliki peluang untuk terus memperbaiki citra tersebut diantaranya dengan analisa kecepatan pada interval CMP yang lebih rapat, analisa pada super gather, dll.

Informasi yang harus kita miliki untuk melakukan elevation statics adalah elevasi sumber-penerima relatif terhadap datum serta kecepatan sedimen di bawah sumber-penerima sehingga diperoleh waktu tempuh gelombang dari elevasi yang bersangkutan terhadap datum.

Untuk data ini, waktu tempuh sumber dan penerima (elevation statics) telah dilakukan dan sudah diselipkan ke dalam trace header. Sehingga, jika kita melakukan surange diperoleh hasil sbb:

Dari hasil surange di atas, terlihat elevation statics dari sumber (selev) dan penerima (gelev) telah berada pada trace header.

Jika kita tidak memiliki informasi gelev dan selev, maka kita harus menghitungnya melalui first break picking dari gelombang refraksi.

Elevation statics dilakukan dengan perintah:

sustatic < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom.su > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_elev_stat.su hdrs=1

sustatic adalah perintah untuk elevation static, hdrs=1 menunjukkan bahwa data statics berada pada trace header. Untuk mempelajari options atau parameter yang digunakan pada sustatic, ketiklah sustatic pada terminal linux.

Gambar dibawah ini adalah shot gather sebelum elevation static:

Gambar dibawah ini adalah shot gather setelah elevation statics yang ditampilkan dengan perintah:

suwind < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_elev_stat.su key=ep min=32 max=32 | suxwigb perc=80 & Perbandingannya (klik untuk memperbesar):

Gambar dibawah ini menunjukkan stack sebelum elevation statics yang kita peroleh sebelumnya:

Gambar di bawah ini adalah stack setelah elevation statics dengan dengan mengunakan model kecepatan dari analisa kecepatan sebelumnya. Lakukan proses yang sama (Tahap 16), dengan mengganti data input dan output.

Saya melakukan sorting dari shot ke cmp dengan perintah:

susort cdp offset < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_elev_stat.su > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev.su

Lalu, pada nmo.sh:
Input: Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev.su
Output: Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2.su

Run dengan mengetikkan sh nmo.sh

Lakukan stacking:
sustack < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2.su > Line001_stack_vel2_elev.su

Tampilkan:
suximage < Line001_stack_vel2_elev.su cmap=hsv17 perc=90 title='After Elevation Statics' & Berikut perbandingannya:

Tahap 17:

Pada Seismic Unix, residual static dilakukan dengan perintah suresstat dimana metoda yang digunakan mengacu pada Ronen dan Claerbout, Geophysics 50, 2759-2767 (1985).

Berbeda dengan elevation statics, residual statics dilakukan setelah koreksi NMO, akan tetapi harus dilakukan pada domain shot gather dengan key=fldr.

Berikut adalah tahapan untuk melakukan residual statics:

Lakukan sorting:
susort < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2.su > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2_fldr.su fldr offset

Hitung residual statics:
suresstat < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2_fldr.su ssol=sstats rsol=rstats ntraces=70782 ntpick=50 niter=5 nshot=481 nr=282 nc=70782 sfold=282 rfold=282 cfold=284

ntraces=70782 (jumlah seluruh trace pada data)
ntpick=50 (banyaknya sample maksimum untuk melakuan shifting dalam korelasi) niter=5 (jumlah iterasi, pada mesin 32bit dengan Intel Core Duo, 1.5GB Mem memerlukan waktu 12-15 menit)
nshot=481 (fldr maksimum, lakukan surange untuk melihat semua key)
nr=282 (jumlah receiver maksimum pada shot)
nc=70782 (harus sama dengan banyaknya seluruh trace) sfold=282 (harus sama dengan nr)
rfold=282 (maksimum ep)
cfold=284 (maksimum cdpt)

Perintah di atas akan menghasilkan dua file i.e. sstats and rstats yang masing-masing berisikan source and receiver statics. Terapkankanlah statics tersebut dengan perintah:

sustatic < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2.su > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2_rstat.su \
hdrs=3 sou_file=sstats rec_file=rstats ns=481 nr=1282

hdrs=3 (Angka 3 memerinta sustatic untuk membaca statics dari file)
sou_file=sstats (file source statics)
rec_file=rstats (file receiver statics)
ns=481 (harus sama dengan nshot pada command sebelumnya)
nr=1282 (banyaknya cdp didalam stack data)

Gambar di bawah ini adalah stack sebelum elevation dan residual statics.
Gambar di bawah ini adalah penampang setelah elevation dan residual statics.
Lakukan perintah-perintah sbb:

Stacking:
sustack < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_elev_nmo_v2_rstat.su > Line001_stack_vel2_elev_rstat.su

Dikarenakan kedua koreksi statics tersebut mengakibatkan time shifting, maka untuk kenyamanan dalam membandingkan hasilnya, saya lakukan shiting dengan perintah:

suresamp < Line001_stack_vel2_elev_rstat.su > Line001_stack_vel2_elev_rstat_shift.su ns=1501 dt=0.002 tmin=0.25

ns=1501 (banyaknya sample)
dt=0.002 (samping interval dalam detik)
tmin=0.25 (di shifting ke atas 0.25 detik)

Tampilkan:
suximage < Line001_stack_vel2_elev_rstat_shift.su cmap=hsv17 perc=90 title='After Elevation and Residual Statics' &
Dan perbandingannya:

Dari perbandingan di atas, terlihat bahwa residual statics masih belum memberikan hasil yang optimal, walaupun di beberapa tempat memberikan hasil yang lebih baik.

Anda masih bisa memperbaiki hasilnya dengan melakukan perubahan dari parameter-parameter yang saya pilih sebelumnya.

Tahap berikutnya

Tips untuk membuat animasi GIF dan menampilkannya pada blogspot:

Install imagemagick untuk mengkonversi png ke gif format, sudo apt-get install imagemagick (saya perlu ini karena hasil screen capture saya berupa png)
Pada folder dimana file-file png berada ketik: for file in *.png; do echo "convert $file $(basename $file .png).gif"; done
Copy semua print out dari tahap 2 lalu paste pada terminal linux (print outnya akan seperti: convert filename.png finename.gif)
Install gifsicle untuk menghasilkan animasi gif: sudo apt-get install gifsicle
Contoh perintah untuk menganimasi dua file gif: gifsicle --delay 200 --loop=1000 stack1.gif stack2.gif > anim.gif
Buat account photobucket http://photobucket.com/
Upload anim.gif
Pada photobucket, double clicks gambar anim.gif
Click HTML code (di sebelah kanan anim.gif, dalam 'Share this photo' box)
Masuk ke blogspot, new post in HTML Mode lalu paste html code dari tahap 9

Seismic Processing with Seismic Unix - Part7

2010-11-27T11:45:00.000-08:00

Lihat Tahap sebelumnya

Tahap 15:

Pada tahapan ini, saya akan menunjukkan bagaimana melakukan velocity analysis dengan Interactive Velocity Analysis.

Lihat juga contoh penggunaan tool tersebut disini.

Lakukan pemilihan setiap 50 CMP pada data input berikut:
Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp.su, yang kita miliki sebelumnya.

suwind < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp.su key=cdp min=200 max=1000 j=50 > inva_every_50cmp.su

Jika proses picking telah dilakukan untuk seluruh cmp yang dimiliki, maka secara otomatis akan terbentuk file inva_every_50cmp.par

Tahap 16:

Setelah diperoleh model kecepatan, maka kita siap untuk melakukan koreksi NMO untuk seluruh CMP.

sunmo < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp.su > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_nmo_v2.su par=inva_every_50cmp.par

Dengan output
Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_nmo_v2.su

Setelah itu lakukan proses stacking:
sustack < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_nmo_v2.su > Line001_stack_vel2.su

Lalu tampilkan dengan perintah:

suximage < Line001_stack_vel2.su cmap=hsv17 perc=90 title='Setelah Velocity Picking' & Gambar di atas adalah stack dengan menggunakan velocity analysis yang saya pick. Bandingkan hasilnya dengan brute stack yang dihasilkan pada tahapan sebelumnya (gambar di bawah):

Berikut adalah animasi perbandingan antara brute stack (kecepatan tunggal) dengan stack dari velocity pick (click untuk memperbesar).

Tahap berikutnya

Seismic Processing with Seismic Unix - Part6

2010-11-26T04:47:00.000-08:00

Lihat tahapan sebelumnya

Tahap 13:
Pada tahap ini saya akan menunjukkan bagaimana caranya melakukan koreksi NMO dan membuat Brute Stack.

Pada Seismic Unix, koreksi NMO dilakukan dengan perintah sunmo. Informasi penting untuk sunmo yang harus kita berikan adalah pasangan kecepatan(m/s) dan waktu(s).

Untuk kasus data ini, saya melakukan trial and error untuk mencari 3 pasangan waktu dan nilai velocity yang paling optimal (yang ditentukan dengan reflektor yang flat). Test ini saya lakukan pada CMP 1000 s/d 1002 (lihat tahapan sebelumnya untuk melihat CMP yang belum dikoreksi NMO). Pasangan kecepatan dan waktu yang diperoleh adalah vnmo=1700,2750,3000 dan tnmo=0.1,1.0,2.0

Berikut adalah kode untuk koreksi NMO (buat dalam satu baris):

sunmo vnmo=1700,2750,3000 tnmo=0.1,1.0,2.0 < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp.su > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_nmo_v1.su

Lalu saya menampilkan hasilnya:

suwind < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_nmo_v1.su key=cdp min=1000 max=1002 | suxwigb perc=95

Sehingga diperoleh gambar sbb (klik untuk memperbesar):

Pada gambar di atas terlihat bahwa reflektor pada CMP 1000-1002 sudah terlihat cukup flat. Disamping itu proses muting nampaknya telah dilakukan secara otomatis bersamaan dengan nmo. Muting tersebut dilakukan untuk event yang mengalami NMO strech yang signifikan.

Tahap 14:

Setelah kita cukup puas dengan hasil NMO, maka kita siap untuk memproduksi brute stack dengan perintah:

sustack < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp_nmo_v1.su > Line001_brute_stack.su

Lalu saya tampilkan dengan:

suximage < Line001_brute_stack.su cmap=hsv17 perc=90

Sehingga diperoleh penampang seismik dibawah ini (klik untuk memperbesar):

Untuk mengubah skala warna pada brute stack, kita bisa melakukannya dengan mengubah parameter cmap. Berikut adalah contoh jika saya menggunakan cmap=hsv4.

Tahap berikutnya

Seismic Processing with Seismic Unix - Part5

2010-11-25T19:57:00.000-08:00

Lihat tahapan sebelumnya

Tahap 11:

Pada tahapan ini, saya akan menunjukkan bagaimana caranya melakukan edit geometry. Edit geometry sangat penting dilakukan karena bermanfaat untuk sorting dari shot gather ke CMP gather, koreksi statik, regularisasi bin, dll.

Jika kita melakukan surange, berikut adalah key yang kita miliki:

Dari hasil surange di atas terlihat bahwa data tersebut belum memiliki informasi geometri seperti koordinat sumber, penerima, offset, cdp, dll.

Marilah kita lihat kembali 3 file geometri untuk data ini i.e. Line_001.SPS (sumber), Line_001.RPS (receiver), Line_001.XPS (relational).

Jika kita membuka file Line_001.SPS dengan gedit, maka akan diperoleh informasi sbb:

Di dalam file Line_001.SPS, baris ke 1 s/d 20 merupakan header yang menjelaskan lokasi kolom (COLS) untuk setiap parameter sumber: point number (nomor sumber), Static Correction, Map grid easting (X coordinate), Map grid northing (Y coordinate), Surface Elevation, dll.

Data ini memiliki point number (shot number) 701, 703, 705, ...1201. Informasi 1V1 bukanlah bagian dari point number (lihat kembali header, COLS dari point number adalah 18-25), dengan gedit informasi posisi baris (Ln) dan kolom (Col) bisa dilihat di pojok kanan bawah. Dari sini kita mengetahui bahwa jumlah sumber (shot) adalah (1201-701)/2+1 = 251 shots. Koreksi statik untuk shot pertama dan seterusnya: -50, -50, -51, dst. Koordinat X shot pertama: 688081.8, koordinat Y shot pertama: 3838302.1, dan elevasi shot pertama: 46.0.

Demikian juga dengan file Line_001.RPS (informasi receiver):

Point number (receiver number) adalah 561, 562, ..., 1342. Koreksi statik untuk receiver 561: -48, Koordinat X untuk receiver 561:684590.2, koordinat Y receiver 561:3837867.6, dan elevasi receiver 561:41.0.

Sedangkan file Line_001.XPS, memuat informasi hubungan sumber penerima:

Perhatikan informasi di atas dengan sebaik-baiknya:
Point Number (COLS 30-37): 701, 703, ...1201 merupakan penomoran untuk shot pertama, kedua, dst.

From receiver untuk shot 701: 561 (receiver pertama untuk shot 701)
To receiver untuk shot 701:842 (receiver terakhir untuk shot 701)
From Channel : 1 dan To Channel: 282 untuk semua shot adalah sama artinya setiap shot memiliki jumlah trace 282.

Setelah kita memahami konfigurasi file di atas, saya akan membuat sebuah matrix dengan jumlah kolom 10 i.e. [sx,xy,selev,sstat,gx,gy,gelev,gstat,cdp,offset] dengan jumlah baris sebanyak jumlah trace yang saya miliki (70782).

Dengan menggunakan gedit copy-lah kode berikut, save, lalu beri nama geom

#!/bin/sh
# skip header dengan (NR>20), remove 1V1
# ekstrak source number (col2), sx (col8),sy (col9), selev (col10),sstat (col3)
awk ' gsub(/1V1/,"") {if (NR > 20) {print $2,$8, $9, $10,$3 }}' Line_001.SPS > sps.txt

# skip header dengan (NR>20), remove 1G1
# ekstrak receiver number (col2), gx (col8),gy (col9), gelev (col10),gstat (col3)
awk ' gsub(/1G1/,"") {if (NR > 20) {print $2,$8, $9, $10,$3 }}' Line_001.RPS > rps.txt

Pada terminal linux ketik sh geom

Lalu dengan gedit copy-lah kode berikut, save, lalu beri nama geomoctave.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc
load sps.txt
load rps.txt

% Ingat no of shots: 251,no of receiver in each shot: 282,total number of traces is 251*282=70782

%%%menyusun matrix sps untuk seluruh trace (70782)%%%%%%
for i=1:251
sps_for_traces_in_each_shot{i}=repmat(sps(i,:),282,1);
end
sps_all_traces=cell2mat(sps_for_traces_in_each_shot');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%menyusun matrix rps untuk seluruh trace (70782)%%%%%%
for i=1:251
rps_for_traces_in_each_shot{i}=rps([(i*2)-1:281+(i*2)-1]',:);
end
rps_all_traces=cell2mat(rps_for_traces_in_each_shot');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% menghitung offset
sx=sps_all_traces(:,2);
sy=sps_all_traces(:,3);
selev=sps_all_traces(:,4);
sstat=sps_all_traces(:,5);

gx=rps_all_traces(:,2);
gy=rps_all_traces(:,3);
gelev=rps_all_traces(:,4);
gstat=rps_all_traces(:,5);
ox=gx-sx;
oy=gy-sy;
offset=sqrt(ox.^2+oy.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%menghitung cdp untuk masing-masing trace
%%%dikarenakan interval geophone 25m dan interval sumber 50m, maka:
%%%cdp untuk shot pertama adalah 1 s/d 282
%%%cdp untuk shot kedua adalah 5 s/d 286
%%%cdp untuk shot ketiga adalah 9 s/d 290 dst....(lihat ilustrasi stacking diagram dibawah untuk memahaminya)
for i=1:251
cdp_each_shot{i}=[(4*i-3):281+(4*i-3)]';
end
cdp_all_traces=cell2mat(cdp_each_shot');

geom_header=[sx,sy,selev,sstat,gx,gy,gelev,gstat,cdp_all_traces,offset];

save -ascii geom_header.txt geom_header

%%%%plot koordinat sumber dan penerima
plot(sx,sy,'r*'); hold on
plot(gx,gy,'b^')
%%%%akhir dari kode

Lalu ketik octave untuk memasuki terminal octave. Jika belum memiliki octave (octave adalah open source yang mirip dengan Matlab), anda bisa menginstall dengan mengetik sudo apt-get install octave3.2

Pada terminal octave ketik geomoctave untuk mengeksekusi kode di atas sehingga kita akan memiliki matriks geom_header.txt dan plot koordinat sumber dan penerima.

Keluarlah dari octave dengan menekan ctrl+z

Gambar dibawah ini adalah ilustrasi stacking chart untuk interval sumber 2X interval receiver. Perhatikan CDP yang pertama untuk shot kedua terletak pada CDP ke 5 dari shot pertama, CDP yang pertama untuk shot ketiga terletak pada CDP ke 9 untuk shot kedua, dst.

courtesy freeusp.org

Gambar di bawah ini adalah koordinat sumber (merah) dan penerima (biru). Perhatikan bahwa posisi sumber dan penerima tidak benar-benar berada dalam satu garis lurus, hal demikian terjadi karena kondisi medan yang bersangkutan. Jadi, walaupun secara teoritik kita menghendaki group interval 25m dan shot interval 50m, tetapi pada kenyataannya sangat sulit untuk diwujudkan.

Selanjutnya, ubahlah format dari geom_header.txt menjadi binary dengan mengetikkan:

a2b < geom_header.txt n1=10 > geom_header.bin

a2b adalah perintah untuk mengubah format ascii ke binary, n1=10 adalah jumlah kolom dalam matriks geom_header.txt

Check-lah dengan mengetikkan ls -l, lihat jika ukuran file geom_header.bin adalah 2831280 byte atau 4*10(kolom)*70782(jumlah baris = jumlah trace).

Setelah itu anda siap untuk menuliskan geom_header.bin ke dalam trace header dari data anda dengan perintah di bawah ini (buat dalam satu baris):

sushw < Line_001_kill_vibro_agc_d2_fk_flt_decon.su infile=geom_header.bin key=sx,sy,selev,sstat,gx,gy,gelev,gstat,cdp,offset > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom.su

Dengan perintah

surange < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom.su

Perhatikan bahwa informasi geometry dan cdp sudah berada dalam trace header data seismik anda.

Untuk keperluan sorting dari shot gather ke CMP gather, dua informasi penting yang harus dimiliki adalah cdp dan offset.

Data ini memiliki cdp: 1 1282 (1 - 1282), jika angka yang berada di luar kurung sama dengan yang berada di dalam kurung maka penomoran cdp data ini benar, sedangkan jika kita lihat nilai offset: 12 3525 (3518 - 3509), rentang angkanya berbeda. Hal ini terjadi karena angka offset yang exact sangat sulit diperoleh akibat medan akuisi (lihat plot x-y coordinate sebelumnya). Oleh mengantisipasi hal ini, saya akan melakukan 'regularisasi' offset.

Kita mengetahui dari shot gather bahwa geometri akuisisi data tersebut adalah split-spread dengan interval geophone = 25m dan jarak dari sumber ke geophone pertama=25m.

Ketik gedit copy, paste kode di bawah ini dan save dengan nama regoff.m

clear; clc
%%%offset regularization...
right=[25:25:3525]';
left=[-3525:25:-25]';
offset_each_shot=[left;right];
offset_all_shot=repmat(offset_each_shot,251,1);

%%%mengganti offset pada geom_header.txt dengan offset yang baru
load geom_header.txt
geom_header_reg_offset=[geom_header(:,[1:9]),offset_all_shot];
save -ascii geom_header_reg_offset.txt geom_header_reg_offset
%%%%%

Ketik octave, lalu pada terminal octave ketik regoff

Keluar dari octave dengan ctrl+z, lalu pada terminal linux ketik:

a2b < geom_header_reg_offset.txt n1=10 > geom_header.bin

Ulangi perintah sushw di atas, lalu dengan surange diperoleh:

Dari hasil surange di atas, terlihat bahwa cdp dan offset memiliki nilai yang sama antara di luar kurung dan dalam kurung.

Tahap 12:

Pada tahap ini, kita akan melakukan sorting dari shot ke cmp gather dengan perintah susort berikut (buat dalam satu baris):

susort cdp offset <Line_001_kill_vibro_agc_d2_fk_flt_decon_geom.su > Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp.su

Lalu pilihlah cdp tertentu dan tampilkan:

suwind < Line_001_kill_vibro_agc_d2_fk_flt_decon_geom_cdp.su key=cdp min=1000 max=1002 | suxwigb perc=95

Sehingga diperoleh gambar sbb (klik untuk memperbesar):
Tahap berikutnya

Seismic Processing with Seismic Unix - Part4

2010-11-24T11:23:00.000-08:00

Lihat tahapan sebelumnya

Tahap 9:

Setelah kita melakukan eliminasi groundroll dengan F-K filtering, pada tahapan ini kita akan melakukan eliminasi noise-noise yang lainnya dengan menggunakan bandpass filter (sufilter).

Didalam bandpass filter terdapat empat parameter frekuensi yang harus kita tentukan f1,f2,f3 dan f4. Pemilihan rentang frekuensi yang akan diloloskan merupakan hal yang sangat krusial, jangan sampai proses ini menghilangkan data reflektor yang anda miliki. Untuk menghindari hal ini anda harus memperhatikan spektrum frekuensi serta rekaman reflektor yang dilihat pada shot gather.

Didalam terminal linux, ketiklah gedit & lalu copy, paste dan save kode berikut lalu beri nama filter

#!/bin/sh

#memilih shot dengan ep=80
suwind < Line_001_kill_vibro_agc_d2_fk.su key=ep min=80 max=80 > Line_001_kill_vibro_agc_d2_fk_ep80.su

#menampilkan shot dengan ep=80
suximage < Line_001_kill_vibro_agc_d2_fk_ep80.su perc=80 &

#memilih shot dengan ep=80 dan tracf 100-105
suwind < Line_001_kill_vibro_agc_d2_fk_ep80.su key=tracf min=100 max=105 > Line_001_kill_vibro_agc_d2_fk_ep80_tracf100_105.su

#spectrum untuk shot dengan ep=80 dengan tracf 100-105
suspecfx < Line_001_kill_vibro_agc_d2_fk_ep80_tracf100_105.su | suxwigb &

#melakukan filter untuk shot dengan ep=80
sufilter < Line_001_kill_vibro_agc_d2_fk_ep80.su f=10,15,50,60 > Line_001_kill_vibro_agc_d2_fk_ep80_flt.su

#menampilkan shot yang telah difilter
suximage < Line_001_kill_vibro_agc_d2_fk_ep80_flt.su perc=80 &

#memilih shot yang telah difilter dengan tracf 100-105
suwind < Line_001_kill_vibro_agc_d2_fk_ep80_flt.su key=tracf min=100 max=105 > Line_001_kill_vibro_agc_d2_fk_ep80_flt_tracf100_105.su

#spectrum untuk shot dengan ep=80 dengan tracf 100-105 setelah filter
suspecfx < Line_001_kill_vibro_agc_d2_fk_ep80_flt_tracf100_105.su | suxwigb &

Lalu ketik sh filter sehingga kita memperoleh gambar berikut:

Gambar di atas adalah shot gather dengan ep=80 (kiri) serta spektrum frekuensi untuk shot gather tersebut dengan tracf ke 100 sampai 105 (kanan). Perhatikan noise yang terdapat pada data, serta rentang frekuensi yang masih full bandwidth dari 8 sampai 100Hz (sumbu vertikal pada gambar sebelah kanan adalah frekuensi dalam Hz). Alasan saya memilih tracf ke 100 sampai 105 karena trace tersebut adalah mid offset, dari surange pada tahapan-tahapan sebelumnya kita mengetahui bahwa tracf yang kita miliki dimulai dari 1 s/d 282.

Gambar di bawah ini adalah hasil setelah difilter dengan bandpass filter 10,15,50,60Hz, yang berarti saya hanya meloloskan frekuensi dari 10 sampai 60Hz. Anda mungkin berargumen bahwa pemilihan cut off frekuensi 10Hz sangat membahayakan data seismik yang kita miliki, akan tetapi saya beralasan bahwa untuk kasus data ini, walaupun cut off 10Hz, reflektor yang saya miliki masih bisa terselamatkan (lihat shot gather). Tentu saja sebelumnya saya melakukan test dengan berbagai kombinasi frekuensi dari mulai 3,6,50,60Hz, 4,8,50,60Hz, 5,10,50,60Hz dan seterusnya.

Lakukanlah proses di atas untuk beberapa lokasi ep, jika anda sudah memperoleh hasil yang paling baik, maka anda dapat menerapkannya untuk seluruh data yang kita miliki dengan perintah sbb:

sufilter < Line_001_kill_vibro_agc_d2_fk.su f=10,15,50,60 > Line_001_kill_vibro_agc_d2_fk_flt.su

Tahap 10:
Pada tahapan ini kita akan menerapkan proses deconvolusi yang bertujuan untuk meningkatkan resolusi temporal dari reflektor serta menekan multiple. Namun sebelum melakukan deconvolusi, saya akan menerapkan autocorrelation terlebih dahulu yang sangat membantu mempelajari perilaku multiple.

Kode di bawah ini adalah kode untuk melakukan test autocorrelation dan deconvolusi. Parameter utama yang harus kita perhatikan adalah minlag dan maxlag, sedangkan ntout adalah jumlah sampel hasil autocorrelasi yang akan dihasilkan. Anda bisa melakukan test dengan pnoise yang berbeda.

Pada terminal linux ketik gedit, copy paste kode berikut lalu save dengan nama decon

#!/bin/sh

minlag=0.02
maxlag=0.1
pnoise=0.001
ntout=120

#memilih data dengan ep=150
suwind < Line_001_kill_vibro_agc_d2_fk_flt.su key=ep min=150 max=150 > Line_001_kill_vibro_agc_d2_fk_flt_ep150.su

#menampilkan data dengan ep=180
suximage < Line_001_kill_vibro_agc_d2_fk_flt_ep150.su perc=80 &

#melakukan autocorrelation dan menampilkan autocorrelation-nya
suacor < Line_001_kill_vibro_agc_d2_fk_flt_ep150.su suacor ntout=$ntout | suximage perc=80 &

#melakukan deconvolusi
supef < Line_001_kill_vibro_agc_d2_fk_flt_ep150.su > Line_001_kill_vibro_agc_d2_fk_flt_ep150_decon.su minlag=$minlag maxlag=$maxlag pnoise=$pnoise

#menampilkan hasil deconvolusi
suximage < Line_001_kill_vibro_agc_d2_fk_flt_ep150_decon.su perc=80 &

#melakukan autocorrelation dari data yang telah di-deconvolusi
suacor < Line_001_kill_vibro_agc_d2_fk_flt_ep150_decon.su suacor ntout=$ntout | suximage perc=80 &

Lalu pada terminal linux ketik sh decon

Gambar dibawah ini adalah gather sebelum deconvolusi (kiri) serta autocorrelation (kanan)

Untuk menentukan minlag dan maxlag, lihatlah penampang autocorrelation. minlag dihitung dari peak amplitude sampai zero crosing yang kedua. Hasil pengamatan saya menunjukkan bahwa peak amplitude data ini adalah sekitar 0.12s dan zero crossing yang kedua sekitar 0.14s, dengan demikian minlag=0.14-0.12=0.02, sedangkan maxlag 0.1 dipilih karena reverberasi masih terulang sampai sekitar 0.22s. sehingga maxlag=0.22-0.12=0.1s.

Gambar di bawah ini adalah hasil deconvolusi serta autocorrelation-nya. Perhatikan bentuk wavelet setelah deconvolusi yang lebih ramping (meningkat resolusi temporal), serta reverberasi yang sudah tereliminasi. Hal ini bisa kita lihat baik pada shot gather walaupun pada autocorrelation.

Setelah anda puas dengan hasil test, maka terapkanlah deconvolusi untuk seluruh data, dengan perintah sbb (satu baris):

supef < Line_001_kill_vibro_agc_d2_fk_flt.su > Line_001_kill_vibro_agc_d2_fk_flt_decon.su minlag=0.02 maxlag=0.1 pnoise=0.001

Tahap berikutnya

Seismic Processing with Seismic Unix - Part3

2010-11-23T19:20:00.000-08:00

Lihat tahapan sebelumnya

Tahap 8:

Sebagaimana yang kita perhatikan pada tahapan sebelumnya, shot gather yang kita miliki masih mengandung berbagai macam noise, diantara noise yang paling dominan adalah ground roll.

Untuk mengeliminasi ground roll, pada tahapan ini kita akan melakukan F-K filtering. Dimana F-K merupakan spectrum frekuensi (F) terhadap bilangan gelombang (K). Fungsi Seismic Unix untuk melakuan F-K filtering adalah sudipfilt.

Sebelum melakukan F-K filtering, data yang kita miliki harus memiliki sampling spatial (d2) yang pada hakikatnya merupakan jarak antara trace atau geophone interval (dalam hal ini 25 meter atau 0.025km).

sushw < Line_001_kill_vibro_agc.su key=d2 a=0.025 > Line_001_kill_vibro_agc_d2.su

Setelah memasukkan d2 terhadap trace header (perintah di atas), marilah kita melakukan test dengan memilih shot gather 32 (ep=32).

suwind < Line_001_kill_vibro_agc_d2.su key=ep min=32 max=32 > Line_001_kill_vibro_agc_d2_shot32.su

Didalam terminal ketiklah gedit & lalu copy-paste dan save kode berikut lalu berilah nama testfk

#!/bin/sh
slopes=-0.5,-0.3,0.3,0.5
amps=1,1,1,1
bias=0.0
dx=0.025
dt=0.002

sudipfilt < Line_001_kill_vibro_agc_d2_shot32.su dt=$dt dx=$dx \
slopes=$slopes amps=$amps bias=$bias |
suximage title="slope=$slopes amps=$amps bias=$bias" \
windowtitle="Shot 32" \
label1="Samples" label2="Trace" f1=1 d1=1 f2=1 d2=1 perc=80&

sudipfilt < Line_001_kill_vibro_agc_d2_shot32.su dt=$dt dx=$dx \
slopes=$slopes amps=$amps bias=$bias | suspecfk | suximage title="F-K Spectrum of Shot 32" \
windowtitle="F-K" label1="Frequency" label2="K" legend=1 cmap=hsv1 perc=97&

Lalu ketik pa da terminal linux sh testfk

Maka anda akan memperoleh gambar seperti di bawah ini. Data seismik sebelum F-K filter (kiri) dan Spektrum F-K (kanan). Klik untuk memperbesar gambar.

Parameter utama dari operasi F-K filtering adalah slopes, amps, bias. Slopes adalah kemiringan dari F-K spectrum, amps=1,1,1,1 (tidak melakukan filter), amps=1,0,0,1 (reject filter), amps=0,1,1,0 (pass filter). dx adalah geophone interval dalam km dan dt adalah time sampling interval dalam detik. Sedangkan bias adalah parameter aliasing. Bias=0 artinya data yang kita miliki tidak memiliki aliasing. Coba anda lakukan test jika bias=0,25 lalu lihat perbedaannya.

Gambar dibawah ini adalah hasil jika saya menggunakan:
slopes=-0.5,-0.3,0.3,0.5, amps=0,1,1,0, bias=0.0, dx=0.025, dt=0.002

Bandingkan shot gather diatas (setelah F-K) dengan shot gather sebelumnya (sebelum F-K), perhatikan ground roll yang sudah mulai tereliminasi.

Gambar dibawah ini adalah jika saya memilih amps=1,0,0,1 (reject filter). Kita melihat bahwa inilah noise yang kita kehendaki untuk di-reject.

Untuk analisa F-K, anda bisa melakukan test dengan berbagai macam slopes katakanlah slopes=-1,-0.5,0.5,1. Lihat hasilnya, lalu anda putuskan mana yang paling baik. Lakukanlah untuk beberapa lokasi ep, jika sudah memperoleh hasil yang terbaik dengan parameter tertentu, maka anda bisa melakukan F-K filter untuk seluruh shot gather yang kita miliki dengan mengetikkan kode berikut (buat satu baris).

sudipfilt < Line_001_kill_vibro_agc_d2.su dt=0.002 dx=0.025 slopes=-0.5,-0.3,0.3,0.5 amps=0,1,1,0 bias=0 > Line_001_kill_vibro_agc_d2_fk.su

Lihat berikutnya

Seismic Processing with Seismic Unix - Part2

2010-11-22T20:27:00.000-08:00

Lihat tahapan-tahapan sebelumnya

Tahap 6:

Pada tahapan ini, saya akan menghilangkan (kill), trace vibroseis dengan key=tracf -1 dan 0. Dari surange, kita mengetahui bahwa rentang tracf dimulai dari -1,0,1,2,3...,282. Perintah berikut adalah cara untuk tidak melibatkan trace vibroseis yang tertanam pada shot gather.

suwind < Line_001.su key=tracf min=1 max=282 > Line_001_kill_vibro.su

Coba tampilkan dengan perintah berikut lalu zoom in dengan meng-klik left button, tahan dan geser untuk memastikan trace vibroseis telah hilang (untuk meng-unzoom, klik pada window xwigb).

suwind < Line_001_kill_vibro.su key=ep min=32 max=32 | suxwigb perc=80

Lalu saya ingin menampilkan shot 32 dalam bentuk image,

suwind < Line_001_kill_vibro.su key=ep min=32 max=32 | suximage perc=80

Dari gambar di atas terlihat sebuah rekaman yang masih penuh dengan noise seperti ground roll, air blast, direct wave, coherent noise, trace yang tidak koheren (time shift), amplitudo yang tidak sama antara zona dangkal dan dalam (akibat geometrical spreading), dll.

Tahap 7:
Perintah berikut ini adalah cara untuk mengkompensasi penurunan amplitudo dengan AGC (Automatic Gain Control). Perlu diingat bahwa AGC merupakan operasi trace by trace, bisa digunakan hanya untuk interpretasi bukan untuk analisis AVO. Untuk analisis AVO sebaiknya digunakan dB/sec gain setelah koreksi NMO.

sugain < Line_001_kill_vibro.su agc=1 wagc=0.2 > Line_001_kill_vibro_agc.su

Tampilkan:

suwind < Line_001_kill_vibro_agc.su key=ep min=32 max=32 | suximage perc=80

Gambar di atas adalah shot gather AGC

Tahapan berikutnya

Seismic Processing with Seismic Unix - Part1

2010-11-20T07:28:00.000-08:00

Sebagaimana yang kita ketahui, Seismic Unix memiliki ratusan fungsi yang berguna untuk mengolah data seismic.

Untuk melakukan latihan pengolahan data seismik, berikut adalah tahapan yang bisa diikuti:

Tahap1:
Download data seismik 2 dimensi dengan format segy (388.9MB) di sini

Tahap2:
Simpanlah file tersebut dalam folder tertentu, katakanlah processing, lalu ekstrak sehingga kita akan memperoleh beberapa file berikut: Line_001.sgy, Line_001.SPS, Line001.RPS, Line_001.XPS, dan Line_001.TXT. File dengan ektensi SPS dan RPS berisikan informasi sumber-penerima seperti indeks (nomor), koordinat, elevasi, statik, dll. File XPS berisikan informasi hubungan sumber-penerima (relational) dan TXT berisikan informasi tentang parameter survey. Berikut adalah isi dari Line_001.TXT:

Tahap 3:
Untuk melakukan processing dengan Seismic Unix, terlebih dahulu kita melakukan konversi format data seismik dari segy ke su.

segyread tape=Line_001.sgy verbose=1 endian=0 > Line_001.su

Pada perintah di atas, saya mendefiniskan endian=0 karena saya menggunakan mesin little endian (Laptop/PC). Jika anda menggunakan workstation anda harus medefinisikan endian=1.

Tahap 4:
Analisa header data su dengan perintah

surange < Line_001.su

Sehingga diperoleh:

Dari informasi di atas terlihat bahwa data seismik ini memiliki jumlah trace 71284, ep adalah shot number dari 32 sampai 282 = 251 shots, serta indeks trace dalam setiap shot tracf dari -1 sampai 282. Dengan kata lain jumlah trace dalam setiap shot adalah 284.

Tahap 5:
Untuk mengevaluasi data yang kita miliki lakukan peritah berikut:

suwind < Line_001.su key=ep min=32 max=32 | suxwigb perc=80

Perintah di atas, kita memilih data dengan suwind untuk shot ke 32, lalu ditampilkan sebagai wiggle dengan suxwigb dengan persentase amplitudo 80.

Gambar di bawah ini (kiri) adalah hasil dari perintah di atas, serta zoom in (kanan) dari trace-trace awal. Terlihat bahwa 2 trace pertama adalah data source signature yang didefinisikan dengan tracf=-1 dan tracf=0.

Lakukanlah proses ini untuk beberapa tempat dengan ep berbeda. Kesimpulan yang saya peroleh adalah semua shot memiliki 284 trace dengan 2 trace pertama sebagai source signature yang ditanamkan pada setiap shot record.

Tahapan berikutnya

Madagascar

2010-11-18T15:24:00.000-08:00

Madagascar adalah sebuah open source yang berisi ratusan tools yang berguna untuk melakukan data analisis multi dimensi seperti perilaku gelombang seismik, ray tracing, waveform modeling, pengolahan citra, termasuk didalamnya fungsi-fungsi pengolahan data seismik.

Berikut adalah tahapan Instalasi Madagascar dengan platform Linux Ubuntu 10.10.

Download Madagascar untuk Linux di sini
Simpanlah file yang telah didownload tersebut pada directory yang bernama madagascar (contoh).
Masuk ke directory madagascar lalu ketik tar -xvvf madagascar-1.0.tar.gz
Semua file yang telah diekstrak akan disimpan di dalam directory baru yang bernama madagascar-1.0.
Pada terminal linux ketik:
sudo apt-get install freeglut3-dev g++ gfortran libc6-dev libg2-dev libglew1.5-dev \
libjpeg62-dev libx11-dev libxaw7-dev libnetpbm10-dev swig python-dev python-scipy \
python-numpy scons units
Masih pada terminal linux ketik
export RSFROOT=/home/agus/madagascar/madagascar-1.0
export PYTHONPATH=$PYTHONPATH:$RSFROOT/lib
export PATH=$PATH:$RSFROOT/bin
Masuk ke directory madagascar, lalu ketik:
./configure lalu enter scons; scons install lalu enter source env.sh lalu enter
Test program dengan mengetikkan pada terminal (satu baris)
sfspike n1=1000 k1=300 title="\s200 Welcome to \c2 RSF" | sfbandpass fhi=2 phase=1 | sfwiggle | sfpen
Jika instalasi dan konfigurasi anda benar, maka anda akan memperoleh gambar seperti di bawah ini.

Setelah reboot atau keluar terminal, kadang-kadang command madagascar tidak dikenal lagi oleh sistem dengan munculnya error message: "....command not found" atau error message yang lainnya, untuk mengantisipasi hal ini, coba ketik lagi kode pada tahap 6 di ./bashrc dengan mengetikkan gedit ~/.bashrc

Contoh di bawah ini saya melakukan wavefront modeling 2D dengan madagascar, dimana kodenya bisa dilihat di bawah artikel ini.

Berikut adalah tahapan untuk menjalankan kode tersebut.

Buatlah directory baru katakanlah /demo/2dmodel
Jika saya ketik pwd saya akan memiliki path /home/agus/madagascar/madagascar-1.0/demo/2dmodel
Pada path tersebut ketik gedit
Lalu copy paste kode pemodelan di bawah dan save dengan nama SConstruct
Ketik scons pada terminal linux, sehingga beberapa file akan digenerate
Bukalah directory Fig i.e. /home/agus/madagascar/madagascar-1.0/demo/2dmodel/Fig
Pada directory tersebut saya ketik sfpen < mod1.vpl

Jika saya ketik sfpen < wave1.vpl
Maka saya akan memiliki movie penjalaran gelombang dimana snapshotnya seperti ini:

Berikut adalah animasinya (klik untuk memperbesar):

Untuk mempelajari lebih jauh tentang Madagascar, anda bisa mempelajarinya disini

Jika anda bekerja pada ubuntu 11.04 anda dapat menginstall madagascar 1.2 .

Kode Pemodelan Wavefront:

##########################################################
# Created by Kyle Shalek, slightly modified by Agus Abdullah #
# Last Modified: 9-4-09 #
# Creates a 2D model with velocity and density inputs. #
# Acoustic FD modeling with awefd2d (Psava) #
# Single trace output #
##########################################################

from rsf.proj import *
import string
########################################################
#1) create velocity model
# Depth in km
xmax = 4.0
zmax = 2.0

layers = ((0.1,0.1,0.1,0.1,0.1),
(0.1,0.4,0.7,1.0,1.3),
(1.5,1.5,1.5,1.5,1.5))

layer1=0.1
layer2=1.3-0.1
layer3=1.5-1.3
layer4=2-1.5

# Velocity (Vp) in km/s
vp = (0.0,
6,
4.6657,
6)
#Velocity (Vs) in km/s
vs = (0.0,
2,
4,
2)

# Density in g/cc, converted to kg/ckm
densities = (0.1*1000000000000,
2.8*1000000000000,
2.5*1000000000000,
2.8*1000000000000)

def arr2str(array,sep=' '):
return string.join(map(str,array),sep)

n1 = len(layers[0])
n2 = len(layers)

Flow('layers',None,
'''
echo %s
n1=%d n2=%d o1=0 d1=%g
data_format=ascii_float in=$TARGET
''' % (string.join(map(arr2str,layers),' '),
n1,n2,xmax/(n1-1)))

d = 0.01 # sampling steps

Flow('mod1','layers',
'''
spline o1=0 d1=%g n1=%d |
unif2 d1=%g n1=%d v00=%s |
dd form=native
''' % (d,int(1.5+xmax/d),
d,int(1.5+zmax/d),
arr2str(vp,','),))
Flow('rho','layers',
'''
spline o1=0 d1=%g n1=%d |
unif2 d1=%g n1=%d v00=%s |
dd form=native
''' % (d,int(1.5+xmax/d),
d,int(1.5+zmax/d),
arr2str(densities,','),))

Result('mod1',
'''
grey color=j title="Velocity Model"
allpos=y titlesz=8 labelsz=6 screenratio=0.5
scalebar=y barlabel='velocity (km/s)' barlabelsz=6
label1="Depth (km)"
label2="Distance (km)"
''')

# Contour the layers
Plot('modline','mod1','contour title=" " label1='' label2='' wantaxis=n screenratio=0.5')
##############################################################

#2) 2D Acoustic wave FD Modeling
# Source location
sx = 2.5
sz = 0.1
Flow('source',None,
'spike n1=3 nsp=3 k1=1,2,3 mag=%g,%g,1 o1=0 o2=0' % (sx,sz))
# Receiver location
rx=2
rz=0.1
Flow ('receiver',None,
'spike n1=3 nsp=3 k1=1,2,3 mag=%g,%g,1 o1=0 o2=0' % (rx,rz))

########################
time=0.5
timesteps=1500 # timesteps*0.001=sec
frequency=30

##########################

#2)a) Source wavelet
Flow('wavelet',None,
'''
spike nsp=1 n1=%g d1=%g k1=%d |
ricker1 frequency=%g |
transp
''' % (timesteps,0.001,200,frequency))

#2)b) Awefd modeling
Flow('data wave','wavelet mod1 source receiver rho',
'''
awefd2d verb=y free=y expl=y snap=y jsnap=10 dabc=y
db=0 dbx=0 o1=0 o2=0
vel=${SOURCES[1]}
sou=${SOURCES[2]}
rec=${SOURCES[3]}
den=${SOURCES[4]}
wfl=${TARGETS[1]}
''')

#################################################################

#2)c) Movie of wave
Plot('wave',
'''
grey gainpanel=all title=Wave label1=Depth unit1=km label2=Lateral unit2=km
''',view=1)
Result('wave1','wave',
'''
grey color=j scalebar=y barlabel='Amplitude' barlabelsz=6
gainpanel=all title=Wave
label1=Depth unit1=km label2=Lateral unit2=km screenratio=0.5
''')

# Wavefield snapshot at 'time' value
Plot('snap','wave',
'''
window n1=201 n2=401 n3=1 min3=%g |
grey color=j gainpanel=all title="Wave Snapshot"
label1=Depth unit1=km label2=Lateral unit2=km screenratio=0.5
scalebar=y barlabel='Amplitude' barlabelsz=6
'''%(time))
##################################################################

#2)d) Overlay wave and layer contours
Result('snap1',['snap','modline'],'Overlay')
#2)e) Single Trace
Result('trace','data',
'''
transp |
spline o1=0 d1=0.001 n1=%g |
wiggle label1=Time unit1=s label2=Amplitude unit2='' title='Vertical Trace with Source'
'''% (timesteps))
Result('trace_window','data',
'''
transp |
spline o1=0 d1=0.001 n1=%g |
window min1=0.3 |
wiggle label1=Time unit1=s label2=Amplitude unit2='' title='Vertical Trace without Source'
'''% (timesteps))
End()
################################################################################

Crooked Line

2010-11-15T18:47:00.001-08:00

Crooked line adalah lintasan seismik dimana posisi sumber dan penerima tidak berada disepanjang garis lurus. Sebuah lintasan seismik dikatakan crooked jika lintasan tersebut membentuk sudut diatas 10 derajat.

Sebagaimana yang kita ketahui, jika lintasan seismik lurus-lurus saja, maka midpoints (titik semu di bawah permukaan) akan berada tepat di bawah lintasan seismik, sedangkan untuk crooked line, midpoints tersebut akan tersebar secara acak di sekitar garis yang membengkok (crooked).

Gambar di bawah ini mengilustrasikan lokasi midpoints yang tersebar secara acak pada crooked seismic line (wilayah abu-abu).

Courtesy Western Geco

Jika kita melakukan sorting berdasarkan CMP (Common Mid Point), maka titik-titik semu di bawah permukaan tersebut akan membentuk sebuah wilayah tertentu daripada sebuah titik seperti yang lazimnya kita peroleh dari lintasan seismik yang lurus. Hal inilah akan menghasilkan efek yang populer disebut dengan CMP smearing. Untuk mengurangi efek CMP smearing, harus dilakukan pendefinisian lintasan baru melalui midpoints yang tersebar tadi, sehingga kita memperoleh crooked profile.

Crooked profile dibuat melalui midpoints dengan interval sumber-penerima (offset) tertentu. Pada praktiknya, crooked profile dikelompokkan menjadi beberapa interval berdasarkan jarak sumber-penerima (offset). Ujung akhir dari kelompok tersebut menjadi pusat CMP yang baru yang mewakili populasi midpoints tadi.

Courtesy Western Geco

Gambar di atas mengilustrasikan proses penyambungan midpoints dengan offset 100 meter untuk menghasilkan crooked profile dengan threshold offset 100m. Masing-masing offset akan menghasilkan crooked profile yang berbeda. Dengan menetapkan rentang offset tertentu, sebuah crooked profile akhir akan diperoleh melalui perata-rataan semua crooked profile dari semua rentang offset.

Gambar dibawah ini mengilustrasikan penampang crooked line serta penampang tiga dimensi dari sebaran midpoints.

Courtesy Nedimovic M.R. and West G.F., GEOPHYSICS, VOL. 68, NO. 2003

Gambar di bawah ini menunjukkan distribusi midpoints (garis yang tegak lurus dengan lintasan) untuk sebuah crooked line. Banyaknya midpoints pada setiap lokasi ditentukan oleh lebar garis yang tegak lurus tersebut. Perhatikan efek crooked line pada lintasan yang menyudut (elips merah), panjangnya garis lurus tersebut umumnya lebih pendek dan tidak merata. Jika kita melakukan koreksi crooked line (seperti prosedur di atas), maka distribusi midpoints tersebut akan lebih merata di seluruh lintasan seismik (tentu saja setelah melewati zona roll-on roll-off).

Modified from Hall K.W. et al., Crewes

Seismik Refraksi

2010-11-07T08:17:00.000-08:00

Studi seismik refraksi ditujukan untuk memetakan karakteristik lapisan dekat permukaan (near surface) seperti kedalaman lapisan lapuk (weathering), bed rocks, pemetaan air tanah, lingkungan, dll. Informasi geofisika yang diperoleh dari studi ini adalah model kecepatan serta kedalaman lapisan bawah permukaan. Informasi tersebut diturunkan dari first break serta geometri sumber-penerima.

Peralatan yang digunakan didalam survey seismik refraksi, biasanya terdiri dari 12 sampai 24 channel geophone dengan interval 2-5 meter dan frekuensi 8-14Hz, dengan sumber gelombang berupa palu ataupun dinamit serta perekam yang biasanya jauh lebih portable daripada peralatan seismik refleksi. Akan tetapi pada sebuah survey seismik refleksi, rekaman refraksi bisa diperoleh seiring dengan perekaman gelombang refleksi sendiri.

Gambar dibawah ini menunjukkan perekam yang dikhususkan untuk survey seismik refraksi.

Courtesy Geometrics

Gambar di bawah ini adalah layout perekaman seismik refraksi. Geophone diletakkan disepanjang lintasan survey, dimana offset (bentangan kabel) harus 3-5 kali lebih panjang dari kedalaman target. Jadi jika panjang offset nya adalah 600 meter, maka kedalaman maksimum yang akan terdeteksi adalah 200 meter.

Courtesy Geometrics

Sumber gelombang pada sebuah offset dari survey refraksi, sedikitnya dua sumber S1 dan S2 (lihat gambar di bawah) yang biasanya diletakkan di sisi kiri dan kanan, dengan jarak dari geophone ½ dari group interval. Ada baiknya juga dilakukan penembakan S3. Jika crossover distance diketahui (lihat gambar anonymous (GEOL 335.3)). Ada baiknya dilakukan pengukuran S4 dan S5, tergantung tujuannya S6 dan S7 kadang-kadang diperlukan.

Courtesy Geometrics

Pada rekaman seismik (shot gathers), first break merupakan sinyal yang pertama kali terekam oleh penerima. Sinyal tersebut berasal dari direct wave dan head wave. Direct wave adalah gelombang yang merambat dari sumber langsung ke penerima melewati lapisan pertama, Sedangkan head wave adalah gelombang yang melewati lapisan pertama lalu merambat disepanjang lapisan kedua. Syarat terjadinya head wave adalah sudut tembak gelombang harus melewati critical angle dan lapisan kecepatan lapisan tersebut harus lebih cepat dari lapisan sebelumnya.

Berikut adalah ilustrasi jejak sinar, kurva serta persamaan waktu tempuh dari direct wave (merah), head wave (biru) dan refleksi (hijau).

Courtesy anonymous (GEOL 335.3)

Gambar di bawah ini menunjukkan rekaman (shot gather) serta interpretasi first break untuk direct wave (merah), head wave yang merambat melewati lapisan pertama dan disepanjang lapisan kedua (biru), serta head wave yang melewati lapisan pertama, kedua dan disepanjang lapisan ketiga (hijau). Kedalaman dan kecepatan lapisan pertama dapat dianalisis dari kurva warna merah, lapisan kedua dari kurva warna biru dan lapisan ketiga dari kurva warna hijau. Perhatikan, banyaknya perlapisan ditunjukkan dengan berapa banyak kurva tersebut saling memotong (crossover).

Courtesy Mitchell J.F. and R.J. Bolander, 1986

Gambar dibawah ini menunjukkan perhitungan kecepatan dan kedalaman perlapisan dari kurva waktu tempuh terhadap offset untuk model 3 lapisan bumi berlapis horizontal.

Courtesy Exploration Geophysics: Refraction Seismic Notes 06/20/02

Kecepatan lapisan V1, V2, dan V3 merupakan slope dari masing-masing kurva, sedangkan kedalaman lapisan ditentukan oleh persamaan di bawah ini:

Courtesy Exploration Geophysics: Refraction Seismic Notes 06/20/02

Karakteristik kemiringan perlapisan batuan, dapat dianalisis melalui kesimetrisan kurva waktu tempuh dari beberapa shot, setidak-tidak nya S1 dan S2 pada desain survey di atas.

Gambar di bawah ini merupupakan plot waktu tempuh untuk semua shot. Terlihat bahwa kurva yang dari arah kiri simetris dengan kurva dari arah kanan, sehingga kita berkesimpulan bahwa perlapisan pada survey ini adalah horizontal.

Courtesy Geometrics

Gambar di bawah ini menunjukkan karikatur survey seismic refraksi untuk model 3 perlapisan yang miring. Sumber gelombang diletakkan pada shot a dan shot b.

Courtesy Exploration Geophysics: Refraction Seismic Notes 06/20/02

Masing-masing kurva waktu tempuh dari shot a dan shot b ditunjukkan pada gambar di bawah ini:

Courtesy Exploration Geophysics: Refraction Seismic Notes 06/20/02

Kecepatan lapisan pertama V1 dihitung dari slope berwarna hijau, dimana kedua shot akan memiliki nilai yang sama. Sedangkan untuk V2 (merah) dan V3 (pink), masing-masing merupakan rata-rata dari slope shot a dan slope shot b. Kedalaman lapisan dibawah masing-masing shot dihitung dengan rumus di atas, yang tentunya akan menghasilkan kedalaman yang berbeda untuk shot a dan shot b. Kemiringan lapisan ditentukan dari perbedaan kedalaman tersebut.

Static Correction

2010-11-05T20:58:00.001-07:00

Didalam pengolahan data seismik, terdapat dua jenis static correction (koreksi statik) yang harus dilakukan yakni elevation (field) statics dan near surface (weathering) correction.

Elevation Statics
Elevation statics adalah koreksi karena perbedaan elevasi source dan receiver. Elevation statics dilakukan dengan meletakkan source dan receiver pada posisi virtual dengan elevasi yang sama (datum) yang biasanya sedikit dibawah elevasi source dan receiver yang terendah. Untuk proses ini, diperlukan informasi replacement velocity dari material antara datum dengan masing-masing source dan receiver. Replacement velocity biasanya diperoleh dari pengetahuan sebelumnya dari daerah yang bersangkutan atau dari pengukuran uphole time.

Persamaan elevation statics, diberikan oleh:

t_D = [(E_S – Z_S - E_D) + (E_R – Z_R - E_D)]/V_r

Dimana E_S elevasi dari source (di atas permukaan laut), Z_S kedalaman dari source ( 0 untuk vibroseis), E_R elevasi dari receiver, Z_R kedalaman dari receiver, E_D elevasi datum, dan V_radalah replacement velocity.

Gambar di bawah ini mengilustrasikan parameter-parameter yang digunakan dalam elevation statics.
Selanjutnya dilakukan perhitungan:

TWTes=TWT-t_D

TWTes adalah waktu tempuh (TWT) dari trace pasangan sumber penerima setelah koreksi statik, TWT adalah waktu tempuh trace sebelum koreksi statik.

Near surface (weathering) correction
Setelah elevation statics, selanjutnya dilakukan near surface correction untuk mengkompensasi variasi ketebalan dan perubahan lateral kecepatan dari lapisan lapuk (weathering zone).

Metoda yang biasa dilakukan untuk near surface correction adalah: uphole surveys, refractions statics, residual statics dan tomo statics.

Uphole survey dilakukan dengan meletakkan rangkaian geophone pada lubang bor yang menembus kedalaman weathering layer dan sub weathering layer, sumber gelombang diletakkan di permukaan di dekat lubang bor, lalu waktu tempuh gelombang dari sumber ke masing-masing penerima diplot untuk menghitung kecepatan. Uphole survey dilakukan pada beberapa lokasi, dimana kecepatan lapisan akan diinterpolasi dari satu lokasi ke lokasi lainnya.

Refraction statics dilakukan untuk mengestimasi long wavelength statics. Wavelength statics mengacu pada lebarnya perubahan lateral kecepatan dan ketebalan dari lapisan lapuk relatif terhadap offset maksimum.

Refraction statics dilakukan untuk mengestimasi ketebalan dan kecepatan lapisan lapuk pada lintasan survey. Terdapat beberapa metoda refraction statics diantaranya: Delay-time, Generalized Reciprocal Method (GRM), dan Least-squares method. Dua metoda pertama memerlukan proses first break picking dan geometri jejak sinar, sedangkan least square method memiliki konsep yang sama dengan residual statics, hanya saja least square menggunakan gelombang refraksi daripada refleksi.

Gambar di bawah ini menunjukkan hasil pekerjaan first break picking, profil weathering dan sub weathering zone dari studi seismik refraksi serta profil untuk elevation statics dan weathering statics.

Courtesy Lawton, D.C., Geophysics Vol. 54 NO. 10, 1989.

Dari hasil studi refraksi di atas terlihat bahwa kecepatan weathering zone (Vw) = 520m/s sedangkan sub weathering zone (Vb)=3103m/s. Demikian juga kita memperoleh informasi kedalaman weathering zone (Zw) pada setiap titik pengamatan.

Koreksi statik untuk mengkompensasi weathering zone terhadap datum D dapat dirumuskan sbb:

t_D=(-2Zw/Vw)+(2(E_D-E_S+Zw)/Vb)

Dimana E_Dadalah elevasi dari datum dan E_Selevasi dari source. Selanjutnya t_D akan dikurangkan dari TWT trace untuk shot yang bersangkutan.

Gambar di bawah ini menunjukkan perbandingan antara data awal (a), elevation statics (b), dan elevation statics+weathering statics (c). Walaupun pada contoh ini kita tidak melihat perbedaan yang sangat mencolok, tetapi kalau kita perhatikan lebih teliti koreksi statics memberikan peningkatan koherensi dari satu trace ke trace yang lainnya.

Courtesy Lawton, D.C., Geophysics Vol. 54 NO. 10, 1989.

Residual statics atau dikenal juga dengan reflection statics dilakukan untuk mengkoreksi short wavelength statics, dimana metoda yang paling banyak digunakan adalah surface consistent method.

Residual statics dirumuskan dengan persamaan sbb:

t_ijk=r_i+s_j+G_k+M_kx_ij²

dimana t_ijk adalah koreksi residual statics total untuk receiver ke-i, source ke - dan CMP ke-k, r_i adalah residual statics untuk receiver ke-i, s_j adalah residual statics untuk source ke-j, G_k perbedaan TWT antara CMP ke-k dengan CMP referensi yang biasanya CMP pertama, M_kx_ij² adalah residual moveout.

Mengacu pada persamaan di atas, komputasi least square dilakukan untuk memperoleh parameter r_i, s_j, G_k dan M_k sebagai unknown parameters dengan input t_ijk dan x_ij (offset) sebagai known parameters.

Pada praktiknya residual statics dilakukan dengan beberapa tahap: menghitung t_ijk, mendekomposisi persamaan di atas untuk memperoleh unknown parameters, menerapkan r_i dan s_j pada gathers sebelum koreksi NMO.

Cara yang umum digunakan untuk menghitung t_ijk adalah pilot trace method, dimana tahapannya sbb:

Sebuah CMP dengan S/N ratio yang bagus di-gain lalu di-NMO dengan menggunakan model kecepatan awal.
Memilih event refleksi tertentu.
CMP gather di stack untuk menghasilkan sebuah stack trace.
Masing-masing trace pada gather di kros-korelasi dengan stack trace.
Time shifts t_ijk¹ merupakan time shift dengan kros-korelasi tertinggi.
Melakukan shifting trace original dengan t_ijk¹ .
Pilot trace awal diperoleh dengan melakukan stacking trace yang di time shifting pada gather tersebut.
Pilot trace ini di kros-korelasi dengan trace yang di shifting pada gather untuk memperoleh t_ijk².
Shift masing-masing trace yang telah di shift dengan menggunakan t_ijk².
Time shift total saat ini adalah t_ijk = t_ijk¹+t_ijk²
Pilot trace akhir dibangun lagi dengan melakukan stacking time shift total.
Pilot trace akhir ini di kros korelasi dengan trace-trace pada gather berikutnya untuk membuat pilot trace pada gather yang bersangkutan.
Proses ini dilakukan pada semua CMP dengan menggerakkan ke-kanan atau ke kiri dari CMP referensi.

Gambar di bawah ini menunjukkan stack sebelum dan setelah residual statics dengan referensi events pada 3.5-4.0s . Perhatikan perbedaan antara keduanya.

Sebelum reflection statics, courtesy Ronen J. dan Claerbout, SEG Atlanta, 1984

Setelah reflection statics, courtesy Ronen J. dan Claerbout, SEG Atlanta, 1984

Pada hasil studi seismik refraksi di atas (Lawton, D.C), kita melihat bahwa profil model kecepatan untuk weathering dan sub weathering zone memiliki resolusi lateral yang kurang baik. Hal seperti ini memberikan hasil yang kurang akurat untuk mengkoreksi elevation statics maupun refraction statics.

Untuk mengantisipasi kekurangan tersebut, teknik tomo statics digunakan untuk memberikan hasil yang lebih baik karena menghasilkan model kecepatan yang lebih detail.

Hal yang terpenting pada tomo statics adalah determinasi model kecepatan tomography yang secara garis besar diperoleh dengan tahapan sbb: penelusuran jejak sinar refraksi (ray tracing), mendiskritisasi model bawah permukaan menjadi sel-sel tomography, perhitungan waktu tempuh gelombang refraksi (first break picking), pembuatan model kecepatan awal dari studi refraksi dan inversi untuk memperoleh model kecepatan tomography.

Keakuratan model kecepatan tomography terutama ditentukan oleh densitas sinar dan arah sinar yang merata ke semua arah pada masing-masing sel tomography.

Gambar di bawah ini menunjukkan densitas sinar serta model kecepatan refraction tomography.

Densitas Sinar, courtesy Doll et al., 2010

Model kecepatan refraction tomography, courtesy Doll et al., 2010