SOCIS 2019: Orbit Determination with SatNOGS

noeljanes · July 17, 2019, 6:15pm

Thanks for those corrections.
I ran the orbit determination with the following settings on date and p-v in the maxvalier.in file:
orbit.date = 2019-06-11T16:35:29.149
orbit.cartesian.px = -5273194
orbit.cartesian.py = -4450295
orbit.cartesian.pz = 110075
orbit.cartesian.vx = -585
orbit.cartesian.vy = 790
orbit.cartesian.vz = 7610
and drag.cd.estimated = false

The px,py,pz coordinates correspond to an orbital position of 503.4km and vx,vy,vz correspond to an orbital velocity of 7673m/s (compared to the 475km and 7643m/s used in the generator program.)

When run through the orbit determination program it yielded enormous residuals (up to 7million) and an orbital altitude of that appeared much too large. Upon some investigation I found that the generated measurements were in m/s rather than the km/s used by the orbit determination program, so I changed this to the correct units.

Upon running it again, with no other changes, I got the same error as before (although with smaller RMS values, so that’s something):

I have tried to alter the orbital parameters to: px = -5423194 , py = -4550295 , pz = 250075 , vx = -735 , vy = 740 , vz = 7760, but got the same error:

I have attached the generated files for the initial run with doppler measurements in m/s.
orbit-determination-range-rate-residuals-generated.out (79.7 KB)
orbit-determination-log-generated.out (5.1 KB)

luc · July 18, 2019, 7:52am

The error about altitude below 120km is related to the atmosphere model used for drag force computation. The orbit determination tutorial has been set up using the DTM2000 model which does not extends down to ground. Such types of limitations are frequent in atmosphere models, as they are intended for satellites. At 120km, you are already very low, at reentry level. We could use another atmosphere model (and in fact a much better one as DTM2000 is really not recommended anymore). The NRLMSISE2000 model is now the one we recommend, and it has the additional benefit of extending down to ground.

This, however, is only a sign of something else wreaking havoc: we should not go from 500km to 120km with a “normal” drag coefficient and no maneuver. In fact, the OD process get lost right at the beginning. Looking at your first iteration, it tries to move the satellite by almost 14000km. This can be avoided by reducing the initial step bound factor (parameter estimator.Levenberg.Marquardt.initial.step.bound.factor). This parameter has some black magic included: there are no perfect rules to set it, its mostly trial and error. When close to Earth, as atmosphere density changes exponentially, some test orbits the optimizer wants to evaluate may end up triggering the exception you have seen, so we try to prevent it testing orbits too far, even if it implies doing more iterations. 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.

So I tried with this:

orbit.date = 2019-06-11T16:35:29.149
orbit.cartesian.px = -5273194
orbit.cartesian.py = -4450295
orbit.cartesian.pz = 110075
orbit.cartesian.vx = -585
orbit.cartesian.vy =  790
orbit.cartesian.vz = 7610
...
estimator.Levenberg.Marquardt.initial.step.bound.factor = 100.0

Then the process converges with the following results:

iteration evaluations      ΔP(m)        ΔV(m/s)           RMS          nb Range    nb Range-rate  nb Angular     nb PV
     0          1                                 502.559212955049       0/0          713/713        0/0          0/0   
     1          2       41897.042311  5.372745063 498.755565145088       0/0          713/713        0/0          0/0   
     2          3       76995.153909 11.272244311 491.320402999840       0/0          713/713        0/0          0/0   
     3          4      104784.168766 24.667473984 477.497649253095       0/0          713/713        0/0          0/0   
     4          5       53358.435567 51.472432766 442.321488682624       0/0          713/713        0/0          0/0   
     5          6       12346.634942 78.272402581 309.656257555581       0/0          713/713        0/0          0/0   
     6          7       10743.579822 92.299454042  85.332350236302       0/0          713/713        0/0          0/0   
     7          8      171563.213336 170.204556144  63.152158757981       0/0          713/713        0/0          0/0   
     8          9       60045.019913 178.367553223  34.203601645514       0/0          713/713        0/0          0/0   
     9         10       75703.814940 61.910993424   9.203202160007       0/0          713/713        0/0          0/0   
    10         11       12612.507287 25.106258191   1.093589206022       0/0          713/713        0/0          0/0   
    11         12         456.426950  1.202983056   1.049911494077       0/0          713/713        0/0          0/0   
    12         13           1.698701  0.001659690   1.049911167960       0/0          713/713        0/0          0/0   
    13         14           0.002574  0.000002804   1.049911167960       0/0          713/713        0/0          0/0   
Estimated orbit: Cartesian parameters: {2019-06-11T16:35:29.149, P(-5253009.438816775, -4400520.679405428, 80427.86229386796), V(-558.9409310057042, 764.267511436562, 7584.967970138002), A(6.505523654578354, 5.449769642672016, -0.09960487730587128)}
Estimated orbital parameters changes: 
   1 Px                             +20184.561183225363  (final value:  -5253009.438816775000)
   2 Py                             +49774.320594572460  (final value:  -4400520.679405428000)
   3 Pz                             -29647.137706132038  (final value:   80427.862293867960)
   4 Vx                             +26.059068994296  (final value:  -558.940931005704)
   5 Vy                             -25.732488563438  (final value:   764.267511436562)
   6 Vz                             -25.032029861998  (final value:   7584.967970138002)
Estimated propagator parameters changes: 
Estimated measurements parameters changes: 
   1 39-CGBSAT-VHF/range rate bias  -0.009368765028  (final value:   374.990631234972)
Number of iterations: 13
Number of evaluations: 14
Measurements type: range-rate
   number of measurements: 713
   residuals min  value  : -48.68390007557355
   residuals max  value  : 45.32828574650921
   residuals mean value  : -2.8188214470925566E-9
   residuals σ           : 15.75972309642003
wall clock run time (s): 80.059

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

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?

We are improving, though.

noeljanes · July 18, 2019, 1:56pm

Hi Luc, thanks for the corrections.

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:

What do you mean by 457m in position?

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 · July 18, 2019, 2:29pm

Good!

Yes, this is this model.

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.

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).

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 · July 22, 2019, 9:41am

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 · July 22, 2019, 11:56am

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 · July 22, 2019, 12:11pm

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 · July 23, 2019, 7:53am

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

noeljanes · July 24, 2019, 7:35pm

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 · July 25, 2019, 1:54pm

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 · July 29, 2019, 11:19am

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

noeljanes · July 31, 2019, 6:48pm

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 · July 31, 2019, 7:56pm

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 · July 31, 2019, 9:33pm

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

noeljanes · August 1, 2019, 4:47pm

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?

How would I best go about doing this?

luc · August 1, 2019, 6:40pm

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 · August 1, 2019, 6:51pm

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 · August 2, 2019, 5:39pm

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

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.

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 · August 6, 2019, 3:23pm

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 · August 6, 2019, 4:59pm

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.