GPS says ~5m accuracy but the fixes agree to 9 cm. Anyone else seen this?

[manankharwar]

So I’ve been chasing something for a couple of days and I figured I’d just post it, partly in case it helps someone and partly because I’m not totally sure I’m reading it right.

I’ve got a u-blox on a small rover doing 1 Hz fixes. The receiver puts ~5m horizontal 1-sigma in position_covariance, which honestly seemed believable for where I was driving. But then I looked at how much the fixes actually move relative to each other, and the median second difference is around ~0.1 m. The worst one in the whole run was about 1 m. If the noise were really white at ~5m I’d have expected that number to be somewhere around 14 m.

So there’s roughly a factor of 160 between what it says and how it behaves.

My guess, and this is the bit I’d like a sanity check on, is that the receiver smooths internally. So the 5m is the absolute error, mostly multipath, and it’s an honest number. But the fix-to-fix noise is centimetres because consecutive fixes are now heavily correlated. Which would mean the covariance is true, it’s just not the thing a Kalman filter is asking for, since the filter assumes the noise is white and independent.

Anyway, if that’s right then it explains my gate. My NIS sits around 0.03 where it should be near 3, so the chi-squared threshold ends up miles above anything that ever actually happens and nothing can trip it.

Two things I genuinely don’t know:

1. Is this a u-blox thing specifically, or do other receivers do the same? I’ve only got the one so I can’t really tell.
2. And what do people actually do about it in practice? Shrinking R feels wrong to me, because then you’d just track the multipath bias and confidently report centimetre accuracy while sitting metres off. I’ve seen a slowly varying GPS bias state mentioned as the proper fix but I haven’t tried it yet.

Here’s the script I used if anyone wants to check theirs. Works on any bag with a NavSatFix, nothing else needed:

#!/usr/bin/env python3
"""python3 gps_cov_check.py <bag_dir>"""
import sys, math, statistics, rosbag2_py
from rclpy.serialization import deserialize_message
from rosidl_runtime_py.utilities import get_message

bag = sys.argv[1]
r = rosbag2_py.SequentialReader()
r.open(rosbag2_py.StorageOptions(uri=bag, storage_id=""), rosbag2_py.ConverterOptions("", ""))
topic = next(t.name for t in r.get_all_topics_and_types()
if t.type == "sensor_msgs/msg/NavSatFix")
r.set_filter(rosbag2_py.StorageFilter(topics=[topic]))
M = get_message("sensor_msgs/msg/NavSatFix")

lat, lon, sig = [], [], []
while r.has_next():
_, data, _ = r.read_next()
m = deserialize_message(data, M)
if m.position_covariance_type == 0:
continue
lat.append(m.latitude); lon.append(m.longitude)
sig.append(math.sqrt(max(m.position_covariance[0], 0.0)))

