How to#

Installation#

Installation of relMT is easy. The basic installation requires:

  • The FFTW library (often provided via the system package manager)

  • The Python packages NumPy, SciPy, PyYAML

Note

relMT is currently tested only on Linux systems.

Windows users, please use the Linux Subsystem for Windows (WSL) and proceed with the instructions below.

Mac users will need to install the fortran compiler gfortran on your system. Please consult the Fortran documentation for details. Then proceed below.

Pre-requisite#

This installation instruction assumes installation with Conda, pip and git. Experienced users may use uv or other tools instead.

As a pre-requisite, please install Miniconda.

Based on Conda and Pip#

We recommend to create a conda environment or to install into an existing one. Choose any Python version greater or equal 3.10

# Create and activate the environment
conda create -n relmt python
conda activate relmt

# Install the Fastest Fourier Transform in the West
conda install -c conda-forge fftw

If not already present, install pip

conda install pip

If not already present, install git

conda install git

Now install relMT locally

git clone https://github.com/wasjabloch/relMT
cd relMT
pip install .

Additional dependencies#

For plotting, we require:

  • Matplotlib for all plotting

  • networkx to visualize connections of equations in the linear system

  • Pyrocko to plot moment tensors

Consider installing these packages using the plot optional dependency:

pip install .[plot]

Some additional functionality requires community packages:

  • Import of waveforms and station inventories via ObsPy

  • Computation of spectra with Multitaper

  • Conversion to and from Cartesian coordinates with UTM

Consider installing these packages using the extra optional dependency:

pip install .[extra]

If you are working in IPython, or Jupyter, install the package in the same Conda environment to avoid version conflicts

conda install ipython

or

conda install jupyter

If you consider contributing to relMT, please install the development version

pip install .[dev]

Synopsis of processing steps#

This guide gives a brief overview over the steps necessary to compute relative moment tensors.

The algorithm is broken down in 5 steps:

  1. Initialize a project directory (relmt init) and create the input files

  2. Iteratively align, exclude and re-align event waveforms (relmt align and relmt exclude)

  3. Determine optimal bandpass windows and measure relative amplitudes (relmt amplitude)

  4. Admit a subset of amplitudes into the linear system of equations (relmt admit)

  5. Create the linear system of equations with a reference MT, apply weighting, and solve the equation system (relmt solve)

See the detailed description below to follow them or try out an example.

Set up a project#

1. Initialize directory#

In the folder where the files of your new project should reside, call from a terminal

relmt init

Alternatively, to create a new, empty project folder, e.g. myproject/, call

relmt init myproject

This creates the following structure of directories and template files:

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
|   +-- default-hdr.yaml
+-- align1/
+-- align2/
+-- amplitude/
+-- result/

2. Create the additional text files#

relMT needs to know, where the seismic stations and events are located, at what angle the seismic rays take off and which are the values of the reference moment tensors. This information is stored in the station, event, phase and reference MT files in the data/ subfolder:

data/
+-- stations.txt
+-- events.txt
+-- phases.txt
+-- reference_mt.txt

The file names are arbitrary and must correspond to the respective entries in config.yaml:

config.yaml#
# Path to the seismic event catalog, e.g. 'data/events.txt'
# (str)
event_file:

# Path to the station location file, e.g. 'data/stations.txt'
# (str)
station_file:

# Path to the phase file, e.g. 'data/phases.txt'
# (str)
phase_file:

# Path to the reference moment tensor file, e.g. 'data/reference_mt.txt'
# (str)
reference_mt_file:

The files obey a simple, whitespace-separated text file format. For details, see:

Tip

The following functions may be useful when creating the text files from external resources.

3. Create the waveform files#

For each station and both (P, S) phases, gather all event waveforms and store them as a 3-dimensional NumPy array as waveform files. Note that the approximate wave arrival (“pick”) is assumed at the center sample.

