SOCIS 2019: Orbit Determination with SatNOGS

[noeljanes]

Hi Luc, thanks for the corrections.

@luc

The second attempt with vx set to -735 is clearly already too far, 160m/s difference with the real orbit is a tremendous initial error.

I had a feeling this would be the case, but felt that it was worth trying it with a significantly larger orbit to see if that would get it to converge (not that it worked out.)

I initially ran the orbit determination with the same parameters that you set and got the same results. I then changed the atmospheric model from DTM2000 to the NRLMSISE2000 model, (this is called NRLMSISE00 in Orekit correct?,) both in the generator and the orbit determination program, and reran the orbit determination which yielded very similar results (which isn’t particularly surprising considering the atmospheric density at around 500km.)

Here are the residuals that I got for both models:

@luc

This is much better than earlier, but the Doppler noise prevents getting a good orbit. We are at 457m in position.

What do you mean by 457m in position?

@luc

We estimated the range rate bias at almost 0, though, there is perhaps a problem there, it should be estimated at 375m/s as it is what is configured in the generation program. Did you remove the bieas before generating the data?

Both the bias in the maxvalier.in file and in the measurement generation program are set to 375m/s with a sigma of 15 in both.
Below are the variables in the maxvalier.in file:

ground.station.range.rate.sigma                [0] = 15
ground.station.range.rate.bias                 [0] = 375.0
ground.station.range.rate.bias.min             [0] = -0.0
ground.station.range.rate.bias.max             [0] = +500.0
ground.station.range.rate.bias.estimated       [0] = true

[luc]

@noeljanes

I initially ran the orbit determination with the same parameters that you set and got the same results.

Good!

@noeljanes

this is called NRLMSISE00 in Orekit correct?

Yes, this is this model.

@noeljanes

which yielded very similar results (which isn’t particularly surprising considering the atmospheric density at around 500km.)

Well, it’s not that obvious. Atmosphere model may be very different from each other: atmosphere is very difficult to model accurately. It is also very time-dependent due to solar wind that injects high energy particles during solar storms (which creates beautiful auroras but also changes atmosphere density). The density can change by a factor of 2 in less than 3 hours! This density is low (we are in near vacuum) but a factor 2 in density is a factor 2 in drag force.

@noeljanes

What do you mean by 457m in position?

This is the difference between the position estimated by the orbit determination (look at the final values fro PX, Py, Pz, Vx, Vy, Vz) and the values used for measurements generation which correspond to the “true” orbit in this case (because we generate measurements from there).

@noeljanes

Both the bias in the maxvalier.in file and in the measurement generation program are set to 375m/s with a sigma of 15 in both.

My bad, I read the -0.0093 m/s in the output as the final value and not as the change between initial and final value. You are right, we end up estimating a bias at 374.99067 which is a good estimate.

Waiting for new real measurements, could you look at the influence of:

• estimating drag coefficient (assuming we start OD with a slightly wrong estimate wrt the one used for generation)
• standard deviation of noise in measurements generation
• adding one or two stations?

[cgbsat]

@noeljanes

I’ve seen on the SatNOGS website that you appear to have a UHF antenna as well as a VHF one, is that correct? If so are you able to run your own software in the UHF frequency range?

Yes, I can operate it in the UHF band. I also have an S-band antenna and can use it there as well. Using S-band is a bit easier as it does not require taking down a SatNOGS station.

Do you have any preference for specific satellites?

[yzokras]

In S-band you could listen to TechnoSat (NORAD ID 42829) at 2263 MHz, I can ask the operations team to turn on the S-band transmitter at the times you want.

Besides, as mentioned in this thread, TechnoSat has precise orbit determination from laser ranging.

[luc]

This would be great to gather some validation data for Orekit, but we should not forget that the goal of the SOCIS project is to check if we can do OD with the SatNOGS antennas.

[cgbsat]

What is TechnoSat transmitting on S-band? Is it just a carrier or does the signal contain data?

[noeljanes]

@cgbsat

Do you have any preference for specific satellites?

I was thinking we could look at the TechnoSat operated by TU Berlin, but if there are other more promising candidates then that would be better

[cgbsat]

