Orbit Determination of GEO satellite

Orekit Python Wrapper usage
1

Hi,
I’m trying to implement an OD procedure through batch LSE for a real GEO measurements.
Here the pseudocode
`

Set time scale to UTC

utc = TimeScalesFactory.getUTC()

Define reference system and earth, moon and sun model

inertialFrame = FramesFactory.getEME2000()   
teme = FramesFactory.getTEME()
gcrf = FramesFactory.getGCRF()
itrf = FramesFactory.getITRF(ITRFVersion.ITRF_2020,IERSConventions.IERS_2010, False) # boolean if true, tidal effects are ignored when interpolating EOP
eci  = gcrf
ecef = itrf
wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(ecef)
moon  = CelestialBodyFactory.getMoon()
sun   = CelestialBodyFactory.getSun()
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf)

Read TDM files

raw_tdm_file_path  = "/home/tanozar/Research/Orbit_Determinartion/TDM_prova_OD/TDM_38091_2023_02_14_new.kvn"
tdm = createTdm(raw_tdm_file_path)

Read TLE

tleLine1 = "1 38091U 12008A   23045.13109083  .00000066  00000-0  00000-0 0  9993"
tleLine2 = "2 38091   2.1185  71.6844 0013244 308.9480 229.1262  1.00270740 40252"
tleReference = TLE(tleLine1, tleLine2)

name_obs     = TdmMetadata.cast_(Segment.cast_(tdm.getSegments().get(0)).getMetadata()).getParticipants().get(1)
sat_name     = TdmMetadata.cast_(Segment.cast_(tdm.getSegments().get(0)).getMetadata()).getParticipants().get(2)
startDateTdm = TdmMetadata.cast_(Segment.cast_(tdm.getSegments().get(0)).getMetadata()).getStartTime()
finalDateTdm = TdmMetadata.cast_(Segment.cast_(tdm.getSegments().get(tdm.segments.size()-1)).getMetadata()).getStopTime()

Set observatory location on Orekit Frames

longitude = radians(41.765278)
latitude  = radians(13.375)
altitude  = float(550.0)

station_frame = TopocentricFrame(earth, GeodeticPoint(latitude, longitude, altitude), name_obs)
station = GroundStation(station_frame)
station.getPrimeMeridianOffsetDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH)
station.getPolarOffsetXDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH)
station.getPolarOffsetYDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH)

Retrieve measurements from List Observvation

angle_1_ra, angle_2_dec, epoch, epoch_dt= [], [], [], []

listSegments_tot = ArrayList()
listSegments     = tdm.getSegments()
listSegments_tot.addAll(listSegments)
listObservationsBlock = [ObservationsBlock.cast_(Segment.cast_(segment).getData()) for segment in listSegments_tot]

for obsBlock in listObservationsBlock:
    angle_1_ra.append(obsBlock.getObservations().get(0).getMeasurement())
    angle_2_dec.append(obsBlock.getObservations().get(1).getMeasurement())
    epoch_dt.append(absolutedate_to_datetime(obsBlock.getObservations().get(0).getEpoch()))
    epoch.append(obsBlock.getObservations().get(0).getEpoch())

SGP4 TLE Propagator

nadirPointing           = NadirPointing(eci, wgs84Ellipsoid)
sgp4Propagator          = SGP4(tleReference, nadirPointing, 2000.)
tleInitialState         = sgp4Propagator.getInitialState()                       # tle initial state
tleEpoch                = tleInitialState.getDate()                              # tle release epoch
tleOrbit_TEME           = tleInitialState.getOrbit()                             # orbit from tle in teme ref
tlePV_ECI               = tleOrbit_TEME.getPVCoordinates(epoch[0],eci)           # state in eci ref at first measurement time
tleOrbit_ECI            = CartesianOrbit(tlePV_ECI, eci, wgs84Ellipsoid.getGM()) # cartesian orbit from eci state
tleOrbit_keplerian_conv = OrbitType.KEPLERIAN.convertType(tleOrbit_ECI)          # conversion in keplerian orbit
tleOrbit_keplerian      = KeplerianOrbit(tleOrbit_keplerian_conv.getPVCoordinates(), inertialFrame, epoch[0], Constants.EIGEN5C_EARTH_MU)
tle_at_epoch = TLE(tleReference.getSatelliteNumber(), tleReference.getClassification(), tleReference.getLaunchYear(), 
    tleReference.getLaunchNumber(), tleReference.getLaunchPiece(), tleReference.getEphemerisType(), tleReference.getElementNumber(), 
    tleOrbit_keplerian.date, tleOrbit_keplerian.getKeplerianMeanMotion(), tleReference.getMeanMotionFirstDerivative(), tleReference.getMeanMotionSecondDerivative(), 
    tleOrbit_keplerian.e, tleOrbit_keplerian.i, tleOrbit_keplerian.perigeeArgument, tleOrbit_keplerian.rightAscensionOfAscendingNode, tleOrbit_keplerian.meanAnomaly, 
    tleReference.getRevolutionNumberAtEpoch(), tleReference.getBStar(), TimeScalesFactory.getUTC())

Read ephemeris and interpolate to measurements time

pod1 = '/home/tanozar/Research/Orbit_Determinartion/Effemeridi/ephem_38091.SP3'#/2023_02_14/WUM0MGXFIN_20230450000_01D_05M_ORB.SP3'
eph1 = ephem_read(pod1,'PC05',epoch_dt,station)
ra_ephem = eph1.loc[:,'Ra'].values
dec_ephem = eph1.loc[:,'Dec'].values

INITIAL ORBIT DETERMINATION (CIRCULAR HYPOTESIS)

raDec1  = AngularRaDec(station, inertialFrame, epoch[0], [radians(angle_1_ra[0]),radians(angle_2_dec[0])], JArray_double([1e-5,1e-5]), JArray_double([1.,1.]), ObservableSatellite(0))
raDec2  = AngularRaDec(station, inertialFrame, epoch[1], [radians(angle_1_ra[1]),radians(angle_2_dec[1])], JArray_double([1e-5,1e-5]), JArray_double([1.,1.]), ObservableSatellite(0))

iod_laplace_cartesian_circ_orbit = IOD_2points(raDec1,raDec2,station,gcrf,wgs84Ellipsoid,inertialFrame)

# ) Keplerian orbit (Laplace circular)
iod_laplace_keplerian_circ_orbit = KeplerianOrbit(iod_laplace_cartesian_circ_orbit.getPVCoordinates(),inertialFrame,raDec1.getDate(),Constants.EIGEN5C_EARTH_MU)

# ) Laplace circular TLE
iod_laplace_circ_TLE = TLE(tleReference.getSatelliteNumber(), tleReference.getClassification(), tleReference.getLaunchYear(), 
        tleReference.getLaunchNumber(), tleReference.getLaunchPiece(), tleReference.getEphemerisType(), tleReference.getElementNumber(), 
        iod_laplace_keplerian_circ_orbit.date, iod_laplace_keplerian_circ_orbit.getKeplerianMeanMotion(), tleReference.getMeanMotionFirstDerivative(), tleReference.getMeanMotionSecondDerivative(), 
        iod_laplace_keplerian_circ_orbit.e, iod_laplace_keplerian_circ_orbit.i, iod_laplace_keplerian_circ_orbit.perigeeArgument, iod_laplace_keplerian_circ_orbit.rightAscensionOfAscendingNode, iod_laplace_keplerian_circ_orbit.meanAnomaly, 
        tleReference.getRevolutionNumberAtEpoch(), tleReference.getBStar(), TimeScalesFactory.getUTC())

ORBIT DETERMINATION

# ) Estimator parameters
estimator_position_scale = 100.0 # m
estimator_convergence_thres = 1e-3
estimator_max_iterations = 50
estimator_max_evaluations = 35

# ) Orbit propagator parameters
prop_min_step = 0.001 # s
prop_max_step = 30.0 # sdd
prop_position_error = 10.0 # m

integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
#integratorBuilder = DormandPrince54IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)
propagatorBuilder_tmp = NumericalPropagatorBuilder(iod_laplace_cartesian_circ_orbit, integratorBuilder, PositionAngle.MEAN, estimator_position_scale)
propagatorBuilder_tmp.setAttitudeProvider(nadirPointing)
propagatorBuilder_tmp.setMass(2000.0) #mass of satellite

# ) Add Perturbations
propagatorBuilder = addPerturbations(propagatorBuilder_tmp, ecef, moon, sun, wgs84Ellipsoid)

# ) Set estimator properties
matrix_decomposer = QRDecomposer(1e-11)
#optimizer = GaussNewtonOptimizer(matrix_decomposer, False)
optimizer = LevenbergMarquardtOptimizer()
estimator = BatchLSEstimator(optimizer, propagatorBuilder)
estimator.setParametersConvergenceThreshold(estimator_convergence_thres)
estimator.setMaxIterations(estimator_max_iterations)
estimator.setMaxEvaluations(estimator_max_evaluations)

# ) OD parameters
weight = JArray_double([1., 1.])
sigma = JArray_double([4e-5, 4e-5])

# ) Add measurements to batch estimator
for i in range(len(angle_1_ra)):
    currentRaDec = [radians(angle_1_ra[i]), radians(angle_2_dec[i])]
    measurement  = AngularRaDec(station,inertialFrame, epoch[i], currentRaDec, sigma, weight, ObservableSatellite(0))
    estimator.addMeasurement(measurement)
    
# ) Estimator
estimator.setObserver(my_observer())
estimatedPropagatorArray = estimator.estimate()
estimatedPropagator = estimatedPropagatorArray[0]
estimatedInitialState = estimatedPropagator.getInitialState()

# ) Retrieve keplerian orbit
batch_orb = OrbitType.KEPLERIAN.convertType(estimatedInitialState.getOrbit())
batch_keplerian_orbit = KeplerianOrbit(batch_orb.getPVCoordinates(), inertialFrame, epoch[0], Constants.EIGEN5C_EARTH_MU)
batch_tle =  TLE(tleReference.getSatelliteNumber(), tleReference.getClassification(), tleReference.getLaunchYear(), 
    tleReference.getLaunchNumber(), tleReference.getLaunchPiece(), tleReference.getEphemerisType(), tleReference.getElementNumber(), 
    batch_keplerian_orbit.date, batch_keplerian_orbit.getKeplerianMeanMotion(), tleReference.getMeanMotionFirstDerivative(), tleReference.getMeanMotionSecondDerivative(), 
    batch_keplerian_orbit.e, batch_keplerian_orbit.i, batch_keplerian_orbit.perigeeArgument, batch_keplerian_orbit.rightAscensionOfAscendingNode, batch_keplerian_orbit.meanAnomaly, 
    tleReference.getRevolutionNumberAtEpoch(), tleReference.getBStar(), TimeScalesFactory.getUTC())

Propagate and compare real measurements and TLE wrt the precise ephemeris

ref_tle = TLE(tleLine1, tleLine2)
new_tle = TLE(batch_tle.getLine1(), batch_tle.getLine2())
sgp4Propagator1 = SGP4(ref_tle, nadirPointing, 2000.)
sgp4Propagator2 = SGP4(new_tle, nadirPointing, 2000.)

# start_time_dt = absolutedate_to_datetime(startDateTdm)
# end_time_dt   = absolutedate_to_datetime(startDateTdm.shiftedBy(3600.0 *24*2))
# time_n        = np.arange(start_time_dt, end_time_dt, dt.timedelta(seconds=60)).astype(dt.datetime)

norm1,norm2, norm3 = [], [], []
x1, y1, z1, x2, y2, z2 = [], [], [], [], [], []
ra1, ra2, ra3, dec1, dec2, dec3 = [], [], [], [], [], []

for i in range(len(epoch_dt)):
#for i in range(int(len(time_n))):
    
    time_n_tmp = epoch_dt[i]
    #sec = float(time_n_tmp.second) + float(time_n_tmp.microsecond)/1e6
    startDate = datetime_to_absolutedate(time_n_tmp)#AbsoluteDate(int(time_n_tmp.year), int(time_n_tmp.month), int(time_n_tmp.day), int(time_n_tmp.hour), int(time_n_tmp.minute), sec, TimeScalesFactory.getUTC())

    obs_state = station_frame.getPVCoordinates(startDate, inertialFrame)
    obs_tmp = obs_state.getPosition()

    # Ref TLE state
    sc_state1 = sgp4Propagator1.getPVCoordinates(startDate, inertialFrame)
    pos_tmp1 = sc_state1.getPosition() 
    x1.append(pos_tmp1.x) 
    y1.append(pos_tmp1.y)
    z1.append(pos_tmp1.z)
    norm1.append(np.linalg.norm(pos_tmp1.toArray()))
    ra_tmp1, dec_tmp1 = rECI2TopoRaDec(pos_tmp1, obs_tmp) 
    ra1.append(ra_tmp1)
    dec1.append(dec_tmp1)
    
    # OD state
    sc_state2 = sgp4Propagator2.getPVCoordinates(startDate, inertialFrame)
    pos_tmp2 = sc_state2.getPosition()  
    norm2.append(np.linalg.norm(pos_tmp2.toArray()))
    x2.append(pos_tmp2.x) 
    y2.append(pos_tmp2.y)
    z2.append(pos_tmp2.z)
    ra_tmp2, dec_tmp2 = rECI2TopoRaDec(pos_tmp2, obs_tmp) 
    ra2.append(ra_tmp2)
    dec2.append(dec_tmp2)

    # sc_state3 = sgp4Propagator3.getPVCoordinates(startDate, inertialFrame)
    # pos_tmp3 = sc_state3.getPosition()  
    # norm3.append(linalg.norm(pos_tmp3.toArray()))
    # ra_tmp3, dec_tmp3 = transform.rECI2TopoRaDec(pos_tmp3, obs_tmp) 
    # ra3.append(ra_tmp3)
    # dec3.append(dec_tmp3)


####################################
############### PLOT ###############
####################################

err_ra_tle_ephem  = [np.abs(i-j) for i,j in zip(ra1, ra_ephem)]
err_dec__tle_ephem = [np.abs(i-j) for i,j in zip(dec1, dec_ephem)]
err_ra_OD_ephem  = [np.abs(i-j) for i,j in zip(ra_ephem, ra2)]
err_dec_OD_ephem = [np.abs(i-j) for i,j in zip(dec_ephem, dec2)]
plt.figure()
plt.subplot(211)
plt.plot(epoch_dt,err_ra_tle_ephem,'r*',label='RA error TLE/Ephem')
plt.plot(epoch_dt,err_ra_OD_ephem,'g*',label='RA error OD/Ephem')
plt.legend()
plt.grid('minor')
plt.xlabel('DATE [UTC]')
plt.ylabel('Angular error [deg]')
plt.subplot(212)
plt.plot(epoch_dt,err_dec__tle_ephem,'r*',label='Dec error TLE/Ephem')
plt.plot(epoch_dt,err_dec_OD_ephem,'g*',label='Dec error OD/Ephem')
plt.legend()
plt.grid('minor')
plt.xlabel('DATE [UTC]')
plt.ylabel('Angular error [deg]')
plt.show()`

