[magnitudes] Add Mwpd duration-amplitude moment magnitude plugin#125
Open
comoglu wants to merge 8 commits into
Open
[magnitudes] Add Mwpd duration-amplitude moment magnitude plugin#125comoglu wants to merge 8 commits into
comoglu wants to merge 8 commits into
Conversation
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.
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 |
…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.
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 |
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.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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). UnlikeMwp, which uses the peak of the displacement integral,Mwpdintegrates the P-displacement over the apparent source durationT0, so it does not saturate for great events.The change is purely additive: a new plugin directory
plugins/magnitudes/mwpd/plus oneSUBDIRS(mwpd)line inplugins/magnitudes/CMakeLists.txt. No existing behaviour changes. The plugin registers an amplitude processor and a magnitude processor, both of typeMwpd(vertical broad-band, amplitude unitnm*s).Method
T0from the high-frequency (1–5 Hz) envelope via its 90/80/50/20 %-of-peak last-crossing times.min(T0, S−P, MAX_MWPD_DUR)— theS−Pcap (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).setUpdateEnabled): the value updates as data stream in and finalises onceT0resolves.Validation
Calibrated against GCMT moment magnitudes for 972 events (Mw 6.0–8.8, 2015–2026):
MwpdtracksMw(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 indescriptions/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