M7-T1 — The Position Fix: Pseudoranges, Least Squares and Dilution of Precision

Learning objectives

After this lesson, you should be able to:

• Write the pseudorange equations and linearize them about an estimate.
• Solve for position and clock bias by iterated least squares and recognize the convergence.
• Compute the dilution of precision from the beacon geometry and convert it into a position error.
• Choose a beacon position that improves the worst direction.

The intuition

UAV M2-T3: GNSS receivers and positioning describes how a receiver turns signals into a position, and weighted least squares gives the algebra of the solution. This lesson reduces the problem to its core and uses a local beacon system instead of satellites, so that the geometry can be changed and the consequences seen: Hopper's range has beacons around it, a receiver on the vehicle, and an unsynchronized clock.

Why four unknowns

The receiver measures the arrival time of each beacon's signal and multiplies by the speed of light. If its clock were exact that would be the true distance. It is not, so every measurement carries the same unknown offset, the clock bias expressed in metres. The position (three numbers) and the bias (one) are four unknowns, so four beacons are the minimum, and more give a least-squares solution. With six beacons there are two spare equations, and they will become the redundancy that the integrity lesson of this module uses.

Linearize and iterate

Each pseudorange is a nonlinear function of the position, the distance to the beacon, plus the bias. Linearized about a guess, it becomes a linear equation in the correction: the row of the matrix is the unit vector from the beacon to the guess, plus a one for the clock. Solve for the correction, move the guess, repeat. Starting from the origin with a clock bias of 120 m the iteration converges in six steps: the corrections are 143, 79, 38, 2.1 m, then 0.006 m, then 10⁻⁷ m, the signature of Newton convergence.

Geometry decides the accuracy

The measurement error does not map one-to-one onto the position error. The map is the matrix $(\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}$, whose diagonal gives the dilution of precision. When the lines of sight are spread around the sky, the matrix is well conditioned. When they are nearly parallel, as the ground beacons around Hopper are when seen from 60 m up, the height is hardly constrained: moving up by a metre changes every range by a few millimetres. Hopper's six beacons give a horizontal DOP of 0.90 and a vertical DOP of 59: with 1 m of range noise, 0.9 m horizontally and 59 m vertically.

A fix is a linearized least-squares solution, and its accuracy is the measurement noise times the dilution of precision. Beacons on the ground cannot see the height: add one above.

Definitions

Pseudorange $\rho_i$

The measured range to beacon i including the receiver clock error: $\rho_i = \|\mathbf{p} - \mathbf{b}_i\| + c\,\Delta t + \nu_i$.

Clock bias $c\,\Delta t$

The receiver clock error times the speed of light, in metres, common to all pseudoranges.

Geometry matrix $\mathbf{H}$

Rows $[\mathbf{u}_i^{\mathsf{T}},\ 1]$ with $\mathbf{u}_i$ the unit vector from beacon i to the estimated position.

Dilution of precision

The factor between the range noise and the position error: $\mathrm{HDOP}$, $\mathrm{VDOP}$, $\mathrm{TDOP}$, $\mathrm{GDOP}$.

Derivation: linearization and the covariance of the fix

From the nonlinear measurement equation to the normal equations and the DOP.

1. $$\rho_i = \|\mathbf{p} - \mathbf{b}_i\| + \beta + \nu_i, \qquad \beta = c\,\Delta t$$

   The pseudorange: true distance, common clock bias, and measurement noise of variance $\sigma_\rho^2$.

2. $$\|\mathbf{p} - \mathbf{b}_i\| \approx \|\hat{\mathbf{p}} - \mathbf{b}_i\| + \mathbf{u}_i^{\mathsf{T}}(\mathbf{p} - \hat{\mathbf{p}}), \quad \mathbf{u}_i = \dfrac{\hat{\mathbf{p}} - \mathbf{b}_i}{\|\hat{\mathbf{p}} - \mathbf{b}_i\|}$$

   First-order Taylor expansion of the distance about the estimate; the gradient of a norm is the unit vector.

3. $$\rho_i - \|\hat{\mathbf{p}} - \mathbf{b}_i\| - \hat\beta = \mathbf{u}_i^{\mathsf{T}}\Delta\mathbf{p} + \Delta\beta + \nu_i$$

   Move the known terms to the left. The right side is linear in the corrections to position and bias, with the row $[\mathbf{u}_i^{\mathsf{T}},\ 1]$.

