Trilateration, clock bias, and the geometry of space-based positioning
2026-02-26
On 2 May 2000, the United States government switched off Selective Availability — the deliberate degradation of GPS accuracy that had been limiting civilian receivers to about 100 metres. Overnight, accuracy improved to roughly 20 metres. Within two years, GPS receivers were in cars, boats, and hiking packs across the world, and within fifteen years they were in every smartphone. The underlying technology — announced by the Department of Defense in the 1970s, operationally complete by 1995 — is one of the most consequential engineering achievements of the twentieth century, and it rests on a piece of mathematics you can derive in an afternoon.
The principle is deceptively simple: if you know your distance from three points in space whose positions are known, you can compute your position. GPS executes this in three dimensions using satellites instead of fixed ground beacons, measuring distances from the travel time of radio signals at the speed of light. The complication — and the mathematical interest — is that the receiver’s clock is not synchronised with the satellites’ atomic clocks. This introduces an unknown, the clock bias, which requires a fourth satellite to solve. The resulting system of four non-linear equations in four unknowns (three spatial coordinates, one clock bias) is solved iteratively on cheap hardware a dozen times per second. This model derives those equations from first principles and shows why the geometry of the satellite constellation — the Dilution of Precision — matters as much as the signal quality.
How does a handheld GPS receiver know where it is on Earth?
Your phone reports your location to within a few meters. It does this by listening to radio signals from satellites 20,000 km overhead. Each satellite broadcasts: - Its position in space - The exact time the signal was transmitted
Your receiver measures when the signal arrives. From time and the speed of light, you can calculate distance. Multiple satellites give multiple distances, which should intersect at your location.
The mathematical question: How do we solve for position from distance measurements? And why do we need four satellites instead of three?
Radio signals travel at the speed of light: c = 299,792,458 m/s.
If a satellite transmits a signal at time t_{\text{transmit}} and your receiver detects it at time t_{\text{receive}}, the signal traveled for:
\Delta t = t_{\text{receive}} - t_{\text{transmit}}
Distance:
d = c \cdot \Delta t
Example: If \Delta t = 0.07 seconds (70 milliseconds), then:
d = 299,792,458 \times 0.07 = 20,985,472 \text{ m} \approx 21,000 \text{ km}
This places you on a sphere of radius 21,000 km centered at the satellite.
With three satellites, you have three spheres. Your position is where all three intersect.
Triangulation (used in surveying): Measure angles to known points, solve for position using geometry of triangles.
Trilateration (used in GPS): Measure distances to known points, solve for position as the intersection of spheres.
Difference: GPS receivers measure time (which gives distance), not angles.
GPS satellites broadcast their positions in ECEF coordinates (Earth-Centered, Earth-Fixed): - Origin at Earth’s center - Z-axis through North Pole - X-axis through (0° lat, 0° lon) — the intersection of equator and prime meridian - Y-axis completes the right-handed system (through 0° lat, 90° E lon)
Satellite i has position (x_i, y_i, z_i) in ECEF.
Receiver position: (x, y, z) in ECEF.
Distance from receiver to satellite i:
r_i = \sqrt{(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2}
This is the geometric range — the straight-line distance through space.
The pseudorange \rho_i is the measured distance based on signal travel time:
\rho_i = c \cdot \Delta t_i
In a perfect world:
r_i = \rho_i
Each satellite gives one equation. Three satellites → three equations → solve for three unknowns (x, y, z).
But the world is not perfect.
Your receiver’s clock is not atomic. It’s a cheap quartz oscillator that drifts.
Let b = receiver clock error (seconds). If your clock is 0.001 seconds fast, you’ll measure arrival times 0.001 seconds too early.
Effect on pseudorange:
\rho_i = c(\Delta t_i + b)
The measured travel time includes both the true travel time and your clock error.
This adds an unknown — the clock bias b — to the system.
Four unknowns: (x, y, z, b)
Four equations required: Need signals from four satellites minimum.
A 1-microsecond clock error causes a 300-meter position error:
\text{Error} = c \cdot b = 299,792,458 \times 0.000001 = 300 \text{ m}
GPS timing is synchronized to nanoseconds (billionths of a second). A 1-nanosecond error = 30 cm.
Satellite clocks are atomic (rubidium or cesium) and synchronized to GPS time. Your phone’s clock is not. So b must be solved alongside position.
For each satellite i, we have:
\sqrt{(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2} + cb = \rho_i
Where: - (x, y, z) is the receiver position (unknown) - b is the clock bias in seconds (unknown) - (x_i, y_i, z_i) is satellite i’s position (known) - \rho_i is the measured pseudorange (observed) - c is the speed of light (known constant)
With m satellites (typically 6–12 visible), we have m equations and 4 unknowns. This is an overdetermined system — more equations than unknowns. We solve in a least-squares sense.
The square root makes this a nonlinear system. You can’t just write it as a matrix equation \mathbf{Ax} = \mathbf{b} and invert.
Solution method: Linearize around an initial guess, then iterate.
Start with a guess: (x_0, y_0, z_0, b_0)
A crude guess might be: - (x_0, y_0, z_0) = center of Earth (0, 0, 0) or the centroid of satellite positions - b_0 = 0 (assume no clock bias initially)
At the current guess, compute:
Predicted geometric range to satellite i:
r_{i,0} = \sqrt{(x_0 - x_i)^2 + (y_0 - y_i)^2 + (z_0 - z_i)^2}
Predicted pseudorange (including clock bias):
\hat{\rho}_i = r_{i,0} + c b_0
Residual (difference between observed and predicted):
f_i = \hat{\rho}_i - \rho_i = (r_{i,0} + c b_0) - \rho_i
If the guess is perfect, all residuals are zero. Otherwise, we need to adjust the guess.
Approximate each equation with a linear function near the guess.
The Jacobian matrix \mathbf{H} has one row per satellite:
\mathbf{H}_i = \left[ \frac{\partial r_i}{\partial x}\bigg|_0, \; \frac{\partial r_i}{\partial y}\bigg|_0, \; \frac{\partial r_i}{\partial z}\bigg|_0, \; c \right]
Partial derivatives:
\frac{\partial r_i}{\partial x} = \frac{x_0 - x_i}{r_{i,0}}, \quad \frac{\partial r_i}{\partial y} = \frac{y_0 - y_i}{r_{i,0}}, \quad \frac{\partial r_i}{\partial z} = \frac{z_0 - z_i}{r_{i,0}}
Interpretation: These are the direction cosines — the components of the unit vector pointing from the satellite to the receiver.
The last column is c because \frac{\partial(cb)}{\partial b} = c.
We want to find \Delta \mathbf{x} = (\Delta x, \Delta y, \Delta z, \Delta b) such that:
\mathbf{H} \Delta \mathbf{x} \approx -\mathbf{f}
This is a linear least-squares problem. If \mathbf{H} has more rows than columns (more satellites than unknowns), solve:
\Delta \mathbf{x} = -(\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T \mathbf{f}
x_1 = x_0 + \Delta x y_1 = y_0 + \Delta y z_1 = z_0 + \Delta z b_1 = b_0 + \Delta b
Repeat steps 2–5 until |\Delta \mathbf{x}| is smaller than some threshold (e.g., 1 cm for position, 1 ns for time).
Typically converges in 3–5 iterations.
A full 3D GPS calculation is tedious by hand. Here’s a 2D analogue to show the mechanics.
Assume signal speed c = 1 (toy units).
Three transmitters at known positions: - S_1 = (0, 0) - S_2 = (10, 0) - S_3 = (0, 10)
True receiver position: (x, y) = (2, 3) with clock bias b = 0.5.
Measured pseudoranges:
Calculate true ranges: - r_1 = \sqrt{2^2 + 3^2} = \sqrt{13} = 3.606 - r_2 = \sqrt{(2-10)^2 + 3^2} = \sqrt{73} = 8.544 - r_3 = \sqrt{2^2 + (3-10)^2} = \sqrt{53} = 7.280
Add clock bias (c = 1): - \rho_1 = 3.606 + 0.5 = 4.106 - \rho_2 = 8.544 + 0.5 = 9.044 - \rho_3 = 7.280 + 0.5 = 7.780
Initial guess: (x_0, y_0, b_0) = (1, 1, 0)
Predicted ranges: - r_{1,0} = \sqrt{1^2 + 1^2} = 1.414 - r_{2,0} = \sqrt{(1-10)^2 + 1^2} = \sqrt{82} = 9.055 - r_{3,0} = \sqrt{1^2 + (1-10)^2} = \sqrt{82} = 9.055
Predicted pseudoranges (with b_0 = 0): - \hat{\rho}_1 = 1.414 - \hat{\rho}_2 = 9.055 - \hat{\rho}_3 = 9.055
Residuals: - f_1 = 1.414 - 4.106 = -2.692 - f_2 = 9.055 - 9.044 = 0.011 - f_3 = 9.055 - 7.780 = 1.275
Jacobian rows:
For S_1 = (0, 0): \left[ \frac{1-0}{1.414}, \frac{1-0}{1.414}, 1 \right] = [0.707, 0.707, 1]
For S_2 = (10, 0): \left[ \frac{1-10}{9.055}, \frac{1-0}{9.055}, 1 \right] = [-0.994, 0.110, 1]
For S_3 = (0, 10): \left[ \frac{1-0}{9.055}, \frac{1-10}{9.055}, 1 \right] = [0.110, -0.994, 1]
\mathbf{H} = \begin{bmatrix} 0.707 & 0.707 & 1 \\ -0.994 & 0.110 & 1 \\ 0.110 & -0.994 & 1 \end{bmatrix}, \quad \mathbf{f} = \begin{bmatrix} -2.692 \\ 0.011 \\ 1.275 \end{bmatrix}
Solve \mathbf{H} \Delta \mathbf{x} = -\mathbf{f} (using least squares or matrix inversion for the square case):
\Delta \mathbf{x} \approx \begin{bmatrix} 1.0 \\ 2.0 \\ 0.5 \end{bmatrix}
Updated guess: (x_1, y_1, b_1) = (1, 1, 0) + (1.0, 2.0, 0.5) = (2.0, 3.0, 0.5)
This is the true solution! In practice, you’d verify by computing residuals again (they should be near zero).
Below is an interactive GPS solver.
<label>
Number of satellites:
<input type="range" id="num-sats-slider" min="4" max="12" step="1" value="6">
<span id="num-sats-value">6</span>
</label>
<label>
Measurement noise (meters):
<input type="range" id="noise-slider" min="0" max="20" step="1" value="5">
<span id="noise-value">5</span> m
</label>
<label>
Clock bias (nanoseconds):
<input type="range" id="bias-slider" min="-100" max="100" step="10" value="50">
<span id="bias-value">50</span> ns
</label>
<div class="button-group">
<button id="solve-gps">Solve Position</button>
<button id="randomize-sats">Randomize Satellites</button>
</div><h4>Solution:</h4>
<p><strong>True position:</strong> <span id="true-pos"></span></p>
<p><strong>Estimated position:</strong> <span id="est-pos"></span></p>
<p><strong>Position error:</strong> <span id="pos-error"></span> m</p>
<p><strong>True clock bias:</strong> <span id="true-bias"></span> ns</p>
<p><strong>Estimated clock bias:</strong> <span id="est-bias"></span> ns</p>
<p><strong>Iterations:</strong> <span id="iterations"></span></p>
<p><strong>PDOP:</strong> <span id="pdop-value"></span></p>Try this: - Increase satellites: More satellites reduce error and improve PDOP - Add noise: See how measurement errors affect position accuracy - Change clock bias: Even large biases (100 ns) are solved correctly - Randomize satellites: Different geometries produce different PDOP values - Notice: The estimated position (pink) converges to true position (green) despite clock bias
Key insight: The fourth unknown (clock bias) is why GPS needs four satellites minimum. The system simultaneously solves for position and time.
The receiver doesn’t need an atomic clock because: 1. All satellites share the same time reference (GPS time) 2. The clock bias affects all pseudoranges equally (adds the same error to each) 3. This creates a geometric constraint that can be solved
Analogy: If four people tell you distances to a landmark, but your measuring tape is miscalibrated by the same amount for all measurements, you can still figure out both: - Where the landmark is - How miscalibrated your tape is
Position Dilution of Precision quantifies how satellite geometry affects accuracy.
Good geometry (low PDOP < 3): - Satellites spread out in all directions - Some high in sky, some near horizon - Measured distances intersect at steep angles
Bad geometry (high PDOP > 6): - Satellites clustered in one part of sky - All at similar elevations - Measured distances intersect at shallow angles
Formula: PDOP relates measurement error to position error:
\sigma_{\text{position}} = \text{PDOP} \times \sigma_{\text{range}}
Example: If ranging error is 5 m and PDOP = 2:
\sigma_{\text{position}} = 2 \times 5 = 10 \text{ m}
The 31 satellites in the GPS constellation are arranged so that: - 6 orbital planes at 55° inclination - 4–5 satellites per plane - At least 4 satellites visible from anywhere on Earth (usually 6–12)
This ensures good PDOP globally at all times.
Radio signals slow down in the ionosphere (charged particles) and troposphere (water vapor).
Effect: Signals take longer than expected → pseudoranges are too large → position error.
Mitigation: - Dual-frequency receivers: Use two frequencies (L1 and L2) to estimate ionospheric delay - Atmospheric models: Predict tropospheric delay from temperature/pressure - Differential GPS: Use corrections from a nearby base station
Signals bounce off buildings, terrain, or water before reaching the antenna.
Effect: Reflected signals travel farther → pseudoranges too large → position error (typically 5–10 m in urban areas).
Mitigation: - Antenna design: Reject signals from below (likely reflected) - Signal processing: Correlate to detect direct vs. reflected signals - Carrier phase tracking: More resistant to multipath than code phase
If all visible satellites are in one part of the sky (high PDOP): - Position error is amplified - Vertical accuracy suffers most (satellites overhead improve HDOP but not VDOP)
Mitigation: Wait for better satellite configuration, or use multi-GNSS (GPS + GLONASS + Galileo + BeiDou).
Satellites are moving at 3.9 km/s in weaker gravity than Earth’s surface.
Special relativity: Moving clocks run slower (−7 μs/day)
General relativity: Clocks in weaker gravity run faster (+45 μs/day)
Net effect: Satellite clocks run fast by +38 μs/day
Without correction, this would cause 10 km/day position error!
Solution: Satellite clocks are pre-adjusted to run at 10.22999999543 MHz instead of 10.23 MHz exactly.
Standard GPS uses code phase (timing of pseudorandom noise pattern) → 1–5 m accuracy.
RTK (Real-Time Kinematic) uses carrier phase (actual radio wave oscillations) → 1–2 cm accuracy.
Carrier phase ambiguity:
The carrier phase measurement is:
\Phi_i = \frac{1}{\lambda}(r_i + cb) + N_i
Where: - \lambda is the carrier wavelength (19 cm for GPS L1) - N_i is an integer number of full wavelengths (ambiguity)
Challenge: Solve for both position and integer ambiguities simultaneously.
Method: Use a nearby base station with known position to form differential measurements that eliminate common errors, then resolve integer ambiguities using statistical methods (LAMBDA algorithm).
With more equations than unknowns (m > 4), the system is overdetermined. There’s no exact solution that satisfies all equations perfectly (due to measurement noise).
Least squares finds the solution that minimizes the sum of squared residuals:
\min_{\mathbf{x}} \sum_{i=1}^{m} f_i^2(\mathbf{x})
For linear system \mathbf{Hx} = \mathbf{b}, the least-squares solution is:
\mathbf{x} = (\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T \mathbf{b}
This minimizes \|\mathbf{Hx} - \mathbf{b}\|^2.
Covariance: If measurement errors have variance \sigma^2, then:
\text{Cov}(\mathbf{x}) = \sigma^2 (\mathbf{H}^T \mathbf{H})^{-1}
The diagonal elements are the variances of x, y, z, and b.
Next steps: Further exploration of this subject could include: - Differential GPS and SBAS (Satellite-Based Augmentation Systems) - Kalman filtering for continuous position tracking - Integration with INS (Inertial Navigation Systems) - Multi-GNSS constellation usage