Angle oly IOD with simulated measurements getting horrible results

I’ve been building playing with a sequential UKF and generating synthetic angle only measurement for sample data. The filter seems to work great which gives me confidence that the measurements are being generated correctly. I wanted to produce the initial guess via IOD and the results have been absolutely horrible, even for long arcs.

Here is how I create measurements

def simluate_orbit_measurements(  # type: ignore
    reference_orbit: Orbit,
    topo_frame: TopocentricFrame,
    t0: AbsoluteDate,
    tf: AbsoluteDate,
    step: float = 1.0,
    az_cov: float = 0.001,
    el_cov: float = 0.001,
    small: float = 0.00001,
    sigma_az: float = 0.001,
    sigma_el: float = 0.001,
):
    """
    Create simulated measurements.

    Args:
        reference_orbit (Orbit):        reference orbit
        topo_frame (TopocentricFrame):  topo frame of the observer
        t0 (AbsoluteDate):              meas period start time
        tf (AbsoluteDate):              meas period end time
        step (float):                   measurement time step
        az_cov (float):                 process noise covariance matrix index [0, 0].
                                        Azimuth and elevation are being treated as
                                        independent variable for simplicity’s sake so
                                        az_cov and el_cov are used to populate a
                                        diagonal matrix
        el_cov (float):                 process noise covariance matrix index [1, 1].
        small (float):                  Diagonal elements threshold under which column
                                        are considered to be dependent on previous ones
                                        and are discarded. This is  straight from
                                        Orekit definition
        sigma_az (float):
        sigma_el (float):

    Returns
    -------
        Set[AngularAzEl]
    """
    # create propagator for reference orbit
    prop_kep = KeplerianPropagator(reference_orbit, EARTH_MU)
    # create ground station that's measuring the simulated measurements
    ground_station: GroundStation = GroundStation(topo_frame)
    # create measurement generator
    meas_gen = Generator()  # type: ignore
    obs_sat = meas_gen.addPropagator(prop_kep)
    subscriber = GatheringSubscriber()  # type: ignore
    subscriber.init(t0, tf)
    meas_gen.addSubscriber(subscriber)

    # this adds noise to the measurements
    vec_gen = CorrelatedRandomVectorGenerator(
        DiagonalMatrix([az_cov, el_cov]),
        small,
        GaussianRandomGenerator(RandomDataGenerator()),
    )

    # give equal weights to az and el hence [1.0, 1.0] magic numbers. No need to make
    # this function infinitely configurable
    azel_builder = AngularAzElBuilder(
        vec_gen, ground_station, [sigma_az, sigma_el], [1.0, 1.0], obs_sat
    )

    # generate measurements only when visible, that way we don't get atmospheric
    # corrections erroring out or some other unforeseen behavior
    elevation_detector = ElevationDetector(topo_frame)
    event_scheduler = EventBasedScheduler(
        azel_builder,
        FixedStepSelector(step, UTC),
        prop_kep,
        elevation_detector,
        SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE,
    )

    meas_gen.addScheduler(event_scheduler)

    meas_gen.generate(t0, tf)
    measurements = subscriber.getGeneratedMeasurements()

    return measurements

This generates about 2000 measurements from which I pick three and try Gauss, Laplace, and Gooding IODs. I won’t go into my implementation of how I call these IOD methods because I’ve run them with other data and had great results (great for angles iod at least) so I want to point the finger to either an ill conditioned scenario or something to do with how I generated the measurements. The pass is rather low, with peak elevation angle < 40 deg

These are the results I get with 3 measurements that cover 400 seconds of separation

Reference Orbit:
Keplerian parameters: {a: 2.6871E7; e: 0.1; i: 55.0; pa: 0.0; raan: -20.0; v: 180.0;}

Gauss IOD
Keplerian parameters: {a: 4.237043009921982E7; e: 0.8878027310562274; i: 105.71985914321517; pa: 136.77671702324247; raan: -18.417826165868302; v: -68.4322943654454;}

Laplace IOD
Keplerian parameters: {a: 5.568142916102404E7; e: 0.9156744453115059; i: 105.49433163980962; pa: 137.15157457590072; raan: -18.964216651890524; v: -69.15886823807563;}

Gooding IOD
Keplerian parameters: {a: 3511316.136873923; e: 0.9998490378536055; i: 39.65962372833863; pa: 303.40046926246333; raan: 45.5261157625205; v: -179.99955679702893;}

