From Least Squares to Fault Exclusion
A GNSS receiver never observes its position directly. All it measures are pseudoranges, which are distances to satellites contaminated by clock error and noise. How well we can turn those measurements into a position depends on three things: the estimator, how we treat measurements of different quality, and whether we notice when one measurement is simply wrong.
In this article we build a concrete example and attack it step by step. The scenario has nine satellites and elevation-dependent noise, and one low-elevation satellite is corrupted by a 30 m multipath-like fault. We compare three estimators over 20,000 Monte Carlo trials:
- OLS: ordinary least squares
- WLS: weighted least squares
- WLS + FDE: weighted least squares with fault detection and exclusion
1. Problem Setup
Work in a local East-North-Up (ENU) frame with the true receiver position at the origin. For satellite $i$ at position $\mathbf{s}_i$, the pseudorange is
$$
\rho_i = \lVert \mathbf{s}_i - \mathbf{x} \rVert + b + \varepsilon_i + \delta_i ,
$$
where $\mathbf{x}=(E,N,U)^\top$ is the receiver position, $b = c,\delta t_r$ is the receiver clock bias expressed in meters, $\varepsilon_i$ is random noise, and $\delta_i$ is a fault bias that is nonzero for only one satellite.
The unknown vector is $\mathbf{p} = (E, N, U, b)^\top$. Linearizing around an approximate position gives
$$
\Delta\boldsymbol{\rho} = H,\Delta\mathbf{p} + \boldsymbol{\varepsilon}, \qquad
H = \begin{bmatrix} -\mathbf{u}_1^\top & 1 \ \vdots & \vdots \ -\mathbf{u}_n^\top & 1 \end{bmatrix}, \qquad
\mathbf{u}_i = \frac{\mathbf{s}_i - \mathbf{x}}{\lVert \mathbf{s}_i - \mathbf{x} \rVert}.
$$
Low-elevation signals travel through more atmosphere and are more prone to multipath, so we model the noise as elevation dependent:
$$
\sigma_i^2 = \sigma_a^2 + \frac{\sigma_b^2}{\sin^2\theta_i}, \qquad \sigma_a = 0.3\ \text{m},\quad \sigma_b = 1.5\ \text{m},
$$
where $\theta_i$ is the elevation angle of satellite $i$.
2. Three Estimators
OLS treats every satellite equally:
$$
\Delta\hat{\mathbf{p}}_{\text{OLS}} = (H^\top H)^{-1} H^\top \Delta\boldsymbol{\rho}.
$$
WLS gives noisy satellites less influence through $W=\mathrm{diag}(1/\sigma_1^2,\dots,1/\sigma_n^2)$:
The geometry quality is summarized by the dilution of precision. With $Q=(H^\top H)^{-1}$:
$$
\mathrm{GDOP}=\sqrt{\operatorname{tr}Q},\quad
\mathrm{PDOP}=\sqrt{Q_{11}+Q_{22}+Q_{33}},\quad
\mathrm{HDOP}=\sqrt{Q_{11}+Q_{22}},\quad
\mathrm{VDOP}=\sqrt{Q_{33}}.
$$
WLS + FDE adds a consistency check. With $K=(H^\top W H)^{-1}H^\top W$, the weighted residual sum of squares
$$
T = \lVert W^{1/2}(I - HK),\Delta\boldsymbol{\rho} \rVert^2
$$
follows a $\chi^2_{n-4}$ distribution when all measurements are healthy. If $T$ exceeds the threshold $\chi^2_{n-4,,1-P_{FA}}$, a fault is declared. We then recompute the solution with each satellite left out in turn and keep the subset whose residual statistic $T_{(-j)}$ is smallest:
$$
j^\ast = \arg\min_{j}, T_{(-j)} .
$$
Finally, to visualize how the WLS solution is found, we use the cost function with the vertical and clock components profiled out:
$$
J(E,N)=\min_{U,,b};(\Delta\boldsymbol{\rho}-H\Delta\mathbf{p})^\top W(\Delta\boldsymbol{\rho}-H\Delta\mathbf{p}).
$$
3. Full Source Code
1 | import time |
4. Code Walkthrough
Section 1: Scenario definition
Instead of full orbital mechanics, the receiver sits at the origin of an ENU frame and each satellite is placed 22,000 km away in the direction given by its azimuth and elevation. This keeps the geometry realistic enough to study estimation error without any ephemeris handling.
The nine satellites are deliberately unevenly spread. G01, G02, G05 and G07 are high, while G06 and G08 sit at 15° and 12°, which is typical of a real sky. The array u holds the unit line-of-sight vectors, and sigma implements the noise model $\sigma_i^2=\sigma_a^2+\sigma_b^2/\sin^2\theta_i$. The weights w are simply $1/\sigma_i^2$.
The faulty satellite is G04 (index 3, elevation 22°). Its 30 m bias is about seven times its own standard deviation, which is large but plausible for a strong multipath reflection.
Section 2: Design matrix and DOP
H is built by stacking $-\mathbf{u}_i^\top$ with a column of ones for the clock. The function dop_values inverts $H^\top H$ (or $H^\top W H$ if weights are given) and extracts the DOP components. Because the frame is ENU, HDOP and VDOP fall directly out of the diagonal of $Q$ with no rotation needed.
Section 3: Gauss-Newton
This is the true nonlinear solver. At each iteration it recomputes the ranges $\rho_i$ from the current estimate, rebuilds $H$ from the updated line-of-sight vectors, and solves the weighted normal equations for the step dx. The loop stops when the step is smaller than 0.1 mm.
The initial guess is deliberately poor, roughly 99 km off in position and 300 m off in clock. This shows that the method does not need a good starting point. The single-epoch measurement pr_one contains noise and the 30 m fault, so the converged answer will not be exactly zero error. It will be the WLS answer for that noisy, faulty epoch.
Section 4: Vectorized Monte Carlo
A naive implementation would loop over 20,000 trials and, for the FDE method, solve nine leave-one-out least squares problems per trial, which is 180,000 small linear solves in a Python loop. That takes seconds to minutes depending on the machine.
The speedup rests on one observation. Every estimator here is linear in the measurement residual vector $\Delta\boldsymbol{\rho}$. So the operators can be computed once, outside the trial loop:
K_olsandK_wlsmap residuals to state corrections, so all 20,000 solutions are a single matrix productdy @ K.T.K_loo[j]is the estimator with satellite $j$ removed, stored as a $4\times n$ matrix whose $j$-th column is zero.M_loo[j]is the matching whitened residual operator, with row $j$ zeroed so the excluded satellite does not contribute to the test statistic.np.einsumapplies all nine leave-one-out operators to all 20,000 trials at once, givingsol_loowith shape $(N, n, 4)$ andT_loowith shape $(N, n)$.
The decision logic is also vectorized. detected compares each trial’s full-set statistic with the $\chi^2$ threshold. j_star picks the best satellite to exclude for every trial, and np.where selects either the excluded solution or the plain WLS solution. The printed timing is the time for this whole block.
Linearizing around the true position for the Monte Carlo is legitimate here. With errors of tens of meters and satellites 22,000 km away, the neglected second-order term is about $\lVert\Delta\mathbf{x}\rVert^2/(2R)\approx 30^2/(2\times 2.2\times10^{7})\approx 2\times10^{-5}$ m, which is negligible.
Section 5: Cost surface
profile_cost evaluates $J(E,N)$ on a grid. For each horizontal point $(E,N)$ it subtracts that point’s contribution from the residuals and then solves the remaining two-parameter problem $(U,b)$ in closed form, again for all grid points at once. The result is a smooth bowl over the horizontal plane, and its lowest point is the WLS horizontal estimate. The grid is scaled from the actual WLS error so the truth and the minimum both fit in the picture.
Section 6: Visualization
All six panels go into one figure with constrained_layout, which handles the polar plot, the two 3D axes and the color bar without overlaps. In the 3D cost surface, computed_zorder = False together with explicit zorder values keeps the markers drawn on top of the semi-transparent surface. Only 1,500 of the 20,000 points are drawn in the 3D scatter to keep it readable, while the bar chart and the CDF use all 20,000.
5. Execution Results
=== Satellite geometry and noise model === Sat Az[deg] El[deg] sigma[m] weight G01 20.0 78.0 1.56 0.410 G02 75.0 55.0 1.86 0.290 G03 130.0 35.0 2.63 0.144 G04 185.0 22.0 4.02 0.062 <-- faulty G05 240.0 48.0 2.04 0.240 G06 300.0 15.0 5.80 0.030 G07 340.0 62.0 1.73 0.336 G08 100.0 12.0 7.22 0.019 G09 210.0 30.0 3.01 0.110 === DOP (all satellites, unweighted) === GDOP: 1.872 PDOP: 1.640 HDOP: 0.944 VDOP: 1.341 TDOP: 0.904 === Gauss-Newton iterations (weighted, single epoch) === iter 3D pos err [m] clock err [m] 0 98994.9494 300.0000 1 196.5270 39.6186 2 8.8030 6.6882 3 8.8020 6.6870 4 8.8020 6.6870 === Monte Carlo (20000 trials, vectorized: 0.133 s) === Test threshold (chi-square, dof=5, P_FA=0.001): 20.52 Fault detection rate : 98.95 % Correct exclusion rate : 98.93 % Method RMS 3D RMS H RMS V 95% 3D mean error (E, N, U) [m] OLS 14.14 11.47 8.28 19.34 [ 1.965 10.395 4.459] WLS 10.27 5.12 8.90 16.59 [-0.019 3.923 6.725] WLS + FDE 6.84 3.36 5.95 12.49 [0.003 0.015 0.098] 3D RMS reduction (OLS -> WLS + FDE): 51.7 %
6. Reading the Results
The geometry (panel 1)
The sky plot shows where the satellites sit and how noisy each one is. High-elevation satellites near the center (G01, G07, G02) are light yellow-green, meaning $\sigma$ below 2 m. Satellites near the horizon (G08 at 12° and G06 at 15°) are dark, with $\sigma$ of 7.2 m and 5.8 m. The red circle marks G04, the faulty satellite. It is a low-elevation satellite with a fairly large $\sigma$ of about 4 m, which matters for detection. A fault on a noisy satellite is harder to notice than one on a clean satellite.
The unweighted DOP values printed in the console tell the same story. The horizontal geometry is good, with HDOP below 1, while the vertical is weaker with VDOP around 1.3. That is the usual GNSS pattern, because all satellites are above the receiver and none are below.
Nonlinear convergence (panel 2)
Starting about 99 km from the truth, Gauss-Newton drops the position error to roughly 200 m after one iteration and to under 10 m after two. Iterations three and onward change nothing visible. The clock error follows the same pattern. The curve flattens at a level of a few meters because that is the noise and fault floor of this single epoch, not a convergence failure. This is also why the linearized Monte Carlo in the next step is justified. Once you are within tens of meters, one linear step is essentially the final answer.
The 3D error clouds (panel 3)
Each cloud is the distribution of position errors in East, North and Up, and the black star is the truth. The three clouds are different in both center and spread.
The OLS cloud (red) is displaced from the origin. The fault at G04, which lies to the south at low elevation, pushes the estimate away from the satellite, mostly toward the north. OLS gives G04 the same weight as a satellite at 78°, so it absorbs the full bias. The mean error from the console confirms this, with a North component of about 10 m.
The WLS cloud (orange) is closer to the truth in the horizontal plane because G04’s weight is small, but it is still shifted, particularly in the vertical.
The WLS + FDE cloud (blue) is centered on the star and tighter overall. Once the faulty satellite is excluded, the bias disappears, and the mean error drops to a few centimeters or less on every axis.
The cost surface (panel 4)
This panel shows why a least squares solution exists and where it lands. The surface is a smooth bowl, and the weighted cost $J$ has a single minimum, marked by the red dot. The black star is the true position. They do not coincide. The gap between them is the bias introduced by the fault. The bowl is also elongated, being steeper in one horizontal direction than the other. That elongation is the geometry at work, since directions with well-spread satellites are tightly constrained and directions with poor coverage are shallow, allowing larger error. The filled contour map on the floor is the same information seen from above.
RMS by method (panel 5)
The bar chart condenses the Monte Carlo into numbers. In my run the 3D RMS error fell from about 14 m with OLS to about 10 m with WLS and about 7 m with WLS + FDE, a reduction of roughly half. Two details are worth noting.
First, the horizontal improvement is much larger than the vertical one. Horizontal RMS drops from about 11.5 m to 5 m with weighting alone and to about 3.4 m with FDE. Vertical RMS barely moves with weighting and is even slightly worse for WLS than for OLS. The reason is that low-elevation satellites are the main source of vertical information, so downweighting them, while sensible for noise, weakens the vertical geometry. WLS is not a free lunch, and it trades vertical strength for lower noise.
Second, FDE reduces error in every category. Removing the biased measurement fixes the bias that weighting could only dampen.
The error distribution (panel 6)
The cumulative distribution shows that the improvement is not just in the average. The blue curve (WLS + FDE) sits far to the left of the orange (WLS) and red (OLS) curves across the entire probability range. At the 95 % line, the 3D error falls from roughly 19 m for OLS to roughly 12.5 m for WLS + FDE. The OLS curve also has essentially zero probability of an error below 5 m, which is the fingerprint of a systematic bias, since a zero-mean estimator would have plenty of small errors.
The console reports how the test behaved. With a false alarm probability of $10^{-3}$ and 30 m of fault on a satellite with $\sigma\approx 4$ m, the fault was detected in about 99 % of trials and the right satellite was excluded in nearly all of those. The small remaining fraction is where the noise happened to cancel the fault and the test statistic stayed below the threshold.
7. Conclusion
Starting from the same nine pseudoranges, three levels of processing produced very different accuracy:
- Weighting by elevation-dependent noise removes much of the horizontal error but weakens vertical geometry.
- A residual test with leave-one-out exclusion removes the systematic bias that neither OLS nor WLS can handle, and it did so with a correct identification rate near 99 %.
- Because every estimator is linear in the residuals, precomputing the operators turns a 180,000-solve loop into a handful of matrix products, so the full Monte Carlo runs in a fraction of a second.
The same structure extends naturally to more realistic setups. The elevation model can be replaced by C/N0-based weights, the single fault by several candidates, and the single-epoch solution by a Kalman filter. The core idea stays the same: model the noise honestly, check the measurements for consisten
