ODTS for GPS satellite

[JordiS]

Hello all,

I am trying to use Orekitv10.1 for determining the Orbit and Clock offset of a GPS satellite. As an input I have Rinex files from 15 (not all of them have visibility of the satellite) different stations.

So far I have succeeded in parsing the Rinex files, extracting the Range, and adding modifiers (Tropo, and Iono) to the observations. My satellite and station configuration file is set-up as follows:

• Orbit estimation (set-up as a driver)
• Satellite clock offset estimated (set-up as a driver)
• Holmes Featherstone attraction model
• Station clock offsets estimated (set-up as a driver)

I also have build a propagator based on the Dormand-Prince propagator and build my Kalman filter.

The first problem I’ve found is that I have little idea of the process noise matrix values for the clock parameters? Do you have some guesses regarding these values for the clock ?

My second question is regarding the algorithm to propagate the clock states for the Kalman filtering. How is this done ? When I am initializing the covariance matrix I am using what is provided in the test AbstractDetermination.java file and that is synthetized in the following lines:

// Orbital covariance matrix initialization // Jacobian of the orbital parameters w/r to Cartesian final double[][] dYdC = new double[6][6]; initialGuess.getJacobianWrtCartesian(propagatorBuilder.getPositionAngle(), dYdC); final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC); RealMatrix orbitalP = Jac.multiply(cartesianOrbitalP.multiply(Jac.transpose()));

Should I also compute the Jacobian with respect the parameters ?

Another question is regarding the estimated parameters. Before going to the kalman filter I parse all Rinex files and multiplex theme according to their epoch as per the example above. Then I insert the all drivers. I was wondering that when the Kalman filter is executed at every iteration is trying to estimate all the drivers or only those that are observed ? When a new station is visible how is this included into the Kalman ?

My final point is regarding the performance of the Kalman filtering. I am executing 43200 measurement. The first 10000 are fairly quick to process but afterwards it slows down and it takes more than one day to complete. I’ve seen in the forum that you had identified a bug regarding this point which you solved by adding the NewtonianAttraction model in the force models. I am using version 10.1 (with the python wrapper) which supposedly resolves this problem. Could it be that the problem is not fully resolved ?

Thanks in advance for your help.

Jordi

[MaximeJ]

Hi @JordiS,

What you’re doing looks like a large and advanced usage of Orekit!
Congrats!

Sorry I can’t help you on your first question regarding how to compute the process noise matrix for the clock offset. I hope someone else here can.

@JordiS

Satellite clock offset estimated (set-up as a driver)
Station clock offsets estimated (set-up as a driver)

I’m afraid if you’re estimating all the clock offsets you might have an undetermination since the Kalman will not be able to decide where to put all the clock offsets.
But again I may be wrong.

@JordiS

Should I also compute the Jacobian with respect the parameters ?

I don’t think so, but you need to initialize the covariance matrix for all the parameters though. For propagation or measurements parameters it is usually a diagonal matrix
with the square values of the initial errors of your parameters. These initial errors you will have to guess though, it depends on your simulation and knowledge of the parameters at the start of your run.

@JordiS

I was wondering that when the Kalman filter is executed at every iteration is trying to estimate all the drivers or only those that are observed ?

Only those that are observed in your current measurements.

@JordiS

When a new station is visible how is this included into the Kalman ?

The new station’s clock offset should be a modifier of the measurement. Thus the derivatives of the measurement with respect to this clock offset will be computed and will populate the measurement matrix (usually called H) in the Kalman.

@JordiS

My final point is regarding the performance of the Kalman filtering.
[…]
Could it be that the problem is not fully resolved ?

It does look like it indeed
In Java we use some tools to evaluate the performances and the different calls to functions. I don’t know which equivalent tool could be used in Python.
If you could share your code we could translate it in Java and see what we get.
Some of us have been wanting to test GPS OD with larger data sets for a while but we lack time…

I hope some of my answers will help you.

Maxime

[JordiS]

