[magnitudes] Add Mwpd duration-amplitude moment magnitude plugin

[comoglu]

Summary

Adds Mwpd — a per-station duration–amplitude moment magnitude for large earthquakes from teleseismic P waveforms, after Lomax & Michelini (2009) (a port of the Mwpd implementation in Early-est). Unlike Mwp, which uses the peak of the displacement integral, Mwpd integrates the P-displacement over the apparent source duration T0, so it does not saturate for great events.

The change is purely additive: a new plugin directory plugins/magnitudes/mwpd/ plus one SUBDIRS(mwpd) line in plugins/magnitudes/CMakeLists.txt. No existing behaviour changes. The plugin registers an amplitude processor and a magnitude processor, both of type Mwpd (vertical broad-band, amplitude unit nm*s).

Method

• Restitute counts → ground velocity, high-pass at 0.005 Hz (Butterworth, order 2 — the Early-est BRB-HP), single-integrate to displacement.
• Accumulate the running double integral into positive/negative displacement lobes; the larger lobe is the amplitude (eq. 3 of the paper).
• Estimate the apparent source duration T0 from the high-frequency (1–5 Hz) envelope via its 90/80/50/20 %-of-peak last-crossing times.
• Integrate over min(T0, S−P, MAX_MWPD_DUR) — the S−P cap (travel-time table) keeps the integral free of S/surface-wave energy.
• M0 = C_M · ∫u dt · Δ (Tsuboi constant), Mwpd = ⅔(log10 M0 − 9.1) minus the INGV/Early-est distance correction, plus the PREM step depth correction and the duration-gated moment-scaling ramp (T0 > 90 s) for great/slow events (eq. 5a).
• Computation is progressive (setUpdateEnabled): the value updates as data stream in and finalises once T0 resolves.

Validation

Calibrated against GCMT moment magnitudes for 972 events (Mw 6.0–8.8, 2015–2026): Mwpd tracks Mw(GCMT) along the 1:1 line with slope ≈ 1.02, mean difference +0.05, σ ≈ 0.15, no saturation and no significant depth dependence — consistent with operational Early-est and with the paper (±0.2). Per-magnitude and per-depth tables and the scatter figure are in descriptions/global_mwpd.rst.

Documentation

• descriptions/global_mwpd.rst — method, processing, calibration (figure + tables), configuration, references.
• descriptions/global_mwpd.xml — per-parameter binding descriptions for scconfig.

References

• Lomax, A. & Michelini, A. (2009), Geophys. J. Int. 176, 200–214, doi:10.1111/j.1365-246X.2008.03974.x
• Tsuboi, S. et al. (1995), Bull. Seismol. Soc. Am. 85, 606–613.

[gempa-jabe]

Interesting addition, thank you. From the technical side I prefer to implement as much as possible in private namespace to not expose symbol which won't be used outside the plugin and which are cluttering the symbol table. I still have no final recipe to do that with multiple object files, but visibility(hidden) and some compiler flags could help in that regard.

[comoglu]

Thanks! Addressed in 5fd3a1b: the plugin's classes and helpers used across object files now get hidden visibility via a small MWPD_LOCAL

[jsaul]

This looks pretty good. I haven't tried it but I trust you. 👍 Certainly a very valuable addition, thanks a lot!

In your comparison with Mw(GCMT) it would seem to me that Mw is slightly but systematically overestimated for Mw around 6 - 6.5 - is this just a visual impression? Generally speaking, even small trends can be corrected for by defining a mapping Mwpd -> Mw(Mwpd) just like we have it for other Mw proxies. Please evaluate if that would make sense. Can also be added later, of course.

Technically there are some minor code formatting suggestions especially avoid one-liners like

if ( ...) { ... }

but write them as

if ( ... ) {
        ...
}

[comoglu]

Thanks, @jsaul !

You are right on the low end , it's a real bias but small. It is about +0.07 on average for Mw 6-7, and it's inside the ~0.14 scatter. It's not really a clean slope either (7-7.5 actually sits at −0.05). Since Mwpd already comes out as a moment magnitude, I didn't want to bake in a network-specific fit, so I left it uncorrected by default. But anyone who wants to flatten it can do it from config with the usual scmag knob parameter magnitudes. Mwpd.multiplier/.offset for a linear station correction, or magnitudes.Mwpd.MwMapping for a proper Mw(Mwpd) mapping like the other proxies. I added a note about that to the docs.

