RouteMatching

Estimating Clock Drift with Linear Regression

The timestamp synchronisation topic establishes that a device clock has two independent errors: a constant offset present from the first fix, and a drift that accumulates through the session. A single median offset removes the first and leaves the second, which is why a fleet that “corrected its timestamps” still has half a second of error by the end of a shift.

Both terms come out of one linear fit. The complications are that the fit must be robust to a handful of bad pairs, that a power cycle resets the offset and invalidates a fit spanning it, and that the residuals need checking — a clean-looking slope on a structured residual is a model that is wrong in a way the number does not show.

Compatibility and Configuration Requirements

Requirement Value Notes
Reference time GNSS time from the fix Sub-millisecond, arrives with the observation
Device time The device’s own clock, unmodified If the device already corrects itself, there is nothing to fit
Session detection Power-cycle events, or a step detector A fit across a reset is meaningless
Fit Theil-Sen or Huber Ordinary least squares is rotated by a few bad pairs
Minimum pairs ~200 per session Fewer and the slope is dominated by noise

The Fit

from __future__ import annotations

import numpy as np
from sklearn.linear_model import TheilSenRegressor


def fit_drift(t_ref_s: np.ndarray, t_dev_s: np.ndarray) -> dict:
    """Fit device_time - reference_time = offset + drift * elapsed.

    Both inputs are seconds since an arbitrary epoch. Returns the offset in
    seconds, the drift in parts per million, and the residual scatter.
    """
    if len(t_ref_s) < 200:
        raise ValueError(f"{len(t_ref_s)} pairs is too few for a stable slope")

    elapsed = t_ref_s - t_ref_s[0]
    error = t_dev_s - t_ref_s

    model = TheilSenRegressor(random_state=0).fit(elapsed.reshape(-1, 1), error)
    offset_s = float(model.intercept_)
    drift_ppm = float(model.coef_[0]) * 1e6

    resid = error - model.predict(elapsed.reshape(-1, 1))
    return {
        "offset_s": offset_s,
        "drift_ppm": drift_ppm,
        "resid_mad_s": float(np.median(np.abs(resid - np.median(resid))) * 1.4826),
        "n_pairs": int(len(elapsed)),
    }


def apply_drift(t_dev_s: np.ndarray, t0_ref_s: float, fit: dict) -> np.ndarray:
    """Undo the fitted offset and drift, returning corrected times."""
    elapsed = t_dev_s - t0_ref_s
    return t_dev_s - fit["offset_s"] - fit["drift_ppm"] * 1e-6 * elapsed

Theil-Sen rather than ordinary least squares because the input contains outliers by construction. A handful of fixes whose GNSS time was momentarily wrong, or whose device timestamp was written during a flush, will rotate a least-squares line noticeably — and a rotated line applies a wrong correction to every fix in the session, not just to the bad ones.

The residual scatter is returned because it is the diagnostic. A well-behaved device gives a residual MAD in the low tens of milliseconds; a value in the hundreds means either the reference is not what you think or the session contains a reset.

One session fits; two sessions do not Clock error against elapsed time. Within one powered session the error is a clean line with an intercept of 1.8 seconds and a slope of 11 parts per million. Fitted across a power cycle that reset the offset, the same points produce a slope four times larger that describes neither session. Device clock error against elapsed time +3 s +1.5 s 0 session 1 — 11 ppm power cycle session 2 — 9 ppm fitted across both — 42 ppm, describes neither A power cycle resets the offset but not the crystal, so the slope is preserved and the intercept is not. Fitting across the reset produces a slope that is really the step divided by the window length.

Detecting Session Boundaries

Where the device reports a power-cycle event, use it. Where it does not — which is most fleets — detect the step directly. A reset shows up as a jump in the error series with no corresponding gap in the data, and it is straightforward to find with a rolling median.

def session_breaks(t_ref_s: np.ndarray, t_dev_s: np.ndarray,
                   step_s: float = 2.0, window: int = 31) -> np.ndarray:
    """Indices where the clock error steps, indicating a reset."""
    error = t_dev_s - t_ref_s
    med = np.array([
        np.median(error[max(0, i - window) : i + 1]) for i in range(len(error))
    ])
    return np.flatnonzero(np.abs(np.diff(med)) > step_s) + 1

The threshold should sit above the drift the clock can accumulate within the window and below the smallest reset worth splitting on. At 11 parts per million a 31-sample window at 1 Hz accumulates under a millisecond of drift, so a two-second threshold is comfortably clear of the noise floor.