Hello Maxime,

Thanks for your reply.

@MaximeJ

I’m afraid if you’re estimating all the clock offsets you might have an undetermination since the Kalman will not be able to decide where to put all the clock offsets.
But again I may be wrong.

If I don’t estimate then how is supposed Orekit to decorrelate the station clock from the Orbit ? If I remove it from the estimation parameters, the measurement will be corrected for the station bias ? In that case when I initialise the ground station I am obliged to set up the getClockOffsetDriver().setValue() to a very accurate value.

By the way what is the difference between setValue and setReference ?

[quote=“MaximeJ, post:2, topic:775”]
I don’t think so, but you need to initialize the covariance matrix for all the parameters though. For propagation or measurements parameters it is usually a diagonal matrix
with the square values of the initial errors of your parameters. These initial errors you will have to guess though, it depends on your simulation and knowledge of the parameters at the start of your run.
[/quote

Here I put you some code I’ve wrote regarding this matter:

Process Noise matrix Orbit part

cartesianOrbitalQ = MatrixUtils.createRealDiagonalMatrix([1.0e-4, 1.0e-4, 1.0e-4, 1.0e-10, 1.0e-10, 1.0e-10])

Coovariance Matrix orbit part

cartesianOrbitalP = MatrixUtils.createRealDiagonalMatrix([1.0e4, 4e3, 1.0, 5e-3, 6e-5, 1.0e-4])

Propagation parameters part

propagationP = MatrixUtils.createRealDiagonalMatrix([0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.001])
propagationQ = MatrixUtils.createRealDiagonalMatrix([1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-12])

Coovariance of the measurements per each station

measurementP = MatrixUtils.createRealDiagonalMatrix([1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0,
1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0])

#measurementP = None

Process noise matrix of the measurement model

measQ = 1e-6*1e-6
measurementQ = MatrixUtils.createRealIdentityMatrix(measurementP.getRowDimension()).scalarMultiply(measQ)
measurementQ = None
measurementP = None

Jacobian of the orbital parameters w/r to Cartesian

dYdC = JArray_double2D(6,6)
propagatorBuilder.getPositionAngle()
initialOrbit.getJacobianWrtCartesian(propagatorBuilder.getPositionAngle(),dYdC)
Jac = MatrixUtils.createRealMatrix(dYdC)
orbitalP = Jac.multiply(cartesianOrbitalP.multiply(Jac.transpose()))

Orbital Process noise matrix

orbitalQ = Jac.multiply(cartesianOrbitalQ.multiply(Jac.transpose()))

Build the full covariance matrix and process noise matrix

nbPropag = propagationP.getRowDimension() if propagationP != None else 0
nbMeas = measurementP.getRowDimension() if measurementP != None else 0
print(nbPropag, nbMeas)
initialP = MatrixUtils.createRealMatrix(6 + nbPropag + nbMeas, 6 + nbPropag + nbMeas)

Q = MatrixUtils.createRealMatrix(6 + nbPropag + nbMeas, 6 + nbPropag+nbMeas)

Orbital part

initialP.setSubMatrix(orbitalP.getData(), 0, 0)
Q.setSubMatrix(orbitalQ.getData(), 0, 0)

Propagation part

if propagationP != None:
initialP.setSubMatrix(propagationP.getData(), 6, 6)
Q.setSubMatrix(propagationQ.getData(), 6, 6)

Measurement part

if measurementP != None:
initialP.setSubMatrix(measurementP.getData(), 6 + nbPropag, 6 + nbPropag)
Q.setSubMatrix(measurementQ.getData(), 6 + nbPropag, 6 + nbPropag)

Notice that the measurementP and Q I change it to empty because when I am calling the:

Configure propagator

kalman = KalmanEstimatorBuilder().addPropagationConfiguration(propagatorBuilder,
ConstantProcessNoise(initialP, Q)).estimatedMeasurementsParameters(estimatedMeasurementsParameters).build()

Was giving me some errors regarding the dimensions. If you know the problem I would appreciate it.

@MaximeJ

It does look like it indeed
In Java we use some tools to evaluate the performances and the different calls to functions. I don’t know which equivalent tool could be used in Python.
If you could share your code we could translate it in Java and see what we get.
Some of us have been wanting to test GPS OD with larger data sets for a while but we lack time…

The code is large. If you tell me which part to share I can copy paste it. The Kalman filter is being called as:

from org.orekit.estimation.sequential import KalmanObserver
from org.orekit.propagation.integration import AbstractIntegratedPropagator
from org.orekit.python import PythonKalmanObserver
import MyObserver

Add observer

kalman.setObserver(MyObserver.MyObserver())
#aux = ObservedMeasurement(2)

Process the list of measurements

finalMeas = ArrayList().of_(ObservedMeasurement)

for meas in multiplexedMeasurements:
finalMeas.add(meas)

estimated = kalman.processMeasurements(finalMeas)[0].getInitialState().getOrbit()

The error I get in this line is:

JavaError Traceback (most recent call last)
in
14 finalMeas.add(meas)
15
—> 16 estimated = kalman.processMeasurements(finalMeas)[0].getInitialState().getOrbit()

JavaError: <super: <class ‘JavaError’>, >
Java stacktrace:
org.orekit.errors.OrekitException: minimal step size (1.00E-03) reached, integration needs 9.26E-04
at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:494)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:414)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:397)
at org.orekit.estimation.sequential.KalmanModel.predictState(KalmanModel.java:1028)
at org.orekit.estimation.sequential.KalmanModel.getEvolution(KalmanModel.java:897)
at org.orekit.estimation.sequential.KalmanModel.getEvolution(KalmanModel.java:58)
at org.hipparchus.filtering.kalman.extended.ExtendedKalmanFilter.estimationStep(ExtendedKalmanFilter.java:56)
at org.orekit.estimation.sequential.KalmanEstimator.estimationStep(KalmanEstimator.java:218)
at org.orekit.estimation.sequential.KalmanEstimator.processMeasurements(KalmanEstimator.java:236)
Caused by: org.hipparchus.exception.MathIllegalArgumentException: minimal step size (1.00E-03) reached, integration needs 9.26E-04
at org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator.filterStep(AdaptiveStepsizeIntegrator.java:330)
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:281)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:468)

I know there is a lot of information and questions but I would appreciate at least some guidance whether I am in the good direction and what can I do to improve the process and reach convergence. My priority at this point is to compute the Orbit, then the satellite clock and then the stations clocks.

Thanks in advance
Jordi

[MaximeJ]

Hello Jordi,

@JordiS

If I don’t estimate then how is supposed Orekit to decorrelate the station clock from the Orbit ? If I remove it from the estimation parameters, the measurement will be corrected for the station bias ? In that case when I initialise the ground station I am obliged to set up the getClockOffsetDriver().setValue() to a very accurate value.

Well I’m not an expert in GNSS OD but I’m afraid if you select all the clock offsets (satellite + all stations) the optimization process may be under-determined.
Example: the process could add 10000s to the satellite offset and subtract 10000s to all stations’ clock offsets without changing the value of the range measurement.
But I’m saying this given my knowledge of range biases in measurements, maybe the clock offsets work differently and this would not cause any issue.

@JordiS

By the way what is the difference between setValue and setReference ?

Reference value is usually the initial value (for batch least-square OD) or predicted value (for Kalman filter). It is used to normalize the parameter drivers with:
value_{N} = (value - value_{ref}) / scale
Normalization is used to avoid numerical errors in matrix calculations.

@JordiS

Notice that the measurementP and Q I change it to empty because when I am calling the:

Configure propagator

kalman = KalmanEstimatorBuilder().addPropagationConfiguration(propagatorBuilder,
ConstantProcessNoise(initialP, Q)).estimatedMeasurementsParameters(estimatedMeasurementsParameters).build()

