Earthquake Location

MATLAB/Octave Code

A simple earthquake location program was written and can be run under MATLAB/Octave. This discussion will focus on Octave runnding on LINUX/OSX.

1. Download the following files:

2. Create a directory in your login directory called OctaveCodes
            mkdir ~/OctaveCodes
3. Create a sub-directory Location
            mkdir ~/OctaveCodes/Location
4. Unpack the code base into the Location sub-directory
            unzip location.zip
5. Start octave. Tell octave where to find the source code
            octave
            addpath('~/OctaveCodes/Location')
            savepath
            quit
    The result of these two commands is that when you start octave again, the codes will be found.

Data Set

Now create a directory for this test. For example,
            mkdir Location
            cd Location
            unzip locationdata.zip

Test Runs

cd Location
cd ELOCATE01

You will see an elocate.txt which was created by the FORTRAN program elocate. Look at this file, and at the bottom you will see the velocity model used, e.g., WUS or CUS. In this case the WUS velocity model was used. If you look at what is in the directory, you see some other files:
WUS.txt CUS.txt HALF.txt SCM.txt
These are the velocity model files which have a very simple format. For example, CUS.txt is

     H       P      S     Pg        Lg
  1.0000  5.0000  2.8900  6.1      3.55
  9.0000  6.1000  3.5200  6.1      3.55
 10.0000  6.4000  3.7000  6.1      3.55
 20.0000  6.7000  3.8700  6.1      3.55
  0.0000  8.1500  4.7000  6.1      3.55

The first line is an important identifier, and the following lines give layer thickness and velocity of that phase in the layer. You will see that this model supports P, S, Pg and Sg phases. The travel time routine computes the first arrival time of each phase. If you have a picked phase in your data set, e.g., PKP, and if that phase is not in the header line, then that phase is not used for the location.

The routine getinput.m knows the names of these files and provide a menu to select a velocity model for use. I mention this in case you wish to use another velocity model.

You will also see the file elocate.dat. This file is created by the program
            sac2eloc sacfiles
            sac2eloc *SAC
                   ...
sac2eloc reads the A, T0 and T1 times from the headers of the sacfiles. I usually set the times using the ppk regional command of gsac.

The elocate.dat has entries like the following:

BRLDA  Z  2016 01 12 18 27 35.279        1 i S         0 X         1   54.0920 -117.4038     1224.        0  9999999
BRLDA  Z  2016 01 12 18 27 30.026        2 e P         2 -         2   54.0920 -117.4038     1224.        0  9999999
NBC5   Z  2016 01 12 18 28 29.705        3 e P         2 +         3   57.5230 -122.6700     1154.        0  9999999
STPRA  Z  2016 01 12 18 28 10.986        4 e Lg        2 X         4   55.6606 -115.8323      761.        0  9999999
STPRA  Z  2016 01 12 18 27 51.063        5 i P         0 D         5   55.6606 -115.8323      761.        0  9999999
SWHSA  Z  2016 01 12 18 27 43.119        6 e S         2 X         6   54.8994 -116.7518      914.        0  9999999

The entries are station name, component, year, month, day, hour, minute, second of the phase arrival time, and ID (not used), the quakity of the phase (i = sharp, e = emergent),, the phase name, the phase weight (0=weight of 1.0, 2 is a weight of 0.5, and 4 is a weight of 0.0; these are HYPO71 weights), the P-wave first motion polarity, another ID (not used), station latitude, station longitude, station elevation (meters), ondate (not used), offdate (not used). These entries provide everything required to locate the event.

Sample Run

GNU Octave, version 3.4.0
Copyright (C) 2011 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. For details, type `warranty'.

Octave was configured for "x86_64-apple-darwin10.7.3".

Additional information about Octave is available at http://www.octave.org.

Please contribute if you find this software useful.
For more information, visit http://www.octave.org/help-wanted.html

Read http://www.octave.org/bugs.html to learn how to submit bug reports.

For information about changes from previous versions, type `news'.

octave-3.4.0:1> elocate

Select Velocity model
     CUS [0]
     WUS [1] 
     SCM [2] 
    HALF [3] 