Splitting too eagerly is cheap — two shorter fits of the same crystal give nearly the same slope — while splitting too rarely is expensive, because one missed reset corrupts the whole session’s correction. When in doubt, split.

Resets are steps, not drift Rolling median of clock error across a twelve-hour shift. The error climbs gently within each powered session and steps sharply at two points where the device lost power. The steps are two and three seconds, far larger than the drift accumulated within the detection window. Rolling median of clock error, one 12-hour shift reset −2.1 s reset +3.0 s session 1 session 2 session 3 Each session drifts at roughly the same rate — the crystal is unchanged — and starts from a new offset. Fitting one line to all three would report a slope four times the true one and correct nothing properly.

Reading the Residuals

A slope and an intercept are two numbers, and two numbers cannot tell you whether the model fits. The residuals can, and three patterns are worth recognising.

Flat and small. The model is right. Residual MAD in the low tens of milliseconds is what a healthy device gives, and nothing more is needed.

Flat and large. The model is right and the measurement is noisy — usually a device that timestamps on queue flush rather than on observation. The correction still helps, but the residual is a floor on how good the timestamps can get, and it should be published alongside the fit.

Structured. A residual with a visible curve or a periodic component means the linear model is wrong. A curve usually means temperature: crystals drift faster when warm, so a vehicle that heats through the morning has a slope that changes. A periodic component usually means the reference is being sampled at a rate that beats against the device’s own.

def residual_shape(resid: np.ndarray) -> str:
    """Crude classifier: is what is left over structured?"""
    half = len(resid) // 2
    drift_between_halves = abs(np.median(resid[:half]) - np.median(resid[half:]))
    scatter = np.median(np.abs(resid - np.median(resid))) * 1.4826
    return "structured" if drift_between_halves > 2 * scatter else "flat"

Where the residual is structured and temperature is the suspected cause, a piecewise fit over shorter windows is usually enough — an hour at a time keeps the linear approximation good without needing a thermal model.

Applying and Recording the Correction

Apply the correction to produce a new column rather than overwriting the original. The device time is evidence about the device, and the corrected time is a derived value; conflating them makes the fit impossible to re-derive or to revise.

fixes = fixes.with_columns(
    ts_corrected_utc=apply_drift(fixes["ts_device_s"], t0, fit),
    clock_offset_s=pl.lit(fit["offset_s"]),
    clock_drift_ppm=pl.lit(fit["drift_ppm"]),
    clock_resid_mad_s=pl.lit(fit["resid_mad_s"]),
)

Storing the fit parameters on every row looks redundant and is not. It makes the correction auditable without joining to a separate table, it makes a bad session identifiable by filtering on the residual, and it lets a downstream consumer decide that a session with a 400-millisecond residual is not good enough for its purposes.

Aggregate the drift figures per device family and watch them. A crystal’s drift is a physical property and should be stable; a family whose median drift moves has had a hardware or firmware change, and that is worth knowing before it shows up as a matching problem.

Reading what is left over Three residual patterns after fitting offset and drift. A flat, tight residual means the model is right. A flat but wide residual means the timestamps themselves are noisy. A curved residual means the linear model is wrong, usually because the crystal warmed through the session. Residual after subtracting the fitted offset and drift flat and tight MAD 18 ms — model is right flat and wide MAD 210 ms — noisy timestamps, publish it as a floor curved model is wrong — usually the crystal warming through the shift Only the first case needs no further thought. The second is honest and bounded; the third needs shorter fit windows. An hour-long piecewise fit keeps the linear approximation good without needing a thermal model. Store the residual MAD per session so a consumer can reject a fit that is not good enough for its purposes.

Common Pitfalls Specific to This Technique

Fitting across a power cycle. Produces a slope that is really the reset step divided by the window length, and applies it to every fix in both sessions.

Using ingest time as the reference. Network latency and queue buffering are not clock error, and a fit against arrival time mostly measures the mobile network.

Correcting a device that already disciplines its own clock. The fit then finds a slope near zero with a large residual, and applying it adds noise. Check for a near-zero slope with a poor residual before assuming a correction is needed.

Overwriting the raw timestamp. The original device time is the only evidence from which the fit can be re-derived if the method changes.

Related

Up to section Timestamp Synchronization for Multi-Device GPS Logs Related How to align GPS timestamps across mixed OBD-II and mobile devices Related Interpolating GPS Gaps During Tunnel Signal Loss