{'fraction_non_alphanumeric': 0.02878613417368644, 'fraction_numerical': 0.0075752984667595905, 'mean_word_length': 4.851507092198582, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Automatic maneuver detection and tracking of space objects in optical survey scenarios based on stochastic hybrid systems formulation

The state space representation of active resident space objects can be posed in the form of a stochastic hybrid system. Satellite maneuvers may be accounted for according to control cost or heuristical considerations, yet it is possible to jointly consider a combination of both. In this work, Sequential Monte Carlo filtering techniques are applied to the maneuvering target tracking problem in an optical survey scenario, where the maneuver control inputs are characterized in a Bayesian inference process. Due to the scarcity of data inherent to space surveillance and tracking, model switching probabilities are not estimated but derived from the ability of the state representation to fit incoming measurements. A Markov Chain Monte Carlo sampling scheme is used to explore the region assumed accessible to the object in terms of the hypothesized post-maneuver observation and a novel and efficient control distance metric. Results are obtained for a simulated optical survey scenario, and comparisons are drawn with respect to a moving horizon least-squares estimator. The proposed framework is proved to allow for a capable implementation of an automated online maneuver detection algorithm, thus contributing to the reduction of uncertainty in the state of active space objects.


Introduction
The interest in Space Surveillance and Tracking (SST) has been steadily growing in the last two decades, as the effects of Earth orbital congestion become more evident. Means to characterize the Earth's orbital population are then required in order to ensure a safe and orderly growth of space activities. Accordingly, sensor data retrieved by surveillance radars and optical telescopes is used to build up and maintain catalogs of Earth orbiting objects. SST cataloging systems are in charge of processing incoming measurements, thereby updating the state of cataloged objects or detecting new ones. These measurements are usually obtained as a time series of observations that stem from a common object, also known as tracks. If the distance between a well-established orbit of a known object and a given track is sufficiently small, then such set of observations is assigned (correlated ) to the object and its corresponding orbital state is updated. On the contrary, the track is categorized as uncorrelated and compared against other uncorrelated tracks. Different combinations of uncorrelated tracks are evaluated until a sufficiently high number (typically ∼ 3-4, see Hill et al. (2012)) is compliant with a specific orbit, resulting in a new object identification. Operational or active spacecraft perform frequent maneuvers, usually dictated by mission requirements or collision avoidance, potentially hindering the data association problem. In the absence of capable maneuver detection and estimation methods, post-maneuver measurements may trigger the identification of a new (duplicated) object as part of the catalog build-up and maintenance process. To overcome this limitation, the space accessible to an active spacecraft can be expressed in terms of a control effort (Holzinger (2011)), with a post-maneuver state probability distribution based on the control effort itself or otherwise previously characterized maneuvers. Uncertainty regarding the state of active spacecraft poses a threat to data association and tracking of space object catalogs. Incorporating knowledge regarding the maneuver history of this active population can aid in providing better predictions for their state, partially contributing to the reduction of epistemic uncertainty at a Space Traffic Management (STM) system level (see Delande et al. (2018), Hilton et al. (2019)).
Various efforts have been directed in this regard, initially to identify maneuvers performed by known objects. Kelecy et al. (2007) devised a method to detect maneuvers based on changes in Two Line Element (TLE) data: by comparing orbits reported at different epochs one could identify non-environmental perturbations, i.e. thrust arcs. The interest in maneuver detection then moved towards solving the data association problem. The former work by Holzinger et al. (2012), but also Lubey (2015), reduces the combined maneuver detection and data association problem to solving the optimal control input required to be compliant with a certain observation given a prior orbit. The aim is not to characterize the maneuver in terms of control input, but to provide a control distance metric that can be used for hypothesis testing, i.e. decide whether an incoming observation is triggered by a maneuvered object. A thorough implementation of this approach, with certain modifications, has been recently explored by Serra et al. (2021), whose main focus is on the data association problem at a SST system level considering multiple targets. Additional efforts have been directed towards jointly solving the data association and maneuver estimation, as proposed by Pastor et al. (2020). Therein, maneuvers are estimated in the form of impulsive burns performed within a discrete temporal grid, where the aim is to minimize a cost function weighting the distance to the observations and the required control effort in terms of ∆V magnitude. In parallel, an increasing interest was given to the so-called patterns of life (cf. Cox et al. (2016)). In general, spacecraft maneuver to maintain their location within certain orbital slots, so the use of previously characterized maneuvers may result in an increased performance in terms of maneuver detection and state estimation. Machine Learning (ML) techniques are well suited to this type of applications, and implementations can be found in the works of Shabarekh et al. (2016), Singh et al. (2016) and Abay et al. (2018). Not only ML techniques have been applied to this end, but also former statistical methodologies. Siminski et al. (2017) propose the use of a Kernel Density Estimator (KDE) to characterize maneuvers in terms of the relative change in orbital elements and the pre-maneuver orbital state. In this approach, data association is based on the information contained within the KDE and a control distance metric ∆V . An admissible region is defined based on a maximum threshold ∆V max so the generated post-maneuver state hypotheses are tested to lie within such region.
In this paper, the authors propose the use of a Stochastic Hybrid System (SHS) (cf. Hu et al. (2000)) formulation as it conforms a natural statistical framework for the modelling of systems for which there are multiple accessible dynamical models, as is the case for maneuvering targets. The flexibility inherent to this formulation allows fusing information coming from very different sources, e.g. knowledge derived from historical data and optimal control. State space filtering of SHS is usually solved by means of multiple model gaussian (cf. Jia et al. (2015), Goff et al. (2015)) and non-gaussian (cf. Andrieu et al. (2003)) algorithms. Due to the non-linearity of the space environment and the special characteristics of optical survey scenarios (McCabe and DeMars (2014)), a Sequential Monte Carlo (SMC) scheme is used to estimate or infer the state and active mode of the system conditioned on a sequence of observations. Our proposal takes the form of a bi-modal system, in which either the non-maneuvering or maneuvering mode is active prior to the last hypothesized observation. In the non-maneuver mode, deterministic dynamical models usually provide sufficiently accurate approximations. However, in the presence of maneuvers, the control input to the system is indirectly inferred through its effects on the state. To this end, Markov Chain Monte Carlo (MCMC) techniques are used to explore the posterior state distribution conditioned on different prior assumptions, effectively conforming alternative hypotheses. These hypotheses are then pruned, promoted and merged, ultimately converging to the true association sequence as more information is available. This procedure is illustrated in Fig. 1, which depicts the post-maneuver measurement association sequence in the state space, clearly distinguishing between the ballistic and maneuvered accessible regions. The accessible space in the event of a maneuver, termed admissible control region, defines the boundaries for MCMC sampling. The latter can be regarded as a revision of the admissible region (cf. Milani et al. (2004)) based on the control effort required to transfer from an initial to a final orbit. The requirements imposed by the definition of the admissible control region and MCMC sampling lead to the derivation of a novel and efficient control distance metric, P , which allows for a fast computation of the control cost separating two orbits within a specific time of flight. Under the assumption of optimal maneuvers in a control cost sense, one can elaborate a candidate post-maneuver state distribution function with higher probability densities in lower P regions. However, note that the orbital information contained within a few optical observations is limited, and so the optimal transfer assumption may not be sufficient to provide an estimate of the post-maneuver orbital state. The use of prior knowledge, i.e patterns of life, can improve the estimation performance given the target follows a repetitive maneuvering plan. Contributions of this work are thus not limited to a combined optimal-heuristic maneuver detection but also a realistic uncertainty characterization of the state of active space objects following a Bayesian inference process.
The content of the paper is structured as follows. The problem and its mathematical formulation is discussed in Section 2. Different methods for estimating the state of stochastic hybrid systems are discussed in Section 3. Special emphasis is given to non-linear non-gaussian filtering both in the presence and absence of a proper definition for the underlying dynamical model. Section 4 deals with the generation of maneuver hypotheses, i.e. post-maneuver state estimation based on alternative prior assumptions. In this regard, the space accessible to the system is bounded in terms of a maximum expected control effort, derived by means of a novel control distance metric. Hypotheses are then generated via MCMC sampling based on optimal control or conditioned on the patterns of life of the target. Results for simulated optical observations are presented and analyzed in Section 5, where the filter performance is compared against different maneuvering target tracking filter implementations. Finally, conclusions are drawn in Section 6 together with guidelines for future lines of research. to the object is given by the admissible control region.

Problem formulation
The problem of maneuver detection and data association can be regarded as a maneuvering target tracking problem if required to be solved in an automated manner, whose mathematical description can be naturally framed under the stochastic hybrid systems formulation. Herebelow, the general description of a SHS is first presented, and then particularized to the case of maneuvering space object tracking based on optical observations.

Stochastic hybrid systems
The general stochastic hybrid state space filtering problem is governed by the following system of stochastic difference equations (cf. Blom and Bloem (2004); Lygeros and Prandini (2010)): where x k ∈ R n represents the physical state of the system, or continuous valued variables being estimated, u k is the control input sequence from t k−1 to t k , and m k ∈ M refers to the discrete valued system mode active at time t k . The mapping function f (·) represents the underlying dynamical model, which typically follows a system of first order differential equations. g(·), usually referred to as diffusion coefficient, is included in order to account for any kind of colored noise for the Wiener process µ k , and is commonly approximated in a scalar or matrix form. Mode transitions are given by the mapping π(·), where ν k is an independent random variable representing the stochastic process, analogous to µ k in Eq. (1). The observations available to the system are indicated by z k ∈ R p , in which h(·) is the deterministic non-linear mapping between the system internal and observed states, and q(·) is again introduced to allow for any type of noise model realized through the random variable γ k .

Maneuvering space object tracking
For the intended application, the state of the system and observed quantities read being r k the Cartesian representation of the position vector,ṙ k its temporal derivative at time t k and B k the solar radiation pressure coefficient of the target object. The latter is included in order to account for variations in the area to mass ratio (exposed to the sun) that stem from changes in attitude (orientation) and mass expenditure due to maneuvers, among others. Note the focus is on optical observations, which are mainly used to track high altitude objects. Thereafter, no atmospheric effects are considered as the atmospheric density is assumed to have a negligible impact on the system dynamics.
Observation data from optical sensors is retrieved as a time series, i.e. track, of correlated right ascension α and declination δ pairs, which are expressed in a topocentric reference frame centered at the sensor location. These correlated pairs conform the so-called optical tracks, whose typical duration spans from 2 to 10 minutes depending on the survey strategy. The orbital period of the observed objects is usually of the order of days or tens of hours, so the orbital arc described within a track is relatively small (< 1%). Information contained within such short tracks is therefore limited so they are usually approximated by a linear regression at the mean epoch. The output of this linear regression is commonly referred to as attributable, and was first proposed in Milani et al. (2004). An attributable (α, δ,α,δ) can then be used to define the line-of-sight of an object and its temporal derivativė The line-of-sight and its temporal derivative are related to the state variables via where r s andṙ s are the position and velocity of the observing site, ρ is the range or positional distance between the object and the sensor, andρ its time derivative. With regard to the system dynamics, the mapping f (·) in Eq. (1) is generally approximated by a dynamical model that takes the form of a perturbed restricted two body problem, with the Earth as central body. Therein, the system is assumed to be a point particle upon which different forces and perturbations act. Besides the central gravity field exerted by the Earth, the effect of Sun, Moon and planets are also considered, together with a non-spherical Earth model and solar radiation pressure: recall atmospheric effects are disregarded due to the orbital regions of interest in optical surveys. Approximating reality in a purely deterministic fashion requires an extensive modelling effort, including but not limited to the roto-translational effects (c.f. Früh et al. (2013), Misra et al. (2016)). To partially mitigate the error introduced by mismodeled dynamics, it is common practice to introduce process noise, which in this case is represented by the second term in the right-hand side of Eq. (1).
The subset u k represents the control input to the system and is related to the active mode m k by Due to the scarcity of data inherent to optical survey scenarios, it is convenient to estimate the state after the maneuver x k instead of the control sequence u k , which involves inferring the active mode m k at time t k . Accordingly, the proposed filtering problem considers maneuver detection and tracking but not the estimation of the control input to the system.

State estimation of stochastic hybrid systems
Within this paper we propose a method to perform state space filtering for the SHS described in Section 2. Hereunder, we briefly introduce the concept of state space filtering and discuss on different filtering techniques commonly applied to state estimation of non-linear systems. Bi-modality of the system is considered at an upper filtering level, resulting in different schemes for estimating the state in the non-maneuvering and maneuvering modes. In the former case we propose a SMC method following the known or ballistic system dynamics. In the maneuver mode we generate samples from two different post-maneuver state distributions, conditioned on 1) the control distance from the post-maneuver state to the pre-maneuver orbit and 2) heuristics derived from the maneuver history, as detailed in Section 4. The latter procedure results in different maneuver hypotheses whose sampled distributions are determined from the application of an advanced MCMC technique, i.e. the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm.

State space filtering
State space filtering refers to the process of estimating the optimal sequence of states x 0:Nz according to a set of system observations z 1:Nz = z(T ) at times T = {t 1 , t 2 , ..., t Nz }. Numerous techniques can be applied to the aforementioned problem attending to the mathematical behavior of the system of interest. Under linear dynamics and Gaussian unbiased noise assumptions, an analytical solution exists, as is the well-known Kalman Filter. This scheme has been successfully applied to non-linear systems under appropriate transformations, see the work in Julier and Uhlmann (1997) and Arasaratnam and Haykin (2009). Based on the assumption of Markovian dynamics, the Kalman filter adopts a sequential scheme so that the estimation is assumed to be optimal at each observation and the sequence of state estimates is not refined as new observations arrive. Another class of estimation methods, termed Batch Least-Squares (BLS), are aimed at minimizing the squared distance to the observations. Under the assumption of a Gaussian state probability distribution, the estimation is performed on a batch of measurements. These methods do not target to solve the state estimation problem in a sequential manner but are rather aimed at fitting a dynamical model to a set of observations. Note however, under certain modifications they can be readily applied to sequential estimation, as suggested in Robertson et al. (1996) for the moving horizon estimator. Whenever the probability distribution function (pdf) of the state is not known a priori (or it is subject to non-linear transformations), the complexity of the estimation problem significantly increases. This is to say, if the Gaussian assumption x(t) ∼ p(x(t), t) ≈ N x(t);x(t), Σ(t) , wherê x and Σ are the expected state estimate and co-variance, is dropped, then the evolution of the probability distribution function p(x(t), t) with time is governed by the Fokker-Planck-Kolmogorov equation, Eq. (10), which makes reference to the work by Fokker (1914) and Planck (1917).
The term f (·) in Eq. (10) follows the definition given in Eq.
(1), whereas w(t) corresponds to the second term of the right hand side of Eq. (1). This complicated partial differential equation is often approximated using Monte Carlo integration (cf. Metropolis (1987)), e.g. via SMC or MCMC techniques. Disregarding the assumption of a Gaussian state distribution leads to a re-definition of how the measurement information is introduced into the system. Kalman and Batch Least-Squares filters use the observed quantities to reduce the state co-variance and update the expected value (or state estimate). However, if the state distribution is allowed to take any realization, the update needs to be performed on the entire distribution.

