Removing Water-Column Noise from Multibeam Data #

Point-cloud cleaning for a multibeam survey is not only about spikes on the seabed — the harder problem sits in the water column above it. Fish schools, entrained bubbles under the hull, thermocline layering, vessel wake, and multipath reflections all return acoustic energy that the sounder logs as valid soundings floating above the true bed. Feed those directly into a gridder and the damage depends on the rule: a minimum-depth (shoal-biased) navigation surface will lock onto the shallowest false return and report a hazard that is not there, while a mean surface will pull the whole cell upward. This page sits inside the multibeam point-cloud filtering workflow and shows how to strip mid-water returns with a modal-seafloor gate, density-based 3D filtering, and a keep-below-surface rule — while deliberately protecting real shoal features that a naive filter would delete.

Why Naive Min-Z Gridding Is Dangerous #

The physics is unavoidable: sound scatters off anything with an acoustic impedance contrast, so a herring swim-bladder, a bubble curtain, or the sharp density step of a thermocline all generate returns long before the ping reaches the seabed. In a well-behaved cell the depth distribution is a tight cluster at the bed; with water-column noise it grows a long shallow tail. The shoal-bias principle that keeps hydrography conservative — when uncertain, report the shallower depth — becomes a liability here, because the shallowest point in the cell is now a fish, not a rock. The reproduction below builds one cell’s depth distribution and shows the mid-water cloud sitting well above the seabed mode:

import numpy as np

rng = np.random.default_rng(7)
# True seabed near -42 m: a tight, dense cluster of returns.
seabed = rng.normal(-42.0, 0.15, size=900)
# Water-column noise: a sparse cloud 5-25 m above the bed.
noise = rng.uniform(-37.0, -18.0, size=70)
cell_z = np.concatenate([seabed, noise])

hist, edges = np.histogram(cell_z, bins=24)
modal_depth = edges[hist.argmax()]          # near -42 m, the real bed
naive_min = cell_z.max()                     # a fish return, ~ -18 m
print(round(float(naive_min - modal_depth), 1))  # tens of metres too shoal

A minimum-depth grid of that cell would publish a depth roughly twenty metres shallower than the real seabed — a fabricated shoal directly in the navigable water column. The task is to remove the tail before any shoal-bias rule runs, without touching the dense seabed cluster or any genuine feature that rises from it.

Step-by-Step Fix #

1. Estimate the modal seafloor and apply a seabed-tracking gate #

The dense mode of the depth distribution in each cell is a robust estimate of the local seabed. Reject returns that sit more than a tolerance above that mode; keep everything at or below it, so real rock never gets clipped:

import logging
import numpy as np

logger = logging.getLogger(__name__)


def seabed_gate(depths: np.ndarray,
                above_tolerance_m: float,
                bins: int = 32) -> np.ndarray:
    """Boolean keep-mask rejecting returns far above the modal seafloor.

    Args:
        depths: 1-D array of Z values for one cell (negative = down).
        above_tolerance_m: metres above the modal bed still kept as valid.
        bins: histogram bins used to locate the seabed mode.

    Returns:
        Boolean mask, True for points to keep.

    Raises:
        ValueError: if the cell is empty or the tolerance is negative.
    """
    if depths.size == 0:
        raise ValueError("empty cell passed to seabed_gate")
    if above_tolerance_m < 0:
        raise ValueError("above_tolerance_m must be >= 0")

    hist, edges = np.histogram(depths, bins=bins)
    modal_depth = 0.5 * (edges[hist.argmax()] + edges[hist.argmax() + 1])
    keep = depths <= modal_depth + above_tolerance_m
    logger.debug("cell mode=%.2f m, kept %d/%d returns",
                 modal_depth, int(keep.sum()), depths.size)
    return keep

The tolerance is the one safety-critical knob: set it from the expected feature relief plus the sounder’s vertical uncertainty, never smaller than the tallest real feature you must preserve, so a wreck rising ten metres proud of the bed survives while a fish twenty metres up does not. This is the same conservative posture as automated spike removal in sonar datasets, applied to the water column instead of the bed.

2. Refine with density-based 3D filtering #

A per-cell gate is fast but 1-D. Water-column noise is physically sparse in 3D, whereas the seabed and any real feature are dense and spatially continuous. A statistical outlier removal (SOR) pass, or a DBSCAN clustering that keeps only the dominant cluster, catches diffuse returns the vertical gate lets slip and respects continuity across cell boundaries:

import logging
import numpy as np
from sklearn.cluster import DBSCAN

logger = logging.getLogger(__name__)


def density_filter(xyz: np.ndarray,
                   eps_m: float,
                   min_samples: int) -> np.ndarray:
    """Keep only points belonging to the dominant spatial cluster.

    Uses DBSCAN in 3D: dense seabed and feature returns form one large
    cluster; sparse mid-water returns fall into noise (label -1) or small
    satellite clusters that are dropped.

    Args:
        xyz: (N, 3) array of easting, northing, depth in metres.
        eps_m: neighbourhood radius; ~ 1.5-3x the point spacing.
        min_samples: minimum neighbours for a core point.

    Returns:
        Boolean keep-mask of length N.

    Raises:
        ValueError: if xyz is not (N, 3) or eps_m is non-positive.
    """
    if xyz.ndim != 2 or xyz.shape[1] != 3:
        raise ValueError("xyz must be an (N, 3) array")
    if eps_m <= 0:
        raise ValueError("eps_m must be > 0")

    labels = DBSCAN(eps=eps_m, min_samples=min_samples,
                    n_jobs=-1).fit_predict(xyz)
    if (labels >= 0).sum() == 0:
        raise ValueError("DBSCAN found no dense cluster; check eps_m")

    counts = np.bincount(labels[labels >= 0])
    dominant = counts.argmax()
    keep = labels == dominant
    logger.info("density filter kept %d/%d points (%.1f%%)",
                int(keep.sum()), labels.size, 100.0 * keep.mean())
    return keep

For survey-scale clouds that will not fit in memory, run the same logic through a windowed PDAL pipeline rather than a single in-memory array — the streaming, chunked approach and its stages are covered in using PDAL for bathymetric point-cloud cleaning. PDAL’s filters.outlier (SOR) and filters.elm handle the mid-water case natively at scale.

3. Keep-below-surface and re-grid #

Apply both masks in order — vertical gate first to remove the obvious tail cheaply, density filter second to clean the diffuse remainder — then grid the survivors with the shoal-bias rule now that it can be trusted. Because every rejected point is recorded, the decision is auditable rather than a silent drop.

Verification and Acceptance Test #

Gate on the physical expectation: after filtering, the vertical spread of returns below the modal seabed in each populated cell should be within the sounder’s uncertainty, and the retained-point ratio should sit in a plausible band (a filter that keeps 20% has eaten the seabed; one that keeps 99.9% did nothing):

import logging
import numpy as np

logger = logging.getLogger(__name__)


def verify_filtered(depths_before: np.ndarray,
                    keep_mask: np.ndarray,
                    max_spread_m: float) -> dict:
    """Assert per-cell tightness and log the retained-point ratio."""
    kept = depths_before[keep_mask]
    if kept.size == 0:
        raise AssertionError("filter removed every return in the cell")
    spread = float(np.percentile(kept, 95) - np.percentile(kept, 5))
    if spread > max_spread_m:
        raise AssertionError(
            f"kept spread {spread:.2f} m exceeds {max_spread_m} m tolerance")
    ratio = float(keep_mask.mean())
    logger.info("retained %.2f%% of returns, spread %.2f m",
                100.0 * ratio, spread)
    return {"retained_ratio": ratio, "kept_spread_m": spread}

Difference the filtered surface against the raw minimum-depth surface to visualise exactly which shoal returns were removed, and confirm every removed cell was water-column noise rather than a genuine hazard before publishing.

Edge Cases and Gotchas #

• Real shoal features versus noise. A rock pinnacle, a wreck, or a boulder field rises from the bed as a dense, spatially continuous cluster, whereas a fish school is dense but disconnected from the seabed. Keep the density filter’s min_samples and eps_m tuned so a continuous feature stays part of the dominant cluster; when a shoal is safety-critical, apply the IHO shoal-bias rule and preserve the shallowest validated sounding rather than trusting the automatic filter. A hydrographer should review flagged Class-B (charting) features before they are deleted.
• Dense fish schools. A large, tight school can itself form a DBSCAN cluster big enough to rival the seabed, defeating the “keep the largest cluster” rule. Anchor the decision to depth as well as density — the dominant cluster must also be the deepest coherent surface — and fall back to the modal-seafloor gate, which is immune to how dense the mid-water cloud is.
• Wreck and obstruction shoal-biasing. Over known wrecks and obstructions the convention inverts: you deliberately retain the shoalest reliable return to protect mariners. Maintain a feature exclusion layer so the water-column filter is suppressed inside charted obstruction polygons, and route those cells through the least-depth workflow in removing bathymetric artifacts and noise rather than the automatic gate.

Related #

• Using PDAL for Bathymetric Point-Cloud Cleaning — running SOR and ELM filters at survey scale through a streaming pipeline.
• Automated Spike Removal in Sonar Datasets — the seabed-side counterpart to water-column cleaning.
• Point Cloud Filtering for Multibeam Sonar — parent workflow: ingest, classify, filter, and grid.

Up: Point Cloud Filtering for Multibeam Sonar