Seismic Tomography

Seismic Tomography adalah teknik pencitraan untuk mendeterminasi variasi sifat fisis batuan seperti kecepatan gelombang P, gelombang S, atenuasi, dll.

Terdapat beberapa teknik tomography diantaranya refraction traveltime tomography, finite-frequency traveltime tomography, reflection traveltime tomography dan waveform tomography.

Seismic tomography biasanya dirumuskan sebagai sebuah inverse problem. Dalam refraction traveltime tomography, data yang diobservasi adalah first-arrival traveltimes t dan parameter model slowness s.

Hubungan antara t dan s dirumuskan sbb:

t = Ls

dimana L adalah kernel matrix berupa jarak jejak sinar dalam sistem sel tomography.

Selain inverse problem, forward problem merupakan hal yang sangat penting dalam kesuksesan tomography. Forward problem bertujuan untuk menelusuri jejak gelombang (ray tracing) dari sumber ke penerima di dalam medium. Medium dalam wilayah study tomography selanjutnya didiskritisasi menjadi sistem sel tomography.

Gambar di bawah ini merupakan contoh ray tracing dengan menggunakan teknik Bresenham. Dalam hal ini, diasumsikan jejak sinar dari sumber ke penerima berupa garis lurus.

Gambar di bawah ini adalah contoh hasil tomography untuk model checkerboard dengan konfigurasi sumber penerima di atas.

Model Checkerboard

Hasil rekonstruksi Tomography untuk model checkerboard di atas

Berikut adalah kode Matlab untuk Tomography di atas:
%%%Traveltime Tomography Code
%%%Developed by Agus Abdullah, PhD
%%%Prepared for ensiklopediseismik.blogspot.com
%%%You are allowed to use and copy this as long as you mention i source: 'ensiklopediseismik'.

clear; clc
a1=[0:9]*100;
a2=a1*0;
A=[a1,a2]';
B=[a2,a1]';

c1=a2+900;
C=2*[a1,c1]';
D=2*[c1,a1]';

for i=1:20
for j=1:20
coord{i,j}=brlinexya(B(i),A(i),C(j),D(j));
end
end

coord2=reshape(coord,[],1);
coord3=cell2mat(coord2);
figure
plot(coord3(:,1), coord3(:,2));


vel=1000*(1+repmat(eye([2,2]),5,5));
vel=vel(1:9,1:9);
vel=reshape(vel,[],1);
save vel.txt vel -ascii

h = waitbar(0,'ray tracing...');
for i = 1:400 ;
waitbar(i/400)
lon=coord2{i}(:,1);
lat=coord2{i}(:,2);
jejak=[lon,lat];
Listpath = dir('*.path');
Npath = length(Listpath);
newpath=Npath+1;
filepath = sprintf('jejak%i.path', newpath);
save(filepath,'jejak','-ASCII')
end
close(h)

% clear; clc
aa1=[0:200:1600]';
aa2=[200:200:1800]';
aa=[aa1,aa2];
latminmax=repmat(aa,9,1);
latMIN=latminmax(:,1);
latMAX=latminmax(:,2);