For each waveform file, populate a corresponding header file with at least the following attributes

STATION_PHASE-hdr.yaml#
# Station code
station:

# Seismic phase type to consider ('P' or 'S')
phase:

# Event numbers corresponding to the first dimension of the waveform array.
events_:

The events_ parameter is a list of integer numbers. The position in the list corresponds to the position of the waveform along the first dimentsion of the waveform array, while the value corresponds to the event number (first row) in the event file.

Default parameters that are equal for multiple stations and phases may be declared only once in default-hdr.yaml. Any values found in a specific STATION_PHASE-hdr.yaml will overwrite the values defined here:

default-hdr.yaml#
# One-character component names ordered as in the waveform array, as one string
# (e.g. 'ZNE')
# (str)
components:

# Sampling rate of the seismic waveform (Hertz)
# (float)
sampling_rate:

# Time window symmetric about the phase pick (i.e. pick is near the central
# sample) (seconds)
# (float)
data_window:

# Start of the phase window before the arrival time pick (negative seconds before
# pick).
# (float)
phase_start:

# End of the phase window after the arrival time pick (seconds after pick).
# (float)
phase_end:

# Combined length of taper that is applied at both ends beyond the phase window.
# (seconds)
# (float)
taper_length:

# Common high-pass filter corner of the waveform (Hertz)
# (float)
highpass:

# Common low-pass filter corner of the waveform (Hertz)
# (float)
lowpass:

Example

A waveform array containing P wavetrains recorded at station “BSTA” of events 1, 3, 7, 11 and 4 (in that order), recorded on three chanels with the names “Z”, “N”, “E” (in that order) with 500 samples (e.g. 5 seconds length with 100 samples per seconds) would have a shape of (5, 3, 500). The header file would have the fields:

BSTA_P-hdr.yaml#
station: BSTA
phase: P
components: "ZNE"
sampling_rate: 100
data_window: 5
events_: [1, 3, 7, 11, 4]

Tip

To create waveform arrays from an ObsPy Stream, have a look at:

4. Example data structure#

For the case of three stations ‘ASTA’, ‘BSTA’, and ‘CSTA’, all of which have P- and S-wave observations, the resulting file structure would look like this:

data/
+-- stations.txt
+-- events.txt
+-- phases.txt
+-- reference_mt.txt
+-- default-hdr.yaml
+-- ASTA_P-hdr.yaml
+-- ASTA_P-wvarr.npy
+-- ASTA_S-hdr.yaml
+-- ASTA_S-wvarr.npy
+-- BSTA_P-hdr.yaml
+-- BSTA_P-wvarr.npy
+-- BSTA_S-hdr.yaml
+-- BSTA_S-wvarr.npy
+-- CSTA_P-hdr.yaml
+-- CSTA_P-wvarr.npy
+-- CSTA_S-hdr.yaml
+-- CSTA_S-wvarr.npy

Congratulations

You are now ready to align the waveforms!

Align waveforms#

1. Prerequisites#

We assume a project has been set up. For simplicity, we here only align P and S on one station, ‘ASTA’. The corresponding minimal project directory looks like this:

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- stations.txt
    +-- default-hdr.yaml
    +-- ASTA_P-hdr.yaml
    +-- ASTA_P-wvarr.npy
    +-- ASTA_S-hdr.yaml
    +-- ASTA_S-wvarr.npy
+-- align1/
+-- align2/

Note that the event file, phase file, and reference MT file are not required at this point.

2. Synopsis#

The aim is to have waveforms aligned to sub-sample accuracy, so that as much energy as possible is projected onto the first principal component of the waveform matrix for P-waves, or the first two principal components for S-waves.

The objective function to be minimized for P-waves is:

\[\epsilon_P = 1 - C_{ij}^2\]

and for S-waves

\[\epsilon_S = 1 + 2 C_{ij} C_{jk} C_{ik} - C_{ij}^2 - C_{jk}^2 - C_{ik}^2\]

