diff --git a/docs/ri8507_compliance_curve.md b/docs/ri8507_compliance_curve.md new file mode 100644 index 0000000..befb1a3 --- /dev/null +++ b/docs/ri8507_compliance_curve.md @@ -0,0 +1,135 @@ +# USBM RI8507 / OSMRE Blasting Compliance Curve — Reference + +Reference for the **velocity-vs-frequency blasting compliance chart** Blastware +draws on its Event Report ("USBM RI8507 And OSMRE"), and how seismo-relay +reproduces it. Implemented in [`sfm/compliance.py`](../sfm/compliance.py); the +spectral (FFT) side lives in [`waveform_fft.py`](../waveform_fft.py). + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each +with a Blastware Event Report + FFT Report as ground truth. Curve values from +USBM RI8507 Appendix B and 30 CFR 816.67. + +--- + +## What it is + +Two closely-related sources for the same limit curve: + +- **USBM RI8507** — Bureau of Mines *Report of Investigations 8507* (Siskind + et al., 1980), *"Structure Response and Damage Produced by Ground Vibration + From Surface Mine Blasting."* The curve is **Figure B-1**, Appendix B + ("Alternative Blasting Level Criteria"), p.73–74. +- **OSMRE / OSM** — the Office of Surface Mining Reclamation and Enforcement + codified it as **30 CFR 816.67, Figure 1**. "CFR" = the U.S. Code of Federal + Regulations. Same curve, regulatory force. + +The chart plots each geophone channel's significant vibration cycles as +`(frequency, peak velocity)` points against this limit. A point **below** the +line passes; **above** fails. + +--- + +## The limit curve + +A structure has a resonance band (~4–12 Hz for whole structures) where it is +most vulnerable, so the safe velocity is **lower** at those frequencies and +**higher** away from them. The curve captures this by alternating two kinds of +bound: + +- **Constant-velocity** segments — a flat horizontal line at a fixed PPV. +- **Constant-displacement** segments — a fixed peak *displacement* `d`. For + simple harmonic motion, peak velocity `v = 2πf·d`, so on a velocity-vs- + frequency **log-log** plot this is a straight line of slope +1 (velocity rises + with frequency). This is why the low- and high-frequency bounds are sloped. + +### Two lines — structure type + +RI8507 gives two lines for two interior-wall constructions (Table 13, p.67): + +| line | construction | plateau PPV | +|---|---|---| +| **Drywall** (solid) | modern gypsum wallboard | **0.75 in/s** | +| **Plaster** (dashed) | older plaster on wood lath | **0.50 in/s** | + +Plaster-on-lath is more damage-prone, hence the lower limit. You apply **one** +line depending on the monitored structure. + +### The four segments (Figure B-1, p.74) + +Going low → high frequency, each line is: + +1. **Ultimate low-frequency bound** — constant displacement **0.030 in** + (`v = 2πf·0.030`). Only relevant below ~4 Hz. +2. **Plateau** — constant velocity **0.75** (Drywall) / **0.50** (plaster) in/s. +3. **Rising diagonal** — constant displacement **0.008 in** (`v = 2πf·0.008`), + climbing from the plateau up to the high-frequency cap. +4. **High-frequency cap** — constant velocity **2.0 in/s** above ~40 Hz. + +The segments are drawn **continuous**: each bound is used over the frequency +range where it is the binding (lowest) limit, and consecutive bounds meet where +they are equal — so there are no vertical steps. Transition frequencies come +straight from the values (`f = V / (2π·d)`): + +| transition | formula | Drywall | Plaster | +|---|---|---|---| +| 0.030 in → plateau | `V_mid / (2π·0.030)` | 3.98 Hz | 2.65 Hz | +| plateau → 0.008 in | `V_mid / (2π·0.008)` | 14.92 Hz | 9.95 Hz | +| 0.008 in → 2.0 in/s | `2.0 / (2π·0.008)` | 39.79 Hz | 39.79 Hz | + +Because both lines share the same **0.008 in** rising diagonal, above ~15 Hz +they lie on the *same* line (both reach 2.0 in/s at ~40 Hz) — RI8507's literal +construction merges them there. Blastware renders the dashed line as a separate +parallel diagonal, but that is cosmetic: above ~15 Hz both structure types carry +the identical limit, so compliance is unaffected. + +> ⚠ RI8507's *Table 13* is a simpler two-range criterion with a **sharp +> discontinuity at 40 Hz** (flat plateau, then a jump to 2.0). Figure B-1 is the +> **smoothed** version that adds the 0.008 in transition — that is the one drawn +> on reports and implemented here. + +--- + +## The compliance scatter (the points) + +The cloud is **not** the FFT spectrum. It is a per-cycle, time-domain measure by +the **zero-crossing method** (`channel_compliance_points`): + +- Split the channel's waveform at its zero crossings. +- Each half-cycle contributes one point: **frequency** `= 1 / (2 · half-period)` + (from the samples between the two crossings), **velocity** `= peak |amplitude|` + in that half-cycle. + +This yields ~90–110 points per channel, and — by construction — each channel's +**highest** point equals that channel's PPV. Verified against Blastware: the +cloud shape, density, and ceiling all match. + +### Why not the FFT? + +A broadband blast spreads its energy across many FFT bins, so no single bin +reaches the time-domain peak — the FFT amplitudes come out ~10× below the +compliance-chart velocities. The compliance chart is a *per-cycle peak* view; +the **FFT** is a separate analysis (Blastware's *FFT Report*), reproduced by +[`waveform_fft.py`](../waveform_fft.py) and used for the dominant-frequency +readout and the #10 FFT view — not for this scatter. + +--- + +## Implementation + +- `sfm/compliance.py` + - `limit_at(freq, curve)` — the limit PPV at a frequency (`curve` = `"Drywall"` + or `"Plaster"`); curves are data in `_CURVES`, so more standards can be added. + - `channel_compliance_points(samples, sps)` — the zero-crossing scatter. + - `draw_compliance_chart(ax, channels, sps)` — matplotlib rendering (both + limit lines + per-channel scatter, Blastware's tick scales and channel + markers: Tran `+` red, Vert `×` green, Long `o` blue). +- Tests: `tests/test_compliance.py`. + +--- + +## Sources + +- USBM **RI8507** (Siskind, Stagg, Kopp, Dowding, 1980), Appendix B / Figure B-1, + p.73–74; Table 13, p.67. (`ref-stuff/usbm-ri8507-ground_vibration.pdf`.) +- **30 CFR 816.67**, "Use of explosives: Control of adverse effects," Figure 1 — + diff --git a/sfm/compliance.py b/sfm/compliance.py new file mode 100644 index 0000000..2b21585 --- /dev/null +++ b/sfm/compliance.py @@ -0,0 +1,132 @@ +"""USBM RI8507 / OSMRE blasting compliance chart. + +Renders the velocity-vs-frequency compliance scatter Blastware draws on its Event +Report: each channel's significant waveform cycles as ``(frequency, peak +velocity)`` points on log-log axes against the regulatory limit curve(s). A point +below the curve passes; above fails. + +Two pieces, kept separate so both can be reused/extended: + * ``limit_at`` / ``limit_curve`` — the regulatory limit curve(s), as data. + * ``channel_compliance_points`` — the per-cycle (freq, velocity) scatter, by + the zero-crossing method (matches Blastware: each channel's cloud tops out + at that channel's PPV). + +Limit curves (USBM RI8507 Figure B-1 / OSM 30 CFR 816.67), drawn CONTINUOUS — a +constant-displacement bound (sloped, ``v = 2πf·d``) meets a constant-velocity +plateau at the frequency where they're equal, so there are no vertical steps +(matching how Blastware draws it). Two lines: + * **Drywall** (modern gypsum board) — 0.75 in/s plateau (solid). + * **Plaster** on wood lath (older homes) — 0.50 in/s plateau (dashed). +Both use a 0.030 in low-frequency displacement bound and rise through a 0.010 in +displacement bound to a 2.0 in/s high-frequency plateau. Values from USBM RI8507 +(Appendix B) / 30 CFR 816.67; ⚠ confirm the exact shape against a Blastware +report before trusting for compliance. +""" +from __future__ import annotations + +import math +from typing import Dict, Sequence, Tuple + +import numpy as np +from matplotlib.ticker import FixedLocator, NullLocator + +# curve name → (low-freq "ultimate" displacement in, mid velocity plateau in/s, +# high-freq displacement in, high-freq velocity plateau in/s). +# RI8507 Fig B-1 (p.74): ultimate max displacement 0.030 in (< ~4 Hz), plateau +# 0.75 (Drywall) / 0.50 (plaster), rising diagonal at 0.008 in displacement up to +# a 2.0 in/s plateau reached at ~40 Hz. +_CURVES: Dict[str, Tuple[float, float, float, float]] = { + "Drywall": (0.030, 0.75, 0.008, 2.00), + "Plaster": (0.030, 0.50, 0.008, 2.00), +} +# how each curve is stroked on the chart +_CURVE_STYLE = {"Drywall": {"ls": "-", "lw": 1.0}, "Plaster": {"ls": "--", "lw": 0.9}} + +STANDARDS = tuple(_CURVES) + +# Blastware's channel markers/colours on the compliance chart. +_CHANNEL_STYLE = { + "Tran": ("+", "#d62728"), # red + + "Vert": ("x", "#2ca02c"), # green x + "Long": ("o", "#1f77b4"), # blue o +} + + +def limit_at(freq_hz: float, curve: str = "Drywall") -> float: + """Max allowed PPV (in/s) at ``freq_hz`` for ``curve`` (continuous).""" + d_low, v_mid, d_high, v_high = _CURVES[curve] + f = max(freq_hz, 1.0) + f_a = v_mid / (2.0 * math.pi * d_low) # disp_low → vel_mid + f_b = v_mid / (2.0 * math.pi * d_high) # vel_mid → disp_high + f_c = v_high / (2.0 * math.pi * d_high) # disp_high → vel_high + if f <= f_a: + return 2.0 * math.pi * f * d_low + if f <= f_b: + return v_mid + if f <= f_c: + return 2.0 * math.pi * f * d_high + return v_high + + +def limit_curve(curve: str = "Drywall", fmin: float = 1.0, fmax: float = 100.0, n: int = 400): + """(freqs, limits) sampled across the band for plotting one curve.""" + freqs = np.logspace(np.log10(fmin), np.log10(fmax), n) + return freqs, np.array([limit_at(f, curve) for f in freqs]) + + +def channel_compliance_points( + samples: Sequence[float], sps: float, fmin: float = 1.0, fmax: float = 100.0, + vmin: float = 0.0, +) -> Tuple[np.ndarray, np.ndarray]: + """Per-cycle (frequency, peak velocity) scatter for one channel. + + Zero-crossing method: split the trace at sign changes; each half-cycle + contributes one point at ``(1/(2·half_period), max|amplitude|)``. Matches + Blastware — the cloud's ceiling is the channel PPV. ``samples`` must be in the + velocity unit you want plotted (in/s). Points outside ``[fmin, fmax]`` or at + or below ``vmin`` are dropped. + """ + x = np.asarray(samples, dtype=float) + if x.size < 3: + return np.empty(0), np.empty(0) + zc = np.where(np.diff(np.signbit(x)))[0] + freqs, vels = [], [] + for a, b in zip(zc[:-1], zc[1:]): + half_period = (b - a) / sps + if half_period <= 0: + continue + freqs.append(1.0 / (2.0 * half_period)) + vels.append(float(np.abs(x[a:b + 1]).max())) + f = np.array(freqs) + v = np.array(vels) + keep = (f >= fmin) & (f <= fmax) & (v > vmin) + return f[keep], v[keep] + + +def draw_compliance_chart(ax, channels: Dict[str, Sequence[float]], sps: float) -> None: + """Draw the compliance chart (both limit curves + per-channel scatter).""" + for name, style in _CURVE_STYLE.items(): + cf, cv = limit_curve(name) + ax.plot(cf, cv, color="#333", zorder=3, **style) + + for ch, (marker, color) in _CHANNEL_STYLE.items(): + samples = channels.get(ch) + if samples is None or len(samples) == 0: + continue + f, v = channel_compliance_points(samples, sps) + ax.scatter(f, v, marker=marker, s=12, c=color, linewidths=0.7, zorder=4, label=ch) + + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlim(1, 100) + ax.set_ylim(0.0394, 10) + xt = [1, 2, 5, 10, 20, 50, 100] + yt = [0.0394, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10] + ax.xaxis.set_major_locator(FixedLocator(xt)); ax.xaxis.set_minor_locator(NullLocator()) + ax.yaxis.set_major_locator(FixedLocator(yt)); ax.yaxis.set_minor_locator(NullLocator()) + ax.set_xticklabels([str(v) for v in xt]) + ax.set_yticklabels([("%g" % v) for v in yt]) + ax.set_xlabel("Frequency (Hz)", fontsize=7) + ax.set_ylabel("Velocity (in/s)", fontsize=7) + ax.tick_params(labelsize=6) + ax.grid(True, which="both", ls=":", lw=0.4, color="#ccc") diff --git a/tests/test_compliance.py b/tests/test_compliance.py new file mode 100644 index 0000000..c102e9c --- /dev/null +++ b/tests/test_compliance.py @@ -0,0 +1,35 @@ +"""USBM/OSMRE compliance curve + scatter logic (sfm.compliance). +Rendering is verified visually against Blastware reports.""" +import math + +import numpy as np + +from sfm.compliance import limit_at, channel_compliance_points + + +def test_osmre_velocity_segments(): + assert abs(limit_at(6.0) - 0.75) < 1e-9 # 3.5–12 Hz flat + assert abs(limit_at(50.0) - 2.00) < 1e-9 # 30–100 Hz flat + + +def test_displacement_segments(): + assert abs(limit_at(2.0) - 2 * math.pi * 2.0 * 0.030) < 1e-9 # low-freq 0.030 in + assert abs(limit_at(20.0) - 2 * math.pi * 20.0 * 0.008) < 1e-9 # rising diagonal 0.008 in + + +def test_limit_clamps_below_1hz(): + assert limit_at(0.1) == limit_at(1.0) + + +def test_scatter_ceiling_is_ppv_at_dominant_freq(): + # ~27 Hz blast-like trace whose energy peaks mid-record (inside full cycles, + # as a real event does): the scatter cloud's ceiling is the trace PPV and the + # top point sits near the dominant frequency. + sps, n = 1024.0, 3328 + t = np.arange(n) / sps + env = np.exp(-((t - 1.5) ** 2) / (2 * 0.3 ** 2)) + x = 0.9 * env * np.sin(2 * np.pi * 27.0 * t) + f, v = channel_compliance_points(x, sps) + assert len(f) > 20 + assert v.max() >= 0.99 * np.abs(x).max() + assert 20.0 < f[int(np.argmax(v))] < 35.0