Covariance matrix before last measurements

Hello to everyone !

I am currently working on the calculation of covariance matrix through orbit determination with Least square and Kalman filtering. In the context of space debris, I need to make orbit determination over a period of 5 or 10 days and to know the covariance matrix at different epoch and not only at the last state. For now, I am able to get the covariance matrix at the last measurement, I would like to know if it is possible to get the covariance matrix at a previous state after the orbit determination for the two methods.

I have other questions too, could someone explains to me why the standard deviations in the covariance matrix for positions and velocities are greater for a period of 5 days than for a period of 10 days. At first, I thought than making an orbit determination over a longer period would increase the incertitudes at the last measurement, but it seems not.

And a last one, why the overall standard deviations are much greater for Kalman than for LS ?

Thanks for your help !

Sacha

Hi @Sacha

Welcome to the Orekit forum!

Yes it is. For the Kalman Filter, you can easily access the covariance matrix, at each measurement epoch by setting a KalmanObserver to your KalmanEstimator using

kalmanEstimator.setObserver(new MyKalmanObserver());

Here an example for the implementation of MyKalmanObserver() class

public class MyKalmanObserver implements KalmanObserver {

    @Override
    public void evaluationPerformed(final KalmanEstimation estimation) {
        // Access the physical covariance matrix
        final RealMatrix covariance = estimation.getPhysicalEstimatedCovarianceMatrix();
        // Do what you want with the covariance
    }

}

Please note that the evaluationPerformed(final KalmanEstimation estimation) method is called each time a measurement is processed by the Kalman filter.

For the batch least squares covariance matrix, the method is a bit different. Indeed, the algorithm optimize the least squares problem using the whole package of measurements. Therefore, you can only access the covariance at the orbit determination epoch. However, it is possible to propagate this estimated covariance, at the end of the orbit determination process, using Orekit. Indeed the covariance at a given epoch can be expressed as a function of the covariance at the orbit determination epoch: P(t) = STM(t,t0) * P(t0) * transpose(STM(t,t0)). Where STM is the state transition matrix.
I will give you an example considering that only orbital parameters are estimated.
Let estimator by your BatchLSEstimator object, first you need to access the estimated orbit propagator and the estimated covariance using

final NumericalPropagator estimatedPropagator = estimator.estimate()[0];
final RealMatrix estimatedCovariance = estimator.getPhysicalCovariances(1.0e-15)

Then configure your estimated orbit propagator to propagate the covariance

// Additional equations name
final String equationName = "partial-derivatives";
// Initialize spacecraft state partial derivatives object
final PartialDerivativesEquations partials = new PartialDerivativesEquations(equationName, estimatedPropagator);
// Get initial state
final SpacecraftState state = estimatedPropagator.getInitialState();
// Add derivatives to initial state
final SpacecraftState stateWithDerivatives = partials.setInitialJacobians(state);
estimatedPropagator.resetInitialState(stateWithDerivatives);
// Initialise the Jacobians mapper
final JacobianMapper jacobiansMapper = partials.getMapper();

Finally, you need to create and add a step handler to your orbit propagator to perform the covariance matrix propagation.

// Set master mode and step handler
estimatedPropagator.setMasterMode(stepSize, new MyStepHandler(jacobiansMapper, estimatedCovariance));

Where stepSize is the step at which you want the computation of the propagated covariance matrix. Here an example of MyStepHandler class:

public class MyStepHandler implements OrekitFixedStepHandler {

    private JacobianMapper mapper;
    private RealMatrix estimatedCovariance;

    MyStepHandler(final JacobianMapper mapper, final RealMatrix estimatedCovariance) {
          this.mapper = mapper;
          this.estimatedCovariance = estimatedCovariance;
    }

    @Override
    public void handleStep(SpacecraftState currentState, boolean isLast) {

            // Current state transition matrix
            final double[][] aYY0 = new double[6][6];
            mapper.getStateJacobian(currentState, aYY0);
            final RealMatrix dYdY0 = new Array2DRowRealMatrix(aYY0, false);

            // Current covariance
            final RealMatrix currentCovariance = dYdY0.multiply(estimatedCovariance.multiplyTransposed(dYdY0))

            // Do what you want with the covariance

    }
}

