Programming,
Prosesing Data Seismik menggunakan Matlab (Part 1) - Load data, otak-atik dan cek geometri
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)')
>>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
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
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
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')
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
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)
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
mas bisa request upload ebooknya ga?/
ReplyDelete