where the indices \({i, j, k}\) indicate the event indices of the combined pairs (for P-waves) and triplets (for S-waves). The objective functions imply that a polarity reversal between any two waveforms is expected. For this reason, the choice of an accurate time window is critical and waveform alignment is most of the times a two-step process.

We will first exclude corrupt or noisy traces. We then choose relatively wide time windows and pre-align waveforms. We may then choose narrower time windows and different filter passbands. We then compare measures that indicate success of alignment, before we ultimately align the waveforms with high accuracy.

3. Data quality control#

For a data set of local seismicity with a magnitude of 1 to 2, with wavetrains about 1 to 2 seconds long and a pick accuracy better than 1 second, the following configuration could be a reasonable one:

data/default-hdr.yaml#
# Sampling rate of the seismic waveform (Hertz)
sampling_rate: 100

# Time window symmetric about the phase pick (i.e. pick is near the central
# sample) (seconds)
data_window: 8

# Start of the phase window before the arrival time pick (negative seconds before
# pick)
phase_start: -1

# End of the phase window after the arrival time pick (seconds after pick).
phase_end: 3

# Combined length of taper that is applied at both ends beyond the phase window.
# (seconds)
taper_length: 0.5

# Common high-pass filter corner of the waveform (Hertz)
highpass: 0.2

# Common low-pass filter corner of the waveform (Hertz)
lowpass: 10

Corrupt traces may take the form of nan values within the trace or traces that are all 0. Such traces can be excluded using relmt exclude --nodata. Very small values, e.g. 1e-3, can also be interpreted as corrupt using the null_threshold parameter.

Noisy traces can be excluded via their signal-to-noise ratio using relmt exclude --snr. The signal amplitude is meassured within the time window given by phase_start and phase_end and the noise amplitude before phase_start within the frequency band enclosed by highpass and lowpass.

To apply the no data and the too noisy criteria at once, one can define:

data/default-hdr.yaml#
# Regard absolute amplitudes at and below this value as null
null_threshold: 0.001

# Minimum allowed signal-to-noise ratio (dB) of signals for event exclusion
min_signal_noise_ratio: 0

and call on the command line, from within myproject/

shell#
relmt exclude --nodata --snr

Assuming the P-phase of event 1 on station ‘ASTA’ had a poor signal-to-noise ratio, and the S-phase of event 2 on station ‘ASTA’ was corrupt, this would generate the following entries in exclude.yaml:

exclude.yaml#
phase_auto_snr:
- 1_ASTA_P

phase_auto_nodata:
- 2_ASTA_S

These phases will be excluded from further processing.

4. Pre-alignment#

Waveforms can be aligned using multi-channel cross-correlation (MCCC) and principal component analysis (PCA). MCCC can accommodate time-shifts as long as data_window and is limited to sample accuracy, while PCA assumes that waveforms are already aligned to 1/4 wavelength and can achieve sub-sample accuracy. For pre-alignment, MCCC suffices.

shell#
relmt align --mccc

Hint

To align many stations in parallel, consider granting access to various CPUs. Some internal NumPy routines will still use all available resources, regardless of this setting.

config.yaml#
# Number of threads to use in some parallel computations
ncpu:

This call reads waveforms from data/ and writes aligned waveforms to align1/, so that the project directory now looks like this:

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- stations.txt
    +-- default-hdr.yaml
    +-- ASTA_P-hdr.yaml
    +-- ASTA_P-wvarr.npy
    +-- ASTA_S-hdr.yaml
    +-- ASTA_S-wvarr.npy
+-- align1/
    +-- ASTA_P-hdr.yaml
    +-- ASTA_P-wvarr.npy
    +-- ASTA_S-hdr.yaml
    +-- ASTA_S-wvarr.npy
+-- align2/

5. Alignment control#

One can check visually if the alignment did succeed by inspecting the aligned waveforms:

