Initial orbit determination and noisy measurements

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:

  1. I take a known satellite TLE
  2. I generate positions using the Skyfield library at various time steps
  3. 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;}

Actually, it is even worse – perturbing the times by one millisecond is enough to produce a wrong result:

# Perturbed case
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 1 millisecond
dates_isot = ['2019-11-02T12:37:18.999Z', '2019-11-02T12:37:24.999Z', '2019-11-02T12:37:31.001Z']
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: 5.118234929517522E7; e: 0.35231887307405213; 
# i: 15.782594092433643; pa: 71.05876538242521; raan: -8.241984222407815; 
# v: 80.19334429060021;}

Hi @JulienPeloton

Welcome to the Orekit forum!

Results are still wrong, but there are better.
I don’t think that introducing a bias in the measurement epoch is a good idea. Generally, playing with time is a bad idea. Especially with IOD methods which are not accurate and sometimes erratic if the measurements are very noisy :slight_smile:

Did you already try to use the measurement generation functionalities? They are available in the org.orekit.estimation.measurements.generation package of Orekit. I think those functionalities could be very useful to generate noisy measurements. Also, you can find examples on those functionalities in different places (sorry it’s in Java but it can be adapted to Python):

  1. In the class tests.
  2. In the tutorials

Maybe you can also try to directly perform an orbit determination.

Best regards,
Bryan

1 Like

Hi @bcazabonne – thank you for your reply and the pointers.

Actually, I perturbed the time in my example because I wanted to illustrate my real-life case on simulated data (I have accurate spatial measurements of a satellite passing by, but with inaccurate timings). But eventually perturbing time or space is equivalent, and leads to inaccurate results for small deviation:

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

# perturb ra/dec by up to one arcsecond
np.random.seed(100)
ra_ = [i.radians + np.random.rand() * 1./3600 * np.pi/180 for i in ra]
dec_ = [i.radians + np.random.rand() * 1./3600 * np.pi/180 for i in dec]

sat = MovingObj(dates_isot, ra_, dec_)
sat.locate_observer('Palomar')
sat.fit_orbit()
print(sat.initialOrbit)
# Keplerian parameters: {a: -1496705.2568740027; e: 3.7268028259516863; 
# i: 46.115980795547564; pa: 97.76292281464372; raan: -29.889960544006804; 
# v: -69.58166186191569;}

I am not sure I completely understand how to use them, especially from the Python wrapper. My java skills are quite low, and any help would be appreciated here :slight_smile:

I am not sure I understand what you suggest – how this would be different from the Laplace orbit determination I am using here?

Ok, I understand now.
So, generating measurements using the measurement generation functionalities will not be useful in you case.

How many measurements do you have?
If you have only three measurements, I understand the use of an IOD method. However, IOD methods are not precise. As a results, if your measurements are not precise too (i.e., very noisy) it is “normal” to have wrong solution when using IOD.

IOD are generally used to calculate an initial guess for an orbit determination process. Indeed, IOD methods use empirical equations made with simplification and approximations. The fact that they generally assume a Keplerian motion is one of the most important simplification.

if you have enough measurements, I strongly recommend you to use an estimation algorithm to determine your orbit. At the opposite of IOD methods, estimation algorithms optimize the problem to find the better estimate of the orbit using the observed measurements. In addition, thanks to orbit determination, you can have a better representation of the dynamic of the satellite by adding orbital pertubations (i.e., atmospheric drag, third body attraction, earth gravity perturbations, …). Orekit proposes two algorithms: a batch least squares and a Kalman Filter. For your example, I think a batch least squares algorithm is the best candidate.

If you are Python user, I recommend you to look at the tutorial developed by Clément Jonglez about batch least squares orbit determination using Orekit. I think that you can adapt it to your needs.

Best regards,
Bryan

3 Likes

Hi @bcazabonne – thank you for the detailed explanations! I will go through the tutorial, and see if I can get out a better result. I think the biggest difficulty will be to guess an initial orbit, but I can probably feed the output of the IOD.

Hi @bcazabonne

I followed the tutorial you mentioned, but unfortunately I encounter a Jacobian matrix for type KEPLERIAN is singular with current orbit error. I found some posts in the forum dealing with this, but I couldn’t get rid of this error, even in a perfect measurement case (see below). I am wondering if the data I have is just hopeless: single small arc (up to one degree) spanning 30 seconds maximum (exposure time)?

Here is an example, where I first use IOD to guess a first orbit, and then use it in a batch least squares orbit determination. I briefly highlight important steps here, but the full script is attached below.

I have measurements corresponding to simulations of the ANIK C3 (TELESAT-5)

dates_isot = [
    '2019-11-02T12:37:19Z', '2019-11-02T12:37:25Z', 
    '2019-11-02T12:37:30Z', '2019-11-02T12:37:34Z',
    '2019-11-02T12:37:40Z', '2019-11-02T12:37:49Z'
]

ra = [
    2.5834605698194246, 2.583886840258935, 
    2.5842420571513363, 2.584526225381548, 
    2.584952463685652, 2.5855918049260747
]

