Moment-tensor
inversion
Seismogram Preparation

I haven't had time to check these macros yet, so be careful.

The codes use SAC formatted data so naturally the data preparation is primarily performed using SAC. An overview of the procedure is shown in the flow chart to the right. The specific way you prepare the observations also depends on the type of observations you have - you would treat three-component regional signals different than teleseismic P-waveforms.

Setting Seismogram Header Values

The header information consists of details like the station and event locations, the origin time, etc. Station information differs for each seismometer and some values like the component azimuth differ for each seismogram. Often station information is already set in the header, and you only need to set the event information, which is common to all the data from a single event. You set those values using SAC's "changeheader" command. Here's a macro to read in signals that end with the suffixes bh[z,n,e] for vertical, north, and east.

do file wild *bhz
setbb base ( BEFORE "bhz" $file )
r %base%bhz %base%bhn %base%bhe
synch
ch evla 34.944
ch evlo  -98.508
ch evdp 5
ch kevnm 98_oklahoma
ch o gmt 1998 118 14 13 00 300
wh
enddo

The remaining steps depend on the type of signals you are working with, I'll describe an example applicable to regional seismograms.

Windowing the Signals

In SAC you window signals using the cut command. This command operates when the data are read in, so first you set the cut limits, then you read in the signals, then you save (write) the signals to keep the window. For regional seismograms you can use the origin time to define one edge of the window and a group velocity specification to define the other. For example, if your data aren't too far away, a group velocity of 1.0 to 1.5 km/s will usually be adequate to capture the complete signal. You may want to keep more if you are going to use a narrow-band filter (which has a long time-domain response).

If you defined the event origin time in the SAC header (as described above) you can use the header value "o" (the letter o) to reference it in the cut command. To get the other window, you can use the marktime command (again, provided that the origin time is set). Or you can make the calculation yourself, as shown in the macro below. I recommend that you save the original data before overwriting them as done in the following macro.

do file wild *bhz
setbb base ( BEFORE "bhz" $file )
r %base%bhz %base%bhn %base%bhe
setbb c0 &1,o
evaluate to c1 &1,dist / vgmin + %c0
*
cut %c0 %c1
r %base%bhz %base%bhn %base%bhe
cut off
*
w over
enddo

Instrument Equalization

To model the signals you have to account for the effects of instrument response (amplitude, group delay or phase) which distorts the ground motion. You can either put the response into your Green's functions or remove the response from the observations. Usually with broadband data we remove the instrument response from the seismograms. The SAC command to do this is called "transfer" and it performs a frequency domain deconvolution. Deconvolutions are often unstable so it is often necessary to limit the bandwidth of the resulting displacement seismogram using a frequency-domain taper. Technically, if you apply the taper to the displacements, you should include the same taper in your synthetics, but as long as you stay away from the regions where the taper is not unity, you can ignore it, or at least most people do.

As with all digital processing it is a very good idea, if not necessary, to remove the mean and taper the signals and this is a convenient time to do that to our data so I'll include it in the macro that does the instrument equalization. For SAC, instrument response is summarized in an ASCII "polezero" file that can be written out using RDSEED, or is usually supplied with the seismograms. I usually keep the names of each component in the polezero file, for example, the polezero file for ccm.bhz is named ccm.bhz.pz. I'll assume that is the convention in the example macro to remove instrument responses. The first two variables are the low-frequency corners for the frequency taper (this is set up for a small eastern US earthquake which does not have strong long-period signals). You should try to keep as wide a bandwidth as possible so you may need to experiment with these numbers. The high frequencies are usually much less of a problem, so I set the taper values to large values, which results in no high-frequency taper.

