Satellite maneuver detection and estimation with optical survey observations

Maneuver detection and estimation is deemed crucial for maintaining catalogs of Resident Space Objects (RSOs) as it helps to avoid sets of duplicated objects and track correlation issues. In fact, maneuvers, along with launches and break-up events, are the main source of potential new object detections during RSOs cataloging activities. For the continuous and reliable provision of Space Situational Awareness (SSA) and Space Traffic Management (STM) services, a challenging trade-off between detection time and characterization accuracy of maneuvers needs to be performed. In this paper, two novel and operationally feasible methodologies are proposed for maneuver detection and estimation. The first, a track-to-orbit methodology, uses a pre-maneuver orbit to linearize the dynamics and estimate the single burn that minimizes the residuals of the post-maneuver tracks. The second, an orbit-to-orbit methodology, estimates the double burn that solves a minimization problem between the pre-maneuver and post-maneuver orbits. Both methods, based on an optimal control approach, are not only proposed to tackle the maneuver estimation problem but also to be integrated on operational and robust association frameworks. Results are presented for optical scenarios with both simulated and real data, providing insightful conclusions on the capabilities, performance and limitations of the proposed methods. Particular emphasis is given to the importance of the track association, since a single track is usually not enough to perform a reliable estimation of the maneuver. Besides, the capability of the methods to provide a solution to the association problem, even when not perfectly characterizing the true maneuver, is discussed.