Was giving me some errors regarding the dimensions. If you know the problem I would appreciate it.

The size of your initial matrices must match the number of estimated measurements (orbital, propagation and measurements).
From what I’ve understood from your code you have 14 propagation parameters and 13 measurements parameters:

• What are the 14 propagation parameters ? Propagation parameters are parameters from the force models (thrust of maneuvers, drag coefficient etc.). From what you said in your first post you should have none. I am right ?
• Measurements parameters are the parameters linked to measurements. There we should find your 15 stations’ clock offsets plus the satellite clock offset (if you estimate them all).

@JordiS

The code is large. If you tell me which part to share I can copy paste it.

I understand.
We would actually need all the configuration (stations + measurements + initial orbit etc.) but it may be too large or too much work to convert to Java. We will try to come up to something on our side to test with a very large number of measurements. However I cannot tell you when this will be done.
Could you open up an issue on the Gitlab forge related to your performance problem so we don’t forget to look into it ?

@JordiS

org.orekit.errors.OrekitException: minimal step size (1.00E-03) reached, integration needs 9.26E-04
at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)

This means that somewhere you put a minimum step size of 1e-3 seconds for your integrator and it wants to get below this (9.26e-4) to be able to integrate to the next step.
Try with a smaller minimum step size (say 1e-4s to start with).

@JordiS

I know there is a lot of information and questions but I would appreciate at least some guidance whether I am in the good direction and what can I do to improve the process and reach convergence. My priority at this point is to compute the Orbit, then the satellite clock and then the stations clocks.

Ok then do it in 3 steps, 1. estimate the orbital parameters only (no propagation or measurements’ parameters), 2. add the satellite clock offset as a measurement parameter, 3. add one and then as many as you want of the stations’ clock offsets to the measurements parameters.

I hope I helped you at least a little.
Maxime

[JordiS]

Hello Maxime,

@MaximeJ

Well I’m not an expert in GNSS OD but I’m afraid if you select all the clock offsets (satellite + all stations) the optimization process may be under-determined.
Example: the process could add 10000s to the satellite offset and subtract 10000s to all stations’ clock offsets without changing the value of the range measurement.
But I’m saying this given my knowledge of range biases in measurements, maybe the clock offsets work differently and this would not cause any issue.

Yes you are right, if I do select all clock offsets, the system will be undetermined. I guess that one possible solution is to fix one station clock so the clock offsets are determined with respect this clock. Do you know if this is possible in Orekit ?

@MaximeJ

• What are the 14 propagation parameters ? Propagation parameters are parameters from the force models (thrust of maneuvers, drag coefficient etc.). From what you said in your first post you should have none. I am right ?
• Measurements parameters are the parameters linked to measurements. There we should find your 15 stations’ clock offsets plus the satellite clock offset (if you estimate them all).

Ok I understand. SO I think this part of the code is incorrect since I was setting the measurement parameters to the propagator parameters. For what you say my propagation Q and covariance should be 0 since no force other than the gravitational are added.

@MaximeJ

I understand.
We would actually need all the configuration (stations + measurements + initial orbit etc.) but it may be too large or too much work to convert to Java. We will try to come up to something on our side to test with a very large number of measurements. However I cannot tell you when this will be done.
Could you open up an issue on the Gitlab forge related to your performance problem so we don’t forget to look into it ?

Yes I can open an issue. Just a small update regarding this topic. I changed the minimal step and run again I the following error was given to me (which points from my point of view to some sort of memory issue in the Java side).

JavaError Traceback (most recent call last)
in
14 finalMeas.add(meas)
15
—> 16 estimated = kalman.processMeasurements(finalMeas)[0].getInitialState().getOrbit()