I’ll have a look, but I will wait for a few days for the European heatwave to pass as I don’t want to risk overheating the PC and SDR.

[yzokras]

TechnoSat’s S-band transmitter does not support carrier-only, when it transmits it is on a 1.5 MHz wide band.

[noeljanes]

@luc

Waiting for new real measurements, could you look at the influence of:

• estimating drag coefficient (assuming we start OD with a slightly wrong estimate wrt the one used for generation)
• standard deviation of noise in measurements generation
• adding one or two stations?

I decided to use 7 different stations to look at the influence of these factors. The idea behind this was to confirm that the trends carried over between stations and weren’t as a result of other factors (for example potential errors due to low elevation above the horizon and such.)

The stations used were:

• 39-CGBSAT-VHF (52.834°, 6.379°, altitude 10m)
• ZL1WJQ (-37.188°, 174.875°, altitude 65m)
• F4KLD-UNIVERSITE-PAUL-SABATIER-Toulouse-III (43.562°, 1.469°, altitude 180m)
• Apomahon (38.048°, 23.739°, altitude 119m)
• Marcs (48.644°, 2.231°, altitude 110m)
• GU-virt (55.861°, -4.251°, altitude 49m)
• Kir-virt (67.840°, 20.409°, altitude 401m)

All bar the stations with the -virt suffix are part of the SatNOGS network, whereas I created the other two to simulate measurements from a wider range of latitudes and longitudes.

Changing Standard Deviation of noise in the measurement generation

I changed the standard deviation of the noise in both the measurement generation and the maxvalier.in files.
I noticed a general trend of a smaller difference between the estimated Cartesian parameters and those used in the measurement generation as the standard deviation decreased (difference in position is in the top row, and velocity is in the bottom row):

In general the various stations followed a similar trend with smaller Δ’s for smaller standard deviations, but the estimated position and velocity would sometimes shoot off to a much larger value (most noticeably at really small values of σ.)
I’m unsure why the OD sometimes converges on an incorrect result and sometimes doesn’t (are there much higher orbits that give a similar radial velocity to those in LEO?)

A sigma of 5m/s seems to be the best compromise between accuracy in the estimated parameters and reproducibility between stations.

Multiple Stations

For the multiple stations I set a range.rate.bias=375m/s and range.rate.sigma=15 and I got the following results:

.

I noticed that all the major outliers stemmed from combinations including the GU-virt station, so I decided to discard this station (and all combinations with it) got the following plot:

This shows that increasing the number of stations used in the OD does, in general, improve the accuracy of the estimated position. That being said, it did appear that a careful selection of stations could yield better results with fewer stations.
An example of this is the combination of 39-CGBSAT-VHF and Marcs which gives an estimated position that is 342m away from that used in the generation (which isn’t amazing, but we are most likely limited by the Doppler noise.)

Changing the coefficient of drag

When I tried to run the OD with the estimated CD flag set to true, they all failed with the following error:

When I set the CD to be an incorrect value (with drag.cd.estimated = false) I got the following series of results (a Cd of 2 is used in the measurement generation):

This made me wonder if this nearly flat behaviour continued to higher drag coefficients (I suspected that it didn’t, but wanted to confirm that this was the case,) so I added in another couple measurements to give the following:

These plots seem to indicate that a slight error in the drag coefficient used in the input file is tolerable, and doesn’t have much of an effect on the estimated parameters, as long as it isn’t wildly inaccurate.

Next steps
I was thinking of looking at the influence of adjusting the range-rate bias next, potentially with a variety of standard deviations as well, to see what impact that would have on the accuracy of the estimated parameters.

PS. I have uploaded all the data files used for these plots to the git repository

[luc]

Thanks for the results, which are promising.
I would say that 3-4 stations and a Doppler noise at 5 m/s or below is a target we could propose to SatNOGS network.
I have concerns about the drag coefficient. The error with hyperbolic eccentricity is probably due to an orbit that is changed too much by the process attempting to test a huge drag coefficient and the satellite diving in atmosphere (perhaps performing a reentry) and the shape of the orbit changing drastically (velocity becoming quasi vertical, hence eccentricity increasing a lot with a perigee below Earth surface, and an hyperbolic orbit). Perhaps you could prevent the coefficient to reach to high value by limiting it between 0 and 8 for example. However, the OD may want to increase it, so at the end you will see it is stuck at its upper limit. There is something that does not work properly here, but I don’t know what. The Orekit OD is used in several project and works well and in this case we even can’t blame the atmosphere since the measurements are simulated using the same models as used in OD (at least I think the models are the same, maybe you should double check we did not change anything between the two).