Introduction
The increasing number of Resident Space Objects (RSOs) and congestion of the orbital debris environment renders the space cataloging activities more challenging year after year. Currently, there are over 500 operational satellites only in Geostationary Earth Orbit (GEO) [1], most of which perform maneuvers every one or two weeks and for which Space Surveillance and Tracking (SST) systems predominantly have optical observations.
The main issue for SST systems related to maneuvering RSOs is the challenging correlation of observations. In the case of survey sensors, correlation is usually performed by comparing real and simulated observations, generated from predicted orbits of already cataloged RSOs. Unless these predicted orbits take into account the maneuver plan followed by the operators, correlation analysis of the first track after a maneuver will fail for sufficiently large magnitude maneuvers and after enough time to impact the orbit. Even for tracking sensors this is a problem, since the RSOs might not be located where expected, therefore leading to an observability issue. This may result in a loss of the RSO or ambiguous correlation situations, where the identity of the observed object is not clear. In those cases, although the correlation could be achieved, maneuver estimation is still required to properly update the state of the RSO through orbit determination. In fact, maneuvers represent currently he primary contribution to potential new detections (more than 500 operational satellites only in GEO [2], most of which perform maneuvers every one or two weeks in the case of chemical propulsion or even daily in the case of electric propulsion), exceeding those related to satellite launches (less than 400 spacecraft launched per year [2]) and break-up events (less than ten events per year, 98% of which involves less than 300 debris cataloged [3]), so detecting them is crucial for maintaining catalogs of RSOs. Capable maneuver detection and estimation methods are a must since otherwise these potential new detections would be promoted to actual new objects leading to sets of duplicated RSOs, thereby polluting the catalog and hampering the provision of Space Situational Awareness (SSA) and Space Traffic Management (STM) services.
At this point it may be convenient to clarify some terms that are extensively used along this paper and whose definitions may depend on the particular field of study. On the one hand, we use the term track to refer to a set of observations taken by a single sensor usually over a short time period, originated from the same RSO and frequently not enough to reliably estimate an orbit. Tracks (or tracklets) obtained from optical and radar surveillance sensors are usually referred to as uncorrelated optical observations (UCOs) and uncorrelated tracks (UCTs) when they cannot be associated to any cataloged RSO, respectively.. We will refer to them as tracks, regardless of the sensor type. On the other hand, a well-established estimation of the trajectory of an RSO in the catalog is called orbit. Accordingly, we refer to track-totrack (T2T), track-to-orbit (T2O) and orbit-to-orbit (O2O) as the association or correlation of tracks, tracks and orbits, and orbits, respectively.
Maneuver detection and estimation can be tackled as part of the association or correlation problem and should be integrated within the cataloging chain herein used, as depicted in Fig. 1, which shows a view of the different processes involved. This paper focuses on two methods for maneuver detection and estimation (located inside branches 1 and 2). In order to update the cataloged orbits and solve the T2O and O2O association problems, these two methods must be integrated in an association framework. The very first UCT received after a maneuver will most likely not be reliably correlated against any RSO in the catalog due to the velocity change and its effects on the dynamics, so it enters into the T2T association algorithm, intended for the detection of new objects. Initially, this first post-maneuver UCT cannot be associated to any other track in the T2T process. However, as more post-maneuver UCTs are obtained, these can be associated. If nothing is done to detect the maneuvers, the set of correlated tracks would coalesce in a new object. In order to prevent this, depending on the complexity of the maneuver, two main possibilities arise (branches 1 and 2 in Fig. 1): -Maneuver detection and estimation from uncorrelated tracks: the new UCT is first associated with the corresponding orbit of the RSO via T2O correlation considering a single-burn maneuver (branch 1 in Fig. 1). This allows to establish a first and preliminary link -association or hypothesis in the Multiple Hypothesis Tracking (MHT) framework -between an orbit and a single UCT, although the maneuver cannot be yet confirmed nor estimated reliably due to the scarce information available. It is important to note that not every RSO in the catalog is considered, but only a subset of candidates, such as those identified as active satellites that have not been recently updated. The maneuver detection and estimation should be performed in the observations space by means of a T2O methodology including maneuvers as described below. The rationale behind this is that the orbital estimation derived from the new UCT (or the few associated UCTs) is still not reliable enough to be directly used. As more UCTs after the maneuver become available to the system, new associations of more tracks arise until there is enough information to promote, i.e. confirm, the hypothesized maneuver. The number of tracks for the association required to properly confirm and estimate a single burn maneuver (i.e., promote a hypotheses) is expected to be around two or three, as suggested in [4]. In any case, less than the four required for a full and nominal RSO initialization in the catalog [5,6]. The proposed T2O methodology is able to estimate maneuvers based on residuals between the estimated orbit before the maneuver and observations afterwards. Figure 2 shows the residuals of observations from four telescopes of the International Scientific Optical Network (ISON) and a GEO satellite before and after a North-South maneuver (vertical dashed line). The more time after the maneuver elapses, the greater the divergence of the residuals. This is an indication of the footprint of the maneuver on the residuals of the post-maneuver tracks when considering the pre-maneuver orbit. -Maneuver detection and estimation from potential new object: the information contained on a small number of tracks might not be enough to estimate the parameters characterizing a maneuver of two burns, since the T2O methodology (branch 1 in Fig. 1) would not be able to associate the orbit with the postmaneuver UCTs. Therefore, it is required to associate a higher number of postmaneuver tracks to obtain a new and reliable orbit estimation without the use of prior information, i.e. perform a potential new RSO detection by means of T2T association. The number of tracks required is greater than in the previous situation, since a full RSO initialization needs to be performed. Once accurate postmaneuver orbital information is available, the maneuver detection and estimation can be done in the orbit space by means of an O2O methodology considering maneuvers (branch 2 in Fig. 1). This problem corresponds to the estimation of two maneuvers capable of linking two already well-established orbits. The case of low-thrust maneuvers may pose challenges to this approach. Nonetheless, the methodology is expected to remain applicable as long as the low-thrust maneu- Fig. 2 Residuals between the estimated orbit before the maneuver and measurements after the maneuver ver has finished, although the estimated impulsive burns would not be very representative of the actual low-thrust maneuver performed.
Multi Target Tracking (MTT) methods, traditionally applied to sensing, guidance, navigation and air traffic control, among others [7], have been also used to tackle the T2T association problem. One can find the work of Stauch et al. [8] based on Joint Probabilistic Data Association (JPDA) [9] and Constrained Admissible Regions (CAR) [10], and the one by Pirovano et al. [11], still based on JPDA and CAR but with a different treatment for the uncertainty following Differential Algebra (DA). JPDA consists in the joint association of objects and observations disregarding the one-to-one association constraint, thus alleviating the decision-making process. Alternatively, Aristoff et al. [12] propose a MHT approach, in which only one-to-one associations are allowed, yet maintaining multiple parallel hypotheses to be tested as new observations arrive. Jones and Vo [13] suggest utilizing a purely statistical framework for RSO catalog maintenance based on Random Finite Sets (RFSs), thus accommodating uncertainty in the number of objects. Beside decision making and statistical association methods, Siminski et al. [14] discuss on the particular application to GEO objects and optical observations. Finally, Furfaro et al. [15] apply machine learning methods to study the variability in space object catalogs and identify feasible trends. However, none of these methods is capable of dealing with maneuvering objects, which may lead to object duplication or false associations if the maneuver exceeds the association uncertainty or likelihood.
The automatic detection of maneuvering RSOs can be framed within the general multi-target tracking-association problem, in which the tracked objects are allowed to maneuver. Some approaches [16,17] take digested information from Two Line Elements (TLEs) to isolate, and thus detect, maneuver intervals. Moreover, a proper treatment of historical data may be used to integrate maneuver detection capabilities within the association problem, i.e. in an on-line fashion. In this regard, Siminski et al. [18] present an implementation based on actual maneuver data reported by the operator, which employs a kernel density estimator to cluster maneuvers of different nature and define admissible regions. Alternatively, Shabarekh et al. [19] developed a machine learning approach aimed at determining the Patterns of Life (PoL) in order to predict and characterize maneuvers. In fact, these techniques [18,19] can be readily integrated in an operational framework as they are devised to inherit real measurement data.
Maneuver estimation has also been assessed by means of optimal control methods. For instance, Holzinger and Scheeres [20] and Holzinger et al. [21] developed an alternative method to JPDA based on distance control metrics as opposed to the usual Mahalanobis distance [22], which is aimed at determining the minimum required control effort to fit a given observation. On the same topic, Lubey [23] emphasizes on the ability to jointly consider maneuver detection and data association, contributing with a thorough uncertainty characterization. These methods characterize the maneuver as a continuous perturbing acceleration of varying magnitude, whereas an impulsive-based approach may also be embraced. The work in [24,25] suggests utilizing the concept of the state transition matrix to evaluate the impact of an impulsive maneuver in the form of a perturbation. A minimization problem can then be formulated in order to derive the parameters defining the maneuver. Note however, the authors only consider an O2O scenario and use a simplified maneuver model: either one or two burns in the direction of the orbital velocity. An extension to [24,25] is the work by Yang et al. [26], which discuss on the uncertainty propagation of maneuvering objects by means of state transition tensors of order four, as opposed to the first order approach based on the use of transition matrices. However, this first order approach has proved to be accurate enough to estimate large and small maneuvers [4,27] and have been also applied to relative motion dynamics [28].
Most of the aforecited methods must be tested and tuned for each specific application. The characteristics of the maneuver detection and estimation problem, especially considering scarcity of data and large time intervals between tracks makes the tuning demanding. In previous years, a series of works have tried to address the gaps of the classical approaches with methods developed adhoc for the SST problem [18,23,27]. In those works, RSO maneuvers are characterized a-priori, in order to incorporate more information to the problem, and make it more tractable. Alternatively, thrust Fourier coefficients have been used to estimate equivalent maneuvers with the same secular behavior and without the need of any a-priori information [29]. Nevertheless, most of them provide results for particular and isolated test cases, rather than representative enough scenarios, analogous to an operational SST system. Scalability concerns, as in the track association problem [5], is of a major importance, given the dimensions of the problem and the huge number of a-priori possible combinations to evaluate.
To address the above-mentioned gaps, we propose two maneuver detection and estimation methods to be used in operational and robust cataloging chains. They do not require any a-priori information of the maneuver and are able to provide an estimation to update the cataloged orbit of the involved RSO. To do so, uncorrelated tracks are associated among them (T2T) and then cataloged orbits are correlated against these post-maneuver track associations (T2O including candidate maneuvers of both single and double burns). In this way, associations of orbits and tracks, or hypotheses, are generated, scored, pruned and promoted in such a way that the involved cataloged orbit and sensing data belongs to a common maneuvered RSO. For instance, maneuvers involving excessive control effort are discarded since they are not realistic but solutions to the orbit linkage problem. Note that in this case the maneuver detection and track correlation problem are coupled. Solving this coupled problem yields both the object identification for which the sensor data is generated as well as the maneuver time, direction and size. This would allow to considerably reduce the maneuver detection time, since T2O association, as opposed to O2O, does not require a new object detection and initiation. However, an alternative formulation is proposed to consider O2O correlation scenarios in which two established objects may be reduced to a single maneuvered one. Opposed to the single impulsive burn of the T2O correlation, this last O2O scenario assumes a double impulsive burn to allow for a transfer orbit capable of approximating maneuvers in a more robust manner.
The present paper is structured in five sections. Section 1 has introduced the role of maneuver detection and estimation on cataloging activities, a state-of-the-art review and the framework of our methodologies. Sections 2 and 3 present the T2O and O2O methodologies for single and double burns maneuver estimation, respectively. Section 4 shows and analyzes the results with both simulated and real observations. Finally, Section 5 gathers the conclusions of the paper and discusses the current status of the work.

