README file for focal mechanism determination software ampinv
ampinv ver. 4.3 (2003/08/04）
Satoshi Ide
Department of Earth and Planetary Science, University of Tokyo
Introduction
ampinv is a program that determines the seismic moment and the focal mechanism (or moment tensor) of an earthquake using wave amplitudes and polarities of P and S waves. It uses Fourier spectral amplitude which is estimated assuming an omega-square model and path attenuation of constant Q. If site amplification function is available, you use it for amplitude correction. There are several ways to solve the problem, linear, linearlized nonlinear, and grid search methods. You can select the method depending on how good data and spatial coverage you have. There are various ways to determine the weight of each data. Since this calculation use linear scale and log scale in mixed manner, there is no ‘best’ way of weighting. You can try various weighting and if the results change too much for different weighting, the data set is not good enough to solve the problem unfortunately. We can make some figures using output files of this program and Generic Mapping Tool ver. 3.4 (Wessel and Smith, 1998).
At some stages of analysis, we use somehow ad-hoc assumption, for example, automatic picking, overall corner frequency determination, and weighting. Some of these assumptions are used to compromise the mixture of linear and log scales. Spectrum amplitude is measured in log scale but seismic moment is estimated by linear inversion. Ambiguity of crustal structure affects multiplicatively while ambiguity of radiation pattern makes difference by factors. Therefore, we should use the solution realizing that it is not a unique best solution from statistical viewpoint.
Data Preparation
First, we need to prepare SAC waveform files with proper information of arrival time and polarities of P and S waves. The unit of waveform should be SI unit (m/s, m, or m/s/s). Vertical component is used for P and two horizontal components (radial and transverse) are used for SV and SH waves. We need to rotate horizontal wave records properly in advance. Although SV can be used in case vertical wavenumber of P at surface is real, the program does not change S-window not to include P wave and special care is necessary to get reliable result from SV waves. In each SAC file the followings headers should be set (OPT is potional).
real headers
DELTA: data sampling interval
B: start time of the record relative to the reference time (RT)
O: hypocentral time of the event.
T0(*): P arrival time (manually picked)
T1(*): S arrival time (manually picked)
T5(*): P arrival time (autopicked)
T6(*): S arrival time (autopicked)
RESP0(OPT): value of 1LSB
RESP1(OPT): dynamic range (bit)
RESP2(OPT): pendulum frequency (Hz)
RESP3(OPT): damping constant
RESP4(OPT): record type (-1:disp 1:acc otherwise: vel)
STLA: station latitude (degree)
STLO: station longitude (degree)
STEL: station height (m)
EVLA: event latitude
EVLO: event longitude
EVDP: event depth(km)
CMPAZ: component azimuth (degree crockwise from N)
CMPINC: component incident angle (degree from vertical U=0)
integer headers
NZYEAR: year of reference time (RT)
NZJDAY: Julius day of RT
NZHOUR: hour of RT
NZMIN: minute of RT
NZSEC: second of RT
NZMSEC: mili-second of RT
NPTS: numbers of data sample
character headers
KSTNM: station name (up to 8 chars)
KT0(OPT): P wave polarity (P, PU, PD, PX)
KT1(OPT): S wave polarity (S, H, V, HU, HD, HX, etc)
(*): optional but at least one must be set. See pick information.
Pick information (arrival times and polarities)
Pick information is read from SAC header of each waveform file. The program uses the following headers to calculate arrival times: B, O, T0, T1, T5, T6.
There are several possibilities that give different arrival times so the code uses information as the following order.
P-wave:1: T02: T5
S-wave:1: T12: T63: O+(T0-O)*1.734: O+(T5-O)*1.73
Polarity information is also supplied in SAC headers:
KT0: P wave polarity (P, PU, PD, PX)
KT1: S wave polarity (S, H, V, HU, HD, HX, etc)
If there is ‘D’ letter in these headers, the program realizes that the polarity of this phase is down. It is same for the letter ‘U’ and polarity up. If there is ‘X’, the program will not try autodetermination of polarity. Otherwise, it determines polarities using a subroutine. So far, this routine is not sophisticated and I recommend to supply manual information for reliable estimation.
Input Files
You have to prepare i_ampinv file before executing the program. The following is a typical example:
File: i_ampinv
4 1imethod iweight (fstr fdip fslp if imethod=5) (see below)
1.0 4.0tahead tlen (see below)
0.5 30.0 20freqmin freqmax nlog (see below)
1(if positive, verbose output)
4.0(Pwave velocity at stations)
/dat/etc/struct.tbl(structure file name)
/dat/etc/ssqinv(directory name where site information are stored)
Control parameters ([ ] is a recommended value)
imethod: 1 moment tensor inversion (linear), 2 moment tensor inversion without isotropic component (linear), 3 double couple inversion (non linear), 4 double couple grid search, 5 only seismic moment (in this case the DC mechanism is given in the first line as 3rd-5th parameters) [4]
iweight: 0 no weighting, 1 scale with the difference from average, 2 scale with the difference from median [2]
tahead: time length (s) before pick [1.0 for M4]
tlen: the length (s) of time window [4.0 for M4]
freqmin: minimum frequency (Hz) for corner frequency search [0.5 for M4]
freqmax: maximum frequency (Hz) for corner frequency search [30.0]
nlog: the number of corner frequency grids in each decade [20]
The structure information is used to calculate incident angle. This is standard structure file for win system (Japanese standard seismic data analysis system).
File: struct.tbl
35.500 139.500 30.0 (Reference location: NOT USED)
6 ERI1 (nlayer) (comment)
5.5 5.51 6.1 6.11 6.7 6.71 8.0 (VP)
8.2 (VP)
4.00 0.01 10.60 0.01 16.90 0.01 600.0 (Width)
5.0 100.0 100.0 30.0 (UNUSED)
This is a FORTRAN formatted file. The first line can be ignored in this analysis, but something must be there in (3F10) format. The second line is (I5,2X,A3,2f10). These are parameter is the number of layer (nlayer, actually this is smaller by one), comment, and dummy numbers not used in this analysis. The following lines have (7F10) format and you have to supply nlayer+2 P velocities and nlayer+1 layer width. Velocity is assumed to change linearly within each layer from top to bottom velocities. The last line is not used either but you need 4 real numbers.
If any, you can use site information file such as TEST.P.tbl.
File: TEST.P.tbl
-0.350000 –0.0 0.1
-0.300000 –0.5 0.1
-0.250000 –0.3 0.1
-0.200000 –0.0 0.1
-0.150000 –0.0 0.1
-0.100000 0.2 0.1
(continue)
Each line contains log frequency, log amplitude, and its error. If this information isprovided, this amplitude is subtracted from each observed spectrum amplitude.
Execution and making outputs
Toexecute ampinv, you type
% ampinv < LIST
LIST is a file that contains the list of SAC file that is used for analysis. If you have ????.U and ????.T files for upward components and transverse components, you can type
% ls ????.[UT] | ampinv
to execute ampinv.
After execution, we get the following output files:
o_mech, o_mech.d, o_mech.p, o_mech.q: everytime
(STA).(CMP).am: made for verbose output
o_grid: made in grid search
o_mech: summary of mechanism solution
1999 04 04 08 32 55 961 35.2774 141.1939 27.08 3.68 4.170 4.467 3.548 26
P -4.170 113 16 N 0.000 -153 11 T 4.170 -31 69 VR 95 6.707 2.706
EX14 1.173 3.098 -1.746 1.606 3.335 0.000 -173 30 67 32 62 102
ERR 0.000 0.000 0.000 0.000 0.000 0.000 0 0 0 0 0 0
line 1: year, month, day, hour, minute, second, millisecond, latitude, longitude, depth, Mw, moment(1), fcp, fcs, number of waveforms
line 2: magnitude, azimuth, and plunge of P, N, and T axes, variance reduction, Vp at source, stress drop
line 3: moment(2) mrr, mtt, mff, mrt, mrf, mtf, strike, dip, and rake angles for two nodal planes
line 4 error value for line 3 (if grid search not shown)
o_mech.d: amplitudes of stations
YMZ PU 335.3 57.0 0.5053E-01 0.1627E+00 1.00 0.1929E+06 0.119E+00
YMZ HX 335.3 57.0 0.2090E-01 0.1041E+00 0.01 0.2037E+05 0.633E-01
OMZ PU 275.5 57.0 0.3496E-01 0.4400E-01 1.00 0.1855E+06 0.435E-01
(continue)
each line: station name, component (with up, down, unknown), azimuth, incident angle, observed amplitude, calculated amplitude, weight, amplitude correction, toq
o_mech.p: pick information
YMZ PU Auto 0 1 Fc 3.162 Wei 1.00
YMZ HX Manual 1 0 Fc 1.122 Wei 0.01
OMZ PU Auto 0 1 Fc 4.467 Wei 1.00
(continue)
each line: station name, component (with up, down, unknown), picker (auto or manual) , picker (if auto 0), polarity (1: up, 0:unknown, -1 down), Fc for each station, weight
o_mech.q: information for Q
YMZ PU 0.2915E+02 0.1193E+00 0.5913E-02 0.994 244
YMZ HX 0.5043E+02 0.7109E-01 0.2604E-01 1.000 796
OMZ PU 0.2803E+02 0.4349E-01 0.4123E-02 1.000 644
(continue)
each line: station name, component (with up, down, unknown), arrival time, toq (arrival time over Q), error of toq, fitting information (output of fit subroutine, see numerical recipes), Q
(STA).(CMP).am:
-0.350000 -7.46216 0.874329 -7.63813
-0.300000 -7.36535 0.681458 -7.58939
-0.250000 -7.24197 0.508539 -7.52509
-0.200000 -7.12107 0.384136 -7.46039
-0.150000 -7.10696 0.382799 -7.44953
-0.100000 -7.08593 0.368236 -7.43568
-4.999999E-02 -7.05936 0.310387 -7.44698
0.000000 -7.04939 0.251931 -7.48352
(continue)
each line: log frequency, amplitude (observed), error of amplitude, noise amplitude
o_grid: result of grid search (imethod = 4)
storing strike, dip, slip, seismic moment value and residual for those values in each line
If you are using GMT, it is easy to plot focal mechanism and spectrum. Just type
% ss.plmech o_mech
This produces mech.eps file that shows the mechanism, comparing observed amplitudes and polarities with calculated values. If you add –E option as
% ss.plmech –E o_mech
It uses o_grid file to show possible range of the solution in grid search procedure.
% ss.plfreq (STA).(CMP).am
makes (STA).(CMP).eps file that shows observed spectrum and fitted omega-square model.
Programs and compilation
ampinv: main program for mechanism determination
ss.plmech: plot mechanism, polarity, and amplitude (C-shell script)
ss.plfreq: plot spectra comparison (C-shell; script)
psmt: handmaid GMT-like routine to plot MT and DC.
ampinv can be compiled using attached Makefile and g77. It uses subroutines of lapack (sub.lapack.f), minpack (sub.minpack.f), fftpack (sub.fftpack), and Numerical Recipes (sub.recipes.f). If you have those, use them instead. psmt can be compiled with gcc. sub.trv.f contains subroutines from Japanese standard phase picking tool, ‘win’. sub.filter.f is digital filter source code written by M. Saito. Source codes used for psmt except for psmt.c are those originally used by H. Kawakatsu. Most sources are written using basic grammars (I think), so most compilers may be available as well.