Orbit Determination using PVT input

Hello,

I was wondering if there’s a way to modify the W3B.aer file to handle the input as position, velocity, and time instead. The measurement data I have are the GPS data (ECEF) from the satellite. Initially, I tried converting the PVT into AER, but I’m curious if there’s a better way Orekit can handle this issue.

Thank you,

Non Meeboon

Hi @Non, welcome

Yes, it is possible and quite straightforward, with one minor challenge.
In the Orekit tutorials, the various orbit determination tutorials extend AbstractOrbitDetermination which contains the readMeasurements method that parses the measurements file. When a line is read from this measurements file, the second field in the line specifies the measurement type. Depending on this type, a specific line parser is used to parse the other fields and create the measurement. Supported types are RANGE, RANGE_RATE, AZ_EL and PV. You can look at the switch statement in this method and at the various line parsing classes (RangeParser, RangeRateParser, AzElParser and PVParser). If you want to add a POS type, you just have to create a new alternative in the readMeasurements switch statement and a new PosParser class. The parser will create a Position measurement, which already exists in Orekit (since version 9.3).

Now comes the challenge. The Position measurement expects the position to be expressed in the inertial frame used for orbit determination. If your coordinates are in a terrestrial frame, you will first have to convert it to inertial frame during parsing. It may be wise to have the frame used for measurements parsing customizable, perhaps as a field in the measurement line to have lines of the form:

2010-11-02T03:21:48.2653   POSITION      ITRF-2010      12345.678   34567.89    -653.23433

Thank you very much, Luc. I really appreciate your help. I will do what you suggested and let you know the result very soon.

Hello Luc,

Thank you for your suggestion. I was able to define the input as PV in the measurement file.

Currently, I am trying to compare the result of OD from Orekit and from the software I have. The difference of OD result between the two system is a bit off approximately up to 100 metres on position xyz axis. I was wondering which input parameters to adjust to get the results which will be closed to each other. I tried turning on/off the drag and solar radiation pressure but the change is only slightly a few centimetres. Since my measurement files are only PV type, I did not adjust any parameters related to azimuth, elevation, and range. I tried changing PV.measurements.position.sigma and PV.measurements.velocity.sigma, but I still could not get the right number. Can you help me point out which input parameters I should be looking into?

This is the input of Orekit OD:
px = -3929923.912622450 m
py = 5175319.620766900 m
pz = 3105571.431586890 m
vx = -933.112766490 m/s
vy = 3275.756285987 m/s
vz = -6616.034855312 m/s

This is the results from Orekit OD:
px = -3929941.971366900 m
py = 5175333.548900160 m
pz = 3105591.133806370 m
vx = -933.229532925 m/s
vy = 3275.784087278 m/s
vz = -6616.046099862 m/s

Here is the OD results from the software I have:
px = -3929912.022938370 m
py = 5175261.944422130 m
pz = 3105654.409526630 m
vx = -933.173860502 m/s
vy = 3275.832151294 m/s
vz = -6616.002801524 m/s

Input Parameters I’ve used:
Drag cd = 2.3
Drag area = 5.0 m^2
Solar radiation pressure cr = 1.0
Solar radiation pressure area = 1.0 m^2
Central body degree = 9
Central body order = 9
PV outlier rejection multiplier = 2
PV outlier rejection starting iteration = 2
PV measurements base weight = 2
PV measurements position sigma = 10
PV measurements velocity sigma = 0.01

Sorry for the long post. Thank you very much.

Yes, 100m is indeed very high. I am surprised drag only has centimeter level effects, how long is the measurement time span ? Does it covers at least one full orbit?
Another point may be to check the frames, mainly at measurement conversion stage. You should check that your orekit-data does include proper EOP for the measurements dates.

The measurement time span covers from 2020-08-22T04:54:45.890 to 2020-08-24T02:53:09.898. The last point is what I input as an orbital parameter. One orbital period is around 101 minutes. The time step for each point is 128 seconds. I’ve attached the measurement data. W3B_skptest.txt (68.3 KB)

Yes, with/without drag changes only a few centimeter. I was using NRLMSISE00 drag model. May be modified it wrongly? Here is the code:

For the frame transformation, I followed an example you discussed in the other Orekit forum. Here is the code:

I used the EOP data from Orekit 10.0 sources. Is there a new update in EOP? Where can I find the most updated ones?

