Next: 21.5 PLOTSPEC Up: 21. Automatic routines in SEISAN Previous: 21.3 AUTOSIG Contents Index
Program to make automatic readings of WA amplitudes and make spectral fitting for channels for one or several events in an S-file. This is a simplified version of AUTOSIG and the program does not require any parameter file. The program will, by default, process Z channels for which the stations has a P or S reading and the window following the S-phase will be used. Optionally P can be used for spectra and N and E channels can also be used. A check is made for P-window length so if the window will include the S-waves, no spectrum is made. If only a P-phase available, the S-phase time will be calculated from the P-phase arrival time. If the epicentral distance is more than 500 km, the S-phase is assumed to be the Lg phase and an Lg velocity of 3.5 km is used. However if the depth is more then 50 km, S is assumed. In any case, if an S is present, the observed S-time is used.
The program can be used stand alone or in EEV with command automag(am) or command ami (input of parameters).
***************** The agnecy used is taken from SEISAN.DEF *******************************
A special option is to do a grid search for the best Q-parmeters that fit the Brune spectrum using one or several events.
If the S-file has distances, both Ml and Mw can be calcculated. The parameter for Ml is taken from STATION0.HYP in local directory or in DAT. If not found, the standard Hutton and Boore(1987) values are used. The additional filter limits for Ml (if any) are taken from MULPLT.DEF, in local directory or DAT. If not filters in MULPLT.DEF, no filters are set. If a filter is used, the amplitudes read are corrected for the gain of the filter at the determined peiod.
The average Ml and Mw are shown on screen after one run. In order to get the average Ml and Mw in the S-file after using AUTOMAG, an update must be made. However each value is printed out and can also be shown with program PLOTSPEC (see later).
The program has a lot of output on screen when running. This can mostly be used to check why a particualr spectrum or amplitude operation failed.
Program input
Interactive input: Give s-file name, spec window and wa window, if zero window length, the corresponding operation is not done, give frequecy range for spectral fit, give channel names. If station is called ALL, all Z-channels from all stations with readings (P or S) are used. Alterentively other components can be used for all stations. If not using option ALL, the user must give each channel in which case any channel can be used. Give if P or S-spectrum.
If the spectral window is given, data outside this window will not be used. However, the signal to noise ratio is still check so the actual spactral range might be smaller than the spcified range.
Input from prompt, options are:
S-file name: If only S-file name, all stations are used and default values of window lengths are applied.
s xxx : where xxx is spectral window
w xxx : where xxx is wa window
p : P-spectrum (S is default)
n : use north-south component
e : use east-west component
overwrite : If given, input file will be overwritten, but automag.out also made. Option can only be used with one event. Used with EEV.
grid : Grid search, see below
Examples, input file is sfile:
Automag sfile s 10 : Makes S-spectra only, with 10 s window length
Automag sfile s 10 p : Makes P-spectra only, with 10 s window length
Automag sfile w 30 : Makes wa operation only with window length of 50
Automag sfile s 10 w 30 : Do both
Automag sfile s 10 w 30 n : Do both but use NS channel
Automag sfile s 10 w 30 overwrite : Do both and overwrite sfile
Automag grid : Do grid search, other options cannot be used, see below
Channels to process: If interactive option, each channel must be giveni o rthe ALL option used. If non interactive: Stations with P and/or S readings are selected if station has a reading on one of the channels.
Before running AUTOMAG, it might be an advantage to remeove all old amplitude and spectral values from the input file. This can be done with program DELS or command DELS in EEV.
Location of waveform data: AUTOMAG will look for waveform files in the standard SEISAN way: Local directory, WAV, default data base under WAV (file names must be standard SEISAN) or in any other directory specified in SEISAN.DEF. NOTE that to use another WAV data base than the the default data base when AUTOMAG is running outside EEV, the data base name must be given in SEISAN.DEF, paramter WAVEFORM_BASE.
Program flow:
Read RESET TEST value for Ml in STATION0.HYP in local directory or DAT
Read filter values for Ml
Read S-file with readings and location
Get waveform files from S-file
Enter or find station and component for instrument correction
Find S/Lg-time or P-time from readings
Select out a time window from waveform channel starting a bit before the P or S-time (see below)
Find channel number in waveform file(s) corresponding to desired channel
Read the P/S-time window from waveform channel, default 50 s, starting 4 s before P/S. If S/Lg has been calcualted from P, extra time is allowed before S calculated as distance(km)/100.
Read response file
Prepare response removal, different filters and poles and zeros used for Wood-Anderson response
Correct for instrument response
Automatically pick maximum amplitude and corresponding period. The amplitude is checked for s/n by taking the spectral level around the pick period and time in the S-signal and comparing it to the spectral level in same frequency band at the start of the trace. s/n must be at least 2.
Correct for additional filter if chosen
Write out results in internal s-file, overwrite old results
Do the automatic spectral fitting using a time window of 20 s starting 4 s before the P/S. Grid search is used
Write results in internal S-file, overwrite old results
Write results in automag.out
Write channels processed in automag.list, used for plotting with plotspec
Grid search
This option allows the spectral fitting to do a grid search for the best fitting Q0, Qalpha and kappa. The user will be asked for input file name, grid parameters, spectral window, frequency range and stations. For each group of events in input file, the average residual will be calculated and the results output in file automag_grid.out in order of increasing average residual, see example below. The parameter used in one run is output in automag.par, see example later. The fit can be made with many traces from many events; up to 4000 spectra can be used for each attenuation combination in the grid search. Both P and S-waves can be used. In order to not include very bad results, the results from spectra resulting in a stress drop not in the range 1-150 bar, will not be included in the average. This means that the only data where sterss drop can be calculated will be used and the events should therefore be updated before using automag.
Q0, qalpha and kappa will play off each other, therefore a comparable fit might be obtained with a different combinations and it might not be easy find the one to use. Low Q0 and high qalpha as well as high Q0 and low qalpha might be combinations giving similar fit. However, the lower Q0 will generally give a higher spectral level and therefore a higher magnitude. Similarly, kappa also plays off with Q0, but to a lesser degree. Hence, for any given Q0 and qalpha, there is generally a best kappa in the range 0.01 to 0.05. In order to limit number of parameters, an initial value of kappa= 0.02 can be used. For the qalpha parameter, the best fits are usually obtained for qalpha in the range 0.3 to 0.9. For Q0, the values with a good fit are often in the range 100-400, depending on qalpha. A good fit to the Brune spectrum does not automatically indicate correct attenuation parameters, since the spectrum could have another shape than the theoretical shape, but at least it gives an indication of possible values. equency range can in most cases be default. However, sometimes the antialias filter cut off the signal in the high end (not corrected by SEISAN) which will give a bad fit so a high cut should be selected. Check first with MULPLT how it looks or run for one event and plot all spectra with PLOTSPEC.
The frequecy below 1 Hz at which Q becomes constant cannot be searched for and is a fixed parameter set in SEISAN.DEF. The default value is 1.0 Hz.
There is thus many combinations that give a similar residual and it might be difficult to chose. One option is to average all Q-results with the lowest residual. This can be done with the AVQ program. The program will read all Q-results from automag_grid.out and average the results over a given user specified number of results. This process is then done for all the results in a running average so it is possible to get an idea of how the parameters change with increasing residual. For details on how AVQ calculates the average, see the AVQ program under section CODAQ. An example of running the program is shown below. It is important that curves that are averaged using the same kappa since the averaging does not take kappa into account. AVQ will also make a plot showing the Q-relations for the Q-curves averaged (the first avarage) as well as the average curve.
Example of running the AUTOMAG program:
automag grid
Give input S-file
dels.out
Give spec window, enter for def (20)
25
Give spectral range, enter for 0.05-srate/2.5
give max distance, enter for all
Start q0, step Q0 and number of steps
50 50 10
Start qalpha, step qalpha and number of steps
0.3 0.1 9
Start kappa, step kappa and number of steps
0.03 0.0 1
Number of tests is: 90
Continue(n/y=enter)
Station code and component (e.g. BER BHZ), end with blank
Station code and component (e.g. BER BHZ), end with blank Giving station code ALL : Use all Z channels with readings Giving station code ALL C : Use all C channels with readings
ALL
P or S-spectrum=default (enter)
Output file for the above input:
q0= 50.0 qa= 0.500 ka= 0.030 nf= 5 re= 0.081 mw= 3.46 f0= 4.19 st= 42.0 q0= 50.0 qa= 0.600 ka= 0.030 nf= 8 re= 0.092 mw= 3.41 f0= 4.18 st= 31.8 q0= 100.0 qa= 0.900 ka= 0.030 nf= 50 re= 0.095 mw= 3.24 f0= 5.11 st= 26.4 q0= 150.0 qa= 0.900 ka= 0.030 nf= 49 re= 0.095 mw= 3.17 f0= 4.01 st= 9.3 q0= 100.0 qa= 1.000 ka= 0.030 nf= 51 re= 0.096 mw= 3.25 f0= 4.18 st= 13.1 q0= 50.0 qa= 0.300 ka= 0.030 nf= 4 re= 0.096 mw= 3.48 f0= 5.80 st= 45.3 q0= 200.0 qa= 0.800 ka= 0.030 nf= 47 re= 0.097 mw= 3.11 f0= 4.23 st= 9.0 q0= 300.0 qa= 1.100 ka= 0.030 nf= 33 re= 0.097 mw= 3.17 f0= 3.42 st= 5.4 q0= 400.0 qa= 1.000 ka= 0.030 nf= 33 re= 0.097 mw= 3.15 f0= 3.51 st= 5.4 q0= 500.0 qa= 0.900 ka= 0.030 nf= 33 re= 0.097 mw= 3.14 f0= 3.58 st= 5.5 q0= 150.0 qa= 0.800 ka= 0.030 nf= 49 re= 0.097 mw= 3.16 f0= 4.90 st= 16.4 q0= 100.0 qa= 0.800 ka= 0.030 nf= 43 re= 0.097 mw= 3.22 f0= 5.94 st= 34.9 q0= 150.0 qa= 0.600 ka= 0.030 nf= 43 re= 0.097 mw= 3.20 f0= 5.82 st= 35.6 q0= 350.0 qa= 1.100 ka= 0.030 nf= 33 re= 0.097 mw= 3.16 f0= 3.38 st= 5.2 q0= 450.0 qa= 1.000 ka= 0.030 nf= 33 re= 0.097 mw= 3.15 f0= 3.47 st= 5.2 q0= 50.0 qa= 1.000 ka= 0.030 nf= 37 re= 0.097 mw= 3.35 f0= 6.22 st= 55.2 q0= 300.0 qa= 0.800 ka= 0.030 nf= 41 re= 0.097 mw= 3.12 f0= 3.71 st= 6.3 .... .... q0= 400.0 qa= 0.400 ka= 0.030 nf= 45 re= 0.103 mw= 3.05 f0= 4.36 st= 9.0 q0= 400.0 qa= 0.500 ka= 0.030 nf= 42 re= 0.103 mw= 3.08 f0= 4.15 st= 8.0 q0= 350.0 qa= 0.600 ka= 0.030 nf= 42 re= 0.103 mw= 3.09 f0= 4.27 st= 8.3 q0= 500.0 qa= 0.500 ka= 0.030 nf= 41 re= 0.104 mw= 3.09 f0= 3.97 st= 6.9 q0= 450.0 qa= 0.500 ka= 0.030 nf= 42 re= 0.104 mw= 3.08 f0= 4.17 st= 7.6 q0= 150.0 qa= 0.300 ka= 0.030 nf= 24 re= 0.104 mw= 3.34 f0= 4.19 st= 21.1 q0= 500.0 qa= 0.400 ka= 0.030 nf= 41 re= 0.105 mw= 3.08 f0= 4.12 st= 7.6 q0= 200.0 qa= 0.300 ka= 0.030 nf= 37 re= 0.106 mw= 3.18 f0= 5.38 st= 29.6 q0= 100.0 qa= 0.400 ka= 0.030 nf= 19 re= 0.108 mw= 3.34 f0= 4.16 st= 27.0 q0= 100.0 qa= 0.300 ka= 0.030 nf= 13 re= 0.109 mw= 3.32 f0= 4.54 st= 28.4 q0= 100.0 qa= 0.500 ka= 0.030 nf= 20 re= 0.112 mw= 3.33 f0= 4.69 st= 30.1 q0= 50.0 qa= 0.800 ka= 0.030 nf= 13 re= 0.112 mw= 3.26 f0= 6.09 st= 58.0 q0= 100.0 qa= 0.600 ka= 0.030 nf= 25 re= 0.112 mw= 3.26 f0= 5.31 st= 32.5 q0= 50.0 qa= 0.400 ka= 0.030 nf= 4 re= 0.114 mw= 3.47 f0= 6.04 st= 51.2 nf is the number of fits for each test, mw is average Mw, f0 is average corner frequency and st the average stress drop.
How to find the best Q-parmeters, use of program AVQ
avq File name, enter for automag_grid.out Min number of fits to use in average, enter for 1 q0= 50.0 qa= 0.500 ka= 0.030 nf= 5 re= 0.081 mw= 3.46 f0= 4.19 st= 42.0 q0= 50.0 qa= 0.600 ka= 0.030 nf= 8 re= 0.092 mw= 3.41 f0= 4.18 st= 31.8 q0= 100.0 qa= 0.900 ka= 0.030 nf= 50 re= 0.095 mw= 3.24 f0= 5.11 st= 26.4 q0= 150.0 qa= 0.900 ka= 0.030 nf= 49 re= 0.095 mw= 3.17 f0= 4.01 st= 9.3 q0= 100.0 qa= 1.000 ka= 0.030 nf= 51 re= 0.096 mw= 3.25 f0= 4.18 st= 13.1 q0= 50.0 qa= 0.300 ka= 0.030 nf= 4 re= 0.096 mw= 3.48 f0= 5.80 st= 45.3 q0= 200.0 qa= 0.800 ka= 0.030 nf= 47 re= 0.097 mw= 3.11 f0= 4.23 st= 9.0 q0= 300.0 qa= 1.100 ka= 0.030 nf= 33 re= 0.097 mw= 3.17 f0= 3.42 st= 5.4 q0= 400.0 qa= 1.000 ka= 0.030 nf= 33 re= 0.097 mw= 3.15 f0= 3.51 st= 5.4 q0= 500.0 qa= 0.900 ka= 0.030 nf= 33 re= 0.097 mw= 3.14 f0= 3.58 st= 5.5 ... ... q0= 500.0 qa= 0.400 ka= 0.030 nf= 41 re= 0.105 mw= 3.08 f0= 4.12 st= 7.6 q0= 200.0 qa= 0.300 ka= 0.030 nf= 37 re= 0.106 mw= 3.18 f0= 5.38 st= 29.6 q0= 100.0 qa= 0.400 ka= 0.030 nf= 19 re= 0.108 mw= 3.34 f0= 4.16 st= 27.0 q0= 100.0 qa= 0.300 ka= 0.030 nf= 13 re= 0.109 mw= 3.32 f0= 4.54 st= 28.4 q0= 100.0 qa= 0.500 ka= 0.030 nf= 20 re= 0.112 mw= 3.33 f0= 4.69 st= 30.1 q0= 50.0 qa= 0.800 ka= 0.030 nf= 13 re= 0.112 mw= 3.26 f0= 6.09 st= 58.0 q0= 100.0 qa= 0.600 ka= 0.030 nf= 25 re= 0.112 mw= 3.26 f0= 5.31 st= 32.5 q0= 50.0 qa= 0.400 ka= 0.030 nf= 4 re= 0.114 mw= 3.47 f0= 6.04 st= 51.2 Number of curves to average: 90 Running average over how many, enter for average of all? 20 Q0= 155.1 qalpha= 0.83 res= 0.096 mw= 3.22 f0= 4.31 avst= 18.4 Q0= 161.2 qalpha= 0.85 res= 0.097 mw= 3.20 f0= 4.28 avst= 16.6 Q0= 164.2 qalpha= 0.87 res= 0.097 mw= 3.19 f0= 4.26 avst= 15.4 Q0= 176.7 qalpha= 0.86 res= 0.097 mw= 3.18 f0= 4.18 avst= 14.4 Q0= 183.2 qalpha= 0.86 res= 0.097 mw= 3.18 f0= 4.16 avst= 14.2 Q0= 196.3 qalpha= 0.81 res= 0.097 mw= 3.18 f0= 4.22 avst= 14.7 Q0= 202.8 qalpha= 0.84 res= 0.097 mw= 3.16 f0= 4.10 avst= 12.7 ... ... Q0= 243.4 qalpha= 0.55 res= 0.102 mw= 3.13 f0= 4.52 avst= 14.2 Q0= 243.2 qalpha= 0.54 res= 0.102 mw= 3.13 f0= 4.54 avst= 14.2 Q0= 238.8 qalpha= 0.53 res= 0.102 mw= 3.14 f0= 4.50 avst= 14.4 Q0= 241.0 qalpha= 0.54 res= 0.102 mw= 3.13 f0= 4.47 avst= 14.0 Q0= 234.1 qalpha= 0.53 res= 0.103 mw= 3.14 f0= 4.53 avst= 15.0 Q0= 222.7 qalpha= 0.52 res= 0.103 mw= 3.15 f0= 4.52 avst= 15.9 Q0= 216.9 qalpha= 0.53 res= 0.103 mw= 3.16 f0= 4.49 avst= 16.2 Q0= 206.6 qalpha= 0.52 res= 0.104 mw= 3.17 f0= 4.51 avst= 17.3 Q0= 192.4 qalpha= 0.55 res= 0.104 mw= 3.18 f0= 4.56 avst= 19.2 Q0= 182.1 qalpha= 0.56 res= 0.105 mw= 3.18 f0= 4.58 avst= 20.0 Q0= 186.4 qalpha= 0.54 res= 0.106 mw= 3.20 f0= 4.59 avst= 20.8
It is seen that the average parameters, including residual, mw and corner frequecy now is a slowly varying function of the average residual. It is also seen that a high Q0 and a high qalpha give a better fit than a low q0 and a low qalpha. If some test have few good fits, they can be eliminated by usin the parameter " Min number of fits to use in average". In this example, 20 could have been used.
Example of the output automag.par file Ml constants: 1.000000 1.11 Component used: Z Spectrum type: S Spec dur. wa. window: 20.0 50.0 Filter: 0.05 0.00 Max distance: 10000.0 Start Q0, delQ, nQ 100.0 50.0 15 Start Qa, delQa, nQa 0.10 0.10 11 Start ka, delka, nka 0.030 0.000 1
Hardwired parameters set in start of program
Default WA window: 50 s
Start 4 s before the S-time
Initial frequency range for spectral analysis: 0.05 to (sample rate)/2.5
WA period must be less than 5 s for results to be saved
Corner frequency must be larger than 0.05 for data to be used
The residual of the spectral fit must be smaller than 0.3 for spectral fit to be used
The Lg-velocy is 3.5 km/s
When doing grid search, the stress drop must be in the range 1-150 for the event to be used.
Parameters from parameter files
Q, kappa, density and velocity come from the spec model in SEISAN.DEF. If no spec model available, the model is taken from MULPLT.DEF. Notice that the frequecy dependence of Q from frequency lower than 1 Hz can only be set with the spectral model parameters in in SEISAN.DEF. See also MULPLT and SEISAN.DFF. if the parameter form MULPLT.DEF are used, Q is assumed constant=q0 below 1 Hz.
For explantion of the output of spectral parameters in the S-file, see MULPLT, spectral section.
The Ml constant are taken form STATION0.HYP
The Ml filter constants are taken from MULPLT.DEF
The spectra for one event can be plotted after running AUTOMAG, see program PLOSPEC.