Batch LS wouldn't give correct result for OD...(Code and data attached)

Orekit usage
1

Hi experts, I’m new to Orekit.

I’m trying to develop a pipeline to do IOD from optical observations, then use all the thousands of data points for refining the Keplarian elements, then output usable TLEs. For the current stage, it’s a traditional single station, angle-only observation.

My testing camera is quite simple and cheap, gives me 8 arcsec resolution. With manual control, I can cover an arc of over 10 minutes for a single pass (> 120 deg distance on the sky). I suppose it is enough for the testing of the orbit determination pipeline and for generating an OD result with not too bad precision?

I searched the Orekit forum and found some example scripts to do orbit determination. Based on them, I produced my script, attached here: Dropbox - Forum - Simplify your life
(as you can see, my data is of quite high quality, with very low scattering, on the scale of one pixel.)

Now the script can do IodGooding and gave reasonable results. (I have a small function to estimate the range of the satellite. But with or without it, IodGooding always works to different extent.)

But in the batch LS part, it never succeeded.

I have tried different parameters, feeding different numbers of data points. (from just using a few data points, to a few dozens, to hundreds, till all of the 3000+ points.)

But it’s either a “true anomaly -109.956 out of hyperbolic range (e = 1, -3.142 < v < 3.142)” or a “Jacobian matrix for type KEPLERIAN is singular with current orbit”. (This particular satellite pass had a highest altitude of ~40 deg, so I’m not in the orbit plane, there shouldn’t be a singularity problem as I understood.)

I have tried to switch on/off the perturbations one by one but it does not affect the faulty result.

Is there anything obvious (or not so obvious) got wrong in my code?

Thank you for any help and comments!

==============
BTW, though I’m not using the Atmospheric Drag part of the code now, but when I tried to do so, there will be a bug report on this line:

msafe = MarshallSolarActivityFutureEstimation(‘DEFAULT_SUPPORTED_NAMES’,MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)

It seems maybe I’m writing ‘DEFAULT_SUPPORTED_NAMES’ in a wrong way. But I cannot find an example showing the correct way. Any suggestions on this?

===============
BTW, what do people normally do when the size and shape of a satellite is unknown but they still need to propagate the orbit? Any empirical formula? Maybe a correlation with brightness?

2

Hi @spto

I don’t see anything obvious.
Maybe look at yours dates, as julian day starts at Noon (i.e. 0s in the day is noon).
Another point is that your float variables should probably be double.

Concerning DEFAULT_SUPPORTED_NAMES, it is a predefined constant in the Java file MarshallSolarActivityFutureEstimation.java. It is a regular expression representing the MSAFE files taht can be loaded. It’s value is "\\p{Alpha}\\p{Lower}\\p{Lower}\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}(?:f|F)10(?:_prd)?\\.(?:txt|TXT)"

What people do when they don’t know the satellite is to guess. If they are looking at a cubesat or at a commercial heavy satellite, they adapt their input accordingly. In any case, errors will end up in the estimated drag coefficient.

3

Hi @spto,

What is the frame of your RADEC observations ?
Is it really the stationFrame (so a Topocentric frame, rotating with the Earth) or is it an inertial frame (like ICRF) centered at the station ?

A good way to debug the batch LS is to place a break point in class MeasurementHandler after the definition of the EstimatedMeasurement.
That way you can see what happens with the calculation of the theoretical measurement and check the residuals one by one.
I’m not sure you can do this with the Python wrapper though.

Maxime

4

Thank you @luc ! I have tried out your suggestions and not solving the problem yet:

======

I have checked my jd values, they are correct. The satellite passed at 22:21 (UTC+8), my jd values ends with 0.0985…, or about 2.364 jd hours, i.e. 14.364 UTC => 22.364 (UTC+8).

======
Seems that python does not has double type. Float is a double precision type for python.

======
After changing the MSAFE line to:

MarshallSolarActivityFutureEstimation("\p{Alpha}\p{Lower}\p{Lower}\p{Digit}\p{Digit}\p{Digit}\p{Digit}(?:f|F)10(?:_prd)?\.(?:txt|TXT)",MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)

there is an error:
DM.feed(msafe.getSupportedNames(), msafe)
TypeError: descriptor ‘feed’ for ‘DataProvidersManager’ objects doesn’t apply to a ‘str’ object

======

5

Hi @MaximeJ , thank you!

I have tried to change the stationFrame on line #161 and #292 to eme (FramesFactory.getEME2000()), (EME is said to be the other name for ICRF?), the results are:
Line #161 change makes no impact on the IodGooding result using 3 points. (exactly the same results)
Line #292 change will give the following error: Jacobian matrix for type KEPLERIAN is singular with current orbit

======
Do you mean I should check the result of AngularRaDec?

6

Do you know the frame in which your camera or post-processing returns the RADEC angles ?