I know eccentricities are tough with angles only but I have gotten much better results in the past with real data. Even the orbit plane is pretty bad here.

Is there anything that immediately jumps out to you?

Thank you in advance.

Two wild guesses:

  • make sure your topocentric frame is based on a real Earth frame (we have seen inertial/terrestrial confusions before)
  • try to pick up the measurements over a longer arc, 400s seems short to me for a Middle Earth Orbit

Hello,

One small correction maybe, but I doubt it will drastically change the results.
az_cov and el_cov should be (0.001)² otherwise you’re generating measurements with a 1.8° error.

Cheers,
Maxime

@MaximeJ just so I understand how the measurement simulation works:

The az_cov and el_cov that are used for the vector generator, do those represent the noise of the actual orbital motion and the simga_az and sigma_el represent the measurement noise when creating the AngularAzElBuilder?

azel_builder = AngularAzElBuilder(
        vec_gen, ground_station, [sigma_az, sigma_el], [1.0, 1.0], obs_sat
    )

So in Kalman terms the covariance matrix is related to Q and the az/el sigmas are related to R?

##############################################################

@luc
I agree with what you are saying about 400s being a short arc. My intent is to do this for LEO where the passes would be shorter but I made the orbit MEO so it would be easier to see above the horizon, it was due to laziness of not wanting to design the orbit to be visible from my sensor location. This actually generated 2000-2500 seconds of measurements but I actually got better IOD results using measurements from the middle 400 seconds rather than the full arc. I’m guessing these angles only IOD methods perform better at higher elevation angles, at least in terms of resolving the orbit plane.

I’m not quite sure what you mean by making sure the TopocentricFrame is based on a real Earth frame. This is how I built it and have been doing it this way for other things that seem to work fine. Does this look right?

ITRF = FramesFactory.getITRF(IERS_2010, True)
EARTH_ELLIPSOID = OneAxisEllipsoid(
    Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ITRF
)
station_gp = GeodeticPoint(radians(lat), radians(lon), alt)
station_topo = TopocentricFrame(EARTH_ELLIPSOID, station_gp, "station")

It is fine, your frame is correct as it is based on an ITRF.

  • az/el sigmas are related to R indeed, these are the known, or theoretical standard deviations of your measurement system.
    It’s used in batch least squares and Kalman to weight the measurements
  • The covariance matrix represents the noise you’re going to simulate in your MeasurementBuilder, so the noise actually applied to your measurements when you generate them. And this should contain sigmas squared, hence the (.)² I proposed in my previous message.
    But they have no link to the process noise matrix Q in Kalman.

OK sorry just to make sure I really get it.

The AzEl sigmas are the sigmas that get passed to the actual AngularAzEl constructor and is our best knowledge of measurement noise.

The covariance matrix is the stochastic noise actually being added to the geometric calculated measurements by the meas simulator? Bad analogy on my part with Q matrix.

So they both represent measurement noise, the former actually being directly applied to the R matrix via OD, and the latter being the noise actually applied by the simulator.

Yes that’s it

@MaximeJ @luc
I had some good results implementing a double-r IOD, I’ve only tested with a few simulated data sets that had low noise but it looks promising. As far as I know Orekit does not have anything similar. If I can find some time, would you guys be interest in a contribution? I could do python or java though my code is not going to be good as you guys write.

Here’s some sample results for 3 measurements across 400 seconds

Seed Orbit
Keplerian parameters: {a: 7371000.0; e: 0.1; i: 55.0; pa: 0.0; raan: 150.0; v: 42.29855896568464;}

IOD Orbit
Keplerian parameters: {a: 7372447.925164646; e: 0.09784620596101012; i: 54.97802124475939; pa: 1.2418487586662612; raan: 149.9627465987682; v: 54.27369646986524;}

Yes, we would be delighted to have another type of IOD. We prefer Java so it could be included into the core library. Don’t worry about code quality, we can deal with that while reviewing the contribution.

Thanks in advance!

I agree with Luc, we’re always delighted to have new features and new contributors!
Thanks!

Hi Luc,

Would there be any interest in having an IOD method that uses more than three measurements? Some techniques provide highly accurate and robust estimation. Typical IOD method for space based iod seems to have a lot of trouble converging in case of noisy measurement or complicated geometry.

Sure, all these kind of improvements are welcome!