Sequential Monte Carlo
SMC methods are referred to as Particle Filters (PFs) when applied to filtering problems. The term particle stems from their approach to managing non-linear transformations of statistical distributions, i.e. approximating Eq. (10). Given an initial distribution, a statistically significant number of samples N are randomly drawn, approximating the initial pdf as a sampled distribution where ω i is the weight associated to sample x i and δ(·) is the Dirac delta function. Each of these samples, or particles, is then propagated following the system dynamics so the non-linear transformation can be fully characterized at least in the statistical region of interest. This procedure results in an approximation of the a posteriori distribution of the state that converges to the true solution as the number of particles is increased, i.e. N → ∞. Conceptually, PFs are appealing since at a higher computational cost, it is possible to perform an arbitrarily accurate non-linear propagation of the state uncertainty.

Bootstrap particle filter
The combination of Sequential Importance Sampling and Resampling (SISR) conforms the Bootstrap Particle Filter (BPF), which is in fact the simplest realization of a capable SMC method for state estimation. Sequential importance sampling refers to the process of updating the sampled state distribution according to the information conveyed in the observation sequence. This update is reflected on the weight associated to each particle, so that at some point, a given particle may have negligible weight. The latter problem can be mitigated through the use of resampling. Sampling from an already sampled distribution results in a population with a lower weight variability, eliminating particles with smaller weights and duplicating those with higher probability. Algorithm 1 summarizes the main steps required by the BPF scheme, which is particularly simple and intuitive from a statistical perspective. These type of filters usually require a significant number of particles N M C to operate and, as such, they suffer from the curse of dimensionality (see Bengtsson et al. (2008)): the sample size is required to grow exponentially with the number of state variables to avoid collapse, which is expected to occur as log(N M C )/n → 0. Applications to high-dimensional systems also feature faster particle depletion rates, so naive resampling techniques as the one in Algorithm 1 may lead to a population that is concentrated on a single particle for relatively low dynamical noise levels.

Algorithm 1: Bootstrap Particle Filter
Apply the measurement update to obtain p(x k |z 1:k ): Resampling 3) Compute the Effective Sample Size: and update the particles according to:

Regularized particle filter
To deal with the particle degeneracy problem, there exist multiple resampling procedures of diverse complexity. One of these is given by the regularized particle filter, in which the resampling process no longer consists in duplicating particles. The discrete particle population is regularized to adopt a continuous form. Each individual particle is assigned a kernel, usually a Gaussian function, with a given bandwidth h. The population then takes the form of a KDE being K(·) a non-negative window function, e.g. a Gaussian. Musso et al. (2001) discuss on the implementation of regularized particle filters and the design of the kernel estimator. The authors propose to use Gaussian kernels and set the bandwidth according to the sample co-variance of the particle population and the dimension of the state space (cf. Silverman (2018)). This results in an adaptive tuning of the bandwidth, allowing for a more robust resampling step.

Markov Chain Monte Carlo
For some applications, there is not a clear definition of the underlying dynamical model, i.e the function f (·) in Eq. (1) cannot be properly characterized. In fact, this is the approach followed for the maneuvering mode of the stochastic hybrid system defined in Eq. (1), as stated in Section 2.2. In these cases, an alternative procedure is given by Markov Chain Monte Carlo (MCMC) methods, which focus on the exploration of a target probability distribution function, disregarding the underlying physical processes. The aim of MCMC algorithms is to sample from the posterior distribution rather than actually solving the Bayesian inference problem. According to Bayes' rule (cf. Bayes (1763)), the following assumption holds meaning that the probability of the state x conditioned on the observation z is directly proportional to the likelihood of a given state realization p(z k |x k ) multiplied by the probability of the state realization itself p(x k |z 1:k−1 ). Note that in this case the problem is not related to solving Eq. (10) since, in principle, we are not capable of approximating p(x k |z 1:k−1 ) following physical process assumptions, i.e. f (·) is not defined. Instead, one can elaborate an arbitrary prior, or proposal, distribution p(x k |z 1:k−1 ) from which to generate samples that are then retained or discarded according to their likelihood. MCMC simulation consists in generating a sequence, or chain, of samples from such proposal distribution. As the length of the chain increases, the sampled distribution converges to the posterior distribution p(x k |z 1:k ).