setbb c1 0.01
setbb c2 0.02
*
do file wild *bhz
setbb base ( BEFORE ".bhz" $file )
r %base%.bhz
rmean
taper w 0.1
trans from polezero subtype %base%.bhz.pz to none freq %c1 %c2 200 201
w %base%_dsp.z
*
r %base%.bhn
rmean
taper w 0.1
trans from polezero subtype %base%.bhn.pz to none freq %c1 %c2 200 201
w %base%_dsp.n
*
r %base%.bhe
rmean
taper w 0.1
trans from polezero subtype %base%.bhe.pz to none freq %c1 %c2 200 201
w %base%_dsp.e
*
enddo

The above macro assumes that you have three-component observations for each vertical component and will halt execution if it encounters a read error.

Rotating the Signals

Once you have the displacements you can convert the two horizontal components into radial and transverse signals. This might be a good place to generate a plot of the signals too, so I'll include that in the macro.

do file wild *_dsp.z
setbb base ( BEFORE "_dsp.z" $file )
*
r %base%_dsp.n %base%_dsp.e
rotate
w %base%_dsp.r %base%_dsp.t
*
r %base%_dsp.r
ch kcmpnm Radial
wh
r %base%_dsp.t
ch kcmpnm Trnsvrs
wh
r %base%_dsp.z
ch kcmpnm Vertical
wh
*
r %base%_dsp.z %base%_dsp.r %base%_dsp.t
p1
*
pause
*
enddo

The macro pauses between stations so that you can inspect the signals.

Marking the First Arrival Time

For the inversions I align the observed and synthetic seismograms using the first arrival so the the last step in the data preparation is to mark the first arrival time. For regional data I often find it convenient to convolve the signals with a short-period WWSSN seismometer response before picking the time. This filters away the microseism noise and generally cleans up the signal. Here's a macro to help automate the procedure. In SAC you can pick times using the PlotPK command and the "a" option. I pick the time on the vertical seismogram and then copy the value into the horizontals.

The Macro

Comments

do file wild *_dsp.z
*
setbb base ( BEFORE "_dsp.z" $file )
*
evaluate to c0  &1,dist / 8.0 + 5.0
evaluate to c1 %c0 + 5.0
evaluate to c2 %c1 + &1,o
ch t8 %c2
ch kt8 "Pn"
* 
r %base%_dsp.z
trans from none to wwsp
*
ppk
setbb a0 &1,a
*
r %base%bh? %base%_dsp.?
ch a %a0
wh
*
enddo
Loop over vertical components
         
Build a "base" name
         
Compute a rough estimate of the 
Pn arrival time for a
10 km deep source (40 km crust). 
This is not saved in the data header.
         
Read and install a WWSP response
         
         
Pick and store the time in a blackboard
variable.
         
Read in the displacements and store
the time pick in the header variable "a".
         
         
end the loop

This macro requires interaction - for each station it will read in the vertical and drop you into the "ppk" command. You pick the time using the cursor and the "a" option, then hit a "q" and the macro reads in the three components and sets the header values before moving on to the next vertical component.

Preparing for An Inversion - Alignment and Sample Rates

That completes the steps that I perform to prepare the displacement seismograms but there are several more steps that I do before the inversion. First I back up these signals so that I maintain a set of broadband displacement seismograms. Then I cut the the observations and the Green's functions so that they begin at the same time relative to the first arrival. For example, if you followed the procedure outlined above then the SAC header value "a" stores the first arrival time and the commands

cut a -2 e +0
r *_dsp.?
w over

will cut the data two seconds before the first arrival. Apply the same cut window to the Green's functions and you will have the synthetics and the data aligned for the inversion.

The last step is to decimate the observations so that they have the same sample as the Green's functions. For example, if the sythnetic seismograms are computer with a sample rate of 1 sample per second, and the original data had a sample rate of 20 samples per second, then the commands

cut off
r *_dsp.?
decimate 4; decimate 5
w over

would reduce the data rate to that of the synthetic seismogram rate.


Overview | Seismogram Prep | "Green's Functions" | Input File Format
Program Output | References | Back to Ammon's Home

Prepared by: Charles J. Ammon

Made on a Mac