Single burn maneuver detection and estimation via T2O
The T2O methodology is proposed to solve the maneuver detection and estimation problem when a pre-maneuver orbit and post-maneuver tracks are available. On the one hand, we assume a pre-maneuver orbit (subscript A), an orbit estimated before the maneuver, is available on the catalog. This means that an extended state vector: T ∈ ℝ 6 and p A ∈ ℝ n p represent the state vector and dynamical parameters, respectively, being r A ∈ ℝ 3 and v A ∈ ℝ 3 the corresponding position and velocity vectors. On the other hand, a set of N optical observations, z t i ∈ ℝ 2 for i = 1, … , N , has been received by the sensor network. Each observation contains a pair of right ascension and declination measurements referred to t i , packed in tracks so that each track contains only observations from a common sensor over a short time period. These N observations may come from n T ≤ N optical tracks, which shall have been previously associated with other methods, such as the association framework presented in [5].
The maneuver can be detected by inspection of the residuals of the post-maneuver tracks and the pre-maneuver orbit (see Fig. 2). The divergence of the residuals can be found by setting a threshold on the residuals weighted with the expected measurement noise of the involved sensor. In this way, the maneuver detection and estimation is triggered whenever absolute or relative threshold criteria are met.
Let t, t 0 be the full transition matrix, allowing to propagate perturbations of the extended state vector from t 0 to t under the linear dynamics assumption, i.e.: Note that t, t 0 contains the so-called state transition, and sensitivity matrices [30]: Subsets of these matrices will be referred to as (⋅) , corresponding to (t)∕ t 0 . At a certain epoch, t M , an impulsive maneuver takes place causing a sudden change in the velocity, i.e.: u = v B t M − v A t M , as shown in Fig. 3. Since the position of the two orbits is intersecting at t M , i.e.; r A t M = r B t M , the post-maneuver state (subscript B) can be obtained by considering the maneuver a perturbation at time t i as follows: Then, the maneuver magnitude, u ∈ ℝ 3 , for a given t M , is estimated so that the residuals of the observations, i = z t i − h t i , y t i , difference between actual measurements and measurements reconstructed from the post-maneuver orbit, perturbed with the solve-for maneuver, is minimized, i.e.: being W the weighting matrix, which contains the inverse of the expected variance of each measurement. The solution can be obtained via a weighted non-linear leastsquares method. To do so, the problem is linearized around a reference point, u 0 , and the corresponding correction, u , is obtained by solving the following linear system: where G is the Jacobian, i.e.: partials of the measurements with respect to the estimated parameters, and z is the difference between the actual observations and the observations predicted from the reference trajectory. This expression is similar to the classical orbit determination problem: H T WH y = H T W z [30], in which the residuals are minimized, being H the Jacobian with respect to the state vector in this case (partials of the measurements with respect to the state vector). In the problem at hand the Jacobian G is computed with respect to the parameters of the single burn maneuver u as opposed to y in H. Note that the state vector at each observation epoch, corresponding to the postmaneuver orbit and required for the evaluation of the Jacobian and residuals, can be linearly propagated by means of Eq. 5.
Starting from a null initial solution and iteratively solving Eq. 7 it is possible to obtain an estimation for u . The Jacobian should be updated at each iteration to ensure convergence. Then, the contribution to G of the i th observation is: In order to detect the maneuver, the least-squares problem must be solved for a range of t M values. In principle, the maneuver is assumed to have occurred after t + A , the epoch of the last observation considered for the estimation of the pre-maneuver orbit, and t 1 , the epoch of the first available new observation. This thus makes it desirable to have a very computationally efficient estimation method, which is the rationale behind using a linearized post-maneuver orbit propagation, rather than a fully numerical one.
Finally, a set of estimations {u k } and corresponding objective function values {J k } for each t k ∈ T for which the problem could be solved are to be obtained. Since the final goal of the maneuver estimation problem presented is the T2O association, we suggest taking every local minima of J and then selecting the solution as the one leading to minimum |u| (where |⋅| denotes the Euclidean norm when applied to a vector). Additional constraints can be included, such as introducing a u max value to avoid the consideration of unrealistic maneuvers. Moreover, the n solutions leading to minimum |u| can be retained for a MHT approach. Note that the methodology provides a set of compatible solutions and their corresponding score (i.e. maneuver magnitude), and as such can be readily integrated in an association framework to solve the T2O correlation problem. The complete method is summarized in Algorithm 1, where hat denotes estimated values.

