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

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

1 Like

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.

2 Likes

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.

1 Like

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?

1 Like

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.

Nice write-up. We build receivers on Septentrio modules (Eview GNSS, so I’m biased), and a few things from our side may help with the “which mode is it in” problem.

  1. Gate R on the solution type, not only on the covariance. The receiver reports the mode on every epoch (standalone, SBAS, DGNSS, RTK float, RTK fixed, PPP). In ROS 2 the septentrio_gnss_driver exposes it in the PVTGeodetic message and in gps_msgs/GPSFix, while NavSatStatus collapses most of these into a few values. A per-mode floor in metres, as you ended up with, is the right shape. I’d also reset your correlation estimate whenever the mode changes, because the error process is different on each side of the switch.

  2. Expect the correlation term to fire on Septentrio too. Multipath and residual atmospheric errors are slow processes. In RTK fixed the magnitude is centimetres, but it is still correlated over minutes. In float it can wander by decimetres while the reported sigma keeps shrinking. Your effective-sample-size correction is the honest way to handle it.

  3. For a cheap calibration, park the antenna on a fixed point with open sky for a few hours, log the raw output and compute the lag-1 autocorrelation per mode. That gives you a prior for r instead of learning it only in stationary windows.

  4. Check where your covariance comes from. The Septentrio driver fills it from the receiver’s PosCovGeodetic block, while some NMEA-based setups derive it from HDOP or GST, which is coarser.

Happy to compare numbers if you run the idle test on a Septentrio unit.

2 Likes

Ahh! @Eviewgnss, thanks, that is the piece I was missing.

I implemented something I’m naming as the “mode-change reset”. It was a real gap: the solution type was already attached to every fix but my window was simply ignoring it. Now with the restart disabled the same sequence reports 1.10 m, still dominated by the mode the receiver had already left. So it was quietly blending two error processes and I had no way to see it.

That also explains the thing that confused me most in the original post, where the same receiver was 160x too pessimistic on one pass and 13.6x too optimistic later the same day. I had treated those as two separate oddities rather than as one receiver in two modes.

The per-mode autocorrelation prior is the part I want to do next, and I cannot without hardware. I have one u-blox. If you are willing to send a unit I will run the idle test on it, publish the per-mode numbers either way, and add a certified config so your customers get tested settings and a one-command check they can run on their own bag. I have done that for the Xsens MTi-680G, a Microstrain and the BNO085 already, they are public in the repo under docs/certified-configs. DM me or create an issue. I’m very active on github hehe

Either way I will post whatever the Septentrio numbers look like if I get hold of one.

Manan

in addition to everything else said, afaik, the explanation is also partially “government regulation”
basically, the GPS satellites don’t give precise enough data to allow the GPS module to know the true position to that degree of accuracy
for example, the position it has could be consistently off by like 1m, even though the deviation is ~0.1m

@solonovamax I have a bag that separates them. Rover parked on one spot for 17.5 minutes, all four wheel encoders reporting zero distance travelled so the antenna really is static, u-blox M9N, no corrections, 1053 fixes, 8 satellites:

reported 1-sigma horizontal              median 2.64 m
displacement between consecutive fixes   median 0.062 m
largest separation between any two fixes        7.49 m

So it is repeatable to about 6 cm from one epoch to the next, and the solution still walks 7.5 m across the run while the antenna does not move at all.

The part I found most useful for telling the two apart: average the fixes in 60 second blocks, about 58 samples each, and the 18 block means still span 6.69 m. If the epoch-to-epoch scatter were white noise, those means would collapse to the order of a centimetre. They do not, so what is left is slowly varying and common mode, which is what @peci1 and @Eviewgnss were pointing at from two different directions earlier in the thread.

On the cause, I do not think Selective Availability is in play any more. It was set to zero in May 2000 and the capability was removed from the GPS III satellites in 2007, and gps.gov still lists that as current policy. What is left is propagation and broadcast data: residual ionospheric delay on a single-frequency receiver, troposphere, and satellite orbit and clock error from the broadcast ephemeris. That is also why differential corrections fix it without anything changing at the satellite, since those errors are shared by receivers that are near each other and cancel in the difference.

Which gets to the same place you did.