Noise Injection & Stochastic Drift for Synthetic Trajectories

Deterministic path generation produces geometrically clean curves, but a clean curve is a tell: real telemetry carries sensor degradation, multipath jitter, mean-reverting bias, and clock skew that idealized paths never show. Part of Trajectory & Movement Simulation, this page covers the specific sub-problem of perturbing a faithful baseline path into a statistically representative one — turning the smooth output of a routing or kinematics stage into data that trains robust models, exercises QA suites against the failure cases that actually occur, and still reconstructs byte-for-byte from a seed. Where the parent area frames the fidelity-versus-privacy trade-off across every movement primitive, here the contract is narrower and sharper: the noise must look like a named sensor’s error distribution, must be applied in the right coordinate space, and must never push a coordinate into a physically impossible state.

Problem Framing

The failure this stage prevents is plausible-but-wrong synthetic data. Three concrete defects recur when noise injection is done naively:

• Coordinate-space corruption. Adding metres of noise directly to latitude/longitude degrees applies a latitude-dependent, anisotropic distortion that nobody intended. A 5 m standard deviation interpreted as degrees is meaningless; interpreted as degrees-times-a-fixed-scale it is correct at the equator and 40% too small at 60° N. This is the same class of bug as silent coordinate reference system drift — the numbers stay finite and the map still renders, so the error survives review.
• Statistically wrong error structure. Independent Gaussian noise per fix has no temporal autocorrelation, so it averages out over a window and a Kalman filter erases it instantly. Real GNSS error is correlated — it drifts, reverts, and occasionally jumps — and a model trained against white noise learns nothing about the failure modes that matter in production.
• Physically impossible motion. Perturbing a path that was carefully constrained by physics-based path generation re-introduces the very violations that stage removed: acceleration spikes, negative speeds, and sub-mechanical turning radii. Noise applied after and never re-checked against kinematics defeats the point of using a physically grounded baseline at all.

The fix is an ordered, reproducible transform: project to a local metric frame, perturb with a correlated stochastic process aligned to the heading, re-validate kinematics, perturb the timestamps, then inverse-project to geodetic coordinates. Order is load-bearing — get it wrong and the errors correlate in ways no downstream consumer expects.

Prerequisites & Toolchain

Noise injection sits between routing/kinematics upstream and serialization downstream, so it needs the same projection and array stack as the rest of the pipeline. Pin versions explicitly — pyproj ships its own PROJ data, and a mismatched transform pipeline is a silent reproducibility hazard.

Copy
python>=3.10
numpy>=1.26,<2.0
scipy>=1.11,<2.0
pyproj>=3.6,<4.0      # bundles PROJ 9.x; do not rely on a system PROJ_LIB
geopandas>=0.14,<0.15
shapely>=2.0,<3.0

Set the environment so transforms and PRNG draws are deterministic across hosts:

bashCopy
export PYTHONHASHSEED=0
export PROJ_NETWORK=OFF        # never fetch grid shifts at runtime; pin them
python -c "import pyproj; print(pyproj.proj_version_str)"   # record this in the seed registry

Treat the noise stage as a stateless transform with one explicit dependency — the seed. Everything else (distribution family, parameters, CRS, PROJ version) is configuration that must be logged alongside the output so a run reproduces exactly, the same discipline that CI/CD integration for spatial data expects from every pipeline stage.

Core Concept: Projection-Aware, Correlated, Heading-Aligned Perturbation

Three ideas distinguish production noise injection from a one-line + np.random.normal().

1. Perturb in a local metric frame. Convert the baseline to a local tangent plane — East-North-Up (ENU) or the appropriate UTM zone — so a metre is a metre everywhere on the path. The OGC WKT Coordinate Reference Systems Standard formalizes this projection contract; in practice you pick a metric CRS at ingestion, do all perturbation there, and inverse-project only at export.

2. Use a correlated process, not white noise. The error model should match the physics of the sensor class:

• Gaussian white noise captures thermal and quantization error — a reasonable floor, but it has zero autocorrelation and washes out under filtering.
• The Ornstein-Uhlenbeck (OU) process captures mean-reverting drift, the dominant signature of low-cost GNSS and IMU bias. It is defined by $dx_t = \theta(\mu - x_t)\,dt + \sigma\,dW_t$, where $\theta$ is the reversion rate, $\mu$ the long-term mean, and $\sigma$ the volatility. Larger $\theta$ pulls the error back toward $\mu$ faster; smaller $\theta$ lets bias wander for longer, which is what makes it survive a downstream smoother.
• Lévy flights / α-stable distributions model heavy-tailed jumps — urban-canyon multipath and satellite-handover spikes. They preserve scale invariance and reproduce real outlier frequencies that Gaussian assumptions badly under-count.

