This notebook demonstrates an orbit determination from position data using a Keplerian propagator and a Batch-Least-Squares (BLS) estimator.
It also includes an IOD (Initial Orbit Determination) to avoid needing to know a first guess of the orbit.
Parameters
sigma_position = 100e3 # Noise (in terms of standard deviation of gaussian distribution) of input position data in meters
sigma_velocity = 100.0 # Noise of input velocity data in meters per second
# Estimator parameters
estimator_position_scale = 1.0 # m
estimator_convergence_thres = 1e-2
estimator_max_iterations = 25
estimator_max_evaluations = 35
# Orbit propagator parameters
prop_min_step = 0.001 # s
prop_max_step = 300.0 # s
prop_position_error = 10.0 # m
Importing generated position/velocity data
import pandas as pd
points_df = pd.read_csv('pos_vel_data_gcrf_with_noise.csv', index_col=0, parse_dates=True)
points_true_df = pd.read_csv('pos_vel_data_gcrf_without_noise.csv', index_col=0, parse_dates=True)
display(points_df)
x | y | z | vx | vy | vz | |
---|---|---|---|---|---|---|
2020-01-01 00:00:00 | 8.056982e+06 | 1.468446e+04 | 3.100794e+04 | -41.522207 | 3885.734219 | 6677.381032 |
2020-01-01 00:08:20 | 7.154730e+06 | 1.852806e+06 | 3.170712e+06 | -2908.078162 | 3434.859916 | 6098.294914 |
2020-01-01 00:16:40 | 5.006926e+06 | 3.606465e+06 | 5.901557e+06 | -5315.241457 | 2799.273333 | 4452.587151 |
2020-01-01 00:25:00 | 2.352366e+06 | 4.434280e+06 | 7.593867e+06 | -6313.109418 | 1530.104869 | 2432.364617 |
2020-01-01 00:33:20 | -1.007620e+06 | 5.033059e+06 | 8.323170e+06 | -6278.175584 | 296.658983 | 754.782685 |
2020-01-01 00:41:40 | -4.067582e+06 | 4.962326e+06 | 8.267137e+06 | -5918.025799 | -525.722689 | -1004.132238 |
2020-01-01 00:50:00 | -6.908941e+06 | 4.381691e+06 | 7.569041e+06 | -5101.186402 | -1353.391861 | -2258.428265 |
2020-01-01 00:58:20 | -8.900902e+06 | 3.333955e+06 | 5.922577e+06 | -3865.802334 | -1957.304014 | -3284.611687 |
2020-01-01 01:06:40 | -1.072952e+07 | 2.449153e+06 | 4.063725e+06 | -2690.374787 | -2389.261654 | -3898.051694 |
2020-01-01 01:15:00 | -1.170893e+07 | 1.477098e+06 | 2.010125e+06 | -1299.566595 | -2455.605459 | -4262.773492 |
2020-01-01 01:23:20 | -1.207381e+07 | -2.286726e+05 | -7.278299e+04 | 71.787359 | -2567.253309 | -4433.738988 |
2020-01-01 01:31:40 | -1.151708e+07 | -1.445310e+06 | -2.439107e+06 | 1461.300633 | -2516.046662 | -4368.852580 |
2020-01-01 01:40:00 | -1.046759e+07 | -2.608557e+06 | -4.459756e+06 | 2907.517698 | -2242.412883 | -3912.112590 |
2020-01-01 01:48:20 | -8.869487e+06 | -3.508512e+06 | -6.282518e+06 | 4138.495196 | -1696.232518 | -3211.194173 |
2020-01-01 01:56:40 | -6.532619e+06 | -4.101071e+06 | -7.719091e+06 | 5268.872249 | -1257.629013 | -2265.282454 |
2020-01-01 02:05:00 | -3.706500e+06 | -4.769382e+06 | -8.363096e+06 | 5737.178909 | -481.579749 | -954.340504 |
2020-01-01 02:13:20 | -6.329865e+05 | -5.049763e+06 | -8.386323e+06 | 6528.516821 | 562.970395 | 664.530033 |
2020-01-01 02:21:40 | 2.787555e+06 | -4.219518e+06 | -7.385754e+06 | 6229.740090 | 1492.106617 | 2848.933959 |
2020-01-01 02:30:00 | 5.307667e+06 | -3.143930e+06 | -5.741309e+06 | 4766.119029 | 2638.319578 | 4667.686305 |
2020-01-01 02:38:20 | 7.327372e+06 | -1.768722e+06 | -3.058312e+06 | 2596.747747 | 3467.515705 | 6077.323966 |
Firing up a JVM for Orekit
import orekit
orekit.initVM()
<jcc.JCCEnv at 0x7f420d36d150>
Downloading and importing the Orekit data ZIP
from orekit.pyhelpers import download_orekit_data_curdir, setup_orekit_curdir
download_orekit_data_curdir()
setup_orekit_curdir()
Downloading file from: https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip
from org.orekit.frames import FramesFactory
gcrf = FramesFactory.getGCRF()
from org.orekit.time import TimeScalesFactory
utc = TimeScalesFactory.getUTC()
Picking 3 vectors for orbit determination. These vectors are selected evenly-spaced in the dataframe.
import math
i = math.ceil(len(points_df.index) / 3)
points_for_iod = points_df.iloc[::i, :]
display(points_for_iod)
from org.hipparchus.geometry.euclidean.threed import Vector3D
from orekit.pyhelpers import datetime_to_absolutedate
pos_1 = points_for_iod.iloc[0]
vector_1 = Vector3D(pos_1[['x', 'y', 'z']].to_list())
date_1 = datetime_to_absolutedate(points_for_iod.index[0])
pos_2 = points_for_iod.iloc[1]
vector_2 = Vector3D(pos_2[['x', 'y', 'z']].to_list())
date_2 = datetime_to_absolutedate(points_for_iod.index[1])
pos_3 = points_for_iod.iloc[2]
vector_3 = Vector3D(pos_3[['x', 'y', 'z']].to_list())
date_3 = datetime_to_absolutedate(points_for_iod.index[2])
x | y | z | vx | vy | vz | |
---|---|---|---|---|---|---|
2020-01-01 00:00:00 | 8.056982e+06 | 1.468446e+04 | 3.100794e+04 | -41.522207 | 3885.734219 | 6677.381032 |
2020-01-01 00:58:20 | -8.900902e+06 | 3.333955e+06 | 5.922577e+06 | -3865.802334 | -1957.304014 | -3284.611687 |
2020-01-01 01:56:40 | -6.532619e+06 | -4.101071e+06 | -7.719091e+06 | 5268.872249 | -1257.629013 | -2265.282454 |
Performing the Initial Orbit Determination using Gibb's method. It assumes that at least 3 data points are available
from org.orekit.estimation.iod import IodGibbs
from org.orekit.utils import Constants as orekit_constants
iod_gibbs = IodGibbs(orekit_constants.EIGEN5C_EARTH_MU)
orbit_first_guess = iod_gibbs.estimate(gcrf,
vector_1, date_1,
vector_2, date_2,
vector_3, date_3)
from org.orekit.propagation.analytical import KeplerianPropagator
kepler_propagator_iod = KeplerianPropagator(orbit_first_guess)
display(orbit_first_guess)
<KeplerianOrbit: Keplerian parameters: {a: 9941292.946431652; e: 0.18973160346429957; i: 61.29160398126308; pa: 3.4221202861263906; raan: -0.5802427038357567; v: 139.49521421280545;}>
Setting up a numerical propagator. It is not possible in Orekit to perform orbit determination with a Keplerian propagator.
from org.orekit.propagation.conversion import DormandPrince853IntegratorBuilder
integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
from org.orekit.propagation.conversion import NumericalPropagatorBuilder
from org.orekit.orbits import PositionAngle
propagatorBuilder = NumericalPropagatorBuilder(orbit_first_guess,
integratorBuilder, PositionAngle.TRUE, estimator_position_scale)
from org.hipparchus.linear import QRDecomposer
matrix_decomposer = QRDecomposer(1e-11)
from org.hipparchus.optim.nonlinear.vector.leastsquares import GaussNewtonOptimizer
optimizer = GaussNewtonOptimizer(matrix_decomposer, False)
from org.orekit.estimation.leastsquares import BatchLSEstimator
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)
Feeding position measurements to the esimator
from orekit.pyhelpers import datetime_to_absolutedate
from org.orekit.estimation.measurements import Position, ObservableSatellite
observableSatellite = ObservableSatellite(0) # Propagator index = 0
for timestamp, pv_gcrf in points_df.iterrows():
orekit_position = Position(
datetime_to_absolutedate(timestamp),
Vector3D(pv_gcrf[['x', 'y', 'z']].to_list()),
sigma_position,
1.0, # Base weight
observableSatellite
)
estimator.addMeasurement(orekit_position)
Performing the orbit determination
estimatedPropagatorArray = estimator.estimate()
dt = 60.0
date_start = datetime_to_absolutedate(points_df.index[0])
date_end = datetime_to_absolutedate(points_df.index[-1])
# First propagating in ephemeris mode
estimatedPropagator = estimatedPropagatorArray[0]
estimatedInitialState = estimatedPropagator.getInitialState()
display(estimatedInitialState.getOrbit())
estimatedPropagator.resetInitialState(estimatedInitialState)
estimatedPropagator.setEphemerisMode()
estimatedPropagator.propagate(date_start, date_end)
bounded_propagator = estimatedPropagator.getGeneratedEphemeris()
<Orbit: Keplerian parameters: {a: 1.000312013771908E7; e: 0.19990357878201737; i: 59.84680615787619; pa: 0.5589419916700391; raan: -0.16034361509430223; v: 142.22691242122573;}>
Propagating the bounded propagator to retrieve the intermediate states
# BLS = batch least squares
import numpy as np
from orekit.pyhelpers import absolutedate_to_datetime
pv_bls_df = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
pv_iod_df = pd.DataFrame(columns=['x', 'y', 'z', 'vx', 'vy', 'vz'])
date_current = date_start
while date_current.compareTo(date_end) <= 0:
datetime_current = absolutedate_to_datetime(date_current)
spacecraftState = bounded_propagator.propagate(date_current)
pv_bls = spacecraftState.getPVCoordinates(gcrf)
pos_bls = np.array(pv_bls.getPosition().toArray())
pv_bls_df.loc[datetime_current] = np.concatenate(
(pos_bls,
np.array(pv_bls.getVelocity().toArray()))
)
pv_iod = kepler_propagator_iod.getPVCoordinates(date_current, gcrf)
pos_iod_gcrf = np.array(pv_iod.getPosition().toArray())
pv_iod_df.loc[datetime_current] = np.concatenate(
(pos_iod_gcrf,
np.array(pv_iod.getVelocity().toArray()))
)
date_current = date_current.shiftedBy(dt)
Plotting orbit in 3D.
import plotly.io as pio
pio.renderers.default = 'jupyterlab+notebook+png' # Uncomment for interactive plots
import plotly.graph_objects as go
fig_data = data=[go.Scatter3d(x=pv_iod_df['x'], y=pv_iod_df['y'], z=pv_iod_df['z'],
mode='lines',
name='IOD solution'),
go.Scatter3d(x=pv_bls_df['x'], y=pv_bls_df['y'], z=pv_bls_df['z'],
mode='lines',
name='Batch least squares solution'),
go.Scatter3d(x=points_df['x'], y=points_df['y'], z=points_df['z'],
mode='markers',
name='Measurements'),
go.Scatter3d(x=points_for_iod['x'], y=points_for_iod['y'], z=points_for_iod['z'],
mode='markers',
name='Measurements used for IOD')]
scene=dict(aspectmode='data', #this string can be 'data', 'cube', 'auto', 'manual'
)
layout = dict(
scene=scene
)
fig = go.Figure(data=fig_data,
layout=layout)
fig.show()
Computing residuals
residuals_df = pd.DataFrame(columns=['bls_minus_measurements_norm', 'iod_minus_measurements_norm', 'bls_minus_truth_norm', 'iod_minus_truth_norm'])
for timestamp, pv_gcrf in points_df.iterrows():
date_current = datetime_to_absolutedate(timestamp)
pv_bls = bounded_propagator.getPVCoordinates(date_current, gcrf)
pos_bls = np.array(pv_bls.getPosition().toArray())
pv_iod = kepler_propagator_iod.getPVCoordinates(date_current, gcrf)
pos_iod = np.array(pv_iod.getPosition().toArray())
pv_measurements = points_df.loc[timestamp]
pos_measurements = pv_measurements[['x', 'y', 'z']]
pv_true = points_true_df.loc[timestamp]
pos_true = pv_true[['x', 'y', 'z']]
bls_minus_measurements = np.linalg.norm(pos_bls - pos_measurements)
iod_minus_measurements = np.linalg.norm(pos_iod - pos_measurements)
bls_minus_truth = np.linalg.norm(pos_bls - pos_true)
iod_minus_truth = np.linalg.norm(pos_iod - pos_true)
residuals_df.loc[timestamp] = [
np.linalg.norm(pos_bls - pos_measurements),
np.linalg.norm(pos_iod - pos_measurements),
np.linalg.norm(pos_bls - pos_true),
np.linalg.norm(pos_iod - pos_true),
]
display(residuals_df)
bls_minus_measurements_norm | iod_minus_measurements_norm | bls_minus_truth_norm | iod_minus_truth_norm | |
---|---|---|---|---|
2020-01-01 00:00:00 | 56284.541463 | 4.406553e+05 | 35197.554591 | 410477.798023 |
2020-01-01 00:08:20 | 124848.068179 | 4.512250e+05 | 29525.788027 | 459466.806069 |
2020-01-01 00:16:40 | 229709.397325 | 7.092156e+05 | 32147.889070 | 469140.661240 |
2020-01-01 00:25:00 | 122488.157125 | 3.941168e+05 | 41267.597273 | 440068.128587 |
2020-01-01 00:33:20 | 177989.312521 | 5.735803e+05 | 51704.258637 | 383718.574137 |
2020-01-01 00:41:40 | 145384.854773 | 4.714415e+05 | 60720.042314 | 315624.280127 |
2020-01-01 00:50:00 | 199651.855773 | 4.289431e+05 | 67340.500439 | 254940.221995 |
2020-01-01 00:58:20 | 236697.521045 | 1.041250e-09 | 71344.845271 | 226318.598419 |
2020-01-01 01:06:40 | 68850.073865 | 2.212800e+05 | 72834.435338 | 246095.676748 |
2020-01-01 01:15:00 | 263136.928530 | 4.124667e+05 | 72060.792997 | 301602.606451 |
2020-01-01 01:23:20 | 203790.425959 | 4.508294e+05 | 69360.588820 | 370368.068233 |
2020-01-01 01:31:40 | 148366.500635 | 2.920698e+05 | 65139.093232 | 437821.516601 |
2020-01-01 01:40:00 | 97167.162597 | 3.922434e+05 | 59876.916042 | 495578.656280 |
2020-01-01 01:48:20 | 111680.778305 | 4.879147e+05 | 54142.044591 | 538013.368693 |
2020-01-01 01:56:40 | 318131.193403 | 4.935799e+05 | 48576.224429 | 560694.189462 |
2020-01-01 02:05:00 | 83406.147972 | 5.204905e+05 | 43791.551469 | 559942.593999 |
2020-01-01 02:13:20 | 180707.939038 | 6.596178e+05 | 40101.479799 | 533204.774547 |
2020-01-01 02:21:40 | 285624.069241 | 2.754855e+05 | 37148.959456 | 480520.331691 |
2020-01-01 02:30:00 | 218996.522755 | 4.825182e+05 | 33764.598152 | 407782.535158 |
2020-01-01 02:38:20 | 137352.096042 | 4.670766e+05 | 28364.897320 | 331800.248974 |
Showing the residuals, i.e. the distance between the measurement points and the estimated positions (by the IOD and the BLS respectively).
This tells how well the estimation fits the measurements, but not how well it fits the "true" orbit, because the measurements contain significant noise here.
fig = go.Figure(data=[
go.Scatter(x=residuals_df.index, y=residuals_df['bls_minus_measurements_norm'],
name='BLS - measurements',
mode='markers+lines'),
go.Scatter(x=residuals_df.index, y=residuals_df['iod_minus_measurements_norm'],
name='IOD - measurements',
mode='markers+lines')
])
fig.show()
Showing the "estimation error", i.e. the difference between the truth (used to generate the measurement data) and the estimation.
The batch least squares estimation gives better results than the measurement noise: it is able to "filter out" the noise a little bit as it fits the ellipse.
fig = go.Figure(data=[
go.Scatter(x=residuals_df.index, y=residuals_df['bls_minus_truth_norm'], name='BLS - truth', mode='markers+lines'),
go.Scatter(x=residuals_df.index, y=residuals_df['iod_minus_truth_norm'], name='IOD - truth', mode='markers+lines'),
])
fig.show()