Respond 0, 1, 2 or 3: 1
Enter depth: 10
Fix depth? y/n: n
NID   STA   PHASE    DIST(km)   AZIMUTH        BAZ          T        DTDR       DTDH        WT        RES  IO(P) (1)
  1  BRLDA      S      12.89     210.36      30.36       5.29        0.24       0.12       0.71      -0.04
  2  BRLDA      P      12.89     210.36      30.36       3.14        0.14       0.07       0.71      -3.14   117
  6  SWHSA      S      86.34      24.35     204.35      24.77        0.27       0.01       0.50     -11.68
  7  SWHSA      P      86.34      24.35     204.35      14.76        0.16       0.00       1.00     -10.08    92
 19  TD09A      P     153.28     156.80     336.80      25.43        0.16       0.00       1.00      -2.80    91
 14  TD010      P     184.32     159.57     339.57      30.38        0.16       0.00       1.00      -3.32    90
  5  STPRA      P     188.43      29.93     209.93      30.94        0.13      -0.10       1.00      -9.90    53
 11  TD008      S     197.86     141.26     321.26      54.58        0.22      -0.16       0.50       0.11
 12  TD008      P     197.86     141.26     321.26      32.13        0.13      -0.10       1.00      -3.23    53
 22  WTMTA      P     207.82     323.49     143.49      33.39        0.13      -0.10       1.00      -9.38    53
  9  TD002      P     208.93     113.61     293.61      33.53        0.13      -0.10       1.00      -5.75    53
 15  TD011      P     217.79     147.00     327.00      34.65        0.13      -0.10       1.00      -3.60    53
 13  TD009      P     218.07     162.60     342.60      34.69        0.13      -0.10       0.71      -3.40    53
 17  TD013      P     239.91     140.89     320.89      37.45        0.13      -0.10       1.00      -3.90    53
 18  TD029      P     260.49     147.46     327.46      40.06        0.13      -0.10       1.00      -3.59    53
 16  TD012      P     262.03     150.98     330.98      40.25        0.13      -0.10       1.00      -3.59    53
 20  TD13A      P     296.06     145.12     325.12      44.56        0.13      -0.10       0.71      -4.09    53
  3   NBC5      P     499.34     317.88     137.88      70.29        0.13      -0.10       0.71     -10.61    53
 21   WALA      P     616.97     157.69     337.69      85.18        0.13      -0.10       1.00      -5.18    53
------------------------------------------------------------------------------------------
     Lat        Lon          H        Epoch       Year Mo Dy  Doy  Hr Mm Sc Msc       SD     (2)
------------------------------------------------------------------------------------------
   54.4319  -117.3221      10.15   1452623243.301 2016 01 12 (012) 18 27 23 301      5.639
   54.4277  -117.3740      10.17   1452623243.156 2016 01 12 (012) 18 27 23 156      0.551
   54.4291  -117.3771      11.66   1452623243.222 2016 01 12 (012) 18 27 23 222      0.548
   54.4282  -117.3663      10.58   1452623243.249 2016 01 12 (012) 18 27 23 249      0.360
   54.4272  -117.3609      10.20   1452623243.265 2016 01 12 (012) 18 27 23 265      0.326
   54.4283  -117.3584      11.06   1452623243.322 2016 01 12 (012) 18 27 23 322      0.315
   54.4271  -117.3576      10.27   1452623243.283 2016 01 12 (012) 18 27 23 283      0.323
   54.4268  -117.3572      10.11   1452623243.279 2016 01 12 (012) 18 27 23 279      0.312
   54.4282  -117.3572      11.05   1452623243.328 2016 01 12 (012) 18 27 23 328      0.310
   54.4270  -117.3572      10.26   1452623243.285 2016 01 12 (012) 18 27 23 285      0.322
   54.4268  -117.3570      10.11   1452623243.280 2016 01 12 (012) 18 27 23 280      0.312
   54.4282  -117.3572      11.04   1452623243.328 2016 01 12 (012) 18 27 23 328      0.310
   54.4270  -117.3571      10.26   1452623243.285 2016 01 12 (012) 18 27 23 285      0.322
   54.4268  -117.3570      10.11   1452623243.280 2016 01 12 (012) 18 27 23 280      0.312
   54.4282  -117.3572      11.04   1452623243.328 2016 01 12 (012) 18 27 23 328      0.310
   54.4270  -117.3571      10.26   1452623243.285 2016 01 12 (012) 18 27 23 285      0.322
   54.4268  -117.3570      10.11   1452623243.280 2016 01 12 (012) 18 27 23 280      0.312
   54.4282  -117.3572      11.04   1452623243.328 2016 01 12 (012) 18 27 23 328      0.310
