Magnitude calculation

[mmesim]

Ideas for magnitude calculation based on what is already available and used in the literature (etc).

• SVD method
  https://pubs.geoscienceworld.org/ssa/bssa/article/100/5A/1952/325099/Precise-Estimation-of-Repeating-Earthquake-Moment

• Amplitude ratio between template and detected event
  e.g https://www.nature.com/articles/ngeo697

• Schaff & Richards, 2004.
  https://academic.oup.com/gji/article/197/1/335/682318?login=true

• Template based magnitude-amplitude regression
  https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JB017468

[damb]

Thanks for opening this discussion, @mmesim.

Previously, I listed some of the alternatives already in #51.

As already discussed, I'm focusing on https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JB017468 while I'm going to prepare the overall implementation for alternative magnitude calculation strategies. Particularly, the approach from Schaff & Richards, 2004 seems interesting while the SVD method seems to be only of use for very high CCC (i.e. >= 0.85).

[mmesim]

Let's try something similar to:

https://github.com/mmesim/qks_detection/blob/46ff5982c29b5fd952214823562b80f0707b79b4/src/get_magnitudes.m#L1-L41

The method is based on

Amplitude ratio between template and detected event
e.g https://www.nature.com/articles/ngeo697

Finally, we compute the magnitude of the detected event based on the median value of the maximum amplitude ratios for all channels between the template and detected events, assuming that a tenfold increase in amplitude corresponds to one unit increase in magnitude.

And also here:
https://www.science.org/doi/10.1126/sciadv.1601946

Magnitudes were calculated by using only phases with measured differential times. For these phases, peak amplitude ratios were calculated between the detected phase and the template phase. The final magnitude estimate was formed by taking the median over all amplitude ratios and assuming that a factor of 10 in amplitude corresponds to one magnitude unit difference (50).

[damb]

Perfect. Thanks for the references.

[damb]

Finally, we compute the magnitude of the detected event based on the median value of the maximum amplitude ratios for all channels between the template and detected events, assuming that a tenfold increase in amplitude corresponds to one unit increase in magnitude.

@mmesim, does this include also the channels which are not used for detection? E.g. given a detector is configured for NET.STA.LOC.HHZ, only. Is the magnitude computed including both NET.STA.LOC.HHN and NET.STA.LOC.HHE or based on NET.STA.LOC.HHZ, exclusively?

[damb]

And another question (for clarification): the approach seems to be used to compute an overall event magnitude for a new detection. Is this correct? Is this approach usable for computing station magnitudes, as well?

[damb]

Also: what kind of preprocessing is required? From https://github.com/mmesim/qks_detection/blob/46ff5982c29b5fd952214823562b80f0707b79b4/src/get_templates.m#L8-L10 and https://github.com/mmesim/qks_detection/blob/46ff5982c29b5fd952214823562b80f0707b79b4/src/get_templates.m#L36-L46 I can see that:

• the linear trend is removed
• the mean is removed
• a filter is applied (are there specific requirements which cannot be covered with the SeisComP filter grammar
• finally the waveform is trimmed.

Does the waveform around the detection require the same preprocessing?

[mmesim]

@mmesim, does this include also the channels which are not used for detection? E.g. given a detector is configured for NET.STA.LOC.HHZ, only. Is the magnitude computed including both NET.STA.LOC.HHN and NET.STA.LOC.HHE or based on NET.STA.LOC.HHZ, exclusively?

In my implementation I used all the available phases for the earthquake template. The magnitude is computed on the channels that are available and there is a P or S phase picked.

And another question (for clarification): the approach seems to be used to compute an overall event magnitude for a new detection. Is this correct? Is this approach usable for computing station magnitudes, as well?

Yes, this method calculates an event magnitude. The process is to find the maximum amplitude for each template and the corresponding detection for each station and channel. Then compute the ratio for all these pairs, and take the median. So, perhaps, you can do that per station.

Also: what kind of preprocessing is required? From https://github.com/mmesim/qks_detection/blob/46ff5982c29b5fd952214823562b80f0707b79b4/src/get_templates.m#L8-L10 and https://github.com/mmesim/qks_detection/blob/46ff5982c29b5fd952214823562b80f0707b79b4/src/get_templates.m#L36-L46 I can see that:

the linear trend is removed
the mean is removed
a filter is applied (are there specific requirements which cannot be covered with the SeisComP filter grammar
finally the waveform is trimmed.
Does the waveform around the detection require the same preprocessing?

Yes, this is the typical processing, and realtime data and templates should be similarly processed. No special filtering is required. In my case I applied a bandpass 2 pole butterworth filter (1.5-15 Hz).

[damb]

The magnitude is computed on the channels that are available and there is a P or S phase picked.

So should the time window include both the P and the S phase? Alternatively, we could rely on the time window(s) of the template waveform specified by the user within the templates.json template configuration. IMO, this would be more transparent to the user.

So, perhaps, you can do that per station.

This was my plan, anyway.

No special filtering is required. In my case I applied a bandpass 2 pole butterworth filter (1.5-15 Hz).

So I'd leave it up to the user which filter to apply. As a default, we could start with the filter you were using.

[mmesim]

Alternatively, we could rely on the time window(s) of the template waveform specified by the user within the templates.json template configuration. IMO, this would be more transparent to the user.

Exactly.

So I'd leave it up to the user which filter to apply. As a default, we could start with the filter you were using.

Also for the filters the default could be the one that it is defined in the .json for template configuration.

[damb]

Also for the filters the default could be the one that it is defined in the .json for template configuration.

Is it sensible to use same filter for both processing and amplitude calculation? If yes, this can be implemented, of course.

[damb]

Alternatively, we could rely on the time window(s) of the template waveform specified by the user within the templates.json template configuration. IMO, this would be more transparent to the user.

Exactly.

So, we're computing the amplitudes only from those streams which are used for detection, right?

[mmesim]

Is it sensible to use same filter for both processing and amplitude calculation?

I think so. In my code I do the prepossessing once and save waveforms and templates in a structure. Then instead of reprocess everything I used the already processed data. That's all.

[damb]

I think so. In my code I do the prepossessing once and save waveforms and templates in a structure. Then instead of reprocess everything I used the already processed data. That's all.

Unfortunately, in real-time this approach is not feasible.

[damb]

I think, I understand the problem now. Let's start with the implementation.

[damb]

See #70

[mmesim]

So, we're computing the amplitudes only from those streams which are used for detection, right?

Yes

[damb]

Thanks for all the details, @mmesim.

[mmesim]

To my understanding the current implementation returns only StationMagnitudes. Is it easy to calculate network magnitudes? For example mean or median of all the StationMagnitudes.

Network magnitudes

The network magnitude is a magnitude value summarizing several station magnitudes values of one origin. Different methods are available for summarizing the station magnitudes:

mean: the usual mean value

trimmed mean value: To stabilize the network magnitudes the smallest and the largest 12.5% of the station magnitude values are removed before computing the mean.

median: the usual median value

median trimmed mean: Removing all station magnitudes with a distance greater than 0.5 (default) from the median of all station magnitudes and computing the mean of all remaining station magnitudes.

[damb]

Is it easy to calculate network magnitudes?

Usually this is done with scmag.

[mmesim]

Ok. I tried that and seems that scmag prefers

Primary magnitudes Currently the following magnitude types are implemented:

ML
Local magnitude calculated on the vertical component using a correction term to fit with the standard ML

MLv
Local magnitude calculated on the vertical component using a correction term to fit with the standard ML

MLh
Local magnitude calculated on the horizontal components to SED specifications.

MLr
Local magnitude calculated from MLv amplitudes based on GNS/GEONET specifications for New Zealand.

MN
Canadian Nuttli magnitude.

mb
Narrow band body wave magnitude measured on a WWSSN-SP filtered trace

mBc
Cumulative body wave magnitude

mB
Broad band body wave magnitude

Mwp
The body wave magnitude of Tsuboi et al. (1995)

Mjma
Mjma is computed on displacement data using body waves of period < 30s

Ms_20
Surface-wave magnitude at 20 s period

Ms(BB)
Broad band surface-wave magnitude

Md
Duration magnitude as described in https://earthquake.usgs.gov/research/software/#HYPOINVERSE

It doesn't recognize MRelative or MLx.