Formatting is fixed too, expanding the one-liner ifs. Will push shortly.

[jsaul]

Awesome job, thanks again! @gempa-jabe do you have any objections to merging?

[gempa-jabe]

Actually this is not the decision of the AmplitudeProcessor. This is application controlled. And applications expect this value to be set to false initially. If we need update in scamp then it should be added there.

[jsaul]

Important point. As scamp only computes amplitudes, an amplitude update would be required for each incoming data record until the AmplitudeProcessor_Mwpd decides that it has finished its job, which may take several minutes in case of very large EQs. These updates would be needed in order to allow the downstream (especially scmag) to capture all intermediate amplitudes to actually make use of them. It will then also have to generate intermediate station magnitudes to continuously update the network magnitude.

My concern would be that this might stress the processing/messaging at a time (soon after a potentially catastrophic EQ) where smooth processing is absolutely crucial.

[gempa-jabe]

As usual, I prefer private namespaces for plugin classes to not expose symbols into the global namespaces which are actually not required. Using private linkage for the plugin is my preference.

[comoglu]

A correction to my earlier change. I removed setUpdateEnabled(true) after your review, but testing it in a real playback showed it stops Mwpd from working at all.

With updates disabled, the base only calls computeAmplitude once the full signalEnd window (1200 s) has streamed. Mwpd is a duration amplitude, it has to watch the data to find where the source stops (T0), so it needs progressive evaluation. Without it, T0 never resolves in time, and no Mwpd is produced.

@jsaul, on your concern about messaging load: I tested this on a real playback (Japan M7.7) with updates enabled. computeAmplitude returns InProgress on every record until T0 resolves and then emits the amplitude exactly once, so there is no per-record spam. scamp stayed at 3.9-6.6% CPU the whole run and produced a clean network Mwpd (~7.4, 10 stations, 6.97-7.65 range) with no messaging issues. So I think the concern does not apply here in practice, but let me know if you would rather see it differently.

I have put updates back on (8451f49) with a comment in the code explaining why.

[gempa-jabe]

With updates disabled, the base only calls computeAmplitude once the full signalEnd window (1200 s) has streamed. Mwpd is a duration amplitude, it has to watch the data to find where the source stops (T0), so it needs progressive evaluation. Without it, T0 never resolves in time, and no Mwpd is produced.

That explanation does not make much sense to me. Updates just allow computation before the full time window has been streamed. But that is not guaranteed, also the time steps are not guaranteed and rely on the size of the records. It could be possible that you will receive one huge record covering the entire 1200s window. Not common but possible. So the internal logic should be able to handle that. Would record size of 4096 by also cause issues by e.g. computing a different T0? What would be the difference in feeding data sample by sample?

[jsaul]

Amplitude updates for large events with long rupture duration are extremely useful of course. But with the extreme example by @gempa-jabe in mind, what if records are really short and updates are triggered by every short record? That may indeed be more updates than ever needed by the downstream.

IMHO there needs to be a parameter that controls the minimum interval between amplitude updates sent to the messaging. Probably once every 30 seconds would be a good default value. Also the downstream must be able to know whether an amplitude computation is finished or not yet.

[comoglu]

Good point. didn't want to just move on, so I had Claude walk the actual code path instead of assuming.
Here they are.
Two things settle it:

TimeWindowProcessor::fill() (common/libs/seiscomp/processing/timewindowprocessor.cpp:84-86) just does _data.append(n, samples) — strictly append-only, no truncation, no ring buffer, no eviction. Whatever chunking arrives on the wire, the data/i2 handed to computeAmplitude is always the complete buffer from the start of the signal window, on every call.

AmplitudeProcessor_Mwpd carries no cross-call state — its only members are the config, the output array, the tt-table pointer, and the cached _maxDur. Nothing stores T0 or level-crossing progress between calls. Every computeAmplitude call rebuilds velocity, the HP/BP filters, and rescans the 90/80/50/20% crossing loop from iPick to n from scratch, with fresh filter instances and purely local state.

