2011-10-30 20:13
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 Vj
and 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
Execute the shell script:
DOIT or ./DOIT
All processing results are placed in a sub-directory RESULT:
ls RESULT elocate.sum.1 elocate.sum.14 elocate.sum.5 elocate.txt.1 elocate.txt.14 elocate.txt.5 elocate.sum.10 elocate.sum.15 elocate.sum.6 elocate.txt.10 elocate.txt.15 elocate.txt.6 elocate.sum.11 elocate.sum.2 elocate.sum.7 elocate.txt.11 elocate.txt.2 elocate.txt.7 elocate.sum.12 elocate.sum.3 elocate.sum.8 elocate.txt.12 elocate.txt.3 elocate.txt.8 elocate.sum.13 elocate.sum.4 elocate.sum.9 elocate.txt.13 elocate.txt.4 elocate.txt.9
The files starting with elocate.txt are the complete output of
the elocate program run in batch mode.
The files starting with elocate.sum are the simple summary of
the location:
cat RESULT/elocate.sum.10
0.035 37.9515 -77.9603 5.69 20111026223658.022 1319668618.022
The entries are the RMS error, latitude, longitude, depth, origin time as YYYYMMDDHHmmSS.MSEC and the epoch time in seconds after January 1, 1970.
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.987
where 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