Side note:
ICRF is the heliocentric reference frame, GCRF has the same axis but is geocentric.
And EME2000 is geocentric and has axes that are very close to GCRF (less than a µrad difference between axes)

I meant that you should check the value of the Estimated<AngularRADEC>, they are the theoretical RADEC that are computed through the propagation based on you model and the observations’ dates.
But again it’s better to be able to do this in debug mode.

7

I used Astrometry.net for plate solving and it returns J2000 RA/Dec. For astronomy purposes it doesn’t matter where the coordinate system centers, but for satellite observing it is centered on my camera (of course :smile:)

8

I’m not familiar with Astrometry.net, is it able to give you the RA/Dec of your object in J2000 by analyzing your pictures ?

Since Orekit’s EME2000 is centered on the Earth it won’t work as is. I think you may need to translate the EME2000 at the ground point of your station using a Transform or better, a TransformProvider.
And use this new frame as the reference frame when building your AngularRaDec objects.

9

Yes, Astrometry.net compares the stars in the picture with known catalogue and outputs the J2000 coordinates of the picture. (so I can calculate what coordinate is each pixel on)

Is there a standard example on how to set up the observer’s ground frame correctly? I guess it’s a very common situation that many will face. (I’m too new to Orekit to develop code without any existing examples… :joy:)

10

Hi @spto and @luc @MaximeJ ,

About the MSAFE, you can use it like this in python:

        msafe = MarshallSolarActivityFutureEstimation(
            MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
            MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
        # Feeding the F10.7 bulletins to Orekit's data manager
        DM = DataContext.getDefault().getDataProvidersManager()
        DM.feed(msafe.getSupportedNames(), msafe)

While I also tried the regular expression string directly, it is OK too.

11

Hi there,

I am actually curious as to how outputs from astrometry.net should be processed with Orekit.
@MaximeJ I am a bit confused: the ground station is passed to AngularRaDec and in theoreticalObservation the relative position from station to spacecraft is computed, so how can the origin of the frame matter?

Cheers,
Romain.

12

Thank you @lirw1984 ! This solved the Atmospheric Drag bug.

=======

The following discussion is based on my current version of the program. For the convenience of everyone, I copied some key lines here (about the definition of frames):

sitella = GeodeticPoint(site_lat, site_lon, site_alt)
stationFrame = TopocentricFrame(wgs84Ellipsoid, sitella, “Ground Station”)
station = GroundStation(stationFrame)
eme = FramesFactory.getEME2000()
radec = AngularRaDec(station, eme, JDobj, observedRaDec, sigma, baseweight, sat)
estimate_Gooding = IOD_Gooding.estimate(eme, obs1, obs2, obs3, range_1, range_3)

For the overall behaviour of the program, my progresses are:

1, I realized that I’m suffering from the Mars Coordinate problem in China (which casues a few hundred meters dispalcement of the ground station). But correcting the coordinates only brings a small difference on the results.

2, I tried different combinations of inputs for IOD_Gooding. The reaction of IOD_Gooding is interesting: 1) it gave Semimajor axis (a) < Earth radius and Eccentricity (e) ~0.5 (which should be ~0) when I only provided the initial or last part of the observation arc;
2) it gave the best estimation of a and e when the input observation arc is centered on the highest elevation point during this pass (data row number ~1600);
3) the total temporal length of the observation arc shouldn’t be too long or too short – e reaches the smallest (~0.004) when the total temporal length is ~6 min.

3, BLS only works when the additional measurements added is very so little few: typically 2-3 points. And it’s very unstable, with the selected added data shifting for several rows in time (0.1s difference per row), BLS may change from working to failing. When I tried to add 4 or 5 points, it never succeeded. But – I have 3000+ samplings on the satellite passing arc, it is a great waste if I cannot use all or most of them. And apparently it would not reach optimal precision if most data are thrown.

4, I have attempted to play with the thresholds, sigma values, position error, iteration / evaluation numbers, and other parameters. They seem not affecting the outcome very much.

5, I used .getObservedValue() to extract the (ra,dec) values stored in the output of AngularRaDec function. The values are exactly the same with what I have input. Is that normal or abnormal? Should I expect AngularRaDec to change the (ra,dec) values by frame transformations?

I have attached my current version of code here (though not much has changed):

Any ideas on what should be done next to improve the pipeline to include more date into the BLS?

13

Hi @Serrof,

You are right, sorry for the confusing comment.
I got mixed up with another code I’m using.
I should have remembered since I wrote the AngularRaDec code…
In my defense it was a while ago :wink:

14

@MaximeJ no worries, thanks for the clarification

@spto For 5/ I think you want to use getEstimatedValue() rather than getObservedValue(), as the latter is indeed the one you passed as real measurements.

Cheers,
Romain.

15