Yes, they must be updated. In fact, the date that is included in the source tree is only intended for tests purposes, not for real computation. some of the data are updated regularly (yearly, monthly, weekly depending on data). As an example, operational system use rapid data for EOP, which are often retrieved weekly from BulletinA files. One way to get slightly less outdated data is explained in the data configuration page here: https://www.orekit.org/static/data/default-configuration.html#Quick_recommended_setup. I have just updated the data there, but the update is not very frequent. The orekit-data project however includes an update.sh script that can be used to update the data directly from the source sites. It is users responsibility to update this data set, the one linked above is only a convenience and no guarantee that it will be updated (or even maintained).

Hi Luc,

I’ve updated the EOR data using the update.sh script. The result; however, did not change much from before (only centimeters level). Can you suggest which parameters should I be looking into next? What would be the appropriate values for PV weight and sigma? How many orbit should the measurement data cover?

Thank you

Here are a few hints, but I’m not sure it will help, some are even minor details.

When you build the OneAxisEllipsoid, don’t get the equatorial radius from the gravity field, but rather use the constant from the same place you get the flattening (i.e. Constants.WGS84_EARTH_EQUATORIAL_RADIUS). The radius in a gravity field, is basically a scaling factor related to dynamics. It could have been chose by gravity field designers to any value, even one not really related to geometry. The radius used to build OneAxisEllipsoid on the other had is a geometric model that precisely defines a reference shape.

Instead of an empty InputParameters for atmosphere model, you could switch to Orekit version 10.2 and use the new CssiSpaceWeatherData that was contributed by Clément in this version, it will use up to date flux (and of course you should also update the flux data somehow by yourself).

In the frame transform, you could void calling getTransformTo twice, just save the transform (or better save the full transormed PV) rather than recomputing everything and extracting once position and then velocity. This is just an optimization, the result will not change at all and your implementation of frame transform seems fine to me.

I am wondering were the velocity in vxList, vyList and vzList comes from? Are they already present in the file you loaded? Do you compute them by finite differences from positions in the file? If they are present in the file, there may be a glitch as velocity in ECEF may be computed either as velocity wrt inertial frame projected to ECEF by just using a rotation or velocity wrt inertial frame transformed to ECEF by applying both rotation and velocity composition so they are really velocity wrt Earth frame. Here you assume the velocity are expressed wrt Earth frame and you convert them back. If they were in fact inertial velocities just rotated, then to get the inertial velocities back you should do it differently (i.e. call transform.transformPosition on the position part and transform.transformVector on the velocity part to prevent Orekit to include velocity composition). I don’t think there is a problem here however, you most probably did the transform the right way because if you didn’t, you would have really inconsistent results. If you computed the velocities by finite differences from position-only ECEF data, then this would probably introduce some inaccuracies. The velocity data would really represent velocity wrt ECEF (which is good and corresponds to what you transform back to inertial frame), but they would include interpolation errors that would skew the orbit slightly. So if in your file you have only position (or if you don’t trust the velocity that is present), then you should probably better ignore the velocity and feed Orekit only with Position masurements, not PV measurements (then back to the initial posts in this thread: you’ll have to add support for Position in the tutorial).

If your data include only PV measurements and you don’t trust some of them more than the other ones, the weight is irrelevant, you could scale all of them up or down by the same factor without any effect. Weight is only important when you mix data with different trust levels, this factors allows you to hint the optimizer to use more some measurements than other ones. It is rarely used. In most cases everybody will just use the default weight set to 1.0 and only use the measurements sigma as a way to normalize measurements.

The sigma should (roughly) represent the standard deviation of the measurement errors. If your data comes from a GNSS receiver, you probably have a data sheet from the manufacturer. You should typically expect meter-level sigma if it uses only code measurements, and centimeter-level if it uses phase measurements as well. If you don’t have these information, you can retrieve the residuals from the output of the tutorial printed at the end of a run (look at “residuals σ” for each measurement file in the log).

In the list above, I think the CssiSpaceWeatherData and the Position only measurements are the most important ones.

Hello Luc,

Thank you very much for your suggestion. I think my previous measurement data was not enough to perform OD. I’ve tried adding more measurement data (date back 4 days before the reference epoch) as my input and the result was really impressive! The difference came down to only a few meter on each xyz axis. I will look into the new weather data and try to modify the drag coefficient and area in order to optimize the result.

Thank you very much. Have a nice day!

Non