[noeljanes]

@luc

The Orekit OD is used in several project and works well and in this case we even can’t blame the atmosphere since the measurements are simulated using the same models as used in OD (at least I think the models are the same, maybe you should double check we did not change anything between the two).

Is this the atmospheric models, such as NRLMSISE2000, that we discussed a couple of weeks ago?

[noeljanes]

As far as I can tell we are mostly using the same models in both the measurement generation and OD.
I did notice one difference though:

• The measurement generation programme uses the Vienna3 model as it’s tropospheric model whereas the OD uses the SaastamoinenModel.

Would this be able to have the effect that we are looking for?

@luc

Perhaps you could prevent the coefficient to reach to high value by limiting it between 0 and 8 for example.

How would I best go about doing this?

[luc]

@noeljanes

Is this the atmospheric models, such as NRLMSISE2000, that we discussed a couple of weeks ago?

Yes, this is the model (or DTM2000), which is associated with a satellite shape model, either a BoxAndSolarArray with several facet or an IsotropicDrag. We should use the same model for
both programs, and I would prefer NRLMSISE2000 to DTM2000.

[luc]

@noeljanes

• The measurement generation programme uses the Vienna3 model as it’s tropospheric model whereas the OD uses the SaastamoinenModel.

I don’t think the discrepancy would be large, but we should fix it to be sure. I propose to select Vienna 3 for both programs.

I thought we had a configuration parameter for this, but it seems we forget ti put one. So fast hack could be to change lines 637-639 in OrbiyDetermination.java from

         if (driver.getName().equals(DragSensitive.DRAG_COEFFICIENT)) {
              driver.setSelected(true);
          }

to

         if (driver.getName().equals(DragSensitive.DRAG_COEFFICIENT)) {
              driver.setSelected(true);
              driver.setMinValue(0.0);
              driver.setMaxValue(8.0);
          }

A better way would be to add keys ParameterKey.DRAG_CD_MIN and ParameterKey.DRAG_CD_MAX in the enumerate and to parse them, so you can set the values in the configuration file instead of changing the source code. You can find an example around lines 1283-1288 for the on board clock offset whose min and max values can be configured by user.

[noeljanes]

@luc

Yes, this is the model (or DTM2000), which is associated with a satellite shape model, either a BoxAndSolarArray with several facet or an IsotropicDrag . We should use the same model for
both programs, and I would prefer NRLMSISE2000 to DTM2000 .

I just doublechecked, and both have been set to use the NRLMSISE2000 model

@luc

I don’t think the discrepancy would be large, but we should fix it to be sure. I propose to select Vienna 3 for both programs.

Looking through the OD program it doesn’t seem to correct for the troposphere for range-rate measurements at all, so I am currently working on adding it in.

@luc

A better way would be to add keys ParameterKey.DRAG_CD_MIN and ParameterKey.DRAG_CD_MAX in the enumerate and to parse them, so you can set the values in the configuration file instead of changing the source code. You can find an example around lines 1283-1288 for the on board clock offset whose min and max values can be configured by user.

Ok, thank you very much. I had a feeling that it would be done like that, but I didn’t want to begin poking at it without asking first.

[noeljanes]

With the Vienna model and the ParameterKey.DRAG_CD_MIN and ParameterKey.DRAG_CD_MAX added I reran the OD, but got the same hyperbolic error as before. I ran it at several different DRAG_CD_MAX, but none of them managed to finish an estimation run.

With DRAG_CD_MAX set to 8 in the config file:

7:

6:

5:

4:

3:

I added the configurable values using the following snippet of code:

if (cdEstimated) {
                for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
                    if (driver.getName().equals(DragSensitive.DRAG_COEFFICIENT)) {
                        driver.setSelected(true);
                    }
                    if (parser.containsKey(ParameterKey.DRAG_CD_MIN)) {
                        driver.setMinValue(parser.getDouble(ParameterKey.DRAG_CD_MIN));
                    }

                    if (parser.containsKey(ParameterKey.DRAG_CD_MAX)) {
                        driver.setMinValue(parser.getDouble(ParameterKey.DRAG_CD_MAX));
                    }
                }
            }

I’m a bit unsure of where to go from here. It would appear to be something to do with the generation of measurements as the OD finishes without issue if the measurements provided by @cgbsat are used instead of the generated measurements.
I’m unsure why this is as the two programs use the same models

Should the range rate measurements have a refraction correction in the OD? I can’t see it being used at all once created

[luc]

@noeljanes

if (parser.containsKey(ParameterKey.DRAG_CD_MAX)) { driver.setMinValue(parser.getDouble(ParameterKey.DRAG_CD_MAX)); }

It should probably be driver.setMaxValue. with this probable copy/paste error, you first set a min value, then overwrite it with a new larger value and do not set any max.

[noeljanes]

Oh, that’s embarrassing . I must have mistyped it, and then not seen it when double checking.
I fixed it, and now don’t get the hyperbolic orbit error. Instead it seems to fail to converge on a solution and hits a max number of iterations (600)

This was still the case when the max number of iterations was raised to 1000, but it seems to be converging (albeit very slowly)

Would it be worth trying to make it converge by raising the maximal count to something much higher? Or would it be better if I altered the output to provide it’s previous estimate when it hits the max count?

[luc]

@noeljanes

Would it be worth trying to make it converge by raising the maximal count to something much higher?

No, an orbit determination process that works should return a value in a few iterations (I would say 5 to 20). We often put the max to values like 20, 25 or 50 just to avoid wasting time for nothing. When the convergence is too bad, there is a deeper problem somewhere. The infamous
estimator.Levenberg.Marquardt.initial.step.bound.factor is sometimes a clue, but I don’t think it is the culprit here because it changes the start of the process (hence the initial in the name), and here the problem is final convergence.

I don’t know what happens, the log does not contain sufficient information. I guess the optimizer is somehow stuck in a loop, or a very slowly decreasing spiral arount optimal point, despite Levenberg-Marquardt is an algorithm considered to be very robust.

What you could do is try to understand what are the various points the optimizer uses. Currently, the log was designed for classical cases (i.e. convergence in a few iterations), where we clearly see the orbit move and converge, with a simple value representing the Euclidean distance between previous and current orbit reference point (i.e. the position at the initial orbit date) and the RMS of the residuals. You can look how this synthetic information is extracted during the process. This is done using an anonymous class implementing BatchLSObserver and passing it to the estimator by calling estimator.setObserver(new BatchLSObserver() { ... }).

The evaluationPerformed method of the observer is called at the end of each iteration and can extract data. currently, it calls the underlying loggers for measurements residuals and uses orbits[0].getPVCoordinates() to get the position and compute the distance with respect to the value stored from previous iterations. So the synthetic log ignores all the information about all estimated parameters changes. We should probably display them and look at their current values. When the method is called, it gets the estimated orbital parameters, the estimated propagator parameters (the drag coefficient belongs to this list) and the estimated measurements parameters (the range rate bias for example belongs to this list). Perhaps you could print the values for these parameters too? Each ParameterDriver hold both its name and its current values (and also min, max, scale factor and other things but we don’t care about them). So you could print big lines with the iteration number and a list of all parameters in name=value pattern. Also note that the current anonymous BatchLSObserver prints the data both on standard output (which corresponds to the screenshot you use) and also print the in a file if the output.base.name key has been set up in the configuration file. Looking in the file may be easier than looking on the console output if the lines are too long.

With such a modification, the log could help us understand what the optimization tries to do.

Another thing you could try easily would be to change the optimizer by replacing

estimator.optimization.engine = Levenberg-Marquardt

by

estimator.optimization.engine = Gauss-Newton

We currently have only these two alternatives: Levenberg-Marquardt or Gauss-Newton.
The estimator.Levenberg.Marquardt.initial.step.bound.factor is of course ignored if we choose Gauss-Newton.