Estimate usage? (from ObservedMeasurement)

Orekit MATLAB usage
1

I am having problems to understand how this method works, as it is throwing errors at evaluation (I am using Orekit in Matlab):

Incorrect number or types of inputs or outputs for function 'estimate'.

From looking at the documentation I see it needs three things:

estimate(int iteration, int evaluation, SpacecraftState[] states)

I am not sure what the first two inputs are, why are there iterations and evaluations involved? Is it resampling the measurement using the noise generation part of the measurement action?

But I don’t have a problem with those inputs, given that they are just integers (I am using 0 for both). The problem is that I am not sure how I am supposed to get the states part. The measurements come from a Generator configured to simulate a radar, so in the process I have to create a Propagator object, that I .addPropagator to it.

After that, when I already have the measurement I don’t have easy access to the spacecraft state at the instant of the given measurement, so what I am doing right now is to use the same propagator, and .propagate to the date of the measurement, thus giving the supposed spacecraft state at that moment. This seems to me a bad way of doing it as the propagation was already performed, and I don’t want to re-do that, but I can’t see another way of accesing the state.

So by doing that I still get the error mentioned, and I don’t understand what I am doing wrong. Reading the documentation it says that states is:

`states` - orbital states corresponding to `getSatellites()` at measurement date

But doing getSatellites to the measurement gives an ObservableSatellite object, which from looking at the documentation doesn’t have information of the spaceraft state in the methods available, only the index of the propagator. So again, I don’t see any other way of geting the states info, aside from re-computing the propagation with the propagator…

My main question is why is the estimate method giving me the error, and if possible I would like to understand if I am geting the state wrongly.

These are the code lines in question:

state = myPropagator.propagate(myMeasurements(i).getDate);

EstimationObject = myMeasurements(i).estimate(0,0,state);

Where myPropagator is the same propagator object that I fed into the measurement Generation object.

Thanks in advance, and sorry if there are not enough details here, as I am writing it in a rush. If needed, this is a script with the complete example of what I am trying to do:

clear all; clc; close all;
% TEST

import org.orekit.orbits.*
import org.orekit.frames.*
import org.orekit.time.*
import org.orekit.utils.*
import org.orekit.propagation.analytical.*
import org.orekit.estimation.*
import org.hipparchus.geometry.euclidean.threed.*
import org.orekit.bodies.*
import org.orekit.geometry.fov.*
import org.orekit.propagation.events.*

% IC
x0  = 1e6*[1.459975626800790   0.436989307528660  -6.916264503127049  -0.003895232294893  -0.006282079270069  -0.001219414996904];
t0  = datetime(2022,05,01);
at0 = AbsoluteDate(2022, 05, 01, 0, 0, 0, TimeScalesFactory.getUTC());
% endDate = at0.shiftedBy(2*86400); 
endDate = at0.shiftedBy(10*86400);
% orbit

mu              = Constants.WGS84_EARTH_MU;
inertialFrame   = FramesFactory.getEME2000();
myOrbit         = EquinoctialOrbit(PVCoordinates(Vector3D(x0(1:3)),Vector3D(x0(4:6))), inertialFrame, at0, mu);
myPropagator    = KeplerianPropagator(myOrbit);

% Radar
ITRF            = FramesFactory.getITRF(ITRFVersion.ITRF_2000,IERSConventions.IERS_2010, 1);
earth           = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ITRF);

sensor          = GeodeticPoint(0.64867, -0.0976, 142.33);
stationFrame    = TopocentricFrame(earth, sensor,'myRadar');
myRadar         = measurements.GroundStation(stationFrame);
% fov
center  = [0 0 1]';
ax1     = [1 0 0]';     hA1 = 35*pi/180;
ax2     = [0 1 0]';     hA2 = 35*pi/180;

% radar_FOV = DoubleDihedraFieldOfView(Vector3D(center),Vector3D(ax1),hA1,Vector3D(ax2),hA2,0);
radar_FOV = CircularFieldOfView(Vector3D(center),35*pi/180,0);