Hi luc,
I am learning hwo to operate Orbit Determination, I am modifying the configure file: .\orekit-10.3\src\test\resources\orbit-determination\W3B\od_test_W3.in . But I donot know how to config third body force beside sun and moon. By the way, is there any more detail description for this configure file?
The another question is: there are 60-100m difference between Orekit and my own Orbit Determination program when processing the same observe data, I wonder if it is because of different frame. So, I want test by modifying “body.frame” keywords, is there other value except “CIO/2010-based ITRF simple EOP” ? Thank you so much!

Hi @Tanggeshi

That’s because you are using the Orekit tests. For the Orekit tests, only the Sun and Moon attractions can be added to the dynamical model. I would recommend you to use instead the orekit-tutorials. Indeed, for the Orekit tutorials, you can add the third body attraction for the main bodies of the Solar System (i.e., Saturn, Jupiter, etc.). You can look at for instance the input file of LaserRangingOrbitDetermination tutorial for which we add Venus attraction in the calculations. The file is located in src/main/resources/slr-od folder.

Please note that Orekit orbit determination accuracy has been strongly validated against centimeter level reference orbits. You can check for instance if the models used to define the satellite dynamic and the measurements are the same. Furthermore, I would recommend you to make sure that you are comparing orbits in the same orbital frames.

You can access the list of predefined frames in the Predefined.java enumerate.

I hope my answers will help you.

Best regards,
Bryan

Hi Bryan,
Your answers are very helpful, I can learn more detail by your help!
Thank you so much!

Best regards,
TangGeshi

Hi Bryan,
After the new Space Weather file downloaded, I find the format is little different, and orekit runs error as follows:
Exception in thread “main” org.orekit.errors.OrekitException: unable to parse line 2,029 of file D:\ODSOFT\PODsoft\callOrekit\orekit-data\CSSI-Space-Weather-Data\SpaceWeather-All-v1.2.txt:
2021 06 28 2562 26 13 13 13 13 13 13 13 13 104 5 5 5 5 5 5 5 5 5 83.0 78.1 77.2 80.3 75.7 75.5
at org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherDataLoader.loadData(CssiSpaceWeatherDataLoader.java:502)
at org.orekit.data.DataProvidersManager$MonitoringWrapper.loadData(DataProvidersManager.java:404)

Would you please give me some suggestion to solve this problem? Thank you!

Here is the screen:

Hi @Tanggeshi

Could you try to update again these data?

I just imported them (i.e. 28th of June at 7.00 am UTC), ran a test using them and everything work correctly.

Bryan

I had updated the latest file, it is the same error.
SpaceWeather-All-v1.2.txt (281.9 KB)
I had upload the file, would you please check it for me? I put this file in the directory: D:\ODSOFT\PODsoft\callOrekit\orekit-data\CSSI-Space-Weather-Data, and make the root directory as follow: File orekitData = new File(“D:\ODSOFT\PODsoft\callOrekit\orekit-data”);
DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
manager.addProvider(new DirectoryCrawler(orekitData));

Thank you!

Tanggeshi

I see.

Below, is a screenshot of your file.

Now, a screenshot of my file.

As you can see, they are missing data in your file. I think that’s why the parser is not able to read your file.

How do you update the space weather data file? Do you use the update.sh script available in the orekit-data folder? Using this script make very easy to update all your orekit-data in few seconds. That’s what I used yesterday to update my data.
If not, I can recommend you to use this script available in the orekit-data project.

Otherwise, you can use my own file:
SpaceWeather-All-v1.2.txt (3.0 MB)

Bryan

Hi,

The file given on Celestrak website (Tanggeshi’s one): https://celestrak.com/SpaceData/SW-Last5Years.txt is different from the file on AGI ftp (Bryan’s one): ftp://ftp.agi.com/pub/DynamicEarthData/SpaceWeather-All-v1.2.txt
That explains why the Orekit parser fails on the Celestrak’s file; but I don’t know why there are such differences between the two files.

Maxime

Hi,

When I wrote the parser for CSSI space weather files, I automatically assumed that the files from both the AGI FTP and Celestrak were identical, which turns out not to be the case as you pointed out.

The MONTHLY_PREDICTED and MONTHLY_FIT values in Celestrak’s file don’t contain any Kp nor Ap values, so it makes it unusable for both NRLMSISE00InputParameters and DTM2000InputParameters.

We should stress in the documentation that only the AGI FTP file should be used. We could also consider improving the error handling in the file parser.

Cheers
Clément