DiffeRential Evolution Adaptive Metropolis (DREAM)
The Metropolis-Hastings (MH) algorithm, first proposed in Hastings (1970), is one of the simplest and most used Markov Chain Monte Carlo methods. Given a prior jumping distribution π(x |x) and an initial point x 0 , candidate samples x k+1 are subsequently drawn from the jumping distribution according to The likelihood ratio is then compared to a random number u ∼ U(0, 1) so that if α ≥ u the candidate x k+1 = x k+1 is accepted, otherwise x k+1 = x k . As the number of samples increases, their distribution converges to the desired posterior p(x|y). This is due to the fact that samples are drawn from a jumping distribution, simulating a random walk whose movements are dictated by the ratio of likelihoods. The only tunable parameter within this method is the design of the jumping distribution, which in the most naive approximation may be a Gaussian centered on x k featuring a user-defined covariance matrix. There exist multiple alternatives for the definition of the jumping distribution in the MH scheme, tailored to different types of inference problems and target distributions. Among these, one can find the DREAM algorithm, developed by Vrugt et al. (2009). DREAM is a multi-chain algorithm, meaning that multiple x k are updated in parallel to improve convergence and efficiency, especially when exploring multi-modal distributions. Moreover, the parameters of the jumping distribution are dynamically adapted to avoid sampling from outside of the posterior support. Successful applications of the DREAM algorithm can be found in the works by Minasny et al. (2011) and Lu et al. (2017), where Bayesian inference and parameter uncertainty analyses are applied to complex environmental problems.

Maneuver hypothesis generation
This section discusses the generation of maneuver hypotheses, which is analogous to characterizing the posterior state distribution p(x k |z k , m k = 1) when the maneuvering mode is active in the interval t k−1:k . Exploration of this distribution is performed by means of MCMC techniques, for which the prior is uniformly distributed and bounded by the admissible control region developed in Section 4.2. Based on the belief of fuel optimal transfers, an approximation for the posterior can be written as where κ is a constant parameter and P (x k ) is the control distance metric defined in Section 4.1. The parameter κ controls the relative importance of the control cost compared to the measurement likelihood p(z k |x k ), so it can be set a priori, e.g. based on the control distance corresponding to the centroid of the admissible control region P (x ) defined in Section 4.2.1. Another option is to use previously characterized maneuvers (patterns of life) to approximate the posterior state distribution in the event of a maneuver, i.e.
M(x k , t k ) in Eq. 17 is a KDE that contains statistics of the maneuver sequence prior to t k , and κ h is a constant used to control the relative importance of both contributions. These statistics are derived from the relative variation in mean orbital elements implied by the maneuver, as well as the pre-maneuver orbital parameters, as discussed in Section 4.3.
Maneuver hypotheses are generated whenever the likelihood of an incoming track p(z k |x k ) ≤ p th falls below certain threshold, and they are processed in parallel by means of a regularized particle filter. Whenever the Maximum Likelihood Estimate (MLE) corresponding to certain maneuver hypothesis features a sufficiently low measurement likelihood, such hypothesis is pruned and so its assigned particles are eliminated from the sampled distribution representing the state of the target. At some point, the surviving maneuver hypothesis is required to replace the ballistic hypothesis so the filter can detect future maneuvers autonomously. To this end, we keep track of the measurement likelihood of each individual maneuver hypothesis r up to time t k as where φ ≈ 0.95 is a fading memory factor used to favour more recent tracks. The active ballistic hypothesis at t k is then indicated by the maximum L(r), providing a prior orbit to be used in the generation of subsequent maneuver hypotheses.

Control distance metric
We aim to develop an inexpensive metric to characterize the admissible control region as defined in Section 4.2, which can also be used in MCMC methods in the form of a log-likelihood. Former definitions of a control distance metric applied to maneuvering space objects can be found in the works by Holzinger et al. (2012) and Lubey (2015), which require to determine the solution to an optimal control problem. Therein, the authors emphasize on the scarcity of data for the particular problem of maneuver detection in optical survey scenarios. In the presence of long re-observation times and hidden (unobserved) states, it is convenient to derive a computationally efficient control distance metric that can be used e.g. in statistical sampling methods without significantly increasing the computational cost. Through a careful simplification of the dynamical model, it is possible to derive metrics that provide a sufficiently close estimation of the control effort required to acquire the post-maneuver orbit from the pre-maneuver one. To this end, we propose to approximate maneuvers as instantaneous velocity changes, or impulsive burns. Let the position and velocity vectors, expressed in Modified Equinoctial Elements, be governed by the set of first order differential equations given in the work by Walker et al. (1985), Eq. (9). Hereunder, we will refer to such system of equations as where the linear dependency of the state derivatives on the perturbing accelerations a p has been made explicit for convenience. Under the assumptions that 1) the control input is an impulsive burn, i.e. a p ∼ δ(t − t M ), and 2) the sensitivity matrix of the state with respect to the control input, A, remains approximately constant between the pre-and post-maneuver orbits, the relative change in orbital elements stemming from an instantaneous velocity variation can be approximated by the linear system ∆oe ≈ A(oe 0 )∆V.
The reader is referred to the former work by Escribano et al. (2021) for the complete mathematical derivation. Transfers between two general orbits require at least two impulsive burns, since single-burn maneuvers lead to a final trajectory that intersects the prior one. At this point, an additional dynamical approximation is introduced: the only accelerating perturbations acting on the dynamical system of interest are those due to the control effort. This assumption of Keplerian motion implies that the spatial geometry of the orbit is invariant in the absence of maneuvers, such that only the location of the object within such orbit (the true longitude L) varies in time. The effect of a sequence of n burns in the time-invariant orbital elements oe i can then be approximated as where the dependency of the sensitivity matrix A on time is implicit through the true longitude at the maneuver epoch L j . Note, however, the intention is to solve the inverse problem, i.e.: given an initial (oe i 0 , L 0 ) and final state (oe i f , L f ), determine the associated control cost. This boundary value problem is posed in the form of an optimization one, since it is of interest to determine the optimum maneuver sequence in a control cost sense. To simplify the problem, only two impulsive burns are considered, whose joint contribution to the change in orbital elements may be estimated as Note the couple of burns (∆V 1 , ∆V 2 ) can be readily solved in a least-squares sense. Let the target change in orbital elements be ∆oe i t , and the maneuver sequence ∆V tot = (∆V 1 , ∆V 2 ), then we can define a cost function being c 1 a cost index used to express the relative importance of the control cost with respect to the injection error. Note the subscripts t and e are used to indicate the target and estimated changes in orbital elements, respectively. The optimal maneuver sequence ∆V * tot that minimizes the cost function can then be obtained by setting the partial derivative of J with respect to ∆V tot to zero, thereby leading to in which A is the augmented sensitivity matrix of Eq. (22), i.e.
Despite the solution given by Eq. (24) simply requires a matrix inversion, the following constrained non-linear optimization problem needs to be solved to determine the optimum pair of true longitudes (L 1 , L 2 ) Compared to the control distance metric proposed by Holzinger et al. (2012), the size of the Non-Linear Programming (NLP) problem is dramatically reduced in this case. A practical approach to solve the optimization in Eq. (26) may consider the use of a gradient descent method, which can be applied to an appropriate set of initial conditions to improve convergence to a global minimum.
Finally, the definition of the control distance metric P reads where || · || is used to indicate the Euclidean norm of a vector, and (∆V * 1 , ∆V * 2 ) are the first and second burns that stem from evaluating Eq. (24) at the optimum set of true longitudes (L * 1 , L * 2 ). The metric P has therefore units of velocity, as it is the result of solving the fuel optimal two-burn transfer between two Keplerian orbits under certain linear approximations.