Double burn maneuver detection and estimation via O2O
The O2O methodology is proposed to solve the maneuver estimation problem when a pre-maneuver and post-maneuver orbits are available. In this case, we assume both the pre-maneuver orbit, y A , and the post-maneuver orbit, y B , are given. The estimation process can be triggered when a residuals divergence is identified, analogously as in Section 2. Besides, this methodology should be applied before publishing the second orbit in the catalog to avoid duplicated objects. The problem is then to find the transfer orbit, x T , that connects the pre-maneuver and post-maneuver orbits, i.e.: where u i represents the i th burn, i.e., delta-V, to be estimated.
The problem, illustrated in Fig. 4 can be solved for t M1 , t M2 ∈ T such that t + A < t M1 (i.e.: the first burn occurs after the last observation used in the estimation of the premaneuver orbit) and t M1 < t M2 < t − B (i.e., the second burn takes place after the first burn and prior to the first observation used in the estimation of the post-maneuver orbit).
Although the solution could be determined by solving Lambert's problem between r A t M1 and r B t M2 and then recovering the maneuver magnitudes directly from Eqs. 9 and 10, we propose a method that uses linearized dynamics to propagate the effect of the maneuvers (as perturbations) in the orbits. This way, we avoid the limitation of the two-body motion dynamical model (i.e. Kepler's law) by considering relevant perturbations while keeping the computational effort low. Two equations arise from the propagation of the pre-maneuver and post-maneuver orbits, respectively: Although either Eqs. 13 or 16 could be directly solved to obtain u 1 and u 2 , we suggest combining the two linear systems to obtain the following problem: The solution of this overdetermined system (12 equations, 6 unknowns) can be obtained via least-squares, i.e., by solving the following linear system: where X is the left-hand-side matrix and the right-hand-side vector in Eq. 18. u = {u 1 , u 2 } T can be solved with any factorization method such as Cholesky decomposition.
Since this method is based on a linearization of the dynamics of the two orbits, there is an inherent applicability limitation related to the magnitude of the perturbations above which the linear dynamics assumption is expected to fail.
The solution of Eq. 18 is expected to be a smooth combination of the two dynamics. Moreover, note that the contribution of each orbit could be weighted if required with a confidence level or covariance, for instance. The complete method is summarized in Algorithm 2.

