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

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

# Fetch the current TLE for INTELSAT 19 (NORAD 38356) directly from CelesTrak.
# IS-19 is a geostationary communications satellite at 166°E longitude,
# visible from the SST.
iss19_url = "https://celestrak.org/NORAD/elements/gp.php?CATNR=38356&FORMAT=TLE"
tle = sk.TLE.from_url(iss19_url)
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()

# Measurements taken once every 5 minutes over a 7-hour window starting
# 9 hours after the TLE epoch (approximating a single night of observation).
sample_start = epoch + sk.duration(hours=9)
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

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

# Fetch the current TLE for INTELSAT 19 (NORAD 38356) directly from CelesTrak.
# IS-19 is a geostationary communications satellite at 166°E longitude,
# visible from the SST.
iss19_url = "https://celestrak.org/NORAD/elements/gp.php?CATNR=38356&FORMAT=TLE"
tle = sk.TLE.from_url(iss19_url)
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()

# Measurements taken once every 5 minutes over a 7-hour window starting
# 9 hours after the TLE epoch (approximating a single night of observation).
sample_start = epoch + sk.duration(hours=9)
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-04-05T13:40:12.728064Z
Iteration 0
Iteration 1
Iteration 2
Initial State Error: 5.319 km, 0.133 m/s
Final State Error: 2.521 km, 0.261 m/s