Admissible control region
In order to apply MCMC methods to estimate the state transition in maneuvering intervals, it is convenient to define proper bounds for a prior distribution of the state. In the following, we take advantage of the work by DeMars et al. (2012). who propose a revisit of the admissible region developed by Milani et al. (2004) tailored to the needs of Resident Space Object (RSO) cataloging. Therein, bounds for the expected range and range-rate values (ρ,ρ) are developed based on different orbital regimes. This is, given an observation in the attributable format (α, δ,α,δ), the position and velocity of a target object is completely defined by the position and velocity of the observing site (r s ,ṙ s ) and the range and range-rate (ρ,ρ) according to the expressions in Eqs. (7-8).
It is then possible to define bounds on the range and range-rate given admissible sets for the semi-major axis a and eccentricity e. This method is applicable to the association of two observations, where there is not a clear definition of the orbital state of the observed object. In the case of maneuver detection, the association problem is defined between a clearly established orbit (of a cataloged object) and a hypothetical post-maneuver observation. Direct application of the constrained admissible region as defined by DeMars et al. (2012) can be used to reduce the search space for post-maneuver state estimation, but does not consider the underlying physical process undergone by the object.
Hereunder, we embrace an alternative definition, termed the admissible control region, in which the space accessible to the object is not bound in terms of semi-major axis and eccentricity, but based on the control effort required to reach certain range and range-rate values. A similar approach has already been explored by Serra et al. (2021), who propose a convexification of the admissible region based on a control-related energy metric. Still, bounds for the accessible space are given in terms of the maximum variation in this energy metric, largely based on the former control distance metric by Lubey (2015). These bounds are then translated into maximum and minimum semi-major axis and eccentricity values, thus disregarding the relative geometry of the initial and final orbits in terms of angular distance. This approach is beneficial from the association standpoint as there is no need to solve a NLP problem to accept or reject the maneuver hypothesis in a preliminary step.
The interest of the present work is not to tackle the maneuver detection and data association considering multiple-targets as in Serra et al. (2021), but to emphasize on the single-target case. Thereafter, we propose a target-based admissible control region in which an admissible set for the post-maneuver range and range-rate is determined in terms of a maximum expected control effort P adm . Attending to the topology of the control distance metric derived in Section 4.1, upper and lower bounds for the range and range-rate are elaborated.

Centroid of the admissible region
In contrast to the former definition of the admissible region, which yields a symmetric set in the range-rate space, our proposal requires a starting point. This point can be thought of as the centroid of the admissible set. An approximation to the latter may be given by the state compliant with an attributable that is closest to the pre-maneuver orbit x b in terms of the metric P , and is defined as The determination of this point involves solving an optimization problem over the range and range-rate in order to find the global minimum of the control distance metric. Note, in general, such point may not provide an accurate estimate of the true post-maneuver state since the actual purpose of a maneuver is to acquire a final orbit; not an optical track. Thereafter, and in the interest of computational efficiency, we propose to use an approximation to this global minimum x ≈ x opt solely based on geometrical considerations. It is, in essence, a valid initial guess to be used in the former optimization problem. Figure 2 illustrates the definition of x , which can also be expressed as being i and θ the inclination and true anomaly of the orbit, respectively. This centroid is a lowdimensional transformation of the pre-maneuver orbit, involving the minimum number of orbital parameters required to match the observation. Fig. 3 compares the post-maneuver estimation accuracy for the optimal state in terms of the control distance metric x opt and the proposed centroid x for two different maneuver types. Albeit small, the fuel optimal estimate shows an increase in accuracy when estimating east-west maneuvers since they usually imply changes in semi-major axis and/or eccentricity. These are not accounted for in the definition of the centroid, and so errors in range and range rate are mostly due to variations in those orbital parameters. On the contrary, the estimate given by the centroid for north-south maneuvers is significantly closer to the true post-maneuver state than that of x opt . The latter supports the fact that a fuel optimal transfer from an initial orbit to a single optical observation does not necessarily provide a good estimate. In this case, since out-of-plane maneuvers require higher impulses than in-plane ones, x opt overestimates the relative change in semi-major axis and eccentricity, with the aim of minimizing the required inclination change to match the observation. Based on this analysis, approximating the centroid of the admissible region as x seems a good trade off between computational efficiency and accuracy, and it is indeed compliant with an orbit preserving assumption, in the sense that one only expects phasing (∆θ) and inclination (∆i) change maneuvers.

