Covariance Propagation¶
This tutorial demonstrates how to propagate state covariance alongside position and velocity using satkit's satstate class. Covariance propagation is essential for:
- Uncertainty quantification: Understanding how errors in an initial state estimate grow over time
- Conjunction assessment: Determining the probability of collision between two objects requires knowing the uncertainty in both objects' positions
- Sensor tasking: Deciding when and where to point a sensor to re-acquire a tracked object
Under the hood, satkit propagates the 6x6 state transition matrix $\Phi(t, t_0)$ alongside the equations of motion. The covariance at time $t$ is then computed as:
$$P(t) = \Phi(t, t_0) \, P(t_0) \, \Phi(t, t_0)^T$$
where $P(t_0)$ is the initial covariance matrix (Montenbruck & Gill 2000, §8.1; Tapley, Schutz & Born 2004, §4.2).
import satkit as sk
import numpy as np
import math
import matplotlib.pyplot as plt
import scienceplots # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']
Setting Up a Satellite State¶
We create a satellite in a near-circular LEO orbit at roughly 400 km altitude. The position and velocity are specified in the GCRF (Geocentric Celestial Reference Frame) inertial frame.
# Epoch
t0 = sk.time(2024, 6, 15, 12, 0, 0)
# LEO orbit: ~400 km altitude, near-circular
r = 6.778e6 # meters (Earth radius + 400 km)
v = math.sqrt(sk.consts.mu_earth / r) # circular orbit speed
pos = np.array([r, 0.0, 0.0]) # meters, GCRF
vel = np.array([0.0, v, 0.0]) # m/s, GCRF
state = sk.satstate(t0, pos, vel)
print(f"Position: {state.pos} m")
print(f"Velocity: {state.vel} m/s")
print(f"Altitude: {np.linalg.norm(state.pos) - 6.378e6:.0f} m")
Position: [6778000. 0. 0.] m Velocity: [ 0. 7668.6356752 0. ] m/s Altitude: 400000 m
Setting Initial Covariance¶
Uncertainty is most naturally expressed in a satellite-local frame whose axes follow the orbit geometry. set_pos_uncertainty(sigma, frame) and set_vel_uncertainty(sigma, frame) take 1-sigma values along the axes of the given frame — RTN (a.k.a. RSW / RIC), NTW, LVLH, or GCRF — build a covariance that is diagonal in that frame, and rotate it into GCRF for storage. Calling them in sequence builds up the full 6×6 covariance; each call preserves the block it isn't updating.
This tutorial uses RTN (Radial / Transverse / Normal), the frame CCSDS uses for covariance in Orbit Data Messages and Conjunction Data Messages. Its axes are defined exactly from the position and angular momentum:
- $\hat{R} = \mathbf{r} / |\mathbf{r}|$ — radial, outward
- $\hat{N} = \mathbf{h} / |\mathbf{h}|$, with $\mathbf{h} = \mathbf{r} \times \mathbf{v}$ — orbit normal (cross-track)
- $\hat{T} = \hat{N} \times \hat{R}$ — transverse (in-track), in the orbit plane and perpendicular to $\mathbf{r}$
$\hat{T}$ is perpendicular to the position, not parallel to the velocity: the two differ by the flight-path angle $\gamma$, which is zero only on a circular orbit. A frame can be tied exactly to $\mathbf{r}$ or to $\mathbf{v}$, but not both — see RTN vs. NTW on an Eccentric Orbit below for when that matters.
# 1-sigma position uncertainty in RTN: [radial, transverse, normal], meters
sigma_pos_rtn = np.array([10.0, 50.0, 20.0])
state.set_pos_uncertainty(sigma_pos_rtn, frame=sk.frame.RTN)
# 1-sigma velocity uncertainty in RTN, m/s
sigma_vel_rtn = np.array([0.01, 0.05, 0.02])
state.set_vel_uncertainty(sigma_vel_rtn, frame=sk.frame.RTN)
print("Position covariance block (GCRF, stored internally), meters^2:")
print(state.cov[0:3, 0:3])
print("\nFull 6x6 covariance matrix (GCRF):")
print(state.cov)
Position covariance block (GCRF, stored internally), meters^2: [[ 100. 0. 0.] [ 0. 2500. 0.] [ 0. 0. 400.]] Full 6x6 covariance matrix (GCRF): [[1.0e+02 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00] [0.0e+00 2.5e+03 0.0e+00 0.0e+00 0.0e+00 0.0e+00] [0.0e+00 0.0e+00 4.0e+02 0.0e+00 0.0e+00 0.0e+00] [0.0e+00 0.0e+00 0.0e+00 1.0e-04 0.0e+00 0.0e+00] [0.0e+00 0.0e+00 0.0e+00 0.0e+00 2.5e-03 0.0e+00] [0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 4.0e-04]]
Propagating with Covariance¶
When a satstate has a covariance matrix set, calling propagate will propagate both the orbit and the covariance using the state transition matrix. Let's propagate forward by one orbital period and examine how the covariance changes.
# Orbital period
a = np.linalg.norm(state.pos) # semi-major axis (circular orbit)
period_sec = 2.0 * math.pi * math.sqrt(a**3 / sk.consts.mu_earth)
print(f"Orbital period: {period_sec / 60:.1f} minutes")
# Propagate forward by one orbit
state_1orbit = state.propagate(sk.duration.from_seconds(period_sec))
print("\nInitial 1-sigma position uncertainty (GCRF diagonal), meters:")
print(f" {np.sqrt(np.diag(state.cov[0:3, 0:3]))}")
print("\nAfter 1 orbit, 1-sigma position uncertainty (GCRF diagonal), meters:")
print(f" {np.sqrt(np.diag(state_1orbit.cov[0:3, 0:3]))}")
Orbital period: 92.6 minutes Initial 1-sigma position uncertainty (GCRF diagonal), meters: [10. 50. 20.] After 1 orbit, 1-sigma position uncertainty (GCRF diagonal), meters: [ 20.18767281 854.68339281 19.99871386]
Examining Uncertainty Growth Over Time¶
To understand how uncertainty evolves, we propagate the state at regular intervals over 24 hours and extract the total position uncertainty (RSS of the 1-sigma values) at each step.
# Propagate at 5-minute intervals over 24 hours
dt_minutes = 5
num_steps = int(24 * 60 / dt_minutes)
times_hr = np.zeros(num_steps)
sigma_pos_rss = np.zeros(num_steps) # RSS position uncertainty
sigma_x = np.zeros(num_steps)
sigma_y = np.zeros(num_steps)
sigma_z = np.zeros(num_steps)
for i in range(num_steps):
dt = sk.duration.from_minutes(dt_minutes * (i + 1))
s = state.propagate(dt)
times_hr[i] = dt_minutes * (i + 1) / 60.0
pos_cov = s.cov[0:3, 0:3]
sigmas = np.sqrt(np.diag(pos_cov))
sigma_x[i] = sigmas[0]
sigma_y[i] = sigmas[1]
sigma_z[i] = sigmas[2]
sigma_pos_rss[i] = np.sqrt(np.trace(pos_cov))
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
# Total RSS position uncertainty
ax1.plot(times_hr, sigma_pos_rss / 1e3, linewidth=1.5)
ax1.set_ylabel("RSS Position 1-sigma (km)")
ax1.set_title("Position Uncertainty Growth Over 24 Hours")
ax1.grid(True, alpha=0.3)
# Per-component uncertainty (GCRF)
ax2.plot(times_hr, sigma_x / 1e3, label="X (GCRF)", linewidth=1.0)
ax2.plot(times_hr, sigma_y / 1e3, label="Y (GCRF)", linewidth=1.0)
ax2.plot(times_hr, sigma_z / 1e3, label="Z (GCRF)", linewidth=1.0)
ax2.set_xlabel("Time (hours)")
ax2.set_ylabel("Position 1-sigma (km)")
ax2.set_title("Per-Component Position Uncertainty (GCRF)")
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Rotating into RTN¶
sk.frametransform.from_gcrf(frame, pos, vel) returns the 3×3 rotation matrix $C$ taking GCRF vectors into a satellite-local frame at a given state. Its rows are the frame's unit axes expressed in GCRF, so we can check the RTN definition directly: $\hat{r}$ maps to $[1, 0, 0]$, $\hat{h}$ to $[0, 0, 1]$, and $\hat{v}$ to $[\sin\gamma, \cos\gamma, 0]$.
C = sk.frametransform.from_gcrf(sk.frame.RTN, state.pos, state.vel)
r_hat = state.pos / np.linalg.norm(state.pos)
h = np.cross(state.pos, state.vel)
h_hat = h / np.linalg.norm(h)
v_hat = state.vel / np.linalg.norm(state.vel)
print(f"r_hat in RTN (expect [1, 0, 0]): {C @ r_hat}")
print(f"h_hat in RTN (expect [0, 0, 1]): {C @ h_hat}")
print(f"v_hat in RTN (expect [sin γ, cos γ, 0]; γ = 0 here): {C @ v_hat}")
r_hat in RTN (expect [1, 0, 0]): [1. 0. 0.] h_hat in RTN (expect [0, 0, 1]): [0. 0. 1.] v_hat in RTN (expect [sin γ, cos γ, 0]; γ = 0 here): [0. 1. 0.]
Covariance in the RTN Frame¶
satkit stores and propagates the covariance in GCRF. To view it in RTN at any time, rotate the position block with the rotation at that time's state:
$$P_{RTN} = C \, P_{GCRF} \, C^T$$
def pos_cov_in(sat, frame):
"""Rotate the 3x3 position covariance from GCRF into a satellite-local frame."""
C = sk.frametransform.from_gcrf(frame, sat.pos, sat.vel)
return C @ sat.cov[0:3, 0:3] @ C.T
# Propagate and extract RTN uncertainties
sigma_radial = np.zeros(num_steps)
sigma_transverse = np.zeros(num_steps)
sigma_normal = np.zeros(num_steps)
for i in range(num_steps):
dt = sk.duration.from_minutes(dt_minutes * (i + 1))
s = state.propagate(dt)
sigmas = np.sqrt(np.diag(pos_cov_in(s, sk.frame.RTN)))
sigma_radial[i] = sigmas[0]
sigma_transverse[i] = sigmas[1]
sigma_normal[i] = sigmas[2]
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(times_hr, sigma_radial / 1e3, label="Radial (R)", linewidth=1.5)
ax.plot(times_hr, sigma_transverse / 1e3, label="Transverse (T)", linewidth=1.5)
ax.plot(times_hr, sigma_normal / 1e3, label="Normal (N)", linewidth=1.5)
ax.set_xlabel("Time (hours)")
ax.set_ylabel("Position 1-sigma (km)")
ax.set_title("Position Uncertainty in RTN Over 24 Hours")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
The RTN plot shows the characteristic pattern of orbital uncertainty growth: transverse (in-track) uncertainty dominates and grows roughly linearly with time, because a small error in velocity (and hence orbital period) makes the satellite drift ahead of or behind its predicted position. Radial and normal uncertainties stay comparatively bounded.
RTN vs. NTW on an Eccentric Orbit¶
On the circular orbit above, $\gamma = 0$ everywhere and $\hat{T}$ happens to lie along the velocity. On an eccentric orbit it does not, and the choice of frame changes what the sigmas mean.
The NTW frame is tied to the velocity instead of the position:
- $\hat{T} = \mathbf{v} / |\mathbf{v}|$ — tangent, exactly along the velocity
- $\hat{W} = \mathbf{h} / |\mathbf{h}|$ — orbit normal (same as RTN's $\hat{N}$)
- $\hat{N} = \hat{T} \times \hat{W}$ — in the orbit plane, perpendicular to the velocity (outward)
NTW is RTN rotated about the orbit normal by $\gamma$. The dominant error — the timing error from period uncertainty, $\delta\mathbf{r} \approx \mathbf{v}\,\delta t$ — points along $\mathbf{v}$, so in NTW it lands entirely on $\hat{T}$. In RTN it is split between $\hat{T}$ and $\hat{R}$, and the radial sigma picks up a component of roughly $\sigma_{\text{along}} \sin\gamma$. The largest flight-path angle on an orbit of eccentricity $e$ is $\arcsin e$.
Below we propagate an orbit with $e = 0.2$ starting at perigee, where $\gamma = 0$ and the two frames coincide, so both views start from the same physical covariance.
# Eccentric orbit: 500 km perigee altitude, e = 0.2, 51.6 deg inclination
ecc = 0.2
rp = sk.consts.earth_radius + 500e3
a_ecc = rp / (1.0 - ecc)
vp = math.sqrt(sk.consts.mu_earth * (2.0 / rp - 1.0 / a_ecc))
inc = math.radians(51.6)
# Start at perigee, where RTN and NTW coincide
ecc_state = sk.satstate(
t0,
np.array([rp, 0.0, 0.0]),
np.array([0.0, vp * math.cos(inc), vp * math.sin(inc)]),
)
ecc_state.set_pos_uncertainty(np.array([10.0, 50.0, 20.0]), frame=sk.frame.RTN)
ecc_state.set_vel_uncertainty(np.array([0.01, 0.05, 0.02]), frame=sk.frame.RTN)
period_ecc = 2.0 * math.pi * math.sqrt(a_ecc**3 / sk.consts.mu_earth)
print(f"Period: {period_ecc / 60:.1f} min; max flight-path angle: {math.degrees(math.asin(ecc)):.1f} deg")
# Propagate at 2-minute intervals over 12 hours
t_ecc_hr = np.arange(2, 12 * 60 + 1, 2) / 60.0
sig_rtn = np.zeros((len(t_ecc_hr), 3))
sig_ntw = np.zeros((len(t_ecc_hr), 3))
gamma_deg = np.zeros(len(t_ecc_hr))
for i, hr in enumerate(t_ecc_hr):
s = ecc_state.propagate(sk.duration.from_hours(hr))
sig_rtn[i] = np.sqrt(np.diag(pos_cov_in(s, sk.frame.RTN)))
sig_ntw[i] = np.sqrt(np.diag(pos_cov_in(s, sk.frame.NTW)))
# Flight-path angle: sin(gamma) = r_hat . v_hat
gamma_deg[i] = math.degrees(math.asin(
np.dot(s.pos, s.vel) / (np.linalg.norm(s.pos) * np.linalg.norm(s.vel))
))
Period: 132.2 min; max flight-path angle: 11.5 deg
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
ax1.plot(t_ecc_hr, sig_rtn[:, 1] / 1e3, label="RTN T ($\\perp$ position)", linewidth=1.5)
ax1.plot(t_ecc_hr, sig_ntw[:, 1] / 1e3, "--", label="NTW T ($\\parallel$ velocity)", linewidth=1.5)
ax1.set_ylabel("1-sigma (km)")
ax1.set_title("Along-track")
ax1.legend()
ax1.grid(True, alpha=0.3)
ax2.plot(t_ecc_hr, sig_rtn[:, 0] / 1e3, label="RTN R ($\\parallel$ position)", linewidth=1.5)
ax2.plot(t_ecc_hr, sig_ntw[:, 0] / 1e3, label="NTW N ($\\perp$ velocity)", linewidth=1.5)
ax2.set_ylabel("1-sigma (km)")
ax2.set_title("In-plane, across track")
ax2.legend()
ax2.grid(True, alpha=0.3)
ax3.plot(t_ecc_hr, gamma_deg, linewidth=1.5)
ax3.set_xlabel("Time (hours)")
ax3.set_ylabel("γ (deg)")
ax3.set_title("Flight-path angle")
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
The along-track sigmas are nearly identical, since $\cos\gamma \approx 1$ even at the largest $\gamma$. The in-plane cross-track sigmas are not: RTN's radial sigma peaks twice per orbit, where $|\gamma|$ is largest, and the peaks grow with the along-track sigma: after 12 hours, about 8 km along-track × $\sin 11.5°$ ≈ 1.6 km appears as "radial" uncertainty. NTW's $\hat{N}$ holds only the uncertainty that is actually perpendicular to the motion, and stays bounded at a few hundred meters. At perigee and apogee ($\gamma = 0$) the two agree.
Choosing a frame:
- RTN — to exchange covariance with other tools (CCSDS OEM / CDM), or when the radial component is what you care about (e.g. radial miss distance in conjunction assessment). On eccentric orbits, expect the radial sigma to include part of the along-track timing error.
- NTW — to separate the timing error (entirely along $\hat{T}$) from the rest of the uncertainty on an eccentric orbit.
The frame also matters on input. set_pos_uncertainty(sigma, frame) builds a covariance that is diagonal in the given frame, and away from $\gamma = 0$ a covariance diagonal in RTN is not diagonal in NTW. Pass the frame your sigmas were actually specified in.
LVLH (z = nadir, y = −orbit normal) is also accepted. It spans the same axes as RTN with different labels and signs, and is intended for porting GN&C code written in that convention.