Optical Observations of Satellites

This provides an example of using optical observations of satellites (e.g., from a telescope on the ground) to update the state of the satellite. This is commonly used to update satellite ephemeris, in particular for satellites in MEO and higher orbits that lack easy access to the service volume of the GPS contellation and can be observed for long durations at night without being shadowed by the Earth. The US government and several commercial companies such as Exoanalytics make use of these techniques to maintain the ephemeris on these satellites.

This example uses a batch least-squares technique for the fitting. This can also be performed with sequential filters such as an Extended or Unscented Kalman Filter.

This example updates the initial state vector to mimimizes the difference between observed and predicted line of site to the satellite for all observations:

\[\hat{\mathbf x}_0 = \arg\min_{\mathbf x_0\in\mathbb R^6} \sum_{k=1}^{N} \left\| E_k\Big(\hat{\mathbf u}_m(\alpha_k,\delta_k)-\hat{\mathbf u}(\phi(t_k,t_0;\mathbf x_0),\mathbf r_{0}(t_k))\Big) \right\|_{R^{-1}}^{2}\]

where

\[\begin{split}\hat{\mathbf u}_m(\alpha,\delta)= \begin{bmatrix} \cos\delta\cos\alpha\\ \cos\delta\sin\alpha\\ \sin\delta \end{bmatrix}\end{split}\]

is the measurement (\(\alpha\) is right ascension and \(\delta\) is declination),

\[\hat{\mathbf u}(\mathbf x_k,\mathbf r_0)= \frac{\mathbf r_k-\mathbf r_0}{\|\mathbf r_k-\mathbf r_0\|}\]

is the predicted unit vector pointing from the observer at \(\mathbf r_0\) to the satellite at \(\mathbf r_k\)

\[\mathbf x_k=\phi(t_k,t_0;\mathbf x_0)\]

is the propagated state at time \(t_k\), given initial state \(\mathbf x_0\)

\(E\) is just a rotation into the local tangent plane for each observation (to avoid singularities and keep measurement errors constant in local measurement plane)

import satkit as sk
import numpy as np
import math as m
import datetime

# Take our observer to be the Space Surveillance Telescope (SST)
# in Exmouth, Western Australia
sst_lat = -22.675
sst_lon = 114.235
sst_alt = 0.5 # km
sst = sk.itrfcoord(latitude_deg=sst_lat, longitude_deg=sst_lon, alt_m=sst_alt*1e3)

# Lets look at INTELSAT 19 (IS-19), a geostationary communications satellite.
# It is currently located at 166E longitude, so it should be visible from the SST.
# We can get the TLE for INTELSAT 19 from Celestrak:
# https://celestrak.com/NORAD/elements/geo.txt
tle = sk.TLE.from_lines([
    "INTELSAT 19 (IS-19)",
    "1 38356U 12030A   26045.51850411 -.00000051  00000+0  00000+0 0  9998",
    "2 38356   0.0150  97.8441 0002877 233.5096 165.8257  1.00272325 49221"
])
epoch = tle.epoch

# Lets make the initial "true" state the tle-predicted state at epoch
pteme, vteme = sk.sgp4(tle, epoch)
q = sk.frametransform.qteme2gcrf(epoch)
pgcrf = q * pteme
vgcrf = q * vteme
state0 = np.hstack((pgcrf, vgcrf))
print(epoch)
begin_time = epoch
end_time = epoch + sk.duration(days=2)
settings = sk.propsettings()
settings.precompute_terms(begin_time, end_time)
# True state over 1 day
truth = sk.propagate(state0, begin_time, end_time, propsettings=settings)

# Assume position knowledge to 10 km, velocity to 5 cm/s
pos_noise = 10e3
vel_noise = 0.05

# initial state estimate
state0_est = state0 + np.random.normal(scale=[pos_noise]*3 + [vel_noise]*3)
state0_est_prior = state0_est.copy()