JavaError: <super: <class ‘JavaError’>, >
Java stacktrace:
java.lang.OutOfMemoryError: Java heap space
at org.hipparchus.analysis.differentiation.DerivativeStructure.(DerivativeStructure.java:91)
at org.hipparchus.analysis.differentiation.DSFactory.build(DSFactory.java:128)
at org.hipparchus.analysis.differentiation.DerivativeStructure.add(DerivativeStructure.java:329)
at org.hipparchus.analysis.differentiation.DerivativeStructure.add(DerivativeStructure.java:65)
at org.hipparchus.analysis.interpolation.FieldHermiteInterpolator.value(FieldHermiteInterpolator.java:155)
at org.orekit.frames.EOPHistory.interpolate(EOPHistory.java:762)
at org.orekit.frames.EOPHistory.getPoleCorrection(EOPHistory.java:494)
at org.orekit.frames.ITRFProvider.getTransform(ITRFProvider.java:100)
at org.orekit.frames.Frame.getTransformTo(Frame.java:291)
at org.orekit.estimation.measurements.GroundStation.getOffsetToInertial(GroundStation.java:532)
at org.orekit.estimation.measurements.GroundStation.getOffsetToInertial(GroundStation.java:481)
at org.orekit.estimation.measurements.Range.theoreticalEvaluation(Range.java:165)
at org.orekit.estimation.measurements.AbstractMeasurement.estimate(AbstractMeasurement.java:204)
at org.orekit.estimation.measurements.MultiplexedMeasurement.theoreticalEvaluation(MultiplexedMeasurement.java:143)
at org.orekit.estimation.measurements.AbstractMeasurement.estimate(AbstractMeasurement.java:204)
at org.orekit.estimation.sequential.KalmanModel.getEvolution(KalmanModel.java:911)
at org.orekit.estimation.sequential.KalmanModel.getEvolution(KalmanModel.java:58)
at org.hipparchus.filtering.kalman.extended.ExtendedKalmanFilter.estimationStep(ExtendedKalmanFilter.java:56)
at org.orekit.estimation.sequential.KalmanEstimator.estimationStep(KalmanEstimator.java:218)
at org.orekit.estimation.sequential.KalmanEstimator.processMeasurements(KalmanEstimator.java:236)

@MaximeJ

Ok then do it in 3 steps, 1. estimate the orbital parameters only (no propagation or measurements’ parameters), 2. add the satellite clock offset as a measurement parameter, 3. add one and then as many as you want of the stations’ clock offsets to the measurements parameters.

Ok I can do that and try, but I guess that since my clock offsets are no estimated and evolve over time the orbit error will be degrading won’t be ?

Thanks a lot

Jordi

[MaximeJ]

Hello JordiS,

@JordiS

I guess that one possible solution is to fix one station clock so the clock offsets are determined with respect this clock. Do you know if this is possible in Orekit ?

I think you just need to “not select” the ParameterDriver involved in that station clock offset and it will work.
Somewhere in your code you must have something like ParameterDriver.setSelected(true) to indicate the Kalman it must estimate the clock offsets. You just need to replace true with false for that station. Or do nothing since the default selected value is false.

@JordiS

For what you say my propagation Q and covariance should be 0 since no force other than the gravitational are added.

Yes I think so.

@JordiS

Yes I can open an issue.

I’ve seen the issue, thank you for this

@JordiS

Just a small update regarding this topic. I changed the minimal step and run again I the following error was given to me (which points from my point of view to some sort of memory issue in the Java side).

Ok…
That one looks tough. I think we will really need to get our hands into it to fix it.
And that will take time, sorry for this…
Just to know, do you use a lot of Rinex files ? And are they “big” files?
If not maybe you could put some configuration info in the issue, like:

• Configuration used for the Kalman: Initial orbit, mass, force model, propagator used etc.
• Stations used and their configs
• Rinex files used
• … Anything that could help us reproduce your bug and fix it faster

@JordiS

Ok I can do that and try, but I guess that since my clock offsets are no estimated and evolve over time the orbit error will be degrading won’t be ?

Yes indeed. The idea is to first have something that works properly (in terms of code) with a “bad” precision and then slowly look for precision by adding corrections and compare the results.
We would actually be very interested in the results !

Maxime