LAMPIRAN I
Lampiran I. Kode RelocationGeigerMethod.m, Kode RayTracingShootingMethod.m, Kode MatrixJacobi.m, Kode EventGather.m, Format penulisan Kecepatan, Waktu tiba gelombang, Stasiun, Hasil pengolahan dengan GAD dan script program Mathlab. 1. Kode RelocationGeigerMethod.m clear,clc %% Simple Geiger Method %% Load DataCatalog [FileName,PathName] = uigetfile(‘*.txt’,'Select the data catalog’); Data = load([PathName FileName]); %% Check Event Gather [EventID EventHypo First Last] = EventGather(Data); %% Velocity Model [FileName,PathName] = uigetfile(‘*.txt’,'Velocity 1-D Model’); Velocity = load([PathName FileName]); VelocityLayer = Velocity(:,2); TopOfLayer = Velocity(:,1); VelocityModel = struct(‘VelocityLayer’,VelocityLayer,’TopOfLayer’,TopOfLayer’); %% Calculate the Time Deviation of Initial Hypocenter for i=1:size(Data,1) Sumber = Data(i,1:3); Penerima = Data(i,4:6); [RayPath TravelTime] = RayTracingShootingMethod(Sumber, Penerima, VelocityModel); DtAwal(i) = TravelTime-Data(i,7); end NewData = Data; NewData =[ NewData zeros(size(Data,1),1)]; WB = waitbar(0,’Please wait…’); for i =1:EventID waitbar(i/EventID) InitialHypo(i,:) = EventHypo(i,:); nReceiver = Last(i) – First(i) + 1; RelativeRms = 10; ite = 0; while RelativeRms > 0.00001 ite = ite+1; MatrixAllJacobi = []; for n = 1:nReceiver
if ite ==1 Sumber = InitialHypo(i,:); ObservationTime(n) = Data((i-1)*nReceiver+n,7); else Sumber = Update; ObservationTime(n) = Data((i-1)*nReceiver+n,7)+DeltaModel(4); end Penerima = Data((i-1)*nReceiver+n,4:6); %Ray Tracing using shooting method based on Snell’s Law, only %direct wave [RayPath TravelTime] = RayTracingShootingMethod(Sumber, Penerima, VelocityModel); NewData((i-1)*nReceiver+n,8) = TravelTime; DelayTime(n,:) = TravelTime-ObservationTime(n); DelayTimeAll((i-1)*nReceiver+n,:) = TravelTime-ObservationTime(n); %Build Matrix Jacobi [Jacobi(n,:)] = MatrixJacobi(RayPath,VelocityModel); end %Add Damping Alpha = 0.00001; DampedJacobi = [Jacobi; Alpha*eye(size(Jacobi,2),size(Jacobi,2))]; DTzeros = [DelayTime; zeros(size(Jacobi,2),1)]; %Inversion Solver DeltaModel=lsqr(DampedJacobi,DTzeros); %Update Model (Hypo & Origin Time) Update = Sumber+DeltaModel(1:3)’; ErrorRms(ite) = sqrt(mean(DelayTime.^2)); if ite >1 RelativeRms = ErrorRms(ite)-ErrorRms(ite-1); end for m=1:nReceiver NewData((i-1)*nReceiver+m,1:3) = Update; NewData((i-1)*nReceiver+m,7) = NewData((i-1)*nReceiver+m,7)+DeltaModel(4); end end UpdatedHypo(i,:) = Update; end close(WB)
figure(1) plot3(EventHypo(:,1),EventHypo(:,2),EventHypo(:,3),’.r’) hold on plot3(UpdatedHypo(:,1),UpdatedHypo(:,2),UpdatedHypo(:,3),’.k’) hold off set(gca,’Zdir’,'Reverse’) grid on figure(2) subplot(211) hist(DtAwal,-2:0.2:2) subplot(212) hist(DelayTimeAll,-2:0.2:2) % NewData = [NewData DelayTimeAll]; % [filename, pathname] = uiputfile(‘*.txt’, ‘Save as data catalog’); % save(filename,’NewData’,'-ASCII’)
2. Kode RayTracingShootingMethod.m function [RayPath TravelTime] = RayTracingShootingMethod (Sumber, Penerima, VelocityModel) Source = struct(‘X’,Sumber(1),’Y',Sumber(2),’Z',Sumber(3)); Receiver = struct(‘X’,Penerima(1),’Y',Penerima(2),’Z',Penerima(3)); epicenter = sqrt((Source.X-Receiver.X)^2+(Source.Y-Receiver.Y)^2); %% NumberLayer [NumberLayerOfSource] = CheckLayer(Source.Z,VelocityModel.TopOfLayer); [NumberLayerOfReceiver] = CheckLayer(Receiver.Z,VelocityModel.TopOfLayer); %% ShootingMethod %InitialAngle ArrivalAngle = [0.000000000001 89.9999999999999 (89.99999-0.00001)/2]; Error = 10; while Error>0.001 ArrivalAngle(3) = (ArrivalAngle(2)+ArrivalAngle(1))/2; for i = 1:length(ArrivalAngle) %Start from Lower Layer RayTracingX = Source.X; RayTracingY = Source.Y; RayTracingZ = Source.Z; if NumberLayerOfSource == NumberLayerOfReceiver RayTracingX = [RayTracingX Receiver.X]; RayTracingY = [RayTracingY Receiver.Y];
RayTracingZ = [RayTracingZ Receiver.Z]; DeltaError = zeros(1,3); else for nLay = NumberLayerOfSource:-1:NumberLayerOfReceiver if nLay ~= NumberLayerOfReceiver if nLay == NumberLayerOfSource deltaZ = Source.Z – VelocityModel.TopOfLayer(nLay); deltaRI = tan(deg2rad(ArrivalAngle(i)))*deltaZ; TransmissionAngle = asin(sin(deg2rad(ArrivalAngle(i)))*VelocityModel.VelocityLayer(nLay1)/VelocityModel.VelocityLayer(nLay)); RayX = Source.X; RayY = Source.Y; else deltaZ = VelocityModel.TopOfLayer(nLay+1) – VelocityModel.TopOfLayer(nLay); deltaRI = tan(TransmissionAngle)*deltaZ; TransmissionAngle = asin(sin(TransmissionAngle)*VelocityModel.VelocityLayer(nLay1)/VelocityModel.VelocityLayer(nLay)); end RayZ = VelocityModel.TopOfLayer(nLay); else deltaZ = VelocityModel.TopOfLayer(NumberLayerOfReceiver+1) – Receiver.Z; deltaRI = tan(TransmissionAngle)*deltaZ; RayZ = Receiver.Z; end deltaY = (abs(Source.Y-Receiver.Y)/epicenter)*deltaRI; deltaX = sqrt(deltaRI^2-deltaY^2); if Source.X<=Receiver.X RayX = RayX + deltaX; else RayX = RayX – deltaX ; end if Source.Y<=Receiver.Y RayY = RayY + deltaY; else RayY = RayY – deltaY; end RayTracingX = [RayTracingX RayX]; RayTracingY = [RayTracingY RayY]; RayTracingZ = [RayTracingZ RayZ]; end
ErrorDistance(i) = sqrt(((RayTracingX(length(RayTracingX))Source.X))^2+((RayTracingY(length(RayTracingY))-Source.Y))^2); DeltaError(i) = epicenter-ErrorDistance(i); Error = abs(DeltaError(i)); end end if (DeltaError(1) * DeltaError(3) < 0) ArrivalAngle(2) = ArrivalAngle(3); elseif (DeltaError(1) * DeltaError(3) > 0) ArrivalAngle(1) = ArrivalAngle(3); else break; end end %plot3(RayTracingX,RayTracingY,RayTracingZ,’-k’) RayPath = struct(‘X’,RayTracingX,’Y',RayTracingY,’Z',RayTracingZ); TravelTime = TravelTimeCalculation (RayPath, VelocityModel); function TravelTime = TravelTimeCalculation (RayPath, VelocityModel) TravelTime =0; for i=1:length(RayPath.X)-1 MidPoint(i) = (RayPath.Z(i)+RayPath.Z(i+1))/2; [NumberLayer] = CheckLayer(MidPoint(i),VelocityModel.TopOfLayer); LengthSegment(i) = sqrt(((RayPath.X(i+1)-RayPath.X(i))^2)+((RayPath.Y(i+1)RayPath.Y(i))^2)+((RayPath.Z(i+1)-RayPath.Z(i))^2)); TravelTime = TravelTime + LengthSegment(i)/VelocityModel.VelocityLayer(NumberLayer); end function [NumberOfLayer] = CheckLayer(Depth,TopOfLayer) NumberOfLayer = 1; for j = 1:length(TopOfLayer)-1 if Depth >= TopOfLayer(j) && Depth < TopOfLayer(j+1) NumberOfLayer = j; end if Depth == TopOfLayer(j) && j~=1 NumberOfLayer = j – 1 ; end end if Depth > TopOfLayer(length(TopOfLayer)) NumberOfLayer = length(TopOfLayer); end
3. Kode MatrixJacobi.m
function [RowOfJacobi] = MatrixJacobi(RayPath,VelocityModel) dtdx=0; dtdy=0; dtdz=0; for i =1:length(RayPath.X)-1 NumberLayer = CheckLayer((RayPath.Z(i+1)+RayPath.Z(i))/2,VelocityModel.TopOfLayer); Divider = VelocityModel.VelocityLayer(NumberLayer)*sqrt((RayPath.X(i+1)RayPath.X(i))^2+(RayPath.Y(i+1)-RayPath.Y(i))^2+(RayPath.Z(i+1)RayPath.Z(i))^2); dtdx = dtdx + ((RayPath.X(i+1)-RayPath.X(i))/Divider); dtdy = dtdy + ((RayPath.Y(i+1)-RayPath.Y(i))/Divider); dtdz = dtdz + ((RayPath.Z(i+1)-RayPath.Z(i))/Divider); end RowOfJacobi = [dtdx dtdy dtdz 1];
4. Kode EventGather.m function [EventID EventHypo First Last] = EventGather(Data) First = 1; Last = []; EventHypo =[]; for i=1:size(Data,1)-1 if Data(i,1)~= Data(i+1,1) First =[ First i+1]; Last =[Last i]; EventHypo = [EventHypo; Data(i,1:3)]; end end Last = [Last size(Data,1)]; EventHypo = [EventHypo; Data(size(Data,1),1:3)]; EventID = length(Last);
5. Kecepatan (Velocity.dat) a. Format penulisan file Velocity.dat, yaitu sebagai berikut: Baris 1: jumlah layer Baris 2: ketidak seragaman koordinat Z dalam kasus jumlah layer= 1. Baris ini dilewati oleh program. Baris 3: nilai kecepatan gelompang P, dituliskan dalam format f8.3 untuk setiap layer. Baris 4: nilai kecepatan gelombang S, dituliskan dalam format f8.3 untuk setiap layer. b. Data file velocity.dat yang digunakan dalam perhitungan hiposenter: nzLayer : 2 (15x,i4) zLayer : 0.0 (15x,5f8.3) VpL(nzLayer+1): 2.8 3.9 (15x,6f8.3) VsL(nzLayer+1): 1.556 2.167 (15x,6f8.3)
6. Waktu tiba (Arrival.dat) a. Format penulisan file Arrival.dat, yaitu sebagai berikut: Kolom 1-10: YYMMDDHHmm Kolom 11: spasi atau “,” Kolom 12-14: kode stasiun (3 karakter) Kolom 15: spasi atau “,” Kolom 16-21: waktu tiba gelombang P dalam format f6.3. Jika tidak ada data maka dituliskan “99.990” Kolom 22: spasi atau “,” Kolom 23: polarisasi gelombang P, dituliskan “+” atau “-“ atau spasi. Kolom 24: spasi atau “,” Kolom 25: jika awalan dari gelombang P jelas, dituliskan “I”, sebaliknya jika tidak jelas dituliskan “E”. Kolom 26: spasi atau “,” Kolom 27-32: waktu tiba gelombang S, ditulis dalam format f6.3. Jika tidak ada data maka dituliskan “99.990”. Kolom 33: spasi atau “,” Kolom 34: jika awalan dari gelombang S jelas, dituliskan “I”, sebaliknya jika tidak jelas dituliskan “E”. b. Contoh data file Arrival.dat yang digunakan dalam perhitungan hiposenter gempa Gunung Lokon, yaitu sebagai berikut: 1204011440 EMP 22.701 + I 24.805 I 1204011440 KIN 22.904 + I 25.611 I 1204011440 WLN 23.367 + E 25.536 I 1204011440 TTW 22.972 + E 25.632 I 7. Stasiun (Station.dat) a. Format penulisan file Station.dat, yaitu sebagai berikut: Baris 1: Jumlah stasiun Baris 2: Kode stasiun (3 karakter) dengan koordinat stasiun (X,Y,Z) dalam meter. b. Data yang digunakan dalam file Station.dat, yaitu sebagai berikut: 5 (i2) EMP 0.102 0.398 -1.143 (a3,2x,3f10.3) KIN 1.907 0.253 -0.914 WLN 0.306 -1.668 -1.072 SEA 0.049 0.620 -1.162 MHW 6.693 -1.412 -1.410
8. Hasil (Results.dat) Contoh tampilan hasil pengolahan GAD gempa vulkanik Gunung Lokon. nst : 5 Station List EMP .102 .398 -1.143 KIN 1.907 .253 -.914 WLN .306 -1.668 -1.072 SEA .049 .620 -1.162 MHW 6.693 -1.412 -1.410 nZLayer: 2 zLayer : .000 Vp : 2.800 3.900 Vs : 1.556 2.167 Hypocenter Date 12 5 6 Time 8:56 Focal Element Probable Error X -.372 3.614 Y 5.643 2.564 Z -.003 2.632 T 36.157 .251 Travel time residual rms= .144sec. ST P S Cal (Obs-Cal) EMP 37.580 37.792 -.212 WLN 38.530 38.306 .224 SEA 37.727 37.739 -.012 EMP 39.013 39.099 -.086 WLN 40.013 40.026 -.013 SEA 39.103 39.004 .099 9. Script program hiposenter Mathlab Script program hiposenter 3D Mathlab, yaitu sebagai berikut: close all;clear all;clc; data=load('kontur.txt'); x=data(:,1);y=data(:,2);z=data(:,3); dx=0.1; dy=0.1; xi=-10:dx:10; yi=-10:dy:10; [XI,YI]=meshgrid(xi,yi); ZI=griddata(x,y,z,XI,YI); mesh(XI,YI,ZI) hold on axis([-6 6 -6 6])
hipo=load('hipo.txt'); hipo(:,3)=abs(hipo(:,3)).*-2; plot3(hipo(:,1),hipo(:,2),hipo(:,3),'ok','MarkerFaceColor',' r') daspect([0.5 0.5 0.4])
