Anomaly Detection in AIS Trajectories #

Operational deployment of anomaly detection within the AIS Vessel Tracking & Route Automation pipeline demands deterministic scoring, strict memory boundaries, and geodesically accurate kinematic profiling. Heuristic thresholding is systematically replaced by state-driven kinematic analysis to ensure reproducible maritime analytics across coastal and offshore jurisdictions. This pipeline ingests raw AIS message streams, computes speed and heading derivatives, flags spatial outliers, and segments behavioral states without materializing full trajectory datasets in memory. The design prioritizes streaming parity, projection-agnostic distance calculations, and deterministic rollback pathways for production maritime systems.

Reference Configuration #

Detection Flow #

Memory-Constrained Ingestion and CRS Validation #

Raw AIS data arrives as high-frequency, irregularly sampled point streams in WGS84 (EPSG:4326). Loading multi-terabyte historical archives or high-throughput real-time feeds into monolithic DataFrames triggers out-of-memory failures and introduces projection distortion during distance calculations. Production workflows must stream partitioned Parquet files, validate coordinate bounds, and enforce strict temporal ordering per MMSI before any kinematic computation begins.

The ingestion layer uses Polars lazy evaluation to maintain sub-2 GB RAM footprints per worker, applying WGS84 bounds validation and monotonic timestamp sorting. This pattern integrates directly with the real-time AIS stream ingestion pipeline by consuming identical partitioned schemas while adding pre-computation guards for downstream scoring. Coordinate validation drops corrupted sensor fills — for example, lat > 90.0 or lon == 0.0 default fills — before they propagate into derivative calculations.

Geodesic Kinematic Derivative Profiling #

Speed and heading anomalies manifest as instantaneous jumps exceeding physical vessel maneuverability constraints. Rather than applying static thresholds, the pipeline computes rolling derivatives using geodesic distance and spherical trigonometry, eliminating false positives caused by coordinate projection artifacts, particularly in high-latitude coastal zones where planar approximations fail.

Using pyproj.Geod, the pipeline calculates forward azimuths and orthodromic distances between consecutive pings. The IMO AIS Performance Standards mandate update rates that vary by vessel class and speed, meaning derivative windows must dynamically adjust to expected reporting intervals. Vectorized computation ensures batch-level throughput without Python loop overhead. Time deltas are computed in seconds, with division-by-zero guards applied before converting orthodromic meters to nautical miles per hour. This kinematic foundation also powers the speed and heading profiling workflow, which handles SOG and COG quality filtering upstream of anomaly scoring.

Production-Grade Implementation #

The following implementation demonstrates a memory-safe, streaming-compatible pipeline. It uses Polars lazy evaluation for partitioned ingestion, vectorized pyproj operations for geodesic profiling, and rolling MAD for deterministic anomaly scoring.

import logging
import polars as pl
from pyproj import Geod
import numpy as np
from typing import Optional

logger = logging.getLogger(__name__)

# Geodesic calculator for WGS84 — avoids projection overhead and planar distortion
GEOD = Geod(ellps="WGS84")

# Anomaly scoring constants
MAD_MULTIPLIER: float = 3.5
ROLLING_WINDOW: str = "24h"
GAP_THRESHOLD_SECONDS: float = 900.0  # 15-minute gap mask boundary


def load_ais_trajectory_batch(
    parquet_path: str,
    mmsi_filter: Optional[list[int]] = None,
) -> pl.LazyFrame:
    """Stream AIS Parquet partitions with strict memory bounds.

    Applies WGS84 coordinate validation and temporal sorting per MMSI
    before any kinematic computation. Drops corrupted sensor fills
    (lat > 90.0, lon == 0.0) that would propagate into derivative calculations.
    """
    lf = pl.scan_parquet(parquet_path)
    if mmsi_filter:
        lf = lf.filter(pl.col("mmsi").is_in(mmsi_filter))
        logger.info("Filtered to %d MMSI identifiers.", len(mmsi_filter))

    # Enforce WGS84 bounds and reject null timestamps
    lf = lf.filter(
        (pl.col("lat").is_between(-90.0, 90.0))
        & (pl.col("lon").is_between(-180.0, 180.0))
        & pl.col("timestamp").is_not_null()
    )
    return lf.sort(["mmsi", "timestamp"])


def _compute_kinematics_chunk(df: pl.DataFrame) -> pl.DataFrame:
    """Vectorized geodesic distance and heading computation per MMSI partition.

    Computes orthodromic distance (meters) and forward azimuth between consecutive
    pings using pyproj.Geod.inv(), then converts to knots and degrees.
    Prepends NaN for the first ping of each partition (no previous position).
    """
    lat: np.ndarray = df["lat"].to_numpy()
    lon: np.ndarray = df["lon"].to_numpy()
    ts: np.ndarray = df["timestamp"].cast(pl.Int64).to_numpy()  # microseconds since epoch

    # Forward azimuth and orthodromic distance (meters) between consecutive pings
    fwd_azim: np.ndarray
    dist_m: np.ndarray
    fwd_azim, _, dist_m = GEOD.inv(lon[:-1], lat[:-1], lon[1:], lat[1:])

    dt_sec: np.ndarray = np.diff(ts) / 1_000_000.0  # microseconds → seconds

    # Guard zero/negative time deltas — they indicate duplicate or out-of-order records
    dt_sec = np.where(dt_sec <= 0, np.nan, dt_sec)

    speed_kn: np.ndarray = (dist_m / dt_sec) * 1.94384  # m/s → knots
    heading_deg: np.ndarray = fwd_azim % 360.0

    # Tag gaps exceeding the 15-minute threshold; derivative accuracy is undefined across them
    gap_mask: np.ndarray = np.concatenate([[False], dt_sec > GAP_THRESHOLD_SECONDS])

    return df.with_columns([
        pl.Series("speed_kn", np.concatenate([[np.nan], speed_kn])),
        pl.Series("heading_deg", np.concatenate([[np.nan], heading_deg])),
        pl.Series("dt_sec", np.concatenate([[np.nan], dt_sec])),
        pl.Series("gap_mask", gap_mask),
    ])


# Scale factor that makes MAD a consistent estimator of the standard deviation
# under a normal distribution: 1 / Φ⁻¹(0.75) ≈ 1.4826. Using scaled MAD keeps the
# 3.5× multiplier interpretable as "robust sigmas" while staying outlier-resistant.
MAD_SCALE: float = 1.4826


def build_anomaly_pipeline(lf: pl.LazyFrame) -> pl.LazyFrame:
    """Chain ingestion, kinematic profiling, and MAD-based scoring.

    Rolling median and Median Absolute Deviation are computed over a 24-hour
    window anchored to UTC timestamps. Points where
    |speed_kn - rolling_median| > 3.5 * (1.4826 * rolling_MAD) are flagged
    ANOMALOUS. Unlike mean/standard-deviation baselines, the median and MAD
    are not pulled toward the very spikes we are trying to detect, so a single
    erroneous fix cannot inflate the baseline and mask neighbouring anomalies.
    Gap-masked rows are excluded from scoring and receive a NULL anomaly_flag
    to prevent stale baselines contaminating post-gap statistics.
    """
    # Apply geodesic kinematics per vessel partition — maintains MMSI ordering
    scored: pl.DataFrame = lf.group_by("mmsi", maintain_order=True).map_groups(
        _compute_kinematics_chunk
    )

    # Rolling median over the 24-hour window keyed to UTC timestamp
    med: pl.Expr = (
        pl.col("speed_kn")
          .rolling_median(window_size=ROLLING_WINDOW, by="timestamp")
          .over("mmsi")
    )

    windowed: pl.DataFrame = scored.with_columns(
        med.alias("speed_rolling_med")
    ).with_columns(
        # Absolute deviation from the rolling median — the per-row input to MAD
        (pl.col("speed_kn") - pl.col("speed_rolling_med")).abs().alias("abs_dev")
    ).with_columns(
        # MAD = median of the absolute deviations within the same rolling window
        pl.col("abs_dev")
          .rolling_median(window_size=ROLLING_WINDOW, by="timestamp")
          .over("mmsi")
          .alias("speed_rolling_mad")
    ).with_columns(
        # Scaled MAD is the robust analogue of the standard deviation
        (pl.col("speed_rolling_mad") * MAD_SCALE).alias("speed_robust_sigma")
    )

    # MAD-based anomaly flag; gap-masked rows receive NULL to preserve auditability
    return windowed.with_columns([
        pl.when(pl.col("gap_mask"))
          .then(pl.lit(None))
          .when(
              (pl.col("speed_kn").is_not_null())
              & (pl.col("speed_robust_sigma") > 0)
              & (
                  pl.col("abs_dev")
                  > MAD_MULTIPLIER * pl.col("speed_robust_sigma")
              )
          )
          .then(pl.lit("ANOMALOUS"))
          .otherwise(pl.lit("NOMINAL"))
          .alias("anomaly_flag"),
        pl.when(
            (pl.col("speed_robust_sigma") > 0) & (~pl.col("gap_mask"))
        )
          .then(
              1.0 - (
                  pl.col("abs_dev")
                  / (MAD_MULTIPLIER * pl.col("speed_robust_sigma") + 1e-9)
              ).clip(0.0, 1.0)
          )
          .otherwise(pl.lit(None))
          .alias("kinematic_confidence"),
    ])

Statistical Outlier Flagging and Behavioral Segmentation #

Once kinematic derivatives are computed, anomalies are isolated using robust statistical methods resistant to heavy-tailed maritime noise. Median Absolute Deviation and rolling z-scores replace mean/variance baselines, preventing port congestion, loitering, or fishing patterns from skewing global thresholds. Detected outliers trigger state transitions that feed directly into segmenting vessel routes by behavior, where deterministic state machines classify transit, loitering, or anomalous drift.

The scoring engine applies a rolling window of 12–24 hours (configurable per vessel class) to establish local kinematic baselines. For a windowed speed series $v_t$ with rolling median $\tilde{v}$, a ping is flagged when its absolute deviation exceeds the scaled Median Absolute Deviation by the configured multiplier:

$$\left| v_t - \tilde{v} \right| > k \cdot c \cdot \operatorname{median}\!\left( \left| v_i - \tilde{v} \right| \right), \quad k = 3.5,\; c = 1.4826$$

The constant $c = 1 / \Phi^{-1}(0.75) \approx 1.4826$ makes the scaled MAD a consistent estimator of the standard deviation $\sigma$ under a Gaussian, so the threshold $k$ stays interpretable as “robust sigmas.” Unlike a mean/variance baseline, the median and MAD are not dragged upward by the very spikes they must detect, which keeps precision high across heterogeneous traffic densities and eliminates the need for jurisdiction-specific threshold tuning. The kinematic_confidence column provides a continuous score between 0.0 and 1.0 for downstream ranking and prioritization.

The contrast below is the core reason this stage uses MAD rather than a mean/standard-deviation baseline: a single erroneous fix inflates the mean band wide enough to swallow the genuine spike, whereas the median band — anchored on robust statistics — holds steady and isolates it cleanly.

Signal Gap Handling and Deterministic Recovery #

AIS dropouts occur in coastal shadow zones, during high-traffic interference, or when transponders enter standby mode. Blind linear interpolation across gaps exceeding 15 minutes introduces false kinematic spikes that corrupt downstream scoring. The pipeline implements gap-aware masking rather than imputation: any interval where dt_sec crosses the 15-minute threshold is flagged, the spanning derivative is discarded, and the post-gap rolling baseline is allowed to re-establish from real observations only.

Gap boundaries are explicitly tagged via the gap_mask boolean column in the output schema. Derivative calculations are suspended across masked intervals, and rolling statistics are reset to prevent stale baselines from contaminating post-gap scoring. This deterministic isolation guarantees that missing data never masquerades as anomalous behavior. Operators can audit gap distribution per MMSI by grouping on gap_mask == True before any behavioral classification step — a discipline also applied in the real-time Kafka consumer workflow where message deduplication and timestamp drift are the upstream sources of artificial gaps.

Validation Gates and Quality Control #

Three mandatory checkpoints must pass before anomaly-flagged output is promoted to the downstream segmentation stage:

Kinematic sanity bounds. Speed values must not exceed 50 knots for surface vessels (adjust for hydrofoils or fast patrol craft). Heading values must fall within [0.0, 360.0]. Any partition where more than 5% of pings breach these bounds is quarantined for manual review.

def validate_kinematics(df: pl.DataFrame) -> None:
    """Raise ValueError if kinematic bounds are violated at population level."""
    speed_violations = df.filter(pl.col("speed_kn") > 50.0).height
    total_valid = df.filter(pl.col("speed_kn").is_not_null()).height
    if total_valid > 0 and speed_violations / total_valid > 0.05:
        raise ValueError(
            f"Kinematic bounds violated: {speed_violations}/{total_valid} pings exceed 50 kn. "
            "Partition quarantined for forensic review."
        )
    logger.info("Kinematic validation passed: %d violations in %d pings.", speed_violations, total_valid)

Anomaly rate drift monitoring. When the fleet-level ANOMALOUS rate exceeds 15% across a partition set, the pipeline emits a warning and halts promotion. A sudden spike typically indicates corrupted ingestion batches, clock drift in the AIS transponder, or a major maritime event requiring manual triage. When the root cause traces back to the message layer — deduplication misses or non-monotonic timestamps — diagnose it against the AIS Kafka consumer troubleshooting reference before re-running this stage.

def check_anomaly_rate(df: pl.DataFrame, ceiling: float = 0.15) -> None:
    """Halt pipeline promotion if fleet anomaly rate exceeds the configured ceiling."""
    scored = df.filter(pl.col("anomaly_flag").is_not_null())
    rate = scored.filter(pl.col("anomaly_flag") == "ANOMALOUS").height / max(scored.height, 1)
    if rate > ceiling:
        raise RuntimeError(
            f"Fleet anomaly rate {rate:.1%} exceeds ceiling {ceiling:.1%}. "
            "Revert to last validated checkpoint and inspect quarantined partitions."
        )
    logger.info("Anomaly rate check passed: %.1f%%.", rate * 100)

Schema completeness. Output Parquet must contain mmsi, timestamp, lat, lon, speed_kn, heading_deg, dt_sec, gap_mask, anomaly_flag, and kinematic_confidence columns. Missing columns trigger a raised exception — not a silent default — to preserve lineage auditability across regulatory audits.

Common Failure Modes and Diagnosis #

Projection-induced speed spikes. If planar distance functions (e.g., Shapely distance() in a projected CRS) are used instead of pyproj.Geod.inv(), distances near high-latitude boundaries and antimeridian crossings can be orders of magnitude wrong, generating false ANOMALOUS flags for vessels in Arctic shipping lanes. Symptom: speed_kn values above 200 kn on otherwise normal tracks. Fix: enforce EPSG:4326 throughout and use geodesic distance exclusively.

Rolling window anchoring failures. Polars rolling_median(window_size="24h", by="timestamp") requires the timestamp column to be a true Datetime type with UTC timezone. If timestamps are stored as integer microseconds or naive datetimes, the window anchor silently falls back to row-count-based windows, breaking the time-keyed baseline. Symptom: anomaly rates collapse to near-zero or spike to 100% across partition boundaries. Fix: cast timestamp to pl.Datetime("us", "UTC") at ingestion, before any rolling computation.

Gap mask reset failures. If gap_mask is not propagated correctly through group_by + map_groups, stale rolling statistics from pre-gap pings bleed across dropout boundaries. Symptom: the first 10–20 pings after a long AIS dropout are systematically flagged ANOMALOUS even when kinematic values are normal. Fix: verify that _compute_kinematics_chunk receives sorted, per-MMSI partitions and that dt_sec > GAP_THRESHOLD_SECONDS correctly identifies the first post-gap ping.

Memory exhaustion from eager collection. Using lf.collect() on a multi-terabyte Parquet corpus without streaming triggers OOM termination. Symptom: worker processes killed at the OS level with no Python traceback. Fix: use lf.sink_parquet(output_path, streaming=True) for final materialisation, or collect in bounded chunks using fetch(n_rows) for development validation.

Pipeline Integration and Downstream Handoff #

The anomaly-flagged output Parquet partitions are written with the following column contract:

These partitions are consumed directly by the behavioral segmentation pipeline, which reads anomaly_flag and kinematic_confidence as inputs to its finite state machine. The gap_mask column prevents gap-boundary pings from triggering spurious state transitions between TRANSIT and LOITERING classifications. Streaming execution is enforced via sink_parquet(streaming=True) to guarantee bounded memory consumption throughout. All derivative windows are anchored to UTC timestamps to prevent timezone-induced scoring artifacts — a requirement that aligns with the schema standards enforced in the AIS NMEA sentence decoding workflow.

Related #

• Real-Time AIS Stream Ingestion Pipelines — upstream ingestion stage that produces the partitioned Parquet this pipeline consumes
• Speed and Heading Profiling for Maritime Analytics — companion kinematic workflow covering SOG/COG quality filtering
• Segmenting Vessel Routes by Behavior — downstream stage that classifies transit, loitering, and anomalous drift using this pipeline’s output
• Building an AIS Kafka Consumer in Python — real-time ingestion troubleshooting reference for gap and timestamp drift issues
• AIS Vessel Tracking & Route Automation — parent section