# sample times for measurements.  Every 5 minues for 6 hours beginning at local 9pm and going to local 4am
sample_start = sk.time(2026, 2, 15, 20, 0, 0) # local 9pm at SST
sample_end = sample_start + sk.duration(hours=7)
sample_times = [sample_start + sk.duration(minutes=5)*i for i in range(7*60//5)]

# Truth at first sample time
state0_truth = truth.interp(sample_times[0])

# the optical measurement is a line-of-sight vector from observer to satellite in the inertial frame
sst_pos = [sk.frametransform.qitrf2gcrf(t) * sst.vector for t in sample_times]
sat_state = [truth.interp(t) for t in sample_times]
los_meas = [sat_state[idx][0:3] - sst_pos[idx] for idx in range(len(sample_times))]
# Normalize the line-of-sight measurements to get unit vectors
los_meas = [los / np.linalg.norm(los) for los in los_meas]

# add noise to the measurements.  Assume 30 microradians of angular noise
# Note: I made up 30 microradians (6 arcseconds).
ang_noise = 30e-6

def add_ang_noise(u, ang_noise):
    # add noise in a random direction perpendicular to the line of sight
    # first get a random vector
    rand_vec = np.random.normal(size=3)
    # make it perpendicular to the line of sight
    rand_vec -= rand_vec.dot(u) * u
    rand_vec /= np.linalg.norm(rand_vec)
    # rotate the line of sight by the noise angle in the direction of the random vector
    theta = np.random.normal(scale=ang_noise)
    u_noisy = u * m.cos(theta) + np.cross(rand_vec, u) * \
        m.sin(theta) + rand_vec * rand_vec.dot(u) * (1 - m.cos(theta))
    return u_noisy

los_meas_noisy = [add_ang_noise(los, ang_noise) for los in los_meas]


# Now, make up an initial state estimate that we are going to refine with measurements
state0_est = truth.interp(sample_times[0]) + np.random.normal(scale=[pos_noise]*3 + [vel_noise]*3)

# Our initial state error .. record for later comparison
state0_est_prior = state0_est.copy()

def unit_vector_jacobian(rho):
    """
    Jacobian of unit vector with respect to the original vector
    """
    rho_norm = np.linalg.norm(rho)
    u = rho / rho_norm
    return (np.eye(3) - np.outer(u, u)) / rho_norm


def tangent_basis(u):
    """
    RA/Dec tangent basis at predicted direction
    """
    ux, uy, uz = u
    alpha = np.arctan2(uy, ux)
    delta = np.arcsin(uz)

    e_alpha = np.array([-np.sin(alpha), np.cos(alpha), 0.0])
    e_delta = np.array([
        -np.sin(delta)*np.cos(alpha),
        -np.sin(delta)*np.sin(alpha),
        np.cos(delta)
    ])
    return np.vstack((e_alpha, e_delta))

# OK, lets do it!
# 3 iterations of Gauss-Newton least squares
for n in range(0,3):
    print(f"Iteration {n}")

    # Propagate state to get state transition matrix at each time
    res = sk.propagate(state0_est, sample_times[0], sample_times[-1], output_phi=True, propsettings=settings)
    [state,phi] = zip(*[res.interp(t, output_phi=True) for t in sample_times])

    residuals = []
    J_list = []

    for idx, t in enumerate(sample_times):
        # Get the predicted line of sight measurement at this time
        p_sat = state[idx][0:3]
        p_obs = sst_pos[idx]
        los_vec = p_sat - p_obs
        los_pred = los_vec / np.linalg.norm(los_vec)

        # Get the tangent basis in predicted direction
        E = tangent_basis(los_pred)

        # Get residuals
        r = E @ (los_meas_noisy[idx] - los_pred)
        residuals.append(r)

        # IMPORTANT: Jacobian must be evaluated on un-normalized LOS vector
        U = unit_vector_jacobian(los_vec)
        # A is our state transition matrix for this time,
        # which maps changes in initial state to changes in state at this time.
        # We only care about the position part of the state
        # since the measurement is only a function of position.
        A = phi[idx][0:3, :]
        Hk = - E @ U @ A
        J_list.append(Hk)

    # OK, now do the least squares solve to get state update
    r = np.hstack(residuals)
    H = np.vstack(J_list)
    # Solve for state update using least squares
    dx = np.linalg.lstsq(H, r, rcond=None)[0]
    # Adjust the initial state estimate by the state update
    state0_est -= dx

print(f"Initial State Error: {np.linalg.norm(state0_est_prior - state0_truth)/1e3:.3f} km, {np.linalg.norm(state0_est_prior[3:6] - state0_truth[3:6]):.3f} m/s")
print(f'Final State Error: {np.linalg.norm(state0_est - state0_truth)/1e3:.3f} km, {np.linalg.norm(state0_est[3:6] - state0_truth[3:6]):.3f} m/s')

2026-02-14T12:26:38.755104Z
Iteration 0
Iteration 1
Iteration 2
Initial State Error: 13.335 km, 0.091 m/s
Final State Error: 2.129 km, 0.228 m/s