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 AbstractOrbitDetermination.java 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
.getLastEstimations();
// 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)) {
@SuppressWarnings("unchecked")
final EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) estimatedMeasurement;
rangeLog.add(evaluation);
// 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)) {
@SuppressWarnings("unchecked")
final EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) estimatedMeasurement;
rangeRateLog.add(evaluation);
} else if (measurementType.equals(AngularAzEl.MEASUREMENT_TYPE)) {
@SuppressWarnings("unchecked")
final EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) estimatedMeasurement;
azimuthLog.add(evaluation);
elevationLog.add(evaluation);
} else if (measurementType.equals(PV.MEASUREMENT_TYPE)) {
@SuppressWarnings("unchecked")
final EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) estimatedMeasurement;
positionLog.add(evaluation);
velocityLog.add(evaluation);
}
}
(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
continue;
}
AbstractOrbitDetermination.java (133.8 KB)
odin_stella_20230101T00.yaml (5.6 KB)
stella_202301.np2 (69.7 KB)
Best regards,
Rongwang