Interesting behavior for Batch Least Squares

Orekit usage
1

Hello Everyone,

I’ve found some strange behavior. The batch least squares is not converging in a I would expect. First let me describe my methodology. I first take an orbit (I’ll refer to this as the true orbit) and propagate it forward. I have a groundstation mimicking a radar. When the true orbit encounters the edge of the radar, the propagation stops takes angular and range observations and propagates forward ten seconds, stops, and takes another set of observations until the true orbit exits the field of view of the groundstation. I then take the first, middle and last observations to produce an orbital estimate using herrick gibbs. I then use a LevenbergMarquardt least squares optimizer to run batch least squares with the herrick gibbs guess as an initial guess and all the observations as input. My comparisons from the BLS converged orbit and the true orbit is done using Mahalanobis Distance (Basically standard deviation distance).

Well I’ve found some strange behavior. First using perfect, no error observations does not cause the BLS to converge on the correct answer, even if I use the truth as oppose to the herrick gibbs as the input. That is strange issue number one.

Secondly, the more observations I feed into the least squares, the greater the number of deviations the converged orbit is from the true orbit. This is counter intuitive.

Finally the last if I take the covariance and propagate it forward, and compare it to the true orbit once again, a weird amplifying sinusoidal error. Check it outimage

The x axis is days and y axis mahalanobis distance (lower is better). The lowest line is 14 observations, the one above 15 and the highest is 19 observations. This graph was done using perfect no error observations.

Now from looking at the covariance matrices, it looks as if they are somewhat close to being singular, and also the main contribution seem to come from the eccentricity and semimajor axis error. In fact from the covariance matrix it seems they are highly correlated (0.9+)

I believe these errors are probably all stemming from the same issue. They are present when using all orbital elements, but seem to be worse in keplarian than equinoctial, and worst in Cartesian. My only theory is that the problem is from numerical inaccuracies stemming from the starting and stopping of the propagation of the true orbit.

Does anyone have any suggestions or input?

2

As I understand, you use simulated measurements.
How are these measurements computed? Are the simple direct computation of line of sight direction and norm of station to satellite vector or do they take into account all physical and practical details (signal trael time, motion of satellite and station during signal fligth, tropospheric and ionospheric effects…)? For perfect measurements to work properly in a later orbit determination, the same model must be used in measurement simulation and measurements use in OD. As measurement use in OD does take into account these time of flight and local motion effects, they must be included in the measurement simulation phase too. There is a measurement simultation feature in Orekit that is consistent with the measurements model used in orbit determination, did you use it?

1 Like
3

Hi Luc, thanks for your quick response. Let me try to address your suggestions as best I can!

For the range, I use the base frame of the groundstation, and use the .getRange function. I don’t think I take into account any physical or practical details .I simply just feed that double, and the time directly into the Range type, with no changes. So I understand that I should be adding the time duration of the signal flight to the date assigned to both the angular and range measurements?

I believe I have done this. I use the same propagator with the same settings.

Ah thank you for pointing me here. I did not use this. I assume I have to rebuild my observation generation using the MeasurementBuilder class.

Thanks again for all your help

4

For the range, I use the base frame of the groundstation, and use the .getRange function. I don’t think I take into account any physical or practical details .

No, this method from TopocentricFrame simply computes the norm of the vector from ground station to
external point at one date, ignoring travel time.

I simply just feed that double, and the time directly into the Range type, with no changes. So I understand that I should be adding the time duration of the signal flight to the date assigned to both the angular and range measurements?

Yes, taking into account that for a two-way measurement, there are three positions to consider, at three different times: signal transmission from ground, transition on spacecraft, signal reception on ground and that everything moves.

I believe I have done this. I use the same propagator with the same settings.

Fine, but measurement model should also be the same, otherwise the more realistic model used in OD will consider that the observed measurement is offset, because it was computed using a different model. Worse, the offset is not a noise that can be filtered out using a sufficient number of measurements, it is a slowly varying deterministic difference, and as such will be combined with the orbit, so it may explain weird differences.

Ah thank you for pointing me here. I did not use this. I assume I have to rebuild my observation generation using the MeasurementBuilder class.

Exactly. There is a RangeBuilder and an AngularAzElBuilder classes that can be useful for you. Note that if you add some EstimationModifier to your OD (say tropospheric and ionospheric effects, or biases), then you should also add them to the builders before generating the measurements.

As your measurement scheduling seems quite specific, you may need to implement your own DateSelector and use it together with an ElevationDetector to build an EventBasedScheduler so you can generate measurements according to your specifications.

The measurements generation feature is quite rich and complex. Look at the documentation here: https://www.orekit.org/static/architecture/estimation.html, there are some static class diagrams that may help sorting out the architecture.

2 Likes
5

Thanks Luc! Let me try those suggestions and Ill let you know what I find!

6

@luc
Hello, again’