shell#
relmt plot-alignment align1/ASTA_P-wvarr.npy
relmt plot-alignment align1/ASTA_S-wvarr.npy

The resulting plot shows the applied time shifts, the waveforms, the signal-to-noise ratio, the averaged cross correlation coeffcient per trace, the cross-correlation matrix and the expansion coefficient norm per trace in an interactive window.

Based on these plots, one can find optimal filter corners and time windows for which waveforms appear similar and characteristic. Let us assume the P-wavetrain on station ‘ASTA’ turns out to be 1.5 seconds long, was shifted so that it starts 0.5 seconds after the central sample and is most pronounced in the 1-5 Hz band. We manipulate align1/ASTA_P-hdr.yaml in the following way:

align1/ASTA_P-hdr.yaml#
# Start of the phase window before the arrival time pick (negative seconds before
# pick)
phase_start: 0.5

# End of the phase window after the arrival time pick (seconds after pick).
phase_end: 2.0

# Common high-pass filter corner of the waveform (Hertz)
highpass: 1

# Common low-pass filter corner of the waveform (Hertz)
lowpass: 5

Re-running the plot function shows the plot with the changed parameters:

shell#
relmt plot-alignment align1/ASTA_P-wvarr.npy

From intuition or reasoning one can now define parameters to exclude waveforms that do not align well. Two parameters are most relevant here:

  • min_correlation: This parameter between 0 and 1 measures how well on average a waveform correlates or anti-correlates with one other (P-waves) or two others (S-waves), where 1 is best and indicates perfect correlation. This parameter is optimized by the MCCC algorithm.

  • min_expansion_coefficient_norm: This parameter between 0 and 1 measures how well a waveform is reprented by the first principal seismogram (P-waves) or the first two principal seismograms (S-waves), where 1 is best and indicates perfect representation. This parameter is optimized by the PCA algotithm.

It is possible and can be meaningful to define different parameters for different stations and/or phases, for example:

align1/ASTA_P-hdr.yaml#
# Minimum allowed absolute averaged correlation coefficient of a waveform for
# event exclusion
min_correlation: 0.8
align1/ASTA_S-hdr.yaml#
# Minimum allowed norm of the principal component expansion coefficients
# contributing to the waveform reconstruction for event exclusion
min_expansion_coefficient_norm: 0.7

When exluding phases, we must now provide the -a 1 argument to indicate that we are reading from the align1/ directory. To exclude based on the average cross-correlation coefficient, we provide the --cc flag. To exclude based on the expansion coefficient norm, we provide the --ecn flag. Note that unset (i.e. None) values in the header files indicate that the criterion is not used for phase exclusion.

shell#
relmt exclude -a 1 --cc --ecn

Assuming the P-wave of event 3 had a cc value lower than 0.8 on station ‘ASTA’ and the S-wave of event 4 an ecn value below 0.7, the following new entries would be written to exclude.yaml:

exclude.yaml#
phase_auto_cc:
- 3_ASTA_P

phase_auto_ecn:
- 4_ASTA_S

One can as well exclude phases manually from further processing, e.g.:

exclude.yaml#
phase_manual:
- 5_ASTA_P

6. Re-align#

To re-align waveforms, we recommend to use the PCA algorithm as a last step. If one is confident that all waveforms are aligned to 1/4 wavelength one can only perform PCA alignment. Remember to read from align1/ using the -a 1 flag.

shell#
relmt align -a 1 --pca

In case cycle skips are still present we recomment to align using MCCC and PCA in succession. This is the default behaviour of relmt align when no method is specified.

shell#
relmt align -a 1

7. Result#

After this step, the align2/ directory is populated with the aligned waveforms.

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- stations.txt
    +-- default-hdr.yaml
    +-- ASTA_P-hdr.yaml
    +-- ASTA_P-wvarr.npy
    +-- ASTA_S-hdr.yaml
    +-- ASTA_S-wvarr.npy