Results
The T2O and O2O methods, presented in Sections 2 and 3, have been applied to scenarios with simulated and real observations of GEO RSOs. These scenarios are briefly introduced next: propulsion. This scenario intends to study the performance of the methodologies when the impulsive maneuver hypothesis does not hold and the duration of the burn is of several hours.
The real scenarios use real observations from telescopes of the ISON sensor network covering maneuvers of GEO operational satellites used for communications and meteorological purposes (semi-major axis, eccentricity and inclination are included in Table 12 for reference). Besides, Figs. 27 and 28 present the distribution of the track duration (time between first and last observation of the track) and observation spacing (time between observations in a track) for all the tracks considered in the real scenarios. A total of 26 telescopes, listed in Table 13 and featuring a measurement error ranging from 0.5 to 1 arcsecond, have been considered in the analysis. Regarding the real maneuvers, reference values have been generated using additional external data (satellite operator maneuver plans and ranging observations with few meter-level errors on a post-maneuver orbit determination) and are considered as truth values.
Regarding the T2O methodology, the pre-maneuver orbit is estimated by using the observations available before the burn via a batch least-squares orbit determination and propagated to the future to cover the post-maneuver tracks. Then, the post-maneuver tracks are associated incrementally via T2T association [5]. Finally, the T2O method, presented in Section 2, is applied to each pair of premaneuver orbit and post-maneuver track association. A linear grid of maneuver epoch, t M , values has been considered in all the results presented, with a time step of 1 hour in the simulated case (Sat#0), 15 minutes in the first real case (Sat#1) and 6 minutes in the rest of real cases (Sat#2, Sat#3 and Sat#4). The sequence of events is depicted in Fig. 5.
In the case of the O2Omethodology, the pre-maneuver and post-maneuver orbits are estimated analogously, by means of an orbit determination with the observations available before and after the burns, respectively. Then, the O2O method, presented in Section 3, is applied to the pair of orbits. The sequence of events is depicted in Fig. 6.
The dynamical model considered for the state propagation, usual in the case of GEO RSOs, consists in a 30x30 Earth gravitational field, Moon and Sun third body perturbations and cannonball model for the Solar Radiation Pressure (SRP). Regarding the orbit determination, the state and the SRP coefficient are estimated.
Finally, the computational time of each case presented in this section is presented in Table 14. Note that in Sat#1, Sat#3 and Sat#4 scenarios the T2O methodology is applied to several associations of one, two, three and four tracks, so an average time is presented. Besides, the cases in which the T2O methodology is applied to all tracks (suffix A) is included for reference, but do not represent a typical case of application. The computational times correspond to single thread executions on a 2.60 GHz Intel(R) Xeon(R) Gold 6142 CPU.

Sat#0 scenario: simulated observations
The first scenario consists in a set of simulated single and double burn maneuvers, intended to validate the T2O and O2O methodologies, respectively. The involved Fig. 5 Sequence of events for applying the T2O methodology object, Sat#0, is located on a GEO orbit (see Table 1 for further details) and representative of a typical GEO satellite.

Single burns, T2O methodology
The first assessment of the T2O methodology is aimed at investigating the impact of the number of tracks and performing a preliminary analysis with the method under a simulated scenario. To do so, the initial state vector (Table 1) was propagated from t 0 up to t 0 + 7 days considering an impulsive burn in the local RIC frame (radial, intrack and cross-track) at the middle of the interval. Four cases have been studied: 1) radial burn, 2) in-track burn, 3) cross-track burn and 4) additional cross-track burn of high magnitude. The magnitude of the maneuvers is |u| = 1.0 m∕s for the three first cases and |u| = 10.0 m∕s for the last case.
On the one hand, a pre-maneuver orbit was estimated with simulated observations before the maneuver and propagated forward in time without considering the maneuver. On the other hand, an optical sensor station was simulated to generate three tracks on the 5 th day (12 h after the maneuver), 6 th day (36 h after the maneuver) and 7 th day (60 h after the maneuver), according to the reference orbit. Each track has a duration of 15 min and contains one observation (pair of right ascension, , and   J that not always coincide, with a time separation of around 1 day, which is the orbital period. This non-linear behavior suggests that a joint estimation of both the maneuver epoch and magnitude may not be a good choice, at least for an initial maneuver estimation approach. The details of each √ J local minima are compiled in Table 2, where apart from the maneuver estimation results, the errors in the semi-major axis, eccentricity and inclination estimations (corresponding to the pre-maneuver orbit after the application of the estimated maneuver) are also shown. Firstly, in the case of the radial burn, the solution with lowest |û| does not exactly correspond to the true solution, although the three local minima differ in less than 0.2% with respect to the true magnitude. Besides, the one corresponding to the true solution has the lowest error in semi-major axis and eccentricity. Secondly, in the case of the in-track burn, there are two solutions on the vicinity of the epoch (one hour before and one hour after). They have the lowest |û| , as well as lowest orbital differences. Thirdly, in the case of the cross-track burn, there are two solutions with lower |û| than the one of the true maneuver and the one corresponding to the true solution has the lowest orbital differences. Finally, in the case of the cross track burn of high magnitude, there are solutions with √ J > 1 , unlike in the rest of the cases, because of the high magnitude of the burn. This last case was included as a limiting one to study potential limitations of the linearization. As opposed to the cross-track burn of 1 m/s, the estimated maneuver magnitude in the in-track direction is not negligible, thus making the error in the semi-major axis to increase. However, the solution with lowest √ J corresponds to the true solution and the total estimated maneuver magnitude is of 9.95 m/s, i.e.: less than 0.5% error with respect to the true value.
It is expected that cases involving three tracks, {1, 2, 3} , are the best conditioned ones. Figure 10 shows the |û| and √ J distribution of the √ J local minima found for each combination of tracks, along the semi-major axis and eccentricity errors in the radial burn case. Although several √ J local minima are found when considering only a single track, even with |û| ∼ 1m∕s and √ J ∼ 1 , more information (tracks) is required to reliably estimate the effect of the maneuver on the orbit. Only when two and three tracks are involved, solutions with |a −â| < 10 m and |e −ê| < 10 −5 are found and the local minima converge to the truth. In other words, although the methodology can be used with a single track, the lack of information in that case may lead to an insufficiently accurate orbit estimation, potentially compromising subsequent cataloging activities. The methodology allows the linkage of pre-maneuver and post-maneuver tracks, thus reducing the number of tracks required after the maneuver if compared to a new detection from scratch (three or four tracks are usually required [4]) that considers only post-maneuver tracks.

Double burns, O2O methodology
The first assessment of the O2O methodology is aimed at performing a preliminary analysis of the method and comparison against Lambert's problem solution.
To do so, an initial state vector was propagated from t 0 up to t 0 + 7 days considering two impulsive burns in the local RIC frame. As in Section 4.1.1, this orbit will be referred to as reference orbit and six cases have been studied: 1) two radial burns (RR), 2) two in-track burns (II), 3) two cross-track burns (CC), 4) radial and in-track burns (RI), 5) radial and cross-track burns (RC) and 6) in-track and cross-track (IC) burns. The simulated burns are u 1 = +0.1 e i m∕s and u 2 = −0.1 e i m∕s , being e i the unitary vector in the R, I or C direction, i.e.: the sense of the first burn is positive, while the one of the second is negative with respect to the RIC frame. The burns are simulated at 12:00 of the 3 rd and 4 th days. As before, the pre-maneuver orbit was generated by performing the same propagation as the reference orbit but without considering the maneuvers. On the other hand, the post-maneuver orbit was generated by performing a back-propagation of the last state vector of the reference orbit without considering the maneuvers. The initial state vector is the same as in the T2O simulated scenario. In this case, the outcome of the maneuver estimation method are the two estimated burns, û 1 and û 2 , for each t M1 and t M2 in the mesh used to sample T . Figure 11 shows the distribution of the total velocity increase, i.e.: | |û1 | | + | |û2 | | , along t M1 ,t M2 ∈ T in the II (left) and RI (right) cases. Note that solutions with |û| > 5 m∕s have been discarded. This representation, known as porkchop plot and typically used for interplanetary transfers, provides the velocity increase required for each combination of t M1 and t M2 that connects the pre-maneuver and post-maneuver orbits. Since there are no metrics to select a candidate estimation besides the maneuver magnitude, the optimal maneuver, i.e. minimum |û| = | |û1 | | + | |û2 | | is accepted as solution. However, in Fig. 11, the presence of multiple local minima confirms the existence of maneuvers with similar control effort that connect the two orbits. In principle, they are equivalent from the association point of view, meaning that an accurate estimation of the maneuver epoch is not ensured.
The three lowest |û| of each case are presented in Table 3. In most cases, the corresponding t M1 and t M2 do not match the true values and |û| < |u| , meaning that a maneuver with lower |u| than the true one has been found. The norm of the overdetermined linear system, , suggests that the linear system has been properly solved. As expected, those cases involving an in-track burn present higher values of (although within acceptable bounds), being this an indicator of the higher nonlinearity in this direction.
In order to evaluate the accuracy of the estimations, the resulting maneuver magnitudes at t M1 and t M2 (true values) for each case are summarized in Table 4. As expected, the estimated values of the burns coincide with the true ones, even if considered independently. This shows how the methodology is able to find a solution that matches the true maneuver, although it may not be selected as lower |u| ones are found. This is expected, since there are more variables driving operator's maneuver plans, such as working hours, holidays or frequency, among others. Given the unavailability of this external data, sticking to lowest |u| solution is suggested, even though in some cases it may not correspond to the true one. In any case, both are expected to solve the linkage problem and as soon as more tracks after the maneuver are received, the maneuver estimation (not only magnitude on each direction but also time) could be refined.

