Single track orbit determination analysis for low Earth orbit with approximated J2 dynamics

Single track orbit determination analysis for low Earth orbit with approximated J2 dynamics Jose M. Montilla a,⇑,1 , Jan A. Siminski b,2 , Rafael Vazquez a,1 a Universidad de Sevilla, Escuela Te ´cnica Superior de Ingenierı ´a, Camino de los Descubrimientos s.n, Sevilla 41092, Spain b ESA/ESOC, Robert-Bosch-Str 5, Darmstadt 64293, Germany Received 19 February 2024; received in revised form 10 September 2024; accepted 19 September 2024 Available online 25 September 2024 Abstract In the domain of Space Situational Awareness (SSA), the challenges related to orbit determination and catalog correlation are notably pronounced, exacerbated by data scarcity. This study introduces an initial orbit determination methodology that relies on data obtained from a single surveillance radar, with the need for fast algorithms within an operational context serving as the main design driver. The result is a linearized least-squares fitting procedure incorporating an analytically formulated approximation of the dynamics under the J2 perturbation, valid for short-term propagation. This algorithm utilizes all available observables, including range-rate, distinguishing it from other similar methods. A significant contribution of this paper is the enhancement of estimation quality by incorporating information about the object’s predicted orbital plane into the methodology, a method denoted as OPOD. The proposed methods are evaluated through a series of simulations against a classical range and angles fitting method (GTDS) to examine the effects of track length and measurement density on the quality of full state estimation, including the impact of using arcs that are too short. The OPOD methodology shows promising results throughout a wide range of scenarios. Ó2024 COSPAR. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/ by-nc-nd/4.0/). Keywords: Initial orbit determination; Maneuver detection; Space situational awareness 1. Introduction Over the last sixty years, human activities in space have significantly increased the number of orbiting objects (ESA, 2023), a concern heightened by the deployment of mega-constellations like Starlink and OneWeb. This congestion in low Earth orbit underscores the need for nations and space agencies to enhance tracking and cataloging efforts. Robust and efficient orbit determination (OD) algorithms are crucial in this context, playing a pivotal role in satellite tracking by swiftly identifying correlated tracks to maintain an accurate satellite database. The reliability of these algorithms is essential for Space Situational Awareness (SSA) and the safety of orbital operations. Classical OD strategies involve an Initial Orbit Determination (IOD) phase to generate and initial estimate, which is later refined (Schutz et al., 2004), encompassing methods like those of Laplace, Gauss, and Gooding (Escobal, 1970; Gooding, 1993). These methods require measurements with at least six independent parameters to generate a sixparameter orbit, with the quality of the initial estimation benefiting from measurements originating from different observation arcs, referred to as tracks 3 . Generally, multiple https://doi.org/10.1016/j.asr.2024.09.035 0273-1177/Ó2024 COSPAR. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ⇑ Corresponding author. E-mail addresses: [email protected] (J.M. Montilla), jan.siminski@esa. int (J.A. Siminski), [email protected] (R. Vazquez). 1 Aerospace Engineering and Fluid Mechanics Department, Spain. 2 Space Debris Specialist, Space Debris Office, ESOC, Germany. 3 A series of sequential observations delayed by seconds, which in turn are composed of a set of measurements (of different nature) at a common epoch. www.elsevier.com/locate/asr Available online at www.sciencedirect.com ScienceDirect Advances in Space Research 74 (2024) 4968–4989 observation arcs lead to a linkage problem (Milani et al., 2004), where information from different tracks is correlated for proper OD (Milani et al. (2005);Ma et al. (2018)). This study aims to develop a reliable IOD method for single track radar measurements in the LEO regime, facilitating the automatic cataloging of uncorrelated tracks and subsequent maneuver detection within the context of the Spanish Space Surveillance and Tracking Surveillance Radar (S3TSR) (Gomez et al., 2019), and addressing the challenge of short arcs where arc length might be insufficient for existing methods (Tao et al., 2020). A significant contribution of this study is incorporating orbital plane data into the IOD process to enhance estimation quality. The well-known range and angles method of the Goddard Trajectory Determination System (GTDS) (Long et al., 1989) is used as a reference for single track IOD in this work. However, this method is limited to inertial position information and cannot utilize the full set of observables provided by modern radars. Related to the method introduced here, Vallado and Carter (1998) apply a classical differential correction method for IOD with dense observation arcs and numerical propagation models. In contrast, this work focuses on leveraging analytical methods for scenarios with low-density radar tracks, where data scarcity is the primary challenge. The inclusion of extra information to aid in the estimation problem is an approach already present in other IOD methods proposed in the literature. The idea of initializing tracks with angles and angle-rates data, as explored in DeMars et al. (2012), is only suitable for optical-based sensors and not applicable to the LEO radar measurements used here. Other studies, such as Shang et al. (2018) and Zhang et al. (2019), involve multiple stations or very short tracks with extra radial measurements, which are beyond the scope of this work. For a more general review, see Pastor et al. (2021), which highlights the trade-off between speed and accuracy in combined Initial Orbit Determination and Orbit Determination methods. Building on prior investigations (Montilla et al., 2023a; Montilla et al., 2023b), this work seeks to utilize all information from a single radar track, including range, rangerate, and line-of-sight data (azimuth and elevation). Consequently, we establish an IOD process based on radaracquired data to estimate the satellite’s full state. The use of attributables (Montilla et al., 2023a; Reihs et al., 2021) is replaced by a procedure that includes dynamics information, similar to the GTDS range and angles method (Siminski, 2016; Vananti et al., 2017). GTDS remains relevant for IOD, and it is a good alternative to classical 3 points methods such as Herrick-Gibbs (Kaushik, 2016), but leveraging all available observables, including rangerate data, should lead to better results. The unexpected need to incorporate Earth’s flattening in short track orbital determination is a notable outcome explored in this study, given the relative importance of the J 2 perturbation compared to the main Keplerian term. For long-term non-Keplerian propagation, while approximated semianalytical (Amato et al., 2019) or analytical methods (Martinusi et al., 2015) could be used, these are not suited for short track IOD applications. Numerical methods for uncertainty propagation (Hernando-Ayuso et al., 2023) also do not meet the computational efficiency requirements for real-time track association (Pastor et al., 2021). The proposed method strikes a balance between error and speed, suitable for short propagations. The generalized equinoctial orbital elements in Bau `et al. (2021) are used due to their explicit inclusion of potential-derived perturbations and slow evolution, making them ideal for short-term propagation through a Taylor expansion of the solution. This study produces the analytical time derivatives of these orbital elements, up to the fourth order, and calculates the error state transition matrix through the derivation of the resulting polynomial’s coefficients. This is then used in a differential correction algorithm (Vallado, 2001) for IOD with the J2zonal harmonic included, along with an extended application involving past orbital data. Thus, the main innovation of this work is the development of a comprehensive state estimation algorithm. This algorithm employs a basic linear least-squares method to directly fit radar observations, with state and uncertainty calculations performed using an analytical (approximated) J 2 propagator. Furthermore, the method is enhanced by incorporating expected orbital plane information—referred to as OPOD—into the process. In the LEO regime, this orbital plane can be reliably predicted based on prior estimates, even when considering potential maneuvers. The newly proposed methodologies are evaluated against a widely recognized orbit determination approach, demonstrating significant improvements in IOD for simulated tracks, while preserving the analytical framework essential for practical use. The remainder of this paper is organized as follows. In Section 2, the problem and notation are defined. Section 3 describes the full state estimation algorithm first as a general linear least-squares method, then particularized to three different alternatives reliant only on data from a single radar track. This section includes an extra methodology that employs past orbital information to improve the estimation. Section 4compiles detailed testing of the estimation methods based on simulated tracks. All methods are then applied to various scenarios that introduce variability in track length and measurement frequency. The document is closed with conclusions in Section 5.Aincludes most of the formulation needed in the presented work. 2. Problem statement and notation This work seeks algorithms for IOD using a dynamical model to fit various satellite state-derived variables and radar observables from a single track (Fig. 1). It defines a radar station with known inertial position (Rr), velocity (Vr), and antenna orientation (given by the rotation matrix J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4969 TR I) for any required instant. The radar has a designated Field of Regard (FoR) characterized by a specific revisit time rt. When a Resident Space Object (RSO) passes through the sensor’s FoR, it generates a track composed of observations (also referred to as plots), which are sets of measurements taken at a common epoch tm. Here, m21;n½, where nis the number of individual plots. Typically, radar observables for LEO objects include range qm, azimuth Azm, elevation elm, and range-rate _ qm. Each measurement has an error characterized by a known standard deviation, included in a covariance matrix Czmfor the set, with potential correlations between different observables. Here the letter zdenotes a general observed data point, typically a radar observable, but it can be any state derived quantity used for fitting. The notation is modified to ^z when said quantity is computed from the state that is being estimated in the fitting process. It is assumed that the observation function rm¼OBSfun qm;Azm;elm ðÞmaps the radar observation of range and angles to an inertial position, which can then be treated as the measurement. This is needed for the GTDS method as it fits positions computed from the radar observables of range, azimuth and elevation. The covariance of such position measurement Crmis calculated from the original radar measurement uncertainty (Czm)usingthe Unscented Transform (UT) (Goff, 2015). The sigma points are reconstructed in inertial space to obtain the covariance matrix Crm, offering a more conservative uncertainty representation than the linear method (Siminski, 2016). At the midpoint t0¼t1þðtnt1Þ=2 of the radar track, the satellite’s state estimation ^ y0¼^ r0^ v0 ½and its covariance Cyare computed. Note that ^ r0and ^ v0are the estimated inertial position and velocity in the Cartesian frame I. This approach requires calculating predicted measurements, ^zm, from the estimate ^ y0, necessitating the propagation of the orbital state to the instant tm. The propagation step is encapsulated within ^ ym¼^ rm^ vm ½¼P ^ y0;Dtm ðÞ, where Dtm¼tmt0is the propagation duration from the reference instant t0, which is omitted 4 . The details of the propagation process will be discussed in Sections 3.1 and 3.3. The function z¼hy;tðÞis defined to convert an inertial state yat an instant tto any observable z, allowing the computation of the predicted measurement ^zm¼h^ ym;tm ðÞ. 3. Iterative linearized least-squares fitting algorithm This Section 3reviews a general linear least-squares algorithm for any observable and propagation model. Section 3.1 adapts this algorithm to the GTDS, fitting position measurements with Keplerian dynamics. In Section 3.2, radar observables are fitted using a Kepler model, enhanced by a Taylor expansion for a more accurate computation of the derivatives. Section 3.3 applies generalized equinoctial orbital elements and a Taylor expansion for the inclusion of J2dynamics in short-track IOD. Finally, Section 3.4 proposes incorporating orbital plane information to enhance the orbit determination algorithm. To maintain a broad scope, the fitting algorithm is derived from a general scalar observable zmat time tmwith known standard deviation rm. The predicted measurement ^zm¼h^ ym;tm ðÞgenerally differs from the estimated state ^ y0, even without considering model inaccuracies in ^ ym¼P^ y0;Dtm ðÞ. Thus, zm¼^zmþnm, where nmis the residual. The goal is to minimize Pn m¼1ðwmnmÞ2, where wm¼1=rm: Jz¼min ^ y0X n m¼1ðwmnmÞ2:ð1Þ A prerequisite is that the number of measurements must exceed the state dimension (n>6). This leads to a nonlinear minimization problem, typically solved iteratively for efficiency. Assuming a good initial estimate of ^ y0, the predicted observable ^zmis approximated at iteration ias a linearization around ^zðiÞ mfor a variation D^ yðiÞ 0: hð^ ym;tmÞhP^ yðiÞ 0;Dtm  ;tm  þ@h @^ y0ð^ yðiÞ 0;tmÞ D^ yðiÞ 0 ¼hð^ yðiÞ m;tmÞþ@h @yð^ yðiÞ m;tmÞ @P @^ y0ð^ yðiÞ 0;DtmÞ D^ yðiÞ 0:ð2Þ Considering all measurements z¼z1;;zn ½ |, and the predicted measurements ^ zðiÞ¼^zðiÞ 1;;^zðiÞ n hi |from ^ yðiÞ 0, the error DzðiÞ¼z^ zðiÞis defined. The derivatives of the predicted value for the last iteration @h=@^ y0cð^ yðiÞ 0;tmÞare stacked in the matrix A to get the linear equation DzðiÞ¼ADyðiÞ 0þn. With this expression, the same minimization problem as in Eq. 1can be solved as a linear least-squares problem: Fig. 1. The estimation problem from a single radar track is shown as iterative fitting of radar observables: range, azimuth, elevation, and rangerate qm;Azm;elm;_ qm ðÞ, creating a track of n observations from the same Resident Space Object. 4 The propagation generally requires to specify the reference epoch t0, but potential-based dynamics that go up to the second zonal harmonic does not depend on the reference instant, so only the total propagation time Dtmis included in the formulation for brevity. J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4970 DzðiÞ¼AD^ yðiÞ 0þn; Jz¼min D^ yðiÞ 0 n|Wn¼min D^ yðiÞ 0ðDzðiÞAD^ yðiÞ 0Þ|WðDzðiÞAD^ yðiÞ 0Þ; D^ yðiÞ 0¼A|WAðÞ 1ðA|WDzðiÞÞ; ð3Þ where Cy¼A|WAðÞ 1is the estimation covariance. Residual normalization is required, with W ¼C1 z(inverted measurement covariance) and Cz¼diag r2 1;;r2 n  , for example. This differential correction approach simplifies the problem to iterations of a linear least-squares problem. While general convergence is not assured due to potential local minima, a close initial estimate offers a favorable speed-accuracy trade-off (Vallado, 2001). Iterations continue until the position displacement norm jD^ rðiÞ 0jis below a threshold. The propagation method Pð^ y0;DtmÞand the relationship between spacecraft state and measurements in hðy;tÞare particularized next with three analytical methods. 3.1. Keplerian fitting of non-weighted Cartesian position - GTDS The GTDS range and angles method (Siminski, 2016; Vananti et al., 2017) is a specific case of the general algorithm from Section 3. This method fits measurements in inertial position, represented by rmat each instant tm,so the vector of all measurements is z¼r| 1;r| 2;;r| n  |. Predicted measurements ^ rmare an explicit function of the propagated estimation ^ ym, given by ^ rm¼h^ ym;tm ðÞ. The propagation function P^ y0;DtðÞensures Keplerian dynamics. This is achieved using the universal variable v, which is applicable to all orbit types and allows computation of the fand gKepler functions, yielding the propagated state as a linear combination of the initial conditions: ^ ym¼P^ y0;Dtm ðÞ¼ ^ rð^ y0;DtmÞ ^ vð^ y0;DtmÞ  ¼fð^ r0;^ v0;DtmÞ^ r0þgð^ r0;^ v0;DtmÞ^ v0 _ fð^ r0;^ v0;DtmÞ^ r0þ_ gð^ r0;^ v0;DtmÞ^ v0  : ð4Þ The details for computing fand gare provided in Algorithm 8 of Vallado (2001). For the iterative method, the predicted position ^ rð^ y0;DtÞmust be linearized around ^ rðiÞ 0 and ^ vðiÞ 0for variations D^ rðiÞ 0and D^ vðiÞ 0: ^ rð^ r0;^ v0Þ ^ rð^ rðiÞ 0;^ vðiÞ 0Þþ@^ r @^ r0cðiÞD^ rðiÞ 0þ@^ r @^ v0cðiÞD^ vðiÞ 0; @^ r @^ r0cðiÞ¼fðiÞ@^ r0 @^ r0þ^ rðiÞ 0 @f @^ r0cðiÞþ^ vðiÞ 0 @g @^ r0cðiÞ; @^ r @^ v0cðiÞ¼gðiÞ@^ v0 @^ v0þ^ rðiÞ 0 @f @^ v0cðiÞþ^ vðiÞ 0 @g @^ v0cðiÞ: ð5Þ The dependence over Dthas been dropped for conciseness. The GTDS method simplifies this by setting @f=@r0¼@f=@v0¼@g=@r0¼@g=@v0¼0|, reducing the derivatives to @^ r=@^ r0cðiÞfðiÞI3and @^ r=@^ v0cðiÞgðiÞI3. The matrix I3is the identity matrix of order 3. Applying this technique to the linearization of the predicted position results in: ^ rð^ r0;^ v0Þ ^ rð^ rðiÞ 0;^ vðiÞ 0Þþ@^ r @^ r0cðiÞD^ rðiÞ 0þ@^ r @^ v0cðiÞD^ vðiÞ 0 fðiÞ^ rðiÞ 0þgðiÞ^ vðiÞ 0þfðiÞD^ rðiÞ 0þgðiÞD^ vðiÞ 0 ¼fðiÞð^ rðiÞ 0þD^ rðiÞ 0ÞþgðiÞð^ vðiÞ 0þD^ vðiÞ 0Þ¼fðiÞ^ rðiþ1Þ 0þgðiÞ^ vðiþ1Þ 0: ð6Þ Thus, the predicted measurement at tmcan be expressed as: ^ rm¼^ rð^ r0;^ v0;DtmÞ fðiÞ mI3^ rðiþ1Þ 0þgðiÞ mI3^ vðiþ1Þ 0 ¼ fm00gm00 0fm00gm0 00fm00gm 2 6 43 7 5 ðiÞ ^ rðiþ1Þ 0 ^ vðiþ1Þ 0 "# :ð7Þ Using the form in Eq. (7), the vector of predicted positions (^ z) can be expressed by linearly combining the estimation ^ yðiþ1Þ 0with the Kepler functions at the corresponding instants. The position residuals in z¼A^ yðiþ1Þ 0þnare thus an explicit linear function of the estimation ^ yðiþ1Þ 0. By solving the linear least-squares in Eq. (3) with Wbeing the identity matrix (non-weighted), the solution is given by: ^ yðiþ1Þ 0¼A|AðÞ 1A|z¼Hz; Cy¼HCrH|;whereCr¼diag Cr1;;Crn ðÞ;ð8Þ is that of the GTDS method. Emphasis has to be placed on the choice of not weighting position residuals despite the availability of an uncertainty characterization. The result is a fit that minimizes errors in all directions equally. The covariance of the estimation is computed from the matrix H¼A|AðÞ 1A|as a linear approximation after the estimation is computed (in the last iteration). For a more compact form of the solution to Eq. (8), see Siminski (2016). 3.2. Keplerian fitting of radar observables - KEP This method modifies GTDS by fitting measurements directly with Keplerian dynamics, incorporating rangerate information. The measurement vector includes all observables directly, denoted as z¼q1;Az1;el1;_ q1;...;½ qn;Azn;eln;_ qn|. The function hðy;tÞfor these measurements, as well as the derivatives @h=@y, is found in the A.1. Incorporating range-rate requires predicted position and velocity, computed via _ fand _ gfunctions, as shown in Eq. (4). The computation of @h=@^ y0, see Eq. (2), requires the predicted full state derivatives, @P=@^ y0. The position part is found in Section 3.1, Eq. (5). Predicted velocity derivatives are also required, but the GTDS approximation neglecting Kepler function derivatives is unsuitable here. An analytical approximation for f;g;_ f, and _ gderivatives has been employed using a fourth order Taylor expansion around ^ r0and ^ v0. Detailed expressions are omitted for brevity but can be found in Vallado (2001). Thus, the element @P=@^ y0can be computed at each Dtmwith the necessary precision. With @h=@^ y0computed, all elements to solve the linear least-squares problem in Eq. (3) are available. Errors DzðiÞare computed using the propagated estimates ^ yðiÞ m, J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4971 via Eq. (4), evaluated in the observation function (A.1). This is the second analytical method, denoted simply as KEP. 3.3. J2fitting of radar observables - KEP + J2 A third algorithm incorporating J2perturbation has been developed. While analytical methods exist for oblate planet perturbations, they often sacrifice precision and are valid only for specific eccentricity ranges (Martinusi et al., 2015). This work uses an solution valid for short time spans, based on the generalized equinoctial orbital elements (GEqOE) described in Bau `et al. (2021). 3.3.1. GEqOE formulation Given a perturbation force derivable from a potential energy U, the GEqOE incorporate the perturbing potential in their definition, generalizing the classical equinoctial elements (Broucke and Cefola, 1972). This formulation yields the generalized semi-major axis (a) and Laplace vector (lg), defining a non-osculating ellipse in the orbital plane. The elements p1and p2are projections of galong in-plane axes, while Kepler’s equation is written using the generalized mean longitude L: L¼Kþp1cos Kp2sin K;ð9Þ where Kis the generalized eccentric longitude. The generalized mean motion mderives from the total energy, which includes U. Elements q1and q2(equal to the classical p and q) complete the GEqOE set v¼½m;p1;p2;q1;q2;L. The first-order time derivatives of these elements are explicit functions of Uand other forces contained in P(drag, luni-solar perturbations,...). This general formulation can be simplified to include the J2perturbation only: _ m¼0; _ p1¼p2hc r2I^z  1 c X aþ2p2  U; _ p2¼p1I^zhc r2  þ1 c Y aþ2p1  U; _ q1¼IY r; _ q2¼IX r; _ L¼mþhc r2I^z1 c 1 aþa1r a  U: ð10Þ Here, mis constant since the only perturbation is conservative and therefore the total energy is conserved. Some definitions needed for the use of the GEqOE can be found in the A.4.2.Algorithms 1,2 in the A.7 convert between inertial position/velocity and GEqOE. Numerical integration of Eq. (10) requires the Xand Y values. Given initial GEqOE values, the generalized Kepler’s equation is solved, see Eq. (9), with the output Kused then to compute r;X, and Y: r¼a1p1sin Kp2cos KðÞ; X¼aap1p2sin Kþð1ap2 1Þcos Kp2  ; Y¼aap1p2cos Kþð1ap2 2Þsin Kp1  :ð11Þ In Bau `et al. (2021), GEqOE’s slower evolution is reported to reduce numerical integration error. Applying a Taylor expansion in this state representation is thus an interesting alternative. This section aims to derive an analytical solution for IOD. 3.3.2. Taylor expansion in GEqOE A special notation is introduced to differentiate the ODEs in Eq. (10). The function fs¼1=ssimplifies the application of the chain rule in derivation when an element is dividing. For example, _ r¼l=cp 2sin Lp1cos LðÞ, expressed with the true longitude L, can be rewritten as _ r¼lfcrpl. The time derivative becomes € r¼l_ fcrplþ  fc_ rplÞ, requiring _ fs¼_ s=s2. A function that computes time derivatives of fshas been implemented, taking the vector s;_ s;€ s;½as input (see A.2). The function gs¼s_ sis also introduced, appearing when differentiating s2(see A.3 for time derivatives of gs). First-order derivatives of the Taylor expansion involve evaluating the ODEs. For higher orders, elements in Eq. (10) are split into easily derivable parts, as shown in Eq. (A.9) for ^z. This method compacts the equations of motion, see Eq. A.6. Full derivation up to fourth order and symbol definitions are given in A.4. After computing the derivatives at the epoch, a Taylor expansion propagates the state vector vv 0;Dtm ðÞ¼ m;p1;p2;q1;q2;L½ |to all Dtm, which is then converted to Cartesian using Algorithm 2. The analytical derivative of vv 0;Dtm ðÞwith respect to v0is computed to obtain @v @v0. This involves differentiating the functions fsand gs(and their time derivatives) with respect to a generic orbital element (see A.5). The resulting analytical derivatives have been validated using Taylor differential algebra with Hipparchus. To convert between GEqOE and Cartesian states (see Algorithms 1 and 2) and derivatives, the Jacobians of both transformations, RV2GEQOE @v @X and GEQOE2RV @X @v  , are used (Bau `et al., 2021). Cartesian derivatives are obtained by multiplying the Jacobians evaluated at the correct instant: @X @X0¼@X @vt @v @v0  @v @Xt0 ;ð12Þ validated by numerically integrating the variational equations with the J2perturbation added. This completes the method, referred to as KEP + J2. 3.4. Adding extra information: Orbital Plane-based Orbit Determination - OPOD Precise initial orbit determination using a single radar track is challenging due to the scarcity of data. This section presents an alternative approach, Orbital Plane-based Orbit Determination (OPOD), to improve the quality and J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4972 reliability of the estimation by exploiting information from a known past state. The key assumption is that for LEO objects maneuvers are usually performed to maintain altitude and do not significantly affect the orbital plane, as changing the inclination (i) or right ascension of the ascending node (RAAN, X) is expensive. In OPOD, the orbital plane is predicted from the known past orbit with relatively high precision, assuming that process noise (primarily drag uncertainty) does not directly affect the iand Xprediction, even in the presence of maneuvers. The predicted orbital plane information is then added as virtual measurements at the radar epoch (t0), independent of the measurements, to constrain the region of inertial space that the spacecraft could occupy at that instant. This reduces the uncertainty of the estimation, with the level of improvement depending on the uncertainty assigned to these virtual measurements. It is important to note that OPOD does not decouple plane attitude estimation from the in-plane estimation, nor does it assume the plane to be deterministic or even constant. The state is still estimated using the same leastsquares fitting approach, considering all measurements (radar track and virtual orbital plane) through the dynamical model. The addition of orbital plane information at the radar epoch simply provides extra constraints to improve the estimation accuracy. In fact, when the fitting dynamics include the J 2 perturbation, even if the orbital plane is assumed partially known at a given instant it does change along the track. The measurement function hðy;tÞfor these new observables, as well as the analytical derivative with respect to the Cartesian state @h=@yare in the A.6. The addition of these virtual measurements can be applied to either the KEP or KEP + J2 fitting of radar measurements, but not the GTDS method as it only works with position information. 4. Simulated results For this section, simulations of the satellite’s trajectory and radar measurements are used to evaluate IOD methods’performance. Orekit (Maisonobe et al., 2010) handles all high-fidelity satellite dynamics and radar simulations, see Section 4.1. Section 4.2 presents preliminary testing of radar-only IOD methods, establishing baseline performance and track length dependency. Section 4.3 tests the orbital-plane enhanced method (OPOD), showing the impact of adding orbital plane information and uncertainty even in presence of maneuvers. Section 4.4 details the simulated testing methodology and scenarios, generating track data for Section 4.5, where a k2metric tests uncertainty characterization of estimation errors. 4.1. Dynamical modeling and radar characterization For IOD methods relying on sensor data, the testing involves simulating ‘‘real”satellite trajectories and generating radar tracks. The OPOD algorithm, reliant upon predicted orbital plane, requires two distinct dynamics. Table 1 lists the dynamical models for numerical simulations. The main difference between True and Prediction dynamics is the atmosphere model and harmonics degree. NRLMSISE-00 (Picone et al., 2002) and DTM-2000 (Bruinsma et al., 2003) models account for tabulated solar weather data, differing in order to reflect atmospheric uncertainty in practice. High-fidelity (HF) and low-fidelity (LF) satellite models are defined in Tables 2 and 3, respectively. HF satellites are defined by its dimensions, solar array area (As) and orientation in local axis (Ax), drag (CD), absorption (ac) and reflection (rc) coefficients. LF satellites use a simple spherical drag model. HF models pair with True dynamics for simulating trajectories and radar measurements, while LF models use Prediction dynamics for orbital plane predictions. Satellites named for real-world equivalents vary in orbit type: Sentinel-1A (700 km), Swarm-C (430 km), Starlink-1 and 2 (540 km). Table 4 provides initial conditions for simulated scenarios, generating diverse radar track lengths. Only Sentinel-1A and Swarm-C have corresponding LF models for OPOD testing. Other IOD methods (Sections 3.1, 3.2, 3.3) are tested on all satellites. Simulations use Orekit (Maisonobe et al., 2010) for dynamics and measurement generation. Orekit handles sensor measurements, taking into account reflection time and ionospheric delay. Measurement generation adds noise with constant covariance and no biases, a simple model suitable for limited orbital regimes. The radar’s uncertainty characterization is assumed known for estimation. Table 5 Table 2 Parameters defining the realistic satellites Sentinel-1A, Swarm-C and Starlink-1/2 using Orekit’s BoxAndSolarArraySpacecraft force model class. HF Satellite m(kg)xl(m)yl(m)zl(m)As(m2)Ax[-] CD[-] ac[-] rc[-] Sentinel-1A 2270 1.02 1.34 3.2 25.46 y1.5 0.8 0.5 Swarm-C 440 1.3 8 0.9 0.1 y0.95 0.9 0.5 Starlink-1 1250 2.8 2.3 0.3 25.2 z0.9 0.9 0.5 Starlink-2 800 2.5 2 0.2 20 z0.92 0.9 0.5 Table 1 The two dynamical models defined, one for the true trajectory and the other for the prediction. SRP stands for Solar Radiation pressure. Model Earth harmonics Atmosphere Sun Moon SRP True ½150;150DTM-2000 yes yes yes Prediction ½40;40NRLMSISE-00 yes yes yes J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4973 specifies radar parameters for three different stations used separately in the simulations, enriching track length variability and data density across scenarios. Radar location is given by its geodetic coordinates longitude (k), latitude (/) and height (h). The radar antenna orientation is given by 3 consecutive rotations ðux;uy;uzÞaround its local x;y, and zaxes. Initially, it is parallel to the XY plane of the topocentric frame. The field of regard (FoR) is defined by scanning angles U and V, limited in both positive and negative directions, forming a volume outside two intersecting concentric cones. Radar 1 is a simulation of S3TSR (Gomez et al., 2019), and some of its characteristics have been omitted for confidentiality. Radars 2 and 3 have been defined for testing purposes, with similar characteristics to the S3TSR simulation. In all cases radar characterization is assumed perfectly known. 4.2. Preliminary testing of the IOD algorithms With the simulation and fitting methods in place (see Section 3), a controlled test is needed to compare them. Understanding each method’s limitations helps to interpret the estimation statistics (see Section 4.5.1). The approach involves applying the three IOD algorithms to a single radar track using only radar data. To study the effect of track length on estimation accuracy and uncertainty, a long track is shortened incrementally by removing measurements from the ends and repeating the fitting. The longest track from Radar 3, simulating Sentinel-1A over 4 days, is chosen. It consists of 72 measurements spanning 284 s. 4.2.1. Fitting of noiseless measurements First, it is essential to consider the true values of the measurements, rather than the noisy ones, to emphasize the base limitations of each method. Propagating the true position and velocity at the epoch of the test track with Kepler dynamics results in position errors on the order of 100 meters and range errors over 50 m (Fig. 2, Left). This trajectory coincides with the true one at the epoch, yet the Kepler model incurs range errors an order of magnitude greater than the radar measurements’uncertainty. Applying the GTDS algorithm to the true position measurements results in the estimation with errors shown in Fig. 2 (Center). The square sum of errors (n|n) is minimized, and without weights, the position errors are lower on average along the complete trajectory, but not necessarily the range errors. In Fig. 2 (Right), the position fitting is done by minimizing n|Wn, with W ¼C1 r. The position covariance computed with an UT has greater variance in the range direction than the original measurement space uncertainty, but it is still a thin disk in Cartesian space. The estimation error is on the order of a kilometer, while still fitting exact position values. The increased weight of the error in the range direction forces lower residuals in that direction, while residuals in perpendicular directions are less relevant, especially for measurements further from the radar. The problem lies in pairing Kepler dynamics with a strongly directional weighting, as the model cannot maintain low range residuals for such a long pass (Fig. 2, Left), so the fitting algorithm has to change the trajectory considerably. When fitting true radar observables directly with a Kepler model (see Section 3.2), the problem exacerbates (Fig. 3, Left), with an epoch position error of 1.5 km. This is due to the range being even more relevant to the weighting and the addition of a range-rate constraint with very low uncertainty. The approximated J2propagator (Section 3.3.2) is tested up to the fourth term of the Taylor expansion (see Table 6 for constants). Absolute position errors compared to numerical propagation are below 1 cm for a 100-s propagation (Fig. 4). Fig. 3 (Center) shows the errors when the true state is propagated from the epoch using the J2propagator, with range residuals barely over 1 meter at worst. Consequently, the radar observables fitting is much more similar to the original trajectory (Fig. 3, Right), indicating this might be superior to GTDS if only position errors are considered. Fig. 5 compares the methods for different track lengths with noiseless measurements, representing the error at the epoch when the fitted trajectory gets shorter as measurements are extracted in pairs from the sides of the test track. A shorter trajectory can be better fitted with simplified dynamics. The KEP fit of measurements is the worst of Table 4 Initial condition (I.C.) for the reference orbit defined for the scenarios in Table 8, in Cartesian coordinates (Earth EME2000 inertial frame). The initial epoch t0is at 00:00:00 UT. I.C. x(km)y(km)z(km)vx(km=s)vy(km=s)vz(km=s)t0 S1A-1 1459.975 436.989 6916.264 3.8952 6.282 1.219 01/05/2022 SWC-1 254.597 241.108 6792.396 7.151 2.655 0.1728 06/07/2020 SRL-1 4970.721 3941.433 602.664 4.458 4.973 4.244 27/08/2023 SRL-2 5030.115 3838.607 754.660 4.040 4.249 5.315 27/08/2023 Table 3 Parameters defining the low fidelity satellites Sentinel-1A 1 and Swarm-C 1 using Orekit’s isotropic drag and radiation force model classes. These parameters are: mass, drag cross section, drag coefficient, solar radiation cross section and reflection coefficient. LF Satellite m(kg)S(m2)CD[-] SRPs(m2)rc[-] Sentinel-1A 1 2200 9.5 2 10 1.4 Swarm-C 1 420 4 2.56 4 1.4 J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4974 the three, but very short tracks do not degrade as much in terms of estimated position bias. The KEP + J2 fit with J2 perturbation seems to be the best in this regard, but GTDS has a small advantage in velocity errors, which can only occur under the condition of perfect position measurements. Since the GTDS method assigns equal weight to all residuals, errors stemming from the basic dynamics are effectively averaged out. 4.2.2. Fitting of noisy measurements When noise is introduced (Fig. 6) and results are averaged over multiple fittings, there is little difference between the methods for tracks shorter than 100 s in terms of average position error. GTDS has a slight disadvantage in Fig. 3. Position and range errors using true radar measurements, including range-rate, for the same sample track as in Fig. 2. Left: Kepler dynamics fitting. Center: Analytical J2propagation from the true state at t0. Right: Analytical fitting with J2. Fig. 2. Kepler-based trajectories vs. true positions for a 284-s radar track. Range error is the position error projected along the satellite-radar vector. Left: Trajectory starting from the true state at epoch. Center: Fitting of true positions without weights. Right: Fitting of true positions with weights. Table 5 Complete characterization of the radar stations used for the simulated testing. Name Radar 1 Radar 2 Radar 3 ðk;/Þ()ð5:59;37:17Þð22:3;78:49Þð175:68;37:86Þ h(m) 142.32 22.14 36.2 ðux;uy;uzÞ()- ð40;0;0Þð30;0;0Þ ðUmin;UmaxÞ()- ð30;30Þð20;20Þ ðVmin;VmaxÞ()- ð20;20Þð45;45Þ rq(m) 7 6:56:5 r_ q(m/s) 0:40:35 0:35 rAz ()0:30:25 0:25 rel ()0:20:15 0:15 nAz;el 0:043 0:043 0:043 rt(s) 7 4 4 J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4975 velocity errors for the first half of the interval, due to the lack of range-rate information. As the track lengthens, the Keplerian fitting of measurements shows its fundamental flaw, with added information deteriorating the estimation due to the pairing of non–homogeneous data (diskshaped covariance), with an unsuitable dynamical model. GTDS suffers less from the Kepler approximation thanks to the lack of weighting, enabling error averaging. The approximated J2method performs best overall, particularly for the longest radar tracks. Fig. 7 shows the uncertainty evolution for different track lengths. GTDS is more conservative in uncertainty estimation, partly due to the absence of range-rate data and the indirect use of measurement uncertainty. The information in Fig. 7 does not consider correlations, providing an incomplete picture of the covariance. The evolution is compatible with the errors, suggesting correct behavior of this estimation method. The Keplerian fit of radar observables suffers from overconfidence for tracks longer than 100 s, Fig. 4. GEqOE Taylor expansion relative error improves by an order of magnitude with each added term for the first 200 s. The normalizing value (for relative error computation) is the numerical propagation output with 1014 tolerance (4a). Cartesian error is computed with the fourth-order approximation (4b). Fig. 5. Fitting errors (at t0) for decreasing radar track durations using the three main IOD algorithms from Section 3. The same radar track as in Figs. 2 and 3 is used. All fittings use true measurement values. Table 6 Constants used for the analytical J2propagator developed in Section 3.3.2. l(km3=s2)Re(km) J 2 398600:4418 6378:137 1:082626683553 103 J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4976 4.5.2. k2statistics on OPOD The SEN-1A-2 scenario has been used to test the estimation algorithms’performance, now including predicted orbital plane information. For each sampled radar track, the propagated initial condition information (randomly chosen from the I.C. sampling) was used to include predicted inclination and RAAN in the fitting. This was only done with the radar observables fitting methods, as GTDS is unweighted and does not allow for non-position measurements. Fig. 18 shows that including orbital plane information has little effect on KEP fittings, only solving the issue with very short tracks but not addressing the degradation for longer ones. However, KEP+J2 sees its only drawback fixed and now behaves more appropriately for all track lengths. Although the obtained k2values are lower than expected (Fig. 20), a conservative covariance is generally more desirable to avoid false correlations. This is due to the virtual measurements being more consistent than indicated by the used uncertainty value (ri;X¼5e3in all Fig. 17. Estimation only k2percentage of anomalous estimations for Radar 2 (rt¼4 seconds). Fig. 18. Estimation only k2mean values for Radar 2 (rt¼4 seconds), including the information of the orbital plane with ri;X¼5e3in the KEP and KEP+J2methods. Fig. 19. Estimation only k2variance values for Radar 2 (rt¼4 seconds), including the information of the orbital plane with ri;X¼5e3in the KEP and KEP+J2methods. J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4983 cases), suggesting there is room for reducing the estimation uncertainty. GTDS results are unchanged as no extra information is added to the fit. (see Fig. 19) Another considerably improved aspect not present in the k2statistics is the improved accuracy of the estimations compared to the GTDS method. This is not surprising, as the augmented methodology uses range-rate and orbital plane data to generate a better estimation in all aspects, including uncertainty realism. 5. Conclusions This work introduces a novel approach to initial orbit determination using data from a single radar track combined with the predicted orbital plane of the object. Focusing on the need for rapid algorithms in operational scenarios, a least-squares fitting procedure incorporating an analytically formulated approximation of the dynamics under the J2perturbation is developed. The algorithm includes range-rate observables, distinguishing it from similar methods. A comprehensive evaluation is conducted through simulated tests, comparing its performance against a classical range and angles fitting method (GTDS), exploring the impact of track length and measurement density on full state estimation. Proper full state estimation using radar observables directly does not work as intended with Keplerian dynamics due to the magnitude of propagation errors being greater than precise range measurements modeled after modern radar technology. Improving the fidelity of the fitting dynamics by including the J2perturbation solves this issue while maintaining the fully analytical aspect of the methodology. The resulting fast IOD method is applicable to maneuver detection and finding correlations between cataloged objects and tracks. The proposed method behaves as expected for tracks longer than 40 s, especially with high data density. However, there is a problem for shorter arcs, where the nonlinear relation between the estimate and radar observables affects the velocity part of the estimation. Covariance inflation post-processing could be applied to take advantage of the higher quality estimation. The developed method admits the inclusion of chosen parts of the predicted state as virtual measurements. Given the orbital plane’s stability in LEO and relatively short integration windows, the precision of predicted iand Xjustifies their inclusion in the estimation process. The resulting method, denoted as OPOD, shows improved accuracy and reliability for all track lengths. The analytical time derivatives of the GEqOE equations of motion developed in this work could be employed in an enhanced numerical propagator that makes use of higher order derivatives, potentially improving accuracy and reducing the total number of steps needed for a given propagation length and tolerance. This line of research is currently being pursued as an extended application of the work presented here. Another promising avenue for future research is the direct fitting of GEqOEs. This approach would allow for the incorporation of range-rate information into the estimation process and could potentially improve uncertainty representation for the linear approximation of the covariance in the context of very short tracks. By leveraging the inherent structure of GEqOEs, where two elements are functions of inclination and longitude of the node, this method could naturally integrate orbital plane orientation information, as demonstrated in our OPOD method. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments Jose M. Montilla and R. Vazquez gratefully acknowledge support by grant TED2021-132099B-C33 funded by MICIU/AEI/10.13039/501100011033 and by ‘‘European Union NextGenerationEU/PRTR”. Special thanks are due to Guillermo Escribano for his help in developing the ideas and methods in this work. During the preparation of this work the author(s) used ChatGPT in order to revise and correct grammar and Fig. 20. Estimation only k2percentage of anomalous estimations for Radar 2 (rt¼4 seconds), including the information of the orbital plane with ri;X¼5e3in the KEP and KEP+J2methods. J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4984 refine the language. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication. Appendix A. Auxiliary functions A.1. Radar measurement function derivatives The measurements provided by the radar in LEO are range (q), azimuth (Az), elevation (el) and range-rate (_ q). The observation function q;Az;el;_ q½¼hðy;tÞof these observables is given by: q¼ðrPRÞ|ðrPRÞ½ 1=2¼q|q½ 1=2; Az ¼atan2 e| 1qcR;e| 2qcR ðÞ¼atan2 e| 1TR Iq;e| 2TR Iq  ; el ¼arcsin e| 3qcR q  ¼arcsin e| 3TR Iq q  ; _ q¼1 qðrPRÞ|ðvVRÞ¼q|u q; ðA:1Þ where q¼rPRand u¼vVR. The unit vectors e1;e2;e3 ðÞform a Cartesian frame base. The fitting algorithms in Sections 3.2 and 3.3 require the derivatives of these measurements with respect to position (r) and velocity (v): @h @y¼ qrqv AzrAzv elrelv _ qr_ qv: 2 6 6 6 4 3 7 7 7 5 :ðA:2Þ The explicit functions of these derivatives are: qr¼q| q;qv¼0|; Azr¼1 1þtan2ðAzÞe| 1tanðAzÞe| 2 ½ TR I e| 2qcR;Azv¼0|; elr¼1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1sin2ðelÞ pe| 3TR IsinðelÞqr  1 q;elv¼0|; _ qr¼ðu_ qqrÞ1 q;_ qv¼q q; ðA:3Þ which requires the knowledge of radar’s inertial position (PR), velocity (VR) and orientation (TR I). A.2. Time derivatives of the inverse function In Section 3.3, the linear least-squares algorithm is applied to a dynamical model that includes the J 2 perturbation using the GEqOE from Bau `et al. (2021). To maintain the analytical aspect of the fitting method, a Taylor expansion has been developed. The coefficients of this expansion are detailed in A.4, which extensively uses the inverse function fs¼1 sand its time derivatives: _ fs¼ _ s s2; € fs¼2_ s2 s3€ s s2; f ... s¼6€ s_ s s3s ... s26_ s3 s4; f ð4Þ s¼8s ð3Þ_s s3þ6€s2 s336 _s2€s s4s ð4Þ s224 _s4 s5: ðA:4Þ A.3. Time derivatives of the function of s times _ s As with the inverse function in A.2, the function that multiplies sby its time derivative _ s;gs¼s_ s, is repeatedly used in the derivation of the J 2 Taylor propagator (see A.4). The time derivatives of gsare included here: _gs¼_ s2þs€ s; €gs¼3_ s€ sþss ... g ... s¼3€ s2þ4_ ss ð3Þþss ð4Þ; g ð4Þ s¼10€ ss ð3Þþ5_ ss ð4Þþss ð5Þ: ðA:5Þ A.4. Time derivatives of the J2GEqOE equations of motion To develop an efficient J 2 perturbation propagator, this work uses a Taylor expansion of the solution in terms of GEqOE (Bau `et al., 2021). First, the equations of motion are simplified into a compact form (Eq. A.6), utilizing the inverse function (see A.2). This section includes the time derivatives of Eq. A.6 up to fourth order. All simplifications are defined in A.4.2, facilitating a systematic application of the chain rule. From A.4.1–A.4.4, all steps for computing the Taylor J 2 propagator coefficients are included. A.4.1. First order derivatives The first order time derivatives of the GEqOE can be computed with: _ m¼0;_p1¼p2dwh ðÞfcn1U;_p2¼p1whdðÞþfcn2U; _ q1¼IsL;_ q2¼IcL;_ L¼mþdwhfcCU; ðA:6Þ where sL¼sin L;cL¼cos L. The rest of definitions for the elements in Eq. (A.6) can be found next. Higher order derivatives will require the computation of: cos L¼X r; sin L¼Y r; _ X¼_ rcos Lh rsin L; _ Y¼_ rsin Lþh rcos L: ðA:7Þ A.4.2. Second order derivatives The second order time derivatives of the GEqOE can be computed with: J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4985 € m¼0;€q1¼_ IsLI_sL; €p1¼_p2dwh ðÞþp2_ d_wh  _ fcn1þfc _ n1  Ufcn1_ U;€q2¼_ IcLI_cL; € p2¼_ p1whdðÞþp1_ wh_ d  þ_ fcn2þfc _ n2  Uþfcn1_ U;€ L¼ _ d_ wh_ fcCþfc _ C  UfcC_ U: ðA:8Þ All the elements that appear in Eq. (A.8) are computed with the following formulation. The first time derivative of ^zis computed from: ^z¼2Yq2Xq1 ðÞ rð1þq2 1þq2 2Þ¼C D¼CfD;_ ^z¼_ CfDþC_ fD; C¼2Yq2Xq1 ðÞ;_ C¼2_ Yq 2þY_ q2_ Xq 1X_ q1  ; D¼rð1þq2 1þq2 2Þ¼rqs;_ D¼ _ rqsþr_ qs; qs¼ð1þq2 1þq2 2Þ;_ qs¼2q1_ q1þq2_ q2 ðÞ; ðA:9Þ where _ fDis calculated from _ D. The first time derivative of Uand Uzare: U¼ A r3ð13^z2Þ¼AUzfr3;_ U¼A_ Uzfr3þUz_ fr3  ; Uz¼ð13^z2Þ;_ Uz¼6^z_ ^z: ðA:10Þ The value _ fr3in Eq. (A.10) is computed from the derivative of r3: dðr3Þ dt ¼_ rc¼3r2_ r:ðA:11Þ Given that the generalized semi-major axis is a¼ðl=m2Þ1=3, the derivative of Cand b, which in turn allows to compute _ f1þb;_ fband _ fa, are given by: C¼faþa1r=aðÞ;_ C¼_ faþ_ a1r a  a_r a; a¼1 1þb¼f1þb;_ a¼_ f1þb; b¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1p12p2 2 q;_ b¼_ ps 2b¼1 2_ psfb; ps¼1þp2 1þp2 2;_ ps¼2p1_ p1þp2_ p2 ðÞ: ðA:12Þ The first time derivative of Iis given in: I¼3A hr3^zð1q2 1q2 2Þ¼3A^zd h¼3A^zdfh;_ I¼3A_ ^zdþ^z_ d  fhþ^zd_ fh hi ; d¼1q2 1q2 2;_ d¼_ qs; h¼hr3¼hrc;_ h¼_ hrcþh_ rc; h¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c22r2U p;_ h¼gc2grUr2_ U  fh; c¼l2 m  1=3ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1p2 1p2 2 p¼l2 m  1=3 b;_ c¼l2 m  1=3_ b: ðA:13Þ This includes the computation of _ fhfrom _ h,and_ fhis computed from _ hfor later use. The same goes for _ fcand gc, calculated from _ c. The first time derivative of dis: d¼hc r2¼hcðÞfr2; _ d¼_ h_ c  fr2þhcðÞ _ fr2; dðr2Þ dt ¼_ rs¼2r_ r¼2gr; ðA:14Þ where _ fr2is computed from the derivation of r2. Finally, the time derivatives of wh;n1;n2;sLand cLare: wh¼I^z;_ wh¼_ I^zþI_ ^z; n1¼X aþ2p2;_ n1¼_ X aþ2_ p2; n2¼Y aþ2p1;_ n2¼_ Y aþ2_ p1; sL¼Y r¼Yfr;_sL¼_ YfrþY_ fr; cL¼X r¼Xfr;_cL¼_ XfrþX_ fr: ðA:15Þ In this section the computation of every fstime derivative has been explicitly written. In the following, these comments are omitted for the sake of brevity, and it will be assumed all instances can be computed with the corresponding derivatives of s(to the needed order). A.4.3. Third order derivatives The third order time derivatives of the GEqOE can be computed with: m ... ¼0; q ... 1¼ € IsL2_ I_sLI€sL; q ... 2¼ € IcL2_ I_cLI€cL; p ... 1¼€ p2dwh ðÞþ2_ p2_ d_ wh  þp2€ d€ wh  € fcn1þ2_ fc_ n1þfc € n1  U2_ fcn1þfc _ n1  _ Ufcn1€ U; p ... 2¼€ p1whdðÞþ2_ p1_ wh_ d  þp1€ wh€ d  þ€ fcn2þ2_ fc_ n2þfc € n2  Uþ2_ fcn2þfc _ n2  _ Uþfcn2€ U; L ... ¼€ d€ wh€ fcCþ2_ fc_ Cþfc € C  U2_ fcCþfc _ C  _ UfcC€ U: ðA:16Þ All the elements that appear in Eq. (A.16) are computed with the following formulation. The second time derivative of ^zcomes from: € ^z¼€ CfDþ2_ C_ fDþC€ fD;€ C¼2€ Yq 2þ2_ Y_ q2þY_ q2€ Xq12_ X_ q1X€ q1  ; € D¼ € rqsþ2_ r_ qsþr€ qs!€ fD;€ qs¼2_ q2 1þq1€ q1þ_ q2 2þq2€ q2  :ðA:17Þ The second order time derivatives of Xand Yare computed with: _ X¼_rcos Lh rsin L¼_rcLwsL;€ X¼€rcLþ_r_cL_wsLw_sL; _ Y¼_rsin Lþh rcos L¼_rsLþwcL;€ Y¼€rsLþ_r_sLþ_wcLþw_cL; w¼h r¼hfr;_w¼_ hfrþh_ fr; _ r¼l cp2sin Lp1cos LðÞ¼lfcrpl;€ r¼l_ fcrpl þfc_ rpl  ; rpl ¼p2sLp1cL;_ rpl ¼_ p2sLþp2_sL_ p1cLp1_cL: ðA:18Þ The second time derivative of Uand Uzcan be computed with: € U¼A€ Uzfr3þ2_ Uz_ fr3þUz€ fr3  ;€ Uz¼6_ ^z2þ^z€ ^z  ;d2ðr3Þ dt2¼€rc¼32r_r2þr2€rðÞ: ðA:19Þ The second order derivatives of Cand bare given by: € C¼€ faþ€ a1r a  2_ a_ r aa€ r a;€ a¼€ f1þb; € b¼1 2€psfbþ_ps_ fb  !€ f1þb;€ fb;€ps¼2_p2 1þp1€p1þ_p2 2þp2€p2  : ðA:20Þ J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4986 The second time derivative of Iis completely defined in: € I¼3A€ ^zdþ2_ ^z_ dþ^z€ d  fhþ2_ ^zdþ^z_ d  _ fhþ^zd€ fh hi ; € h¼€ hrcþ2_ h_ rcþh€ rc; € h¼_gc2_grU4gr _ Ur2€ U  fhþgc2grUr2_ U  _ fh; € c¼l2 m  1=3€ b; € d¼ € qs: ðA:21Þ The second time derivative of dresults from: € d¼€ h€ c  fr2þ2_ h_ c  _ fr2þhcðÞ € fr2; d2ðr2Þ dt2¼€ rs¼2_gr:ðA:22Þ The second order time derivatives of wh;n1;n2;sLand cL are: € wh¼€ I^zþ2_ I_ ^zþI€ ^z; € n1¼€ X aþ2€ p2; € n2¼€ Y aþ2€ p1; €sL¼€ Yfrþ2_ Y_ frþY€ fr; €cL¼€ Xfrþ2_ X_ frþX€ fr: ðA:23Þ A.4.4. Fourth order derivatives The fourth order time derivatives of the GEqOE can be computed with: m ð4Þ¼0; q ð4Þ 1¼I ... sL3€ I_sL3_ I€sLIs ... L; q ð4Þ 2¼I ... cL3€ I_cL3_ I€cLIc ... L; p ð4Þ 1¼p ... 2 dwh ðÞþ3€p2_ d_wh  þ3_p2€ d€wh  þp2d ... w ... h  f ... c n1þ3€ fc_ n1þ3_ fc€ n1þfcn ... 1  U3€ fcn1þ2_ fc_ n1þfc € n1  _ U3_ fcn1þfc _ n1  € Ufcn1U ... ; p ð4Þ 2¼p ... 1 whdðÞþ3€ p1_ wh_ d  þ3_ p1€ wh€ d  þp1w ... hd ...  þf ... c n2þ3€ fc_ n2þ3_ fc€ n2þfcn ... 2  Uþ3€ fcn2þ2_ fc_ n2þfc € n2  _ Uþ3_ fcn2þfc _ n2  € Uþfcn2U ... ; L ð4Þ¼d ... w ... hf ... c Cþ3€ fc_ Cþ3_ fc€ CþfcC ...  U3€ fcCþ2_ fc_ Cþfc € C  _ U3_ fcCþfc _ C  € UfcCU ... : ðA:24Þ All the elements that appear in Eq. (A.24) are computed with the following formulation. The third time derivative of ^zis: ^z ... ¼C ... fDþ3€ C_ fDþ3_ C€ fDþCf ... D ;C ... ¼2Y ... q2þ3€ Y_q2þ3_ Y€q2þYq ... 2X ... q13€ X_q13_ X€q1Xq ... 1  ; D ... ¼r ... qsþ3€ r_ qsþ3_ r€ qsþrq ... s ;q ... s¼23 _ q1€ q1þq1q ... 1þ3_ q2€ q2þq2q ... 2  : ðA:25Þ The third order time derivatives of Xand Yare computed in: X ... ¼r ... cLþ2€ r_cLþ_ r€cL€wsL2_w_sLw€sL; Y ... ¼r ... sLþ2€ r_sLþ_ r€sLþ€wcLþ2_w_cLþw€cL; €w¼€ hfrþ2_ h_ frþh€ fr; r ... ¼l€ fcrpl þ2_ fc_ rpl þfc€ rpl  ; € rpl ¼€ p2sLþ2_ p2_sLþp2€sL€ p1cL2_ p1_cLp1€cL: ðA:26Þ The third order time derivative of Uand Uzcan be computed with: U ... ¼AU ... zfr3þ3€ Uz_ fr3þ3_ Uz€ fr3þUzf ... r3 ! ; U ... z¼63 _ ^z€ ^zþ^z^z ...  ; d3ðr3Þ dt3¼r ... c¼32 _ r3þ6r_ r€ rþr2r ... ðÞ: ðA:27Þ The third order derivatives of Cand bare: C ... ¼f ... aþ€ a1r a  3€ a_ r a3_ a€ r aar ... a;a ... ¼f ... 1þb ; b ... ¼1 2p ... s fbþ2€ ps_ fbþ_ ps€ fb  ;f ... b ;p ... s¼23 _ p1€ p1þp1p ... 1þ3_ p2€ p2þp2p ... 2  :ðA:28Þ The third order time derivative of Iis computed from the next expressions: I ... ¼3A^z ... dþ3€ ^z_ dþ3_ ^z€ dþ^zd ...  fhþ3€ ^zdþ2_ ^z_ dþ^z€ d  _ fhþ3_ ^zdþ^z_ d  € fhþ^zdf ... h  ; h ... ¼h ... rcþ3€ h_rcþ3_ h€rcþhr ... c; h ... ¼€gc2€grU6_gr_ U6gr € Ur2U ...  fhþ2_gc2_grU4gr _ Ur2€ U  _ fhþgc2grUr2_ U  € fh; c ... ¼l2 m  1=3 b ... ; d ... ¼q ... s : ðA:29Þ The third time derivative of dis: d ... ¼h ... c ...  fr2þ3€ h€ c  _ fr2þ3_ h_ c  € fr2þhcðÞf ... r2 ;d3ðr2Þ dt3¼r ... s¼2€gr:ðA:30Þ The third order time derivatives of wh;n1;n2;sLand cLare given by: w ... h¼I ... ^zþ3€ I_ ^zþ3_ I€ ^zþI^z ... ; n ... 1¼X ... aþ2p ... 2 ; n ... 2¼Y ... aþ2p ... 1 ; s ... L¼Y ... frþ3€ Y_ frþ3_ Y€ frþYf ... r ; c ... L¼X ... frþ3€ X_ frþ3_ X€ frþXf ... r : ðA:31Þ A.5. Derivatives of the fsand gstime derivatives with respect to the element x When f ðkÞ s is derived with respect to x, then all the derivatives up to s ðkÞwith respect to x(with standard notation s ðkÞ x) are needed (or up to s ðkþ1Þin the case of g ðkÞ s ). The notation for the cross derivative used here meets @fs ðkÞ=@x¼fs ðkÞ;x , chose on purpose for the compact expressions developed in this appendix section. J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4987 Thus, the function that computes f ðkÞ;x s needs not only the vector s;_ s;;s ðkÞ  , but also sx;_ sx;;s ðkÞ x  as input. Here are the expressions of the function that implements these derivatives: @fs @x¼f x s¼ sx s2; @_ fs @x¼f ð1Þ;x s¼ _sx s2þ2_ssx s3; @€ fs @x¼f ð2Þ;x s¼ €sx s2þ2€ssxþ2_s_sx s36_s2sx s4; @f ... s @x¼f ð3Þ;x s¼ s ... x s2þ2s ... sxþ3€ sx_ sþ€ s_ sx ðÞ s318€s_ssxþ_s2_sx s4þ24 _s3sx s5; @f ð4Þ s @x¼f ð4Þ;x s¼ s ð4Þ x s2þ2 s ð4Þsxþ6€ s€ sxþ4s ... x _ sþs ð3Þ_ sx  s364s ð3Þ_ssxþ62_s_sx€sþ_s2€sx ðÞ þ3€s2sx s4 þ12 12_s2€ssx8_s3_sx s5þ120 _s4sx s6: ðA:32Þ And similarly for g ðkÞ;x s the derivatives are: @gs @x¼g x s¼sx_ sþs_ sx; @_gs @x¼g ð1Þ;x s¼2_ s_ sxþsx€ sþs€ sx; @€gs @x¼g ð2Þ;x s¼3_ sx€ sþ3_ s€ sxþsxs ... þss ... x; @g ... s @x¼g ð3Þ;x s¼6€ s€ sxþ4_ sxs ð3Þþ4_ ss ð3Þ xþsxs ð4Þþss ð4Þ x; @g ð4Þ s @x¼g ð4Þ;x s¼10€ sxs ð3Þþ10€ ss ð3Þ xþ5_ sxs ð4Þþ5_ ss ð4Þ xþsxs ð5Þþss ð5Þ x: ðA:33Þ A.6. Derivatives of inclination and right ascension of the ascending node Consider the specific angular momentum as the cross product of position and velocity both given in an inertial frame (I), h¼rv. The inclination ican be computed as a function of rand vfrom: h¼h|hðÞ 1=2;hz¼k|h;i¼arccos hz h  :ðA:34Þ The right ascension of the ascending node Xis computed from: n¼kh;n¼n|n ðÞ 1=2; n¼n n;X¼arccos i|nðÞ:ðA:35Þ For the derivation of iand Xwith respect to position and velocity let us first define the cross product of a 3 n matrix M by a 3 1 vector das the cross product of each column in M by d. This means that the i-eth column of D¼Mdmeets Di¼Mid. In other words, the cross product of matrix and vector acts column-wise for the next derivation. The differentiation with respect to the state for the inclination observable results in the following: @h @r¼@r @rvþr@v @r¼Iv; @h @v¼@r @vvþr@v @v¼Ir; @h @y¼@h @r;@h @v  ; @hz @y¼e| 3@h @y; @h @y¼h| h @h @y; @i @y¼ 1 ffiffiffiffiffiffiffiffiffiffiffiffi 1cos2i p@hz @ycos i@h @y hi 1 h: For the right ascension of the ascending node on the other hand: @n @y¼e3@h @y¼@h @ye3; @n @y¼n| n@n @y; @n @y¼@n @yn@n @y hi 1 n; @X @y¼ 1 ffiffiffiffiffiffiffiffiffiffiffiffiffi 1þcos2X pe| 1@n @y: A.7. Generalized equinoctial orbital elements auxiliary functions Algorithms 1 and 2 are basic conversion functions needed to work with the generalized equinoctial orbital elements. Algorithm 1. v¼mp1p2q1q2L½ |¼RV2GEQOEfunðrcECI ; vcECI ;l;J2;RÞ 1: r¼krk,v¼kvk 2: h¼rv;h¼khk;_ r¼rv r 3: ^z¼rð3Þ=r 4: A¼lJ2R2  2;U¼A r3ð13^z2Þ;Ueff ¼h2 2r2þU 5: er¼r=r;eh¼h=h;ef¼eher 6: ex¼1;0;0½ |;ey¼0;1;0½ |;ez¼0;0;1½ | 7: Total energy:e¼eKþU¼v2 2l rþU 8: Generalized mean motion:m¼1 l2eðÞ 3=2 9: q1¼eh_ ex 1þeh_ ez;q2¼eh_ ex 1þeh_ ez 10: eX¼1 1þq2 1þq2 2 1q2 1þq2 2;2q1q2;2q1  | 11: eY¼1 1þq2 1þq2 2 2q1q2;1þq2 1q2 2;2q2  | 12: Generalized angular momentum:c¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r2Ueff p 13: Generalized velocity vector:t¼_ rerþc ref 14: Generalized eccentricity vector:g¼1 ltrtðÞer 15: p1¼geY;p2¼geX 16: X¼reX;Y¼reY 17: b¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1p12p2 2 q;a¼1 1þb 18: Generalized semi-major axis:a¼l 2e¼l m2  1=3 19: cos K¼p2þ1 ab1ap2 2  Xap1p2Y  20: sin K¼p1þ1 ab1ap2 1  Yap1p2X  21: L¼atan2ðsin K;cos KÞþ 1 abXp1Yp2 ðÞ J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4988 Algorithm 2. XcECI ¼r|;v| ½ |¼GEQOE2RVfunðm;p1; p2;q1;q2;L;l;J2;RÞ 1: Numerically solve Kepler’s equation:K L¼Kþp1cos Kp2sin K 2: eX¼1 1þq2 1þq2 2 1q2 1þq2 2;2q1q2;2q1  | 3: eY¼1 1þq2 1þq2 2 2q1q2;1þq2 1q2 2;2q2  | 4: b¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1p12p2 2 q;a¼1 1þb 5: a¼l m2  1=3 6: X¼aap1p2sin Kþð1ap2 1Þcos Kp2  7: Y¼aap1p2cos Kþð1ap2 2Þsin Kp1  8: r¼XeXþYeY 9: r¼a1p1sin Kp2cos KðÞ 10: ^z¼rð3Þ=r 11: A¼lJ2R2  2;U¼A r3ð13^z2Þ 12: _ r¼ffiffiffiffi la p rp2sin Kp1cos KðÞ 13: cos L¼X r;sin L¼Y r 14: c¼l2 m  1=3ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1q2 1q2 2 q;h¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c22r2U p 15: _ X¼_ rcos Lh rsin L 16: _ Y¼_ rsin Lþh rcos L 17: v¼_ XeXþ_ YeY References Amato, D., Bombardelli, C., Bau `, G., et al., 2019. Non-averaged regularized formulations as an alternative to semi-analytical orbit propagation methods. Celest. Mech. Dynam. Astron. 131, 1–38. Bau `, G., Hernando-Ayuso, J., Bombardelli, C., 2021. A generalization of the equinoctial orbital elements. Celest. Mech. Dynam. Astron. 133, 1– 29. Broucke, R.A., Cefola, P.J., 1972. On the equinoctial orbit elements. Celest. Mech. 5 (3), 303–310. Bruinsma, S., Thuillier, G., Barlier, F., 2003. The DTM-2000 empirical thermosphere model with new data assimilation and constraints at lower boundary: accuracy and properties. J. Atmos. Solar-terrestrial Phys. 65 (9), 1053–1070. DeMars, K.J., Jah, M.K., Schumacher, P.W., 2012. Initial orbit determination using short-arc angle and angle rate data. IEEE Trans. Aerosp. Electron. Syst. 48 (3), 2628–2637. ESA (2023). ESA’s space environment report 2023. URL: https://www. esa.int/Space_Safety/ESA_s_Space_Environment_Report_2023 accessed: 2024-01-10. Escobal, P., 1970. Methods of orbit determination. Krieger Pub Co.. Goff, G.M. (2015). Orbit Estimation of Non-Cooperative Maneuvering Spacecraft. Ph.D. thesis Air Force Institute of Technology. URL: https://scholar.afit.edu/etd/202/. Gomez, R., Salmero ´n, J.M.-V., Besso, P., et al., 2019. Initial operations of the breakthrough Spanish Space Surveillance and Tracking Radar (S3TSR) in the European context. In: In 1st ESA NEO and Debris Detection Conference. Presented paper. Darmstadt, Germany, URL: https://conference.sdo.esoc.esa.int/proceedings/neosst1/paper/479. Gooding, R., 1993. A new procedure for orbit determination based on three lines of sight (angles only). Defence Research Agency Farnborough, UK. Hatten, N., Russell, R.P., 2017. A smooth and robust Harris-Priester atmospheric density model for low Earth orbit applications. Adv. Space Res. 59 (2), 571–586. Hernando-Ayuso, J., Bombardelli, C., Bau `, G., et al., 2023. Near-linear orbit uncertainty propagation using the generalized equinoctial orbital elements. J. Guid., Control, Dynam. 46 (4), 654–665. Kaushik, A.S. (2016). A statistical comparison between Gibbs and Herrick-Gibbs orbit determination methods. Ph.D. thesis. Long, A.C., Cappellari Jr, J., Velez, C., et al., 1989. Goddard trajectory determination system (GTDS) mathematical theory (revision 1). National Aeronautics and Space Administration/Goddard Space Flight Center, FDD/552-89/001 and CSC/TR-89/6001. Ma, H., Bau `, G., Bracali Cioci, D., et al., 2018. Preliminary orbits with line-of-sight correction for leo satellites observed with radar. Celest. Mech. Dynam. Astron. 130 (10), 70. Maisonobe, L., Pommier, V., Parraud, P., 2010. Orekit: An open source library for operational flight dynamics applications. In: In 4th international conference on astrodynamics tools and techniques. European Space Agency Paris, pp. 3–6. Martinusi, V., Dell’Elce, L., Kerschen, G., 2015. Analytic propagation of near-circular satellite orbits in the atmosphere of an oblate planet. Celest. Mech. Dynam. Astron. 123, 85–103. Milani, A., Gronchi, G.F., Knezˇevic ´, Z., et al., 2005. Orbit determination with very short arcs: II. identifications. Icarus 179 (2), 350–374. Milani, A., Gronchi, G.F., Vitturi, M.D.M., et al., 2004. Orbit determination with very short arcs. I admissible regions. Celest. Mech. Dynam. Astron. 90, 57–85. Montilla, J.M., Sanchez, J.C., Vazquez, R., et al., 2023a. Manoeuvre detection in Low Earth Orbit with radar data. Adv. Space Res. 72 (7), 2689–2709. Montilla, J.M., Vazquez, R., & Di Lizia, P. (2023b). Mixture-based cost metrics for maneuver detection using radar track data. In 33rd AAS/ AIAA Space Flight Mechanics Meeting AAS 23–267 (pp. 1–16). Pastor, A., Sanjurjo-Rivo, M., Escobar, D., 2021. Initial orbit determination methods for track-to-track association. Adv. Space Res. 68 (7), 2677–2694. Picone, J., Hedin, A., Drob, D.P., et al., 2002. NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues. J. Geophys. Res.: Space Phys. 107 (A12), SIA-15. Reihs, B., Vananti, A., Schildknecht, T., et al., 2021. Application of attributables to the correlation of surveillance radar measurements. Acta Astronaut. 182, 399–415. Schutz, B., Tapley, B., Born, G.H., 2004. Statistical orbit determination. Elsevier. Shang, H., Chen, D., Cao, H., et al., 2018. Initial orbit determination using very short arc data based on double-station observation. IEEE Trans. Aerosp. Electron. Syst. 55 (4), 1596–1611. Siminski, J. (2016). Techniques for assessing space object cataloguing performance during design of surveillance systems. In 6th International Conference on Astrodynamics Tools and Techniques (ICATT) (pp. 14–17). Tao, X., Li, Z., Gong, Q., et al., 2020. Uncertainty analysis of the shortarc initial orbit determination. IEEE Access 8, 38045–38059. Vallado, D.A., 2001. Fundamentals of astrodynamics and applications, volume 12. Springer Science & Business Media. Vallado, D.A., Carter, S.S., 1998. Accurate orbit determination from short-arc dense observational data. J. Astronaut. Sci. 46, 195–213. Vananti, A., Schildknecht, T., Siminski, J. et al. (2017). Tracklet-tracklet correlation method for radar and angle observations. In 7th European Conference on Space Debris. ESOC Darmstadt volume 4. Zhang, S., Fu, T., Chen, D., et al., 2019. An initial orbit determination method using single-site very short arc radar observations. IEEE Trans. Aerosp. Electron. Syst. 56 (3), 1856–1872. J.M. Montilla et al. Advances in Space Research 74 (2024) 4968–4989 4989