+-- align1/
    +-- ASTA_P-hdr.yaml
    +-- ASTA_P-wvarr.npy
    +-- ASTA_S-hdr.yaml
    +-- ASTA_S-wvarr.npy
+-- align2/
    +-- ASTA_P-hdr.yaml
    +-- ASTA_P-wvarr.npy
    +-- ASTA_S-hdr.yaml
    +-- ASTA_S-wvarr.npy

You can inspect the results with relmt plot-alignment. When re-doing individual stations, delete the corresponding -hdr.yaml and -wvarr.npy files from the target directory. When re-doing all stations, activate the -o flag to overwrite in the target directory, e.g. relmt align -a 1 -o.

Measure amplitudes#

1. Prerequisites#

We assume that we have a directory with aligned waveforms. Following the alignment guide, waveforms are located in align2/

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- stations.txt
    +-- events.txt
    +-- ...
+-- align1/
    +-- ...
+-- align2/
    +-- ASTA_P-hdr.yaml
    +-- ASTA_P-wvarr.npy
    +-- ASTA_S-hdr.yaml
    +-- ASTA_S-wvarr.npy
    +-- BSTA_P-hdr.yaml
    +-- BSTA_P-wvarr.npy
    +-- BSTA_S-hdr.yaml
    +-- BSTA_S-wvarr.npy

2. Choosing a set of amplitude options#

In the configuration file, find the “Amplitude parameters” block that lists the options that control how amplitudes are measured. We start by measuring amplitudes within the frequency band defined in the highpass and lowpass parameters of the individual -hdr.yaml files, in other words choosing the manually determined filter corners.

config.yaml#
# Amplitude parameters
# --------------------
#
# Suffix appended to files, naming the parameters parsed to 'amplitude'
amplitude_suffix:

# Filter method to apply for amplitude measure. One of:
# - 'manual': Use 'highpass' and 'lowpass' of the waveform header files.
# - 'auto': compute filter corners using the 'auto_' options below
amplitude_filter: manual

The auto_ options are ignored in this case so that the next relevant option is the amplitude_measure, with the choices:

  • direct means to consider each event combination (pairs for P-waves, triplets for S-waves) seperately, determine seperate optimal filter passbands and measure relative amplitudes directly on the seismograms. This option allows to bridge larger magnitude differences.

  • indirect means to consider all waveforms (the waveform matrix) jointly, determine one common filter passband for all events in the matrix, and measure relative amplitude as the ratio of the contribution of the principal seismogram to each seismogram in the matrix.

config.yaml#
# Method to meassure relative amplitudes. One of:
# - 'indirect': Estimate relative amplitude as the ratio of principal seismogram
#     contributions to each seismogram.
# - 'direct': Compare each event combination seperatly.
amplitude_measure: direct

3. Measure the amplitudes#

To measure the amplitudes, we call relmt amplitude. We indicate that the aligned waveforms are stored in the align2/ subdirectory by supplying the -a 2 option.

shell#
relmt amplitude -a 2

This creates the amplitdue files in the amplitude/ subdirectory

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- ...
+-- align1/
    +-- ...
+-- align2/
    +-- ...
+-- amplitude/
    +-- P-amplitudes.txt
    +-- S-amplitudes.txt

4. Trying a different parameter set#

You may want to try some of the options to find optimally suited filter bands for the respective event combinations. We here fill in some values that have proven useful for certrain data sets. Note that we parse a different configuration file and set amplitude_suffix, so that the alternative option set becomes apparent through naming. We here choose the suffix auto_amp, and rename the confiugration file accordingly.

config-auto_amp.yaml#
# Suffix appended to files, naming the parameters parsed to 'amplitude'
amplitude_suffix: auto_amp

# Method to meassure relative amplitudes. One of:
amplitude_measure: direct