Sat#1 scenario: real observations and SK maneuvers
Once the performance of the methodologies with simulated data has been assessed, a scenario with real observations is studied to confirm the adequacy of the proposed method. It consists in a real GEO satellite providing coverage to Europe and 4 telescopes from the ISON network. The tracks are distributed along two weeks and the impulsive burns performed by the RSO are: 1) | | u 1 | | ∼ 1.17 m∕s NS and 2) | | u 2 | | ∼ 20 mm∕s EW burn, separated around 62 h. The timeline of tracks and burns is shown in Fig. 12.

Single burns, T2T methodology
In order to understand the details of the T2O methodology performance, a subset of 9 tracks, depicted as red dotted lines in Fig. 12, from one of the telescopes has been used.

NS burn
Five tracks, {1, 2, 3, 4, 5} , were available after the first burn and before the second burn and associations of two and three tracks have been considered for  J local minima half orbital period before and after the reference t M . These correspond to different maneuvers that are compatible with the orbit and tracks considered and in fact they present similar |û| values.
The local minima corresponding to each solution with |û| < 2 m∕s are presented in Table 5. The estimations are consistent along the different track associations and  J ∼ 1 indicated that the maneuver is estimated in such a way that the residuals of the observation matches the expected sensor noise. Note that associating three tracks allows to discard a pair of local minima with √ J ∼ 3.6 , while this is not possible if only associations of two tracks are used. The solution with lowest |û| from the associations of three tracks, that would be identified as the the maneuver in a hard-decision making environment, would be |û| = 1.1512 m∕s at t M = 12-18:30 (from association {1, 3, 5} ), i.e.: 1% of error in magnitude and around 12 h difference in epoch with respect to the reference value. This solution correspond to a burn with opposite sense in the cross-track direction performed on the opposite orbital point (around half period separation).

EW burn
In this case, four tracks, {6, 7, 8, 9} , were available after the second burn and again, associations of up to three tracks have been considered for the estimation of the maneuver. Figure 14 shows the distribution of every |u| local minima found along |u| and √ J . Even though the estimation of this burn is more challenging that the previous NS one due to the lower impulse involved, which is translated into a fainter maneuver footprint on the residuals, associations of two and three tracks are able to estimate a burn within the order of magnitude of the reference value. The estimation obtained with the association of three tracks and lowest |û| is |û| = 44 mm∕s and the corresponding maneuver epoch is estimated with less than one hour of error ( t M = 14-21:15).

Double burns, O2O methodology
The pre-maneuver and post-maneuver orbits, estimated with real observations, have been used to study the performance of the O2O methodology. Figure 15 shows the resulting porkchop plot (zoomed in the vicinity of the reference solution and discarding solutions with |û| > 5 m∕s ) that, as in the simulated scenario, present many local minima corresponding to maneuvers able to link the two orbits. Table 6 presents the four optimal solutions (i.e.: lowest |û| ). Note that | |tM1 − t M1 | | < 15 min and | |tM2 − t M2 | | < 2 h in the case of the optimal solution. Although this maneuver is different than the reference one in terms of | |û1 | | and | |û2 | | , if independently considered, the total velocity increase, |û − u| is estimated with an error lower than 1.5% of the reference value. This means that two potential burns whose |û| is of the same order of magnitude of the real one can be estimated. Besides, Table 6 present the results obtained if Lambert's problem is solved in the same t M1 and t M2 grid than the methodology. The reason behind the unaccurate results provided by Lambert's problem is the orbit dynamics mismodelling of the Lambert's problem solution (two-body motion). This justifies the choice of the linearized orbit model including perturbations over a two-body motion model for a robust approach.
The total estimated maneuver magnitude at t M1 and t M2 (reference epochs) is |û| = 1.29 m∕s , i.e., a relative error of less than 10%. Again, without additional information it is not possible to select this solution since the local minima have similar values of |û| and are able to solve the linkage problem. A projection of this distribution on the t M2 −t M1 plane (Fig. 16) illustrates this fact and also justifies the suitability of defining a threshold u max such that |û| < u max to reduce the number of solutions. In this case, there are a total of 96,141 solutions, which is reduced by 20% with u max = 5 m∕s and 54% with u max = 2 m∕s . Note this is very relevant for the association framework since it provides a mean to prune hypotheses, thereby reducing the computational load. Still, the computational burden associated to the maneuver detection and estimation methodology is not reduced, since these solutions have to be obtained first.

Sat#2 scenario: real observations and SK maneuvers
This scenario focuses on a real GEO satellite and 9 telescopes from the ISON network. The tracks are distributed along one and a half weeks, as shown in Fig. 17, and the impulsive burns performed by the RSO are | | u 1 | | ∼ | | u 2 | | ∼ 0.59 m∕s , in the in-track direction and separated 12 h. Since there are no observations between the two burns, only the O2O methodology has been applied. The pre-maneuver and post-maneuver orbits have been estimated by using all tracks before the first burn and after the second burn, respectively.

