blinkdagger

In this tutorial I want to explain to you how you can use Matlab to solve some real world problems. Our topic today is finding the epicenter (starting location) of an earthquake. Don’t worry if you have little background knowledge about earthquakes. I’ll be giving you enough information to understand this example.

Setting up the Problem

Now let’s pretend you are a seismologist and have deployed 25 seismic stations in Switzerland. These seismic stations are equipped to measure and record the motions of the ground. Your helpful assistants have organized the coordinates of Switzerland boundary in border.xy and data for analysis in file loctim.txt.

Column 1: Latitude (in degrees)
Column 2: Longitude (in degrees)
Column 3: Elevation (in km)
Column 4: Earthquake motion arrival time (in seconds after reference time 16:30)

47.39500 7.23130 0.860 45.969
46.94960 8.24280 0.880 50.980
46.77470 6.96140 0.720 52.146
… data continues

Let’s import these files in Matlab.

% Import Data Files
loc = importdata('loctim.txt');
bound = importdata('border.xy');
lat = loc(:,1);     % latitude
lon = loc(:,2);     % longitude
ele = loc(:,3);     % elevation
tim = loc(:,4);     % time of arrival

Let’s also plot Switzerland and the 25 seismic stations.

% Plot Switzerland
plot(bound(:,1),bound(:,2))
hold on;
% Plot Stations
plot(lon,lat,'*')

You should get

Starting with an Initial Guess

Now to actually find the epicenter, we need to start with an initial guess. What would be a good educated guess? Let’s start by understanding what characteristics the real solution will exhibit.

The true epicenter should be located where the times it takes to travel to each station will be consistent with the ones in loctim.txt. I am going to guess that the epicenter will be in Switzerland and we’ll test to see if the calculated arrival times seem reasonable.

Let’s assume that we are not interested in the elevation of the epicenter and that the earthquake motion travels in Switzerland at the speed of 5.8 km/s.

So I have chosen the initial location R as latitude 46.9 and longitude 9. Let’s see how close it is to the true epicenter.

Let’s choose an arbitrary station G (I chose station 6) out of the 25 stations and calculate the time it takes from our initial guess R to the station G. We can employ the distance formula.

time = distance/velocity
t-t0 = sqrt((x-x0)^2 + (y-y0)^2)/v

In matlab we code

% Initial Guess
lat0 = 46.9; lon0 = 9;
% Plot initial guess and notate with red circle
plot(lon0,lat0,'ro')
% Plot station G and notate with green circle
plot(lon(6),lat(6),'go')
 
% Parameters to change degrees to kilometers
latkm = 111.19; lonkm = 75.82; vp = 5.8;
 
 
% Initial Calculation
t0 = sqrt( ((lon(6)-lon0)*lonkm).^2 + ((lat(6)-lat0)*latkm).^2 )/vp

The code will yield t0 = 7.5216 seconds, meaning it takes 7.5216s for the earthquake to travel from our initial guess R to station G. Now from our loctim.txt I look at the actual recorded arrival time for station G - 53.715 sec with reference time of 16:30. So this implies that the earthquake started at 16:30:46.2 (since 53.715 – 7.5216 = 46.1934 s).

Now if you repeat this process with a different station G, you’ll realize you get different earthquake start times! We know this is impossible since a single earthquake only has a single start time. But this doesn’t mean this is a bad initial guess. It just means that we do not have the true solution yet.

You could painstakingly repeat the process by comparing the values, or you can code an algorithm and make the computer seek the solution for you. The algorithm should minimize the error between the calculated and recorded arrival times. If the error is zero, we know the position is the true epicenter.

Algorithm to Find the Solution

How do we code such an algorithm? It is definitely not an easy task if you are not familiar with linear algebra. But with some high school mathematics background, we can reason out an algorithm.

We know that as we are closer to the true epicenter, our error will also decrease. You can also think of it as we are closer to the true solution, smaller rate of change to the coordinates R is needed to pin down the actual location. Evidently we can represent this through setting the derivatives to zero.

Take the derivatives of aforementioned distance formula with respect to x, y, and t. Or you can just use the following code.

% Initial Guess Stored in Array m
m(1) = t0; m(2)= lon0; m(3) = lat0;
 
% Loop through All Stations
for i=1:length(lat)
        diffx = (lon(i)-m(2))*lonkm;
        diffy = (lat(i)-m(3))*latkm;
        D0(i) = sqrt(diffx^2+diffy^2);          % Distance
        d(i) = tim(i) - D0(i)/vp - m(1);         % Time
        % Store derivative in a matrix G
        G(i,3) = -(diffy/D0(i))/vp; G(i,2) = -(diffx/D0(i))/vp; G(i,1) = 1;
