diff --git a/plugins/magnitudes/CMakeLists.txt b/plugins/magnitudes/CMakeLists.txt index e161ff027..86ac0447f 100644 --- a/plugins/magnitudes/CMakeLists.txt +++ b/plugins/magnitudes/CMakeLists.txt @@ -1 +1,2 @@ -SUBDIRS(nuttli) \ No newline at end of file +SUBDIRS(nuttli) +SUBDIRS(mwpd) diff --git a/plugins/magnitudes/mwpd/CMakeLists.txt b/plugins/magnitudes/mwpd/CMakeLists.txt new file mode 100644 index 000000000..d03f372b1 --- /dev/null +++ b/plugins/magnitudes/mwpd/CMakeLists.txt @@ -0,0 +1,21 @@ +SET(PLUGIN_TARGET mwpd) + +SET( + PLUGIN_SOURCES + config.cpp + amplitude.cpp + magnitude.cpp + plugin.cpp +) + +SET( + PLUGIN_HEADERS + mwpd.h + version.h +) + +SC_ADD_PLUGIN_LIBRARY(PLUGIN ${PLUGIN_TARGET} "") +SC_LINK_LIBRARIES_INTERNAL(${PLUGIN_TARGET} client) + +FILE(GLOB descs "${CMAKE_CURRENT_SOURCE_DIR}/descriptions/*.xml") +INSTALL(FILES ${descs} DESTINATION ${SC3_PACKAGE_APP_DESC_DIR}) diff --git a/plugins/magnitudes/mwpd/amplitude.cpp b/plugins/magnitudes/mwpd/amplitude.cpp new file mode 100644 index 000000000..1ebb3d71e --- /dev/null +++ b/plugins/magnitudes/mwpd/amplitude.cpp @@ -0,0 +1,432 @@ +/*************************************************************************** + * SeisComP duration-amplitude moment-magnitude plugin (mwpd) * + * * + * Amplitude processor. From the raw vertical window it reproduces two * + * Early-est streams: * + * 1. BRB-HP displacement: counts/gain -> velocity, high-passed at * + * MWPD_GAIN_FREQUENCY, single-integrated to displacement; the running * + * double integral is accumulated into separate positive- and * + * negative-displacement lobes (int_int_sum_pos/neg). * + * 2. HF envelope: raw band-passed 1-5 Hz, squared, boxcar-smoothed; its * + * peak/level (90/80/50/20 %) crossings give the source duration T0. * + * The emitted amplitude is the displacement integral at T0 (nm*s); T0 is * + * carried as the amplitude "period". * + * * + * Ports timedomain_processing.c (mwpd + T0 loops) and * + * timedomain_processing_data.c (calculate_Raw_Mwpd_Mag, calculate_duration).* + * * + * GNU Affero General Public License Usage - see LICENSE. * + ***************************************************************************/ + + +#define SEISCOMP_COMPONENT Mwpd + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "mwpd.h" + + +using namespace Seiscomp::Math::Filtering::IIR; + + +namespace Seiscomp { +namespace Magnitudes { +namespace Mwpd { + + +namespace { + +// Centred boxcar smoothing (half-width in samples), matching the T50 envelope +// smoothing. Edges shrink the window to the available samples. +void boxcarSmooth(std::vector &x, int halfWidth) { + if ( halfWidth < 1 || x.size() < 2 ) { + return; + } + const int n = static_cast(x.size()); + + // Prefix sums for an O(n) moving average. + std::vector pre(n + 1, 0.0); + for ( int i = 0; i < n; ++i ) { + pre[i + 1] = pre[i] + x[i]; + } + + std::vector out(n); + for ( int i = 0; i < n; ++i ) { + int a = std::max(0, i - halfWidth); + int b = std::min(n - 1, i + halfWidth); + out[i] = (pre[b + 1] - pre[a]) / (b - a + 1); + } + x.swap(out); +} + +// Source-duration estimate from the level-crossing indices (relative to the +// pick, in samples). Ports calculate_duration(). +double durationFromLevels(int idx90, int idx80, int idx50, int idx20, + double dt, double floorSec) { + const double t9 = idx90 * dt; + const double t8 = idx80 * dt; + const double t5 = idx50 * dt; + const double t2 = idx20 * dt; + + const double durEstimate = (t5 + t8) / 2.0; + double factor = (durEstimate - 20.0) / 40.0; // 0->1 for 20->60 s + if ( factor < 0.0 ) { + factor = 0.0; + } + else if ( factor > 1.0 ) { + factor = 1.0; + } + + double duration = t9 * (1.0 - factor) + t2 * factor; + if ( duration < floorSec ) { + duration = floorSec; + } + return duration; +} + +} // namespace + + +REGISTER_AMPLITUDEPROCESSOR(AmplitudeProcessor_Mwpd, MWPD_TYPE); + + +AmplitudeProcessor_Mwpd::AmplitudeProcessor_Mwpd() +: Processing::AmplitudeProcessor(MWPD_TYPE) { + setUsedComponent(Vertical); + setUnit(MWPD_AMP_UNIT); + // Compute progressively as data streams (Early-est style): the processor + // emits a growing Mwpd and finalizes as soon as T0 resolves/terminates, + // rather than waiting the full (long) signal window. maxDuration is just + // the upper bound on how long we keep looking. + setUpdateEnabled(true); + applyConfig(); +} + + +void AmplitudeProcessor_Mwpd::applyConfig() { + // Pick at signalStart=0; integration starts analysisPreP before it (taken + // from the noise side of the window). Signal end covers the longest T0. + setSignalStart(0.0); + setSignalEnd(_cfg.maxDuration); + setNoiseStart(-(_cfg.analysisPreP + 30.0)); + setNoiseEnd(-_cfg.analysisPreP); + + setMinDist(_cfg.minDistanceDeg); + setMaxDist(_cfg.maxDistanceDeg); + setMinSNR(0); // T0 has its own S/N termination; no scalar SNR gate + computeTimeWindow(); +} + + +bool AmplitudeProcessor_Mwpd::setup(const Processing::Settings &settings) { + if ( !Processing::AmplitudeProcessor::setup(settings) ) { + return false; + } + + _cfg = MwpdConfig(); + if ( !readMwpdConfig(settings, std::string("amplitudes.") + type(), _cfg) ) { + return false; + } + + applyConfig(); + + // Travel-time table for the S-P window cap (graceful if unavailable). + _ttt = nullptr; + if ( _cfg.useSpCap ) { + _ttt = Seiscomp::TravelTimeTableInterface::Create(_cfg.tttBackend.c_str()); + if ( !_ttt || !_ttt->setModel(_cfg.tttModel) ) { + SEISCOMP_WARNING("%s: travel-time table '%s'/'%s' unavailable; " + "S-P window cap disabled", + type().c_str(), _cfg.tttBackend.c_str(), + _cfg.tttModel.c_str()); + _ttt = nullptr; + } + } + + SEISCOMP_DEBUG("%s: HPcorner=%.3fHz HF=%.1f-%.1fHz preP=%.1fs maxDur=%.0fs spCap=%d", + type().c_str(), _cfg.highpassCorner, _cfg.hfFmin, _cfg.hfFmax, + _cfg.analysisPreP, _cfg.maxDuration, _ttt != nullptr); + return true; +} + + +double AmplitudeProcessor_Mwpd::computeSPSeconds() const { + if ( !_ttt || !_environment.hypocenter || !_environment.receiver ) { + return -1.0; + } + try { + const double slat = _environment.hypocenter->latitude().value(); + const double slon = _environment.hypocenter->longitude().value(); + double sdep = 0.0; + try { sdep = _environment.hypocenter->depth().value(); } + catch ( ... ) {} + const double rlat = _environment.receiver->latitude(); + const double rlon = _environment.receiver->longitude(); + double relev = 0.0; + try { relev = _environment.receiver->elevation(); } + catch ( ... ) {} + + const double tP = _ttt->computeTime("P", slat, slon, sdep, rlat, rlon, relev); + const double tS = _ttt->computeTime("S", slat, slon, sdep, rlat, rlon, relev); + if ( tP > 0.0 && tS > 0.0 && tS > tP ) { + return tS - tP; + } + } + catch ( ... ) {} + return -1.0; +} + + +const DoubleArray *AmplitudeProcessor_Mwpd::processedData(Component comp) const { + return comp == targetComponent() ? &_processedData : nullptr; +} + + +bool AmplitudeProcessor_Mwpd::computeAmplitude( + const DoubleArray &data, + size_t i1, size_t i2, + size_t /*si1*/, size_t /*si2*/, + double offset, + AmplitudeIndex *dt, + AmplitudeValue *amplitude, + double *period, double *snr) { + + const double fsamp = _stream.fsamp; + if ( fsamp <= 0.0 ) { + setStatus(Error, 1); + return false; + } + const double deltaTime = 1.0 / fsamp; + + const double gain = std::fabs(_streamConfig[targetComponent()].gain); + if ( gain == 0.0 ) { + setStatus(MissingGain, 1); + return false; + } + + const int n = static_cast(i2); // end of usable window + const int iPick = static_cast(i1); // signalStart = pick + const int preP = static_cast(_cfg.analysisPreP * fsamp + 0.5); + const int iStart = iPick - preP; // integration start + if ( iStart < 0 || n > static_cast(data.size()) || n - iPick < 16 ) { + // Not enough data past the pick yet. With incremental updates this is a + // "keep waiting" state, NOT an error: leave the status InProgress so the + // base keeps feeding. (A terminal status here would kill the processor.) + return false; + } + + const double *raw = static_cast(data.data()); + + // ===================================================================== + // 1) BRB-HP displacement and its positive/negative running integrals. + // ===================================================================== + // Velocity (counts -> m/s), demeaned by the noise offset, high-passed at + // the BRB-HP corner. We keep it as _processedData for inspection. + _processedData.resize(n); + double *vel = _processedData.typedData(); + for ( int i = 0; i < n; ++i ) { + vel[i] = (raw[i] - offset) / gain; + } + + ButterworthHighpass hp(_cfg.hpOrder, _cfg.highpassCorner, fsamp); + hp.apply(n, vel); + + // Running single integral (displacement) and double integral split into + // positive- and negative-displacement lobes (Early-est mwpd loop). + const int nAcc = (n - iStart) + 1; // index 0 == analysis start + std::vector intPos(nAcc, 0.0); + std::vector intNeg(nAcc, 0.0); + double intSum = 0.0, accPos = 0.0, accNeg = 0.0; + int k = 0; + for ( int i = iStart; i < n; ++i ) { + ++k; + intSum += vel[i] * deltaTime; // displacement [m] + const double area = intSum * deltaTime; // [m*s] + if ( intSum >= 0.0 ) { + accPos += area; + } + else { + accNeg += -area; + } + intPos[k] = accPos; + intNeg[k] = accNeg; + } + + // ===================================================================== + // 2) High-frequency envelope and the source duration T0 (durRaw). + // ===================================================================== + double durRaw; // source-duration estimate [s] (used for the ramp) + bool t0Terminated = false; // T0 search hit a termination criterion (final) + double peakRelSec = -1.0; // time of the envelope peak after P (debug) + + if ( _cfg.fixedDuration ) { + durRaw = _cfg.fixedDurationVal; + t0Terminated = true; + } + else { + std::vector env(n); + for ( int i = 0; i < n; ++i ) { + env[i] = raw[i]; + } + ButterworthBandpass bp(_cfg.hfOrder, _cfg.hfFmin, _cfg.hfFmax, fsamp); + bp.apply(n, env.data()); + for ( int i = 0; i < n; ++i ) { + env[i] *= env[i]; // square + } + const int smoothHW = static_cast(_cfg.smoothHalfWidth * fsamp + 0.5); + boxcarSmooth(env, smoothHW); + // NOTE: env stays as the squared (power) envelope -- Early-est's + // data_float_t50 is squared-then-smoothed and is NOT square-rooted; the + // 90/80/50/20% level tracking below runs on power, so a "20% of peak" + // crossing is at 0.2 in power (= ~0.45 in amplitude). + + // Level tracking from the pick onward (Early-est T0 loop). + double ampNoise = env[iPick]; + if ( ampNoise <= 0.0 ) { + ampNoise = DBL_MIN; + } + double ampPeak = -1.0; + int idxPeak = 0; + double amp90 = -1.0, amp80 = -1.0, amp50 = -1.0, amp20 = -1.0; + int idx90 = 0, idx80 = 0, idx50 = 0, idx20 = 0; + const int minDurSamp = static_cast(_cfg.minDuration * fsamp + 0.5); + const int smooth2 = static_cast(2.0 * _cfg.smoothHalfWidth * fsamp + 0.5); + durRaw = -1.0; + + for ( int i = iPick; i < n; ++i ) { + const double amp = env[i]; + const int rel = i - iPick; + + if ( amp > ampPeak ) { // new peak: unset all levels + ampPeak = amp; + idxPeak = rel; + amp90 = amp80 = amp50 = amp20 = -1.0; + } + else { + if ( amp90 < 0.0 ) { + if ( amp <= 0.9 * ampPeak ) { amp90 = amp; idx90 = rel; } + } + else if ( amp > 0.9 * ampPeak ) { + amp90 = amp80 = amp50 = amp20 = -1.0; + } + if ( amp80 < 0.0 ) { + if ( amp <= 0.8 * ampPeak ) { amp80 = amp; idx80 = rel; } + } + else if ( amp > 0.8 * ampPeak ) { + amp80 = amp50 = amp20 = -1.0; + } + if ( amp50 < 0.0 ) { + if ( amp <= 0.5 * ampPeak ) { amp50 = amp; idx50 = rel; } + } + else if ( amp > 0.5 * ampPeak ) { + amp50 = amp20 = -1.0; + } + if ( amp20 < 0.0 ) { + if ( amp <= 0.2 * ampPeak ) { + amp20 = amp; idx20 = rel; + durRaw = durationFromLevels(idx90, idx80, idx50, idx20, + deltaTime, _cfg.durationFloor); + } + } + else if ( amp > 0.2 * ampPeak ) { + amp20 = -1.0; + } + + // Terminate once below the S/N or peak ratio (after min dur). + if ( rel > minDurSamp && + ( amp / ampNoise < _cfg.snT0End || + amp / ampPeak < _cfg.peakRatioT0End ) ) { + t0Terminated = true; + break; + } + } + + // Terminate if well past twice the established duration. + if ( durRaw > 0.0 && amp20 > 0.0 && + rel > static_cast(2.0 * durRaw * fsamp + 0.5) && + rel > smooth2 ) { + t0Terminated = true; + break; + } + } + + peakRelSec = idxPeak * deltaTime; + } + + // ===================================================================== + // 3) S-P-capped integration over the source duration T0. + // ===================================================================== + // Early-est integrates the displacement over min(T0, S-P, MAX_MWPD_DUR) + // (calculate_Raw_Mwpd_Mag / timedomain_processing.c ~1762). Under progressive + // updates we only finalize once T0 is stable -- its search terminated + // (t0Terminated) or the full window has streamed -- otherwise the value would + // be locked on a transient early T0. Until then, keep waiting (InProgress); + // if T0 never resolves by window end the base drops the station (MWPD_INVALID). + if ( durRaw <= 0.0 ) { + return false; // T0 not resolved yet -> wait + } + const double availDur = (n - iPick) * deltaTime; + const bool windowFull = availDur >= _cfg.maxDuration - 1.0; + if ( !t0Terminated && !windowFull ) { + return false; // T0 still provisional -> wait + } + + if ( durRaw > _cfg.maxDuration ) { + durRaw = _cfg.maxDuration; + } + const double sp = computeSPSeconds(); + double integDur = durRaw; + if ( sp > 0.0 && integDur > sp ) { + integDur = sp; // S-P cap on the integration window + } + if ( availDur < integDur - 1.0 ) { + return false; // data must cover the window -> wait + } + + int idxDur = static_cast((integDur + _cfg.analysisPreP) / deltaTime + 0.5); + if ( idxDur >= nAcc ) { + idxDur = nAcc - 1; + } + if ( idxDur < 1 ) { + return false; + } + + double amplitudeIntegral = std::max(intPos[idxDur], intNeg[idxDur]); // [m*s] + if ( amplitudeIntegral <= 0.0 ) { + return false; + } + + // Emit the integral in nm*s (as Mwp emits nm); carry the source-duration T0 + // as the period (the base divides period by fsamp, so multiply here). + amplitude->value = amplitudeIntegral * 1.0e9; + *period = durRaw * fsamp; + *snr = 1000000.0; + + dt->index = iPick; + dt->begin = iStart - iPick; + dt->end = idxDur + iStart - iPick; + + setStatus(Finished, 100.0); // T0 stable -> finalize + SEISCOMP_DEBUG("%s.%s.%s: Mwpd FINAL ampInt=%g nm*s T0raw=%.1fs peak@%.1fs " + "S-P=%.1fs integDur=%.1fs win=%.0fs (pos=%g neg=%g) gain=%g", + _environment.networkCode.c_str(), + _environment.stationCode.c_str(), + _environment.locationCode.c_str(), + amplitude->value, durRaw, peakRelSec, sp, integDur, availDur, + intPos[idxDur] * 1.0e9, intNeg[idxDur] * 1.0e9, gain); + + return true; +} + + +} +} +} diff --git a/plugins/magnitudes/mwpd/config.cpp b/plugins/magnitudes/mwpd/config.cpp new file mode 100644 index 000000000..12c1500a1 --- /dev/null +++ b/plugins/magnitudes/mwpd/config.cpp @@ -0,0 +1,123 @@ +/*************************************************************************** + * SeisComP duration-amplitude moment-magnitude plugin (mwpd) * + * * + * Shared configuration reader and the Early-est constants/corrections * + * (moment constant, INGV_EE distance correction, PREM step depth * + * correction). * + * * + * GNU Affero General Public License Usage - see LICENSE. * + ***************************************************************************/ + + +#define SEISCOMP_COMPONENT Mwpd + +#include + +#include + +#include "mwpd.h" + + +namespace Seiscomp { +namespace Magnitudes { +namespace Mwpd { + + +// PREM depth correction table (depth_corr_mwpd.h), depth[km] -> delta-magnitude, +// 0..600 km. Step lookup: the correction of the largest table depth <= h. +namespace { + +const int NUM_DEPTH_CORR = 11; +const double DEPTH_CORR[NUM_DEPTH_CORR][2] = { + { 0.0, 0.00 }, { 71.0, 0.00 }, { 80.0, 0.00 }, {171.0, -0.01 }, + {220.0, 0.05 }, {271.0, 0.06 }, {371.0, 0.07 }, {400.0, 0.12 }, + {471.0, 0.15 }, {571.0, 0.18 }, {600.0, 0.22 } +}; + +} + + +double mwpdMomentConstant(const MwpdConfig &cfg) { + // Early-est MWP_CONST (Tsuboi 1995): + // 4 pi * rho * Vp^3 * FP * (10000/90 deg->km) * 1000 (km->m) + return 4.0 * M_PI * cfg.rho * cfg.vp * cfg.vp * cfg.vp * cfg.fp + * (10000.0 / 90.0) * 1000.0; +} + + +double mwpdDistanceCorrection(double deltaDeg, double depthKm) { + // calculate_Mwp_correction_INGV_EE (20150312 coefficients). + if ( deltaDeg < 0.0 || depthKm > 100.0 ) { + return 0.0; + } + return -1.32e-6 * deltaDeg * deltaDeg * deltaDeg + + 2.40e-4 * deltaDeg * deltaDeg + - 0.0146 * deltaDeg + + 0.314; +} + + +double mwpdDepthCorrection(double depthKm) { + // get_depth_corr_mwpd_prem: largest table depth <= h; -999 if below row 0. + int index = NUM_DEPTH_CORR - 1; + while ( index >= 0 && DEPTH_CORR[index][0] > depthKm ) { + --index; + } + return ( index < 0 ) ? -999.0 : DEPTH_CORR[index][1]; +} + + +bool readMwpdConfig(const Processing::Settings &settings, + const std::string &prefix, MwpdConfig &out) { + // BRB-HP restitution + displacement integral. + settings.getValue(out.highpassCorner, prefix + ".highpassCorner"); + settings.getValue(out.hpOrder, prefix + ".highpassOrder"); + settings.getValue(out.analysisPreP, prefix + ".analysisPreP"); + + // HF envelope. + settings.getValue(out.hfFmin, prefix + ".hfFmin"); + settings.getValue(out.hfFmax, prefix + ".hfFmax"); + settings.getValue(out.hfOrder, prefix + ".hfOrder"); + settings.getValue(out.smoothHalfWidth, prefix + ".smoothHalfWidth"); + + // T0 estimation. + settings.getValue(out.minDuration, prefix + ".minDuration"); + settings.getValue(out.maxDuration, prefix + ".maxDuration"); + settings.getValue(out.snT0End, prefix + ".snT0End"); + settings.getValue(out.peakRatioT0End, prefix + ".peakRatioT0End"); + settings.getValue(out.fixedDuration, prefix + ".fixedDuration"); + settings.getValue(out.fixedDurationVal,prefix + ".fixedDurationValue"); + + // Magnitude formula + corrections. + settings.getValue(out.rho, prefix + ".density"); + settings.getValue(out.vp, prefix + ".vp"); + settings.getValue(out.fp, prefix + ".radiation"); + settings.getValue(out.useDistanceCorr, prefix + ".distanceCorrection"); + settings.getValue(out.useDepthCorr, prefix + ".depthCorrection"); + settings.getValue(out.useDurationRamp, prefix + ".durationRamp"); + settings.getValue(out.minDistanceDeg, prefix + ".minimumDistance"); + settings.getValue(out.maxDistanceDeg, prefix + ".maximumDistance"); + + // S-P window cap. + settings.getValue(out.useSpCap, prefix + ".spWindowCap"); + settings.getValue(out.tttBackend, prefix + ".tttBackend"); + settings.getValue(out.tttModel, prefix + ".tttModel"); + + if ( out.hfFmax <= out.hfFmin ) { + SEISCOMP_ERROR("%s: hfFmax (%f) must be > hfFmin (%f)", + prefix.c_str(), out.hfFmax, out.hfFmin); + return false; + } + if ( out.maxDuration <= out.minDuration ) { + SEISCOMP_ERROR("%s: maxDuration (%f) must be > minDuration (%f)", + prefix.c_str(), out.maxDuration, out.minDuration); + return false; + } + + return true; +} + + +} +} +} diff --git a/plugins/magnitudes/mwpd/descriptions/global_mwpd.rst b/plugins/magnitudes/mwpd/descriptions/global_mwpd.rst new file mode 100644 index 000000000..989979860 --- /dev/null +++ b/plugins/magnitudes/mwpd/descriptions/global_mwpd.rst @@ -0,0 +1,157 @@ +:math:`M_{wpd}` is a per-station **duration--amplitude moment magnitude** for +large earthquakes, determined from teleseismic *P* waveforms. It integrates the +broad-band *P*-wave ground displacement over the apparent source duration +:math:`T_0`, after Lomax & Michelini (2009), and is a port of the Mwpd +implementation in Early-est (A. Lomax). Because it integrates over the whole +source duration rather than taking the peak (as :math:`M_{wp}` does), it does +**not saturate** for events with long rupture duration, such as tsunami +earthquakes and great earthquakes. + +The plugin contributes both an amplitude processor and a magnitude processor of +type ``Mwpd``. Both use the vertical broad-band component only. + +Method +====== + +Given the far-field *P* displacement :math:`u(t)` for a source of duration +:math:`T_0`, the scalar seismic moment is (Tsuboi et al. 1995; Lomax & +Michelini 2009) + +.. math:: + + M_0 = C_M \int_{t_P}^{t_P + T_0} u(t)\, \mathrm{d}t , + +where :math:`t_P` is the *P*-arrival time. The implementation accumulates the +running displacement integral separately over its **positive and negative +lobes** (to separate the direct *P* wave from later reflected/secondary phases +of opposite polarity, eq. 3 of Lomax & Michelini 2009) and uses the larger of the two at +:math:`T_0`. The moment magnitude follows the standard relation + +.. math:: + + M_{wpd}^{\mathrm{raw}} = \tfrac{2}{3}\left(\log_{10} M_0 - 9.1\right) + - C_\Delta(\Delta, h) , + +with a distance correction :math:`C_\Delta` (INGV/Early-est, subtracted; zero +for :math:`h>100` km), then a PREM step depth correction :math:`C_h(h)` and a +moment-scaling term for great/slow events: + +.. math:: + + M_{wpd} = M_{wpd}^{\mathrm{raw}} + C_h(h) + + r\,(M_{wpd}^{\mathrm{raw}} - 7.2)\cdot 0.45 , + \qquad r = \mathrm{clip}\!\left(\tfrac{T_0-90}{110-90},\,0,\,1\right). + +The duration ramp :math:`r` engages only for :math:`T_0 > 90` s (the moment +scaling for large interplate-thrust / tsunamigenic events, eq. 5a in Lomax & +Michelini 2009). + +Amplitude +========= + +For each *P* pick on the vertical broad-band, the amplitude processor: + +#. restitutes counts to ground velocity (sensitivity gain) and high-passes with + a Butterworth filter at a **0.005 Hz** corner (the Early-est BRB-HP filter); + a low corner is essential to retain the long-period displacement that + carries the moment of great earthquakes; +#. single-integrates to displacement and accumulates the running double + integral into positive/negative lobes; +#. estimates the **apparent source duration** :math:`T_0` from the + high-frequency envelope: the velocity is band-passed 1--5 Hz, squared and + boxcar-smoothed, and :math:`T_0` is derived from the times at which this + envelope last drops below 90/80/50/20 % of its peak; +#. integrates over :math:`\min(T_0,\ t_S-t_P,\ \mathtt{maxDuration})` --- the + :math:`S\!-\!P` cap (from a travel-time table) keeps the integral free of + *S* and surface-wave energy. + +The amplitude carries the displacement integral (unit ``nm*s``); :math:`T_0` is +carried as the amplitude period. In real-time application, the computation is +progressive: the processing is updated as new data arrives and finalised once +:math:`T_0` is resolved. + +Magnitude +========= + +The magnitude processor reads the displacement integral (``nm*s``) and +:math:`T_0` (period) and applies the relations above, with the moment constant + +.. math:: + + C_M = 4\pi\,\rho\,V_p^{3}\,F_p \cdot \tfrac{10000}{90}\cdot 1000 + \approx 4.68\times10^{21} + +(Tsuboi constant; :math:`\rho=3400`, :math:`V_p=7900`, :math:`F_p=2`). Since +:math:`M_{wpd}` is already a moment magnitude, no further :math:`M_w` conversion +is applied (SeisComP's per-station ``estimateMw`` step is the identity here); +the network magnitude is a robust (median / trimmed-mean) average of the station +:math:`M_{wpd}` values. + +Calibration +=========== + +The implementation was validated against GCMT moment magnitudes for 972 events +(:math:`M_w` 6.0--8.8, 2015--2026), processed through this algorithm with +broad-band stations from global networks. :math:`M_{wpd}` tracks +:math:`M_w^{\mathrm{CMT}}` along the 1:1 line with **no saturation** and no +significant depth dependence. + +.. figure:: mwpd_calibration.png + :width: 16cm + :align: center + + :math:`M_{wpd}` (this plugin) versus GCMT :math:`M_w` for 972 events, + coloured by source depth. Overall slope 1.02, mean difference +0.05, + :math:`\sigma = 0.15`, correlation r = 0.95. + +Mean difference :math:`M_{wpd} - M_w^{\mathrm{CMT}}` by magnitude band: + +============= ====== ===== ==== +:math:`M_w` mean sigma N +============= ====== ===== ==== +6.0 -- 7.0 +0.07 0.14 831 +7.0 -- 7.5 -0.05 0.13 94 +7.5 -- 8.0 -0.04 0.20 40 +8.0 -- 9.5 -0.09 0.16 7 +============= ====== ===== ==== + +and by source depth: + +================ ====== ===== ==== +depth (km) mean sigma N +================ ====== ===== ==== +0 -- 70 +0.07 0.16 702 +70 -- 300 -0.02 0.13 175 +300 -- 800 +0.01 0.10 95 +================ ====== ===== ==== + +These results are consistent with operational Early-est (slope ~1.01, +:math:`\pm 0.13`) and with Lomax & Michelini (2009), which reports +:math:`M_{wpd}` matching :math:`M_w^{\mathrm{CMT}}` within :math:`\pm 0.2`. + +Configuration +============= + +Add ``mwpd`` to the ``plugins`` parameter and enable the ``Mwpd`` amplitude and +magnitude in :ref:`scamp` / :ref:`scmag` (or :ref:`scolv`). The defaults +reproduce Early-est; the most important amplitude parameter is +``amplitudes.Mwpd.highpassCorner`` (0.005 Hz --- keep low). See the binding +parameter descriptions for the full list. Because :math:`M_{wpd}` integrates +over the source duration, it requires a long post-*P* window (up to +``maxDuration``), so the final value has higher latency than :math:`M_{wp}`. +For faster tsunami early warning, Lomax & Michelini (2013) describe +:math:`M_{wpd}\mathrm{(RT)}`, a real-time variant that uses a smaller minimum +station distance and avoids some of the duration estimation and re-scaling; +it is not implemented in this plugin. + +References +========== + +* Lomax, A. & Michelini, A. (2009). :math:`M_{wpd}`: a duration--amplitude + procedure for rapid determination of earthquake magnitude and tsunamigenic + potential from *P* waveforms. *Geophys. J. Int.* **176**, 200--214, + doi:10.1111/j.1365-246X.2008.03974.x +* Lomax, A. & Michelini, A. (2013). Tsunami early warning within five minutes. + *Pure Appl. Geophys.* **170**, 1385--1395, doi:10.1007/s00024-012-0512-6 +* Tsuboi, S. et al. (1995). Rapid determination of :math:`M_w` from broadband + *P* waveforms. *Bull. Seismol. Soc. Am.* **85**, 606--613. diff --git a/plugins/magnitudes/mwpd/descriptions/global_mwpd.xml b/plugins/magnitudes/mwpd/descriptions/global_mwpd.xml new file mode 100644 index 000000000..1fa6de7c9 --- /dev/null +++ b/plugins/magnitudes/mwpd/descriptions/global_mwpd.xml @@ -0,0 +1,107 @@ + + + + global + + Duration-amplitude moment magnitude Mwpd (Lomax & Michelini 2009), + ported from Early-est 1.2.9. The vertical broadband is restituted to + displacement, high-passed at the BRB-HP corner and integrated over the + source duration T0; T0 is estimated from the high-frequency envelope. + The amplitude processor emits the displacement integral (nm*s) and + carries T0 as the amplitude period; the magnitude processor forms M0 + and Mwpd with the distance, PREM depth and duration corrections. + + + + + + Measurement parameters for the Mwpd amplitude (displacement + integral over the source duration T0). + + + BRB-HP high-pass corner for the displacement trace (Early-est filter_hp_bu_co_0_005_n_2). Keep low (~0.005 Hz) to preserve the long-period moment of great events. + + + Butterworth order of the BRB-HP high-pass filter. + + + Lead before P at which the displacement integration starts (START_ANALYSIS_BEFORE_P_BRB_HP). + + + Low edge of the high-frequency band used for the T0 envelope. + + + High edge of the high-frequency band used for the T0 envelope. + + + Butterworth order of the high-frequency band-pass. + + + Boxcar half-width that smooths the squared HF trace into the envelope (SMOOTHING_WINDOW_HALF_WIDTH_T50). + + + Earliest the T0 search may terminate (MIN_MWPD_DUR). + + + MAX_MWPD_DUR: cap on T0 and on the post-pick signal/integration window. + + + Envelope signal-to-noise ratio below which the T0 search terminates (SIGNAL_TO_NOISE_RATIO_T0_END). + + + Envelope-to-peak ratio below which the T0 search terminates (SIGNAL_TO_PEAK_RATIO_T0_END). + + + If true, use a fixed source duration instead of the HF-envelope estimate. + + + The fixed source duration used when fixedDuration is true. + + + Cap the displacement-integration window at the theoretical S-P time, as Early-est does, so the integral excludes S/surface-wave energy. The reported source duration T0 (and its ramp correction) is not capped. + + + Travel-time table backend used to compute the S-P cap. + + + Travel-time model used to compute the S-P cap. + + + + + + Moment-magnitude parameters and corrections for Mwpd. + + Density in the Tsuboi moment constant. + + + P velocity in the Tsuboi moment constant. + + + FP average radiation factor in the Tsuboi moment constant. + + + Apply the INGV_EE epicentral-distance correction. + + + Apply the PREM step depth correction (0-600 km). + + + Apply the >90 s duration ramp correction that boosts great/slow events (USE_MPWD_CORR). + + + Reject station magnitudes closer than this epicentral distance. + + + Reject station magnitudes beyond this epicentral distance. + + + Linear station correction a in a*Mwpd+b. Leave 1.0 for none. + + + Constant station correction b in a*Mwpd+b. Leave 0.0 for none. + + + + + diff --git a/plugins/magnitudes/mwpd/descriptions/mwpd_calibration.png b/plugins/magnitudes/mwpd/descriptions/mwpd_calibration.png new file mode 100644 index 000000000..0ef332a17 Binary files /dev/null and b/plugins/magnitudes/mwpd/descriptions/mwpd_calibration.png differ diff --git a/plugins/magnitudes/mwpd/magnitude.cpp b/plugins/magnitudes/mwpd/magnitude.cpp new file mode 100644 index 000000000..40077dc68 --- /dev/null +++ b/plugins/magnitudes/mwpd/magnitude.cpp @@ -0,0 +1,134 @@ +/*************************************************************************** + * SeisComP duration-amplitude moment-magnitude plugin (mwpd) * + * * + * Magnitude processor: displacement integral over T0 -> M0 -> Mwpd, with * + * the INGV_EE distance correction, the PREM step depth correction and the * + * >90 s duration ramp. Faithful to Early-est calculate_Raw_Mwpd_Mag + * + * calculate_corrected_Mwpd_Mag. * + * * + * GNU Affero General Public License Usage - see LICENSE. * + ***************************************************************************/ + + +#define SEISCOMP_COMPONENT Mwpd + +#include + +#include +#include + +#include "mwpd.h" + + +namespace Seiscomp { +namespace Magnitudes { +namespace Mwpd { + + +REGISTER_MAGNITUDEPROCESSOR(MagnitudeProcessor_Mwpd, MWPD_TYPE); + + +MagnitudeProcessor_Mwpd::MagnitudeProcessor_Mwpd() +: Processing::MagnitudeProcessor(MWPD_TYPE) { + _minimumDistanceDeg = _cfg.minDistanceDeg; + _maximumDistanceDeg = _cfg.maxDistanceDeg; +} + + +bool MagnitudeProcessor_Mwpd::setup(const Processing::Settings &settings) { + if ( !Processing::MagnitudeProcessor::setup(settings) ) { + return false; + } + + _cfg = MwpdConfig(); + if ( !readMwpdConfig(settings, std::string("magnitudes.") + type(), _cfg) ) { + return false; + } + + _minimumDistanceDeg = _cfg.minDistanceDeg; + _maximumDistanceDeg = _cfg.maxDistanceDeg; + + SEISCOMP_DEBUG("%s: dist=%.1f-%.1f deg distCorr=%d depthCorr=%d durRamp=%d", + type().c_str(), _cfg.minDistanceDeg, _cfg.maxDistanceDeg, + _cfg.useDistanceCorr, _cfg.useDepthCorr, _cfg.useDurationRamp); + return true; +} + + +MagnitudeProcessor_Mwpd::Status +MagnitudeProcessor_Mwpd::computeMagnitude( + double amplitude, const std::string &unit, + double period, double /*snr*/, + double delta, double depth, + const DataModel::Origin *, + const DataModel::SensorLocation *, + const DataModel::Amplitude *, + const Locale *, + double &value) { + + if ( amplitude <= 0.0 ) { + return AmplitudeOutOfRange; + } + + // The displacement integral is carried in nm*s (as Mwp carries nm). + double ampNmS = amplitude; + if ( !convertAmplitude(ampNmS, unit, MWPD_AMP_UNIT) ) { + return InvalidAmplitudeUnit; + } + + // T0 (source duration) arrives in the amplitude period (seconds). + const double duration = period; + + if ( depth < 0.0 ) { + depth = 0.0; + } + + // --- raw Mwpd (calculate_Raw_Mwpd_Mag) -------------------------------- + const double amplitudeIntegral = ampNmS * 1.0e-9; // nm*s -> m*s + double moment = mwpdMomentConstant(_cfg) * amplitudeIntegral * delta; + if ( moment <= FLT_MIN ) { + return Error; + } + + double raw = (2.0 / 3.0) * (std::log10(moment) - 9.1); + if ( _cfg.useDistanceCorr ) { + raw -= mwpdDistanceCorrection(delta, depth); + } + + // --- corrected Mwpd (calculate_corrected_Mwpd_Mag) -------------------- + double corr = raw; + + if ( _cfg.useDepthCorr ) { + const double dc = mwpdDepthCorrection(depth); + if ( dc > -998.0 ) { // -999 marks an invalid (too-shallow) lookup + corr += dc; + } + } + + // Moment scaling, Early-est calculate_corrected_Mwpd_Mag (USE_MPWD_CORR): + // DURATION-gated ramp -- only for T0 > 90 s (great/slow events), ramped + // linearly 90->110 s, then add (raw - 7.2) * 0.45. (Early-est's pure + // magnitude-gated form, eq 5a, is USE_MOMENT_CORR and is disabled there.) + if ( _cfg.useDurationRamp && duration > _cfg.durRampLow ) { + double ramp = (duration - _cfg.durRampLow) + / (_cfg.durRampHigh - _cfg.durRampLow); + if ( ramp > 1.0 ) { + ramp = 1.0; + } + corr += ramp * (raw - _cfg.magCutoff) * _cfg.rampPow; + } + + value = corr; + + SEISCOMP_DEBUG("%s: ampInt=%g nm*s T0=%.1fs delta=%.2f depth=%.1f " + "-> logM0=%.2f raw=%.2f Mwpd=%.2f", + type().c_str(), ampNmS, duration, delta, depth, + std::log10(moment), raw, value); + + return OK; +} + + +} +} +} diff --git a/plugins/magnitudes/mwpd/mwpd.h b/plugins/magnitudes/mwpd/mwpd.h new file mode 100644 index 000000000..fcbaae05a --- /dev/null +++ b/plugins/magnitudes/mwpd/mwpd.h @@ -0,0 +1,189 @@ +/*************************************************************************** + * SeisComP duration-amplitude moment-magnitude plugin (mwpd) * + * * + * Per-station moment magnitude Mwpd (Lomax & Michelini 2009), porting * + * Early-est 1.2.9. The displacement is integrated over the source * + * duration T0 (estimated from the high-frequency envelope) rather than * + * taking the running-integral peak as Mwp does. * + * * + * GNU Affero General Public License Usage - see LICENSE. * + ***************************************************************************/ + + +#ifndef SEISCOMP_MAGNITUDES_MWPD_PLUGIN_H +#define SEISCOMP_MAGNITUDES_MWPD_PLUGIN_H + + +#include +#include +#include +#include +#include + + +namespace Seiscomp { +namespace Magnitudes { +namespace Mwpd { + + +//! Registered amplitude/magnitude type string. +#define MWPD_TYPE "Mwpd" + +//! Amplitude unit: the displacement integral over T0, carried in nm*s (as Mwp). +#define MWPD_AMP_UNIT "nm*s" + + +// All of this plugin's own symbols are used only inside the plugin (the loader +// needs only createSCPlugin; the processors register via static factories), so +// keep them out of the shared library's dynamic symbol table. Single-TU helpers +// live in anonymous namespaces in the .cpp files; symbols shared across object +// files (the config reader, the corrections, the processor classes) get hidden +// visibility here. +#if defined(__GNUC__) || defined(__clang__) + #define MWPD_LOCAL __attribute__((visibility("hidden"))) +#else + #define MWPD_LOCAL +#endif + + +/** + * Configuration shared by the amplitude and magnitude processors. Defaults + * reproduce the hard-coded Early-est 1.2.9 constants + * (timedomain_processing.c / timedomain_processing_data.c). + */ +struct MwpdConfig { + // --- BRB-HP restitution + displacement integral (amplitude side) ------ + // Early-est high-passes the displacement trace with filter_hp_bu_co_0_005_n_2 + // = Butterworth HP 0.005 Hz, order 2 (mkfilter -Bu -Hp -o 2 -a 2.5e-4). + // (MWPD_GAIN_FREQUENCY=0.15 is only the gain-reference frequency, NOT a + // filter. A 0.15 Hz HP removes periods >7 s and guts great-event moment.) + double highpassCorner = 0.005; //!< [Hz] BRB-HP corner (filter_hp_bu_co_0_005_n_2) + int hpOrder = 2; //!< Butterworth high-pass order for the BRB-HP trace + double analysisPreP = 1.0; //!< [s] START_ANALYSIS_BEFORE_P_BRB_HP + + // --- HF envelope (T0 estimator) --------------------------------------- + // Early-est uses filter_bp_bu_co_1_5_n_4 = Butterworth band-pass 1-5 Hz o4. + double hfFmin = 1.0; //!< [Hz] envelope band low edge + double hfFmax = 5.0; //!< [Hz] envelope band high edge + int hfOrder = 4; //!< envelope band-pass order + double smoothHalfWidth = 5.0; //!< [s] envelope boxcar smoothing half-width + + // --- T0 source-duration estimation ------------------------------------ + double minDuration = 20.0; //!< [s] MIN_MWPD_DUR: earliest the T0 search may terminate + double maxDuration = 1200.0; //!< [s] MAX_MWPD_DUR: T0 search / integration cap + double durationFloor = 0.2; //!< [s] MINIMUN_DURATION + double snT0End = 3.0; //!< SIGNAL_TO_NOISE_RATIO_T0_END + double peakRatioT0End = 0.025; //!< SIGNAL_TO_PEAK_RATIO_T0_END + bool fixedDuration = false; //!< use a fixed T0 instead of the HF-envelope estimate + double fixedDurationVal= 60.0; //!< [s] the fixed T0 when fixedDuration is true + + // --- magnitude formula + corrections ---------------------------------- + double rho = 3400.0; //!< [kg/m^3] density in MWP_CONST + double vp = 7900.0; //!< [m/s] P velocity in MWP_CONST + double fp = 2.0; //!< FP average radiation factor in MWP_CONST + bool useDistanceCorr = true; //!< apply calculate_Mwp_correction (INGV_EE) + bool useDepthCorr = true; //!< apply get_depth_corr_mwpd_prem (PREM step) + bool useDurationRamp = true; //!< USE_MPWD_CORR: >90 s duration ramp + double durRampLow = 90.0; //!< [s] MWPD_MOMENT_CORR_DUR_CUTOFF_LOW + double durRampHigh = 110.0; //!< [s] MWPD_MOMENT_CORR_DUR_CUTOFF_HIGH + double magCutoff = 7.2; //!< MWPD_MOMENT_CORR_MAG_CUTOFF + double rampPow = 0.45; //!< MWPD_MOMENT_CORR_POW + double minDistanceDeg = 5.0; //!< accept station magnitudes from this distance + double maxDistanceDeg = 105.0; //!< ...up to this distance + + // --- S-P window cap (Early-est limits the integral to the S-P time) --- + bool useSpCap = true; //!< cap the displacement integral at the S-P time + std::string tttBackend = "libtau"; //!< travel-time backend for S-P + std::string tttModel = "iasp91"; //!< travel-time model for S-P +}; + + +/** Reads MwpdConfig from a binding. @p prefix is "amplitudes.Mwpd" or + * "magnitudes.Mwpd". Missing keys keep their (Early-est) defaults. */ +MWPD_LOCAL bool readMwpdConfig(const Processing::Settings &settings, + const std::string &prefix, MwpdConfig &out); + + +/** Tsuboi (1995) moment constant, reproduced from Early-est: + * 4 pi * rho * Vp^3 * FP * (10000/90) deg->km * 1000 km->m. */ +MWPD_LOCAL double mwpdMomentConstant(const MwpdConfig &cfg); + +/** INGV_EE distance correction (subtracted from the raw magnitude), + * calculate_Mwp_correction_INGV_EE: 0 if delta<0 or depth>100 km. */ +MWPD_LOCAL double mwpdDistanceCorrection(double deltaDeg, double depthKm); + +/** PREM step depth correction get_depth_corr_mwpd_prem (returns -999 if + * the depth is below the first table row). */ +MWPD_LOCAL double mwpdDepthCorrection(double depthKm); + + +// --------------------------------------------------------------------------- +// Amplitude processor (registered as "Mwpd"). Vertical broadband only. From +// the raw window it builds (a) the BRB-HP displacement and its running +// integral, separated into positive/negative displacement lobes, and (b) the +// high-frequency envelope from which the source duration T0 is estimated. +// Emits the displacement integral over T0 (nm*s) and carries T0 as the +// amplitude "period". +// --------------------------------------------------------------------------- +class MWPD_LOCAL AmplitudeProcessor_Mwpd : public Processing::AmplitudeProcessor { + public: + AmplitudeProcessor_Mwpd(); + + public: + bool setup(const Processing::Settings &settings) override; + const DoubleArray *processedData(Component comp) const override; + + protected: + bool computeAmplitude(const DoubleArray &data, + size_t i1, size_t i2, + size_t si1, size_t si2, + double offset, + AmplitudeIndex *dt, + AmplitudeValue *amplitude, + double *period, double *snr) override; + + private: + void applyConfig(); + //! Theoretical S-P time [s] at the current source/receiver, or -1. + double computeSPSeconds() const; + + private: + MwpdConfig _cfg; + DoubleArray _processedData; + Seiscomp::TravelTimeTableInterfacePtr _ttt; +}; + + +// --------------------------------------------------------------------------- +// Magnitude processor (registered as "Mwpd"). Turns the displacement +// integral into M0 and Mwpd, applying the distance, PREM depth and duration +// corrections. The source duration T0 arrives in the amplitude "period". +// --------------------------------------------------------------------------- +class MWPD_LOCAL MagnitudeProcessor_Mwpd : public Processing::MagnitudeProcessor { + public: + MagnitudeProcessor_Mwpd(); + + public: + void setDefaults() override {} + bool setup(const Processing::Settings &settings) override; + + Status computeMagnitude(double amplitude, const std::string &unit, + double period, double snr, + double delta, double depth, + const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Amplitude *, + const Locale *, + double &value) override; + + private: + MwpdConfig _cfg; +}; + + +} +} +} + + +#endif diff --git a/plugins/magnitudes/mwpd/plugin.cpp b/plugins/magnitudes/mwpd/plugin.cpp new file mode 100644 index 000000000..6da853919 --- /dev/null +++ b/plugins/magnitudes/mwpd/plugin.cpp @@ -0,0 +1,22 @@ +/*************************************************************************** + * SeisComP duration-amplitude moment-magnitude plugin (mwpd) * + * * + * Plugin entry point. Registers the Mwpd amplitude and magnitude * + * processors (a port of Early-est 1.2.9's Mwpd, Lomax & Michelini 2009). * + * * + * GNU Affero General Public License Usage - see LICENSE. * + ***************************************************************************/ + + +#include + +#include "version.h" + + +ADD_SC_PLUGIN( + "Duration-amplitude moment magnitude Mwpd: per-station moment magnitude " + "from the integral of the displacement over the source duration T0 " + "(port of Early-est 1.2.9, Lomax & Michelini 2009).", + "Early-est port", + MWPD_VERSION_MAJOR, MWPD_VERSION_MINOR, MWPD_VERSION_PATCH +) diff --git a/plugins/magnitudes/mwpd/version.h b/plugins/magnitudes/mwpd/version.h new file mode 100644 index 000000000..9d3dd77b0 --- /dev/null +++ b/plugins/magnitudes/mwpd/version.h @@ -0,0 +1,14 @@ +/*************************************************************************** + * SeisComP duration-amplitude moment-magnitude plugin (mwpd) * + ***************************************************************************/ + +#ifndef SEISCOMP_MAGNITUDES_MWPD_VERSION_H +#define SEISCOMP_MAGNITUDES_MWPD_VERSION_H + +#define MWPD_VERSION_MAJOR 0 +#define MWPD_VERSION_MINOR 1 +#define MWPD_VERSION_PATCH 0 + +#define MWPD_VERSION "0.1.0" + +#endif