4. $$\Delta\boldsymbol{\rho} = \mathbf{H}\,\Delta\mathbf{x} + \boldsymbol{\nu} \;\Rightarrow\; \Delta\hat{\mathbf{x}} = (\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}\mathbf{H}^{\mathsf{T}}\Delta\boldsymbol{\rho}$$

   Stack the equations for all beacons and solve in the least-squares sense; the minimum of $\|\Delta\boldsymbol{\rho} - \mathbf{H}\Delta\mathbf{x}\|^2$.

The covariance

1. $$\mathrm{Cov}(\Delta\hat{\mathbf{x}}) = (\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}\mathbf{H}^{\mathsf{T}}\,(\sigma_\rho^2\mathbf{I})\,\mathbf{H}(\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}$$

   Propagate independent noise of variance $\sigma_\rho^2$ through the linear map.

2. $$\mathrm{Cov} = \sigma_\rho^2\,(\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1} = \sigma_\rho^2\,\mathbf{Q}$$

   The factors collapse. The diagonal of Q, square-rooted, are the dilutions of precision.

3. $$Q_{33}\ \text{large} \iff \text{the third column of } \mathbf{H} \text{ is nearly in the span of the others}$$

   The vertical component of the line of sight, the third column, is nearly the same for all ground beacons seen from 60 m, so it is nearly proportional to the clock column of ones: vertical position and clock bias become hard to tell apart.

The accuracy of a fix is the range noise multiplied by the geometry factor $(\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}$. When the unit vectors have little spread in a direction, the matrix is nearly singular in that direction and the error is large there, however good the measurements.

A worked example

Six beacons around Hopper's range (all within 8 m of the ground, spread over ±430 m), Hopper at (40, 20, 60) m, clock bias 120 m, range noise 1 m.

1. Gauss-Newton. From the origin with a zero clock bias the corrections are 143, 78.9, 38.4, 2.14, 0.0064, 9×10⁻⁸ and 4×10⁻¹² m: quadratic convergence in six steps to (40.0000, 20.0000, 60.0000) with a clock bias of 120.0000 m.
2. DOP. HDOP = 0.90, VDOP = 59.06, TDOP = 8.11, GDOP = 59.62. With a range noise of 1 m the 1σ errors are 0.90 m horizontal and 59 m vertical.
3. Check by Monte Carlo. 2000 noisy fixes give a horizontal rms of 0.94 m and a vertical rms of 62.9 m, against the predicted 0.90 and 59.1 m (the difference is from the sample size and the linearization).
4. A seventh beacon. On the ground at the east edge: VDOP 32.2. On a 40 m tower: 8.8. On a 150 m tethered balloon: 1.25, and GDOP 1.55. The horizontal DOP barely changes (0.88, 0.82, 0.83).
5. Which four. With only four beacons, (0, 1, 2, 3) gives VDOP 63 and (0, 2, 4, 5) gives VDOP 90: the choice of subset matters.
6. Reading. The height is the problem. Hopper must land to within 3 m, so a vertical error of 59 m from the fix is useless, and the altitude has to come from the radar altimeter or from the balloon beacon.

import math
import random
import gnckit as g

# six range beacons around Hopper's range (metres, local east-north-up), all on or near the ground
beacons = [(-300.0, -200.0, 2.0), (350.0, -250.0, 5.0), (400.0, 180.0, 0.0), (-250.0, 300.0, 8.0), (50.0, 420.0, 3.0), (60.0, -430.0, 1.0)]
truth = (40.0, 20.0, 60.0)               # Hopper at 60 m altitude
cb_true = 120.0                           # receiver clock bias times the speed of light, metres
sigma = 1.0                               # pseudorange noise, metres

def ranges(p, cb, nz=None):
    out = []
    for i, b in enumerate(beacons):
        r = math.dist(p, b) + cb
        if nz: r += nz[i]
        out.append(r)
    return out

def solve(meas, idx, x0=(0.0, 0.0, 0.0, 0.0), iters=8):
    x = list(x0)
    hist = []
    for _ in range(iters):
        H, res = [], []
        for i in idx:
            b = beacons[i]
            r = math.dist(x[:3], b)
            H.append([(x[0] - b[0]) / r, (x[1] - b[1]) / r, (x[2] - b[2]) / r, 1.0])
            res.append(meas[i] - (r + x[3]))
        HtH = g.mat_mul(g.mat_T(H), H)
        Htr = [sum(H[k][j] * res[k] for k in range(len(idx))) for j in range(4)]
        dx = g.mat_solve(HtH, Htr)
        x = [x[j] + dx[j] for j in range(4)]
        hist.append(math.sqrt(sum(d * d for d in dx)))
    return x, hist

