Batch orbit estimation covariance matrix

Hi All,
I am performig orbit determination following the example of GorgiAstro

laser-orbit-determination/02-orbit-determination-example.ipynb at 6cafcef83dbc03a61d64417d0aeb0977caf0e064 · GorgiAstro/laser-orbit-determination · GitHub

I would like to ask how “to estract” the covariace matrix associated to the OD for ALL the filter estimation time intervall.

I mean I want to plot the error with respect to the ±3-sigma intervall, so I need the time depended covariace matrix associated to the orbit estimation.

In the example, for me, is not clear at which time the covariance is computed

Many thanks for your time and considerations
Kind Regards

Hi @AlessandroV

What type of estimator do you use? Recursive filter or least squares algorithms?

In both case, you can access the covariance matrix for each estimator step (i.e., each measurement for the recursive filter or each iteration for the least squares) by adding an observer to the estimator.
For the batch least squares the observer is represented by the BatchLSObserver class. For the recursive filter it is represented by the KalmanObserver class.
Both class is called by its corresponding estimator just after an estimation step. Therefore, it will provide fresh estimated parameters. Just not that both classes are interface, and Orekit doesn’t provide implementations. As a result, you will have to implement the interface.

An example of BatchLSObserver implementation in Python is available here: Custom BatchLSObserver with Python wrapper - #7 by Anthony-C31 The covariance information for the processed iteration is stored in the lspEval parameter. The covariance is available using lspEval.getCovariances(1.0e-10) (Note that you can access other iteration information like the RMS)

An example of KalmanObserver implementation in Python is available here: KalmanEstimator : KalmanObserver and negative eccentricity with NumericalPropagatorBuilder - #2 by aerow610 The covariance information for the processed measurement is stored in the estimation parameter. The covariance is available by using estimation.getPhysicalEstimatedCovarianceMatrix() (Note that you can access other parameters like the Kalman gain, the value of the estimated parameters, the predicted and corrected measurements, etc.)

Best regards,

Thank you, very much

Hi @bcazabonne , I have a doubt, that I would like to clarify:

With the Batch-type filters, at what state estimated, does the covariance matrix refer to?
I mean the covariance matrix that I get with:

covMat_eci_java = estimator.getPhysicalCovariances(1.0e-10)

I mean this: let’s say I have this scenario:
t1 = start of measurements and start of OD
t2 = end of measurements (but NOT end of OD propagation)
t3 = end of OD propagation

What point does the covacity matrix refer to?
t1, t2, t3, an average of all the points, … ?

Many thanks for your time and help
Kind Regards

This is an excerpt of the code

set estimator

matrix_decomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrix_decomposer, False)

estimator = BatchLSEstimator(optimizer, propagatorBuilder)

observableSatellite = ObservableSatellite(0) # Propagator index = 0

Performing the orbit determination

estimatedPropagatorArray = estimator.estimate()

date_start = datetime_to_absolutedate(posMeas_df_pack.index[0]) # start meas and OD
date_end = datetime_to_absolutedate(posMeas_df_pack.index[-1] + timedelta(seconds= timePropAfterMeas )) end of OD propagation

#Note: datetime_to_absolutedate(posMeas_df_pack.index[-1]) is the time of end meas

First propagating in ephemeris mode

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

generator = estimatedPropagator.getEphemerisGenerator();
estimatedPropagator.propagate(date_start, date_end);
bounded_propagator = generator.getGeneratedEphemeris();

Getting covariance matrix in ECI frame

covMat_eci_java = estimator.getPhysicalCovariances(1.0e-10)

Hi @AlessandroV

The epoch of the covariance given by Orekit batch least squares estimator shall correspond to the epoch of the initial guess of the orbit determination.
Because the purpose of the least squares estimator is to optimize an initial guess according to a set of observed measurements, both the estimated state and the covariance are given at the initial guess epoch.

Best regards,

Perfect, thank you very much

Hi, @bcazabonne ,sorry again,
How can I get the the estimated state , with respect to the orbit guess ? (using Python wrapper)
I am trying with estimator.getLastEstimations() , estimator.getOrbitalParametersDrivers and estimator.getOptimum , but I did not understand how get the array of cartesian/keplerian of the solution (with respect to the initial guess).

Can you help me ?

Many thanks for your time and help
Kind Regards

Hi @AlessandroV,

Method BatchLSEstimator.estimate() returns an array of propagators initialized with the estimated states.
So if you have only one satellite to estimate you will get the estimated orbit with (in Java):

Propagator[] estimatedPropagators = estimator.estimate();
Orbit estimatedOrbit = estimatedPropagators[0].getInitialState().getOrbit();

Here’s an example for using the BatchLSEstimator in Java: AbstractOrbitDetermination, from the tutorials.
And in Python: Python SLR orbit determination example (courtesy of Clément Jonglez).


1 Like

Perfect, thank you very much