Removing Multibeam Refraction Artifacts #

A multibeam echosounder measures range along each beam and bends that ray through the water column using a sound velocity profile (SVP). When the SVP is wrong, the beams bend by the wrong amount and the seafloor is mis-placed — most visibly at the swath edges, where the outer beams travel at steep grazing angles and accumulate the largest error. The result is a swath that curls up at the edges (a “smile”) or droops at the edges (a “frown”), even over perfectly flat seafloor. This page sits inside Removing Bathymetric Artifacts and Noise and treats one specific artifact: refraction curvature from an incorrect SVP. It is distinct from the isolated vertical outliers handled by automated spike removal in sonar datasets — refraction is a systematic, angle-dependent bias, not a scattering of random spikes, and a spike filter will not touch it.

Root Cause: Snell’s-Law Ray Bending With a Wrong SVP #

Sound does not travel in a straight line through the ocean. As it crosses layers of differing temperature, salinity, and pressure, its speed changes and the ray bends according to Snell’s law. Ray tracing enforces the invariant that the ray parameter $p = \sin\theta / c$ is constant along the path, so within a layer of sound speed $c_i$ the beam makes an angle

$$\sin\theta_i = p\, c_i, \qquad p = \frac{\sin\theta_0}{c_0}$$

where $\theta_0$ and $c_0$ are the launch angle and surface sound speed at the transducer. The depth of a beam is the accumulated vertical travel through the layers,

$$z = \sum_i \Delta t_i \, c_i \cos\theta_i$$

and its across-track offset accumulates the horizontal component. If the SVP used in this sum is wrong by $\delta c(z)$, the reconstructed depth is biased. Crucially the bias depends on beam angle: a nadir beam ($\theta_0 \approx 0$) travels nearly straight down and is almost insensitive to $c$, while an outer beam at $\theta_0 = 60^\circ$–$70^\circ$ has a long slanted path where the velocity error integrates into a large vertical displacement. To first order the across-track depth bias grows steeply with the launch angle,

$$\Delta z(\theta_0) \;\approx\; z\,\tan\theta_0 \cdot \frac{\overline{\delta c}}{\bar c}$$

so it is near zero at nadir and maximal at the outer beams, and it is symmetric about nadir when the error is a uniform velocity bias. That symmetric, angle-growing signature is the fingerprint that separates refraction from real terrain.

You can reproduce the signature synthetically. The snippet applies a uniform velocity bias to a flat bottom and plots the across-track residual — the symmetric, angle-growing curve is unmistakable.

import numpy as np

# Flat seafloor at 30 m; symmetric beam angles across the swath.
angles_deg = np.linspace(-65.0, 65.0, 131)
true_depth = np.full_like(angles_deg, 30.0)

# A +2% surface-sound-speed error biases outer beams the most.
dc_over_c = 0.02
biased = true_depth * (1.0 + dc_over_c * np.tan(np.radians(np.abs(angles_deg))))
residual = biased - true_depth

print("nadir residual (m):", round(float(residual[65]), 3))     # ~0.0
print("outer-beam residual (m):", round(float(residual[0]), 3)) # large, symmetric

The nadir residual is essentially zero while the outer beams carry the full error — the defining shape a correction must target.

Step-by-Step Fix With Production Code #

Step 1 — Detect via across-track residual versus beam angle #

De-trend each swath against a smooth reference (a local plane or the nadir-referenced depth) and bin the residuals against beam angle. A refraction artifact shows up as a residual that is small at nadir, grows monotonically toward both edges, and is symmetric. Real slope, by contrast, is asymmetric about nadir.

import logging
from dataclasses import dataclass

import numpy as np
import pandas as pd

logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")
log = logging.getLogger("refraction")


@dataclass(frozen=True)
class RefractionSignature:
    symmetric_gain: float   # residual growth per unit tan(angle)
    asymmetry: float        # port-vs-starboard imbalance; ~0 for pure refraction
    is_refraction: bool


def detect_refraction(swath: pd.DataFrame, sym_tol: float = 0.15) -> RefractionSignature:
    """Test a swath for a symmetric, angle-growing refraction residual.

    Args:
        swath: columns 'beam_angle_deg', 'depth', 'nadir_depth' for one ping set.
        sym_tol: max allowed |port-starboard| gain ratio to call it refraction.

    Raises:
        ValueError: if required columns are missing or the swath is empty.
    """
    required = {"beam_angle_deg", "depth", "nadir_depth"}
    if not required.issubset(swath.columns) or swath.empty:
        raise ValueError(f"swath needs columns {required} and must be non-empty")

    resid = (swath["depth"] - swath["nadir_depth"]).to_numpy()
    t = np.tan(np.radians(swath["beam_angle_deg"].to_numpy()))
    port = t < 0
    stbd = t > 0
    # Slope of residual vs |tan(angle)| on each side.
    gain_p = np.polyfit(np.abs(t[port]), resid[port], 1)[0] if port.sum() > 2 else 0.0
    gain_s = np.polyfit(np.abs(t[stbd]), resid[stbd], 1)[0] if stbd.sum() > 2 else 0.0
    sym_gain = 0.5 * (gain_p + gain_s)
    denom = abs(sym_gain) or 1e-9
    asymmetry = abs(gain_p - gain_s) / denom
    is_refraction = abs(sym_gain) > 1e-3 and asymmetry < sym_tol
    log.info("sym_gain=%.4f asymmetry=%.3f refraction=%s", sym_gain, asymmetry, is_refraction)
    return RefractionSignature(float(sym_gain), float(asymmetry), bool(is_refraction))

Step 2 — Estimate a per-swath correction #

Given a symmetric signature, estimate a correction. The rigorous route re-ray-traces with an adjusted SVP; a fast, defensible field route fits the empirical residual-versus-angle curve and subtracts it. The empirical model below fits the outer-beam de-trend as a function of $\tan\theta$, which matches the first-order physics above.

def estimate_detrend(swath: pd.DataFrame, degree: int = 2) -> np.ndarray:
    """Fit an empirical residual-versus-tan(angle) de-trend for one swath.

    Returns polynomial coefficients (highest degree first) mapping
    tan(beam_angle) -> depth correction to subtract.

    Raises:
        ValueError: if too few beams remain to fit the requested degree.
    """
    resid = (swath["depth"] - swath["nadir_depth"]).to_numpy()
    t = np.tan(np.radians(swath["beam_angle_deg"].to_numpy()))
    if t.size <= degree + 1:
        raise ValueError(f"need > {degree + 1} beams to fit degree {degree}")
    coeffs = np.polyfit(t, resid, degree)
    log.info("de-trend coeffs=%s", np.round(coeffs, 5).tolist())
    return coeffs

Step 3 — Apply the correction per ping and flag the uncorrectable #

Apply the de-trend per ping across the swath, and — critically — flag beams whose residual exceeds the model’s valid envelope rather than force-flattening them. A beam that is 3 m off when the model explains only 0.4 m is not refraction; flattening it would fabricate seafloor.

def apply_refraction_correction(
    swath: pd.DataFrame,
    coeffs: np.ndarray,
    max_residual_m: float = 0.5,
) -> pd.DataFrame:
    """Subtract the fitted de-trend per ping and flag uncorrectable beams.

    Adds columns 'depth_corr' and 'refraction_flag'. Beams whose post-
    correction residual still exceeds max_residual_m are flagged True and
    must be excluded downstream, not trusted as flattened.

    Raises:
        ValueError: if the correction makes agreement worse on average, which
            signals a mis-fit (e.g. real slope mistaken for refraction).
    """
    out = swath.copy()
    t = np.tan(np.radians(out["beam_angle_deg"].to_numpy()))
    correction = np.polyval(coeffs, t)
    out["depth_corr"] = out["depth"].to_numpy() - correction

    before = np.abs(out["depth"].to_numpy() - out["nadir_depth"].to_numpy())
    after = np.abs(out["depth_corr"].to_numpy() - out["nadir_depth"].to_numpy())
    if after.mean() > before.mean():
        raise ValueError("correction increased mean residual; refit — likely real slope")

    out["refraction_flag"] = after > max_residual_m
    log.info("flagged %d/%d beams beyond correctable envelope",
             int(out["refraction_flag"].sum()), len(out))
    return out

Verification and Acceptance Test #

The acceptance criterion is physical: over genuinely flat overlapping terrain, corrected outer-beam depths must agree with nadir depths within a survey-order tolerance. Assert it and refuse to pass a swath that still curls.

def verify_outer_beam_agreement(
    corrected: pd.DataFrame,
    tol_m: float = 0.25,
    outer_angle_deg: float = 45.0,
) -> None:
    """Assert corrected outer beams agree with nadir within tol_m.

    Raises:
        AssertionError: if the RMS outer-vs-nadir disagreement exceeds tol_m,
            meaning residual curvature remains and the swath is not clean.
    """
    outer = corrected[np.abs(corrected["beam_angle_deg"]) >= outer_angle_deg]
    if outer.empty:
        raise AssertionError("no outer beams to verify")
    valid = outer[~outer["refraction_flag"]]
    diff = valid["depth_corr"].to_numpy() - valid["nadir_depth"].to_numpy()
    rms = float(np.sqrt(np.mean(diff ** 2))) if diff.size else float("inf")
    if rms > tol_m:
        raise AssertionError(f"outer-beam RMS {rms:.3f} m exceeds tol {tol_m} m")
    log.info("outer-beam agreement OK: RMS=%.3f m", rms)

A passing swath logs a sub-tolerance RMS; a failing one raises before the data reaches a grid, so a residual smile cannot silently propagate into slope and volume products.

Edge Cases and Gotchas #

• Real slope masquerading as refraction — and vice versa. A uniformly sloping bottom produces an across-track trend too, but an asymmetric one: deeper on the down-slope side, shallower on the up-slope side. Refraction is symmetric about nadir; slope is not. The asymmetry metric in Step 1 is the discriminator — never de-trend a swath whose residual is one-sided, or you will flatten genuine bathymetry. Cross-check against adjacent overlapping lines: real slope is consistent line-to-line, refraction curls independently per line as the SVP drifts.
• Tide confused with refraction. A constant vertical offset between overlapping lines is a tide/datum problem, not refraction — refraction is angle-dependent and near-zero at nadir, whereas a tide error shifts the whole swath uniformly including nadir. If nadir depths themselves disagree line-to-line, fix the vertical reduction first; see Fixing Tide-Induced Offsets Between Survey Lines. Applying a refraction de-trend to a tide error will not remove the offset and can hide it.
• Thermocline variability within a survey. A single cast SVP goes stale as the boat transits into water with a different thermal structure, so refraction can appear and disappear along a line. Fit the correction per swath (or per short time window), not once per survey, and prefer more frequent casts or a moving-vessel profiler over a single morning cast. Beams flagged uncorrectable in Step 3 should feed back into the survey plan as a signal that a fresh SVP is overdue.

Related #

• Automated Spike Removal in Sonar Datasets — remove isolated vertical outliers, distinct from systematic refraction curvature
• Fixing Tide-Induced Offsets Between Survey Lines — the constant-offset artifact refraction is often confused with
• Point Cloud Filtering for Multibeam Sonar — upstream cleaning that normalizes beams before refraction analysis

Up: Removing Bathymetric Artifacts and Noise