Double burns, O2O methodology
The porkchop plot, shown in Fig. 18, shows a clear |û| local minimum on the vicinity of the reference solution (dotted lines). Note that the shape of this porkchop plot is different than the one obtained for the double NS+EW burn (Section 4.2, Fig. 15) mainly because different directions are involved. There are less local minima and higher convexity in the current case, which, in principle, may lead to an easier identification of the solution that is closer to the reference solution. Besides, the distribution of |û| along the time of flight, t M2 −t M1 , is presented in Fig. 19. There are three families of solutions centered at 0.5, 1.5 and 2.5 days (periodicity directly related to the orbital period) that correspond to the three regions around the three local minima in Fig. 18.
The estimated values of the solution with lowest |û| are presented in Table 7, along with the reference values. Both the burns independently (relative error of 0.08% and 1.46%) and the total impulse (relative error of 0.7%), as well as the epochs, can be properly determined. Apart from the global minima (first estimation presented), the two |û| local minima corresponding to the other two families of solutions have been included (second and third estimations presented) for completeness. They are feasible solutions with similar |û| but different t M1 and t M2 . Besides, Table 7 lists also the solution obtained using Lambert's formulation, which exhibits a larger error when compared to the proposed method. In spite of the latter, this error is lower in relative terms than the one presented in Table 6 (Sat#1 scenario), since the dynamics mismodelling becomes more important as the time of flight increases (62 h in Sat#1 and 12 h in Sat#2).

Sat#3 scenario: real observations and large maneuvers
This scenario focuses on a real GEO satellite and 19 telescopes from the ISON network. The tracks are distributed along three days, as shown in Fig. 20

Single burns, T2O methodology
As opposed to Section 4.2 (Sat#1), where the focus was on subsets of tracks, now we have considered every available track and study the distribution of the estimated maneuvers for each association of 1, 2, 3 and 4 tracks. Figure 21 shows the distribution of the track associations along the estimated maneuver magnitude error, |û| , and the estimated maneuver epoch error, Besides, the vertical dashed black line represents the reference epoch of each burn. There is a clear benefit of increasing the number of associated tracks in terms of both epoch and magnitude estimation error. This is expected as more associated tracks imply more observations that are taken into account in the estimation process. Figure 22 presents the histograms of these errors for the particular case of the 3 rd burn. There is a clear improvement in terms of the estimated epoch and magnitude error when moving from one or two tracks to three or four. Note that the O2O methodology could be used when dealing with associations of four tracks since in this case a reliable post-maneuver orbit can be usually obtained. However, the T2O methodology provides not only the estimation of the maneuver (magnitude, direction and epoch) but also √ J , which can be used to prune solutions, of interest for the association problem. For each of the three burns, the estimation obtained from associations of 3 and 4 tracks leading to the minimum maneuver magnitude has been selected and compared against the reference one in Table 8. The maneuver magnitude relative error (selected assoc.) is around 0.6% for the first burn, 1% for the second burn and 19% for the third burn. They are good results taking into account the relatively short time interval of observations considered for the estimation of the pre-maneuver orbit, as well as the short time between the maneuver and the post-maneuver tracks, just a few hours. In this regard, we may assume that, as long as the linearization of the dynamics is valid, the greater the time between the maneuver and the post-maneuver tracks the better, because of the time needed by the dynamics to make noticeable the maneuver effect. Besides, Table 8 includes an additional estimation obtained by considering every available track (a total of 27, 31 and 16 for the first, second and third burns, respectively).  Apart from the expected improvement of the maneuver magnitude estimation, the main benefit of considering more tracks is to discard |û| local minima that, although similar, do not correspond to the true solution (as discussed in Section 4.2.1).

Double burns, O2O methodology
The O2O methodology has been applied to the same three burns but aimed at estimating two double burns: 1) 1 st + 2 nd burns and 2) 2 nd + 3 rd burns. To do so, the corresponding pre-maneuver and post-maneuver orbits have been estimated from the available observations.
The two resulting porkchop plots are shown in Fig. 23, including the reference epochs as black dashed lines. In both cases there is a clear |û| global minima region and, opposed to the previous scenarios, there are no additional regions with local minima, mainly due to the relatively short domain of search, limited by the maneuvers timeline. Although the global minima of the first case ( 1 st and 2 nd burns, left) does not perfectly match the reference epochs, t M1 and t M2 , note that there is a delay of ∼ 2 h between the estimated and reference maneuvers.
Moreover, the |û| global minima are listed in Table 9. As in other scenarios, the error of the estimated maneuver magnitude of the burns, if considered independently (26% and 4% relative error for the first case, 35% and 106% for the second) is higher than the total impulse (11% for the first case, 15% for the second). These results confirm that the O2O methodology is able to deal with high impulses for GEO regime, confirming the validity of the linearization under high perturbations.

Sat#4 scenario: real observations and electric propulsion
This scenario focuses on a real GEO satellite with electric propulsion and 12 telescopes from the ISON network. As opposed to the previous satellites, equipped with chemical propulsion, this one performs continuous maneuvers with a low thrust electric propulsion. The tracks are distributed along several days, as shown in Fig. 24 and the continuous burns performed by the RSO are: 1) | | u 1 | | ∼ 48 mm∕s with a duration of 1.65 min and 2) | | u 2 | | ∼ 135 mm∕s with a duration of 1.75 h, separated around 23 h (from start of the second to end of the first).

