Optical Observations of Satellites¶

Ground-based telescopes can observe sunlit satellites against the night sky, measuring their angular position (right ascension and declination). These line-of-sight measurements can then refine the satellite's orbit. This technique is particularly valuable for satellites in MEO and GEO orbits, which are visible for long periods at night and lack easy access to GPS signals. Organizations like the US Space Force and commercial companies such as ExoAnalytic routinely use optical tracking to maintain satellite catalogs.

Approach¶

This example uses a batch least-squares (Gauss-Newton) technique to update an initial state estimate by minimizing the difference between observed and predicted line-of-sight unit vectors:

$$ \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:

$\hat{\mathbf u}_m(\alpha,\delta) = [\cos\delta\cos\alpha,\; \cos\delta\sin\alpha,\; \sin\delta]^T$ is the measured unit vector from right ascension $\alpha$ and declination $\delta$
$\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 from observer at $\mathbf r_0$ to 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_k$ projects the 3D residual into the local tangent plane (RA/Dec directions) to keep the measurement space 2D and avoid singularities

The state transition matrix $\Phi$ from the high-precision propagator maps initial state perturbations to perturbations at each measurement time, forming the Jacobian of the linearized system.

In [1]:

Copied!





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')
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: 5.706 km, 0.062 m/s
Final State Error: 0.758 km, 0.080 m/s