Adding Realistic GPS Noise to Synthetic Vehicle Trajectories

When a synthetic vehicle path is perturbed with naive isotropic Gaussian noise, it fails downstream validation because real GNSS error is temporally correlated, anisotropic, and satellite-geometry dependent — this page gives the per-receiver error budgets and verification benchmarks that make injected noise indistinguishable from a real receiver’s output.

Part of Noise Injection & Stochastic Drift: the parent area defines the ordered transform — project, perturb, re-validate, inverse-project — and this page fills in the numbers that transform needs: which σ and θ to use for which receiver class, and how to prove the result is statistically faithful before it reaches a training set.

Root Cause: Why White Noise Is Detectable

GNSS positional error is a composite signal driven by atmospheric propagation delay, receiver clock bias, satellite ephemeris inaccuracy, and local multipath reflection. Treating these as independent white noise destroys the autocorrelation structure that downstream tracking filters — Kalman, particle, or spline-based — depend on. A model trained against white noise learns nothing about the failure modes that dominate in production: urban-canyon navigation, tunnel exits, and cold-start drift.

The error must be decomposed into three orthogonal components, each a distinct stochastic process with its own time constant:

Radial GPS error quantiles for a Rayleigh model and a realistic contaminated mixture The horizontal axis is the percentile of fixes, from the first to the ninety-ninth; the vertical axis is radial error in metres. The lower curve is a pure Rayleigh distribution — the radial error you get from independent Gaussian noise in each axis — fitted to a typical reported accuracy. The upper curve is a mixture in which ninety-four per cent of fixes come from that same Rayleigh and six per cent come from a component five and a half times wider, representing the fixes taken with poor satellite geometry, partial obstruction or during reacquisition. Through the first four fifths of the distribution the two curves are almost indistinguishable, which is why a model fitted on the median or the mean passes every check anybody runs. Above the eightieth percentile they separate sharply, and at the ninety-ninth the mixture is more than twice the Rayleigh. Markers give both values at the ninety-fifth and ninety-ninth percentiles. The note underneath states why this matters: the fixes a consumer complains about are exactly the ones in the part of the distribution the Gaussian model gets wrong, and the parameter that fixes it is a contamination fraction — one number, which belongs in the noise manifest alongside the sigma. The two models agree everywhere except where it matters p0 p25 p50 p75 p95 p99 0 10 20 30 40 50 percentile of fixes radial error (m) p95: 16 m vs 11 m p99: 48 m vs 14 m Rayleigh fitted to reported accuracy realistic mixture, 6% degraded component 20,000 fixes per model, σ = 4.5 m per axis, 6% contamination at 5.5σ, fixed seed. The fixes a consumer complains about are exactly the ones the Gaussian model gets wrong, and the parameter that fixes it is a contamination fraction — one number, which belongs in the noise manifest beside σ.
The two models agree through the first four fifths of the distribution and separate exactly where consumers complain.
  1. Slow-varying bias — ionospheric and tropospheric delay plus receiver clock drift. Low-frequency wander over minutes to hours, modeled as a random walk (integrated Wiener process).
  2. Correlated spatial jitter — multipath reflection and signal attenuation in urban canyons, strongly dependent on vehicle heading and surrounding geometry. Best represented by an Ornstein-Uhlenbeck process with environment-dependent reversion speed θ.
  3. High-frequency measurement noise — thermal noise in the receiver front-end and quantization error. Zero-mean Gaussian with variance scaled by satellite geometry (HDOP/VDOP).

Critically, noise must be applied in a local tangent plane (ENU or NED), never directly in latitude/longitude. Applying metric noise to geographic degrees introduces latitude-dependent anisotropy — the same silent-corruption class as coordinate reference system drift, where numbers stay finite and the map still renders so the bug survives review. The injection layer should also operate on the kinematic state vector — position, velocity, heading — so it never violates the constraints enforced by physics-based path generation upstream.

Three-component GNSS noise decomposition summed onto a clean baseline path Top row: a clean, kinematically valid baseline path. Second row adds slow bias drift modeled as a random walk scaled by sigma_b in metres per minute, wandering over minutes to hours. Third row adds correlated multipath jitter modeled as an Ornstein-Uhlenbeck process with reversion rate theta and a 10 to 60 second time constant. Fourth row adds a high-frequency, zero-mean Gaussian measurement floor with standard deviation sigma_base times HDOP, applied per fix. The bottom row sums all four into the final perturbed trace, shown hugging a dashed copy of the clean baseline, which clears the lag-1 autocorrelation and one-over-f power-spectral-density verification gates. Clean baseline path kinematically valid · physics-checked in a local metric frame (ENU/UTM) + Bias drift random walk · σ_b (m/min) τ ≈ minutes → hours · iono/tropo/clock + Multipath jitter Ornstein-Uhlenbeck · reversion θ (s⁻¹) τ ≈ 10–60 s · heading-aligned, correlated + Measurement floor Gaussian · σ = σ_base · HDOP white, per-fix · thermal + quantization = Perturbed trace clears lag-1 autocorr (r₁>0.2) and 1/f PSD shape gates - - - clean baseline (reference)

