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

In [1]:
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

In [2]:
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

In [3]:
import orekit
orekit.initVM()
Out[3]:
<jcc.JCCEnv at 0x7f420d36d150>

Downloading and importing the Orekit data ZIP

In [4]:
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
In [5]:
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.

In [6]:
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

In [7]:
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.

In [8]:
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

In [9]:
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

In [10]:
estimatedPropagatorArray = estimator.estimate()
In [11]:
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

In [12]:
# 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.

In [13]:
import plotly.io as pio
pio.renderers.default = 'jupyterlab+notebook+png'  # Uncomment for interactive plots
In [14]:
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()