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