Introduction

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.

Data Set

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

cd Sac
for i in *Sac
do
saccvt -I < $i > tmp
mv tmp $i
done

Now you can work with the data

Data Preparation

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

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> map r on
r on
Execute using the command: sh map.sh
GSAC> sh map.sh
sh map.sh
pscoast: Working on block #   301
pscoast: Adding Borders...
GSAC>   

The file map.eps is created and you can view it using gs, display or other PostScript viewers. You will see the following: map.png

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.

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> 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:

P001.png

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

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> 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
p002.png

Isolation of Fundamental Mode

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:
mft1.png

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

No 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   ULLBHZ.Sac
CHUBHZ.Sacr DOIT* KANBHZ.Sacs mft96.disp PUSBHZ.Sacr SOGBHZ.Sacr TEJBHZ.Sacr ULLBHZ.Sacr
CHUBHZ.Sacs HDBBHZ.Sac KWJBHZ.Sac mft96.ods PUSBHZ.Sacs SOGBHZ.Sacs TEJBHZ.Sacs ULLBHZ.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
P003.png
Notice how the body waves and scatter surface waves have been removed.

Phase velocity analysis

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

pom.png




 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