Coastal & Marine Spatial Analysis Pipelines

Correcting Heave-Induced Ripple in Bathymetry

A gridded surface shows a regular corrugation running across the survey lines, wavelength a few tens of metres, amplitude ten or twenty centimetres, remarkably even. It looks like a field of sand waves and it is very often the vessel going up and down. Distinguishing the two matters because they call for opposite responses: a bedform is the product, and a motion artefact is a defect that should be corrected at its source rather than filtered out of the surface. This page is the discrimination and the fix, within the artifact removal workflow.

Root cause: heave is compensated, and compensation is never exact

The motion sensor measures heave, roll and pitch, and the acquisition software applies them to place each beam. Three things break that chain, and all three produce a periodic residual at the wave period:

Lever-arm error. The motion sensor and the transducer are at different places on the vessel. Roll and pitch turn that offset into a vertical displacement, and an error of a few centimetres in the measured offset becomes a periodic depth error scaled by the vessel’s rotation.

Latency. Motion samples and ping times are stamped by different clocks. A constant offset of even 20 ms between them applies yesterday’s heave to today’s ping, and at a 6-second swell period with 1 m amplitude that is roughly 2 cm of error — worse in shorter, steeper seas.

Heave filter settling. A heave sensor high-pass filters acceleration to recover displacement, and that filter has a settling time and a period below which it under-reports. In a long swell the reported heave is smaller than the real one, and the residual is a scaled copy of the wave.

import numpy as np
from scipy import signal

def along_track_spectrum(residual: np.ndarray, spacing_m: float, speed_kn: float):
    """Peak period, in seconds, of a residual sampled along the track."""
    if spacing_m <= 0 or speed_kn <= 0:
        raise ValueError("spacing and speed must be positive")
    f, pxx = signal.welch(residual - residual.mean(), fs=1.0 / spacing_m,
                          nperseg=min(1024, len(residual)))
    peak_cycles_per_m = f[np.argmax(pxx)]
    speed_ms = speed_kn * 0.514444
    return float(1.0 / (peak_cycles_per_m * speed_ms)) if peak_cycles_per_m else float("inf")

Convert the spatial wavelength to a temporal period using the vessel speed and compare it against the swell period logged during acquisition. A match within a second or so is the single strongest piece of evidence available, and it takes one function call.

Three Tests, and a Bedform Fails All Three Three diagnostic panels. Panel one compares an along-track spectrum of the depth residual against the swell spectrum recorded by the vessel during acquisition; both peak at about six seconds, which a sand-wave field would not do. Panel two shows a ripple crest drawn across the full width of the swath: for heave the crest is a straight line perpendicular to the track because the whole vessel moved together, while for a bedform the crest is oblique and continues across the neighbouring line. Panel three shows two survey lines run four hours apart over the same patch of seabed: the heave ripples are out of phase between them because the swell was at a different phase, while a bedform would appear identically in both. 1 · period matches the swell 6.1 s solid: depth residual dashed: logged swell 2 · in phase across the swath heave: crests square to the track bedform: oblique, and continues next door 3 · two lines, four hours apart out of phase — the seabed did not move a bedform would match exactly A real bedform fails all three: wrong period, oblique crests, identical between lines. Any one of the three is suggestive. All three together is a diagnosis, and each costs a few lines of code against data the survey already holds. Run them before touching the surface: the fix for a diagnosis is different from the fix for a suspicion.

Step-by-step fix with production code

1. Confirm the phase relationship across the swath

Heave moves the whole vessel, so it displaces every beam of a ping by the same amount. A bedform does not care where the beams are. Comparing the residual at nadir against the residual at the swath edge separates them cleanly:

import logging

import numpy as np
import pandas as pd

logger = logging.getLogger(__name__)


def swath_phase_coherence(df: pd.DataFrame, *, nadir_deg: float = 10.0) -> float:
    """Correlation between the nadir residual and the outer-beam residual, per ping.

    Near 1.0 means the whole swath moved together — a vessel-motion signature.
    Near 0 means the residual is a property of the seabed under each beam.
    """
    inner = (df["angle_deg"].abs() <= nadir_deg)
    outer = (df["angle_deg"].abs() >= 45.0)
    if not inner.any() or not outer.any():
        raise ValueError("need both nadir and outer-beam soundings to compare")

    a = df[inner].groupby("ping_id")["residual_m"].mean()
    b = df[outer].groupby("ping_id")["residual_m"].mean()
    joined = pd.concat([a, b], axis=1, join="inner").dropna()
    if len(joined) < 50:
        raise ValueError(f"only {len(joined)} pings have both — too few to correlate")
    r = float(joined.corr().iloc[0, 1])
    logger.info("nadir/outer residual correlation over %d pings: %.3f", len(joined), r)
    return r

A correlation above about 0.7 is a motion signature. Below about 0.3 the residual is seabed, and the next place to look is the sound-velocity profile rather than the motion sensor.

2. Fix it at the source: lever arm and latency

Once identified, the correction belongs in the reprocess, not in the grid. Both parameters can be solved for from the data itself, because both produce residuals with a distinctive dependence on vessel attitude:

def solve_latency_s(df: pd.DataFrame, candidates=np.arange(-0.100, 0.101, 0.005)) -> float:
    """Time offset between the motion series and the ping series, by minimising
    the variance of the corrected residual."""
    best, best_var = 0.0, float("inf")
    for lag in candidates:
        heave = np.interp(df["ping_time_s"].to_numpy() + lag,
                          df["motion_time_s"].to_numpy(), df["heave_m"].to_numpy())
        var = float(np.var(df["z_uncorrected"].to_numpy() + heave))
        if var < best_var:
            best, best_var = float(lag), var
    if abs(best) >= candidates.max() - 1e-9:
        raise ValueError(
            f"best latency {best:+.3f} s is at the edge of the search range — "
            "widen it, or the residual is not a timing problem")
    logger.info("solved latency %+.3f s (residual variance %.5f m²)", best, best_var)
    return best

The edge-of-range guard is the part worth keeping. An optimiser that returns its own boundary has not found a minimum; it has run out of room, and accepting the value produces a “correction” that makes the next survey worse.

3. Filtering is the last resort, and it is recorded

If the raw motion data is unavailable — a legacy dataset, a third-party deliverable — the ripple can be filtered out along-track with a notch at the known wavelength. This is a cosmetic repair and it must be labelled as one, because it also removes any real bedform at that wavelength, and bedforms at swell wavelength are exactly the ones that a swell-dominated shelf produces.

def notch_along_track(z: np.ndarray, spacing_m: float, wavelength_m: float,
                      q: float = 8.0) -> np.ndarray:
    """Suppress one along-track wavelength. Removes real bedforms at it too."""
    from scipy.signal import iirnotch, filtfilt
    w0 = (1.0 / wavelength_m) / (0.5 / spacing_m)
    if not 0 < w0 < 1:
        raise ValueError(f"wavelength {wavelength_m} m is outside the Nyquist range")
    b, a = iirnotch(w0, q)
    logger.warning("notching %.1f m along-track: this is cosmetic, not a correction — "
                   "record it in the product metadata", wavelength_m)
    return filtfilt(b, a, z, axis=0)

The logger.warning in a function that works correctly is deliberate. The risk here is not that the code fails, it is that a cosmetic repair is mistaken for a correction three months later, and a warning in the run log is the cheapest available defence.

Correct at the Source, or Remove the Bedform Too A profile containing a twenty-metre heave ripple superimposed on a twenty-eight-metre real sand-wave field. In the reprocessed version, the heave ripple is gone because the latency was corrected at the source, and the sand waves remain clearly visible. In the notch-filtered version, the notch centred on twenty metres has a finite width that also covers twenty-eight metres, so both the ripple and the sand waves are gone and the profile is smooth. A caption records that the second surface looks cleaner and contains less of the seabed, and that nothing in the file distinguishes the two. raw: 20 m ripple over 28 m sand waves reprocessed with the solved latency sand waves intact notch filtered at 20 m smoother, and the sand waves are gone The lower surface looks cleaner and holds less seabed. Nothing in the file says which one you have.

Verification and acceptance test 

import pytest

def test_ripple_is_gone_and_the_bedform_is_not(before, after, spacing_m, speed_kn):
    p_before = along_track_spectrum(before, spacing_m, speed_kn)
    p_after = along_track_spectrum(after, spacing_m, speed_kn)
    assert abs(p_before - 6.1) < 1.0, "the input did not contain the expected ripple"
    assert abs(p_after - 6.1) > 1.5, "the ripple period still dominates after correction"

    # Energy at the bedform wavelength must survive the correction.
    from scipy import signal
    def energy_at(z, lam):
        f, pxx = signal.welch(z - z.mean(), fs=1.0 / spacing_m, nperseg=1024)
        return float(np.interp(1.0 / lam, f, pxx))
    assert energy_at(after, 28.0) > 0.6 * energy_at(before, 28.0), (
        "the correction removed the sand-wave energy as well as the ripple")

The second assertion is what separates a correction from a smoothing. A filter that removes the ripple will fail it, which is the intended outcome.

What the ripple costs downstream

A twenty-centimetre corrugation is small against most depth budgets and large against almost every derivative computed from the surface, which is where the cost actually lands.

Slope and curvature. A derivative amplifies short-wavelength content in proportion to how short it is. A 20 cm ripple at 20 m wavelength contributes a slope of roughly 1.1°, which on a shelf whose real slope is 0.3° means the slope raster is dominated by the artefact. Every habitat model keyed on slope inherits it.

Volume differences. Two surveys with ripples at different phases differ by up to twice the amplitude wherever the phases oppose. Over a dredge area of a few hectares that is thousands of cubic metres of apparent change, and it appears as a plausible spatial pattern rather than as noise.

Feature detection. Any detector tuned for bedforms finds the ripple everywhere, at high confidence, because it is more regular than a real bedform field. A pipeline that counts sand waves will report a suspiciously uniform population.

The general shape of this is that the artefact is worst in exactly the products that motivated the survey. A raw depth grid tolerates it; a change-detection or habitat product does not, and those are usually the deliverable.

The Same 20 cm, Priced Against Four Products Four horizontal bars showing the size of a twenty-centimetre heave ripple relative to each product's tolerance. Against a depth grid with a fifty-centimetre budget the bar reaches forty per cent and is annotated as tolerable. Against a slope raster on a shelf of 0.3 degrees real slope, the ripple contributes 1.1 degrees and the bar runs off the scale, annotated as the artefact dominating the signal. Against a volume difference between two surveys, the opposing phases give up to forty centimetres of apparent change over the whole area, annotated as thousands of cubic metres. Against a bedform count, the ripple is more regular than any real bedform field and the detector finds it everywhere, annotated as a uniform false population. a 20 cm ripple, against four tolerances depth grid 40% of a 0.5 m budget — tolerable slope raster 1.1 deg of false slope on a 0.3 deg shelf — the artefact is the signal volume difference up to 40 cm where phases oppose — thousands of m³ bedform count more regular than any real field — found everywhere, at high confidence The artefact is worst in exactly the products the survey was commissioned for. A depth grid tolerates it; a change-detection or habitat product does not.

Edge cases and gotchas

  • Ripple that reverses sign with heading is almost always latency rather than lever arm, because a timing offset flips its effect when the vessel turns around and a geometric offset does not. Splitting the residual by heading before solving is the fastest way to tell them apart.
  • A “ripple” at exactly the ping interval is not heave at all; it is a per-ping bias, usually a bottom-detection threshold interacting with a changing seabed. It has no wave period and will not correlate with the swell log.
  • Real bedforms at swell wavelength exist and are common. On a swell-dominated shelf the wave field and the bedform field share a scale, because one built the other. That is precisely why the cross-line test matters: it is the only one of the three that cannot be confounded by a coincidence of scale.

Up: Removing Bathymetric Artifacts & Noise