Please note that I have just written this code snippet in the forum editor. Therefore, it may contain errors. Furthermore, I think that the given method can be applied to the Kalman estimator by replacing the estimator.getPhysicalCovariances(1.0e-15) instruction by kalmanestimator.getPhysicalEstimatedCovarianceMatrix()

I think that is because you have more measurements. Therefore, the estimation algorithms are able to provide better estimation of the parameters and therefore standard deviations are better for the 10 days fit span.

Good question. Algorithms are very differents. I’m not able to answer this question.

Best regards,
Bryan

1 Like

Hi Bryan,

Thank you very much for you answers ! It is exactly what I need. Although, I have still some difficulties to understand everything about covariance matrix.

In an other example, I have a SP vector at some date, this is my initial state. I want to know what is the covariance matrix 5 days later. I generate theoretical measurements by propagating the initial state until the date I want to have the covariance matrix, with a step of 1 min. Then, for the LS determination, I set up the propagator with the initial state. I set the estimator equal to the BatchLSEstimator(optimizer, propagator) with optimizer = LevenbergMarquardtOptimizer() and the propagator I’ve just initiated. Finally, I add the measurements to the estimator and I run the estimation.

So, is this the right way to make orbit determination to get covariance matrix ?
And at the end, I get the covariance matrix at the date of the initial state, that’s right ? But how it is calculated ?

With this example and with what you just said, if I want the covariance matrix at the last estimated states, I need to propagate the initial covariance, are you agree ?

Thank you again !

Sacha

I forgot to ask you something else. Is it wrong to set the propagator for the estimation with the last measurement and not with the initial state. If I do that, the estimator will return me the initial state as the last measurement and the covariance matrix at this date, so at the date I want. But I’m not convince by doing this way. I’m maybe wrong.

Thanks.

Sacha

At the end of the estimation process, you have a covariance matrix which corresponds to the initial state epoch (i.e. the initial state you used to initialize the orbit determination). If you want the covariance matrix at another epoch, I recommend to propagate the covariance at the needed epoch using my previous recommendations.

The covariance matrix is computed according the the following formula: P = (A^T W A)-1
Where W is the weighting matrix and A is the partial derivatives matrix (i.e. partial derivatives of the measurements with respect to the estimated parameters) and A^T its transpose.
If you want more information about the covariance matrix in a batch least squares orbit determination, I recommend you the paper: Covariance Transformations for Satellite Flight Dynamics Operations It is very useful to understand covariance matrices. It is available here.

it is better to use the same satellite state used to generate the measurements as initial state of the orbit determination.

Best regards,
Bryan

It is more clear for me now.

Thank you very much, it helped me a lot.

I’ll let you know if I have other questions.

Cheers,
Sacha

Here again sorry,
I’ve tried to do what you wrote for the covariance matrix propagation with LS method. I use python btw but I managed to translate your code. So after creating the step handler class, and after the step estimatedPropagator.setMasterMode(stepSize, new MyStepHandler(jacobiansMapper, estimatedCovariance)), I just have to make a estimatedPropagator.propagate(t) with t the date at which I want the covariance matrix ? And I’ll have the matrix at each stepSize ?

Thank you Bryan.

Sacha

Yes, you have just to call the propagate method.

Hello everyone !

I am also working on the calculation of covariance matrix through orbit determination with Least squares method. This discussion is very useful as I try to access the covariance at different epochs.
I followed your instructions @bcazabonne and tried to adapt it to the new version Orekit 12.0. Thus, I replaced PartialDerivativesEquations with EpochDerivativesEquations. However, the method getMapper() does not exist anymore with EpochDerivativesEquations. Therefore, I can’t initialise the Jacobians mapper as you did. I tried to find alternatives but couldn’t find one.

Do you have any clue how to do this with Orekit 12.0 ?

Thank you for your help !

Hugo

Hi @hugo and welcome to the Orekit forum!

Unfortunately, that’s not the correct way. I hope that one day we will remove the EpochDerivativesEquation class :sweat_smile:

You shall use the harvester from the numerical propagator, like that:

String stmAdditionalName = "stm" propagator.setupMatricesComputation(stmAdditionalName, null, null);

My recommendation is to follow the tutorial we developed. It is compatible with Orekit 12.0. It is available here:

Best regards,
Bryan

2 Likes

Hi @bcazabonne,

Thank you very much, the tutorial helps a lot !
I think I got it now :slight_smile:

Best regards,
Hugo