Finding the Breakpoint in Travel Time Data

Plot arrival time against distance for seismic waves and you’ll see a distinctive kink: one slope for surface waves traveling through topsoil, another steeper slope where the wave refracts through bedrock. The breakpoint tells you the depth to the layer boundary. But if you have 30 geophones and noisy data, eyeballing it isn’t enough.

Segmented regression finds that breakpoint algorithmically. You’re fitting two lines with a shared join point, minimizing the squared residuals across both segments. The challenge is that the breakpoint location is itself a parameter you’re estimating—it’s nonlinear in a way that ordinary least squares can’t handle directly.

R’s segmented package iterates: guess a breakpoint, fit two lines, measure the gradient at the join, adjust the breakpoint, repeat until convergence. For refraction data with a clear velocity contrast, it locks on in 3-4 iterations.

library(segmented)
dist <- c(5,10,15,20,25,30,35,40,45,50)
time <- c(0.013,0.026,0.039,0.052,0.061,0.068,0.075,0.082,0.089,0.096)
fit <- lm(time ~ dist)
seg <- segmented(fit, seg.Z = ~dist, psi = 25)
summary(seg)$psi  # Breakpoint estimate
slope(seg)         # Velocities in each layer

In 1985, you’d code the iteration yourself. Pascal’s strong typing and explicit loop control made the convergence logic transparent—you could watch the breakpoint walk toward the minimum as residuals shrank.

program SegmentedFit;
var x, y: array[1..10] of real;
    psi, sse, delta: real;
    i, iter: integer;
begin
  psi := 25.0;  { initial breakpoint guess }
  for iter := 1 to 20 do begin
    sse := FitTwoLines(x, y, psi);  { sum squared errors }
    delta := Gradient(x, y, psi);   { derivative at join }
    psi := psi - 0.5 * delta;
    if abs(delta) < 0.01 then break;
  end;
  writeln('Breakpoint: ', psi:6:2);
end.

The breakpoint moved 2.3 metres between my first and third geophone line this morning. Same algorithm, different answer—the clay lens doesn’t run straight.