# Filter method to apply for amplitude measure.
amplitude_filter: auto

# Method to estimate lowpass filter that eliminates the source time function. One
# of:
# - 'fixed': Use the value 'fixed_lowpass' (not implemented)
# - 'corner': Estimate from apparent corner frequency in event spectrum
# - 'duration': Filter by 1/source duration of event magnitude.
#     Requires 'auto_lowpass_stressdrop_range'
auto_lowpass_method: corner

# When estimating the lowpass frequency of an event as the corner frequency
# (auto_lowpass_method: 'corner'), assume a stressdrop within this range (Pa).
auto_lowpass_stressdrop_range: [1e5, 1e7]  # 0.1 to 10 MPa

# Include frequencies with this signal-to-noise ratio to optimal bandpass filter.
# Respects lowpass constraint. If not supplied, do not attempt to optimize
# passband.
auto_bandpass_snr_target: 0

# Minimum ratio (dB) of low- / highpass filter bandwidth in an amplitude ratio
# measurement. When positive, discard observation outside dynamic range. When
# negative, lower the highpass until the (positive) dynamic range is reached.
min_dynamic_range: 5

We use the -c option to explicitly parse the renamed configuration file.

shell#
relmt amplitude -a 2 -c config-auto_amp.yaml

Each phase observation is now assigned a bandpass, which is stored in amplitude/bandpass.yaml (with the optional amplitude_suffix applied). If you change any of the auto_ options and re-run the above command, remember to parse the -o option to overwrite the bandpass.yaml or change the amplitude_suffix. If a bandpass.yaml with the matching amplitude_suffix is found, we read from that file instead of computing new filter corners.

min_dynamic_range is only applied to waveform combinations and can therfore be changed after computing filter bands.

With the changed naming conventions, we now find an alternative set of amplitudes besides the previously computed.

myproject/
+-- config.yaml
+-- config-auto_amp.yaml
+-- ...
+-- data/
    +-- ...
+-- align1/
    +-- ...
+-- align2/
    +-- ...
+-- amplitude/
    +-- P-amplitudes.txt
    +-- S-amplitudes.txt
    +-- bandpass-auto_amp.yaml
    +-- P-amplitudes-auto_amp.txt
    +-- S-amplitudes-auto_amp.txt
    +-- ...

5. Conclusion#

We showed how to measure amplitudes in two different ways. The amplitude files may grow very large and can contain outliers. When building the linear system it is important to admit a balanced set of well-constrained amplitudes.

Admit amplitudes into the linear system#

1. Prerequisites#

The aim is to admit the “best” relative P and S amplitudes from the potentially very large sets stored in P-amplitudes.txt and S-amplitudes.txt into the linear system of equations that constraints the relative moment tensors. We assume the following file structure:

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- ...
+-- align1/
    +-- ...
+-- align2/
    +-- ...
+-- amplitude/
    +-- P-amplitudes.txt
    +-- S-amplitudes.txt

As relative S amplitudes represent comparission of event triplets the number of possible S combinations is \(\binom{N}{3}\), for \(N\) events, while the maximum number of P pairs is only \(\binom{N}{2}\).

2. Admit parameters#

The following parameters in config.yaml define the critera by which relative amplitude observations should be excluded. As before, different parameter sets are distinguished using admit_suffix, which is appended to the name of the output files.

config.yaml#
# Admission paramters
# -------------------------
#
# Suffix appended to the amplitude suffix, naming the admission parameters
# parsed to 'admit'
# (str)
admit_suffix: admitted

The first parameters pertain to the quality of waveform reconstruction:

  • max_amplitude_misfit discribes how well a waveform is reconstructed, where a value of 0 indicates perfect reconstruction, >1 indicates that the misfit amplitude is larger than the singal amplitude itself. This value should be <1.

  • max_s_amplitude_misfit As max_amplitude_misfit, but sets a different value for S-waves.

  • max_s_sigma1 discribes the linear independence of the two relative S amplitudes. The values of \(B_{abc}\) and \(B_{acb}\) may become arbitrary as this value approaches 1 and the equation may not be linearily independent of others. However, in the case of very similar S waveforms, \(\sigma_1\) may be close to 1 for almost all waveform combinations. One then runs into danger of excluding all observations when choosing too low a value.