% Create the GENERATOR object for the measurements
myGenerator     = measurements.generation.Generator();
mySat           = myGenerator.addPropagator(myPropagator);

% AzEl Builder
sigmas = [0.005235 0.005235];
myAzElBuilder   = measurements.generation.AngularAzElBuilder(null(1), myRadar,sigmas, 1, mySat);

% FixedStepSelector
updateTime          = 7;
fixedStepSelector   = FixedStepSelector(updateTime,null(1));

% EventDetector
handler             = handlers.ContinueOnEvent();

% FOV DETECTOR
% fov_detector        = ElevationDetector(100,1e-8,stationFrame).withConstantElevation(55*pi/180).withHandler(handler);
fov_detector        = GroundFieldOfViewDetector(stationFrame,radar_FOV).withHandler(handler).withMaxCheck(100).withMaxIter(100).withThreshold(1e-8);

fov_detector = fov_detector.withMaxCheck(20).withMaxIter(200).withThreshold(1e-8);


logger              = org.orekit.propagation.events.EventsLogger();
fov_detector        = logger.monitorDetector(fov_detector);

% SCHEDULER
signSemantic        = measurements.generation.SignSemantic.valueOf('FEASIBLE_MEASUREMENT_WHEN_NEGATIVE');
% signSemantic        = measurements.generation.SignSemantic.valueOf('FEASIBLE_MEASUREMENT_WHEN_POSITIVE');
myScheduler         = measurements.generation.EventBasedScheduler(myAzElBuilder, fixedStepSelector, myPropagator, fov_detector, signSemantic);
    
myGenerator.addScheduler(myScheduler);

tic
myMeasurements = myGenerator.generate(at0,endDate);
toc

myMeasurements  = myMeasurements.toArray();
n_mes           = length(myMeasurements);
fprintf('Number of measurements: %i \n',n_mes)

obs_val = zeros(2,n_mes);
dates_radar = datetime.empty(0,n_mes);

for i = 1 : n_mes
    
    azel = myMeasurements(i).getObservedValue;

    dates_radar(i)  = date_conversion(myMeasurements(i).getDate,'orekit2matlab');
    obs_val(:,i) = azel;

    state = myPropagator.propagate(myMeasurements(i).getDate);
        
    EstimationObject = myMeasurements(i).estimate(0,0,state);

end

figure(3); 
polarplot(obs_val(1,:),90-180/pi*obs_val(2,:),'bo'); hold on;
polarplot(ObVal_az_el(1,:),90-180/pi*ObVal_az_el(2,:),'r*')
rlim([0 90]);

myEvents = logger.getLoggedEvents.toArray();
for i = 1 : length(myEvents)
    E = myEvents(i);

    figure(3); hold on
    el = stationFrame.getElevation(E.getState.getPVCoordinates().getPosition(), inertialFrame, E.getDate);
    az = stationFrame.getAzimuth(E.getState.getPVCoordinates().getPosition(), inertialFrame, E.getDate);
    polarplot(az,90-180/pi*el,'m*','MarkerSize',6)
end
2

They are only used as an indication of the progress of the orbit determination. They are generally just passed along in the EstimatedMeasurement that will be created. This is used for example in the {Dynamic}OutlierFilter classes (which are measurements modifiers) that have a warmup index and do not operate on the first few iterations, they just kick in after warmup.

The state will be provided by the estimator itself when doing orbit determination.

If you are using for example batch least square, the state you will get for the first iteration will correspond to the initial guess propagated. Then the least square engine will attempt to modify the estimated parameters like initial orbit and force model parameters, and then run another iteration. This new iteration will therefore have a different state as the orbit has changed. Then the new estimate will change, and hopefully the estimated measurement should be closer to the observed measurement. This continue until convergence and at the end you get orbits that is closest to the set of observed measurements.