Admissible region topology
An analysis of the behavior of the control distance metric in the ρ-ρ space is required in order to define proper bounds for the admissible control region. In general, one can define these bounds  in the form of an upper energy level, as suggested in Eq. (28). Therein, the maximum admissible control effort P adm may be defined in absolute and relative terms as P adm = min(P max , max(P min , k P P (x ))).
This allows to reduce the search space when observations are close to the pre-maneuver orbit, through k P and P min , but also when P (x ) approaches P max . More complex thresholding functions can be elaborated, yet the simplistic approach given by Eq. (31) provides reasonable bounds at a modest tuning effort.  Figure 4 depicts the behavior of the control distance metric P in the range and range-rate space for two different maneuvers, centered at x . It can be seen that x provides a good approximation to the centroid of the admissible region in both cases, and it is found to be close to the optimum in terms of control distance defined in Eq. (29). The topology of the control metric P exhibits a non-smooth behavior between local valleys, potentially hindering the application of gradient-based methods for determining x * . This behavior is thought to be related with the assumptions made in the definition of the control metric P , especially the one in Eq. (20).
The centroid x is uniquely defined in the unobservable range and range-rate space by (ρ ,ρ ), since the observed magnitudes are expected to match the measurements. In order to approximate the complicated outer boundaries of the admissible control region, depicted in Fig. 4, one can define certain search directions in the unobservable space along which to determine the intersections with P adm . For the two cases shown, a polytopic approximation may provide accurate results given a proper selection of search directions. Accordingly, the authors propose to approximate the admissible control region as an orthotope with search lines parallel to the range and range-rate axes and passing through the centroid x , in the form of an axis-aligned minimum bounding box. Moreover, the orthotope that numerically represents the admissible control region is extended in the observable space by 3-σ bounds along each dimension to accommodate measurement uncertainty.
An example of such polytopic approximation can be consulted in Fig. 5, where the admissible control region is depicted in topocentric spherical coordinates for an east-west maneuver. Therein, two-dimensional maps of the control metric are shown for all the possible combinations of states, centered at x . The distribution of P within the α-δ admissible space is nearly isotropic, showing control metric values close to the minimum, thus indicating that the boundaries implied by the measurement uncertainty are much lower that those given by the maximum expected control effort. On the contrary, the behavior of P in theα-δ plane feature higher gradients that eventually result in a over-estimation of the admissible control region if the expected measurement noise is used (determined according to Maruskin et al. (2009), Eq. 8). The majority of projections show a single global minimum with iso-P lines that are either concentric or parallel to one of the axes. That is not the case for the ρ-α plane, showing multiple local minima and some correlation between the range and the right ascension rate. The assumption x opt ≈ x seems to hold or at least x appears to be closer to x opt than the true post-maneuver state. Approximating C(x) as an orthotope following an axis-aligned box may result in an overestimation of the space accessible in terms of P , as indicated by the white regions outside the dotted line boundaries of Fig. 5. Nonetheless, it provides an extremely efficient computation of the expected region wherein to search for the post-maneuver state, which is desirable since C(x) must be determined individually for any post-maneuver observation. In fact, a more precise definition of the boundaries would not yield any increase in estimation accuracy since every candidate point within the admissible region is weighted according to its associated control distance.

Heuristic characterization
It is common practice to approach the maneuver detection and estimation problem using the maneuver history or patterns of life, see for instance Siminski et al. (2017);Shabarekh et al. (2016); Abay et al. (2018). These works propose the use of heuristic methods, based on a statistical characterization of the maneuver history and ML techniques. In general, spacecraft maneuver to be kept inside a specific orbital slot, usually referred to as station-keeping. Dynamical perturbations acting on the spacecraft motion are typically well-characterized, so that the various station-keeping maneuvering modes are already devised at the mission design stage. These usually feature certain temporal frequency and similar control magnitudes, with the aim of simplifying the operational workload. Thereafter, application of heuristics and ML methods seems suitable to this scenario, and has seen successful implementations in the previously cited works.
The proposal herein suggested is largely based on the work by Siminski et al. (2017), in which maneuvers are characterized in terms of relative variations in certain orbital elements ξ = [∆a ∆e ∆i] T . As discussed in Section 3 state space filtering is performed sequentially by means of a regularized particle filter. Post-maneuver observations trigger the generation of maneuver hypotheses, which are tested for correlation with subsequent tracks. At some point, there is a single hypothesis or group of equivalent hypotheses that survive, and then the filter is assumed to have converged in terms of state (and mode) estimation accuracy, i.e. there is a change in the active ballistic hypothesis indicated by L(r). It is not until such convergence is detected that the regularized particle filter is run in reverse mode (Lindsten and Schön (2013)), and the last (earliest) observation that associates with the surviving hypothesis is deemed the first post-maneuver track. Immediate pre-and post-maneuver orbits are then compared in terms of the relative variation in orbital elements. Note that statistics are readily obtained from this comparison since the preand post-maneuver states are given in the form of sampled distributions. Detected maneuvers are thus represented by a KDE M(x k , t k ) that is continuously updated with everyξ k = E[ξ k ] and being n M the number of detected maneuvers up to time t k , and n ξ the dimensionality of the feature vector ξ. In this case we have adopted a multivariate normal kernel estimator due to its simplicity, but there are other alternatives based on e.g. automatic bandwidth selection (see Wand et al. (1994)). The techniques discussed in Section 3.3 are used to explore the post-maneuver state distribution conditioned on both the incoming observations and the maneuver history, i.e. sample from p (x k |z k , M(x k , t k )). In doing this, samples generated from the admissible control region C(x) are evaluated in terms of [∆a ∆e ∆i], and this relative change is compared against the KDE containing the information related to previously characterized maneuvers. Thereafter, changes in orbital elements that are compliant with previous maneuvers are favoured, potentially leading to a more precise post-maneuver state recovery in the presence of a repetitive maneuver plan.

Performance metrics
Two different metrics have been devised to characterize the response of the algorithm in terms of accuracy and consistency. While the former is directly affected by the signal to noise ratio, a statistically consistent (unbiased) filter should be able to adjust its considered uncertainty in order to deliver unbiased estimates.
The objective of the proposed filtering scheme is to jointly solve the maneuver detection and tracking problem. Thereafter, emphasis should be in the post-maneuver state estimation and so the root-mean-square error with respect to the true post-maneuver state x + is given as a function of the number of tracks elapsed after each maneuver, n T , as defined in Eq. (33). Therein, m is the active mode introduced in Eq. (9) and 1 1 = [1 0 ... 0] is the indicator function with dimension n T .
Besides, with the aim of providing a means to analyze the uncertainty realism, or statistical consistency of the filter, the Posterior Cramér Rao Bound (PCRB) is used as proposed in Tichavsky et al. (1998). The PCRB is an analogy of the Cramér Rao Lower Bound (CRLB) (see Rao (1992)) when applied to the estimation process of random parameters. The PCRB can be regarded as a lower bound for the variance of any unbiased estimator, and is given by the inverse of Fisher's information matrix A recursive derivation of the Fisher information matrix, J , for the particular case of Gaussian process and measurement noise g(·)µ k , q(·)γ k with co-variance Q k and R k , is given by Eq. (35) for the non-linear discrete stochastic process described in Eqs. (1) and (3). The expectation operators in Eqs. (37-39) are approximated by means of Monte Carlo averages since the dynamical and measurement models are, in general, non-linear.
Note this derivation is only valid for the particular case of known active mode m k = ∅. In order to simplify the approach followed to compute the PCRB, 1) the active mode between subsequent observations will be considered known, and 2) the dynamical process noise Q k when m k = 1 will be assumed equal to the co-variance of the distribution p (x k |z k , M(x k , T )), with M(x k , T ) a KDE representing the entire set of true maneuvers characterized as described in Section 4.3.
Orbital uncertainty evolution in Cartesian co-ordinates presents a highly non-linear behavior, unlike parametric representations such as classical orbital elements or Modified Equinoctial Elements (MEE), among others (cf. Woodburn and Coppola (2014)). Thereafter, uncertainty is to be characterized in MEE and compared to the aforementioned PCRB to test the statistical efficiency of the proposed filter.
In this regard, we can construct a distance metric analogous to the Mahalanobis distance (see Mahalanobis (1936)) wherex k is the reference or ground truth value for the state estimate and C = J k , Σ −1 k may be the PCRB or the co-variance matrix of the estimates. The above distance d 2 ∼ χ 2 (n x ) follows a chi-square distribution provided the random vector x is normally distributed. An unbiased and statistically efficient method would comply with the latter for C = J k , whereas a consistent uncertainty characterization would feature d 2 ∼ χ 2 (n x ) for C = Σ −1 k . If the distribution of the estimation error distance d 2 features a lower skewness than χ 2 (n x ) then the uncertainty is overestimated and the filter can be considered pessimistic. On the contrary, a higher skewness is an indicative of an underestimated uncertainty (optimistic filter) and possibly a biased estimation.

