Carrier-Phase GNSS: RTK and PPP-RTK

Overview¶

Carrier-phase GNSS reaches mm-cm accuracy once the integer ambiguity of each phase measurement is resolved. This notebook shows the two classic flavours, both as GTSAM factor graphs solved with incremental ISAM2 and fixed with LAMBDA:

Part 1 - RTK (double difference): a nearby base station differences away the clocks and atmosphere. Factors: DoubleDifferencePseudorangeFactor, DoubleDifferenceCarrierPhaseFactor.
Part 2 - PPP-RTK (undifferenced): no base station; the receiver clock, tropospheric ZTD and slant ionosphere become state variables, corrected with QZSS CLAS. Factors: UndifferencedPseudorangeFactor, UndifferencedCarrierPhaseFactor.

All the GNSS plumbing (RINEX / CLAS decoding, satellite states, observation extraction, and the LAMBDA covariance bridge) lives in the helper module gnss_frontend.py so the cells below show only the factor graph.

Authors: Frank Dellaert, et al. (see THANKS for the full author list)

See LICENSE for the license information

# Colab: install GTSAM (with the PPP factors) and the GNSS front-end, and
# download the open-sky QZSS CLAS data used by Part 2.
# (Part 1's rover/base RINEX ships inside the front-end package.)
import os
import urllib.request

try:
    import google.colab  # noqa: F401
    !pip install --quiet numpy gtsam-develop \
        "git+https://github.com/inuex35/cssrlib-numba.git@gtsam-gnss-frontend"
except ImportError:
    pass

DATADIR = os.environ.get("GNSS_DATA", "gnss_data")
_BASE = "https://raw.githubusercontent.com/hirokawa/cssrlib-data/main/data"
_FILES = {
    "doy2025-233/233h_rnx.obs":   f"{_BASE}/doy2025-233/233h_rnx.obs",
    "doy2025-233/233h_rnx.nav":   f"{_BASE}/doy2025-233/233h_rnx.nav",
    "doy2025-233/233h_qzsl6.txt": f"{_BASE}/doy2025-233/233h_qzsl6.txt",
    "clas_grid.def":              f"{_BASE}/clas_grid.def",
    "antex/igs20.atx":            "https://files.igs.org/pub/station/general/igs20.atx",
}
for rel, url in _FILES.items():
    dst = os.path.join(DATADIR, rel)
    if not os.path.exists(dst):
        os.makedirs(os.path.dirname(dst), exist_ok=True)
        print("downloading", rel, "...")
        urllib.request.urlretrieve(url, dst)
print("data ready in", DATADIR)

data ready in /home/ubuntu/cssrlib/cssrlib-data/data

Part 1 - RTK (double difference, base + rover)¶

A static rover and a base at a surveyed location see the same satellites. Differencing between receivers and between satellites cancels the clocks and most of the atmosphere, leaving the rover position and the between-receiver integer ambiguities. gnss_frontend.load_rtk returns the per-epoch double-difference bundles (bundled data, no download).

import numpy as np
import gtsam
from gtsam import symbol
from gtsam.examples import gnss_frontend as gnss

rtk = gnss.load_rtk()
print(f"{len(rtk.frames)} epochs;  reference satellite per system: {rtk.ref_of}")

X = symbol('x', 0)                                  # static rover ECEF position
def AM(sat, f): return symbol('n', sat * 10 + f)   # between-receiver SD ambiguity

59 epochs;  reference satellite per system: {<uGNSS.GPS: 0>: 17, <uGNSS.GAL: 1>: 45}

Build the graph incrementally. For each epoch we add, per double difference, a DoubleDifferencePseudorangeFactor and a DoubleDifferenceCarrierPhaseFactor (each forms the rover-base double difference internally). Only between-satellite differences are observable, so one reference ambiguity per constellation is unobservable -- we gauge-pin it with a prior. ISAM2 uses QR factorization for robust marginals.

params = gtsam.ISAM2Params(); params.setFactorization('QR')
isam = gtsam.ISAM2(params)
seen, pinned = set(), set()