3. Align the noise to the heading. GNSS error is anisotropic: longitudinal (along-track) uncertainty usually exceeds lateral (cross-track) uncertainty because of Doppler smoothing and velocity-estimation bias. Build a diagonal covariance in the vehicle frame, then rotate it into the local ENU frame by the instantaneous bearing $\psi$ before sampling:

$$\Sigma_{\text{ENU}} = R(\psi)\begin{bmatrix}\sigma_\parallel^2 & 0\\ 0 & \sigma_\perp^2\end{bmatrix}R(\psi)^\top,\qquad R(\psi)=\begin{bmatrix}\cos\psi & -\sin\psi\\ \sin\psi & \cos\psi\end{bmatrix}$$

A separate, slowly accumulating drift term sits on top of the per-fix noise and is bounded by a circuit breaker so cumulative bias never diverges over a long run.

Step-by-Step Implementation

The steps below compose into a single deterministic transform. Every stochastic draw comes from one seeded numpy.random.Generator, so identical seeds yield byte-identical output.

Step 1 — Project the baseline to a local metric frame. Load the clean path and transform it from EPSG:4326 to a metric CRS. For city-scale data use the local UTM zone; for a tangent-plane treatment build an ENU transformer around a fixed origin.

pythonCopy
import numpy as np
from pyproj import Transformer

SEED = 20260326
rng = np.random.default_rng(SEED)

# baseline: arrays of lon, lat (deg) and unix-epoch seconds, monotonic in time
lon, lat, t = load_baseline_path()            # shape (N,), (N,), (N,)

to_metric = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)  # UTM 33N
to_geo    = Transformer.from_crs("EPSG:32633", "EPSG:4326", always_xy=True)
east, north = to_metric.transform(lon, lat)   # metres
xy = np.column_stack([east, north])

Step 2 — Derive per-fix heading. The anisotropic covariance needs the bearing at each fix; use the finite-difference direction of travel, falling back to the previous heading where the vehicle is stationary.

pythonCopy
def headings(xy: np.ndarray) -> np.ndarray:
    d = np.diff(xy, axis=0, prepend=xy[:1])
    psi = np.arctan2(d[:, 1], d[:, 0])        # radians, ENU bearing
    psi[0] = psi[1] if len(psi) > 1 else 0.0
    return psi

psi = headings(xy)

Step 3 — Sample a heading-aligned OU process. Integrate the OU SDE with an exact discretization (numerically stable for any step size), then rotate the along/cross-track components into the ENU frame per fix.

pythonCopy
def ou_anisotropic(psi, dt, *, theta, sigma_par, sigma_perp, rng):
    n = len(psi)
    # exact OU update: x_{k+1} = x_k e^{-theta dt} + noise with matched variance
    decay = np.exp(-theta * dt)
    std_par  = sigma_par  * np.sqrt((1 - decay**2) / (2 * theta))
    std_perp = sigma_perp * np.sqrt((1 - decay**2) / (2 * theta))
    par = np.zeros(n)
    perp = np.zeros(n)
    for k in range(1, n):
        par[k]  = par[k-1]  * decay + rng.normal(0.0, std_par)
        perp[k] = perp[k-1] * decay + rng.normal(0.0, std_perp)
    # rotate (along, cross) -> (east, north) by per-fix heading
    e = par * np.cos(psi) - perp * np.sin(psi)
    n_ = par * np.sin(psi) + perp * np.cos(psi)
    return np.column_stack([e, n_])

dt = np.diff(t, prepend=t[0] - 1.0)
ou = ou_anisotropic(psi, dt, theta=0.15, sigma_par=4.0, sigma_perp=2.0, rng=rng)

Step 4 — Add heavy-tailed multipath jumps. Sparingly inject α-stable jumps to reproduce urban-canyon spikes; a Bernoulli mask keeps them rare and the jump magnitude heavy-tailed.

pythonCopy
from scipy.stats import levy_stable

def multipath_jumps(n, *, p_jump, alpha, scale, rng):
    mask = rng.random(n) < p_jump
    mag = levy_stable.rvs(alpha, 0.0, scale=scale,
                          size=n, random_state=rng) * mask
    ang = rng.uniform(0, 2 * np.pi, size=n) * mask
    return np.column_stack([mag * np.cos(ang), mag * np.sin(ang)])

jumps = multipath_jumps(len(xy), p_jump=0.01, alpha=1.6, scale=6.0, rng=rng)

Step 5 — Accumulate bounded drift with a circuit breaker. Keep a separate slow-drift accumulator and clamp it so cumulative bias cannot diverge over a long simulation window.

pythonCopy
DRIFT_MAX_M = 25.0
drift = np.cumsum(rng.normal(0.0, 0.05, size=(len(xy), 2)), axis=0)  # random walk
mag = np.linalg.norm(drift, axis=1, keepdims=True)
over = mag > DRIFT_MAX_M
drift = np.where(over, drift / np.maximum(mag, 1e-9) * DRIFT_MAX_M, drift)

perturbed_xy = xy + ou + jumps + drift

Step 6 — Re-validate kinematics before inverse-projecting. Perturbation can create impossible motion, so differentiate in velocity-acceleration-jerk space and smooth out violations with a Savitzky-Golay pass tuned to the 95th-percentile of real sensor error rather than to hard mechanical limits — over-constraining here silently strips the variance you just added.

pythonCopy
from scipy.signal import savgol_filter

def enforce_kinematics(xy, dt, *, v_max, a_lat_max, smooth_window=7, poly=3):
    v = np.gradient(xy, axis=0) / dt[:, None]
    a = np.gradient(v, axis=0) / dt[:, None]
    speed = np.linalg.norm(v, axis=1)
    a_lat = np.abs(np.cross(v, a)) / np.maximum(speed, 1e-6)
    if (speed > v_max).any() or (a_lat > a_lat_max).any():
        xy = savgol_filter(xy, smooth_window, poly, axis=0)  # preserve topology
    return xy

clean_xy = enforce_kinematics(perturbed_xy, dt, v_max=40.0, a_lat_max=3.5)

Step 7 — Perturb time, then inverse-project. Inject timestamp jitter and a bounded clock-skew random walk, enforce monotonicity so time-series indexing never breaks, simulate packet dropouts, and only then convert back to geodetic coordinates — handing realistic, dropout-aware timing to the temporal synchronization for moving objects stage downstream.

pythonCopy
def perturb_time(t, *, jitter_s, skew_sigma, p_drop, rng):
    jitter = rng.normal(0.0, jitter_s, size=len(t))
    skew = np.cumsum(rng.normal(0.0, skew_sigma, size=len(t)))   # Wiener clock drift
    tp = t + jitter + skew
    tp = np.maximum.accumulate(tp)                  # forbid timestamp inversion
    keep = rng.random(len(t)) >= p_drop             # Bernoulli packet loss
    return tp, keep

tp, keep = perturb_time(t, jitter_s=0.05, skew_sigma=0.01, p_drop=0.02, rng=rng)
out_lon, out_lat = to_geo.transform(clean_xy[keep, 0], clean_xy[keep, 1])

For sensor-class parameter tables and hardware-specific error budgets that calibrate these knobs, the drill-down on adding realistic GPS noise to synthetic vehicle trajectories provides per-receiver $\sigma$, $\theta$, and jump-rate values with validation benchmarks.

Validation & Testing

Noise injection is only trustworthy if it is both correct (right statistics) and reproducible (right bytes). Gate both in CI.

Reproducibility. Hash the full output for a fixed seed and assert it never changes across hosts:

pythonCopy
import hashlib

def trace_hash(lon, lat, t):
    buf = np.concatenate([lon, lat, t]).astype("<f8").tobytes()
    return hashlib.sha256(buf).hexdigest()

assert trace_hash(out_lon, out_lat, tp[keep]) == GOLDEN_HASH

Statistical fidelity. Verify the injected error reproduces the target moments and autocorrelation, and that spatial structure survives — Moran’s I and the rest of the realism metrics and evaluation suite turn “looks noisy” into a number with a tolerance.

pythonCopy
residual = np.column_stack([
    *to_metric.transform(out_lon, out_lat)]) - clean_xy[keep]
err = np.linalg.norm(residual, axis=1)

assert abs(err.mean()) < 1.0                       # near-zero bias after drift clamp
assert 3.0 < err.std() < 6.0                       # matches target sensor sigma band
# lag-1 autocorrelation must be positive — proves OU, not white noise
r1 = np.corrcoef(err[:-1], err[1:])[0, 1]
assert r1 > 0.2

Kinematic and temporal invariants. Assert no impossible motion or time inversion slipped through:

pythonCopy
assert np.all(np.diff(tp[keep]) > 0)               # strictly monotonic timestamps
v = np.linalg.norm(np.gradient(clean_xy[keep], axis=0)
                   / np.diff(tp[keep], prepend=tp[keep][0]-1)[:, None], axis=1)
assert v.max() <= 40.0 + 1e-6                       # speed envelope respected

Pin a golden dataset with a fixed PRNG seed and run these assertions on every pipeline change; a moment that moves outside its tolerance band is a regression, not noise.

Performance & Scale Considerations

• Vectorize the projection, loop only the OU recursion. Transformer.transform is fully vectorized — pass whole arrays, never per-point. The OU update in Step 3 is genuinely sequential, but it is a single O(N) pass over scalars; if it dominates, lift it into Numba or rewrite as a cumulative-product recurrence rather than reaching for multiprocessing.
• Parallelize across agents, not within a trace. A fleet of independent agents is embarrassingly parallel. Derive each agent’s seed deterministically from one base seed (np.random.default_rng(base).integers(...) per agent), so workers run concurrently while the merged fleet stays byte-identical — the same seed-derivation discipline used in Markov-chain routing.
• Stream long traces in chunks. Carry the OU value, drift accumulator, and clock-skew state across chunk boundaries so a multi-hour trace never materializes fully in memory; the exact-discretization OU update makes chunk boundaries seamless.
• Cache the transformer. Construct Transformer objects once and reuse them; building one per call re-parses PROJ pipelines and dominates runtime on short traces.

Failure Modes & Troubleshooting

Noise added in geographic degrees

Symptom: error magnitude varies with latitude; a “5 m” σ produces visibly different spread in the north and south of a dataset. Root cause: sampling noise directly on lon/lat instead of in a metric CRS. Fix: always project to UTM or ENU first (Step 1), perturb in metres, and inverse-project only at export.

White noise that vanishes downstream

Symptom: a Kalman filter or moving average completely removes the injected error and models see effectively clean paths. Root cause: independent per-fix Gaussian draws with zero temporal autocorrelation. Fix: use the OU process (Step 3) and assert lag-1 autocorrelation r1 > 0.2 in CI so white noise can never regress in unnoticed.

Impossible motion after perturbation

Symptom: velocity spikes, negative speeds, or sub-mechanical turning radii in the perturbed output. Root cause: noise applied to a kinematically clean path with no re-validation. Fix: run the velocity-acceleration-jerk filter of Step 6, and calibrate thresholds to the 95th percentile of real sensor error so filtering does not silently erase the variance you injected.

Unbounded drift over long runs

Symptom: cumulative bias grows without limit and late-trace coordinates wander hundreds of metres off the path. Root cause: an OU $\theta$ near zero plus an unclamped random-walk drift term. Fix: keep the drift accumulator separate, clamp it with the circuit breaker of Step 5, and raise $\theta$ so the mean-reverting term actually reverts.

Timestamp inversion breaking replay

Symptom: downstream time-series indexing or causal feature extraction throws on non-monotonic time. Root cause: additive jitter and clock skew can reorder near-simultaneous fixes. Fix: apply np.maximum.accumulate to enforce monotonicity (Step 7) and assert np.diff(t) > 0 before export.

Privacy leakage through under-perturbation

Symptom: an adversarial linkage test re-identifies a distinctive home-to-work trace from the synthetic output. Root cause: perturbation magnitude below the spatial-generalization threshold, so a rare real path remains recognizable. Fix: apply differential privacy mechanisms — ensure perturbation exceeds the minimum safe distance for $k$-anonymity, track the composed $(\varepsilon, \delta)$ budget, and confirm training utility survives with the statistical-fidelity checks above.

Frequently Asked Questions

Related

• Physics-Based Path Generation — the kinematically clean baselines this stage perturbs, and the velocity/jerk envelopes its re-validation enforces.
• Markov Chain Routing Models — the routing layer upstream; keep its transition-matrix perturbation decoupled from per-fix GNSS noise.
• Temporal Synchronization for Moving Objects — consumes the jittered, dropout-aware timestamps this stage emits and aligns them to simulation ticks.
• Adding Realistic GPS Noise to Synthetic Vehicle Trajectories — per-receiver error budgets, σ/θ tables, and validation benchmarks for the parameters above.
• Privacy-Preserving Generation Frameworks — the differential-privacy budgets and $k$-anonymity gates that bound leakage from perturbed traces.