GNSS Velocity from Doppler

Overview¶

Besides the pseudorange and the carrier phase, a GNSS receiver measures the Doppler shift of every signal it tracks: the rate at which the relative motion of satellite and receiver stretches the carrier. One Doppler observation is a range-rate measurement along the line of sight, so a handful of satellites determines all three components of the receiver velocity -- and, as a by-product, the drift of the receiver clock.

DopplerFactor puts that measurement into a factor graph. Its keys are the receiver velocity at epoch $k$ and the receiver clock bias at epochs $k-1$ and $k$ , because the clock drift is not a separate state but the difference of two adjacent clock biases:

e^{T}(v_s - v_r) + c\left(\frac{b_k - b_{k-1}}{\Delta t} - \dot b_s\right) + \text{Sagnac rate} \; - \; (-\lambda D).

(1)

This notebook downloads a compact, prepared real 1 Hz dataset from the same open-sky Septentrio mosaic-X5 receiver as Part 2 of RtkAndPppExample.ipynb. The antenna sits on a surveyed static marker, so its true velocity is exactly zero and every metre per second we estimate is error. The prepared dataset keeps the example independent of the RINEX front-end. To process your own static-receiver files, see Preparing another RINEX dataset.

try:
    import google.colab  # noqa: F401
    %pip install --quiet numpy plotly gtsam-develop
    !mkdir -p data
    !wget -q https://raw.githubusercontent.com/borglab/gtsam/develop/python/gtsam/examples/data/doppler_velocity_example.npz -O data/doppler_velocity_example.npz
except ImportError:
    pass

from pathlib import Path

import numpy as np
import plotly.graph_objects as go

import gtsam
from gtsam import symbol
DATA_FILE = Path("data/doppler_velocity_example.npz")
print("prepared data:", DATA_FILE)

prepared data: data/doppler_velocity_example.npz

The measurements¶

The prepared dataset contains the L1 Doppler of healthy GPS/Galileo/QZSS satellites above a 15° elevation mask. Each record carries only what the factor needs: the satellite ECEF position and velocity, satellite clock drift, carrier wavelength, elevation and epoch interval. The original RINEX observation and broadcast-navigation files are reduced from about 65 MB to a small compressed NumPy file.

with np.load(DATA_FILE) as archive:
    data = {name: archive[name] for name in archive.files}
num_epochs = len(data["epoch_dt"])
num_observations = len(data["doppler"])
mean_satellites = num_observations / (num_epochs - 1)
print(f"{num_epochs} epochs at 1 Hz, "
      f"{mean_satellites:.1f} satellites per factor epoch")
print(f"surveyed marker (true velocity = 0): {data['receiver_ecef']}")

301 epochs at 1 Hz, 18.0 satellites per factor epoch
surveyed marker (true velocity = 0): [-3962108.7007  3381309.5532  3668678.6648]

The factor graph¶

One velocity node and one clock-bias node per epoch. Each Doppler observation becomes a DopplerFactor on [v_k, b_{k-1}, b_k], weighted by $1/\sin(\text{elevation})$ . Only differences of the clock biases are observable, so the chain has one gauge freedom, pinned with a prior on the first bias. The graph is linear in all its states, so a single Levenberg-Marquardt iteration solves it.

V = lambda k: symbol('v', k)   # receiver ECEF velocity at epoch k
B = lambda k: symbol('b', k)   # receiver clock bias at epoch k
SIGMA = 0.05                   # zenith range-rate sigma [m/s]

graph, initial = gtsam.NonlinearFactorGraph(), gtsam.Values()
initial.insert(B(0), 0.0)
graph.addPriorDouble(B(0), 0.0, gtsam.noiseModel.Isotropic.Sigma(1, 1e-9))

factors = []
receiver = gtsam.Point3(*data["receiver_ecef"])
for k in range(1, num_epochs):
    initial.insert(V(k), gtsam.Point3(0, 0, 0))
    initial.insert(B(k), 0.0)
    begin, end = data["epoch_offsets"][k - 1:k + 1]
    for i in range(begin, end):
        weight = 1.0 / max(np.sin(data["elevation"][i]), 0.1)
        factor = gtsam.DopplerFactor(
            V(k), B(k - 1), B(k), data["doppler"][i],
            data["wavelength"][i], gtsam.Point3(*data["satellite_position"][i]),
            gtsam.Point3(*data["satellite_velocity"][i]), receiver,
            data["epoch_dt"][k], data["satellite_clock_drift"][i],
            gtsam.noiseModel.Isotropic.Sigma(1, SIGMA * weight))
        graph.add(factor)
        factors.append((k, factor))

result = gtsam.LevenbergMarquardtOptimizer(graph, initial).optimize()
print(f"{len(factors)} Doppler factors, {result.size()} states")