To exclude the possibility that my data is bad (bad measurement or Mars Coordinate, etc.), I generated a simulated dataset from known TLE, without adding any noise.

But the behaviour is similar:
IOD_Gooding gave result that is very close to the truth. Below is a comparison of the IOD_Gooding result(left) and the known TLE(right):

a: 7481006.339399883_____7482080.055577391
e: 8.523466150913378E-4__0.0003123
i: 123.00476856711326____ 123.0023
pa: 134.92971659662658___83.0462°
raan: 11.670946768399862__28.3223°
v: -17.20631076777036_____-83.08182538373174

a,e,i are almost exactly correct, while pa, raan, v has considerable errors. While pa and v errors are unavoidable for an almost circle orbit, the raan error size is a bit bigger than I would expect.

======
But for BLS, I could still add only 2 more points. 3 additional points will make the pipeline collapse.

I think that means:
1, the way I define the frames should be correct, or the a, e, i result wouldn’t be so close to truth.
2, something is very wrong about my usage of BLS. Since I’m using an artificial perfect dataset now, it would only be the way I have adopted these functions that went wrong.

Any ideas?

Attached is the simulated dataset:

16

Hi @spto

one thing to keep in mind is that the RA-Dec theoretical measurements within the Orekit OD process are computed with light time correction. On the other hand, as far as I am aware, it is not the case for the implementation of Gooding’s algorithlm, that thinks that the angular information is instantaneous.

But now that you are using synthetic data, I would suggest skipping all together the IOD and focus on the OD itself (initializing it with the actual orbit) until you manage to have it converged. Since you are using a lot of data point, maybe use the OutlierFilter? I’ve tried to use it myself recently but without much success, as it seems to either accept all or reject all measurements. Perhaps the Orekit dev team has some tips as on how to set it up properly?

Best,
Romain.

17

Hi @Serrof, @spto,

Filtering outliers strongly depends on the configuration of the OD.
If you suspect your first iterations are far away from the truth you can set the warm up to larger iterations (like 2 or 3) so that the filter activate later in the OD process (that way it won’t remove all the measurements).
If all measurements are accepted, you can try diminishing the multiplier maxSigma so that the filter gets more restrictive when it’s activated.

@spto, some thoughts to help you:

  • In your example you’ve set the RADEC error (or sigma) to 0.2 radians which is about 11.46 degrees.
    This is too big in my opinion, you should try to put something closer to the real resolution of your optical observations (8 arcsec gives 3.88e-5 rad).
  • If your initial orbit is a bit far from the expected solution you can try using a LevenbergMarquardtOptimizer instead of a GaussNewtonOptimizer in your BatchLSEstimator. It has a bigger convergence radius.

Also, one thing you can try is to plot the residuals at iteration 0, that means the residuals of the measurements with respect to your initial orbit and propagation (without any least-square correction).

You could extract them from a BatchLSObserver at iteration 0 of the BatchLSEstimator (see an example of an OrbitDeterminationObserver in the Java tutorials.

Or you could use a specific OrekitStepHandler.
I’ve written these two, they are in Java so you will have to translate them (the first one should be enough):

  1. ResidualsHandler.java (3.8 KB)
    Single satellite version. Once it is used, call the method getEvaluations to have a Map of <ObservedMeasurement, EstimatedMeasurement>.
    ObservedMeasurement are your original measurements, EstimatedMeasurement are the theoretical evaluation of the measurements made by Orekit given your intial orbit and propagation model.
  2. MultiSatResidualsHandler.java (4.1 KB)
    Multi satellite version. Works the same but with a PropagatorParallelizer instead of a classical Propagator.

Both are mere copies of Orekit’s MeasurementHandler.
The difference is that they don’t refer to any batch least-square model.

Usage of the single satellite one (still in Java).

// Suppose you have a sorted list of measurements and a fully configured propagator builder
List<ObservedMeasurement<?>> measurements;
PropagatorBuilder propagatorBuilder;

// Get start/end date
AbsoluteDate startDate = measurements.get(0).getDate();
AbsoluteDate finalDate = measurements.get(measurements.size() - 1).getDate();

// Build the propagator
Propagator propagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());

// Initialize single sat step handler
ResidualsHandler residualhandler = new ResidualsHandler(measurements);

// Add the handler to the propagator
propagator.getMultiplexer().add(residualhandler);

// Propagate from first obs to last obs
propagator.propagate(startDate, finalDate);

// Get the map of results
Map<ObservedMeasurement<?>, EstimatedMeasurement<?> = residualhandler.getEvaluations();

Then, comparing the values of each EstimatedMeasurement.getEstimatedValue() with its EstimatedMeasurement.getObservedValue() will give you the residuals, that you can plot to better understand what happens.

Hope this helps,
Maxime

18

Thank you @MaximeJ for the detailed suggestion!! Let me learn and try :pray: