Z-Score Filtering for Vibration Anomalies in IIoT Telemetry

Vibration telemetry from rotating assets — motors, gearboxes, spindles, pumps — is the primary signal layer for predictive maintenance and real-time OEE optimization. This page is the concrete recipe under outlier detection methods for one specific job: separating genuine mechanical anomalies from baseline operational noise using a rolling Z-score, without drowning the CMMS in false tickets. Static amplitude thresholds fracture under variable load and seasonal process shifts; a dynamic Z-score continuously recalibrates against the local statistical distribution, so only statistically significant deviations propagate downstream to alerting, work-order dispatch, or digital-twin synchronization. The hard part is not the formula — it is making the formula honest across machine states, async batch boundaries, and drifting edge clocks.

The core calculation uses a sliding temporal window to compute the local mean ($\mu_t$) and standard deviation ($\sigma_t$) of acceleration or velocity RMS at each sampling interval. The instantaneous Z-score is:

$Z_t = \frac{x_t - \mu_t}{\sigma_t}$

A score beyond ±3.0 during steady-state operation typically signals incipient bearing degradation, rotor unbalance, or shaft misalignment. The same value during a startup transient, coast-down, or rapid load shift is expected process dynamics — which is why raw Z-score thresholding alone is an anti-pattern on the factory floor.

Variant 1 — State-aware rolling Z-score (the default) #Permalink to this section

The standard production implementation couples the rolling statistic to machine-state context and pauses recalibration outside operational states. Without this masking, transient stress during ramp-up triggers phantom flags that corrupt availability and inflate false-positive maintenance tickets. Machine state should come from canonical PLC tags — enforce PLC tag standardization upstream so RUNNING/STEADY_STATE mean the same thing on every line.

import numpy as np
import pandas as pd


def state_aware_zscore(
    df: pd.DataFrame,
    value_col: str = "vibration_rms",
    state_col: str = "machine_state",
    window: str = "15min",
    z_threshold: float = 3.0,
) -> pd.DataFrame:
    """Rolling Z-score on operational segments only, with defensive denominator."""
    df = df.sort_index()
    if not df.index.is_monotonic_increasing:
        raise ValueError("Timestamp index must be strictly monotonic.")

    # Gap-fill before aggregation so NaN runs do not collapse the denominator.
    df[value_col] = df[value_col].interpolate(method="time", limit=30, limit_direction="both")
    df[value_col] = df[value_col].ffill(limit=60)

    # Mask: compute rolling stats only while the asset is actually turning.
    op_mask = df[state_col].isin(["RUNNING", "STEADY_STATE"])
    roll = df.loc[op_mask, value_col].rolling(window=window, min_periods=8)
    mu = roll.mean()
    sigma = roll.std(ddof=0).clip(lower=1e-4)  # clamp flat-signal std -> no inf

    df["z_score"] = np.nan
    df.loc[op_mask, "z_score"] = (df.loc[op_mask, value_col] - mu) / sigma
    df["is_anomaly"] = df["z_score"].abs() > z_threshold
    return df

Three production details are load-bearing here. The fill step reuses the gap-filling algorithms stage — short interior holes are bridged with linear interpolation for missing values before any statistic runs, because a NaN run silently drives the rolling denominator toward zero. ddof=0 (population standard deviation) matches the streaming assumption that the window is the observed distribution. And sigma.clip(lower=1e-4) stops a perfectly flat signal — common during idle or sensor drift — from producing inf; note that because IEEE 754 doubles cannot represent decimals exactly, this clamp is also why threshold comparisons must use a tolerance, a caveat covered in floating-point drift in sensor readings.

Variant 2 — Robust (modified) Z-score with MAD #Permalink to this section

The classic Z-score has a self-defeating flaw: the mean and standard deviation it depends on are themselves dragged by the very outliers it is trying to catch. A single accelerometer spike inflates $\sigma_t$, masking the next several genuine events. For spiky vibration data, prefer the modified Z-score built on the median and Median Absolute Deviation (MAD), which has a ~50% breakdown point:

$M_t = \frac{0.6745\,(x_t - \tilde{x}_t)}{\text{MAD}_t}$

def modified_zscore(
    df: pd.DataFrame,
    value_col: str = "vibration_rms",
    window: str = "15min",
    m_threshold: float = 3.5,
) -> pd.DataFrame:
    """MAD-based Z-score: resistant to the spikes a mean/std filter masks."""
    x = df[value_col]
    median = x.rolling(window=window, min_periods=8).median()
    abs_dev = (x - median).abs()
    mad = abs_dev.rolling(window=window, min_periods=8).median().clip(lower=1e-4)

    df["mod_z"] = 0.6745 * (x - median) / mad  # 0.6745 = Phi^-1(0.75)
    df["is_anomaly"] = df["mod_z"].abs() > m_threshold
    return df

Use 3.5 as the conventional cutoff (Iglewicz–Hoaglin). MAD costs more per window than mean/std but is the right default when the signal carries impulsive faults — gear-mesh chips, bearing spalls — rather than smooth drift.

Variant 3 — Carrying rolling state across async batches #Permalink to this section

Vibration arrives in chunks (for example 5-minute edge buffers fanned out through async batch processing). If each batch recomputes its window from scratch, $\sigma_t$ is discontinuous at every chunk boundary and the first samples of each batch spike falsely. The fix is to serialize the trailing window and carry it forward into the next worker — a pattern that pairs naturally with a Celery-based MQTT ingestion topology.

from dataclasses import dataclass


@dataclass
class WindowState:
    """Trailing samples persisted between batches (e.g. in Redis)."""
    tail: pd.Series          # last `window` of values, indexed by timestamp
    last_ts: pd.Timestamp


def zscore_batch(batch: pd.Series, state: WindowState, window: str = "15min",
                 z_threshold: float = 3.0) -> tuple[pd.Series, WindowState]:
    if batch.index.min() <= state.last_ts:
        raise ValueError("Batch overlaps prior state; dedupe before scoring.")

    stitched = pd.concat([state.tail, batch]).sort_index()
    roll = stitched.rolling(window=window, min_periods=8)
    z = (stitched - roll.mean()) / roll.std(ddof=0).clip(lower=1e-4)

    z_batch = z.loc[batch.index]               # only emit the new rows
    new_tail = stitched.last(window)           # persist trailing window
    return z_batch.abs() > z_threshold, WindowState(new_tail, batch.index.max())

The overlap guard matters: replays and at-least-once MQTT delivery duplicate samples, and a duplicate inside the stitched window depresses variance and hides the next anomaly.

Gotchas & anti-patterns #Permalink to this section

• Raw thresholding without state masking. Flagging every ±3σ breach counts every startup transient as a fault. Gate on PLC state and add hysteresis to the RUNNING/IDLE boundary so brief dips do not flap the mask.
• Stateless windows across batches. Recomputing per chunk produces a reliable spike on the first samples of every batch. Carry the trailing window forward (Variant 3).
• Mean/std on impulsive signals. One spike inflates σ and masks the next few real events; switch to the MAD variant for bearing/gear faults.
• Scoring before clock sync. Millisecond skew between edge gateways and the historian misaligns window boundaries and inflates σ. Run clock drift correction (NTP/PTP at the edge) before any rolling statistic.
• Treating flags as ground truth for OEE. A flagged sample is a candidate, not downtime. Feed it through OEE formula validation and microstop logic, not straight into availability math.

Quick reference: choosing a variant and threshold #Permalink to this section

Common failure modes and fixes:

For the underlying statistics, the NIST/SEMATECH e-Handbook of Statistical Methods is the standard reference on detecting outliers in measurement data.

Related #Permalink to this section

• Outlier Detection Methods — parent overview of the layered detection stack
• Implementing Linear Interpolation for Missing Sensor Values — the fill stage this filter depends on
• Clock Drift Correction — the prerequisite for stable rolling windows
• Using Celery for High-Throughput MQTT Ingestion — the async topology behind batch state carry-forward
• Adjusting Performance Thresholds for Legacy CNC Machines — turning flags into OEE-grade microstop decisions