Single burns, T2O methodology
The T2O was applied to the estimation of the two burns, independently. Figure 25 shows the distribution of |û| and √ J along t M for one of the associations of 4 tracks. Even though the burns are continuous, the method is able to locate local minima on √ J assuming impulsive ones.  The selected solutions for each case ( √ J local minima with lowest |û| ) are listed in Table 10, where the reference values and the estimation using all tracks are included. Note that there are only 4 available tracks for the 2 nd burn estimation and thus, the selected association contains all tracks. The first burn estimated magnitude relative error is around 29% (22% if considering all tracks) and the estimated epoch differs less than 3 h from the reference. However, the error in the second burn is much higher than the first (around 1400% in magnitude and around 2 h). Note that Table 10 (right) shows data from 04-20:00 for completeness, although the last observation used for the pre-maneuver orbit estimation is on 05-00:50, i.e.: solutions before this epoch should not be considered.  The main difference between the two burns is the duration: 1.65 min and 1.75 h, respectively and according to the reference. The first burn duration is so short that it can be assumed impulsive, while this does not hold for the second burn. On the other hand, the characterization of the second burn is not as accurate as the first one. Despite of this, the association of the pre-maneuver orbit and the post-maneuver tracks, is, in principle possible, since a compatible maneuver has been found with a reasonable matching between the pre-maneuver orbit and the post-maneuver tracks ( √ J ∼ 1).

Double burns, O2O methodology
The double burn maneuver has been estimated with the O2O methodology, obtaining the porkchop plot presented in Fig. 26, where the two burns start and end epochs are depicted with black dashed lines. There are several local minima on the vicinity of the two reference impulsive burns. Table 11 lists the reference and estimated values (global minima) of the two burns. The first and second burn estimated epoch error is of less than 1 h and 9 h, respectively, while the estimated total magnitude relative error is of 40%. Since the two maneuvers are jointly estimated, there is not a significant difference between the two estimated impulsive burns, i.e.: the epoch of the first is estimated with lower error, while the magnitude of the second is estimated more accurately. At this point, it is important to recall that the O2O methodology is intended to estimate the double burn impulsive maneuver that best approximates the real maneuver, which in this case is not a double burn impulsive maneuver. Therefore, this is a an interesting test case with which the limitations of the methodology can be analyzed. In this case, the epochs and magnitude of the two separated burns are presented in Table 11 for completeness, but what is more relevant is the total magnitude (last two rows), which is estimated as 0.2572 m/s. The reference value is 0.1828 m/s, meaning that the obtained maneuver is of the same order of magnitude (around 40% of relative error with respect to the reference value). Consequently, the association of the pre-maneuver and post-maneuver orbits is still possible given a suitable association framework.

Conclusions
Two novel methodologies for the detection and estimation of impulsive maneuvers have been presented. The first, focuses on track-to-orbit associations and is intended to detect and estimate single burn maneuvers, by determining the maneuver that applied to the pre-maneuver orbit minimizes the post-maneuver observations residuals. The second, conceived for the orbit-to-orbit association problem, approximates maneuvers of two burns as linear perturbations over the nominal ballistic motion of the RSOs. They do not require a-priori information nor fine tuning and are able to provide a first estimation of the maneuver that can be improved as more post-maneuver tracks are available. Besides, the two methodologies can be directly applied to observations from other sensors, such as radar or passive ranging stations.
Results show a good performance even when dealing with few tracks over relatively short time periods, including low impulses ( ∼ 10 mm∕s ) and high impulses ( ∼ 10 m∕s ) in GEO. Therefore, the authors propose the use of these methods in operational cataloging chains to increase the flexibility and robustness of the maintenance of the catalogs of RSOs. The computational time, compiled in Table 14 indicate that their integration would not harm the computation load of the association framework. As expected, there is a strong dependency on the size and step of the maneuver time grid considered to estimate the maneuver. Even though this grid was not optimized for reducing the computational burden, the computation time is of the order of seconds for a single maneuver detection and estimation with the T2O methodology and tens of seconds for the O2O methodology.
Moreover, the methods could be easily integrated in association frameworks, such as [5], for the evaluation of hypotheses involving tracks and orbits with maneuvers. Only by doing so, the T2O and O2O association problems under the presence of maneuvers could be solved. The rationale behind this, illustrated in the results, is that a single track is not enough for the T2O methodology to reliably estimate a maneuver and thus, tracks must be associated before applying the T2O methodology. The detection and estimation of the maneuvers is performed using optimal control and simplified dynamical models, allowing to retain local minima for the solutions corresponding to different control efforts. In the case of a MHT framework, multiple local minima may translate into the expansion of the association tree, i.e.: generation of new hypotheses related to the maneuver, whereas in a hard-decision framework, this can be reduced to the optimal maneuver. In order to trim the association tree, or avoid taking a wrong decision, a maximum control effort can be defined such that |û| < u max . The value of u max determines the number of hypotheses that may arise from a maneuver event and also avoid considering unrealistic maneuvers. These methods have been conceived for the association problem, and therefore, their ultimate goal is not to provide the most accurate or realistic maneuver characterization, but one that allows the linkage between tracks and orbits (T2O) and orbits among themselves (O2O), particularly in data-scarce scenarios. Satellite operators may not always perform optimal maneuvers due to experience, safety or even working-hours scheme aspects and therefore optimal control metrics should not be blindly trusted. At the end, during cataloging operations the final goal is to ensure traceability between tracks and orbits so as to optimize SST network sensing data usage and detection of maneuvers, as well as to avoid duplicated objects. To reach this final goal, the two proposed methodologies should be integrated in an association framework. The authors are currently working on this integration and its application to real scenarios, which will provide end-to-end performance metrics, such as association metrics (true positives, false positives and false negatives) and total computational cost.  Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Data availability
The data underlying this article will be shared on reasonable request to the corresponding author.
Code availability The code developed and used for the preparation of this manuscript is commercial software property of GMV.