A Hands-On Example in Python
Introduction
Total Electron Content (TEC) — the number of free electrons integrated along a signal path through the ionosphere — is one of the most important quantities in GNSS positioning, space weather monitoring, and radio propagation studies. A single GPS/GNSS receiver only measures Slant TEC (STEC), the TEC along the actual line-of-sight to a satellite, which is contaminated by satellite/receiver hardware biases and depends heavily on the elevation angle of the signal. To recover a physically meaningful Vertical TEC (VTEC) map of the ionosphere, we need to solve an inverse (optimization) problem that jointly estimates:
- the spatial and temporal shape of the ionosphere,
- unknown receiver and satellite hardware biases,
- and even the effective height of the ionospheric “thin shell” itself.
This is exactly the kind of nonlinear least-squares optimization problem that shows up in real-world GNSS ionospheric modeling (as used by IGS analysis centers such as CODE, JPL, and ESA). In this post, we’ll build a compact but realistic synthetic example, solve it with scipy.optimize.least_squares, and visualize the results — including a 3D reconstruction of the ionosphere.
Mathematical Formulation
1. Background VTEC surface
We model the “quiet-time” background ionosphere over a local region as a quadratic surface in local coordinates $x$ (longitude-like) and $y$ (latitude-like):
$$
V_{bg}(x,y) = c_0 + c_1x + c_2y + c_3x^2 + c_4y^2 + c_5xy
$$
2. Diurnal variation
The ionosphere is driven by solar illumination, producing a strong diurnal cycle that peaks in the early afternoon:
$$
D(t;A) = 1 + A\sin!\left(\frac{2\pi(t-6)}{24}\right)
$$
3. Traveling Ionospheric Disturbance (TID)
To make the problem realistic and visually interesting, we add a traveling wave-like electron density enhancement (a simplified TID) that drifts across the region over time:
$$
B(x,y,t;A_b,v) = A_b\exp!\left[-\frac{(x-x_0(t))^2}{2\sigma_x^2}-\frac{(y-y_0(t))^2}{2\sigma_y^2}\right], \qquad x_0(t) = x_{00} + vt
$$
4. Full VTEC model
$$
V(x,y,t) = V_{bg}(x,y),D(t;A) + B(x,y,t;A_b,v)
$$
5. Single-layer mapping function
STEC and VTEC are related through the classic thin-shell mapping function, parameterized by the effective ionospheric shell height $H$:
$$
\text{STEC} = \text{VTEC}\cdot M(el,H), \qquad
M(el,H) = \frac{1}{\sqrt{1-\left(\dfrac{R_E}{R_E+H}\cos el\right)^2}}
$$
6. Observation equation
Each slant TEC observation from receiver $r$ to satellite $s$ also contains unknown hardware biases:
7. The optimization problem
All unknowns — the surface coefficients, diurnal amplitude, TID amplitude/speed, shell height, and all biases — are stacked into a single parameter vector $\theta$, estimated by nonlinear least squares:
$$
\hat{\theta} = \arg\min_{\theta}\sum_{i=1}^{N}\Big[\text{STEC}^{model}_i(\theta) - \text{STEC}^{obs}_i\Big]^2
$$
Because one satellite bias and one background level are not jointly identifiable (a classic rank-deficiency problem in GNSS bias estimation), we fix one satellite’s bias as the reference datum ($b_{s_0}=0$).
Python Implementation (Google Colab)
1 | import numpy as np |
Code Walkthrough
Section 1–2 (simulation setup and true VTEC field). We define the Earth radius, the “true” ionospheric shell height (350 km, a typical F2-layer peak height), and a 10°×10° local region. The true_vtec() function is our ground truth generator: a smooth quadratic background surface, modulated by a diurnal factor that peaks around 14:00 local time, plus a moving Gaussian “bump” that represents a traveling ionospheric disturbance drifting eastward at 0.35°/hour. This ground truth is what our optimizer will try to recover — it is never given to the estimator directly.
Section 3 (mapping function). mapping_function() implements the standard single-layer thin-shell obliquity factor. Low elevation angles produce large mapping-function values (STEC much larger than VTEC), which is why elevation-dependent weighting matters so much in real ionospheric estimation.
Section 4 (synthetic observations). We simulate 8 ground stations tracking 6 satellites over 24 hourly epochs, each observation carrying a random ionospheric pierce-point (IPP) location, a random elevation angle, and small Gaussian noise (σ = 0.3 TECU). Random receiver and satellite hardware biases are injected, with one satellite fixed as the zero-bias reference to keep the system identifiable — mirroring how real GNSS bias estimation handles rank deficiency.
Section 5 (vectorized residual function). This is the computational core of the optimization. scipy.optimize.least_squares calls the residual function dozens of times per run, so it must be fast. The residuals() function uses pure NumPy broadcasting and fancy indexing (b_r[station_id], b_s[sat_id]) to evaluate all 1,152 observations at once, with no Python-level loop. For comparison, residuals_naive() performs an equivalent per-observation loop. In practice the vectorized version is roughly 30–45× faster than the loop-based version — the printed benchmark shows this difference directly, which matters a lot once you scale up to realistic GNSS networks with tens of thousands of observations.
Section 6 (optimization). All 23 unknowns (6 surface coefficients, diurnal amplitude, TID amplitude and speed, shell height, 8 receiver biases, 5 satellite biases) are packed into one vector theta and estimated jointly with the Trust Region Reflective (trf) algorithm, which supports bound constraints — here used to keep the shell height physically plausible (200–600 km) and the diurnal/TID amplitudes sign-consistent. Notice how the initial guess for the TID drift speed (theta0[idx_bv] = 0.3) is chosen close to a physically reasonable value: nonlinear least squares is sensitive to the starting point because the TID term is a non-convex Gaussian bump, and a poor initial guess can trap the optimizer in a local minimum that mostly “ignores” the disturbance. This is a genuine and important lesson from real ionospheric TID inversion — good priors or multi-start strategies matter.
Section 7–8 (VTEC maps and visualization). We evaluate the true and estimated VTEC models on a 60×60 spatial grid at a single snapshot time (14:00, near the diurnal peak, when the TID bump is also easiest to see) and build one consolidated figure with nine panels.
Understanding the Visualization
The final figure combines everything into a single 3×3 grid:
- Top-left / top-middle (3D surfaces): the true and estimated VTEC fields at 14:00 local time, letting you visually compare the recovered ionospheric shape — including the traveling disturbance bump — against ground truth.
- Top-right (error contour): the spatial difference between estimated and true VTEC, showing where the reconstruction is most and least accurate.
- Middle-left (scatter fit): modeled vs. observed slant TEC for every observation; points hugging the red diagonal indicate a good fit.
- Middle-center (residual histogram): the distribution of fit residuals, which should look roughly Gaussian and centered near zero if the model and noise assumptions are consistent.
- Middle-right (diurnal curve): true vs. estimated diurnal variation curves over a full day.
- Bottom-left / bottom-middle (bias bars): recovered receiver and satellite hardware biases compared to their true (simulated) values — a direct check on whether the bias-separation part of the inversion worked.
- Bottom-right (scorecard): a compact text summary of all key numbers — RMSE, optimizer iterations, and true-vs-estimated parameter values — for quick reference.
Generated 1152 synthetic slant-TEC observations from 8 stations and 6 satellites over 24 epochs. Vectorized residual evaluation: 0.480 ms Naive loop residual evaluation: 27.279 ms Speed-up factor: 56.8x Optimization finished in 0.090 s, cost=47.5828, nfev=8 Estimated shell height H = 348.81 km (true: 350.0 km) Estimated diurnal amplitude = 0.651 (true: 0.65) Estimated bump amplitude = 11.992 (true: 12.0) Estimated bump drift speed = 0.349 deg/h (true: 0.35) Final RMSE of slant TEC residuals: 0.287 TECU VTEC map RMSE at t=14.0h : 0.059 TECU
Discussion
With a reasonable initial guess, this nonlinear least-squares approach recovers the shell height, the diurnal amplitude, the traveling-disturbance amplitude and speed, and every receiver/satellite bias to within a fraction of a TECU of their true simulated values, while the final slant-TEC residual RMSE stays close to the injected noise level — a strong indicator that the model is neither over- nor under-fitting the data. The most instructive part of this exercise, though, is what happens when the optimizer is given a poor starting guess for the traveling-disturbance parameters: because that term is a non-convex Gaussian in space and time, least_squares can converge to a local minimum where the disturbance is essentially averaged away into the background surface, even though the overall cost still looks “reasonably small.” This mirrors a real, well-known difficulty in operational ionospheric monitoring — smooth, large-scale background TEC is comparatively easy to estimate from sparse ground networks, but transient, localized structures like traveling ionospheric disturbances require denser spatial sampling and often benefit from better-informed initial guesses, regularization, or multi-start/global optimization strategies to be reliably recovered.