So for a given cumulative sample count n, the array being scanned is bit-identical no matter how it arrived — one 1200s record or a thousand small ones. Record size changes when we happen to check (updates=true just lets us check before signalEnd), never what the answer is when we do. A 4096-sample record wouldn't give a different T0 than the same samples split up — same array, same scan, same result.

What do you think ? @gempa-jabe

[comoglu]

Amplitude updates for large events with long rupture duration are extremely useful of course. But with the extreme example by @gempa-jabe in mind, what if records are really short and updates are triggered by every short record? That may indeed be more updates than ever needed by the downstream.

IMHO there needs to be a parameter that controls the minimum interval between amplitude updates sent to the messaging. Probably once every 30 seconds would be a good default value. Also the downstream must be able to know whether an amplitude computation is finished or not yet.

Checked the code for this too. emitAmplitude() (the call that sends the message) only fires when computeAmplitude() returns true and Mwpd returns false on every call except the one where T0 actually resolves. So regardless of record size, Mwpd sends at most one Amplitude message per station, ever. No growing sequence of values to throttle, so I don't think the min-interval concern applies here that seems more relevant for amplitude types that do stream intermediate values.

(Side note: the base's Finished status flag lags it's only set once the full 1200s window streams, even though our value went out earlier. Internal only, doesn't affect the message scmag sees.)

I will be checking out in the morning.

[comoglu]

Small correction to my last comment. computeAmplitude does not return false after T0 resolves; it keeps getting called and keeps returning true with the same value. What actually blocks the repeat messages is the base class's own dedup in process(), it only calls emitAmplitude when the new amplitude is strictly greater than the last one, and ours is never larger on the repeat calls since it is a deterministic recompute over unchanged history. So still one message per station, just via the base's dedup, not our own gating.

[gempa-jabe]

You can try to call setStatus(Finished, 100) once you have computed an amplitude. Then your processing should never be called again.

To recap: setting the update flag does not change the way amplitudes are computed or if they are computed at all, it is just a matter of when they are computed, right?

Given the algorithm and nature of the amplitude, I am finde with enabling updates in that particular case.

[comoglu]

One more data point on this. I also went back to the actual Early-est source we ported from (timedomain_processing.c, https://github.com/comoglu/early-est/blob/e69c669d1e7a98a3f5a7aa4b77ef40327bec3f39/timedomain_processing/timedomain_processing.c#L1541-L1691. Turns out Early-est does the opposite of what we do: it is incremental, it keeps a persistent per-station T0 state (amp_peak, amp_90/80/50/20 and their indices) in a struct, and each call only processes the newly arrived samples, updating that running state. It never rescans history.

My port went the other way. No state between calls, full rescan of the whole buffer every time. That works because SeisComP hands us an append-only buffer that never gets truncated, so for a given amount of data we always see the same array no matter how it arrived.

So both designs land on the same guarantee; chunking cannot change the result, just by different means. Early-est does it with careful incremental bookkeeping, and I do it by not keeping state at all and always looking at everything.

Let me know if that matches your understanding of the framework or if we misread something in the SeisComP base class, you know that code better than I do. Your input is highly appreciated.

[comoglu]

You can try to call setStatus(Finished, 100) once you have computed an amplitude. Then your processing should never be called again.

To recap: setting the update flag does not change the way amplitudes are computed or if they are computed at all, it is just a matter of when they are computed, right?

Given the algorithm and nature of the amplitude, I am finde with enabling updates in that particular case.

Right, thanks heaps, updates only change when computeAmplitude gets called, not whether or what it computes.

on setStatus(Finished, 100). Checked the codebase, MLc and MLh already do exactly this once their result is final. I will add it to Mwpd's computeAmplitude once the amplitude is resolved, that stops store() from calling process() again at all instead of relying on the base's dedup to absorb the repeats. Will push shortly.

[comoglu]

Correction, sorry for the back and forth. It turns out computeAmplitude already calls setStatus(Finished, 100.0) the moment it resolves; that has been there since an early commit (de4a13e). So it does not get called again after that, store() rejects every later record before process() runs at all. My last two comments describing the base's dedup as the mechanism were wrong, I mis-checked for this. Your suggestion describes exactly what the code already does, nothing to add.