inversion |
|
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 ValuesThe 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 SignalsIn 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 EqualizationTo 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 SignalsOnce 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 TimeFor 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 RatesThat 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 |