NID   STA   PHASE    DIST(km)   AZIMUTH        BAZ          T        DTDR       DTDH        WT        RES  IO(P) (3)
  1  BRLDA      S      37.50     184.63       4.63      11.75        0.27       0.03       0.69       0.20
  2  BRLDA      P      37.50     184.63       4.63       6.99        0.16       0.02       0.68      -0.29    97
  6  SWHSA      S      65.28      36.61     216.61      19.15        0.27       0.02       0.49       0.64
  7  SWHSA      P      65.28      36.61     216.61      11.41        0.16       0.01       1.15      -0.04    93
  5  STPRA      P     167.98      35.33     215.33      27.77        0.16       0.00       1.15      -0.04    91
 19  TD09A      P     178.88     159.13     339.13      29.51        0.16       0.00       1.14      -0.19    91
 22  WTMTA      P     184.91     319.58     139.58      30.39        0.13      -0.10       1.10       0.32    53
 14  TD010      P     210.19     161.21     341.21      33.59        0.13      -0.10       1.14       0.17    53
 11  TD008      S     220.77     144.88     324.88      59.38        0.22      -0.16       0.26       2.02
 12  TD008      P     220.77     144.88     324.88      34.93        0.13      -0.10       0.96       0.67    53
  9  TD002      P     223.33     119.49     299.49      35.25        0.13      -0.10       0.91      -0.78    53
 15  TD011      P     241.84     149.75     329.75      37.59        0.13      -0.10       1.14       0.15    53
 13  TD009      P     244.18     163.69     343.69      37.89        0.13      -0.10       0.81       0.09    53
 17  TD013      P     262.63     143.98     323.98      40.23        0.13      -0.10       1.15       0.02    53
 18  TD029      P     284.56     149.77     329.77      43.00        0.13      -0.10       1.14       0.17    53
 16  TD012      P     286.72     152.96     332.96      43.28        0.13      -0.10       1.15       0.09    53
 20  TD13A      P     319.61     147.36     327.36      47.44        0.13      -0.10       0.79      -0.27    53
  3   NBC5      P     477.17     316.15     136.15      67.38        0.13      -0.10       0.58      -1.01    53
 21   WALA      P     642.47     158.33     338.33      88.31        0.13      -0.10       0.61      -1.61    53
------------------------------------------------------------------------------------------
     Lat        Lon          H        Epoch       Year Mo Dy  Doy  Hr Mm Sc Msc       SD     (4)
------------------------------------------------------------------------------------------
   54.4270  -117.3571      10.26   1452623243.285 2016 01 12 (012) 18 27 23 285      0.322
    0.0066     0.0237       1.69            0.161 
------------------------------------------------------------------------------------------
Confidence ellipse:                (5)
     Major Axis Length 1.535329
     Minor Axis Length 0.737951
     major Axis Azimuth   269

Comments:
(1) Given an initial location near the station with the earliest arrival time, these are the distances and residuals.
(2) These show the beginning of each iteration and the rms error of the residuals. As iterations continue, the source parameters change and the rms (S) becomes smaller. You will no a bit of instability in the solution
(3) These are all parameters after the last iteration.
(4) These are the event parameters with estimated errors. The confidence can be obtained by multiplying by the value of the Student t-distribution, which requires the number of observations and number of parameters determined.
(5) This defines the error ellipse. To obtain the confidence values, you must multiply these lengths by the F-distribution value for the number of observations and the number of parameters determined. See

Flinn, E. A. (1965). Confidence regions and error determinations for seismic event location, Reviews of Geophysics 3, 157-185

for details.

Your earthquake

In this section you will pick arrival times and then locate the earthquake.

Get the data set DAT.zip

Create the processing directory:

mkdir ELOCATE10
cd ELOCATE10
unzip DAT.zip
cp ../ELOCATE05/CUS.txt ../ELOCATE05/WUS.txt ../ELOCATE05/SCM.txt ../ELOCATE05/HALF.txt .
cd DAT

Now in the DAT directory start gsac to interactively pick arrival times. Since this is an exercise, the Sac files have the distance already assigned:

gsac
GSAC> r *
GSAC> sort up dist
GSAC> hp c 1 n 2
GSAC> ppk regional perplot 3
GSAC> wh
GSAC> q

cd ..
sac2eloc DAT/*

When using the ppk regional you can use the XX, +, -, mouse wheel to adjust the trace. Then select the phase, P or S, and then read the instruction bar at the top.

To go to the next set of traces, hit the 'n' for next page, or 'b' for back to previous page.

The wh is used to finish just to record the arrival time picks. you can see what you did in gsac by using the 'lh' list header command.

The sac2eloc reads the header of the sac files, and sames the arrival time and station information to create the elocate.dat file.

Locate the earthquake

Use the CPS program elocate
```
elocate -VELMOD
elocate 
```
Use the MATLAB/octave program elocate to locate the earthquake

More development

The routine getdistance.m justs uses a local cartesian coordinate system to compute distances and azimuths. This is OK as long as the distances are short and if the location is not near the poles. The next step would be to use spherical coordinates to compute the distances:
http://www.eas.slu.edu/eqc/eqc_course/IntroSeis/Ass07/Ass07.pdf

or ellipsoidal

http://www.eas.slu.edu/eqc/eqc_course/IntroSeis/Ass08/Ass08.pdf