Kalman Filter Covariance Propagation in RIC


I am currently running the Kalman Filter on satellite in low earth orbit and using the Python Wrapper. For this satellite, it would greatly simplify things for me to have the covariance and process noise matrices represented in a local RIC (or TNW) frame. From what I understand, since my initial state and filter output is in the J2000 frame (and kalman propagator is inheriting these in cartesian), both of the aforementioned matrices are represented in this frame.

My current idea for this is to have my getProcessNoiseMatrix function (instance of CovarianceMatrixProvider) take the input J2000 covariance, transform it to the RIC frame, perform the evolution similar to here on p. 5 CovarianceReference, then transform it back to J2000. This method of 2 transforms seems a bit clunky, and also makes any intermediate covariance propagation between measurements difficult as each step will have to do the same. Plus I would guess that this will significantly increase the filter processing time.


Hi @aerow610,

Yes, they have to when they are given to the filter. So you need to transform from RIC to J2000 or backwards at every step of the Kalman.
The transformation has to be done at some point since it is not possible to propagate the orbit and covariance in RIC frame.

The interface CovarianceMatrixProvider was made for the kind of purpose you’re looking for. You need to implement this interface with your own code and perform the transformation from RIC to inertial in there.

Which is what you propose in your second point. Regarding this, I think your idea is a good one.
There is actually something close to what you’re looking for, it is the class UnivariateProcessNoise. Here you can define your process noise matrix with time-dependent univariate functions (like polynomials in your case) in RIC (or TNW… what we call LOF for Local Orbital Frame in Orekit ) frame.

I don’t have much time right now so I cannot read the whole paper you put in reference.
But I think you will encounter two problems:

  1. The sigmas used to compute the QRIC matrix seems to come from the current covariance matrix. In class UnivariateProcessNoise the coefficients of the univariate functions are defined at initialization. It will require a change of interfaces to add a covariance matrix in input.
  2. The UnivariateProcessNoise class assumes the process noise in RIC is a diagonal matrix while in your paper it is obviously not the case.

Tell me if my assumptions regarding the issues you may face are right.
It could be interesting to improve our models so if you need the UnivariateProcessNoise to evolve please open a bug on the Gitlab forge of the Java version of Orekit.

Hi @aerow610,

Edit from my previous post.
I’ve read the paper you mention in your first post.
I found it very interesting and if you have time I will be happy to have a follow up from you on whether you managed or not to reproduce the process noise matrix presented in the paper.

From what I’ve understood the “problem 1” I gave in my previous post is not an issue since the sigmas used in QRIC are computed by a post-processing, apart from the Kalman filter.

You’ll still get the “problem 2” of not being able to define the cross-coupling terms in the process noise matrix QRIC.
If you want this to be corrected for version 10.1 of Orekit, please open an issue on the tracker of the forge.


Hi @MaximeJ,

Thank you for all the info. I have just gotten back into this, so haven’t had time to try out the UnivariateProcessNoise class yet. I imagine that it will be fine using only the diagonal elements since I am working with less accurate orbit determination measurements and don’t need to be too precise at this point. Will let you know if I get it working.

I did however attempt to do the transform for each step, though no longer seem to have a converging filter as a result. Below is a code snippet of my implementation:

class filt_covar_provider(PythonCovarianceMatrixProvider):
    def getInitialCovarianceMatrix(self, initial):
        i_cov = JArray_double(len(init_cov))
        # Initial diagonal covariance into Java
        for x,i in enumerate(init_cov):
            i_cov[x] = i
        return MatrixUtils.createRealDiagonalMatrix(i_cov)
    def getProcessNoiseMatrix(self, prev, curr):
        t_meas = abs(curr.getDate().durationFrom(prev.getDate()))
        if t_meas > 0.0:
            meas_delta_t = t_meas
            t_meas = meas_delta_t
        return xform_vnc_to_j2k(curr, get_proc_noise_mat(t_meas))

vnc_acc_sig = [(2.4e-11)**2, (1.4e-9)**2, (5e-9)**2] # final values from AURA in paper

# Custom Process Noise Matrix
def get_proc_noise_mat(t):
    Q = Array2DRowRealMatrix(6, 6)
    for i in range(3):
        Q.setEntry(i, i, (vnc_acc_sig[i] * t**4 / 3))
        Q.setEntry(i, i+3, (vnc_acc_sig[i] * t**3 / 2))
    for i in range(3,6):
        Q.setEntry(i, i, (vnc_acc_sig[i-3] * t**2))
        Q.setEntry(i, i-3, (vnc_acc_sig[i-3] * t**3 / 2))
    return RealMatrix.cast_(Q)

# Transform the covariance to J2000
def xform_vnc_to_j2k(curr_state, cov_mat_vnc):
    # Get dummy propagator to build VNC frame
    vnc = LocalOrbitalFrame(j2k_frame, LOFType.VNC, propagator_num, 'VNC')
    sc_state = propagator_num.getInitialState()
    # Get frozen orbit and Jacobian
    vnc_to_j2k_frozen = vnc.getTransformTo(j2k_frame, sc_state.date).freeze()
    jacobianDoubleArray = JArray_double2D(6, 6)
    vnc_to_j2k_frozen.getJacobian(CartesianDerivativesFilter.USE_PV, jacobianDoubleArray)
    jacobian = Array2DRowRealMatrix(jacobianDoubleArray)

    cov_mat_j2k_java = jacobian.multiply(cov_mat_vnc.multiply(cov_mat_vnc.transpose()))
    return cov_mat_j2k_java

To start off, I just took the sigmas from Aura listed in the paper (can calibrate later once I get it working).

When I run the filter, it diverges fairly quickly, so I am afraid that my implementation of the transform was not done properly.

Hi @aerow610,

I think there’s an error in your code, right before the last line it should be:

cov_mat_j2k_java = jacobian.multiply(cov_mat_vnc.multiply(jacobian.transpose()))

And I’ll add two remarks:

  1. Not sure if it works in Python but most of the code of xform_vnc_to_j2k function where you get the Jacobian could be replaced with the line (assuming your current state is already in J2000):
LOFType.VNC.transformFromInertial(curr_state.getDate(), curr_state.getPVCoordinates()).getInverse().getJacobian(CartesianDerivativesFilter.USE_PV, jacobianDoubleArray );
  1. I think the sigmas you input may not be the good ones. You did the transform from RIC to VNC by assuming VNC = IRC; while I think it should be VNC = ICR.
    Nevertheless this does not explain the divergence.