Per-Receiver Error Budgets

The parameters below anchor each process to a hardware class. Use these as defaults, then calibrate against a real error trace from the target receiver.

Receiver class Bias σ_b (m/min) OU reversion θ (s⁻¹) Measurement σ_base (m) Multipath jump prob.
Survey / RTK (L1+L2) 0.05 2.0 0.3 <0.1%
Consumer phone (open sky) 0.2 1.0 1.5 0.5%
Consumer phone (urban canyon) 0.4 0.2 3.0 1.0%
Low-cost / IMU-fused 0.5 0.1 3.0 1.5%

The measurement-noise variance is then scaled by simulated satellite geometry, σpos2=σbase2⋅HDOP2\sigma_{pos}^2 = \sigma_{base}^2 \cdot \text{HDOP}^2, so the same receiver shows tighter fixes under good geometry and wider scatter when satellites cluster.

Prerequisite Check

Confirm the projection stack is present and that pyproj carries its own PROJ data — a mismatched transform pipeline is a silent reproducibility hazard.

python
import numpy as np
import pyproj

print("numpy", np.__version__)        # 1.26.x or newer
print("pyproj", pyproj.__version__)   # pin to 3.x
print("PROJ", pyproj.proj_version_str)  # 9.x bundled with pyproj 3.x

# A metre must be a metre before any noise is applied:
assert pyproj.CRS.from_epsg(4326).is_geographic   # degrees — never perturb here

Fix: Three-Component Noise Injection

The injector consumes a clean path in EPSG:4326, a receiver profile, and a deterministic seed, then returns perturbed geodetic coordinates. Each stochastic process is seeded from one base seed so the whole fleet reconstructs byte-for-byte.

Correct and incorrect placements of noise injection in a trajectory pipeline A four-stage pipeline runs left to right: generate the clean path under its physical constraints, project into a local metric frame, inject noise, and resample to the reporting interval before reprojecting out. Two alternative placements are drawn crossed out. Injecting noise before path generation is wrong because the constraint stage that follows will smooth it away: a jerk limiter and a curvature clip exist precisely to remove high-frequency deviation, so the noise is filtered out and the trace ends up cleaner than the one it started from. Injecting noise after resampling is wrong for the opposite reason: resampling interpolates between fixes, so noise added afterwards is applied to points that have already been smoothed, and the correlation structure between consecutive reported fixes is not the one a receiver produces. A note records the third constraint that fixes the position exactly: noise has to be injected in a projected metric frame, because a displacement expressed in degrees is a different distance at every latitude, and the reprojection back to geographic coordinates is therefore the last step rather than an earlier one. Between the constraint stage and the resampler — and nowhere else generate clean path under physical constraints project into a local metric frame inject noise layered, correlated, in metres resample & reproject to the reporting interval, then back to EPSG:4326 too early the constraint stage filters it straight back out too late resampling has already smoothed the trace And it has to happen in metres. A displacement expressed in degrees is a different distance at every latitude, which is why the reprojection back to EPSG:4326 is the last step rather than an earlier one. Both wrong placements produce a trace that has noise in it and the wrong noise spectrum — and a spectrum is not something anybody inspects by eye, which is why this ordering error survives review so reliably.
Both wrong placements produce a trace with noise in it and the wrong noise spectrum — and nobody inspects a spectrum by eye.
python
import numpy as np
from pyproj import Transformer

def inject_gnss_noise(lon, lat, t, profile, seed):
    """Perturb a clean path with correlated GNSS error. Arrays are 1-D, time in seconds."""
    rng = np.random.default_rng(seed)
    dt = np.diff(t, prepend=t[0] - 1.0)              # per-fix interval

    # --- Step 1: project EPSG:4326 -> local UTM so noise is metric ---
    utm = pyproj.database.query_utm_crs_info(
        "WGS 84",
        pyproj.aoi.AreaOfInterest(lon.min(), lat.min(), lon.max(), lat.max()),
    )[0]
    fwd = Transformer.from_crs("EPSG:4326", utm.code, always_xy=True)
    inv = Transformer.from_crs(utm.code, "EPSG:4326", always_xy=True)
    e, n = fwd.transform(lon, lat)

    # --- Step 2: slow bias drift (random walk, m/min -> m/s) ---
    sb = profile["sigma_b"] / 60.0
    bias_e = np.cumsum(sb * rng.standard_normal(len(t)) * np.sqrt(dt))
    bias_n = np.cumsum(sb * rng.standard_normal(len(t)) * np.sqrt(dt))

    # --- Step 3: OU multipath jitter, exact discretization ---
    theta, sigma_j = profile["theta"], profile["sigma_base"]
    J_e = np.zeros(len(t)); J_n = np.zeros(len(t))
    for k in range(1, len(t)):
        a = np.exp(-theta * dt[k])
        sd = sigma_j * np.sqrt((1 - a**2) / (2 * theta))
        J_e[k] = a * J_e[k-1] + sd * rng.standard_normal()
        J_n[k] = a * J_n[k-1] + sd * rng.standard_normal()

    # --- Step 4: HDOP-scaled high-frequency measurement noise ---
    hdop = profile.get("hdop", 1.0)
    meas = profile["sigma_base"] * hdop
    noise_e = meas * rng.standard_normal(len(t))
    noise_n = meas * rng.standard_normal(len(t))

    # --- Step 5: rare heavy-tailed multipath jumps ---
    jump = rng.random(len(t)) < profile["jump_prob"]
    jump_mag = rng.standard_cauchy(len(t)) * meas  # heavy tail
    noise_e += jump * jump_mag
    noise_n += jump * rng.standard_cauchy(len(t)) * meas * jump

    e_p = e + bias_e + J_e + noise_e
    n_p = n + bias_n + J_n + noise_n
    lon_p, lat_p = inv.transform(e_p, n_p)
    return lon_p, lat_p

python
URBAN = {"sigma_b": 0.4, "theta": 0.2, "sigma_base": 3.0,
         "jump_prob": 0.01, "hdop": 1.8}
lon_p, lat_p = inject_gnss_noise(lon, lat, t, URBAN, seed=20260402)

Verification Step

Injected noise must clear three CI gates before it enters any training or routing pipeline. Each gate catches a specific way naive noise betrays itself.

python
def verify_noise(clean_e, pert_e, dt=1.0):
    err = pert_e - clean_e

    # 1. Autocorrelation decay — white noise drops to ~0 at lag 1.
    r1 = np.corrcoef(err[:-1], err[1:])[0, 1]
    assert r1 > 0.2, f"no temporal correlation (r1={r1:.3f}) — white noise leaked in"

    # 2. Spectral shape — real GNSS error is 1/f at low frequency, not flat.
    psd = np.abs(np.fft.rfft(err - err.mean())) ** 2
    lo, hi = psd[1:5].mean(), psd[len(psd)//2:].mean()
    assert lo > 2 * hi, "flat PSD — error has no low-frequency drift"

    # 3. Kinematic continuity — derived accel must stay vehicle-plausible.
    v = np.gradient(pert_e, dt)
    a = np.gradient(v, dt)
    assert np.abs(a).max() < 4.0, "acceleration spike — perturbed beyond vehicle dynamics"
    return {"r1": r1, "psd_ratio": lo / hi, "amax": np.abs(a).max()}

Wire verify_noise into the same gate that runs on every regression build. The autocorrelation assertion is the one that matters most: it is the single check that makes it impossible for a refactor to silently regress the OU process back to independent Gaussian draws. Calibrate the real GNSS time constant against the target hardware — position error typically decays exponentially with a 10–60 s constant.

Edge Cases & Gotchas

  • Antimeridian and UTM-zone straddling. A trajectory that crosses ±180° longitude, or one that spans two UTM zones, will pick a single zone via query_utm_crs_info and accumulate distortion at the far edge. For long-haul routes, perturb in a local ENU tangent plane anchored at the path centroid instead of a fixed UTM zone, or segment the path per zone and stitch the perturbed pieces.
  • Null Island and unset coordinates. A (0.0, 0.0) fix from a failed lock sits off the African coast and will be perturbed as if real, injecting plausible-looking garbage. Drop or interpolate sentinel (0, 0) and NaN fixes before injection — never let the noise layer dignify a dropout with realistic drift.
  • Float precision at low latitude vs. high latitude. Inverse-projecting metre-scale perturbations back to degrees loses precision differently near the poles, where a degree of longitude collapses. Keep all arithmetic in float64 and inverse-project only once, at export — repeated round-trips through EPSG:4326 accumulate sub-metre bias that the PSD gate will eventually flag.

Frequently Asked Questions

Why not just add Gaussian noise to every fix?

Because independent Gaussian draws have zero temporal autocorrelation, so any Kalman filter or moving average downstream erases them instantly and the model sees effectively clean paths. Real GNSS error drifts and reverts; reserve plain Gaussian noise for the thin HDOP-scaled measurement floor layered under the OU drift, never as the whole model.

How do I choose the OU reversion rate θ for a given vehicle?

Match θ to how long the receiver’s bias persists. A survey-grade unit corrects quickly (large θ, error reverts within a second or two); a phone in an urban canyon drifts for tens of seconds (small θ). Fit the OU autocorrelation decay to a real error trace from the target hardware class, then confirm the synthetic lag-1 autocorrelation lands in the same band.

How do I scale measurement noise with satellite geometry?

Multiply the base variance by HDOP², so σpos2=σbase2⋅HDOP2\sigma_{pos}^2 = \sigma_{base}^2 \cdot \text{HDOP}^2. Drive HDOP from a simulated sky-occlusion model — high values when buildings block satellites, low values in open sky — so the same receiver tightens and widens realistically along the route.

Does adding GPS noise make a trajectory dataset private?

No. Perturbation reduces re-identification risk but does not bound it. Pair it with formal differential privacy mechanisms and a spatial k-anonymity gate, and confirm the perturbation magnitude clears the minimum safe distance for every distinctive trace before export.