Coastal & Marine Spatial Analysis Pipelines

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.

Seabed-Tracking Gate Over the Per-Cell Depth Distribution Horizontal axis is return count, vertical axis is depth with shallow at the top and the seabed at the bottom. A tall dense bar cluster marks the modal seafloor near the bottom. Sparse short bars above it represent water-column noise from fish, bubbles, thermoclines and multipath. A dashed gate line sits a tolerance above the modal seafloor: everything above the gate is labelled rejected water-column noise, everything at or below is kept, including the dense bed and a small real feature bump that rises but stays under the gate. depth (shallow ↑) return count per depth bin fish / bubbles / thermocline / multipath seabed gate = mode + tolerance REJECT above gate KEEP at / below gate real feature (kept) modal seafloor

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.

Why Min-Z Picks the Fish and Not the Seabed A histogram of every sounding falling inside one grid cell, plotted against depth. A tall narrow cluster near minus 24 metres is the seabed return. A sparse scatter of returns between minus 6 and minus 18 metres comes from a fish school, a thermocline reflection and bubble wash under the hull. Three candidate estimators are marked: minimum depth, which lands on the shallowest bubble return at minus 6.2 metres and reports a 6-metre shoal; the mean, which is pulled up to minus 21.4 metres by the tail; and the modal or low-percentile estimator, which lands inside the dense cluster at minus 23.9 metres and recovers the true seabed. A caption notes that a navigation-surface convention of "shoalest wins" is exactly what makes the naive minimum dangerous here. depth of returns within one grid cell (m) count bubbles · fish · thermocline seabed min-Z → −6.2 m a shoal that is not there mean → −21.4 m dragged up by the tail mode / p10 of the cluster → −23.9 m, the seabed −6 −14 −22 −28 “Shoalest wins” is the right rule for charting and the wrong rule for gridding raw returns. Apply the shoal-biased rule after the water column is removed, never as the thing that removes it.

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.

Three Kinds of Water-Column Return, Three Different Discriminators Three panels, one per noise source. Panel one, fish school: two overlapping survey lines cross the same patch of seabed at different times, and the shallow cluster appears at different depths and positions in each, so a cross-line consistency test removes it. Panel two, thermocline: a flat layer at a consistent depth spans the whole survey and shows up as a distinct second mode in a survey-wide depth histogram, so a global layer detection removes it. Panel three, bubble wash: returns hug the nadir beam directly beneath the vessel in every ping regardless of position, so a beam-angle and range-gate rule removes it. A footer states that a single statistical outlier filter catches parts of all three and all of none. fish school seabed line 1 (●) and line 2 (○) disagree → cross-line consistency thermocline same depth, survey-wide a second mode in the depth histogram → global layer detection bubble wash always at nadir, every ping → beam-angle + range gate A single statistical outlier filter catches part of each and the whole of none — and the part it catches is the part that happened to be sparse, which on a thermocline is nothing at all. Run the three discriminators in this order: nadir gate, layer detection, cross-line consistency. Each one makes the next one's population cleaner.

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.

Up: Point Cloud Filtering for Multibeam Sonar