These are the TLE that I retrieve from IOD, OD compared wrt the reference TLE

TLE reference
1 38091U 12008A 23045.79524248 .00000066 00000-0 00000-0 0 9999
2 38091 2.2394 72.4721 0012948 306.8422 109.8802 1.00266834 40254
TLE IOD
1 38091U 12008A 23045.79524248 .00000066 00000-0 00000-0 0 9999
2 38091 1.9532 185.4967 0146323 104.0013 202.2751 1.03092309 40258
TLE OD
1 38091U 12008A 23045.79524248 .00000066 00000-0 00000-0 0 9999
2 38091 1.4439 253.1767 0013068 139.2320 98.6334 1.00267526 40256

The figure below shows the RA/DEC error comparison between TLE, OD and precise ephemeris

angular_error
angular_error871×593 32.7 KB

It’s clear that there is something wrong but I’m not sure where the error, or better the errorS are. I’m quite sure that the real measurements are accurate.
Do you have any suggestions?

Cheers
Gaetano

2

Hey Gaetano, that looks like it might be a coordinate system issue to me. When you’re dealing with TLEs you need to be especially careful.

As a tip, you should always start by plotting your measurements against the expected measurements from “ground truth” (here, your TLEs would be OK), to make sure that you understand your data. If the obs and the TLE don’t really line up, then you have to figure out why. Once you know your data, then you can try OD.

3

Hi @markrutten
What you’re saying makes sense.
I’ll try it.
As regard coordinates system, I think that the once used are consistent but I’ll check again

Thanks

Cheers
Gaetano

4

Hi @Tanozar,

I think there’s a misunderstanding when you’re doing this:

Here you’re using osculating values from your estimated orbit as input of a TLE.
TLE values are mean elements in the sense of SGP4/SDP4 theory, they’re not osculating elements.
That may be one of your problems.

When I look at your figure I can’t help but think that there is either a frame issue (like Mark said) or a radian/degree conversion issue somewhere. I would do some extra checks on that.
For example, are you sure that the parsed TDM values are in degree as seem to suggest the lines:

I’m asking because the default behavior of CCSDS ParserBuilder in Orekit is to convert compatible units while reading. In your case, it would convert the angles in the file in degree to radians in the output TDM object.

Hope this helps,
Maxime