If you are using Kalman filter, then the iteration and evaluation are not really meaningful since this kind of estimator does not performs several passes, it just progress always forward in time in one sweep. The states however have similar semantics as batch least square for the last iteration: they correspond to the best estimation of the state the estimator has been able to find up to call point.

If you are generating measurements, then the estimate method is called automatically by the generator, using the initial orbit and force model parameters (that never change in case of measurement generation).

The estimate method is not really intended to be called directly by users. It is designed to be called on the fly by either least square estimator, Kalman filter, or measurement generator and the calling algorithm knows about the current propagators used.

As I understand your code, you first generate the observed measurements in a first pass, and then loop over them in a second pass, calling the estimate method to get the estimated value. So if you did not change the orbit between generation and call to estimate (i.e. you call it directly, not from an orbit determination engine), then the estimated measurement will be the same as the observed measurement plus the noise, but you seem to set up null noise, so they should be perfectly equal in your case.

1 Like
3

Thank you @luc for the fast response. I have been using Orekit for the past two years for my PhD, and between the fantastic documentation and the activity of the forum learn using it is really fast. Still, I have a lot more to learn from the complete toolkit…

Now it makes much more sense. I could not see an easy way to leverage that method. So if it is only meant for inernal use in the orbit determination tools I can understand why it is there.

I’ll explain what I eventually want to compute, so it is easier to understand what I was trying to do. As my PhD is focused in maneuver detection I usually test the metrics by doing my own “real” satellites simulation (with high fidelity modeling) and radar measurements from these simulations. Until now I was using my own implementation of radar pass detection. I did that by leveraging the event detection capabilities of Orekit (using Elevation masks). And aside from the (unnecesary) complexity of the internal logic to register the observables only when inside the mask (at regular time intervals), I was doing a naive aproach to measurement, in that I just got .getRange (or range rate) at the instant of measurement. I now see that this was naive because the instantaneous range is not exactly what a ranging radar sees when it sends the signal and measures the time of return.

So, despite knowing the power of the already implemented tooling of Orekit I wasn’t really using it, with an overly complex radar simulator, and a flawed one to top it. Now I was in the process of completelly changing the workflow of my metric testing, to better use the tools orekit provides, like the measurement generator. This is much more straighforward and easy to use (you just give it the sensor location, FOV, and the scheduler for the observables). What I was trying to do and show in the post example relates to how I test the metrics. I not only simulate the “real” orbit and just do a measurement. What I do is sample each measurement a number of times (say 300), so that I get not just one, but many possible tracks of the same observation event. The reason being that I am using metrics which are distributed in a non trivial way. For example if you use the Mahalanobis Distance to check the statistical distance beetween two states, you are asuming that both are normaly distributed, and thus you know that this distance follows a chi squared distribution. The metric I am testing does not asume the predicted orbit distribution is normally distributed, and it makes the distribution of the metric something not obvious, so I sample measurements to get the actual distribution of a particular case (which is characterized by the maneuver, or lack of it).

Knowing the final goal, it is maybe easier to see that what I was trying to do was to re-sample the measurement or observation action. I thought that estimate was for that. Also, because I was not sure of that, my current approach would be to use a null noise, as seen in the example, to get the perfect/ideal measurement, and then add the noise by myself, so that I can compute the distribution of the detection metric.

The question would then be: once the generator is done, Can I re-do the measurement easily? Ideally to get the same ObservedMeasurement but with a diferent ObservedValue…

Thanks again for the detailed response.

4

Hey there,

I might have missed your point, but can’t you just use the build method of RangeBuilder for instance on a bunch of SpacecraftState, possibly all duplicates of the same guy?

Best,
Romain

5

Hi Serrof. I guess that would be another approach to generate the measurement samples, but still would require to do something previous to know when I have to meassure. So either run a Generator once (with the event detection configured as needed for the station/s), or do as I have been doing for now with my own pass detector (but not just geting the range with .getRange, and instead using the RangeBuilder for realism).