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.

In [1]:

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']

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.

In [2]:

# 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")

# 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 where the axes are aligned with the orbit geometry. satkit provides a unified API for this — set_pos_uncertainty(sigma, frame) and set_vel_uncertainty(sigma, frame) — that accepts any of the supported orbital or inertial frames: LVLH, RIC (a.k.a. RSW/RTN), NTW, or GCRF.

Here we use LVLH (Local Vertical, Local Horizontal), where z = nadir, y = anti-angular-momentum, and x ≈ along-track for a circular orbit. Calling set_pos_uncertainty and set_vel_uncertainty in sequence builds up the full 6×6 covariance — each call preserves the block it isn't updating.

In [3]:

# 1-sigma position uncertainty in LVLH (x ≈ along-track, y = cross-track, z = nadir)
sigma_pos_lvlh = np.array([50.0, 20.0, 10.0])
state.set_pos_uncertainty(sigma_pos_lvlh, frame=sk.frame.LVLH)

# 1-sigma velocity uncertainty in LVLH
sigma_vel_lvlh = np.array([0.05, 0.02, 0.01])  # m/s
state.set_vel_uncertainty(sigma_vel_lvlh, frame=sk.frame.LVLH)

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)

# 1-sigma position uncertainty in LVLH (x ≈ along-track, y = cross-track, z = nadir) sigma_pos_lvlh = np.array([50.0, 20.0, 10.0]) state.set_pos_uncertainty(sigma_pos_lvlh, frame=sk.frame.LVLH) # 1-sigma velocity uncertainty in LVLH sigma_vel_lvlh = np.array([0.05, 0.02, 0.01]) # m/s state.set_vel_uncertainty(sigma_vel_lvlh, frame=sk.frame.LVLH) 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.

In [4]:

# 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 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.18757395 854.68359056  19.99871425]

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.

In [5]:

# 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))

# 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))

In [6]:

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()

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()

The LVLH Frame¶

The LVLH (Local Vertical, Local Horizontal) frame is a coordinate system attached to the orbiting satellite, defined as:

• x-axis: approximately along-track (in the direction of velocity)
• y-axis: cross-track (opposite to the orbital angular momentum vector $\mathbf{h} = \mathbf{r} \times \mathbf{v}$)
• z-axis: nadir (pointing toward Earth center, i.e., $-\hat{r}$)

This frame is useful because orbital uncertainties have a natural physical interpretation in LVLH: along-track errors tend to grow fastest (due to period uncertainty), while radial errors remain relatively bounded.

The qgcrf2lvlh property on satstate provides the quaternion that rotates vectors from the GCRF inertial frame into the LVLH frame.

In [7]:

# Demonstrate the LVLH quaternion
q_lvlh = state.qgcrf2lvlh
print("Quaternion GCRF -> LVLH:", q_lvlh)

# Verify the LVLH axes:
# z-axis should point toward nadir (-r_hat)
r_hat = state.pos / np.linalg.norm(state.pos)
z_lvlh = q_lvlh * r_hat  # Should be [0, 0, -1]
print(f"\nr_hat in LVLH (expect ~[0, 0, -1]): {z_lvlh}")

# y-axis should be along -h_hat
h = np.cross(state.pos, state.vel)
h_hat = h / np.linalg.norm(h)
y_lvlh = q_lvlh * h_hat  # Should be [0, -1, 0]
print(f"h_hat in LVLH (expect ~[0, -1, 0]): {y_lvlh}")

# x-axis should be approximately along velocity
v_hat = state.vel / np.linalg.norm(state.vel)
x_lvlh = q_lvlh * v_hat  # Should be ~[1, 0, 0]
print(f"v_hat in LVLH (expect ~[1, 0, 0]): {x_lvlh}")

# Demonstrate the LVLH quaternion q_lvlh = state.qgcrf2lvlh print("Quaternion GCRF -> LVLH:", q_lvlh) # Verify the LVLH axes: # z-axis should point toward nadir (-r_hat) r_hat = state.pos / np.linalg.norm(state.pos) z_lvlh = q_lvlh * r_hat # Should be [0, 0, -1] print(f"\nr_hat in LVLH (expect ~[0, 0, -1]): {z_lvlh}") # y-axis should be along -h_hat h = np.cross(state.pos, state.vel) h_hat = h / np.linalg.norm(h) y_lvlh = q_lvlh * h_hat # Should be [0, -1, 0] print(f"h_hat in LVLH (expect ~[0, -1, 0]): {y_lvlh}") # x-axis should be approximately along velocity v_hat = state.vel / np.linalg.norm(state.vel) x_lvlh = q_lvlh * v_hat # Should be ~[1, 0, 0] print(f"v_hat in LVLH (expect ~[1, 0, 0]): {x_lvlh}")

Quaternion GCRF -> LVLH: Quaternion(Axis = [0.5774, 0.5774, -0.5774], Angle = 2.0944 rad)

r_hat in LVLH (expect ~[0, 0, -1]): [ 5.55111512e-16  1.11022302e-16 -1.00000000e+00]
h_hat in LVLH (expect ~[0, -1, 0]): [ 0.00000000e+00 -1.00000000e+00  3.33066907e-16]
v_hat in LVLH (expect ~[1, 0, 0]): [1.00000000e+00 4.44089210e-16 1.11022302e-16]

Covariance in the LVLH Frame¶

While satkit stores and propagates the covariance in GCRF, we can rotate it into the LVLH frame at any point to examine the uncertainty in physically meaningful orbit-relative directions.

In [8]:

def cov_gcrf_to_lvlh(sat):
    """Rotate the 3x3 position covariance from GCRF to LVLH."""
    q = sat.qgcrf2lvlh
    # Build rotation matrix by rotating basis vectors
    R = np.column_stack([
        q * np.array([1.0, 0.0, 0.0]),
        q * np.array([0.0, 1.0, 0.0]),
        q * np.array([0.0, 0.0, 1.0]),
    ])
    pos_cov_gcrf = sat.cov[0:3, 0:3]
    return R @ pos_cov_gcrf @ R.T


# Propagate and extract LVLH uncertainties
sigma_along = np.zeros(num_steps)
sigma_cross = np.zeros(num_steps)
sigma_radial = np.zeros(num_steps)

for i in range(num_steps):
    dt = sk.duration.from_minutes(dt_minutes * (i + 1))
    s = state.propagate(dt)
    cov_lvlh = cov_gcrf_to_lvlh(s)
    sigmas = np.sqrt(np.diag(cov_lvlh))
    sigma_along[i] = sigmas[0]   # x = along-track
    sigma_cross[i] = sigmas[1]   # y = cross-track
    sigma_radial[i] = sigmas[2]  # z = radial

def cov_gcrf_to_lvlh(sat): """Rotate the 3x3 position covariance from GCRF to LVLH.""" q = sat.qgcrf2lvlh # Build rotation matrix by rotating basis vectors R = np.column_stack([ q * np.array([1.0, 0.0, 0.0]), q * np.array([0.0, 1.0, 0.0]), q * np.array([0.0, 0.0, 1.0]), ]) pos_cov_gcrf = sat.cov[0:3, 0:3] return R @ pos_cov_gcrf @ R.T # Propagate and extract LVLH uncertainties sigma_along = np.zeros(num_steps) sigma_cross = np.zeros(num_steps) sigma_radial = np.zeros(num_steps) for i in range(num_steps): dt = sk.duration.from_minutes(dt_minutes * (i + 1)) s = state.propagate(dt) cov_lvlh = cov_gcrf_to_lvlh(s) sigmas = np.sqrt(np.diag(cov_lvlh)) sigma_along[i] = sigmas[0] # x = along-track sigma_cross[i] = sigmas[1] # y = cross-track sigma_radial[i] = sigmas[2] # z = radial

In [9]:

fig, ax = plt.subplots(figsize=(10, 5))

ax.plot(times_hr, sigma_along / 1e3, label="Along-track", linewidth=1.5)
ax.plot(times_hr, sigma_cross / 1e3, label="Cross-track", linewidth=1.5)
ax.plot(times_hr, sigma_radial / 1e3, label="Radial", linewidth=1.5)

ax.set_xlabel("Time (hours)")
ax.set_ylabel("Position 1-sigma (km)")
ax.set_title("Position Uncertainty in LVLH Frame Over 24 Hours")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(10, 5)) ax.plot(times_hr, sigma_along / 1e3, label="Along-track", linewidth=1.5) ax.plot(times_hr, sigma_cross / 1e3, label="Cross-track", linewidth=1.5) ax.plot(times_hr, sigma_radial / 1e3, label="Radial", linewidth=1.5) ax.set_xlabel("Time (hours)") ax.set_ylabel("Position 1-sigma (km)") ax.set_title("Position Uncertainty in LVLH Frame Over 24 Hours") ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() plt.show()

The LVLH-frame plot reveals the characteristic pattern of orbital uncertainty growth: along-track uncertainty dominates and grows roughly linearly with time, because a small error in velocity (and hence orbital period) causes the satellite to drift ahead or behind its predicted position. Cross-track and radial uncertainties remain comparatively bounded.