meas = ranges(truth, cb_true)
x, hist = solve(meas, range(6))
print('noise-free fix from the origin with a clock bias of 120 m, size of the correction at each iteration (m):')
print('  ' + ', '.join(f'{h:.3g}' for h in hist))
print(f'  solution ({x[0]:.4f}, {x[1]:.4f}, {x[2]:.4f}), clock bias {x[3]:.4f} m  (truth {truth}, {cb_true})')

def dop(idx, at):
    H = []
    for i in idx:
        b = beacons[i]
        r = math.dist(at, b)
        H.append([(at[0] - b[0]) / r, (at[1] - b[1]) / r, (at[2] - b[2]) / r, 1.0])
    Q = g.mat_inv(g.mat_mul(g.mat_T(H), H))
    return math.sqrt(Q[0][0] + Q[1][1]), math.sqrt(Q[2][2]), math.sqrt(Q[3][3]), math.sqrt(sum(Q[i][i] for i in range(4)))
h_, v_, t_, gd = dop(range(6), truth)
print(f'six beacons: HDOP {h_:.2f}, VDOP {v_:.2f}, TDOP {t_:.2f}, GDOP {gd:.2f}')
print(f'  with sigma = {sigma} m: horizontal error {h_ * sigma:.2f} m, vertical error {v_ * sigma:.2f} m (1 sigma)')

random.seed(3)
errs = []
for _ in range(2000):
    nz = [random.gauss(0, sigma) for _ in range(6)]
    xx, _h = solve(ranges(truth, cb_true, nz), range(6), x0=x)
    errs.append((xx[0] - truth[0], xx[1] - truth[1], xx[2] - truth[2]))
eh = math.sqrt(sum(e[0] ** 2 + e[1] ** 2 for e in errs) / len(errs))
ev = math.sqrt(sum(e[2] ** 2 for e in errs) / len(errs))
print(f'Monte Carlo, 2000 fixes: horizontal rms {eh:.2f} m, vertical rms {ev:.2f} m (the DOP predicts {h_ * sigma:.2f} and {v_ * sigma:.2f})')

print('adding a seventh beacon to improve the vertical accuracy:')
for name, extra in (('on the ground, east edge', (380.0, 10.0, 0.0)), ('on a 40 m tower', (-30.0, 10.0, 40.0)), ('on a 150 m tethered balloon', (-30.0, 10.0, 150.0))):
    beacons.append(extra)
    hh, vv, tt, gg = dop(range(7), truth)
    print(f'  {name:30s} HDOP {hh:.2f}, VDOP {vv:.2f}, GDOP {gg:.2f}')
    beacons.pop()
print('four beacons only (the minimum): ')
for idx in ((0, 1, 2, 3), (0, 2, 4, 5)):
    hh, vv, tt, gg = dop(idx, truth)
    print(f'  beacons {idx}: HDOP {hh:.2f}, VDOP {vv:.2f}, GDOP {gg:.2f}')

noise-free fix from the origin with a clock bias of 120 m, size of the correction at each iteration (m):
  143, 78.9, 38.4, 2.14, 0.00636, 9.12e-08, 4.1e-12, 1.66e-12
  solution (40.0000, 20.0000, 60.0000), clock bias 120.0000 m  (truth (40.0, 20.0, 60.0), 120.0)
six beacons: HDOP 0.90, VDOP 59.06, TDOP 8.11, GDOP 59.62
  with sigma = 1.0 m: horizontal error 0.90 m, vertical error 59.06 m (1 sigma)
Monte Carlo, 2000 fixes: horizontal rms 0.94 m, vertical rms 62.86 m (the DOP predicts 0.90 and 59.06)
adding a seventh beacon to improve the vertical accuracy:
  on the ground, east edge       HDOP 0.88, VDOP 32.24, GDOP 32.58
  on a 40 m tower                HDOP 0.82, VDOP 8.81, GDOP 8.96
  on a 150 m tethered balloon    HDOP 0.83, VDOP 1.25, GDOP 1.55
four beacons only (the minimum): 
  beacons (0, 1, 2, 3): HDOP 1.10, VDOP 62.97, GDOP 63.58
  beacons (0, 2, 4, 5): HDOP 1.26, VDOP 90.22, GDOP 91.12

In robotics

The same equations run in every GNSS receiver, and in every ranging system.

• GNSS. The beacons are satellites moving on known orbits, the clock bias is that of a receiver, and the DOP of a good constellation is about 1 to 2 horizontally and 1.5 to 3 vertically: the vertical is always worse because all the satellites are above the receiver.
• Ranging beacons. Ultra-wideband anchors in a warehouse, radio beacons at a test range or acoustic transponders underwater all give pseudoranges, and the layout of the beacons is a design variable of the system, chosen from the DOP of the volume to be covered.
• Layout and the height. The lesson of Hopper's range is general: place beacons so that the lines of sight spread in every direction, including above. A tower, a tethered balloon or a beacon on a hill buys more than adding ground beacons.
• Weights and outliers. Real receivers weight each pseudorange by its expected variance (elevation, signal-to-noise ratio), and check the residuals, which is the integrity lesson of this module.
• The fix as a measurement. In an aided navigation system the fix, or the pseudoranges themselves, are the measurements of the filter in the next lesson; DOP is the covariance the filter uses for them.

The lab M7-L1: Solve the Fix solves this fix from six pseudoranges and asks which beacon to add. The extra equations of a six-beacon fix are used for fault detection in M7-T6: Navigation Integrity: Fault Detection and Protection Levels. Satellite signals and error sources are covered at UAV M2-T1: GNSS fundamentals.

Common misconceptions

Key takeaways

1. A pseudorange is the true range plus a common clock bias; position and bias are four unknowns, so four beacons are the minimum.
2. Linearize about an estimate: the row of H is the unit vector from the beacon plus a one; iterate the least-squares correction.
3. Convergence from 143 m of initial error took six iterations, quadratically.
4. The covariance of the fix is $\sigma_\rho^2(\mathbf{H}^{\mathsf{T}}\mathbf{H})^{-1}$; its diagonal gives HDOP, VDOP and TDOP.
5. Hopper's six ground beacons give HDOP 0.90 and VDOP 59: 0.9 m horizontal and 59 m vertical for 1 m of range noise.
6. A beacon above the vehicle (a 150 m balloon) brings VDOP to 1.25; a 40 m tower to 8.8.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• P. D. Groves, Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, 2nd edition, Artech House, 2013, chapter 9 on GNSS positioning and DOP. Paid.
• E. D. Kaplan and C. J. Hegarty (eds.), Understanding GPS/GNSS: Principles and Applications, 3rd edition, Artech House, 2017. Paid.
• ESA Navipedia, the GNSS encyclopaedia (use its search for "DOP" and "least squares") , free pages on positioning and dilution of precision. Free to read; the licence is stated on the site.
• MIT 16.322 Stochastic Estimation and Control (OpenCourseWare) , for weighted least squares and its covariance. Free; licence CC BY-NC-SA 4.0.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Lab
  M7-L1: Solve the Fix

  From six pseudoranges to range beacons, solve Hopper's position and clock bias, compute the dilution of precision, and choose which beacon to add to improve the vertical accuracy.

Other courses on this site treat neighbouring topics from their own angle:

• UAV M2-T1: GNSS fundamentals
• UAV M2-T3: GNSS receivers and positioning
• UAV M2-T4: GNSS errors and multipath
• UAV M2-T6: GNSS jamming and spoofing
• UAV M3-T6: Error-State Kalman Filter
• UAV M3-T5: Extended Kalman Filter
• UAV M3-T11: Multi-sensor fusion architectures
• UAV M3-T13: Innovation monitoring
• UAV M3-T16: Sensor fault detection and isolation
• UAV M6-T5: PX4 state estimation and EKF2
• UAV M2-T11: Visual-Inertial Odometry
• SLAM M8-T5: GPS-Denied Localization
• Linear Algebra M5-T6: Weighted Least Squares
• Linear Algebra M6-T5: Singular Value Decomposition
• Probability for Robotics (modules 7, 8 and 12: Kalman, extended and error-state filters)

This course leaves the following to the pages that already teach them:

• Constellations, signals, error sources, multipath, jamming and spoofing
• The Kalman, extended and error-state filter derivations
• Visual, LiDAR and radar odometry, VIO and SLAM as aiding sources
• Innovation tests (NIS, NEES) as such