Hi All,
I am performig orbit determination following the example of GorgiAstro
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,
Bryan
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
Alessandro
This is an excerpt of the code
set estimator
matrix_decomposer = QRDecomposer(1e-11)
optimizer = GaussNewtonOptimizer(matrix_decomposer, False)estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)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()
estimatedPropagator.resetInitialState(estimatedInitialState)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,
Bryan
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).
Cheers,
Maxime
Perfect, thank you very much