Benchmarking methods
Various filter implementations have been evaluated in order to analyze the performance and improvements of the proposal over a standard operational approach. The differences between these methods are mainly related to the filtering scheme, maneuver detection and post-maneuver state estimation, so environmental modelling and measurement association remain the same.
With regard to the former, the dynamical model used for state estimation have the following characteristics: • Non-spherical Earth of degree and order 10. Measurement association is dictated by the Mahalanobis distance from the state to the observation, which is assumed to follow a χ 2 distribution with 4 degrees of freedom. The no-maneuver association threshold, p th is defined as the 3-σ gate of the aforementioned distribution for the current study. Post-maneuver association is subject to a threshold set on the maximum expected control effort, so that if the expected value of the post-maneuver state is farther in terms of P , i.e. P ≥ P max , then the track remains uncorrelated. This maximum expected control effort is highly related with the definition of the admissible control region C(x), whose outer boundaries are set attending to the following thresholds: P max = 10 m/s, P min = 1 m/s and k P = 3. The baseline method, assumed an operational standard, consists in a moving horizon estimator (hereafter MHE) that considers up to 6 subsequent tracks. Measurement association and maneuver detection is based on the aforementioned thresholds, also considering the admissible region defined using the developed control distance metric. Post-maneuver state estimation is conditioned on the sequence of observations after the detected maneuver (no-association), and the loss function is augmented to include the control distance P to the previous orbit, so that the post maneuver estimate corresponds to x opt . As a benchmarking option, the same moving horizon estimator is used, but this time considering the true maneuver sequence in spite of detecting a maneuver in terms of the Mahalanobis distance between the track and the ballistic trajectory. The latter method is termed MHE II, and is assumed to deliver the best estimate in a Bayesian sense when no heuristic or a priori maneuver information is included.
The proposed method consists in a regularized particle filter with a variable population size, wherein multiple maneuver hypotheses may be active. Each of these maneuver hypotheses is associated a fixed number of particles, N H = 1000, which are individually tested for association. A former implementation, termed SHF, discards the use (and automatic generation) of heuristics, therefore solely relying on the control distance metric to approximate maneuvers. The complete proposal is realized through SHF II, where heuristics are automatically generated based on detected maneuvers (in a feedback loop), in principle using all the available information in the estimation process through an alternative sampling conditioned on previous maneuvers.

Synthetic measurements
A simulation is carried out for a Geostationary Earth Orbit (GEO) spacecraft performing stationkeeping maneuvers, and the test scenario is defined as follows: • Subject: GEO Spacecraft equipped with chemical propulsion. Its assigned orbital slot comprises a mean longitude band = −4.8 ± 0.2 • and an inclination band i = 2 ± 0.05 • . The object is simulated for a total duration of 401 days (03/09 -05/10) using a dynamical model including the following perturbations: -Non-spherical Earth of degree and order 70. A total number of 20 detectable maneuvers, consisting of 29 detectable burns, are distributed along 401 days and 362 tracks. 5 of them are single-burn NSSK, while the total number of (doubleburn) EWSK is 16. Note that for these maneuvers (and burns) to be detectable an observation needs to be obtained between them, otherwise two or even three burns may be collapsed into one distinguishable orbit update. Table 1 gathers the results for two alternative methods: the moving horizon estimator (MHE) and the proposed combined optimal-heuristic hybrid filter (SHF II). The latter has shown more accurate in detecting maneuvers (and burns) at the correct epoch, being able to properly identify all detectable maneuvers and an additional burn. The MHE filter is also capable of identifying all the detectable maneuvers but four of them are detected with a delay of one observation. Nonetheless, the higher sensitivity of SHF II triggered two false detections. Both are found to be caused by an overestimation of the eccentricity, which leads to a bad characterization of the mean longitude for an EWSK maneuver and a lower estimated inclination in the case of a NSSK maneuver. On the contrary, since the MHE is capable of performing significant orbit updates, no false detections have been reported.
The latter comes at the expense of a poor state estimation performance, as indicated by the sequential estimation error for the maneuver sequence shown in Fig. 6. Therein, MHE, MHE II and even the (uninformed) SHF present difficulties at recovering the orbit after a maneuver. It is not only that the state estimation error at the first post-maneuver observation is relatively high (of the order of 10 km and 1 m/s), but also two to three additional observations are usually required in order to keep the error within the variance bounds indicated by the PCRB. The proposed method, SHF II, still presents some estimation errors greater than the PCRB but is capable of rapidly converging to the true orbit as observations arrive. This faster convergence is attributed to the lower magnitude of the initial post-maneuver state error, in the order of 1 km and 0.5 m/s.

