Transforming a 7x7 covariance matrix into a local orbital frame

Orekit Python Wrapper usage
1

Dear Orekit community,
I am currently doing orbit determination with Range and BistaticRange measurements and a BatchLSEstimator using the Orekit Python wrapper. As BistaticRange measurements are affected by a clock offset between the ground stations, I utilize the ClockOffsetDriver of the participating ground stations to include and estimate the clock offset.
I would love to look at the covariance matrix in a local orbital frame (LVLH) to learn more about how the clock offset influences the satellite position in this frame. So far I have been using the StateCovariance class to perform this transformation, but as far as I understand the implementation of this class, it only works on 6x6 matrices with the six elements being the kepler elements in my case. According to my current understanding I need to expand the Jacobians using partial derivatives of the coordinates with respect to the clock offset-parameter. I do however not know how to obtain these derivatives in Orekit.
Could you please point me in the right direction how I could implement this transformation?

Cheers,
Pascal

2

@pascal.sauer

  1. So if all you want to do is see the position covariance in the local frame, you don’t actually need to propagate the 7th parameter - taking the 6x6 top-left submatrix and converting it using StateCovariance is fine. However…
  2. The propagation of the 7th factor ClockOffsetDriver has the potential to affect the pos/vel covariance values, but these effects are not calculated in the default 6x6 propagation (something I’d like to fix in future versions). You CAN however pull and propagate them manually using the following code:
// Build the propagator using the NumericalPropagatorBuilder class
parameters = propagatorBuilder.getSelectedNormalizedParameters()
propagator = propagatorBuilder.buildPropagator(parameters)

// Set your parameter ClockOffsetDriver to True
int locationOfYourParameter = X
propagator.getAllForceModels().get(X).getParametersDrivers().get(Y).setSelected(True)

// Propagate the propagator to the covariance start time
propagator.propagate(covTime)

// Initialize the STM matrix harvester
harvester = propagator.setupMatricesComputation("stm", None, None)

// Propagate to next time step
state = propagator.propagate(nextTime)

# Pull STM(N|0)
dYdY0 = harvester.getStateTransitionMatrix(state)

# Pull "extra" STM for extra parameters
if haveExtraParameters
    dYdPp = harvester.getParametersJacobian(state)

stm_N_0 = MatrixUtils.createRealIdentityMatrix(numPropagatedValues)
for i in range(6) :
    for j in range(6) :
        stm_N_0.setEntry(i, j, dYdY0.getEntry(i, j))

        if numPropagatedValues > 6 :
            for j in range(6, numPropagatedValues) :
                stm_N_0.setEntry(i, j, dYdPp.getEntry(i, j-6))

# Update covariance manually
nextCov = stm_N_0.multiply(initCov).multiply(stm_N_0.transpose())

I pulled this from my own code where I had to propagate drag/solar as part of my covariance and modified it a bit to make it a touch more readable, so it might not be completely perfect for your own scenario. But it should give you the right idea.

1 Like
3

Hi!

I’ve something in mind for it. In order to support any type of parameters in the covariance, I would like to add a List<AdditionalCovarianceEquationProvider> in the covariance propagation of Orekit. Thanks to that, users could add additionnal elements to the covariance, and not only the 6 orbital elements. With this, all these elements could be stored in a single StateCovariance matrix.

Bryan

1 Like
4

Dear @baubin,
Thank you for the code. I have played with it a bit and it looks like this would solve the problem if I was estimating ParameterDrivers of the propagator, i.e force models. However the ClockOffsetDriver is a ParameterDriver affecting the measurements, which is in fact not affecting the state transition matrix or the getParametersJacobian method. As far as I can see, the derivatives for measurement ParameterDrivers are not exposed to the user.

Cheers,
Pascal

5

Sigh… sorry @pascal.sauer I didn’t specifically realize that ClockOffsetDriver was a measurement parameter and that it wouldn’t behave the same way as the propagation parameters. @bcazabonne I know we need to add the measurement/propagation parameters to the general covariance propagation in the long term, but in the short term do you know a workaround to solve Pascal’s problem now?

6

Hello everyone,
I might have found a solution to my problem: I manually re-implemented the transformations of the StateCovariance method changeCovarianceFrame and expanded the transformation matrices with an identity row/column (1.0 on the diagonal, 0.0 everywhere else), as a measurement bias should not impact the frame transformations. Is this assumption correct? I also multiplied the covariances in the last column and row with the speed of light to get every unit in meters, and the resulting covariance matrix seems reasonable.

# %% manual frame change for 7x7
params = Initial_parameters("scenario.toml")
lse = params["LSE"]
ini_date = datetime_to_absolutedate(params["Initial_orbit"]["epoch_date"])
lse["delta_t_1_is_estimated"] = True

speed_of_light = 299792458.0

prop, estimator = run_simulation(params, "test2")
cov_matrix = estimator.getPhysicalCovariances(1e-20)
orbit = prop.getInitialState().getOrbit()

# changeCovarianceFrame
## changeFrameAndCreate(orbit, epoch, frame, frameOut, orbitalCovariance, orbitType, angleType);

lof_type = LOFType.LVLH
lof = LOF.cast_(lof_type)
eci_type = orbit.getType()
angle_type = PositionAngleType.MEAN
frame_in = orbit.getFrame()
frame_out = lof

### changeTypeAndCreate
array = JArray_double2D(7, 7)
cart_orbit = OrbitType.CARTESIAN.convertType(orbit)
cart_orbit.getJacobianWrtCartesian(angle_type, array)
# Compute dOutputdCartesian
dOdC = Array2DRowRealMatrix(array, False)
dOdC.setEntry(6, 6, 1.0) # Set the last diagonal element

array = JArray_double2D(7, 7)
cart_orbit = OrbitType.KEPLERIAN.convertType(orbit)
cart_orbit.getJacobianWrtParameters(angle_type, array)
# Compute dCartesiandInput
dCdI = Array2DRowRealMatrix(array, False)
dCdI.setEntry(6, 6, 1.0)

# Compute dOutputdInput
dOdI = dOdC.multiply(dCdI)
# Conversion of the state covariance matrix in target orbit type
cartesian_cov = dOdI.multiply(cov_matrix.multiplyTransposed(dOdI))

## changeFrameAndCreate
# Get the Cartesian covariance matrix converted to frameOut
CtL_transform = KinematicTransform.cast_(frame_out.transformFromInertial(ini_date,
                                                                         orbit.getPVCoordinates(frame_in)))
jacobianFromFrameInToLofOut = MatrixUtils.createRealMatrix(CtL_transform.getPVJacobian())
jacobian_matrix = Array2DRowRealMatrix(JArray_double2D(7, 7), False)
for i in range(7):
    for j in range(7):
        if i==6 and j==6:
            jacobian_matrix.setEntry(i, j, speed_of_light)
        elif i==6 or j==6:
            jacobian_matrix.setEntry(i, j, 0.0)
        else:
            jacobian_matrix.setEntry(i, j, jacobianFromFrameInToLofOut.getEntry(i, j))

covOut = jacobian_matrix.multiply(cartesian_cov.multiplyTransposed(jacobian_matrix))
1 Like