for ei, frame in enumerate(rtk.frames):
    graph, val = gtsam.NonlinearFactorGraph(), gtsam.Values()
    if ei == 0:                                     # one static position + loose prior
        val.insert(X, gtsam.Point3(*rtk.x0))
        graph.add(gtsam.PriorFactorPoint3(
            X, gtsam.Point3(*rtk.x0), gtsam.noiseModel.Isotropic.Sigma(3, 30.0)))
    for o in gnss.dd_observations(frame, rtk):
        if (o.ref, o.f) not in pinned:              # gauge: pin reference SD ambiguity
            val.insert(AM(o.ref, o.f), o.sd_ref)
            graph.addPriorDouble(AM(o.ref, o.f), o.sd_ref,
                                 gtsam.noiseModel.Isotropic.Sigma(1, 0.5))
            pinned.add((o.ref, o.f)); seen.add(AM(o.ref, o.f))
        if AM(o.tgt, o.f) not in seen:
            val.insert(AM(o.tgt, o.f), o.sd_tgt); seen.add(AM(o.tgt, o.f))
        graph.add(gtsam.DoubleDifferencePseudorangeFactor(
            X, *o.pr_args, gtsam.noiseModel.Isotropic.Sigma(1, 0.3 * o.w)))
        graph.add(gtsam.DoubleDifferenceCarrierPhaseFactor(
            X, AM(o.ref, o.f), AM(o.tgt, o.f), *o.cp_args, o.lam,
            gtsam.noiseModel.Isotropic.Sigma(1, 0.01 * o.w)))
    isam.update(graph, val)

res = isam.calculateEstimate()
print("float  2D=%.3f  3D=%.3f m" % gnss.enu_error(rtk.pos_ref, rtk.xyz_ref, res.atPoint3(X)))

float  2D=0.155  3D=0.262 m

Resolve the integers. The ambiguity state N is represented as a continuous double (in cycles) in the GTSAM factor graph, so ISAM2 returns the float ambiguity estimate -- GTSAM never optimizes N as an integer. Integer ambiguity fixing is a separate step, performed here with the LAMBDA method using the float estimate and its covariance. gnss_frontend.resolve_integer_ambiguities writes the float ambiguities and their joint covariance into the front-end nav state and runs the LAMBDA search (resamb_lambda). With a short baseline RTK fixes to the mm-cm level.

last = rtk.frames[-1][2]                            # last epoch's DD bundle
nb, fixed = gnss.resolve_integer_ambiguities(rtk.engine, isam, res, X, AM,
                                             last.sat, last.el, seen, rtk.nf, rtk.syss)
print("float  2D=%.3f  3D=%.3f m" % gnss.enu_error(rtk.pos_ref, rtk.xyz_ref, res.atPoint3(X)))
if fixed is not None:
    e2, e3 = gnss.enu_error(rtk.pos_ref, rtk.xyz_ref, fixed)
    print("FIXED  2D=%.3f  3D=%.3f m  (%d SD ambiguities)" % (e2, e3, nb))

float  2D=0.155  3D=0.262 m
FIXED  2D=0.007  3D=0.015 m  (28 SD ambiguities)

Part 2 - PPP-RTK (undifferenced, single receiver + CLAS)¶

Now drop the base station. Without differencing, the receiver clock, tropospheric zenith wet delay (ZTD) and slant ionosphere no longer cancel, so they become state variables (corrected by QZSS CLAS). load_ppp decodes CLAS and returns the per-epoch undifferenced observations.

ppp = gnss.load_ppp(DATADIR)
print(f"{len(ppp.frames)} epochs;  constellations {ppp.syss}")

X = symbol('x', 0)                          # static receiver ECEF position
def CK(ei, si): return symbol('c', ei * 4 + si)  # per-epoch, per-system clock
def ZT(ei): return symbol('z', ei)               # per-epoch ZTD (random walk)
def IO(s):  return symbol('i', int(s))           # per-satellite slant iono
def AM(s, f): return symbol('n', int(s) * 4 + f) # per-sat/freq ambiguity

def ztd_rw():                               # ZTD random-walk between epochs
    def err(this, values, jac):
        a = values.atDouble(this.keys()[0]); b = values.atDouble(this.keys()[1])
        if jac is not None:
            jac[0] = np.array([[1.0]]); jac[1] = np.array([[-1.0]])
        return np.array([a - b])
    return err

100 epochs;  constellations (<uGNSS.GPS: 0>, <uGNSS.GAL: 1>, <uGNSS.QZS: 2>)

