This exercise shows how to estimate phase velocities beneath a
network from a teleseismic surface wave signal observed by all
stations. The concept is to assume that the Earth structure
beneath the seismic network is uniform and that the incident
wavefront from the distance source follows a great circle path so
that the signal arrivals at two stations at the same distance at the
same time. The analysis essentially projects the observed waveforms
onto a linear array in the center of the network and the p-omega
stacking is then performed.
In order to do this, data
processing must account for the requirements of the stacking program
sacpom96 (which is called by the GUI do_pom). First the
traces must be of the same length, second, the traces must be
resampled to a lower sample rate for efficiency in the stacking (to
make the fast Fourier transforms faster because the number of points
is smaller), the fundamental mode surface wave must be isolated for a
good stack (using the GUI do_mft which calls sacmft96
to determine the group velocity dispersion and then sacmat96
to isolate the fundamental mode. Finally do_pom cat be run.
The data set used for this tutorial is Sac.tgz
Download this file and unpack is using the command
gunzip -c Sac.tgz | tar xvf -
This will create a subdirectory Sac and place the waveform forms in
that directory.
Since these SAC files are coming from an Intel
machine, you must make sure that the files are in the correct byte
order by performing the commands
for i in *Sac
do
saccvt -I < $i > tmp
mv tmp $i
done
You must put the event information into the trace headers. you can
do this using gsac as follows (I assume that you are in the
Sac directory)
gsac
GSAC - Computer Programs in Seismology [V1.1.21 13 SEP 2007]
Copyright 2004, 2005, 2006, 2007 R. B. Herrmann
GSAC> r *Z.Sac
CHUBHZ.Sac HDBBHZ.Sac KANBHZ.Sac KWJBHZ.Sac PUSBHZ.Sac
SEOBHZ.Sac SOGBHZ.Sac TAGBHZ.Sac TEJBHZ.Sac TJNBHZ.Sac ULLBHZ.Sac
GSAC> ch evla 6.898 evlo 126.579 evdp 33
GSAC> ch ocal 2001 01 01 06 57 01 172
GSAC> wh
GSAC> quit
Since I assume that the trace headers already had the station latitude and longitude fields set (STLA and STLO), the trace headers will now have the DIST, GCARC, AZ and BAZ fields set.
If I
have GMT installed, I can do the following from within gsac
r *Z.Sac
CHUBHZ.Sac HDBBHZ.Sac KANBHZ.Sac KWJBHZ.Sac PUSBHZ.Sac
SEOBHZ.Sac SOGBHZ.Sac TAGBHZ.Sac TEJBHZ.Sac TJNBHZ.Sac ULLBHZ.Sac
GSAC> map5 r on
r on
Execute using the command: sh map5.sh
GSAC> sh map.sh
sh map5.sh
pscoast: Working on block # 301
pscoast: Adding Borders...
GSAC>
The file map5.eps is created and you can view it using gs, display or other PostScript viewers. You will see the following: map5.png (Note the image was created using GMT4 and the mao command)
The next step is to view the traces and to examining the header
for the sample interval since the analysis requires that all traces
have the same sample interval, DELTA, and the same number of points.
Note that these traces are the original digital data. The
instrument response has not been removed. Normally deconvolution is
necessary, but if the instruments have identical response, then
deconvolution is not necessary. Such is the case for these data.
SAC> r *Z.Sac
CHUBHZ.Sac HDBBHZ.Sac KANBHZ.Sac KWJBHZ.Sac PUSBHZ.Sac SEOBHZ.Sac
SOGBHZ.Sac TAGBHZ.Sac TEJBHZ.Sac TJNBHZ.Sac ULLBHZ.Sac
GSAC> sort up dist
Sorting on DIST in ascending order
GSAC> lh delta npts
SOGBHZ.Sac (6):
NPTS 6601 DELTA 0.2
KWJBHZ.Sac (3):
NPTS 6601 DELTA 0.2
PUSBHZ.Sac (4):
NPTS 6601 DELTA 0.2
HDBBHZ.Sac (1):
NPTS 80001 DELTA 0.05
TAGBHZ.Sac (7):
NPTS 6601 DELTA 0.2
TEJBHZ.Sac (8):
NPTS 6601 DELTA 0.2
TJNBHZ.Sac (9):
NPTS 80001 DELTA 0.05
SEOBHZ.Sac (5):
NPTS 6601 DELTA 0.2
ULLBHZ.Sac (10):
NPTS 6601 DELTA 0.2
KANBHZ.Sac (2):
NPTS 6601 DELTA 0.2
CHUBHZ.Sac (0):
NPTS 6601 DELTA 0.2
GSAC> p
You will not the different sampling intervals and the different
number of points, The command "p" plots the
traces:
We use this figure to design the cut to isolate the signal and
then we resample to a DELTA=0.5 sec with the following commands
r *Z.Sac
CHUBHZ.Sac HDBBHZ.Sac KANBHZ.Sac KWJBHZ.Sac PUSBHZ.Sac SEOBHZ.Sac
SOGBHZ.Sac TAGBHZ.Sac TEJBHZ.Sac TJNBHZ.Sac ULLBHZ.Sac
GSAC> rtr
GSAC> w
CHUBHZ.Sac HDBBHZ.Sac KANBHZ.Sac KWJBHZ.Sac PUSBHZ.Sac SEOBHZ.Sac
SOGBHZ.Sac TAGBHZ.Sac TEJBHZ.Sac TJNBHZ.Sac ULLBHZ.Sac
GSAC> synchronize o
GSAC> wh
GSAC> cuterr fillz
GSAC> cut o 300 o 1800
cut o 300 o 1800
O 300.000000 O 1800.000000
GSAC> r *Z.Sac
CHUBHZ.Sac HDBBHZ.Sac KANBHZ.Sac KWJBHZ.Sac PUSBHZ.Sac SEOBHZ.Sac
SOGBHZ.Sac TAGBHZ.Sac TEJBHZ.Sac TJNBHZ.Sac ULLBHZ.Sac
GSAC> lp c 1 n 2 p 2
LP: corner fc 1.000000 npoles 2 pass 2 Butterworth
GSAC> interpolate delta 0.5
delta 0.5
GSAC> cd ..
Current directory is
/home/rbh/PROGRAMS.310t/PROGRAMS.330.fixups/TUTORIAL/POMEGA
GSAC> mkdir GOOD
mkdir GOOD
GSAC> cd GOOD
Current directory is
/home/rbh/PROGRAMS.310t/PROGRAMS.330.fixups/TUTORIAL/POMEGA/GOOD
GSAC> w
CHUBHZ.Sac HDBBHZ.Sac KANBHZ.Sac KWJBHZ.Sac PUSBHZ.Sac SEOBHZ.Sac
SOGBHZ.Sac TAGBHZ.Sac TEJBHZ.Sac TJNBHZ.Sac ULLBHZ.Sac
GSAC> pwd
pwd
/home/rbh/PROGRAMS.310t/PROGRAMS.330.fixups/TUTORIAL/POMEGA/GOOD
GSAC>
This is what was done. Read the traces and remove and linear
trend, and then write the traces out. Reset the time reference as the
origin time using synchronize. Then plot the trace and see
that the surface wave is in the window 300 to 1800 seconds after the
origin time. Reread and add zeros to the beginning and end of traces
if necessary. The low pass filter with a zero phase filter at 1.0 Hz,
and finally interpolate to a sample of 0.5 seconds.
Since
I do not want to destroy the original data, I create a new directory
parallel to the Sac directory with the SHELL command cd ..
followed by mkdir GOOD. We then move from the current
directory to that new directory and write the traces there. If
I read these new traces in I see the following with a lh npts
delta
GSAC> lh npts delta
SOGBHZ.Sac (6):
NPTS 3001 DELTA 0.5
KWJBHZ.Sac (3):
NPTS 3001 DELTA 0.5
PUSBHZ.Sac (4):
NPTS 3001 DELTA 0.5
HDBBHZ.Sac (1):
NPTS 3001 DELTA 0.5
and the plot
In the GOOD directory, issue the command
do_mft *Z.Sac
Then select one trace, TAG, for example, and run do_mft and select the dispersion curve. Since this is a large earthquake, I instruct do_mft to look at periods between 4 and 300 seconds. I select dispersion between 20 and 200 seconds:
Now hit the Match button. on the next page, select Match and then Fund for the fundamental mode. Then select No and then Quit
If you look at the
console, you will see the following command when do_mft started
the phase match filter program sacmat96
/home/rbh/PROGRAMS.310t/PROGRAMS.330/bin/sacmat96 -F TAGBHZ.Sac -D disp.d -AUTO
We will use this as a prototype to phase match all traces through a
SHELL script. The reason for doing it this way is that there is no
certainty that all traces will have the
same range of dispersion
periods for and individual phase match filter. So create the
script DOIT. make it executable and run in the GOOD directory
which also has the disp.d dispersion file. Here is DOIT
#!/bin/sh for i in *Z.Sac do sacmat96 -F $i -D disp.d -AUTO done
Now after you execute it you see the following files in the GOOD directory
CHUBHZ.Sac disp.out KANBHZ.Sacr mft96.ctl PUSBHZ.Sac SOGBHZ.Sac TEJBHZ.Sac
CHUBHZ.Sacr DOIT* KANBHZ.Sacs mft96.disp PUSBHZ.Sacr SOGBHZ.Sacr TEJBHZ.Sacr
CHUBHZ.Sacs HDBBHZ.Sac KWJBHZ.Sac mft96.ods PUSBHZ.Sacs SOGBHZ.Sacs TEJBHZ.Sacs
disp.d HDBBHZ.Sacr KWJBHZ.Sacr MFT96.PLT SEOBHZ.Sac TAGBHZ.Sac TJNBHZ.Sac
disp.dp HDBBHZ.Sacs KWJBHZ.Sacs P001.PLT SEOBHZ.Sacr TAGBHZ.Sacr TJNBHZ.Sacr
disp.dv KANBHZ.Sac MFT96CMP P002.eps SEOBHZ.Sacs TAGBHZ.Sacs TJNBHZ.Sacs
The files ending in .Sacs are the isolated
fundamental
mode. The sum of the traces .Sacs and .Sacr
gives the original trace. If you read the *Z.Sacs and plot
under gsac you will see
Notice
how the body waves and scatter surface waves have been removed.
In the directory GOOD run the command
do_pom *Z.Sacs
Hit the commands SelectAll, Do Pom on page 1,
then select 4 and 200 for the period limits, Shade
on, Type Rayleigh, Nray 250, and Length 4
followed by Do Pom on page 2, and wait. You will then be
presented with a editing screen to select the phase velocities. Then
save the dispersion by clicking on Exit and then Yes.
The output dispersion is in the file disp.d which is in a surf96
dispersion format. The dispersion plot with selected dispersion
is
Here are some of the dispersion points
SURF96 R C X 0 86.23 4.21690 0.24200 10.7102
SURF96 R C X 0 80.31 4.22890 0.19760 10.7777
SURF96 R C X 0 84.45 4.22890 0.22800 10.7340
SURF96 R C X 0 81.92 4.24100 0.20980 10.7623
SURF96 R C X 0 128 4.28920 0.49920 10.4922
SURF96 R C X 0 130 4.33730 0.52780 10.4752
SURF96 R C X 0 132.1 4.38550 0.55780 10.4582
SURF96 R C X 0 134.3 4.42170 0.58610 10.4410
SURF96 R C X 0 136.5 4.46990 0.61910 10.4236
SURF96 R C X 0 138.9 4.51810 0.65350 10.4064
SURF96 R C X 0 141.2 4.55420 0.68550 10.3900
You will notice that the dispersion is not perfectly uniform. This
may be due to the fact that the path from the earthquake to each
station may have slightly different propagation which violated our
assumption of circular wavefronts.
Perform this
operation for a numebr of earthquakes, from different regions, so
that you get a good sampling of the dispersion.
Last changed March 11,
2009