Question on State Transition Matrices & Covariance


#1

Hello everyone!

I have some questions on consistency. I am producing a program in orekit which simulates the process of producing and correlating uncorrelated tracks.

I simulate a satellite orbit and take data when it passes though a radar fov. These “tracks” are on the order of 2-5 minutes long and the radar records the range and angular position every ten seconds while in the fov. The data points have random error associated with them ( 0.015 degrees, 40 m of std dev).

The first, last and middle radar data for each track is then plugged into a herrick-gibbs script I wrote to produce an initial orbital guess. That orbit is then cast into the EquinoctialOrbit type, as that suits well for bls.

To refine the guess I use the orekit batch-leasts-quares estimation package. I take that initial herrick-gibbs equinoctial orbit and use it as the initial guess. The rest of the tracks range and angular measurements are then fed into the bls estimator and a refined orbital guess is produced.

For each track from the same satellite, one would expect the bls estimator to produce relatively similar orbits. To check, I am using the mahalanobis distance like in this paper (https://link.springer.com/article/10.1007/s40295-013-0018-1)

To compare two tracks and see if they correlate, I take the covariance and state of the first occurring track and propagate it to the epoch of the second track. This is where my questions lie and where I believe I am having issues.

The propagation of the the epoch 1 state to the second is not an issue. The state vector I use is (a, ex,ey,hx,hy,lv) since those are the orekit elements. To propagate the covariance, however, I need the state transition matrix which I get from the propagation of epoch 1 to epoch 2 shown in the following code block:

  	NumericalPropagatorBuilder numPropBuild = getPropagatorEq(orb1.getOrbit());
  	double[] normalizedParameters = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  	NumericalPropagator numProp = numPropBuild.buildPropagator(normalizedParameters);
  	PartialDerivativesEquations PDE = new PartialDerivativesEquations("dYdY0", numProp);
  	orb1 = PDE.setInitialJacobians(orb1);
  	numProp.setInitialState(orb1);

  	SpacecraftState orbf = (numProp.propagate(dt2));
  	
      Array2DRowRealMatrix dYdY0 = new Array2DRowRealMatrix(6, 6);
      PDE.getMapper().getStateJacobian(orbf, dYdY0.getDataRef());
      //PDE.getMapper().getParametersJacobian(orbf, dYdY0.getDataRef());
  	System.out.println("\n\n Final Jacobian: " + dYdY0);

Then I get the covariance using estimator.getPhysicalCovariance

Then the mahalanobis distance is calculated in the following manor:

Pm(tn) = stm(tn,tm) *Pm( tm) *stm(tn,tm)^T
Pmn(tn) = Pn(tn) + Pm(tn)
z(tn) = xn(tn) - xm(tn), tn>tm
(km,n)^2 = z^T *Pmn(tn)^(-1) *z(tn)

where km,n is the mahalanobis distance. The more related the track the lower this value. tn is the later epoch and tm the earlier. x is the state and n or m mean which track they belong to, and finally stm is the state transition matrix.

Unfortunately, the output I am getting is wrong by orders of magnitude, and points towards an error. I believe that the error arises from the state transition matrix, the covariance or both. From browsing the archived mailing list it seems to me that the stm I’m producing likely might be the cartesian state transition matrix. Additionally, I’ve found both the optimum covariance and the physicalcovariance. I get the feeling the physical covariance may also be for cartesian elements.

Are my hunches correct or is something else happening?

Here I’ve posted the rest of my code, if anyone is interested or thinks something deeper is happening: https://pastebin.com/YfC2v0Cd

Thanks anyone for their help and I hope everyone has a great week!

Sincerely,
Paul


#2

Hi Paul,

Last time I looked in to this only Cartesian elements are supported when computing the STM.[1] Bug report is at [2]. Contributions are welcome.

The linked thread suggests a couple different methods, either using the d(Cartesian)/d(Equinoctial) Jacobian on both ends of the propagation, or using a FieldNumericalPropagator with DerivativeStructure.

[1] https://www.orekit.org/mailing-list-archives/orekit-developers/msg00114.html
[2] https://gitlab.orekit.org/orekit/orekit/issues/351


#3

@evan.ward Ah I think I understand from your links, but I still think I need some help. From the earlier orbit propagated to the later orbit I should use something like:
> double[][] deqdcart = new Array2DRowRealMatrix(6, 6).getData();

  	orbf.getOrbit().getJacobianWrtCartesian(PositionAngle.MEAN, deqdcart);

where orbf is the earlier orbit propegated to the second epoch.
This gives me the d(Eqiu)/d(cart) at later epoch. Then I take the dYd0 after the

line and multiply to get
stm = (d(eq)/d(cart))dYdY0(d(eq)/d(cart))^-1

correct? That’s what I gather from the links but you say

What do you mean by both ends?

Do you mean I should also get d(eq)/d(cart) before propagation like the following:

        double[][] deqdcart = new Array2DRowRealMatrix(6, 6).getData();
		orbi.getOrbit().getJacobianWrtCartesian(PositionAngle.MEAN, deqdcart);

where orbi is the orbit before being propagated to epoch 2?

Sincerely,
Paul


#4

I’ve tried what I described above and i’m still experiencing mahalanobis distances which is enormous. Here is the new calculation area:

NumericalPropagatorBuilder numPropBuild = getPropagatorEq(orb1.getOrbit());
double[] normalizedParameters = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
NumericalPropagator numProp = numPropBuild.buildPropagator(normalizedParameters);
PartialDerivativesEquations PDE = new PartialDerivativesEquations(“dYdY0”, numProp);
orb1 = PDE.setInitialJacobians(orb1);
numProp.setInitialState(orb1);

  	SpacecraftState orbf = (numProp.propagate(dt2));
  	
      Array2DRowRealMatrix dYdY0 = new Array2DRowRealMatrix(6, 6);
      PDE.getMapper().getStateJacobian(orbf, dYdY0.getDataRef());
      //PDE.getMapper().getParametersJacobian(orbf, dYdY0.getDataRef());
  	System.out.println("\n\n Final Jacobian: " + dYdY0);
  	System.out.println("Please be equinoctial:"+orbf.getOrbit().getType());
  	
      double[][] deqdcart = new Array2DRowRealMatrix(6, 6).getData();
  	orbf.getOrbit().getJacobianWrtCartesian(PositionAngle.MEAN, deqdcart);
  	RealMatrix deqdcartRM = new Array2DRowRealMatrix(deqdcart);
  	RealMatrix idedc = MatrixUtils.inverse(deqdcartRM); 
      RealMatrix stmf = deqdcartRM.multiply(dYdY0).multiply(idedc);
      RealMatrix stmi = MatrixUtils.inverse(stmf);

With this regime I get some NaN outputs and outrageously large kmn number.

My secondary hunch was that .getPhysicalCovariance is not giving me the correct number. When I use the optimum.getCovariance the output is much closer, but I’m not just going to change it because it fits expected behavior rather than fits mathematics.

Thanks again for all your help. If you know any additional info let me know!

Sincerely,
Paul


#5

The initial and final Jacobians are different because they depend on the orbital elements. You need to compute:

d(eq_0)/d(cart_0) * d(cart_0)/d(cart_1) * d(cart_1)/d(eq_1)

The subscripts indicate different times. As you can see the first and last terms are separate quantities. So you need to use orb1 to compute the first term instead of orbf.


#6

Ah thanks so much I’ll try this out!

Sincerely,
Paul


#7

Sorry one more question is the d(cart_0) in reference to the initial jacobian dYdy0 and would be identity, and d(cart_1) the dYdY0 after the propegation?

Then this computation is equal to that state transition matrix?


#8

Hi Paul,

Perhaps a text book will explain it better than I can in a message board. Do you have Montenbruck & Gill? It’s a good book that explains the concepts in chapter 7. It’s essentially just applying the chain rule a few times. The middle term is the Cartesian STM. The first and last terms are the Jacobians of the coordinate transformations you do to convert to Cartesian and from Cartesian, respectively.

Regards,
Evan


#9

I just found a copy of that book! I think you also did a good job of explaining it, especially after this last message.

I’m still getting vastly different answers as to which covariance I use between the optimum.getCovariance and the estimator.getPhysicalCovariance. Looking through the declaration I see that the physicalCovariance uses the get******ParametersDrivers.getDrivers to get the scales which it uses to “unormalize” the matrix.

Why is the matrix normalized in the first place?
What are those drivers?

Thanks again for all your help thus far it has been extraordinarily useful and helpful!


#10

Your welcome. Normalization is designed to avoid numerical issues when computing the state updates. It is easier for the algorithm if all quantities are approximately the same magnitude, but in physical units position in meters is many orders of magnitude larger than a drag coefficient, for example. So the the drivers translate back and forth between the normalized values that the numeric algorithms like to see and the physical values that astrodynamicists like to see. To see the difference try setting the referenceValue to zero, scale to 1 and min/max value to -/+ infinity before running the orbit determination.