end

Now I have stored the initial guess in array m, error in time of arrival in d, and the derivatives in G. I also have calculated the values for all stations by putting them in a for loop. So for each station, I can write a single corresponding equation. But instead of listing each of the equations out, I can utilize our variables m, d, and G and write

G*m = d

Consider this as a shorthand for displaying system of equations.

We know that if we improve m (originally stores our initial guess) correctly, we will be closer to our true solution. For the sake of this tutorial let’s assume that we need to update m with the following code

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

Just treat this as a black box for now. The theory behind it is based on least squares solution. The same fundamental theory behind lsqrcurvefit! I’ll be exploring least squares in upcoming posts.

Incorporate the update scheme with your system of equations

% Choose number of iterations
for i = 1:6
    for i=1:length(lat)
        diffx = (lon(i)-m(2))*lonkm;
        diffy = (lat(i)-m(3))*latkm;
        D0(i) = sqrt(diffx^2+diffy^2);
        d(i) = tim(i) - D0(i)/vp - m(1);
        G(i,3) = -(diffy/D0(i))/vp; G(i,2) = -(diffx/D0(i))/vp; G(i,1) = 1;
    end
% Least Squares Solution
dm = inv(G'*G)*G'*d'
% change dm to degrees as our coordinates are in degrees
dm(2) = dm(2)/lonkm;
dm(3) = dm(3)/latkm;
% Update 'm' Vector
m = m+dm'
% Plot Each Iterative Solution as a blue circle
plot(m(2),m(3),'o')
end
 
% Plot Final solution and notate with a black diamond
plot(m(2),m(3),'kd')

The key step here is to create another for loop to carry out the same calculations with each new and improved guess. I’ve set it to do it 6 times. You can experiment how many times you need to arrive at the zero error solution.

Now combine all the code together

% Finding Earthquake Location from Recorded Data
 
% Import Data Files
loc = importdata('loctim.txt');
bound = importdata('border.xy');
lat = loc(:,1);     % latitude
lon = loc(:,2);     % longitude
ele = loc(:,3);     % elevation
tim = loc(:,4);     % time of arrival 
 
% Plot Switzerland
plot(bound(:,1),bound(:,2))
hold on;
% Plot Stations
plot(lon,lat,'*')
% Initial Guess
lat0 = 46.9; lon0 = 9;
plot(lon0,lat0,'ro')
plot(lon(6),lat(6),'go')
 
% Parameters to change degrees to kilometers
latkm = 111.19; lonkm = 75.82; vp = 5.8;
 
 
% Initial Calculation
t0 = sqrt( ((lon(6)-lon0)*lonkm).^2 + ((lat(6)-lat0)*latkm).^2 )/vp
 
 
m(1) = t0; m(2)= lon0; m(3) = lat0;
 
% Choose number of iterations
for i = 1:6
    for i=1:length(lat)
        diffx = (lon(i)-m(2))*lonkm;
        diffy = (lat(i)-m(3))*latkm;
        D0(i) = sqrt(diffx^2+diffy^2);
        d(i) = tim(i) - D0(i)/vp - m(1);
        G(i,3) = -(diffy/D0(i))/vp; G(i,2) = -(diffx/D0(i))/vp; G(i,1) = 1;
    end
% Least Squares Solution
dm = inv(G'*G)*G'*d'
% change dm to degrees as our coordinates are in degrees
dm(2) = dm(2)/lonkm;
dm(3) = dm(3)/latkm;
% Update 'm' Vector
m = m+dm'
% Plot Each Iterative Solution as a blue circle
plot(m(2),m(3),'o')
end
 
% Plot Final solution and notate with a black diamond
plot(m(2),m(3),'kd')

You should get the following

dm =
1.0e-006 *

-0.0262
0.4140
0.9619

m =
40.2578 7.5580 47.2169

dm can be approximated to 0. So our algorithm tells us that the epicenter is at the black diamond

latitude 47.2169 and longtitude 7.5580 and occurred at time 16:30:40.3

Here are the links to download the source code, border.xy, and loctim.txt

Conclusion

I know very few people will be finding epicenters after reading this example, but I hope going through this exercise has reveal to you the basics of approaching a problem through matlab. I know I did not feature a matlab function in this post, but this is to encourage every user that matlab functions are tools to assist our analysis and not as a replacement. I hope this post has also got you interested in linear algebra and least squares solution. I’ll be writing about those topics in detail in future posts.

Was this example helpful? Are there problems in your field of interest that can benefit from this algorithm (i.e. locating whales and foreign submarines)? Please share.

Special thanks to Professor Malcolm Sambridge who taught me Inverse Problems this past quarter.