Maneuvers (Burns) Correct Delayed False
MHE 16 (16) 4 (4) 0 (0) SHF II 20 (21) 0 (0) 2 (2)  Fig. 7 depicts the distribution of the d 2 distance as defined in Eq. (40) for MHE II and SHF II implementations. The conclusions that can be extracted from the pdf of d 2 are twofold: whether the filter is optimal in a statistical sense and how consistent is the estimated uncertainty with respect to the estimation error. The former is dictated by the PCRB, and corresponds to the blue histogram in   the figure. In both cases, the distribution appears to be positively skewed with respect to the target pdf χ 2 (n x ) and so either filter implementation is deemed sub-optimal. Nonetheless, SHF II presents a region of high density near the peak of the theoretical distribution, indicating a higher level of accuracy. Filter consistency is indicated by the orange bars, and can be inferred from the d 2 distance pdf when the estimated filter co-variance Σ k is used. The moving horizon estimator presents a higher skewness, so that the co-variance estimate is even lower than the PCRB. Thereafter, MHE II provides optimistic estimates as information regarding maneuvers is not included in the estimation process: the filter determines the optimal state based on minimizing the observation residuals and the co-variance is determined based on such residuals. The proposed filter, on the contrary, produces consistent state estimates since its estimated co-variance is in line with the estimation error. Regarding the expected estimation performance of the different filters, Fig. 8 shows the Root Mean Square Error (RMSE) as a function of the number of tracks elapsed after every maneuver. Note the estimation error for early post-maneuver tracks is of the order of 10 km in some cases, potentially leading to a wrong track correlation and suggesting an extension of the method to multiple maneuvering targets. The scenario is similar to what was already inferred from Fig. 6: MHE and MHE II show a higher initial error and slower convergence than the proposed method. A possible explanation may be the presence of multiple local minima; or even a global minimum that does not necessarily correspond to the true final orbit. As the number of post-maneuver tracks increases, these local minima collapse to the solution. The proposed filter partially mitigates this problem by 1) considering multiple hypotheses, and 2) introducing prior information to explore solutions that are neither global nor local minima. The former, as shown for SHF, seems to provide a more accurate state characterization from two to five post-maneuver tracks, despite the initial error being similar to that of both MHE implementations. Including heuristics has shown to be beneficial in approximating post maneuver orbits as regions that show no relevant features a priori, are explored as indicated by the patterns of life of the target object. These patterns of life are summarized in Fig. 9. Therein, the actual (detectable) maneuvers performed by the object are compared to the maneuvers identified and characterized by the proposed method. It can be seen that every detectable maneuver is placed within a co-variance (blue) ellipse, which coincides with the number of correctly detected maneuvers reported in Table 1. Note that not all maneuvers are characterized with similar confidence levels, as indicated by the size of the 3-σ co-variance bounds. In particular, there are two NSSK characterizations that show a significant standard deviation in the eccentricity when compared to the rest. One of these is directly related with a false maneuver detection, in particular the one on the top-right corner of the − vs i − plot shown in Fig. 9. This behavior is attributed to the lower control cost related to variations in eccentricity: orbits with higher eccentricity and lower inclination are compliant with the same tracks than higher inclination, more circular ones (if not observed at the antinodes) under certain observability conditions. Regarding EWSK maneuvers, the characterized changes in eccentricity seem to be consistent across all maneuvers, but that is not the case for the semi-major axis. The estimated ∆a values, especially those ∆a > 0 present a significant standard deviation, of the order of 1.3 km. Due to the positioning of the optical telescopes and the special characteristics of the GEO region, information regarding the semi-major axis is mostly conveyed in the time stamp of the observations. It is not until there is a sufficient temporal separation between tracks that the uncertainty in semi-major axis, and also eccentricity, can be reduced to acceptable levels. Accordingly, with the aim of improving the characterization of detected maneuvers, it is recommended to further reduce the uncertainty of the state estimates through a careful smoothing recursion.
In terms of computational cost, the proposed method incurs in a threefold increase with respect to the more efficient moving horizon estimator. This is partially due to the use of a regularized particle filter, which requires propagating a considerable number of state realizations ∼ 10 4 . A trade-off between uncertainty estimation accuracy and computational time may drive the use of, e.g. Gaussian Mixture Filters as in Li et al. (2016). Nonetheless, the intention of the current work is to provide a baseline for the statistical characterization of the state of a maneuvering object and so mode clustering methods as those required in Gaussian Mixture Filtering are dropped in favour of a more general statistical approximation, still at a higher computational cost. Figure 9: Maneuver characterization in terms of relative change in orbital elements (left) and pre-maneuver inclination and mean longitude (right). Blue ellipses are 3-σ co-variance contours for the characterized maneuvers, while black dots correspond to true maneuvers.

Summary and conclusions
A novel approach for the maneuver detection and tracking of space objects has been presented, relying on a stochastic hybrid systems formulation. Due to the scarcity of data inherent to optical space survey scenarios, estimation of the control input in the maneuvering mode is overseen in favour of a post-maneuver state estimation. The definition of an admissible control region based on a novel and efficient control distance metric has proved to be helpful in characterizing the set of feasible maneuvers performed by a target object between two subsequent tracks. Maneuver hypotheses are elaborated based on the control cost and the maneuver history, hence resulting in a combined optimal-heuristic approach. State space filtering is solved via a Regularized Particle Filter implementation, enhanced by a Markov Chain Monte Carlo characterization of post-maneuver state hypotheses.
Results are obtained for synthetic data, and comparisons are drawn against a moving horizon estimator. The proposed framework shows promising performances in post-maneuver tracking and maneuver detection accuracy, being suitable for online maneuver detection and data association purposes. The latter seem to be partially aided by the combined use of optimal control and maneuver heuristics, which helps in 1) identifying non-observable maneuvers, and 2) limiting the uncertainty in the presence of high control cost orbit updates. Moreover, success in the automation of maneuver detection and post-maneuver state estimation yields a proper characterization of the patterns of life, ultimately resulting in an increased predictability of the state of the population of active space objects.
The test scenario is representative in an operational context in terms of measurement uncertainty and environmental modeling, and so it is expected to perform similarly in an operational environment. A trade-off between computational cost and estimation accuracy may drive the use of a less demanding state space filtering technique, e.g. Gaussian Mixture Filtering, yet a careful implementation of the proposed method can enable online tracking of several targets. Certain modifications are required to adapt the methodology to multiple maneuvering target tracking, which could be based on a combination of state of the art data association methods and the proposed admissible control region. A thorough study on the different strategies that can be used to generate and apply heuristics may improve the maneuver detection and post-maneuver state estimation performances, especially if they are set according to known operator decisions or are tailored to specific procedures, e.g. limiting the accessible space to assigned orbital slots. Moreover, under a proper characterization of the different maneuvering modes, the bi-modality of the system can be dropped in favor of a multi-modal approach, ensuring traceability in the detected maneuver type.