Python Time Domain Moment Tensor User's Page
Nera Project « Back to table of contents
NERA: Network of European Research Infrastructures for Earthquake Risk Assessment and Mitigation
NERA (2010-2014) is an EC infrastructure project that integrates key research infrastructures in Europe for monitoring earthquakes and assessing their hazard and risk. NERA integrates and facilitates the use of these infrastructures and access to data for research, provides services and access to earthquake data and parameters, and hazard and risk products and tools. NERA activities are coordinated with other relevant EC-projects and European initiatives, and contribute to the ESFRI EPOS infrastructure and the OECD GEM program. Grants providing access to highly specialised research facilities are also available within NERA.
JRA2: Tools for real time seismology, acquisition and mining « Back to table of contents
Off-line tool-box for Time Domain Moment Tensor Inversion
JRA2 aim is to interface existing research analysis software with the data to the end of improving seismic real-time monitoring. This allows the testing of the software in a semi-productive mode without having to insert it within existing specialised seismic monitoring software packages such as EarthWorm or SeisComP.
Python Time Domain Moment Tensor - pyTDMT « Back to table of contents
Python Time Domain Moment Tensor - pyTDMT - inversion is a python package to invert seismic waveform data for the moment tensor. The inversion is performed into time domain. pyTDMT does not require data-preprocessing and input file formatting. The inversion, from data load to plots, is performed and controlled by arguments and does not require to run sequence of scripts. The inversion is performed using only python modules with exception of the greens computation which is performed using a single fortran code (FKRPROG.F) which is included into the package.
Download « Back to table of contents
pyTDMT package can be downloaded here
Package Requirements and installation « Back to table of contents
Package requirements
- Python2.7 or higher
- Obspy
- numpy
- scipy
- rdseed-5.0 or higher
- Fortran compiler
Note for Mac OS users: numpy and scipy are already included into obspy.dmg installing file
Installation
- 1. Unzip and untar the pyTDMT.tar.gz.
- 2. Compile the greens generator code:
cd src/fkrpog
gfortran -O3 FKRPROG.f -o FKRPROG
or
g77 -ff90 -O3 FKRPROG.f -o FKRPROG
mv FKRPROG ../../bin - 3. Add the path yourlocalpath/tdmt/bin and yourlocalpath/tdmt/py to your $PATH in your profile
Usage « Back to table of contents
>>> pytdmt.py --help or pytdmt.py -h
>>> Usage: pytdmt.py [-h] [--fseed FSEED] [--mseed MSEED] [--sac SAC] [--asc ASC]
[--Paz PAZ] [--Resp RESP] [--stsC STSC] [--deco DECO]
[--ori ORI] [--len LEN] [--pre PRE] [--fsfile FSFILE]
[--dva DVA] [--area AREA] [--epi EPI] [--title TITLE]
[--outdir OUTDIR] [--wsac WSAC] [--model MODEL] [--npts NPTS]
[--delta DELTA] [--cpus CPUS] [--rvel RVEL]
[--bandpass BANDPASS] [--highpass HIGHPASS] [--lowpass LOWPASS]
[--zeroph ZEROPH] [--taper TAPER] [--sim SIM] [--flim FLIM]
[--deci DECI] [--inter INTER] [--war WAR] [--DeltaInv DELTAINV]
[--mti MTI] [--maxw MAXW] [--range RANGE] [--purge PURGE]
[--set SET] [--azi AZI] [--nrIter NRITER] [--iso ISO] [--vr VR]
[--zcor ZCOR] [--plt PLT] [--pltname PLTNAME] [--pltndc PLTNDC]
[--sort SORT] [--verbose VERBOSE]
Options « Back to table of contents
| -h, --help | Show this help message and exit |
| --fseed FSEED | fseed file name inclusive of PATH if PATH different than \. Default=None |
| --mseed MSEED | mseed file name inclusive of PATH if PATH different than \. Default=None |
| --sac SAC | Path for sac binary files. Default=None |
| --asc ASC | Path for sac alpha files. Default=None |
| --paz PAZ | Path for 'PZs' poles and zeros (sac file format) files. Default=None |
| --stsC STSC | Station coordinates file. Mandatory without --fssed option. Default=None. |
| --deco DECO | Apply deconvolution of instrument from observed [Y]/N. |
| --ori ORI | Event Origin Time (e.g.: 2011-07-25T12:30:00). |
| --len LEN | Length of signal in seconds after event origin Time. |
| --pre PRE | Length of signal in seconds before event origin Time. |
| --dva DVA | Data type: 1(displacement); 2(velocity); Default=1 |
| --area AREA | String: earthquake epicenter area. Default="Earthquake" |
| --title TITLE | Title to diplay in the plot. Default="Event". |
| --epi EPI | Epicenter coordinate and depth: Lat Lon Depth. |
| --outdir OUTDIR | Directory for data extraction for seed input file. Default=data |
| --wsac WSAC | Write observed and Greens in binary sac file. Y/[N] |
| --model MODEL | Earth model for greenfunctions. Default="None". |
| --npts NPTS | Number of points for greens. Power of 2. Default="1024" |
| --delta DELTA | Epicenter coordinate and depth: Lat Lon Depth. |
| --cpus CPUS | Number of CPU available for greens computation. Min=1, Max=4. Default="1" |
| --rvel RVEL | Reduction velocity. Recommendet value = 8km/s. Default="8" |
| --bandpass BANDPASS | Bandpass filter "corners fimn fmax". E.g.:"2 0.01 0.1". Default="None" |
| --highpass HIGHPASS | Highpass filter "corners fimn". E.g.:"2 0.01 0.01". Default="None" |
| --lowpass LOWPASS | Lowpass filter "corners fmax". E.g.:"2 0.1 0.01". Default="None" |
| --zeroph ZEROPH | Zerophase for filters. True/[False]. |
| --taper TAPER | Costaper. If taper=-1 no taper is applied. Defaults="0.1" |
| --sim SIM | Remove instrument method: PZs for poles and zeros, RESP, for RESP_ files. Default="PZs" |
| --flim FLIM | Corner frequency for deconvolution filtering. Defaults="0.002 0.005 0.5 1" |
| --deci DECI | Decimation factor for sampling rate. Only integer decimation factor allowed. Default="None" |
| --inter INTER | Interpolate data to correct samplingrate if sampling not correct after decimation [Y]/N. Artifact may occour. Warning about decimation use --war Y. |
| --war WAR | Warnings. Y/[N]. |
| --DeltaInv DELTAINV | Delta in Hz (data and greens) for MT inversion. Default="1.0" |
| --mti MTI | Mxx Myy Mzz Mxy Mxz Myz Mo. Default="1 1 1 0 0 0 1e20" !! Unused !! in this version. |
| --maxw MAXW | Maximum variance pro station allowed. maxw=1 means perfect fit. Default="2.5". !! Unused !! |
| --range RANGE | Min and Max distance range in km for stations to use. Default="None" |
| --purge PURGE | Station list to remove from dataset. Data are unused but not deleted. Default="None" |
| --set SET | Station list to use. Default="All" |
| --nrIter NRITER | Max number of iteration allowed. nrIter=0->No iteration, =-1->auto. Default="0". See tutorial for details. |
| --iso ISO | iso=0 for isotropic component set to 0, else iso=1. Default=0 |
| --vr VR | min VR (variance reduction) pro station allowed.Default="50" |
| --zcor ZCOR | Fix Zcor for all 3 components (iso) or for each component (uniso), Default="iso" |
| --plt PLT | Plot waveformfit [Y]/N |
| --pltname PLTNAME | Plot name file. Default="plot.pdf" |
| --pltndc PLTNDC | Plot non douple couple[Y] / only douple couple [N] focal mechanism. [Y]/N. |
| --sort SORT | Sort stations with respect to C:code D:distance A:azimuth. Default="A D C" |
| --verbose VERBOSE | Verbose option. [0]/1. |
Greens and Eath models « Back to table of contents
Green's functions are generated by the frequency-wavenumber generator FKRPROG for 1D flat layered Earth Model written
by Chandan K. Saikia of URS. References for this code are: 1. DUNKIN J. W. (1965),BSSA,335-358
2. HASKELL N. A. (1964),BSSA,337-393 3. WANG AND HERRMANN (1980),BSSA,1015-1036
4. WATSON T. H. (1970),BSSA,161-166.
Some formatted earth model files are included into the package. The Included models are: 1. PREM, 2. AK135, 3. gil7,
4. Aquila (suitable for central Italy)
See section I/O files for details on the file format.
I/O Files « Back to table of contents
Earth Model
Example for the Aquila Earth Model:
| 0.1500E+01 | 3.7497E+00 | 2.1436E+00 | 2.2753E+00 | 600.00 | 300.00 |
| 0.3000E+01 | 4.9399E+00 | 2.8210E+00 | 2.4858E+00 | 600.00 | 300.00 |
| 0.3000E+01 | 6.0129E+00 | 3.4336E+00 | 2.7058E+00 | 600.00 | 300.00 |
| 0.7000E+01 | 5.5516E+00 | 3.1475E+00 | 2.6093E+00 | 600.00 | 300.00 |
| 1.5000E+01 | 5.8805E+00 | 3.3583E+00 | 2.6770E+00 | 600.00 | 300.00 |
| 0.6000E+01 | 7.1059E+00 | 4.0081E+00 | 3.0002E+00 | 600.00 | 300.00 |
| 0.8000E+01 | 7.1000E+00 | 3.9864E+00 | 3.0120E+00 | 600.00 | 300.00 |
| 1.0000E+01 | 7.9000E+00 | 4.4036E+00 | 3.2060E+00 | 600.00 | 300.00 |
Each line represent a layer of the model. Each line has the following values:
1. Thikness of the layer in [km]
2. p-velocity [km/s]
3. s-velocity [km/s]
4. Density [gr/cc]
5. Q-alpha
6. Q-beta
Station coordinate file
- Using --fseed option: This file is automatic generated from rdseed when using a fseed file as observed waveform input file. No need to write a new one!
-
Using --mseed / --sac / --asc options: You need to generate this file manulally and calling it using --stsC option.
This is a cvs file and has the following format:
StationCode Network Latitude Longitude Elevation
Name format for SAC/ASC and PZs files and PZs file format.
- Using --fseed option: SAC and ASC file are not needed. PZs files automatic generated from rdseed.
- Using --mseed option: SAC and ASC file are not needed. PZs files need to be written.
- Using --sac / --asc options: SAC and ASC files respectively should be placed into some directory. PZs files need to be written.
- PZs file name format: File name must begin with SAC_PZs_NETWORK_CODE_Channel_ (SAC_PZs_IV_SALB_HHZ_)
- SAC/ASC file name format: File name must end with .SAC
PZs file format
PZs file format follows the SAC-PZs file format:
| ZEROS | Nr_of_zeros | |
| zeros_real | zeros_imm | |
| POLES | Nr_of_poles | |
| poles_real | poles_imm | |
| CONSTANT | value |
Example:
| ZEROS | 3 | |
| +0.000000e+00 | +0.000000e+00 | |
| +0.000000e+00 | +0.000000e+00 | |
| +0.000000e+00 | +0.000000e+00 | |
| POLES | 5 | |
| -3.564700e-02 | -3.687900e-02 | |
| -3.564700e-02 | +3.687900e-02 | |
| -2.513300e+02 | +0.000000e+00 | |
| -1.310400e+02 | -4.672900e+02 | |
| -1.310400e+02 | +4.672900e+02 | |
| CONSTANT | +4.634211e+16 |
pyTDMT output
The pyTDMT generates automatically in the sub directory
inv the following output files:
- args_options.log: ASCII run log. List all aptions and their values for the last run
- green.log: ASCII log file generated by FKRPROG
- mt_inv.out: ASCII output file with the MT inversion results. The content is set refreshed at each run.
- mt_inv.log: ASCII output file with the MT inversion results. The content is happened at each run.
plots
The option --plt Y allows to generate waveformfit plots for
each inversion run. The plot opens at the end of the run and automatically saves pdf files for each plot page in subdirectory
inv. The name of the
plot file is set with the option --pltname NameOfThePlot.
For example a two pages plot will generate two files named NameOfThePlot-1.pdf
and NameOfThePlot-2.pdf.
Here an example of waveformfit plot.
The plot shows the following informations (from top to bottom):
- Number of page / Number of total pages
- The three components column: tangential, radial, vertical
- For each station and each component:
- Station Code, Epicentral distance in km, Azimuth, shift between observed and synthetics in samples, Variance reduction
- Observed seismogram (black line) and synthetic seismogram (red dashed line)
- Epicentral Area / Title
- Event origin time and location
- Planes orientation (Strike/dip/Slip)
- MTi Element values
- Scaling factor for MTi elements value, Seismic moment [dyn*cm], Mw
- Weighted reduced variance [0-100] and the quality [0-4] (from min to max).
VR∈[0,20[ -> 0; VR∈[20,40[ -> 1; VR∈[40,60[ -> 2; VR∈[60,80[ -> 3; VR∈[80,100] -> 4; - Nr of statiopn used, filter frequency range, length of the waveform window in seconds
- Beachball. using option --pltndc you may choose to plot the entire MT solution or only the douple couple part.
Performance test « Back to table of contents
We perform a quality and performance test over a catalog of 250 earthquake with magnitude 3.23 <= Mw <=6.3 occurred in the Italian Peninsula and near bordering regions from begin 2006 till end of December 2012. This catalog refers to the tdmt Moment tensor computed at the Centro Nazionale Terremoti (CNT) of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) solutions published at the following page: http://cnt.rm.ingv.it/tdmt (Scognamiglio et al. 2009). The published catalog is manually revised.
We perform a full automatic MT inversion for each event of the catalog. None of the automatic MT is manually revised, since the goal of this test is not to obtain the best solution for each event of the catalog.
The inversions applied the same parameters and earth model for the entire data-set. The full list of the parameters is listed in the table below (where not specified, we used the default parameters of the pytdmt_V2.2 release).
For all our automatic solutions we check the Mw and the orientation of the principal axis with respect the solutions of the INGV catalog
Args options:
| --beg | Origin time from CPT catalog http://cnt.rm.ingv.it |
| --ori | Epicenter coordinates and depth¹ from CPT catalog http://cnt.rm.ingv.it |
| --len | 300 |
| --pre | 50 |
| --model | pytdmt/earth-models/Aquila |
| --lowpass | 2 0.050 |
| --highpass | 2 0.020 |
| --clean | 1 |
| --noise | 0.1 |
| --range | 50 300 |
| --auto | Y² |
| --plt | N; |
¹ :Earthquakes with depth 7.5 km were moved to 8.5 Km since with the earth model used at 7.5 it was not possible to compute the Green's functions since a bug in the code FKRPROG.F. 13 events with depth greather than 50 km where not inverted since our model do not include deep layers.
² : This Option automatically enable --nrIter -1. For details on the other paramenters (Variance reduction, distance range and corner frequency for the filters) we used the default values. Type ./pytdmt.py --help to check the default values.
General overview on the solutions:
We obtain 109 solutions with MT quality larger or equal to 1 (Table 1). The minumum variance reduction imposed using --auto Y will not allow to obtain quality 0 solutions, since all low quality fit traces are removed from data-set during the iterations. The minimum magnitude solution has Mw=3.23 while the maximum magnitude has Mw = 6.3.
We did not obtain solutions for earthquake located not near to the italian peninsula (e.g.: Balkan, Serbia, ...) we did not obtained solutions since the earth-model used does not fit the crustal structure. Often, particularly for the older and smaller events of the catalog, the data obtained were sparsed and not close enought to obtain a solution with quality larger or equal to 1.
Table 1: Performances
| Year | Nr. Events | Nr. MT | No data | Too deep | Rate |
| 2012 | 52 | 45 | 7 | 0 | 85.6% |
| 2011 | 20 | 14 | 6 | 0 | 70.0% |
| 2010 | 22 | 17 | 5 | 0 | 77.3% |
| 2009 | 76 | 65 | 11 | 1 | 84.2% |
| 2008 | 21 | 10 | 7 | 4 | 47.6% |
| 2007 | 38 | 14 | 18 | 6 | 36.8% |
| 2006 | 21 | 15 | 4 | 2 | 71.4% |
| 2006-2012 | 250 | 185 | 41 | 13 | 74.0% |
Table 1: General performance of pyTDMT with respect the INGV tdmt catalog.
Mw:
For the magnitude performance test we check the automatic solutions obtained with at least 1 station for all quality levels. No significative difference on mean and standard difference is observed for solutions obtained with more stations. The table 2 below shows the mean differences between our solutions and the ones of the CNT tdmt catalog. The figure 1 shows our Mw values versus CNT tdmt Mw values.
Table 2: Mw
| Q | Nr. | Min | Max | Mean | Std |
| 1 | 33 | 3.2 | 5.1 | -0.14 | 0.21 |
| 2 | 68 | 3.3 | 4.6 | -0.04 | 0.11 |
| 3 | 72 | 3.5 | 6.3 | -0.04 | 0.10 |
| 4 | 12 | 3.6 | 6.1 | -0.05 | 0.06 |
| 1-4 | 185 | 3.2 | 6.3 | -0.06 | 0.13 |
Table 2: Mw comparison between our values obtained with at least 1 stations and the one of the CNT tdmt catalog.
Figure 1: Mw
Figure 1: Mw comparison between our Mw (vertical axe) obtained with at least 1 stations and the one of the CNT tdmt catalog (horizontal axe). Light grey diamond refers to quality 1 MT solutions; middle grey stars to quality 2 MT solutions; dark grey dots to quality 3 and 4 MT solutions.
Table 1 and figure 1 show that Mw difference are close to 0.05 -/+ 0.1. Only Solutions with quality 1 shows a larger mean difference (-0.12) and standard deviation (0.21). These large values are biased by 3 Mw=2.8-3.5 earthquakes (CNT catalog). The MT solutions was obtained using only one station with very low quality data, that the option --clean and --noise were not able to remove. Generally the our Mw values computed fully automatic using the --auto Y option, show a good agreement with the values of the CNT catalog manually revised. As expected, smaller magnitude events results into MT solutions with lower Variance reduction values (quality 1 and 2), because of the smaller signal-to-noise ratio and a less precise crustal model for the shorter highpass filter used for the inversion.
Principal axis:
For the principal axis performance test we check the automatic solutions obtained with at least 1, 2, 3 and for stations and for all quality levels. The table 3 below shows the mean of the mean angle difference of the principal axis directions (Bernardi et al. 2004) between our solutions and the ones of the CNT tdmt catalog. The figure 3 shows the mean angle difference of the principal axis with respect the number of stations and solution quality. Figure 4 shows the mean angle difference of the principal axis with respect the magnitude Mw and the solution quality.
Table 3: Principal axis
| Q | Stat. | Nr. | Min | Max | Mean | Std | Out |
| =1 | ⟩=2 | 27 | 2 | 60 | 29.9 | 21.3 | 13 |
| =2 | ⟩=2 | 62 | 2 | 64 | 19.1 | 15.6 | 12 |
| =3 | ⟩=2 | 60 | 2 | 62 | 19.2 | 13.4 | 12 |
| =4 | ⟩=2 | 2 | 8 | 26 | 15.6 | 7.2 | 0 |
Table 2: Comparison between the mean diference angle of the principal axis from our values obtained with at least 2 stations and the one of the CNT tdmt catalog.
Figure 2: Principal axis
Figure 2: Mean principal axis angle difference between our MT solutions and the ones on the INGV tdmt catalog with respect to the number of station used to obtain the final solution. Symbols and grey scale refers to figure 1.
Table 2 and figure 2 show the mean difference angle of the principal axis between our MT solutions and the one of the CNT catalog. Only solutions obtained with more than 2 stations are considered. The mean values and the standard deviations shows that generally the mean difference angle is below 30 degree. Solutions with less than 30 degree are well resolved and reliable (Bernardi et al. 2004). Figure 2 shows that generally the mean angle difference decrease increasing the number of station used for the final inversion, although well constrained solution can be obtained also using 2 or more stations.
References:
Bernardi, F., Braunmiller, J., Kradolfer, U., and Giardini, D., 2004. Automatic regional moment tensor inversion in the European-Mediterranean region, Geophys. J. Int., Vol. 156, pp. 1-14
Scognamiglio L., E. Tinti, and A. Michelini, 2009. Real-time determination of seismic moment tensor for italian region, Bull. Seism. Soc. of Am., Vol. 99, No. 4, pp. 2223-2242, doi:10.1785/0120080104.