TLE Fitting¶

Two-Line Element sets (TLEs) are the standard format for distributing satellite orbital elements, but they degrade in accuracy over time. When higher-fidelity state vectors are available (e.g., from GPS or precision orbit determination), it is useful to fit a new TLE to those states.

satkit.TLE.fit_from_states performs this fit using Levenberg-Marquardt non-linear least-squares optimization. It tunes the TLE orbital parameters to minimize the difference between the input state positions and the SGP4-predicted positions. Input states are assumed to be in the GCRF frame and are internally rotated to the TEME frame used by SGP4.

In [1]:

# Imports
import satkit as sk
import numpy as np
import math as m

# Create a high-precision state
# Altitude for circular orbit
altitude = 450e3

# Radius & velocity
r0 = altitude + sk.consts.earth_radius
v0 = m.sqrt(sk.consts.mu_earth / r0)

# Inclination
inclination = 15 * m.pi / 180.0

# Create the state (3D position in meters, 3D velocity in meters / second)
state0 = np.array([r0, 0, 0, 0, v0 * m.cos(inclination), v0 * m.sin(inclination)])
# Make up an epoch
time0 = sk.time(2024, 3, 15, 13, 0, 0)

# Propagate the state forward by a day with high-precision propagator
res = sk.propagate(state0, time0, time0 + sk.duration(days=1.0))

# Get interpolated states every 10 minutes
times = [time0 + sk.duration(minutes=i) for i in range(0, 1440, 10)]
states = [res.interp(t) for t in times]

# Fit the TLE
(tle, fitresults) = sk.TLE.fit_from_states(states, times, time0 + sk.duration(days=0.5)) # type: ignore

# Print the result
print(tle)
print(fitresults['success'])

# Imports import satkit as sk import numpy as np import math as m # Create a high-precision state # Altitude for circular orbit altitude = 450e3 # Radius & velocity r0 = altitude + sk.consts.earth_radius v0 = m.sqrt(sk.consts.mu_earth / r0) # Inclination inclination = 15 * m.pi / 180.0 # Create the state (3D position in meters, 3D velocity in meters / second) state0 = np.array([r0, 0, 0, 0, v0 * m.cos(inclination), v0 * m.sin(inclination)]) # Make up an epoch time0 = sk.time(2024, 3, 15, 13, 0, 0) # Propagate the state forward by a day with high-precision propagator res = sk.propagate(state0, time0, time0 + sk.duration(days=1.0)) # Get interpolated states every 10 minutes times = [time0 + sk.duration(minutes=i) for i in range(0, 1440, 10)] states = [res.interp(t) for t in times] # Fit the TLE (tle, fitresults) = sk.TLE.fit_from_states(states, times, time0 + sk.duration(days=0.5)) # type: ignore # Print the result print(tle) print(fitresults['success'])

            TLE: none
                         NORAD ID: 00000,
                      Launch Year: 2000,
                            Epoch: 2024-03-16T01:00:00.000000Z,
                  Mean Motion Dot: 0 revs / day^2,
              Mean Motion Dot Dot: 0 revs / day^3,
                             Drag: -0.00015321661361392396,
                      Inclination: 14.991373022912788 deg,
                             RAAN: 355.9976365777418 deg,
                            eccen: 0.0013613551572524063,
                   Arg of Perigee: 199.31610033701583 deg,
                     Mean Anomaly: 62.0109822387111 deg,
                      Mean Motion: 15.408691924113539 revs / day
                            Rev #: 0
        
Convergence in parameter value

Generate Test Data¶

To demonstrate the fitting, we create a synthetic truth trajectory by propagating a circular orbit at 450 km altitude with satkit's high-precision propagator. We then sample position and velocity states every 10 minutes over one day, and fit a TLE to those samples.

In [2]:

# Compute position errors (differences between TLE & state)

# Get the positions from sgp4
(pteme, vteme) = sk.sgp4(tle, times)
# Rotate positions from TEME to GCRF frame
pgcrf = [sk.frametransform.qteme2gcrf(t) * p for t, p in zip(times, pteme)]
# Take difference between state vector and SGP4 positions, and compute norm
pdiff = [p - s[0:3] for p, s in zip(pgcrf, states)]
pdiff = np.array([np.linalg.norm(p) for p in pdiff])


# Plot position errors
import plotly.graph_objects as go

fig = go.Figure()
fig.add_trace(go.Scatter(x=[t.datetime() for t in times], y=pdiff, mode='lines', name='Position Error',
                         line=dict(color='black', width=2)))
fig.update_layout(title='TLE Fitting Position Errors',
                  xaxis_title='Time',
                  yaxis_title='Position Error (m)')
fig.update_xaxes(showline=True, linewidth=2, linecolor="black", mirror=True)
fig.update_yaxes(showline=True, linewidth=2, linecolor="black", mirror=True)
fig.update_layout(
    xaxis=dict(
        gridcolor="#dddddd",
        gridwidth=1,
    ),
    yaxis=dict(
        gridcolor="#dddddd",
        gridwidth=1,
    ),
)

fig.show()

# Compute position errors (differences between TLE & state) # Get the positions from sgp4 (pteme, vteme) = sk.sgp4(tle, times) # Rotate positions from TEME to GCRF frame pgcrf = [sk.frametransform.qteme2gcrf(t) * p for t, p in zip(times, pteme)] # Take difference between state vector and SGP4 positions, and compute norm pdiff = [p - s[0:3] for p, s in zip(pgcrf, states)] pdiff = np.array([np.linalg.norm(p) for p in pdiff]) # Plot position errors import plotly.graph_objects as go fig = go.Figure() fig.add_trace(go.Scatter(x=[t.datetime() for t in times], y=pdiff, mode='lines', name='Position Error', line=dict(color='black', width=2))) fig.update_layout(title='TLE Fitting Position Errors', xaxis_title='Time', yaxis_title='Position Error (m)') fig.update_xaxes(showline=True, linewidth=2, linecolor="black", mirror=True) fig.update_yaxes(showline=True, linewidth=2, linecolor="black", mirror=True) fig.update_layout( xaxis=dict( gridcolor="#dddddd", gridwidth=1, ), yaxis=dict( gridcolor="#dddddd", gridwidth=1, ), ) fig.show()

Evaluate Fit Quality¶

Compare the fitted TLE against the original states by propagating the TLE with SGP4, rotating from TEME to GCRF, and computing position differences. Since TLEs are a simplified analytical model, some residual error is expected even with a perfect fit.