Help about estimated propagator by orbit determination

Hi, all,

I am trying to do laser range orbit determination, following orekit-tutorials.
The residual is pretty good (numerical, ~few cm), but the orbit propagated by estimated propagator is ‘wrong’, validated by sp3, especially for leo, like stella.

I dig and dig. Now I am sure that the transit states used in the last iteration of orbit-determination are different with the states propagated by the estimated propagator, mainly caused by gravity force. So the residual is good but the propagated orbit is ‘wrong’.

But why?

I wrote a test, only gravity force, ignore other forces, and station range bias is not estimated. In the case of degree=30, order=30, the differences can be reached tens of meters. Lower order, smaller difference.

Please help me to check it.

orekit in 12.0 and orekit-tutorials in 12.0.

There are some changes in of orekit-tutorials.

(1) After estimation is done, use ephemeris mode of estimted propagator (in the method run(File)).

            final Orbit estimated = estimatedPropagators[0].getInitialState().getOrbit();

            // use Ephemeris mode of estimated propagator.
            final Propagator estimatedProp = estimatedPropagators[0];
            final EphemerisGenerator generator = estimatedProp.getEphemerisGenerator();
            estimatedProp.propagate(estimated.getDate(), estimated.getDate().shiftedBy(4 * Constants.JULIAN_DAY));
            final BoundedPropagator ephProp = generator.getGeneratedEphemeris();

(2) do the comparison of transit state (in the method run(File)).

            // do compare
            // use measurements list instead of Map.entrySet to keep in date-order.
            //    entry.getKey()==>observedMeasurement, entry.getValue()==>estimatedMeasurement
            System.out.println("              transit_date d_pos_m");
            final Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> lastEstimationsMap = estimator
            // Compute some statistics
//            for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
            for (ObservedMeasurement<?> observedMeasurement : measurements) {
                final EstimatedMeasurement estimatedMeasurement = lastEstimationsMap.get(observedMeasurement);
                // Send to dedicated log depending on type
                final String measurementType = observedMeasurement.getMeasurementType();
                if (measurementType.equals(Range.MEASUREMENT_TYPE)) {
                    final EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) estimatedMeasurement;

                    // compare the transit state used in the last iteration of orbit-determination w.r.t propagated by the estimated propagator.
                    // the difference in pos should be small enough (like ~1cm), right? But, many meters. 
                    // degree=30, order=30
                    // 2023-01-01T11:24:21.407144   0.002
                    // 2023-01-01T22:37:52.106740   4.527
                    // 2023-01-02T11:00:12.306801  13.736
                    // 2023-01-02T22:13:11.107262  21.306
                    // 2023-01-03T12:18:12.106309  26.779
                    final SpacecraftState transitState = evaluation.getStates()[0];
                    final AbsoluteDate transitDate = transitState.getDate();
                    final SpacecraftState estimated_transitState = ephProp.propagate(transitDate);
                    final Vector3D pos = transitState.getPosition();
                    final Vector3D estimated_pos = estimated_transitState.getPosition();
                    final double deltaPos = estimated_pos.distance(pos);
                    System.out.println(String.format("%s %7.3f",
                            transitDate.toStringWithoutUtcOffset(TimeScalesFactory.getUTC(), 6), deltaPos));
                } else if (measurementType.equals(RangeRate.MEASUREMENT_TYPE)) {
                    final EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) estimatedMeasurement;
                } else if (measurementType.equals(AngularAzEl.MEASUREMENT_TYPE)) {
                    final EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) estimatedMeasurement;
                } else if (measurementType.equals(PV.MEASUREMENT_TYPE)) {
                    final EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) estimatedMeasurement;

(3) Ignore ‘unkonwn’ station’s observation in the crd file (in the method readCrd()).

            // Station data
            final StationData stationData = stations.get(String.valueOf(header.getSystemIdentifier()));
            if (stationData == null) {
                // ignore 'unkonwn' station's observation
            } (133.8 KB)

odin_stella_20230101T00.yaml (5.6 KB)
stella_202301.np2 (69.7 KB)

Best regards,

Additional, if order=0. The unit is meter.

                    // degree=30, order=0
                    // 2023-01-01T11:24:21.407144   0.004
                    // 2023-01-01T22:37:52.106740   0.005
                    // 2023-01-02T11:00:12.306801   0.112
                    // 2023-01-02T22:13:11.107262   0.322
                    // 2023-01-03T12:18:12.106309   0.300

If degree=2, order=0.

                    // degree=2, order=0
                    // 2023-01-01T11:24:21.407144   0.004
                    // 2023-01-01T22:37:52.106740   0.010
                    // 2023-01-02T11:00:12.306801   0.010
                    // 2023-01-02T22:13:11.107262   0.016
                    // 2023-01-03T12:18:12.106309   0.021

Hi @lirw1984,

I’ve tried to reproduce your results but I’m missing the SINEX files for the station positions.
Could you please share them?

I don’t have an answer yet to your issue but it reminds me of this one: Re-computed estimated measurements.

I don’t think it will change a lot but, in your input file, could you try replacing:

# Propagator definition
  # Numerical integrator (min step (s), max step (s) and position error (m))
    minStep: 0.001
    maxStep: 300
    positionError: 10.0


# Propagator definition
  # Numerical integrator (min step (s), max step (s) and position error (m))
    minStep: 0.001
    maxStep: 60.
    positionError: 0.001

And tell us how it compares with this configuration ?

Hi @MaximeJ ,

I can’t thank you enough.

I tried this configuration, the difference is 0.000 m. That is, it is perfectly solved.

Once I tried reducing positionError as 0.001 for lageos2, the improvement is not significant. I’ve never thought about reducing maxStep.

Yes. I’ve tormented by this question for a long time. I thought there must be something wrong in my code.
And I was focused on drag force.

Though, I uploaded the sinex file for station positions here anyway.
SLRF2020_POS+VEL_2023.03.28.snx (168.0 KB)


You’re very welcome Rongwang, I’m really happy I could help !


Note that the problem you experienced may come from the fact that our tutorials could be configured better. If you’re feeling like opening an issue on that matters you’re welcome.


Hi @MaximeJ ,

Actually, the key point is just the srp force.

@Bryan has pointed that in previous post.

Also, the analysis centers producing Lageos2’s orbit products generally not estimate the reflection coefficient. The value for Lageos2 is well known and equal to 1.13. The reason why we don’t active the solar radiation pressure for the Lageos2 tutorial is because it increases the computation a lot (i.e., ~16sec without SRP, ~69sec with SRP, and ~98sec with estimated SRP) without big improvements in the orbit solution.

While, in v12.0, the performance of srp is improved, so I think that it is acceptable to add srp force.