Hi,
First, thank you for this extremely useful Python wrapper for Orekit!
I was recently playing with it (I never used Orekit directly), and came across some behavior for which I am not sure I completely understand. I am using the 10.3.1 version of Orekit, installed from conda.
I have noisy RA/Dec measurements of satellites at different times, and I would like to extract the orbital parameters. For this I use the Laplace angles-only initial orbit determination, assuming Keplerian motion (src/main/java/org/orekit/estimation/iod/IodLaplace.java · develop · Orekit / Orekit · GitLab). But the output are not as expected.
I illustrate below the problem using simulated data:
- I take a known satellite TLE
- I generate positions using the Skyfield library at various time steps
- From these positions I run the initial orbit determination
At this stage, the output agrees with the reported orbital parameters of the satellite. But now, to mimick the fact that I have noisy measurement, I just slightly perturbed the input times before running the initial orbit determination (hundreds of millisecond scatter at each time step, while time steps are separated by several seconds). And now the output are completely off (sometimes negative semi-major axis, etc.).
I tried to inject uncertainty in the input positions (using the AngularRaDec
), but it does not change anything to the estimated parameters.
So my question is: how to account for the fact that measurements have scatter when estimating the initial orbital parameters?
You will find below a self-contained example to reproduce the findings:
#initialize orekit and JVM
import orekit
orekit.initVM()
# setup the orekit data loading, the file orekit-data.zip shall be in same directory as notebook.
from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()
from orekit import JArray_double
from org.orekit.utils import PVCoordinates
from org.orekit.utils import Constants
from org.orekit.utils import IERSConventions
from org.orekit.time import AbsoluteDate
from org.orekit.time import TimeScalesFactory
from org.orekit.bodies import GeodeticPoint
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.frames import FramesFactory;
from org.orekit.frames import TopocentricFrame
from org.orekit.estimation.iod import IodLaplace
from org.orekit.estimation.measurements import GroundStation;
from org.orekit.estimation.measurements import AngularRaDec;
from org.orekit.estimation.measurements import ObservableSatellite
from org.orekit.orbits import KeplerianOrbit;
from org.hipparchus.geometry.euclidean.threed import Vector3D
import pandas as pd
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord, EarthLocation
from astropy.time import Time
from skyfield.api import EarthSatellite, Topos
from skyfield.api import load
def generate_satellite_positions(satellite, dates_isot, lat='33.3563 N', lon='116.8650 W'):
""" Return RA/Dec Angle of the satellite for the dates.
Default observation site is Palomar
Parameters
----------
satellite: EarthSatellite
EarthSatellite instance
dates_isot: list of str
list of dates as ISO UTC
lat: str
Latitude of the observing site
lon: str
Longitude of the observing site
Returns
----------
ras: list of Angle
decs: list of Angle
"""
# Observations are done from earth
site = Topos(lat, lon)
difference = satellite - site
# wrap dates
ts = load.timescale()
dates = [ts.from_astropy(Time(x)) for x in dates_isot]
# Compute RA/Dec for satellites during the night
positions = [difference.at(t).radec()[0:2] for t in dates]
ras, decs = np.transpose(positions)
return ras, decs
class MovingObj():
""" Moving object container """
def __init__(self, dates, ra, dec):
""" Initialise a moving object
Parameters
----------
dates: list of str
list of dates as ISO UTC
ra: list of float
List of RA measurements (radians)
dec: list of float
List of Dec measurements (radians)
"""
self.dates = dates
self.ra = ra
self.dec = dec
def locate_observer(self, site='Palomar'):
""" Define an observation site
Parameters
----------
site: str
Any astropy known site name (case-insensitive). Default is Palomar
"""
self.site = site
# Observatory
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False);
body = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
location = EarthLocation.of_site('Palomar')
lat, lon, alt = float(location.lat.rad), float(location.lon.rad), float(location.height.value)
self.observer = GroundStation(
TopocentricFrame(body, GeodeticPoint(lat, lon, alt), site)
);
self.observer.getPrimeMeridianOffsetDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
self.observer.getPolarOffsetXDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
self.observer.getPolarOffsetYDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
def fit_orbit(self, scatter=1.0):
""" extract initial orbit parameters using Laplace method
Parameters
----------
scatter: float
sigma theoretical standard deviation. Default is 1.0
"""
gcrf = FramesFactory.getGCRF()
self.obsDate = []
self.raDecs = []
for date_, ra_, dec_ in zip(self.dates, self.ra, self.dec):
d1 = Time(date_).datetime
date1 = AbsoluteDate(
d1.year, d1.month, d1.day,
d1.hour, d1.minute, float(d1.second + d1.microsecond/1e6),
TimeScalesFactory.getUTC()
)
radec1 = JArray_double([ra_, dec_])
sig = JArray_double([scatter, scatter])
w = JArray_double([1.0, 1.0])
raDec1 = AngularRaDec(
self.observer, gcrf, date1,
radec1, # measurement
sig, # sigma theoretical standard deviation
w, # weight
ObservableSatellite(0)
)
self.obsDate.append(date1)
self.raDecs.append(raDec1)
laplace = IodLaplace(Constants.EGM96_EARTH_MU)
self.obsPva = self.observer.getBaseFrame().getPVCoordinates(self.raDecs[1].getDate(), gcrf);
# Use the 3 first measurements
self.orbit = laplace.estimate(
gcrf, self.obsPva,
self.obsDate[0], self.lineOfSight(self.raDecs[0]),
self.obsDate[1], self.lineOfSight(self.raDecs[1]),
self.obsDate[2], self.lineOfSight(self.raDecs[2]),
);
self.initialOrbit = KeplerianOrbit(
self.orbit.getPVCoordinates(),
FramesFactory.getEME2000(), self.obsDate[0],
Constants.EIGEN5C_EARTH_MU
)
@staticmethod
def lineOfSight(radec):
""" compute the line of sight for a given ra/dec
Note: this is included in orekit 11, but not before
Parameters
----------
radec: AngularRaDec
AngularRaDec instance
"""
alpha, delta = radec.getObservedValue()
los = Vector3D(alpha, delta)
return los
and just run
ts = load.timescale()
line1 = '1 13652U 82110C 19305.11481830 -.00000186 +00000-0 +00000-0 0 9999'
line2 = '2 13652 014.8379 353.7256 0003739 341.0436 021.0683 00.99899227022330'
satellite = EarthSatellite(line1, line2, 'ANIK C3 (TELESAT-5)', ts)
print(satellite)
# ANIK C3 (TELESAT-5) catalog #13652 epoch 2019-11-01 02:45:20 UTC
# Perfect case -- all good
dates_isot = ['2019-11-02T12:37:19Z', '2019-11-02T12:37:25Z', '2019-11-02T12:37:31Z']
ra, dec = generate_satellite_positions(satellite, dates_isot)
sat = MovingObj(dates_isot, [i.radians for i in ra], [i.radians for i in dec])
sat.locate_observer('Palomar')
sat.fit_orbit()
print(sat.initialOrbit)
# Keplerian parameters: {a: 4.22820321049187E7; e: 1.5653189820011832E-4;
# i: 14.824427854453244; pa: 10.526587409216345; raan: -6.165885817788712;
# v: 138.722440746431;}
# Perturbed case -- all off!
dates_isot = ['2019-11-02T12:37:19Z', '2019-11-02T12:37:25Z', '2019-11-02T12:37:31Z']
ra, dec = generate_satellite_positions(satellite, dates_isot)
# scatter the times by hundreds of millisecond
dates_isot = ['2019-11-02T12:37:18.9Z', '2019-11-02T12:37:24.8Z', '2019-11-02T12:37:31.1Z']
sat = MovingObj(dates_isot, [i.radians for i in ra], [i.radians for i in dec])
sat.locate_observer('Palomar')
sat.fit_orbit()
print(sat.initialOrbit)
# Keplerian parameters: {a: -9593.956621851428; e: 517.8382836818042;
# i: 43.176695972549574; pa: 85.61033112962203; raan: -27.993882345720124;
# v: 83.43512670532412;}
Note that if I increase the sigma theoretical standard deviation, it does not change the result:
sat.fit_orbit(scatter=10.0)
print(sat.initialOrbit)
# Keplerian parameters: {a: -9593.956621851428; e: 517.8382836818042;
# i: 43.176695972549574; pa: 85.61033112962203; raan: -27.993882345720124;
# v: 83.43512670532412;}