Ray Tracing

Penelusuran jejak sinar (ray tracing) merupakan proses yang sangat penting didalam aktifitas seismik eksplorasi seperti untuk keperluan desain survey, seismic modeling, 4D seismic, seismic tomography, dll.

Teknik ray tracing yang digunakan tergantung kepada kebutuhan dan kompleksitas model bawah permukaan. Untuk model bumi berlapis, ray tracing dapat dilakukan dengan mengikuti Hukum Snellius.

Gambar di bawah ini menjelaskan penelusuran jejak sinar yang melewati lapisan-lapisan bumi pada kedalaman z dengan nilai kecepatan v serta sudut datang dan transmisi sinar.

Courtesy Khanh Duc Nguyen, University of the Faroe Islands

Hubungan antara sudut datang gelombang, sudut transmisi dan kecepatan gelombang ditunjukkan oleh persamaan di atas. Untuk masing-masing sinar akan memiliki ray parameter p tertentu yang sama untuk semua lapisan.

Jarak lateral sinar (offset) dx dan waktu tempuh dt pada masing-masing lapisan ditunjukkan oleh persamaan sbb:

Dengan menjumlahkan seluruh dx dan dt, maka diperoleh offset dan waktu tempuh untuk masing-masing sinar.

Gambar di bawah ini menunjukkan ray tracing untuk model bumi parallel yang terdiri dari 9 lapisan serta kurva waktu tempuh untuk masing-masing lapisan.

Berikut ini adalah code matlab untuk ray tracing di atas.
%%%ray tracing seismik refleksi untuk model bumi berlapis horizontal dengan
%%%by Agus Abdullah, PhD
%%%just for ensiklopediseismik
clear; clc
for z=1:9
lap=[1:9];
nolayers=lap(z);
norays=15; %untuk setting theta di bawah, hati-hati dengan jumlah rays terhadap critical angle!!!
vel=[1500,1800,2200,1850,2400,2000,2700,2000,2900]; %kecepatan setiap lapisan
dz=[300,500,600,250,300,400,120,400,200];%ketebalan setiap lapisan

for i=1:norays
theta(i)=i*2; %sudut tembak sinar
end

for k=1:norays
for i=1:nolayers-1
theta(i+1,k)=(180/pi) * asin(sin(theta(i,k).*pi/180).*(vel(i+1)./vel(i))); %menghitung perubahan sudut sinar di setiap lapisan dengan menggunakan hukum snellius
end
end
for k=1:norays
p(k)=sin(theta(1,k).*pi/180)./vel(1); %menghitung ray parameter
end

for k=1:norays
for i=1:nolayers
dx(i,k)=(p(k)*vel(i).*dz(i))/sqrt(1-p(k)*p(k).*vel(i).*vel(i)); %menghitung jarak lateral di setiap lapisan
dt(i,k)=dz(i)/(vel(i).*sqrt(1-p(k)*p(k).*vel(i).*vel(i))); %menghitung waktu tempuh di setiap lapisan
end
end

for k=1:norays
twt(k)=2*sum(dt(:,k)); %menghitung twt untuk masing-masing sinar
end

%%%memanipulasi offset
dx_down=dx;
dx_up=flipud(dx_down);
dx=[dx_down;dx_up];
dx(1,1)=dx(1,1);
for k=1:norays
for i=2:nolayers*2,
dx(i,k)=dx(i-1,k)+dx(i,k);
end
end
nol=[1:norays]*0;
dx=[nol;dx];

%%memanipukasi kedalaman
dz=dz(1:nolayers);
dz(1)=dz(1);
for i=2:nolayers,
dz(i)=dz(i-1)+dz(i);
end
dz_down=dz';
dz_up=flipud(dz_down);
dz_up=dz_up(2:nolayers);
dz=[0;dz_down;dz_up;0];
offset=dx(nolayers*2+1,:);
% plot hasil
for k=1:norays
subplot(1,2,1)
plot(dx(:,k),dz); hold on
end
xlabel('offset(m)')
ylabel('depth(m)')
title('Jejak Sinar')
state=set(gca,'ydir');
if (strcmp(state,'reverse'))
set(gca,'ydir','reverse')
else
set(gca,'ydir','reverse')
end

a=size(dx);
dx=reshape(dx,a(1,1)*a(1,2),1);
x = [0 max(dx)];
for i=1:nolayers
y = [dz(i) dz(i)];
plot(x,y,'r'); hold on
end
axis([0 max(dx) 0 max(dz)]);
subplot(1,2,2)
plot(offset,twt,'linewidth',3); grid on; hold on
xlabel('offset(m)')
ylabel('twt(s)')
title('Kurva Waktu Tempuh')
state=set(gca,'ydir');
if (strcmp(state,'reverse'))
set(gca,'ydir','reverse')
else
set(gca,'ydir','reverse')
end
clear
end