for i=1:9
bb1a{i}=(0*[1:9]')+(i-1)*200;
bb2a{i}=(200+(0*[1:9]'))+(i-1)*200;
end
bb=[bb1a',bb2a'];
lonminmax=cell2mat(bb);
lonMIN=lonminmax(:,1);
lonMAX=lonminmax(:,2);

latlon=[latminmax,lonminmax];
save latlon.txt latlon -ascii

clear; clc
h = waitbar(0,'menghitung matrix kernel...');
for i = 1:400;
waitbar(i/400)
result=load(sprintf('jejak%i.path',i));
lon=result(:,1);
lat=result(:,2);
cellcorrdinates=load('latlon.txt');
latMIN=cellcorrdinates(:,1);
latMAX=cellcorrdinates(:,2);
lonMIN=cellcorrdinates(:,3);
lonMAX=cellcorrdinates(:,4);

for k=1:length(lat),
for i=1:length(latMIN),
ifz ((lat(k,1)'>'=latMIN(i,1)) & (lat(k,1)'<'latMAX(i,1)) & (lon(k,1)'>'=lonMIN(i,1)) & (lon(k,1)'<'lonMAX(i,1)) );
dlatlondepth(k,i)=1;
end;
end
end

[a,b]=find(dlatlondepth==1);

for i=1:(length(b)-1)
c(i)=b(i+1)-b(i);
end

[d,e]=find(c~=0);
f=[0,e,length(b)];

for i=1:(length(f)-1)
index{i}=a(f(1,i)+1:f(1,i+1)); %%%change a to b
end


for i=1:(length(f)-1)
panj(i)=length(index{i});
end

for i=1:(length(f)-1)
for k=1:panj(i),
latcell3(i,k)=lat(index{i}(k,1));
loncell3(i,k)=lon(index{i}(k,1));
end
end


for i=1:(length(f)-1)
latcell4{i}=latcell3(i,:);
loncell4{i}=loncell3(i,:);
end

for i=1:(length(f)-1)
latcell2{i}=latcell4{i}(1:panj(i));
loncell2{i}=loncell4{i}(1:panj(i));
end


for i=1:(length(f)-1)
if (index{i}(1,1))==1,
indexadd(i)=index{i}(1,1);
else
indexadd(i)=index{i}(1,1)-1;
end
end

for i=1:(length(f)-1)
latadd (i)=lat(indexadd(1,i));
lonadd (i)=lon(indexadd(1,i));
end

for i=1:(length(f)-1)
latcell{i}=[latadd(i),latcell2{i}];
loncell{i}=[lonadd(i),loncell2{i}];
end


for i=1:(length(f)-1)
for k=1:length(latcell{i})-1
latcellSIG(i,k)=(latcell{i}(1,k+1))-(latcell{i}(1,k));
loncellSIG(i,k)=(loncell{i}(1,k+1))-(loncell{i}(1,k));
dist2(i,k)=((latcellSIG(i,k)).^2)+ ( (loncellSIG(i,k)).^2) ;
dist1(i,k)=(dist2(i,k)).^0.5;
end
end

for i=1:(length(f)-1),
dist{i}=dist1(i,:);
dl(i)=sum(dist{i});
end

m=[e,length(b)];
for i=1:length(m),
cell(i)=b(m(1,i));
end
result1=[cell',dl'];

Listpath = dir('*.sell');
Npath = length(Listpath);
newpath=Npath+1;
filepath = sprintf('jejak%i.sell', newpath);
save(filepath,'result1','-ASCII')
clearex('h')
end %end of data
close(h)


clear
ListDir = dir('*.sell');
NFiles = length(ListDir);

for i = 1:NFiles ;
result{i}=load(sprintf('jejak%i.sell',i),'%f');
end

for i = 1:NFiles ;
sel{i}=result{i}(:,1);
dl{i}=result{i}(:,2);
end


clear result
for i=1:NFiles ,
LMAT{i}=[1:81].*0;
LMAT{i}(sel{i})=dl{i};
end

clear sel dl i
G1=cell2mat(LMAT');
clear LMAT
G3=reshape(G1,NFiles ,81);
clear G1

G=sparse(G3);

save GP3D.mat G
clear; clc

load('GP3D.mat'); %load kernel matrix
load vel.txt
slowness_syn=1./vel;
time_syn=G*slowness_syn;
Cd=(1/(0.025*0.025))*sparse(eye(400));
Cm=(1/(0.0001))*sparse(eye(81));
A=G'*Cd*G+0.5*Cm;
b=G'*Cd*time_syn;
slow_inv=lsqr(A,b,0.0001,300);
[X,Y] = meshgrid(0:200:1600);
slow_inv=flipud(slow_inv);
slow_inv=reshape(slow_inv,[],9);

vel_inv=1./slow_inv;
vel=reshape(vel,[],9);
colormap(gray)
ni=3;ri=25; ai=65 ;bi=65;
vel_inv= interp2(vel_inv,ni,'bicubic'); %%vel orig
vel= interp2(vel,ni,'bicubic'); %%vel orig
x1=meshgrid(0: ri:1600,1:bi);
lonA=x1';
latA=meshgrid(0:ri:1600,1:ai);
surface(lonA,latA,vel);
colormap(jet)
shading interp
caxis([1000 2000]);
figure
surface(lonA,latA,vel_inv);
colormap(jet)
shading interp
caxis([1000 2000]);
%%%%%%%%%save function berikut pada .m terpisah
function clearex(varargin)
%% This function clear all workspace
% except for one or more selected variable
%
% File created on Feb 15th 2008
%
% Last modified: Feb 15th 2008
%
% Author: Arnaud Laurent
%
% Inputs: name of variables to keep (e.g. 'a','b','c')
% Note: It is possible to use wildcard (e.g. 'a*')

a = evalin('base','who');

var = cell(size(varargin));

for i=1:nargin
var{i}=varargin{i};
end

assignin('base','ClEaRGsJioU',var);
var = evalin('base','who(ClEaRGsJioU{:})');

clearvar = a(~ismember(a,var));
assignin('base','ClEaRGsJioU',clearvar);

evalin('base','clear(ClEaRGsJioU{:},''ClEaRGsJioU'')')

%%%%%%%%%save function berikut pada .m terpisah
function [Coords]=brlinexya(Sx,Sy,Ex,Ey)
% function [Coords]=brlinexya(Sx,Sy,Ex,Ey)
% Bresenham line algorithm.
% Sx, Sy, Ex, Ey - desired endpoints
% Coords - nx2 ordered list of x,y coords.
% Author: Andrew Diamond;
%
% if(length(M) == 0)
% M = zeros(max([Sx,Sy]),max([Ex,Ey]));
% end
Dx = Ex - Sx;
Dy = Ey - Sy;
% Coords = [];
CoordsX = zeros(2 .* ceil(abs(Dx)+abs(Dy)),1);
CoordsY = zeros(2 .* ceil(abs(Dx)+abs(Dy)),1);
iCoords=0;
if(abs(Dy) <= abs(Dx))
if(Ey >= Sy)
if(Ex >= Sx)
D = 2*Dy - Dx;
IncH = 2*Dy;
IncD = 2*(Dy - Dx);
X = Sx;
Y = Sy;
% M(Y,X) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sx;
CoordsY(iCoords) = Sy;
while(X < Ex)
if(D <= 0)
D = D + IncH;
X = X + 1;
else
D = D + IncD;
X = X + 1;
Y = Y + 1;
end
% M(Y,X) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = X;
CoordsY(iCoords) = Y;
% Coords = [Coords; [X,Y]];
end
else % Ex < Sx
D = -2*Dy - Dx;
IncH = -2*Dy;
IncD = 2*(-Dy - Dx);
X = Sx;
Y = Sy;
% M(Y,X) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sx;
CoordsY(iCoords) = Sy;
while(X > Ex)
if(D >= 0)
D = D + IncH;
X = X - 1;
else
D = D + IncD;
X = X - 1;
Y = Y + 1;
end
% M(Y,X) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = X;
CoordsY(iCoords) = Y;
% Coords = [Coords; [X,Y]];
end
end
else % Ey < Sy
if(Ex >= Sx)
D = 2*Dy + Dx;
IncH = 2*Dy;
IncD = 2*(Dy + Dx);
X = Sx;
Y = Sy;
% M(Y,X) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sx;
CoordsY(iCoords) = Sy;
while(X < Ex)
if(D >= 0)
D = D + IncH;
X = X + 1;
else
D = D + IncD;
X = X + 1;
Y = Y - 1;
end
% M(Y,X) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = X;
CoordsY(iCoords) = Y;
% Coords = [Coords; [X,Y]];
end
else % Ex < Sx
D = -2*Dy + Dx;
IncH = -2*Dy;
IncD = 2*(-Dy + Dx);
X = Sx;
Y = Sy;
% M(Y,X) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sx;
CoordsY(iCoords) = Sy;
while(X > Ex)
if(D <= 0)
D = D + IncH;
X = X - 1;
else
D = D + IncD;
X = X - 1;
Y = Y - 1;
end
% M(Y,X) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = X;
CoordsY(iCoords) = Y;
% Coords = [Coords; [X,Y]];
end
end
end
else % abs(Dy) > abs(Dx)
Tmp = Ex;
Ex = Ey;
Ey = Tmp;
Tmp = Sx;
Sx = Sy;
Sy = Tmp;
Dx = Ex - Sx;
Dy = Ey - Sy;
if(Ey >= Sy)
if(Ex >= Sx)
D = 2*Dy - Dx;
IncH = 2*Dy;
IncD = 2*(Dy - Dx);
X = Sx;
Y = Sy;
% M(X,Y) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sy;
CoordsY(iCoords) = Sx;
while(X < Ex)
if(D <= 0)
D = D + IncH;
X = X + 1;
else
D = D + IncD;
X = X + 1;
Y = Y + 1;
end
% M(X,Y) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = Y;
CoordsY(iCoords) = X;
% Coords = [Coords; [Y,X]];
end
else % Ex < Sx
D = -2*Dy - Dx;
IncH = -2*Dy;
IncD = 2*(-Dy - Dx);
X = Sx;
Y = Sy;
% M(X,Y) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sy;
CoordsY(iCoords) = Sx;
while(X > Ex)
if(D >= 0)
D = D + IncH;
X = X - 1;
else
D = D + IncD;
X = X - 1;
Y = Y + 1;
end
% M(X,Y) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = Y;
CoordsY(iCoords) = X;
% Coords = [Coords; [Y,X]];
end
end
else % Ey < Sy
if(Ex >= Sx)
D = 2*Dy + Dx;
IncH = 2*Dy;
IncD = 2*(Dy + Dx);
X = Sx;
Y = Sy;
% M(X,Y) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sy;
CoordsY(iCoords) = Sx;
while(X < Ex)
if(D >= 0)
D = D + IncH;
X = X + 1;
else
D = D + IncD;
X = X + 1;
Y = Y - 1;
end
% M(X,Y) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = Y;
CoordsY(iCoords) = X;
% Coords = [Coords; [Y,X]];
end
else % Ex < Sx
D = -2*Dy + Dx;
IncH = -2*Dy;
IncD = 2*(-Dy + Dx);
X = Sx;
Y = Sy;
% M(X,Y) = Value;
% Coords = [Sx,Sy];
iCoords = iCoords + 1;
CoordsX(iCoords) = Sy;
CoordsY(iCoords) = Sx;
while(X > Ex)
if(D <= 0)
D = D + IncH;
X = X - 1;
else
D = D + IncD;
X = X - 1;
Y = Y - 1;
end
% M(X,Y) = Value;
iCoords = iCoords + 1;
CoordsX(iCoords) = Y;
CoordsY(iCoords) = X;
% Coords = [Coords; [Y,X]];
end
end
end
end
Coords = [CoordsX(1:iCoords),CoordsY(1:iCoords)];