#!/bin/sh
set -x
#####
# define the location of the deconvolved, QCd waveforms
# of both the large event, DIRBIG and the smaller event used
# for the empirical Green function, e.g.,
#DIRBIG=/home/rbh/PROGRAMS.310t/MOMENT_TENSOR/MECH.TEL/20070815234057/DAT.TEL
#DIRSMALL=./SYN
# DO NOT CHANGE THE LOCATION OF THE SYNTHETICS
#####
DIRBIG=path_to_big_event_waveform_directory
#####
# define the mechanism and depth for the synthetics
# These values can be obtained from the regional
# averages of mechanisms using the new code
# STK = strike
# DIP = dip
# RAKE = rake
# HS = sourrce depth in km
# MW = this is the reference MW, in theory the zero frequency
# level of the derived source time function is the moment ratio
# of the two events, e.g., a ratio of 1000 corresponds to a
# delta Mw of 2.0, so if thw reference Mw=5 then the big event has Mw=7!
#####
STK=171
DIP=55
RAKE=112
HS=30
MW=5.0
#####
# ONLY CHANGE BELOW HERE TO FINE TUNE
#####
DIRSMALL=./SYN
#####
# make the synthetics for all vertical components and
# create synthetic ground velocities at each station in units of m/s
# Place the synthetics in a subdirectory called SYN
#####
DOMKSYN ${DIRBIG} ${STK} ${DIP} ${RAKE} ${HS} ${MW}
#####
# define the CUT window with respect to P
#####
CUTL=-40
CUTH=200
#####
# define the frequency band
#####
FL=0.01
FH=0.20
#####
# define the DELTA for the traces used in the decon
#####
INTDELTA=0.25
#####
# modified 10 AUG 2009 from the script used for Wells and L'Aquila
# Previously:
# ls $DIRBIG $DIRSMALL | sort | uniq > uniq.list
# this will not work since early IRIS data did not use location codes
# So
# 1. Assume thaat only one location code is used
#` 2. build up a list of unique STATION COMPONENT codes
#####
#####
# initialize by cleaning up
#####
rm -f stacmp
rm -f *.decon
#####
# create a synthetic for each of the observed waveforms
#####
if [ -d DECONDIR ]
then
rm -fr DECONDIR
fi
mkdir DECONDIR
for i in ${DIRBIG}/*[ZRT] ${DIRSMALL}/*[ZRT]
do
saclhdr -KSTNM -KCMPNM -NL $i >> stacmp
done
cat stacmp | sort | uniq > uniqstacmp
#####
# do the deconvolution
#####
#####
# select the ALPHA parameter
# The approximate relation between the ALPHA and Gaussian pulse duration is
# ALP Duration Corner Frequency
# 0.5 8 s ~ 0.2 Hz
# 1.0 4 s ~ 0.4 Hz
# 2.5 1.6s ~ 1 Hz
#
# The ALP parameter will affect the resolution of the final fit
# A value of 2.5 may be too extreme for real data
#
# The choice of the ALPHA parameter affects the resolution of the
# resulting empirical source time function
#
# The interative decon program computes a measure of fit, we use the
# LEVEL parameters to guide the acceptance of the results
#
# NOTE: for synthetics the fits will generally be good and the LEVEL
# does little in selecting good traces
#####
for ALP in 0.5 1.0
do
case ${ALP} in
0.5) LEVEL=70. ;;
1.0) LEVEL=60. ;;
2.5) LEVEL=50. ;;
esac
#####
# select the component
# later use a case statement to set the cut, e.g., a vs t1
#####
for COMP in HZ
do
grep ${COMP} uniqstacmp > unuq.list
rm -fr BIG.DIR SML.DIR
mkdir BIG.DIR
mkdir SML.DIR
while read KSTNM KCMPNM
do
if [ -f ${DIRBIG}/${KSTNM}*${KCMPNM} ]
then
if [ -f ${DIRSMALL}/${KSTNM}*${KCMPNM} ]
then
cp ${DIRBIG}/${KSTNM}*${KCMPNM} BIG.DIR/${KSTNM}${KCMPNM}
cp ${DIRSMALL}/${KSTNM}*${KCMPNM} SML.DIR/${KSTNM}${KCMPNM}
gsac << EOF
cuterr fillz
cut a ${CUTL} a ${CUTH}
r BIG.DIR/${KSTNM}${KCMPNM} SML.DIR/${KSTNM}${KCMPNM}
hp c ${FL} n 3
lp c ${FH} n 3
interpolate delta ${INTDELTA}
w N D
q
EOF
saciterd -FN N -FD D -2 -ALP ${ALP} -POS -D 20 -N 100
USER5=`saclhdr -USER5 decon.out`
#####
# only use decon that filt 80%
#####
cat > awkprog << EOF
{
ANS="YES"; if(\$1 < ${LEVEL}) ANS="NO" ; print ANS
}
EOF
ANS=`echo $USER5 | awk -f awkprog`
if [ ${ANS} = "YES" ]
then
mv decon.out ${KSTNM}${KCMPNM}.${ALP}.decon
fi
fi
fi
done < unuq.list
#####
# make the map
#####
gsac << EOF
rh BIG.DIR/* SML.DIR/*
map kstnm on r on
quit
EOF
sh map.sh
gm convert -trim map.eps ${COMP}.${ALP}.map.png
done
for i in Z
do
case $i in
Z) KOLOR=2 ;;
T) KOLOR=4 ;;
esac
gsac << EOF
bg plt
r *${i}.${ALP}.decon
prs shd pos color ${KOLOR} tl -20 200 amp 0.3 vl 0 360 az sa 0 ann sta
q
EOF
if [ -f PRS001.PLT ]
then
plotgif -C16 -F7 -W10 -K < PRS001.PLT > ${i}decon.gif
fi
done
mv Zdecon.gif Zdecon.${ALP}.gif
#mv Tdecon.gif Tdecon.${ALP}.gif
mv *.decon DECONDIR
done
#####
# cleanup
#####
rm -fr BIG.DIR
rm -fr SML.DIR
rm -fr numerator denominator N D observed predicted
rm -fr stacmp PRS001.??? awkprog
rm -fr uniqstacmp unuq.list
rm -fr decon.out
|
This is the most important part since this this tells where
to find the waveforms assiciated with the large earthquake. this
directory has waveform files that end in 'Z'
'R' or 'T' Note we do not need a second directory since we will
create synthetics for the small earthquake.
We must also tell the program the focal mechanism, depth and target Mw
for the small earthquake. Note that the Mw does not affect the duration
of the source time function which has a duration of 4 * 4 *0.20 sec =
3.2 sec Perhaps this is too long
We will place our created small earthquake waveforms in a local
directory called SYN. The script DOMKSYN will create this
directory
Make the synthetics
This says that we will consider the time window 40 seconds before and
200 seconds after the predicted P-wave first arrival
If we filter both sets of waveforms identically, the deconvolution
should yield the same result. However, the high-pass removes some low
frequency noise.
To run the deconvolution, both traces must have the same sample rate.
In addition resampling make the process faster since we do not wish to
fit high frequencies. This value may be adequate when the deconvolved
function is much longer than a few seconds. When the function is
shorter, this must be smaller
Do the deconvolution for these two values of the filter parameter ALP
This is set up for the P-wave. First because we use the vertical
component and because we cut with respect to the P-wave first arrival
time. If you wish to use SH, then the loop should start as
for COMP in HT
and the gsac command cut should be with 't0' instead of 'a'
These two 'if' statements ensure that the station-component is in both
directories, then
copy them to a local work directory
cut both waveforms with respect to P
prefilter both waveforms
resample to a smaller dt
save as N for numerator (big event) and D for denominator (small event)
Do the interative deconvolution and examine the goodness of fit
parameter. To see how 'saciterd' works
saciterd -h
If the fit matches the desired level, save the result with a name that
reflects station, component, ALP. Since we work with synthetics the
fits will be good.
All processing is complete for this value of ALP
Use the 'gsac' command to make a GMT script to plot the stations and
ray paths used for the deconvolution
Use the GraphicsMagick package to convert the map.eps to a PNG
graphics file. If you have ImageMagick installed, the command is
convert -trim map.eps ${COMP}.${ALP}.map.png
All processing is complete
Using 'gsac' plot some aximuthal record sections of the deconvolved
traces.
Convert these directly to a 'gif' image instead of an eps, and then a
PNG
Clean up
|