5400 Doppler factors, 601 states

How fast is a station that does not move?¶

The antenna is static, so the estimated velocity is the error. We express it in the local East/North/Up frame of the marker.

velocity = np.array([result.atPoint3(V(k)) for k in range(1, num_epochs)])
enu = velocity @ data["ecef_R_enu"].T

print("velocity error (truth = 0):")
print("  bias  E/N/U = %+.4f %+.4f %+.4f m/s" % tuple(enu.mean(axis=0)))
print("  RMS   E/N/U =  %.4f  %.4f  %.4f m/s" % tuple(np.sqrt((enu ** 2).mean(axis=0))))
print("  horizontal RMS %.4f m/s, 3D RMS %.4f m/s, worst epoch %.4f m/s"
      % (np.sqrt((enu[:, :2] ** 2).sum(axis=1).mean()),
         np.sqrt((enu ** 2).sum(axis=1).mean()),
         np.linalg.norm(enu, axis=1).max()))

velocity error (truth = 0):
  bias  E/N/U = +0.0011 +0.0034 +0.0040 m/s
  RMS   E/N/U =  0.0119  0.0130  0.0319 m/s
  horizontal RMS 0.0176 m/s, 3D RMS 0.0364 m/s, worst epoch 0.1004 m/s

fig = go.Figure()
for i, (name, color) in enumerate(zip(("East", "North", "Up"),
                                      ("#1f77b4", "#2ca02c", "#d62728"))):
    fig.add_scatter(y=enu[:, i], name=name, line=dict(width=1, color=color))
fig.update_layout(
    title="Velocity error of a static antenna, estimated from Doppler alone",
    xaxis_title="epoch [s]", yaxis_title="velocity error [m/s]",
    height=380, margin=dict(l=60, r=20, t=50, b=50))
fig.show()

A few centimetres per second, with no bias to speak of: the range-rate model -- line of sight, satellite velocity, satellite clock drift and the Earth-rotation (Sagnac) rate -- is consistent with the measurements at the level the receiver can measure them.

The other state the graph estimates is the receiver clock, whose drift is the difference of adjacent bias states.

CLIGHT = 299_792_458.0
drift = CLIGHT * np.diff([result.atDouble(B(k))
                          for k in range(num_epochs)])
print("receiver clock drift: mean %+.4f m/s, std %.4f m/s"
      % (drift.mean(), drift.std()))

residual = np.array([f.evaluateError(result.atPoint3(V(k)),
                                     result.atDouble(B(k - 1)),
                                     result.atDouble(B(k)))[0]
                     for k, f in factors])
print("range-rate residual: RMS %.4f m/s over %d observations"
      % (np.sqrt((residual ** 2).mean()), len(residual)))

receiver clock drift: mean +0.0038 m/s, std 0.0685 m/s
range-rate residual: RMS 0.0290 m/s over 5400 observations

Conclusion¶

Doppler alone, with no position states and no atmosphere model, pins the velocity of a static geodetic antenna to a few centimetres per second and its clock drift to a few tenths of a nanosecond per second. In a real navigation graph these factors sit next to PseudorangeFactor or the double-difference factors, sharing the same clock-bias chain, and give the velocity observability that code measurements alone provide only through differencing.

Sources¶

DopplerFactor.h / DopplerFactor.cpp
DopplerFactor.ipynb -- the factor’s reference documentation
Carrier-Phase GNSS: RTK and PPP-RTK -- same receiver and dataset, carrier-phase positioning
RINEX front-end and prepared-data writer: gnss_frontend.py; open-sky source dataset

Preparing another RINEX dataset¶

The observation file must contain L1 code and Doppler measurements, and the broadcast-navigation file must cover the same epochs. The current helper accepts GPS, Galileo and QZSS from a static receiver. First install the front end:

python -m pip install "git+https://github.com/inuex35/cssrlib-numba.git@gtsam-gnss-frontend"

Then change to python/gtsam/examples in a GTSAM checkout and point the existing support script at the two files:

import os
import tempfile
from pathlib import Path

# Give cssrlib-numba's JIT compiler a writable cache location.
os.environ.setdefault("NUMBA_CACHE_DIR", str(Path(tempfile.gettempdir()) / "numba-cache"))
from gnss_frontend import load_doppler, save_doppler_data

raw = load_doppler(
    "path/to/receiver.obs",
    "path/to/broadcast.nav",
    reference_ecef=[-3962108.7007, 3381309.5532, 3668678.6648],
    n_epochs=301,
)
Path("data").mkdir(exist_ok=True)
save_doppler_data(raw, "data/my_doppler_data.npz")

Replace reference_ecef with the surveyed receiver position in metres. It may be omitted to use the RINEX APPROX POSITION XYZ header value instead.