**brief theory**

Earthquake hypocenter location is an important task in seismology. In the past seismologists had to use geometrical methods (e.g. Circle methods) to find it. It wasn't until the proposal of inversion theory that a new perspective of solving this problem was implemented. The core of this technique assumes that given a set of data we can obtain the model that those data represents based on

d=Gm

Where d is a matrix containing the data with size(nx1)

G is a matrix containing the change in the model with size (nxm)

m is the model matrix with size(mx1)

If G is square, which implies we have a system of n=m equations and n=m unknowns, we could simply multiply both sides of the equation by G ^{-1} the inverse of G

So we would have

G ^{-1} d=m

However, this is not the case in seismology where the number of unknowns is huge. Therefore the matrices are not squares.So, this problem has to be solve by using inversion techniques.Fortunately, nowadays there are several tools available too perform inversions. Here we shall focus on how to find an eathquake hypocenter using matlab.

**Data**

Let's say there are 10 stations to which data(seismograms), you have access

Here is the data for station coordinates and arrival times arrival_times.txt

**Column 1: ** Station name

** Column 2: ** longitude (km)

** Column 3: ** latitude (km)

** Column 4: ** height (m) at station site

** Column 5: ** P wave arrival time (s)

** Column 6: ** S wave arrival time (s)

** Column 7: ** S-P wave arrival time (s)

** Column 8: ** P wave velocity (km/s)

** Column 8: ** S wave velocity (km/s)

Let’s import these files in Matlab.

%assigning values

A = importdata('arrival_times.txt', '\t', 1);

load val.txt

tp=A.data(:,4);

ts=A.data(:,5);

v2=[A.data(:,7);A.data(:,8)];

name = A.textdata(:,1);

xco = [A.data(:,1);A.data(:,1)];

yco = [A.data(:,2);A.data(:,2)];

zco = [A.data(:,3)/1000;A.data(:,3)/1000];

Declare some variables

%Definition of data vector matrix, Kernel matrix, and others

G=zeros(length(xco),4);

dobs=zeros(length(xco),1);

dtdx=zeros(length(xco),1);

dtdy=zeros(length(xco),1);

dtdz=zeros(length(xco),1);

dtdt=zeros(length(xco),1);

arrvt=zeros(2*length(xco),1);

dpred=zeros(length(xco),1);

d=zeros(length(xco),1);

res=zeros(length(xco),1);

mic=zeros(100000,4);

dmn=zeros(4,1);

mestn=zeros(4,1);

OT=187.2;

xg=val(1,1);

yg=val(1,2);

zg=10;

tg=OT;

mi=[xg;yg;zg;tg];

cont=0;

arrvt=[tp;ts];

xg, yg, and zg, corresponds to the intial guess which is a value that(for the first calculation) one have to propose, in this example this value is on the file val.txt

Let’s plot the 10 seismic stations.

% Plotting stations locations

figure(1);

plot(xco(:,1),yco(:,1),'b*')

hold on

Figure 1. Seismic stations location.

Let’s plot our initial guess.

% Plotting the initial guess as a green circle

figure(2);

plot(xco(:,1),yco(:,1),'b*')

hold on

plot(xg,yg,'go')

Figure 2. Green circle represents our initial guess for the earthquake hypocenter location.

**Algorithm to find the solution**

As mentioned before, To solve this problem we have to use the generalized inverse theory, so we need to construct three matrices

First of all the data matrix d=[t1,t2,t3..tn] where t1,t2,t3...tn are arrival times in station 1, station 2 and so on respectively.

then the kernel matrix

G=[dti/dx,dti/dy,dti/dz,dti/dt] where dti/dx,dti/dy,dti/dz,dti/dt are partial derivatives with respect to time.

Finally the model matrix m=[mx,my,mz,mt] where mx,my,mz,mt are for the first iteration the initial guess values this is

mx=xg

my=yg

mz=zg

mt=tg

and after the first iteration the model matrix must be updated with the inverted one each time.

Let’s see the code

% definition of the intial minimum error

E=10;

%begin the iterations

while min(abs(E))>(1*10^0);

cont=cont+1;

if cont<=1

dobs(:)=[tp(:);ts(:)];

for i=1:18

alpha=v2(i);

% calculating the partial derivatives

dtdx(i)=(-(xco(i)-mi(1)));

dtdx(i)=dtdx(i)/(alpha^2*(arrvt(i)-mi(4)));

dtdy(i)=(-(yco(i)-mi(2)));

dtdy(i)=dtdy(i)/(alpha^2*(arrvt(i)-mi(4)));

dtdz(i)=(-(zco(i)-mi(3)));

dtdz(i)=dtdz(i)/(alpha^2*(arrvt(i)-mi(4)));

dtdt(i)=1;

% travel time calculation

ttt(i)=sqrt(((xco(i)-xg)*(xco(i)-xg))+((yco(i)-yg)*(yco(i)-yg))+((zco(i)-zg)*(zco(i)-zg)))/alpha;

dpred(i)=OT+ttt(i);

d(i)=dobs(i)-dpred(i);

end

G=[dtdx dtdy dtdz dtdt];

% Getting Eigenvalues and Eigenvectors of G

[U,S,V] = svd(G);

% Least square solution

dm=inv(G'*G)*G'*d;

% updating the 'mi' vector

mi=mi+dm;

mic(cont,:)=(mi(:))';

dc=G*mi;

E=(sum(d(:)).^2);

% Plot first solution as a cyan circle

figure(3);

plot(xco(:,1),yco(:,1),'b*')

hold on

plot(xg,yg,'go')

hold on

plot(mi(1),mi(2),'co')

After the first iteration (cyan circle) you can notice how the earthquke position has changed in comparison to our intial guess (green circle)

Figure 3. Earthquake hypocenter after the first iteration.

As surely you have noticed the new model matrix has been changed, this assure the next steps in the iteration will have a new "guessed model", but this time it is originated from the inversion

Now, lets look at the second part of the solution, which is pretty much the same as the first part

% Creating a new figure for the next stage

% and plotting first solution as a cyan circle

figure(4);

plot3(mi(1),mi(2),-mi(3),'co')

hold on

% evalauting condition before starting next iterations

elseif cont>1

dobs(:)=[tp(:);ts(:)];

for i=1:18

alpha=v2(i);

% calculating the partial derivatives

dtdx(i)=(-(xco(i)-mi(1,1)));

dtdx(i)=dtdx(i)/(alpha^2*(arrvt(i)-mi(4)));

dtdy(i)=(-(yco(i)-mi(2,1)));

dtdy(i)=dtdy(i)/(alpha^2*(arrvt(i)-mi(4)));

dtdz(i)=(-(zco(i)-mi(3,1)));

dtdz(i)=dtdz(i)/(alpha^2*(arrvt(i)-mi(4)));

dtdt(i)=1;

% travel time calculation

ttt(i)=sqrt(((xco(i)-mi(1,1))*(xco(i)-mi(1,1))+((yco(i)-mi(2,1))*(yco(i)-mi(2,1)))+((zco(i)-mi(3,1))*(zco(i)-mi(3,1))))/alpha;

dpred(i)=mi(4)+ttt(i);

d(i)=dobs(i)-dpred(i);

end

G=[dtdx dtdy dtdz dtdt];

% Getting Eigenvalues and Eigenvectors of G

[U,S,V] = svd(G);

% Least square solution

dm=inv(G'*G)*G'*d;

% updating the 'mi' vector

mi=mi+dm;

mic(cont,:)=(mi(:))';

dc=G*mi;

E=(sum(d(:)).^2);

end

end

%deleting zeros rows in the mic matrix

av=mic;

av(~any(mic,2),:)=[];

mic=av;

clear av

Then,all we need is just to follow small details for graphing(a work around to get different colours in the same graph)

for i=2:cont-1

% Plot each iteratively calculated solution as a blue circle

plot3(mic(i,1),mic(i,2),-mic(i,3),'o','Color',[0 0 1]);

end

In order to know how many iterations took the routine to get over the minimum error constrain

sprintf(' # of iterations %d',cont)

Finally, plot the best solution as a red circle

hold on

plot3(mic(cont,1),mic(cont,2),-mic(cont,3),'o','Color',[1 0 0]);

xlabel('X'),ylabel('Y'), grid;

hold on

plot3(xco(:,:),yco(:,:),-zco(:,:),'*','Color',[0 0 1])

The final graph should look like

Figure 4. Earthquake hypocenter from inversion(Red circle).

Note the change in depth after each iteration(blue circles, Fig. 5). This has been a clear

example of how to use computational tools for solving real life problems.

Here is the script for this example hypocenter.txt after saving the file, just change the

".txt" extension to ".m" and it will work fine.

Finally, a graphical review of the calculations is shown below

Figure 5. Earthquake hypocenter from inversion.Note the change in depth after each iteration.

Red circle represents the best solution.