Skip to content

[magnitudes] Add Mwpd duration-amplitude moment magnitude plugin#125

Open
comoglu wants to merge 8 commits into
SeisComP:mainfrom
comoglu:feature/mwpd-magnitude
Open

[magnitudes] Add Mwpd duration-amplitude moment magnitude plugin#125
comoglu wants to merge 8 commits into
SeisComP:mainfrom
comoglu:feature/mwpd-magnitude

Conversation

@comoglu

@comoglu comoglu commented Jun 24, 2026

Copy link
Copy Markdown
Contributor

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.

comoglu and others added 5 commits June 24, 2026 22:24
Native port of Early-est 1.2.9's Mwpd (Lomax & Michelini 2009): a
per-station moment magnitude from the integral of the vertical
displacement over the source duration T0.

Amplitude processor (type "Mwpd"): from the raw vertical window it
restitutes to velocity, high-passes at the BRB-HP corner (0.15 Hz),
integrates to displacement and accumulates the running double integral
into separate positive/negative displacement lobes. The source duration
T0 is estimated from the high-frequency envelope (1-5 Hz band-pass,
squared, boxcar-smoothed, 90/80/50/20 %-of-peak level crossings). It
emits the displacement integral at T0 in nm*s and carries T0 as the
amplitude period.

Magnitude processor (type "Mwpd"): forms M0 from the Tsuboi constant and
applies the INGV_EE distance correction, the PREM step depth correction
and the >90 s duration ramp.

Builds as plugin mwpd; no component combiner (vertical P only). Bindings
documented in descriptions/global_mwpd.xml.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…progressive updates

Reworks the Mwpd amplitude/magnitude processors after offline + real-time
playback testing against Mwp/mb/CMT:

- Source-duration T0 estimated on the SQUARED (power) HF envelope via the
  90/80/50/20%-of-peak last-crossing logic (Lomax & Michelini 2009, Fig. 2);
  envelope band ~1 Hz.
- Integration over min(maxDuration, S-P): the displacement double-integral is
  accumulated in positive/negative lobes and read at the window end, with the
  window capped at the theoretical S-P time (libtau/iasp91) so it excludes
  S/surface-wave energy. Matches Early-est's MAX_MWPD_DUR window behaviour;
  the eager short-T0 cut under-integrated and read ~0.4 mag low.
- Completeness gate: emit only once data covers the integration window;
  incomplete stations are skipped rather than reported on a partial window.
- Progressive updates (setUpdateEnabled) so the value is computed incrementally
  as data streams; finalizes when the window is complete.
- Moment scaling gated on magnitude (Mwpd += 0.45*(raw-7.2) for raw>7.2),
  eq (5a), instead of T0>90 s.
- S-P cap via a configurable travel-time table (libtau/iasp91).

On a teleseismic M~7.3 test this gives Mwpd 7.2, consistent with Mwp 7.3 /
mb 7.2. Constant remains Early-est's Tsuboi form; absolute calibration vs GCMT
to follow.
…gration

Calibration against GCMT (and Lomax & Michelini 2009, Table 1) showed Mwpd
reading increasingly low with magnitude (saturation). Root cause, confirmed
from the Early-est source: the displacement high-pass is
filter_hp_bu_co_0_005_n_2 = Butterworth 0.005 Hz, order 2 -- not 0.15 Hz
(MWPD_GAIN_FREQUENCY is only the gain-reference frequency, not a filter). A
0.15 Hz high-pass removes periods >7 s and strips the long-period moment that
dominates great-event integrals.

Changes:
- highpassCorner 0.005 Hz / order 2 (was 0.15 Hz / 4); config key renamed
  gainFrequency -> highpassCorner.
- HF envelope band restored to 1-5 Hz (filter_bp_bu_co_1_5_n_4).
- Integrate the displacement over min(T0, S-P, MAX_MWPD_DUR=1200 s) instead of
  a fixed window; finalize once the T0 search terminates or the window is full
  (progressive updates), else keep waiting / drop unresolved (MWPD_INVALID).
- Moment scaling switched to Early-est's active duration-gated form
  (T0>90 s ramp 90->110 s) * (raw-7.2) * 0.45.

Validated against GCMT: big events (M7.6-8.4) mean -0.15 vs CMT with no
saturation (e.g. M8.4 7.33->8.28); full run mean ~-0.05, sigma ~0.16,
slope ~1.0 (no high-M droop).
Method, amplitude/magnitude processing, the 0.005 Hz BRB-HP and T0-based
integration, configuration notes, the GCMT calibration result (slope ~1.02,
sigma ~0.15, no saturation) and references (Lomax & Michelini 2009; Tsuboi 1995).
…e docs

Embeds the validation against GCMT (972 events, Mw 6.0-8.8): Mwpd-vs-Mw scatter
figure plus per-magnitude and per-depth difference tables (slope 1.02, mean
+0.05, sigma 0.15, no saturation), in global_mwpd.rst.
@cla-bot cla-bot Bot added the cla-signed The CLA has been signed by all contributors label Jun 24, 2026
@gempa-jabe

Copy link
Copy Markdown
Contributor

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.

…l table

Per review (gempa-jabe): the plugin's own classes and helper functions were
exported (16 global symbols) although nothing outside the plugin uses them --
the loader needs only createSCPlugin and the processors register via static
factories. Give the cross-object-file symbols (readMwpdConfig, the moment/
distance/depth corrections, and the two processor classes) hidden visibility
via a small MWPD_LOCAL attribute; single-TU helpers were already in anonymous
namespaces. Result: zero exported Mwpd symbols, createSCPlugin unchanged,
plugin still loads and registers.
@comoglu

comoglu commented Jun 24, 2026

Copy link
Copy Markdown
Contributor Author

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

comoglu added 2 commits June 24, 2026 23:24
Drop the now-unused t0Resolved flag (resolution is tested via durRaw>0), which
also clears a -Wunused-but-set-variable warning, and correct the setup debug
label from "gainFreq" to "HPcorner".
- Clarify no-saturation applies to long-rupture events (tsunami / great quakes).
- Cite "Lomax & Michelini 2009" explicitly for eq. 3 and eq. 5a.
- Reword the progressive-computation note for real-time application.
- Replace the SeisComP "estimateMw identity" jargon with a plain statement
  that no further Mw conversion is applied and the network value is a robust
  average of station Mwpd.
- Note Mwpd(RT) (Lomax & Michelini 2013) as a lower-latency real-time variant
  and add the reference.
@jsaul

jsaul commented Jul 15, 2026

Copy link
Copy Markdown
Contributor

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 ( ... ) {
        ...
}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

cla-signed The CLA has been signed by all contributors

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants