Modeling GPS Multipath Error in Urban Canyons

Synthetic urban trajectories perturbed with independent white Gaussian jitter look nothing like real downtown GPS, because the dominant error source between tall buildings is multipath — a slowly wandering, position-dependent bias that white noise cannot reproduce.

Part of Noise Injection & Stochastic Drift: that stage degrades clean generated paths into sensor-realistic ones, and this page addresses the specific failure that appears the moment a trajectory enters a dense urban core — a noise model calibrated on open-sky data undershoots the magnitude, the correlation, and the spatial structure of the error a receiver actually produces in a canyon.

Root Cause: Reflected Signals Produce Correlated, Biased Error

Open-sky GNSS error is close to zero-mean and nearly white: each fix is a fresh, roughly independent draw around the true position, so adding independent Gaussian noise per sample is a fair approximation. In an urban canyon the geometry breaks both assumptions. Tall buildings block the direct line of sight to low-elevation satellites and reflect their signals off glass and stone, so the receiver locks onto a path that travelled farther than the straight-line distance. That extra path length enters the range measurement as a positive bias, and because the reflecting surfaces and blocked satellites change slowly as the vehicle moves down a street, the resulting position error is temporally correlated and position-dependent — it drifts smoothly over seconds and tens of metres rather than resampling every fix.

White Gaussian noise εt∼N(0,σ2)\varepsilon_t \sim \mathcal{N}(0, \sigma^2) with E[εtεt−1]=0\mathbb{E}[\varepsilon_t \varepsilon_{t-1}] = 0 therefore gets three things wrong at once: it is zero-mean (multipath is biased), it is uncorrelated (multipath persists across fixes), and it is spatially uniform (multipath is severe on a narrow street and mild in an open plaza). A first-order Gauss-Markov process captures the first two:

Street-canyon multipath geometry and the systematic position bias it produces A street is drawn in cross-section between two tall buildings, with a receiver on the pavement. Two satellites are shown. The first is high in the sky and its signal reaches the receiver directly; that measurement is good. The second is low toward the horizon and its direct path is blocked by the building, so the only signal that arrives is one that has bounced off the opposite façade. The reflected path is longer than the direct one, and the extra distance is marked; the receiver has no way to know the signal was reflected, so it treats that longer travel time as genuine distance to the satellite. Because every blocked satellite is blocked on the same side and every reflection arrives from the same façade, the errors do not cancel: the position solution is pulled consistently toward the open side of the street, and the displacement is drawn. A panel records the two consequences that follow. First, the error is a bias rather than noise, so averaging more fixes does not reduce it and a longer dwell makes the wrong answer more confident. Second, the error is a function of position and street orientation, so the same receiver passing the same point on a different day reproduces it — which means a synthetic model of it must be keyed on geometry rather than drawn independently per fix. Multipath is a bias keyed on geometry, not noise keyed on time street, ~18 m wide receiver high satellite direct — good low satellite direct path blocked reflection off the façade the reflected path is longer — the receiver cannot tell solved position Two consequences averaging more fixes does not reduce it — a longer dwell makes the wrong answer more confident it reproduces at the same place on a different day, so model it keyed on geometry, not per fix Every blocked satellite is blocked on the same side and every reflection arrives from the same façade, so the errors do not cancel — which is the whole difference between a bias and noise.
Every blocked satellite is blocked on the same side, so the errors do not cancel — which is the whole difference between a bias and noise.

bt=ϕ bt−1+ηt,ηt∼N(0,ση2),0<ϕ<1,b_t = \phi\, b_{t-1} + \eta_t, \qquad \eta_t \sim \mathcal{N}(0, \sigma_\eta^2), \qquad 0 < \phi < 1,

where the bias state btb_t decays toward zero with correlation coefficient ϕ\phi (larger ϕ\phi means longer-lived drift) and is driven by process noise ηt\eta_t. The stationary variance is σb2=ση2/(1−ϕ2)\sigma_b^2 = \sigma_\eta^2 / (1 - \phi^2), and the correlation time is τ=−Δt/ln⁡ϕ\tau = -\Delta t / \ln\phi. Making the driving variance and correlation depend on where the vehicle is — high on a narrow canyon segment, low in the open — adds the spatial structure. This is a strict generalization of the open-sky model that adding realistic GPS noise to synthetic vehicle trajectories establishes: set ϕ→0\phi \to 0 and the spatial term to a constant and the Gauss-Markov process collapses back to white noise.

Urban-canyon multipath: reflected signals and correlated position bias versus white noise Left panel is a street cross-section with two tall buildings forming a canyon. A satellite high overhead sends a direct signal to a vehicle on the street. A second, low-elevation satellite is blocked by a building, and its signal reflects off the opposite facade before reaching the receiver, drawing a longer bent path. A cloud of true positions along the street is offset by a smoothly drifting bias arrow toward the reflecting wall, illustrating position-dependent correlated error. Right panel plots autocorrelation against lag: white Gaussian noise drops to zero after lag zero, while the Gauss-Markov multipath process decays gradually as phi to the power of lag, showing persistence across many fixes. Canyon geometry & drifting bias Error autocorrelation building building street direct low-elev reflected drifting bias bₜ lag (fixes) autocorr white noise → 0 Gauss-Markov φ^lag persists
Reflected low-elevation signals add a slowly drifting, position-dependent bias. Its autocorrelation decays as a power of the lag, whereas white noise is uncorrelated after lag zero.

Minimal Reproducer: Show White Noise Fails the Autocorrelation Test

Confirm the gap by measuring the autocorrelation of a real urban error series and comparing it to what white Gaussian noise produces. Pin the toolchain so the RNG stream is reproducible.

# requirements.txt — pin majors, let patch releases float
numpy==1.26.*
scipy==1.11.*
pyproj==3.*

python
import numpy as np

def lag1_autocorr(err: np.ndarray) -> float:
    """Lag-1 autocorrelation of a scalar error series (e.g. cross-track offset)."""
    e = err - err.mean()
    denom = np.sum(e * e)
    return float(np.sum(e[1:] * e[:-1]) / denom) if denom > 0 else 0.0

rng = np.random.default_rng(7)
white = rng.normal(0.0, 8.0, size=600)          # open-sky-style jitter, sigma 8 m
# Real downtown cross-track error is biased and sticky:
real_urban_lag1 = 0.93                            # measured from a canyon drive log

print(round(lag1_autocorr(white), 3))            # ~ 0.0  -> uncorrelated
print(real_urban_lag1)                            # ~ 0.93 -> strongly correlated
print(round(white.mean(), 2))                     # ~ 0.0  -> zero-mean, but urban error is biased

White noise reports a lag-1 autocorrelation near zero and a mean near zero; a canyon drive log reports autocorrelation above 0.9 and a non-zero mean. Any model whose synthetic error autocorrelation does not approach the measured value will not fool a detector that inspects error dynamics.

Fix: Position-Dependent Gauss-Markov Bias with Occasional Jumps

Model the urban error as the sum of a slowly drifting Gauss-Markov bias whose intensity depends on local canyon severity, plus a small white measurement term, plus rare discrete jumps when the satellite set changes. Work in a metric CRS so the noise is in metres, never in EPSG:4326 degrees.

Expected horizontal GPS error against canyon height-to-width ratio The horizontal axis is the ratio of building height to street width, from a wide boulevard on the left to a narrow canyon between towers on the right. The primary curve, read against the left axis, is the expected horizontal error in metres, derived from the fraction of the sky dome visible between the façades and the dilution of precision that follows from it. It is close to open-sky accuracy at low ratios and rises steeply once the ratio passes about one, because the satellites that remain visible are all high in the sky and therefore give a poor geometric spread. A second curve, read against the right axis, is the number of usable satellites, which falls as the canyon narrows. Three markers label representative street types: a boulevard, a typical European city street, and a financial-district canyon, with the expected error printed for each. The note underneath makes the modelling point: the error profile is a function of the built form and can be computed from a building-height layer the pipeline already has, so it does not need to be a free parameter, and computing it means the same street produces the same error every time an agent passes through it — which is the property that makes the synthetic trace behave like a real receiver. The error profile follows from the built form — so compute it, don't tune it 0.25 0.75 1.25 1.75 2.25 2.75 0 5 10 15 20 25 building height ÷ street width expected horizontal error (m) 0 4 8 12 usable satellites boulevard · 5 m typical city street · 13 m tower canyon · 24 m expected horizontal error (left) usable satellites (right) Derived from the visible sky fraction between two façades and the dilution of precision that follows, with 3.2 m open-sky accuracy. It can be computed from a building-height layer the pipeline already has — and computing it means the same street produces the same error on every pass, which is what makes the trace behave like a receiver.
Computed from the visible sky fraction, so the same street produces the same error on every pass.

Step 1 — A canyon-severity field

Drive the noise from a per-position severity in [0,1][0, 1]: high on narrow streets flanked by tall buildings, low in open space. In production this comes from a building-height raster or sky-view factor; here it is any callable over the path.

python
import numpy as np

def severity_along(path_xy: np.ndarray, sky_view: np.ndarray) -> np.ndarray:
    """Map a per-sample sky-view factor (1=open, 0=fully enclosed) to canyon severity."""
    sky_view = np.clip(sky_view, 0.0, 1.0)
    return 1.0 - sky_view                              # 0 open plaza ... 1 deep canyon

Step 2 — Correlated, position-dependent bias process

Integrate a first-order Gauss-Markov bias whose correlation ϕ\phi and driving standard deviation scale with severity, so the error is stickier and larger in the canyon and relaxes back to open-sky behaviour in the open.

python
import numpy as np

def multipath_bias(severity: np.ndarray, dt: float, *, seed: int,
                   phi_open: float = 0.6, phi_canyon: float = 0.97,
                   sigma_open: float = 1.5, sigma_canyon: float = 12.0
                   ) -> np.ndarray:
    """Return an (N, 2) correlated bias in metres, driven by local canyon severity."""
    rng = np.random.default_rng(seed)
    n = len(severity)
    b = np.zeros((n, 2))
    for t in range(1, n):
        s = severity[t]
        phi = phi_open + (phi_canyon - phi_open) * s        # longer memory in canyon
        sigma_eta = sigma_open + (sigma_canyon - sigma_open) * s
        eta = rng.normal(0.0, sigma_eta, size=2)
        b[t] = phi * b[t - 1] + eta                         # Gauss-Markov update
    # Bias points preferentially away from the reflecting facade: add a mild along-severity
    # offset so deep-canyon segments carry a persistent positive component, not just zero-mean drift.
    b += 0.4 * severity[:, None] * sigma_canyon * np.array([1.0, 0.0])
    return b

Step 3 — Satellite-set jumps and the final observation

Add rare discrete jumps where the visible constellation changes (turning a corner, passing a gap), then combine bias, jumps, and a small white term into the observed position. Keep the timestamps untouched — resampling and clock alignment belong to the temporal synchronization for moving objects stage, not here.

python
import numpy as np

def apply_urban_gps(clean_xy: np.ndarray, severity: np.ndarray, dt: float, *,
                    seed: int, jump_rate: float = 0.01, sigma_white: float = 1.0
                    ) -> np.ndarray:
    """clean_xy: (N,2) metric truth from the physics stage -> observed urban positions."""
    rng = np.random.default_rng(seed)
    bias = multipath_bias(severity, dt, seed=seed)
    jumps = np.zeros_like(clean_xy)
    active = np.zeros(2)
    for t in range(len(clean_xy)):
        if rng.random() < jump_rate * (0.5 + severity[t]):   # more re-locks in canyon
            active = rng.normal(0.0, 6.0 + 10.0 * severity[t], size=2)
        jumps[t] = active                                    # step persists until next re-lock
    white = rng.normal(0.0, sigma_white, size=clean_xy.shape)
    return clean_xy + bias + jumps + white

Because the bias is driven off the clean geometry from the physics-based path generation stage, run this perturbation after the kinematic envelope is enforced — the correlated drift will inflate finite-difference derivatives, so a path that sat on its jerk limit before noise will breach it after.

Verification Step: Gate Autocorrelation and Spatial Structure

Assert that synthetic urban error reproduces both the temporal correlation and the position dependence measured from real drives, within tolerance.

python
import numpy as np

def test_urban_error_is_correlated_and_position_dependent(
        clean_xy, observed_xy, severity, *, target_lag1=0.90, tol=0.08):
    err = observed_xy - clean_xy
    # 1. Temporal correlation must approach the measured urban value, not ~0.
    e = err[:, 0] - err[:, 0].mean()
    lag1 = float(np.sum(e[1:] * e[:-1]) / np.sum(e * e))
    assert lag1 >= target_lag1 - tol, f"error not sticky enough: lag1={lag1:.2f}"

    # 2. Error magnitude must be larger in canyon segments than open ones.
    mag = np.linalg.norm(err, axis=1)
    hi = mag[severity > 0.7].mean()
    lo = mag[severity < 0.3].mean()
    assert hi > 1.8 * lo, f"no position dependence: canyon {hi:.1f} vs open {lo:.1f}"

    # 3. Sanity: not degenerate / not exploding.
    assert np.isfinite(mag).all() and mag.max() < 150.0, "implausible error magnitude"

The autocorrelation assertion is the load-bearing one: white Gaussian noise cannot pass it at any variance, so a green gate proves the model captures error dynamics and not merely error magnitude. The position-dependence assertion prevents a well-tuned but spatially uniform model from sneaking through.

Edge Cases & Gotchas

Degrees are not metres — project first. Applying a metre-valued bias to EPSG:4326 coordinates adds roughly 8 metres of longitude error at the equator but wildly more near the poles, because a degree of longitude shrinks with latitude. Always transform the clean path into a UTM zone or EPSG:6933 before adding noise, then transform back; a bias generated in degrees produces error that is anisotropic and latitude-dependent for entirely the wrong reason.

Correlation time versus sample rate. The correlation coefficient ϕ\phi is tied to the sampling interval through τ=−Δt/ln⁡ϕ\tau = -\Delta t / \ln \phi. A ϕ=0.97\phi = 0.97 calibrated at 1 Hz implies a very different physical memory at 10 Hz, so if you resample the trajectory you must recompute ϕ\phi from the target correlation time, not carry the raw coefficient across. Carrying it unchanged makes the drift look implausibly persistent or implausibly brief after a rate change.

Antimeridian and UTM-zone seams. A vehicle crossing a UTM zone boundary or the antimeridian while a bias state is active will have that bias applied in two different projected frames, injecting a spurious jump at the seam. Keep the whole canyon segment in a single local projected CRS whose central meridian sits inside the region, and only stitch frames at points where the bias state has relaxed near zero.

Frequently Asked Questions

Why is white Gaussian noise inadequate for urban GPS?
White Gaussian noise is zero-mean, uncorrelated, and spatially uniform, but urban multipath is biased, temporally correlated over many fixes, and far worse on narrow streets than in open space. A detector or model that inspects error dynamics sees the difference immediately: real urban error drifts smoothly and sticks, while white noise resamples independently every fix. White noise is only a fair approximation of open-sky conditions.
What does the correlation coefficient phi actually control?
Phi sets how long the bias state remembers its previous value in a first-order Gauss-Markov process. A phi near zero decays instantly and reproduces white noise, while a phi near one drifts slowly and persists across many fixes, which is the urban-canyon regime. It maps to a physical correlation time through tau equals minus delta-t divided by the natural log of phi, so it must be recomputed whenever the sampling rate changes.
How do I obtain a canyon-severity field without ray-tracing every satellite?
A sky-view factor derived from a building-height raster is a good, cheap proxy: it measures the fraction of the sky hemisphere visible from a point, which correlates strongly with multipath severity. You can approximate it from open building footprint and height data with a simple horizon scan, then map one minus sky-view to severity. Full satellite ray-tracing is more accurate but rarely necessary for generating realistic-looking synthetic error.
Should multipath noise be added before or after enforcing kinematic limits?
After. The clean, physically feasible path is produced first, then noise is layered on, because the correlated bias and the discrete jumps both inflate the finite-difference acceleration and jerk reconstructed from the noisy positions. A path that sat exactly on its jerk limit will breach it once noise is applied, so enforce the kinematic envelope on the clean path and give any downstream jerk gate a slack factor to absorb the perturbation.