R = 6371000.0
xs = [math.radians(v - lon[0]) * R * math.cos(math.radians(lat[0])) for v in lon]
ys = [math.radians(v - lat[0]) * R for v in lat]
d2 = sorted(math.hypot(xs[i+1] - 2*xs[i] + xs[i-1], ys[i+1] - 2*ys[i] + ys[i-1])
for i in range(1, len(xs) - 1))
declared = statistics.median(sig)
observed = d2[len(d2) // 2]
expected = math.sqrt(6.0) * 1.1774 * declared

print("declared 1-sigma %.2f m" % declared)
print("median |2nd difference| %.3f m" % observed)
print(" if that were white %.1f m" % expected)
print("ratio %.0fx smoother" % (expected / max(observed, 1e-9)))

Would be curious what ratio other people get and which receiver you’re on, especially if anyone’s running RTK fixed or something that doesn’t smooth its output. Thanks

[peci1]

I think this is just the common problem with GNSS. The error estimates are nowhere close to white noise. They are heavily correlated even on non-RTK sensors - the phase locks themselves introduce a quite large time dependency which breaks the I.I.D. assumption of KF or even factor graphs. As the signal from which you compute the position is digital (the code modulated on the waveforms), the only possible source of a white noise (or a random walk) would be clock drift and multipath.

It seems to me you’re right that the receivers try to express the bias from the absolute pose in the “covariance” outputs (it fits both Emlid and Septentrio receivers we used).

We don’t have a good solution for this. We just use the reported covariance as if it were a Gaussian and that gets fused to our pose estimate covariance (that’s also what robot_localization does).

We’re (very) slowly getting towards a framework that could at least reject unreliable signals when the robot has substantial pitch or roll (because no receiver has IMU-assisted horizon filter). That will hopefully remove some of the multipath stuff we’re fighting with. Similarly, there was also a paper from Norlab a few years ago about rejecting signals whose satellites are hidden behind lidar-detected obstacles. But that doesn’t help with this problem on the theoretical side.

Slow-moving bias might be good, but it really depends on the geometry of the scene which can change abruptly.

I guess the only reliable solution would be to model the whole SPP or RTK problem as separate factors down to the phase lock mechanism, but I have no idea how demanding that would be (both CPU-demanding and implementation-demanding). I have now one student who plays with GNSS++ library which seems like a nice successor to rtklib. But we don’t know yet what the outcome will be. One thing I’m sure with is we don’t want to compete Septentrio/uBlox, so we need to leave the hard parts of the problem on the sensors’ APIs and only process the signals they report (i.e. we don’t aim to do multipath detection, cycle slips, interference etc.).

As a good test whether you handle this problem correctly you could take a problem with a motion model and GNSS. When the robot is idle, the motion model reports zero motion, but the GNSS might be off/jumping. As long as the last position before going idle more or less agrees with the current GNSS “covariance”, the pose estimate should not drift towards the GNSS mean. That I’d consider a holy grail of GNSS integration in localization.

[manankharwar]

Thanks, that’s really helpful. The phase lock thing is the part I didn’t have. I was thinking of it as the receiver smoothing white noise after the fact, which is clearly the wrong picture here I guess, and it also explains why you see the same behaviour on Emlid and Septentrio. Good to know it isn’t just a u-blox quirk.

So I went and ran your idle test. My filter failed it, pretty badly, and chasing down why turned out to be your phase lock point taken one step further. Writing it up here because I think the second half is more general than my setup.

Worst case, rover parked for a minute with the wheel encoders confirming it hadn’t moved:

GNSS reported position moved 9.76 m, and the filter followed it for 10.16 m

So basically 100% tracking of the receiver’s wander, which is exactly the thing you said shouldn’t happen.

The first thing I found is a bit embarrassing and probably obvious to you: ZUPT pins velocity, not position. I fuse [vx=0, vy=0, wz=0] when the wheels say stationary, and that constrains the velocity states and nothing else. Position process noise just carries on growing between fixes exactly as if the robot were driving. So the covariance stays big, the gain stays high, and every fix drags the estimate. “I’m not moving” and “my position hasn’t changed” are the same physical fact but they’re different mathematical statements, and the filter only ever heard the first one.

Holding position process noise down while the wheels say stationary got the average excursion across 11 parked windows from 3.49 m to 1.64 m. Then it stopped helping, and that’s the interesting bit.

This is where your phase locks come back. A KF assumes the measurement noise is white, so N fixes of a fixed point should shrink the position uncertainty like sqrt(N). But if consecutive fixes are largely the same error repeated, then a minute of 1 Hz fixes isn’t 60 independent looks at a stationary point. It’s closer to two. So the filter ends up far more confident than the geometry actually supports, and then it carries that over-confidence into the next leg of the run.

So instead of picking a distrust factor I tried measuring it. When the wheels say stationary, every fix is sampling the same physical point, which as far as I can tell is the one moment in a run where you can check a sensor’s self-reported covariance against evidence instead of against a constant you made up. Two things fall out of those fixes, per axis:

• magnitude: how far they actually spread, divided by the sigma the receiver claims, squared. A receiver that’s as good as it says gets 1 and nothing happens to it.
• correlation: the lag-1 autocorrelation r. Effective sample size is N(1-r)/(1+r), so R needs to carry a factor of (1+r)/(1-r) for the filter’s own posterior to mean what it says.

Multiply the two, cap it, and that’s the whole thing.

What surprised me was which term did the work. On that 57 s window the magnitude term sat at exactly 1.0 the entire time. The receiver was declaring 21 to 45 m while actually spreading 0.4 to 3.3 m, so it wasn’t over-confident at all, it was pessimistic. All of the correction came from the correlation term, which measured 0.63 climbing to 0.985. That’s an effective sample size factor of 5 going up past 100.

That window went from 10.16 m to 0.10 m. Across all 11 windows: 3.49 m with neither of these, 1.64 m with just the process noise one, 0.74 m with both. It helped in 8 of the 11 and hurt in none. Loop closure over the whole run didn’t change either way, and it didn’t reject any extra fixes, which was the thing I was worried about.

One bit worth flagging, since you mentioned Septentrio and Emlid. The correlation term fires on an honest receiver too, and I think it should. An RTK unit correctly reporting 2 cm still has errors correlated over minutes, so a filter that sits for a minute and averages sixty of them is wrong about how much it knows, even though every individual covariance it was handed was accurate. That’s not a bad-receiver problem, it’s a “we’re all assuming the noise is white” problem. I wouldn’t have gone looking for it if you hadn’t named the mechanism.

Also, and this complicates my original post a bit: the same receiver is wrong in both directions depending on what mode it’s in. The 160x I opened with was it being far too pessimistic. Later the same day, after it picked up SBAS, it declared 0.076 m 1-sigma while its own position scattered about 1.03 m over 70 stationary fixes. So 13.6x optimistic that time, and that direction is the more dangerous one, because R is built from that number and the gate is judged against the same R. I ended up putting a floor on it in metres rather than a multiplier, since a multiplier tuned for one mode is wrong the moment the receiver switches to the other.

Usual caveats: one robot, one receiver, 11 windows, and 0.10 m isn’t zero. I’d genuinely be curious whether the correlation term does anything on a Septentrio, because yours may well be better behaved than mine and then it should measure close to 1 and leave you alone. If it doesn’t, that’d be worth knowing.

[peci1]

Great analysis!

Just a note: wheel odometry 0 doesn’t mean the robot has to be stationary. Many robots have non-negligible dynamics, they can skid, slip, jump… If you want to implement some extra behavior for static states, you need a separate estimator that combines multiple factors (at least wheels + IMU).

Regarding the correlation: could you explain a bit more what do you correlate with what? Is it just the subsequent samples, and only in stationary mode?

[manankharwar]

Both fair.

You’re right that wheels at zero isn’t stationary. That bit me for real: This past weekend, an encoder died mid-run during testing and kept reporting zero while the robot drove. The filter believed it was parked and then refused to let GNSS move it, recovering about 7 m of the 20 m actually driven.

My fix is weaker than what you’re describing. Not a better detector, just a cross-check that revokes the parked state, using GNSS since it doesn’t depend on wheels. While the rover is parked, I track displacement from the first fix and path length summed fix to fix, and gate that on the ratio.

So it’s revocation, not detection, it’s useless indoors, and it won’t catch slip or jump, only sustained motion. Wheels plus IMU is the better answer and I don’t have one.

On the correlation: lag-1 autocorrelation of the GNSS fix sequence itself, per axis, parked windows only. Nothing correlated against another sensor.

It only works parked because the true position is constant, so each fix’s deviation from the window mean is the measurement error. That’s what you can’t get while moving without a truth reference. Parked, you get the error sequence free and can just ask whether consecutive errors are independent.

var  = (sum(x^2) - sum(x)*mean) / (n-1)
r    = (sum(x[i]*x[i-1])/(n-1) - mean^2) / var
R   *= (1 + r) / (1 - r)

r clamped to [0, 0.99], applied per axis as R' = D R D so the receiver’s X/Y correlation survives. There’s a magnitude term multiplying it too, but on the window I quoted that sat at exactly 1.0 throughout, so all the correction came from r going 0.63 to 0.985.