I am using the estimation generation, but I have a couple questions.
First a question on usage: Using a observation builder (range or angularAzEl) I can’t seem to get more than one type of observation. I can produced the right number of, and accurate Range measurements or AngularAzEl measurements. The measurements give the correct answer when using bls orbit determination, but it seems like I only generate whichever type of measurement comes first. To fix this I’ve:
changed the visibility handler,
tried using distinct generators
tried using distinct propagators in the scheduler, and I’ve tried changing the propagatorIndex in the Builders.

Additionally I thought maybe it would work if I calculated in parallel, but I have not succeeded. I’m sure the issue is simple to fix and its just a boneheaded mistake. Let me know if the solution jumps out to you!

Additionally, how would I read in real data for observations. I guess obviously a file parser which would allow building range measurements, and then those measurements already take in the considerations of the real world which were causing me issues earlier. Is there some other considerations, or thoughts you would have on reading in real radar data?

Thanks so much!
Paul

7

The proper way to generate several types of measurements is to set up one Generator and one propagator only (adding several propagators is useful only for intersatellites measurements) and to add as many schedulers as there are measurements types. So if you need both range and angular measurements, you should add two schedulers to your generator. Once everything has been set up, you can call the generate method once and get all measurements.

If you did that and you only get one measurement type, then this is a bug because it is the way it is intended to be used.

Concerning measurements reading, you are right, you need a parser for the format that is generated by your stations. The tutorial use a very simple format, but is was custom made and is not a standard by any means, so you can reuse it if you want, but it will not be well suited for exchanging data with other systems. Orekit however does support reading some other standard file formats, like CCSDS TDM files and Rinex observation files. Beware the TDM is not really a complete specification, in the sense that it really needs data producer and consumer to set up an ICD to tune some details.

8

Thanks for the quick response! I unfortunately have tried using multiple schedulers on the same generator. They still only produce one measurement type. I’m using an older versions though, (the development version out before 9.3) so I’ll upgrade and see if it still is an issues. Also thank you for the information on reading in data!

9

Hey @luc

To follow up, I checked with the most recent version of orekit and after changing my code to fit the new constructors, the problem persisted:

I am only able to generate one type of observation, based on what is coded first. Here are some relevant snippets of my code:

final EventDetector sta1Visi = new GroundFieldOfViewDetector(sta1Frame, fov); //
AdapterDetector detector =new AdapterDetector(sta1Visi);
SpacecraftState s0 = new SpacecraftState(initialOrbit);
detector.init(s0, initialDate);

  AdapterDetector det1 = detector;
  AdapterDetector det2 = detector;		
  double[] azElError = new double[]{0.015*Math.PI/180,0.015*Math.PI/180};
  double[] baseweight =new double[]{1.0,1.0};

  double rangeSigma = 40.0;
  double rangeBW = 1;
  ObservableSatellite obs = new ObservableSatellite(0);
  RangeBuilder rB = new RangeBuilder(null, station, false, rangeBW, rangeBW,obs);
  AngularAzElBuilder aAEB = new AngularAzElBuilder(null, station, azElError, baseweight, obs);	
  double  timeToEnd = days*86400; // I think the conversion is 1s/1

  AbsoluteDate finalDate = initialDate.shiftedBy(timeToEnd);
  FixedStepSelector fss = new FixedStepSelector(10., TimeScalesFactory.getUTC());
  EventBasedScheduler<Range> eBS = new EventBasedScheduler<Range>(rB, fss,	 numProp, detector, SignSemantic.FEASIBLE_MEASUREMENT_WHEN_NEGATIVE);
  EventBasedScheduler<AngularAzEl> aeBS = new EventBasedScheduler<AngularAzEl>(aAEB, fss, numProp, detector, SignSemantic.FEASIBLE_MEASUREMENT_WHEN_NEGATIVE);		
  Generator genR = new Generator();
  genR.addPropagator(numProp);
  genR.addScheduler(aeBS);
  genR.addScheduler(eBS);

Again if I just switch the last two lines I get range measurements only, but in this configuration I only get angular measurements. The measurements do converge to the right orbit either was, and they are at the right time and the right number. Any ideas on fixing this?

Thanks,
Paul

10

It looks more and more like a bug to me.
One last check: could you try again but instead of reusing the same instance of detector to build two independent instances before passing them to the EventBasedScheduler constructors? Looking at the core, it should not change anything so I fear you will still get only one type of events. If you confirm this does not solve the problem, then please open a new issue at https://gitlab.orekit.org/orekit/orekit/issues. I consider this issue as important, so we will adress it very soon, before releasing upcoming 10.0.

11

I’ve just tried that and no luck! I’ll open up an issue now!

Thanks
Paul

12

Dear @luc

I believe I have found another problem. When using the generate function, the returned observations only come from the first pass of the satellite into observable space, and none of the subsequent. I’ve tried forcing the generator to restart after the end of first pass, but no dice. In testing I propagated for a time that I calculated to contain 3 passes, but only data from the first was generated. Just in case I made a mistake, I then increased the difference between the first date and second date to ridiculous lengths, and no change. Any suggestions on getting the data from all the passes?

Paul