The Computer Program in Seismology location program elocate is a simple program that can easily be run in a batch mode. The program requires only a velocity model file, VEL.MOD and the data file elocate.dat. The program is described in Chapter 4 of the distributed PROGRAMS.330/DOC/GSAC.pdf/cps330g.pdf.
There are many references to the "leave one out" jackknife procedure for estimating means and variances. A recent reference is
G. A. Prieto, D. J. Thomas, F. L. Vernon, P. M. Shearer and R. L. Parker, (2007). Confidence intervals for earthquake source parameters, Geophys. J. Int 168, 1227-1234. doi:10.1111/j.1365-246X.2006.03257.x
The essence of the the procedure is to systematically leave one observation out of a data set, and then to run the inversion procedure. Then the inverted results are examined. Define K as the original number of observations. The result of each inversion will be a Vj for j=1,K. The mean value is
Vbar = (1/K) SUM Vjand the variance is
Var = (K-1)/K SUM [Vj - Vbar]2
To run this procedure, you must have installed the Computer Program in Seismology location program elocate. You must also have access to the programs awk/gawk and the C compiler gcc.
Download the following:
or t.tgz and then gunzip -c t.tgz | tar xvf - to get the DOIT and elocate.dat files.
Make the shell script executable:
chmod +x DOIT
DOIT
or
./DOIT
All processing results are placed in a sub-directory RESULT:
ls RESULT
elocate.sum.1 elocate.sum.13 elocate.sum.3 elocate.sum.7 elocate.txt.10 elocate.txt.14 elocate.txt.4 elocate.txt.8
elocate.sum.10 elocate.sum.14 elocate.sum.4 elocate.sum.8 elocate.txt.11 elocate.txt.15 elocate.txt.5 elocate.txt.9
elocate.sum.11 elocate.sum.15 elocate.sum.5 elocate.sum.9 elocate.txt.12 elocate.txt.2 elocate.txt.6
elocate.sum.12 elocate.sum.2 elocate.sum.6 elocate.txt.1 elocate.txt.13 elocate.txt.3 elocate.txt.7
cat RESULT/elocate.sum.10The entries are the RMS error, latitude, longitude, depth, origin time as YYYYMMDDHHmmSS.MSEC and the epoch time in seconds after January 1, 1970.
0.035 37.9515 -77.9603 5.69 20111026223658.022 1319668618.022
The listed output of the DOIT script is:
37.9501 0.0040 -77.9569 0.0059 5.6960 0.4665 1319668617.987 0.057 20111026223657.987where the columns are
Latitude Dlat Longitude Dlon Depth Ddep Epoch_time Dtime and origin time in human form.
The location estimate obtained using all data is in the elocate.txt output:
Error Ellipse X= 0.3721 km Y= 0.6225 km Theta = 130.7784 deg
RMS Error : 0.036 sec
Travel_Time_Table: CUS
Latitude : 37.9500 +- 0.0048 N 0.5303 km
Longitude : -77.9570 +- 0.0057 E 0.4946 km
Depth : 5.69 +- 0.88 km
Epoch Time : 1319668617.988 +- 0.11 sec
Event Time : 20111026223657.988 +- 0.11 sec
Event (OCAL) : 2011 10 26 22 36 57 988
HYPO71 Quality : BB
Gap : 102 deg
The jackknife estimated standard error are similar to those for latitude and longitude, and smaller for depth and origin time for this data set.
#!/bin/sh
######
# preserve the original elocate.dat
######
cp elocate.dat elocate.dat.save
#####
# create a work subdirectory and populate it
# if it already exists start anew!
#####
rm -fr TEMP
mkdir TEMP
cp VEL.MOD TEMP
#####
# create a result directory to preserve the
# different run output
#####
rm -fr RESULT
mkdir RESULT
#####
# initialize to get the VEL.MOD file
#####
elocate -VELMOD
#####
# define the velocity model and starting depth
#####
MODEL=4
DEPTH=10
#####
# determine the number of lines in elocate.dat
#####
NLINE=`wc elocate.dat.save | awk '{print $1}' `
#####
# now loop over the lines
#####
COUNT=0
while [ ${COUNT} -lt ${NLINE} ]
do
head -${COUNT} elocate.dat.save > TEMP/elocate.dat
COUNTM1=`echo ${COUNT} ${NLINE} | awk '{print $2 -$1 -1}' `
tail -${COUNTM1} elocate.dat.save >> TEMP/elocate.dat
COUNT=` echo ${COUNT} | awk '{print $1 + 1}' `
#####
# use subshell to do the work
#####
(cp VEL.MOD TEMP; cd TEMP; elocate -M ${MODEL} -D ${DEPTH} -BATCH > elocate.txt)
#####
# preserve the results with a unique name
#####
mv TEMP/elocate.txt RESULT/elocate.txt.${COUNT}
mv TEMP/elocate.sum RESULT/elocate.sum.${COUNT}
done
#####
# Use Jackknife procedure to get location and one sigma
#####
#####
# Reference: Equations (1) - (4)
# G. A. Prieto, D. J. Thomas, F. L. Vernon, P. M. Shearer
# and R. L. Parker, (2007). Confidence intervals for earthquake
# source parameters, Geophys. J. Int 168, 1227-1234.
# doi:10.1111/j.1365-246X.2006.03257.x
#####
# scan the elocate.sum.* files
# to perform the leave one out
#
# the entries in elocate.sum are
# RMS LAT LON DEP OT OT_CSS
# 0.007 37.9478 -77.9593 5.64 20111026223658.017 1319668618.017
#####
# make lists of the variables
#####
rm -f list.lat list.lon list.dep list.otcss
cat RESULT/elocate.sum.* | awk '{print $2}' > list.lat
cat RESULT/elocate.sum.* | awk '{print $3}' > list.lon
cat RESULT/elocate.sum.* | awk '{print $4}' > list.dep
cat RESULT/elocate.sum.* | awk '{print $6}' > list.otcss
#####
# now get the jackknife estimates for each parameter.
# this ignores any possible covariance
#####
for P in LAT LON DEP OT
do
case ${P} in
LAT) IND=1 ;FILE=list.lat ;FMTSD="%9.4f" ;;
LON) IND=2 ;FILE=list.lon ;FMTSD="%9.4f" ;;
DEP) IND=3 ;FILE=list.dep ;FMTSD="%9.4f" ;;
OT) IND=4 ;FILE=list.otcss ;FMTSD="%20.3f" ;;
esac
#####
# M is the number of entries in the file
#####
M=`wc ${FILE} | awk '{print $1}' `
#####
# create and awk scrip to compute the mean
#####
cat > awkscmn << EOF
BEGIN {SUM=0.0}
{SUM+=\$1}
END { printf "${FMTSD}\n", SUM/NR}
EOF
MEAN=`cat ${FILE} | awk -f awkscmn `
#####
# save the mean in an array
#####
MVAL[${IND}]="${MEAN}"
#####
# now computed the standard error from the square root of the variance
#####
cat > awkscsd << EOF
BEGIN {SUM=0.0}
{
SUM+= (\$1 -(${MEAN}) )*(\$1 -(${MEAN}) )
}
END { printf "${FMTSD}\n", sqrt( ( NR - 1 ) * SUM / NR ) }
EOF
SD=`cat ${FILE} | awk -f awkscsd`
SDVAL[${IND}]="${SD}"
done
rm -f awkscmn awkscsd list.dep list.lat list.lon list.otcss
#####
# at this point compile a C program uner LINUX to convert
# epoch time to human time
#####
#####
# now compile a C program to convert the time in seconds from
# January 1, 1970 00:00:00.000 to UTC
#####
cat > c.c << EOF
#include<stdio.h>
#include <time.h>double TIME=${MVAL[4]} ; time_t TIMES; struct tm tm; char ostr[30]; main() { TIMES=(time_t)TIME; gmtime_r( &TIMES, &tm); strftime(ostr,sizeof(ostr),"%Y%m%d%H%M%S", &tm); printf("%s.%3.3d\n",ostr,(int)( 1000*(TIME - (double)TIMES )+0.49)); } EOF gcc c.c TIMEHUMAN=`a.out` ###### # now format the nice summary ##### echo ${MVAL[1]} ${SDVAL[1]} ${MVAL[2]} ${SDVAL[2]} ${MVAL[3]} ${SDVAL[3]} ${MVAL[4]} ${SDVAL[4]} ${TIMEHUMAN} rm -f a.out c.c