config.yaml#
# Discard amplitude measurements with a higher misfit than this. Applies only to P
# amplitudes if 'max_s_amplitude_misfit' is given.
# (float)
max_amplitude_misfit: inf

# If given, discard S amplitude measurements with a higher misfit.
# 'max_amp_misfit' then only applies to P amplitudes.
# (float)
max_s_amplitude_misfit:

# Maximum first normalized singular value to allow for an S-wave reconstruction. A
# value of 1 indicates that S-waveform adheres to rank 1 rather than rank 2 model.
# The relative amplitudes Babc and Bacb are then not linearly independent.
# (float)
max_s_sigma1: 1.0

The next parameters pertain to properties of event combinations.

  • max_magnitude_difference: Large differences in event magnitude may cause problems in the relative amplitude measurement when the lowpass filter corner was chosen too high and the depletion in high frequency energy of the larger-magnitude event becomes significant.

  • max_event_distance: Large inter-event distances may be an indication that the common Green’s function assumption is violated.

config.yaml#
# Maximum difference in magnitude between two events to allow an amplitude
# measurement.
# (float)
max_magnitude_difference:

# Maximum allowed distance (m) between two events.
# (float)
max_event_distance:

The next parameters assure that only observations that contribute to a meaningful solution enter the system of equations.

  • min_equations indicates that each relative moment tensor must be constrained by at least this many equations

  • max_gap is the maximum azimuthal gap allowed for a moment tensor. If this value is exceeded, all ob

If one of the values is exceeded for a moment tensor, all observations pertaining to that MT will be discarded. These criteria are applied iteratively until no equations are left that violate either.

config.yaml#
# Minimum number of equations required to constrain one moment tensor
# (int)
min_equations: 1

# Minimum number of stations required to constrain one moment tensor
# (int)
min_stations: 1

# Maximum azimuthal gap allowed for one moment tensor
# (float)
max_gap: 360.0

The next paramters allow to decimate the number of S observations. There are \(\binom{3}{N}\) possible S triplets, but only \(\binom{N}{2}\) possible P pairs. However, S observations should not necessarily dominate the equation system.

  • two_s_equations determines, if each S amplitude observation contributes two equations (of the two longest polarization vectors \(g\) for the station-event configuration), or only one (of the longest polarization vector)

  • max_s_equations allows to exclude S observations that are redundant in the sense that the events contributing to the triplet have been observed many times and on stations that have many observations. Observations are exluded until the threshold is reached. When set, keep_events is a list of events that should not be considered as redundant and equation_batches controls how often the equations should be re-ranked.

config.yaml#
# Use two equations per S-amplitude observation triplet (`False` only includes the
# one with the highest norm of the polarization vector).
# (bool)
two_s_equations: True

# Maximum number of P-wave equation in the linear system. If more are available,
# iteratively discard those with redundant observations, on stations with many
# observations, and with a higher misfit
# (int)
max_p_equations:

# As `max_p_equations`, but for S equations.
# (int)
max_s_equations:

# When reducing number of equations, increase importance of these events by not
# counting them in the redundancy score. Use to keep more equations e.g. for the
# reference event or specific events of interest.
# (list)
keep_events:

# When reducing the number of S-wave equations, rank observations iteratively this
# many times by redundancy and remove the most redundant ones. A lower number is
# faster, but may result in discarding less-redundant observations.
# (int)
equation_batches: 1

3. Applying admission parameters#

When a set of parameters has been found, run

shell#
relmt admit

which will place the admitted files next to the original ones

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- ...
+-- align1/
    +-- ...
+-- align2/
    +-- ...
+-- amplitude/
    +-- P-amplitudes.txt
    +-- S-amplitudes.txt
    +-- P-amplitudes-admitted.txt
    +-- S-amplitudes-admitted.txt

Solve the linear system of equations#

1. Synopsis#

After amplitude admission, each line of the P amplitude file represents one equation of the linear system and each line of the S amplitude file either one or two (depending on the two_s_equations parameter).

We now combine the amplitudes with the take-off angles stored in the phase file and the reference MT stored in the reference MT file into a linear system of equations. We apply a weighting scheme as to counter-act large amplitude differences and honour variable data quality. We further show how to apply model constraints. Finally, we solve the system in a least square sense and show how to detect outliers trough the bootstrapping method.

2. Defining the reference MT#

As before, the parameter sets can be made distinguishable by defining result_suffix, which will be appended to the filenames created here.

config.yaml#
# Solve parameters
# ----------------
#
# Suffix appended to amplitude and admit suffices indicating the parameter set parsed
# to 'solve'
# (str)
result_suffix:

One or multiple reference MTs can be inserted into the right hand side of the equation system. The event number must correspond to an entry in the event file and the reference MT file. We usually apply a weight of 1000 to the reference MT, which forces that reference event to attain the reference MT. Note that the linear system of equations is normalized so that parameters are in the -1 to 1 range.

config.yaml#
# Event indices of the reference moment tensor(s) to use
# (list)
reference_mts: [0] # For example

# Weight of the reference moment tensor
# (float)
reference_weight: 1000

3. Applying a constraint#

A deviatoric or no constraint (i.e. full MT) can be applied to the solutions, meaning that we either solve for 5 or 6 MT elements. Note that when a deviatoric constraint is applied, we do not consider the isotropic part of the reference MT and the resulting magnitudes will be lower.

config.yaml#
# Constrain the moment tensor. 'none' or 'deviatoric'
# (str)
mt_constraint: none

4. Weighting observations by misfit#

To emphasize better observations, those with a lower misfit are given a larger weight (i.e. unity), which is assigned to all observations with a misfit lower than min_amplitude_misfit. To avoid that observation with a misfit approaching the largest allowed misfit (max_misfit, see admit amplitudes) get a weight close to 0, one can set the min_amplitude_weight.

config.yaml#
# Minimum misfit to assign a full weight of 1. Weights are scaled linearly from
# `min_amplitude misfit` = 1 to `max_amplitude_misfit` = `min_amplitude_weight`
# (float)
min_amplitude_misfit: 0.0

# Lowest weight assigned to the maximum amplitude misfit
# (float)
min_amplitude_weight: 0.0

5. Drawing bootstrap samples#

To detect outlying observation one can draw bootstrap samples. This will create an additional relative MT file with a “-boot” suffix.

config.yaml#
# Number of samples to draw for calculating uncertainties. If not given, do not
# bootstrap.
# (int)
bootstrap_samples: 0

6. Calling the solve routine#

The solution is calculated by calling.

shell#
relmt solve

which will read the amplitude files with the combined amplitude_suffix and admit_suffix defined in the config.py and write the solution to the result/ subdirectory, with the result_suffix appended.

myproject/
+-- config.yaml
+-- exclude.yaml
+-- data/
    +-- ...
+-- align1/
    +-- ...
+-- align2/
    +-- ...
+-- amplitude/
    +-- ...
+-- result/
    +-- relative_mts.txt
    +-- relative_mts-boot.txt  # Optional bootstrapping results

7. Computing synthetic amplitudes#

Amplitude predictions for the found solution are computed using the --predict argument. In that case, we also need to specify which original waveform observations the synthetic amplitudes should be compared to. That is, when the amplitudes were measured on the waveforms stored in align2, we supply the -a 2 argument

shell#
relmt solve --predict -a 2