Build the graph. Each satellite/frequency adds an UndifferencedPseudorangeFactor and an UndifferencedCarrierPhaseFactor, wired to the position, the system clock, the ZTD and the slant iono (and, for carrier, the ambiguity). The clock is loosely-priored, the iono gets a tight CLAS prior, and ZTD is linked across epochs by the random-walk factor.

params = gtsam.ISAM2Params(); params.setFactorization('QR')
isam = gtsam.ISAM2(params)
seen_ck, seen_io, seen_am = set(), set(), set()

for ei, fr in enumerate(ppp.frames):
    graph, val = gtsam.NonlinearFactorGraph(), gtsam.Values()
    if ei == 0:
        val.insert(X, gtsam.Point3(*ppp.x0))
        graph.add(gtsam.PriorFactorPoint3(
            X, gtsam.Point3(*ppp.x0), gtsam.noiseModel.Isotropic.Sigma(3, 30.0)))
    val.insert(ZT(ei), 0.0)
    if ei == 0:
        graph.addPriorDouble(ZT(0), 0.0, gtsam.noiseModel.Isotropic.Sigma(1, ppp.ztd_sig))
    else:
        graph.add(gtsam.CustomFactor(gtsam.noiseModel.Isotropic.Sigma(1, 0.003),
                  gtsam.KeyVector([ZT(ei), ZT(ei - 1)]), ztd_rw()))
    for o in gnss.undiff_observations(fr, ppp):
        ck = CK(ei, ppp.syss.index(o.sys))
        if ck not in seen_ck:
            seen_ck.add(ck); val.insert(ck, o.clock_init)
            graph.addPriorDouble(ck, 0.0, gtsam.noiseModel.Isotropic.Sigma(1, 1e5))
        if IO(o.sat) not in seen_io:
            seen_io.add(IO(o.sat)); val.insert(IO(o.sat), 0.0)
            graph.addPriorDouble(IO(o.sat), 0.0,
                                 gtsam.noiseModel.Isotropic.Sigma(1, o.iono_sig))
        graph.add(gtsam.UndifferencedPseudorangeFactor(
            X, ck, ZT(ei), IO(o.sat), *o.pr_args,
            gtsam.noiseModel.Isotropic.Sigma(1, 0.6 * o.s_el)))
        if AM(o.sat, o.f) not in seen_am:
            seen_am.add(AM(o.sat, o.f)); val.insert(AM(o.sat, o.f), o.amb_init)
            graph.addPriorDouble(AM(o.sat, o.f), o.amb_init,
                                 gtsam.noiseModel.Isotropic.Sigma(1, 5.0))
        graph.add(gtsam.UndifferencedCarrierPhaseFactor(
            X, ck, ZT(ei), IO(o.sat), AM(o.sat, o.f), *o.cp_args,
            gtsam.noiseModel.Isotropic.Sigma(1, 0.006 * o.s_el)))
    isam.update(graph, val)

res = isam.calculateEstimate()
print("float  2D=%.3f  3D=%.3f m" % gnss.enu_error(ppp.pos_ref, ppp.xyz_ref, res.atPoint3(X)))

float  2D=0.203  3D=0.217 m

Resolve the integers, exactly as in RTK: the graph gives the float ambiguities and LAMBDA fixes them afterwards. Unlike RTK, the PPP position is co-estimated with the clock/ZTD/iono and only becomes observable after a few epochs, so resolve_integer_ambiguities takes a convergence gate conv_sigma=1.0 m before attempting integer fixing.

fr = ppp.frames[-1]
nb, fixed = gnss.resolve_integer_ambiguities(ppp.engine, isam, res, X, AM, fr.sat, fr.el,
                                             seen_am, ppp.nf, ppp.syss, conv_sigma=1.0)
print("float  2D=%.3f  3D=%.3f m" % gnss.enu_error(ppp.pos_ref, ppp.xyz_ref, res.atPoint3(X)))
if fixed is not None:
    e2, e3 = gnss.enu_error(ppp.pos_ref, ppp.xyz_ref, fixed)
    print("FIXED  2D=%.3f  3D=%.3f m  (%d SD ambiguities)" % (e2, e3, nb))

float  2D=0.203  3D=0.217 m
FIXED  2D=0.022  3D=0.074 m  (22 SD ambiguities)

Sources¶

PseudorangeFactor.h / CarrierPhaseFactor.h
Differential GNSS Positioning
Front-end helper: gnss_frontend.py; open-sky CLAS dataset