dec = [
    0.054579013415419424, 0.05447025349530273, 
    0.0543796003532739, 0.05430706477282546, 
    0.054198240976653426, 0.05403495544463993
]

If I run the IOD (using the first 3 measurements), I get a rather good estimate, at least this is the good order of magnitude:

Keplerian parameters: {
a: 4.227102759246746E7; e: 0.0014381077442844517; i: 14.820109677796811; 
pa: 254.52333056539018; raan: -6.155828532482599; v: -105.28402914694347;}

Reported values from n2yo:
{a: 4.2270E7, e: 8E-4, i: 14.3}

Then I use this initial orbit for the the full orbit estimation. Let’s first define the orbit propagator parameters. We will take the same values as in the tutorial, except the prop_max_step which is set at 1.0 second here. Indeed, our data is made of 6 measurements spanning only 30 seconds total.

# Orbit propagator parameters
prop_min_step = 0.001 # s
prop_max_step = 1.0 # s
prop_position_error = 1.0 # m

# Estimator parameters
estimator_position_scale = 100.0 # m
estimator_convergence_thres = 1e-3
estimator_max_iterations = 25
estimator_max_evaluations = 35

Let’s set up the propagator and the estimator, by using the previous orbit initialOrbit from the IOD as a first guess:

integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)

# we use the initial orbit as first guess
propagatorBuilder = NumericalPropagatorBuilder(
    initialOrbit,
    integratorBuilder, 
    PositionAngle.MEAN, 
    estimator_position_scale
)
matrixDecomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrixDecomposer, False)

estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)

# we have 6 measurements in total
for radec in raDecs:
    estimator.addMeasurement(radec)

And then estimate:

estimatedPropagatorArray = estimator.estimate()
---------------------------------------------------------------------------
JavaError                                 Traceback (most recent call last)
<ipython-input-7-23bf7ace616e> in <module>
----> 1 estimatedPropagatorArray = estimator.estimate()

JavaError: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
org.orekit.errors.OrekitException: Jacobian matrix for type KEPLERIAN is singular with current orbit
	at org.orekit.propagation.numerical.NumericalPropagator.tolerances(NumericalPropagator.java:710)
	at org.orekit.propagation.numerical.NumericalPropagator.tolerances(NumericalPropagator.java:654)
	at org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder.buildIntegrator(DormandPrince853IntegratorBuilder.java:55)
	at org.orekit.propagation.conversion.NumericalPropagatorBuilder.buildPropagator(NumericalPropagatorBuilder.java:209)
	at org.orekit.propagation.conversion.NumericalPropagatorBuilder.buildPropagator(NumericalPropagatorBuilder.java:47)
	at org.orekit.estimation.leastsquares.BatchLSModel.createPropagators(BatchLSModel.java:323)
	at org.orekit.estimation.leastsquares.BatchLSModel.value(BatchLSModel.java:224)
	at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
	at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:616)
	at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:399)
	at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:436)

Impossible to resolve the system (while the IOD was giving nearly perfect results). Note, I also tried to add perturbation forces to the propagator, but the result remains the same. Is there anything I am doing wrongly?

small_arcs.ipynb (15.7 KB)

Indeed, the problem is difficult to estimate. As you said, observation arc is very small (i.e., 30sec) and you have a small number of measurements (i.e., 6). As a result, IOD is maybe the best solution, but orbit determination can work. Optimizing keplerian parameters increases the difficulty. Since the optimizer changes the values of the parameters to optimize the problem, sometimes the keplerian parameters can have singular values (i.e., inclination or eccentricity close to 0.0). Could you try working in cartesian elements? When initializing the propagatorBuilder, could you convert the initialOrbit in cartesian orbit by using OrbitType.CARTESIAN.convertType(initialOrbit).

Yes, the IOD result is very good and can be used as an initial guess of the orbit determination (by converting in cartesian elements following my remark above).

Best regards,
Bryan

Thanks @bcazabonne, by using cartesian elements, the estimation step does not fail anymore, and the results make sense:

# get an estimate from the IOD
# Keplerian parameters: {
# a: 4.227102759246746E7; e: 0.0014381077442844517; i: 14.820109677796811; 
# pa: 254.52333056539018; raan: -6.155828532482599; v: -105.28402914694347;}

# switch from keplerian to cartesian
initialCartOrbit = OrbitType.CARTESIAN.convertType(initialOrbit)

# set up the propagator and the estimator as above

estimatedPropagatorArray = estimator.estimate() # works now

estimatedPropagator = estimatedPropagatorArray[0]
estimatedInitialState = estimatedPropagator.getInitialState()
estimatedOrbit_init = estimatedInitialState.getOrbit()

OrbitType.KEPLERIAN.convertType(estimatedInitialState.getOrbit())
# <Orbit: Keplerian parameters: 
# {a: 4.228037968321704E7; e: 4.1469834187781324E-4; i: 14.825282314389423; 
# pa: 40.34433277766516; raan: -6.167512678491653; v: 108.88171979828613;}>

Thanks for your patience and help!

Perfect! :slight_smile: