Programming,

Prosesing Data Seismik menggunakan Matlab (Part 1) - Load data, otak-atik dan cek geometri

10:55:00 PM Leo Cahya D 1 Comments

Hello, digital world..

Hari ini mau saya mau iseng-iseng mengolah data seismik menggunakan Matlab setelah melihat buku karangan Wail A Mousa & co., Processing of Seismic Reflection Data Using MATLAB. Kode-kode dan data untuk latihan sudah dibuat beliau dan dapat diunduh di link berikut.

Oke setelah mendapatkan data bernama Book_Seismic_Data.mat tersebut, kita buka dengan perintah :

>>load Book_Seismic_Data.mat

Pada workspace matlab dapat dilihat ada dua buah variabel dari file tersebut, yaitu D untuk matriks rekaman data dengan panjang trace 594 dan panjang sampel rekaman 150, lalu variabel H yang berisi header dari data tiap trace tersebut.

Menampilkan Rekaman
Untuk menampilkan data, pertama-tama harus kita ketahui waktu cuplik (sampling time) rekaman terlebih dahulu dari variabel H dengan menggunakan perintah :

>>dt=H(1).dt/1000/1000;

kenapa perlu dibagi 1e-6 ? karena umumnya header dt itu dalam satuan µsec.
Selanjutnya untuk menampilkan data anda dapat secara manual menggunakan perintah :
>> x=1:size(D,2); z=0:dt:(size(D,1)-1)*dt; 
>> pcolor(x,z,-1*D); shading interp; axis ij; colormap(gray); colorbar;
>> xlabel('trace number'); ylabel('time (s)');
Plot data dengan colormap gray, dapat diketahui jumlah shotnya ada 18

Jika ingin menggunakan warna umum RGB rekaman seismik :
>> x=1:size(D,2); z=0:dt:(size(D,1)-1)*dt; 
>> pcolor(x,z,-1*D); shading interp; axis ij; colormap(seis_colors); colorbar;
>> xlabel('trace number'); ylabel('time (s)');


Dan juga dari program milik Mousa dapat menampilkan wigglenya dengan perintah
>> x=1:size(D,2); z=0:dt:(size(D,1)-1)*dt; 
>>mwigb(D,1,x,z); axis([0 600 0 3]);
>>xlabel('trace number'); ylabel('time (s)') 


Dari gambar di atas dapat diketahui bahwa jumlah shot yang nya ada 18 buah (bisa juga dilihat nomor fldr pada variabel H dengan perintah "H(end).fldr").

Menampilkan Desain Survei Posisi Shot, Receiver, dan Stacking Chart
Pertama-tama kita ektrak terlebih dahulu nilai geometri dari header, anda bisa menggunakan fungsi extract_geometry.m milik Mousa atau run perintah berikut

>>sx=[H.sx]; sy=[H.sy]; gx=[H.gx]; gy=[H.gy]; 
sz=[H.selev]; gz=[H.gelev]; sg=[H.fldr];
count=1; l=1; k=2;
for i=1:length(sg)
    if k<=length(sg) && sg(k-1)==sg(k)
        count=count+1;
    else
        num_trace_per_sg(l)=count;
        num_sg(l)=l;
        count=1;
        l=l+1;
    end
    k=k+1;
end

dimana variabel tersebut merupakan :
sx: sources x-axis locations
sy: sources y-axis locations
gx: receivers x-axis locations
gy: receivers y-axis locations
num_sg: number of shot gathers
num_trace_per_sg: number of traces/shot gather
sz: elevations of the sources
gz: elevations of the receivers

Untuk melihat jumlah trace tiap shot, dapat dilihat menggunakan perintah :
>>figure,stem(shot_gathers,num_trace_per_sg)
xlabel('Shot gather numbers','FontSize',14)
ylabel('Number of traces/shot gather','FontSize',14)
axis([0,max(shot_gathers)+1,0,max(num_trace_per_sg)+2])
set(gca,'YMinorGrid','on')

dapat dilihat bahwa tiap shot memiliki 33 trace. Selanjutnya posisi ketinggian receiver dan shot dapat dilihat dengan perintah :

>>figure,plot(sz,'*-')
xlabel('Number of traces','FontSize',14)
ylabel('Sources elevation (ft)','FontSize',14)
grid

>>figure,plot(gz,'*-')
xlabel('Number of traces','FontSize',14)
ylabel('Receivers elevation (ft)','FontSize',14)
grid

dan stacking chart dapat kita lihat dengan perintah :
>>num_shots=length(shot_gathers);
stacking_chart(sx,gx,num_shots,num_trace_per_sg)


untuk selanjutnya, akan kita bahas tentang editing dan QC data, kill trace seismik yang rekamannya kacau dan sebagainya..

see ya,

L

1 comment: