Fitting GPR Hyperbolas with SciPy to Get Depth and Soil Velocity

A buried cable or pipe shows up in a ground penetrating radar (GPR) scan as a hyperbola. If you pick a handful of points along that curve, you can fit a model and recover three useful numbers at once: where the target is, how deep it is, and how fast the radar wave travels in that soil.

The model

For a point target at horizontal position x0 and depth d, the two-way travel time recorded at antenna position x is:

t(x) = 2 * sqrt((x - x0)^2 + d^2) / v

where v is the wave velocity in the ground (metres per nanosecond). Because depth and velocity both shape the curve, fitting the whole hyperbola separates them, which a single travel time cannot do.

The code

import numpy as np
from scipy.optimize import curve_fit

def hyperbola(x, x0, depth, v):
    """Two-way travel time (ns) for a point target at (x0, depth).
    x and depth in metres, v in m/ns."""
    return 2.0 * np.sqrt((x - x0) ** 2 + depth ** 2) / v

# Synthetic picks: a cable at x0=3.2 m, depth=1.1 m, in soil with v=0.10 m/ns
rng = np.random.default_rng(7)
x = np.linspace(2.0, 4.4, 25)
t = hyperbola(x, 3.2, 1.1, 0.10) + rng.normal(0, 0.4, x.size)  # +/-0.4 ns picking noise

p0 = [x[np.argmin(t)], 1.0, 0.1]                      # apex guess, depth, velocity
bounds = ([x.min(), 0.05, 0.03], [x.max(), 5.0, 0.30])  # 0.03-0.30 m/ns is physical
(x0, depth, v), cov = curve_fit(hyperbola, x, t, p0=p0, bounds=bounds)
err = np.sqrt(np.diag(cov))

eps_r = (0.2998 / v) ** 2   # relative permittivity from velocity (c = 0.2998 m/ns)
print(f"apex position : {x0:.2f} m  (+/- {err[0]:.2f})")
print(f"depth         : {depth:.2f} m  (+/- {err[1]:.2f})")
print(f"velocity      : {v:.3f} m/ns (+/- {err[2]:.3f})")
print(f"permittivity  : {eps_r:.1f}")

Output

apex position : 3.21 m  (+/- 0.01)
depth         : 1.10 m  (+/- 0.01)
velocity      : 0.101 m/ns (+/- 0.001)
permittivity  : 8.9

Why the bounds matter

Radar velocity in real ground sits roughly between 0.03 m/ns (wet clay) and 0.30 m/ns (air). Bounding the fit to that range stops the optimiser from wandering into physically impossible solutions when picks are noisy or the hyperbola is clipped.

The fitted velocity also gives the relative permittivity, eps_r = (c / v)^2, which you can reuse to convert travel times to depth for other targets in the same scan.

In practice

On real surveys, picks come from processed radargrams and are repeated across many targets along a route, which is how depth profiles such as an L-section GPR survey are built. If you want to see how those surveys are scoped and priced, the GPR survey price list is a good reference.