Satellite Manoeuvre Detection with Multistatic Radar

Traditional radar sensors used for surveillance rely on monostatic radar principles. However, recently the use of remote radio frequency telescopes as bistatic receivers represents an interesting way to reuse existing facilities while providing additional information to improve tracking accuracy. In this paper we study the benefits of using such a system for the task of manoeuvre detection in satellites in LEO and MEO. We investigate the conditions in which a multistatic radar is advantageous for this purpose, and show concrete results based on simulated data. Moreover, we propose novel manoeuvre detection methods, and compare their accuracy to methods found in the literature. A more general way of assessing the accuracy of these manoeuvre detection methods is also proposed, with the aim of taking into account that the parameters of the manoeuvre that actually takes place also have an effect on the accuracy. These can be split into optimal control based methods, and statistical methods. We found the addition of multistatic radar to allow considerable improvement in the accuracy of the manoeuvre detection process, an improvement that is shown to be greater the greater the baseline, i.e., the distance of the receiver to the transmitter. Furthermore, the manoeuvre detection methods that accurately model the uncertainty in the measurements were found to be the most accurate.


Introduction
Space Surveillance and Tracking (SST), a subset of Space Situational Awareness, "comprises technologies to detect, catalogue and predict the objects orbiting the Earth, and the derived applications" [1].In this work we focus on a particular aspect of SST, namely manoeuvre detection for satellites in LEO and MEO.Since their 36 Page 2 of 25 inception, ground based SST systems primarily have utilised radar sensors due to their ability to operate in very long ranges and under various atmospheric conditions while also providing very accurate range measurements.Radar systems track objects by providing information about their range, i.e. distance between target and radar antenna, and, with lower accuracy, also on their azimuth and elevation angles.Range measurements considered in this work come with accuracy in the order of 100 m, while the angle measurements come with errors in the order of a degree.Initially used for early missile warning, modern SST radars are designed to monitor targets in Low Earth Orbit (LEO) up to deep space.Having very high power transmissions, in order to improve their efficiency, radar can also operate in tandem with nearby radio telescopes, referred to as bistatic receivers, forming what is known as a multistatic configuration.In this configuration, the reflected signal is received not only by the primary emitter station, but also by distant radio-frequency (RF) telescopes.Such a system measures the bistatic range, i.e. the sum of distances between the target and both receiver and transmitter, with similar level of accuracy.
A prime example of such a multistatic system is the Tracking and Imaging Radar (TIRA) located at Fraunhofer Institute for High Frequency Physics and Radar Techniques (FHR), Germany and the Effelsberg radio telescope which when paired form a bistatic system with a baseline of about 20 km.This can improve the minimum detectable target size from 2 cm at 1000 km down to 1 cm due to the higher sensitivity of Effelsberg.While being technically bistatic, the TIRA transmitter and the Effelsberg radio telescope act more as a quasi-monostatic radar due to the very high target altitudes relative to the baseline.The experienced bistatic angles, BS , are generally small and the perceived radar cross section (RCS) of the target will be very similar to that of a monostatic radar.Recently, the use of long baseline bistatic radar systems for SST has been investigated within the NATO SET-293 RTG. 1 The difference with existing bistatic systems is that the RF telescopes used as receivers are remotely located from the radar emitter allowing larger bistatic angles and essentially viewing the target from different aspect angles.As examples, in this work we consider two long range bistatic pairs, in addition to the previously mentioned.One also has TIRA as the transmitter, but is paired with the radar telescope in Chibolton, England, with a baseline of around 600 km.The other is a trans-continental system where the transmitter is the Millstone Hill Steerable Antenna (MISA), located at the MIT Haystack Observatory, Massachusetts, USA, and the receiver is the Westerbork Synthesis Radio Telescope (WSRT) located in Westerbork, Netherlands.The resulting baseline is around 5600 km.If the baseline is too long, it becomes harder or even impossible to have a line of sight simultaneously from both transmitter and receiver, and for that reason the latter system is not suitable for LEO objects.One of the benefits is that the target might exhibit higher bistatic RCS for some angles thus improving detectability.In addition, the radar system could integrate multiple receiving signals, i.e from multiple transmitters or using multiple receivers, resulting in a higher signal to noise ratio (SNR) [2,3]. 1 https:// www.sto.nato.int/ Lists/ test1/ activ ityde tails.aspx?ID= 16,824.While preliminary analysis on GEO satellites has shown that spaceborne targets can be detected in these bistatic configurations, the processing framework required to make the most use of the data collected by such systems is still being developed.
The use of tracking data from radar systems allows the reconstruction of the motion of space objects beyond the simple orbit determination.If the dynamics of an object being tracked are not fully known, e.g.there is an unknown acceleration, the latter can be reconstructed.If this unknown acceleration is caused by on-board manoeuvring thrusters, its reconstruction can be helpful in recognising patterns or intentions.This process is called behaviour analysis.
In this work we investigate some manoeuvre detection methods in the literature and apply them to our use-case of LEO/MEO spacecraft being observed by bistatic or multistatic radar.We then investigate to what extent the increased data quality coming from bistatic radar observations improves the accuracy of the results of the behaviour analysis process when compared to the monostatic case.We also propose novel methods, and adaptations to existing methods, that show how accurately taking into account the uncertainty in the measurements tends to improve the quality of the results.In general, the metrics studied provide a quantity that is higher the more likely it is that a manoeuvre was performed, and in some cases also a mathematical description of a manoeuvre that explains the observations.
We then test these methods on realistic scenarios to assess their performance relative to each other.We consider in particular two possible scenarios: repositioning to modify its ground track, and an orbit re-positioning to shadow another satellite.For the latter, we consider the specific case of the Kosmos 2542 and 2543 satellites, which, flying in formation, approached the American KH-11 satellite in January 2020.
Because the performance of a metric depends on the manoeuvre that actually took place, a further contribution of this paper is a quantity, here dubbed "quality of metric", which aims to provide a more general description of the accuracy of a manoeuvre detection metric, taking into account how this varies with the actual manoeuvre performed and with the measurement noise.
There are a variety of manoeuvre detection methods in the literature.Some of them are based on Kalman filters.Classic filtering techniques for manoeuvre detection include the work of Castella [4].This filter assumes frequent measurements occurring during manoeuvres, so frequent that the dynamics model is a straight line prediction, i.e., the effect of gravity is ignored in the time frame between measurements.The manoeuvre detection, as is commonly the case, is achieved by analysing the measurement residuals.The method, however, does not explicitly estimate the accelerations, and the noise spectral density is increased based on a function with ad-hoc parameters.
Another well known filter based technique is dynamic model compensation, described in [5].The acceleration is estimated online, assuming exponential correlation between acceleration biases.Lubey and Scheeres [6,7] develop an optimal control based strategy, combining the Pontryagin maximum principle for trajectory optimization with Kalman filters.
While these are very useful for scenarios where the target is observed continuously throughout the manoeuvre, enabling an accurate full reconstruction and 36 Page 4 of 25 tracking during manoeuvre, they have the disadvantage of not being suitable for scenarios where the observations are separated by long periods of time, with the manoeuvre occurring fully within that time.This is because they require the linearisation of the dynamics, and the assumption of Gaussian distributed states.In this work, we consider those scenarios with long stretches of time between observations, which are characteristic of observations from a single multistatic radar setup of a target in LEO/MEO.
Other methods can be more robust to the non-linearity of the dynamics, making them more applicable to our use cases.The work of Serra et al. [8] is one example.It is based on finding the optimal control law linking two observations.Like Lubey and Scheeres [6,7], the optimal control optimises the integral of the squared acceleration, but the dynamics do not need to be linearised.We therefore compare these methods' accuracy for manoeuvre detection with that of optimising the exact deltav, modelled as a sequence of impulse manoeuvres, and modelled as a sequence of thrust and coast arcs [9].
We also test methods that are not based on finding optimal control laws.Just like the navigation filter based methods mentioned previously, the measurement residuals are used to estimate if a manoeuvre occurred, such as described in [10].The residuals are quantified by the Mahalanobis distance (MD) between the expected state, following a no manoeuvre trajectory, and the actual observation.This metric is based on the assumption of a Gaussian distribution over the state variables.When the state is propagated over a large period of time, however, the assumption of Gaussianity loses applicability, due to the nonlinearity of the dynamics, as previously mentioned.Therefore, we propose another approach, where the MD is minimised while allowing the initial and final states to vary, under the constraint that they are linked by a no manoeuvre trajectory.This removes the need for propagating the uncertainty, and so, if the measurement noise is Gaussian, no approximation is being made on the distribution.
We start by describing manoeuvre detection methods found in the literature, along with novel methods proposed by us, in Sect. 2. The use of multistatic radar data is discussed in Sect.3, focusing on investigating the conditions in which the addition of a bistatic receiver most improves the accuracy of manoeuvre detection.A method for assessing the accuracy of manoeuvre detection metrics, the aforementioned "quality of metric", is then discussed in Sect. 4. Finally, test cases based on realistic scenarios, such as a satellite changing its ground track or shadowing another satellite, are described in Sect. 5.The results of applying the methodology in previous sections are then shown.

Manoeuvre Detection
The process of manoeuvre detection consists in determining, based on observations of the state at multiple different times, if there was a manoeuvre in between.Knowledge of the state at one instant in time along with the dynamics allows predicting the state for all time, apart from deviations that can be due to manoeuvres but also sensor error and imperfect knowledge of the dynamics (often modelled as process noise).For this reason, a state observation never exactly corresponds to the expected value.To determine whether this mismatch is due to a manoeuvre or only the other two factors is the task of manoeuvre detection.All of the methods considered here produce a metric which is higher the more likely it is that a manoeuvre occurred.This is then used to obtain Receiver Operator Characteristic (ROC) curves in Sect. 5.
In this Section we describe two manoeuvre detection methods from the literature based on two different approaches, optimal control based and likelihood based.We then follow each with proposed variations or novel methods based on these approaches, intended to improve their accuracy in detecting manoeuvres.
To detect a manoeuvre, a minimum of two observations of the state are required, one before and one after the manoeuvre.This can be seen as a worst case scenario, as in which the least amount of information is available.If more observations are available they can be converted into such two observations by batch estimation or forwards/backwards filtering [11] separately for the observations before and after the suspected manoeuvre.A common example is using the theory of attributables [10], where tracks of observations are converted into single observations.In this work we consider, without loss of generality, only two observations y 0 and y f , of the states x 0 and x f , where the subscripts 0 and f refer, respectively, to observations before and after the putative manoeuvre.Both y 0 and y f are random variables due to the pres- ence of measurement noise.In the absence of noise, the observations would be given by the measurement function h as y = h(x) .The values of the state estimated from measurements y 0 and y f will be written as x0 and xf , and are themselves random variables.An estimate of x f obtained by propagating x0 will be termed xf = F(x 0 ) , where F is the state transition function which propagates the state for a no manoeuvre trajectory.

Optimal Control Based
One can safely assume that any trajectory performed by a spacecraft will be minimizing fuel consumption, as that is a critical resource for any mission.Therefore, a sensible approach is to solve the optimal control problem linking the observed states, and use as metric G, where u(t) is the control law represented as an acceleration over time, the set U, in the most general case, is the set of functions ℝ → ℝ 3 , x is the state as the Cartesian position and velocity, and J is the total delta-v required for control law u, The following subsections contain some examples of such methods. (1) 36 Page 6 of 25

Smooth Cost Function
In the work by Serra et al. [8], instead of minimizing G, they minimise G 2 given as The benefit is that this cost function can be easily minimised using the Pontryagin maximum principle.A disadvantage is that the resulting control law u * 2 will not be optimal with respect to the actual manoeuvre cost, which is measured by J(u) .Even though G 2 does not correspond to a delta-v value, the Cauchy-Schwarz inequality provides an upper bound to the value of J(u * 2 ) , given as Since u * 2 does not optimize J, we can write G ≤ J(u * 2 ) , which, combined with the inequality in Eq. ( 4), means that there is an upper bound on G obtained using G 2 , The solution of Eq. ( 3) will produce a continuous control law, as opposed to a sequence of impulse manoeuvres as defined in Sect.2.3, which are typically the optimal solution to transfer problems without constraints on the acceleration.Thus, while this method is useful for manoeuvre detection and quantification, it is not as useful for manoeuvre reconstruction.In the next sections we will look at other methods that do not have these disadvantages.
Since in [8], the authors consider a case where only angle measurements are available, they replace the constraint on the final position in Eq. ( 3) by where Ω is an admissible region, defined as the set of states that could produce the measurements observed.Because in our work we do have both range and angle measurements, affected by error, we define the Ω as confidence sets around the measurements y , where the expected observation vector for a given state x is given by the observation function h(x) .Therefore, Ω is defined using the Mahalanobis distance: where R is the covariance matrix of the observations, F (−1) is the inverse of the cumulative distribution function (CDF) of the chi-squared 2 k distribution, k is the number of individual observations in y , and l, the p-value, is the minimum ( observation likelihood, i.e. the minimum value of p(y|x) (the probability of obtain- ing measurement y given the target's state x ) for x to be in Ω(y) .The value of l was chosen to be 90% in this work.In this work, we define both the initial state and final state as being in a set Ω , as opposed to Serra et al. [8] who only uses an admissible region in the final state.We then obtain The metric G 2 allows one to explicitly account for the uncertainty in the measure- ments, which will be shown in Sect. 5 to improve the manoeuvre detection accuracy.

Sequence of Impulses
For high thrust, the optimal control law will tend to be composed of relatively short bursts of thrust separated by relatively long periods of coasting.Therefore, a common simplification is to describe the manoeuvre as a sequence of impulses Δv (i) , optimizing the sum of their magnitudes J Δv = ∑ i ‖Δv (i) ‖ , with the same constraints as in Eq. ( 1), resulting in metric G Δv , We optimise this using Matlab's ® function fmincon, with its implementation of the interior point algorithm [12], having as control parameters the Cartesian components of each Δv (i) as well as the time between consecutive manoeuvres Δt i , with i = 1, … , N , where N is the number of impulses considered.The optimization space U Δ can therefore be written as ℝ 3N × ℝ N + .In this work, we used N = 5 , however, the optimiser will often make some of these manoeuvres of zero magnitude, meaning a smaller number of manoeuvres is actually found.
The process of majorization-minimization [13] is also used to deal with the undefined gradient when Δv (i) = 0 .This consists in iteratively optimizing a function Q(u, û(k) ) , called the majoriser: for any û(k) , u , then it can be shown that û(k) will tend to a local optimum as k → ∞ .The same can be said for any function Q a that can be obtained by adding or multiplying a constant to a function Q satisfying the above [13].
G Δv = min {Δv (1) , Δt (1) , …}∈U Δ Our cost function has a similar formulation to the "total variation" in [13], so we obtain our Q following the same process, resulting in which is equivalent to optimising The iterative process in Eq. ( 10) is repeated 10 times, as it was found through testing that the improvement to the delta-v was not significant when using more iterations.An additional implementation detail is that the optimization variable passed to the objective function used by the optimisation algorithm has the delta-vs scaled, i.e., i) .This keeps all parameters in similar scales, which typically improves convergence, and helps with finite difference estimations, thresholds on the first order of optimality, conditioning of hessian matrices, etc.

Continuous Thrust
Alternatively, a manoeuvre can be approximated as being composed of thrust arcs with constant thrust pointing in a constant direction in the radial-transverse-normal frame, which is more suitable for low thrust engines.The tool FABLE [9] efficiently produces propellant cost estimations for these types of trajectories.The references [9,14,15] fully describe the transcription and the algorithm used.A difference in our case is that we do not consider a departure from an hyperbolic orbit around the Earth, so the parameters describing the hyperbolic excess velocity are not used.The cost function G FABLE (x 0 , x f ) is the propellant mass consumption given by FABLE.

Likelihood Based
Our current problem of manoeuvre detection can be viewed as one of null hypothesis testing.The null hypothesis H 0 is that no manoeuvre took place, i.e. x f = F(x 0 ) , and the deviation in the measurement is due to the uncertainty in the state estimation only.Therefore, a suitable metric is the negative log likelihood of the observations assuming the null hypothesis, which, making the common assumption of Gaussian noise and linearising the dynamics, results that the log likelihood grows monotonically with the Mahalanobis distance where Δx = xf − xf is the estimate of Δx , and P is the covariance matrix of this Δx , given by where P 0 and P f are the covariance matrices of x0 and xf , respectively, and is the state transition matrix [16].This procedure is described in [10], where it is referred to as "reachability analysis".Reachability analysis refers to any set of metrics that are based on a distance between the observed final states and the predicted states in the absence of a manoeuvre, and in fact includes the optimal control methods, where the control cost is the aforementioned distance.

Minimal Mahalanobis Distance
The calculation of G MD requires propagating the uncertainty to obtain P .This is done by linearising the dynamics, and assuming a Gaussian distribution for xf .One way to avoid such assumptions, as well as the need for propagating uncertainty altogether, is to use the minimum log likelihood under H 0 , where X is the set of states.In our case, X = ℝ 6 , as the states x are written in Car- tesian coordinates.The optimizing value of x 0 and its propagation F(x 0 ) correspond to the maximum likelihood estimate for those states if a manoeuvre is believed to not have occurred (null hypothesis H 0 ).This is also a batch estimation method.The optimal negative log likelihood value can be used as a manoeuvre estimation metric.This also avoids estimating P 0 and P f , which requires linearising the measurement function h.
We can now formulate the metric G MD by only assuming Gaussian measurement noise with covariance matrix R , which makes Eq. ( 16) equivalent to Like G 2 , the metric G MD explicitly considers the uncertainty in the measurements, and does not require assuming the states have a Gaussian distribution, but it does not require any parameters, whereas G 2 required the parameter l.This can be an advan- tage for G MD if no information is available that suggests a good choice for l, which requires some consideration.If l is too large, the Ω sets in Eq. ( 8) may be empty, whereas if l is too small G 2 may be zero when manoeuvres have actually occurred.On the other hand, having a parameter that can be tuned can also be advantageous.The parameter l allows the user to adjust how conservative this metric is.(16) min 36 Page 10 of 25

Manoeuvre Detection with Multistatic Radar
In monostatic radar, an electromagnetic wave or signal is emitted in a certain direction.If the signal intercepts a target, its reflection is received in the same location where it was transmitted from.We denote the range, azimuth and elevation as observed from the transmitter station as 1 , 1 , 1 , and their respective rates as ρ1 , ̇1 , β1 .Furthermore, 1 is the position vector from the transmitter to the target.
With bistatic radar, the receiver is located far away from the transmitter, at a distance L, the baseline.In this case, the bistatic range and its rate are observed, defined as BS ≜ 1 + 2 − L , with 1 and 2 being the distance from the transmitter to the target; and the distance from the target to the receiver, respectively.We also have the observation angles and their rates as observed from the receiver station, 2 , 2 , ̇2 and β2 .
One disadvantage of bistatic radar is that both stations must observe the target simultaneously.This means that the transmitter and receiver antenna beams must overlap in order to generate a composite antenna beam in which the target is visible by the bistatic system.This constraint reduces the number of observation opportunities.This can be ameliorated by the addition of multiple receivers, in a multistatic configuration.Such a configuration would also provide higher quality data when the target is in the line of sight of multiple receivers.
The accuracy of the manoeuvre detection process depends not only on the metric being used, but also on the quality of the measurements, and in particular for this discussion, on the geometry of the radar transmitter, receiver, and target configuration.
Suppose we have a vector of observations h(x) that allows determining the state per- fectly.The covariance matrix of the state observed in this way can be approximated as where R is the covariance matrix of our observations, which we consider to be a diagonal matrix with the variance of each measurement, and H is the Jacobian of the observation function h(x) .A criterion for a measurement being good at detecting a manoeuvre is having a small variance along the direction Δx , which is the deviation in the observed state caused by the manoeuvre.This occurs when Δx is aligned with the eigenvectors with the smallest eigenvalues (variances).These will be close to the direction of the measurement with the lowest variance, coincident if H is orthogonal as is the case for monostatic observations.In our scenarios, the measurements with the lowest variance will be the range and range rate for the position and velocity respectively.
The position and velocity vectors of the target in an inertial reference frame are r and v respectively.The gradient of BS with respect to r is while the gradient ∇ v ρ is defined similarly.
This provides an intuition for determining when will bistatic radar provide the greatest advantage to the results.We illustrate this with a scenario, described here, and with results in Sect.5.1.Consider a case where the transmitter is at latitude 30° N and the (18 receiver, in the same meridian, is initially in the same position.Suppose a target in an equatorial orbit then performs an in-plane manoeuvre, and is observed along its orbit, such that the Δr and Δv caused by this manoeuvre are in the orbital plane, i.e. the equa- torial plane.As the receiver station is brought closer to the equator, the gradient ∇ will have a larger projection onto the orbital plane, and thus the manoeuvre will have a larger effect on the measured range.This is illustrated for ∇ x  = ∇ v ρ in Fig. 1.In the latter scenario, ∇ x BS can have a larger in-plane component than ∇ x 1 does, which implies a lower variance for an in-plane direction.So we should expect that in these conditions, the bistatic radar measurements on their own will lead to better manoeuvre detection capabilities than monostatic, and indeed in Sect.5.1 this is shown to be the case.
On the other hand, if the baseline is small, the gradient of the bistatic range measurement ∇ BS will have similar direction to the gradient of the monostatic range ∇ 1 , and as such, will contribute mostly to lowering the variance in the direction of ∇ 1 .When the baseline is higher, ∇ BS becomes further from ∇ 1 , leading to lower vari- ances along directions perpendicular to ∇ 1 .Therefore, we can expect an increase in baseline to improve manoeuvre detection results, in the general case where Δx does not align too closely with ∇ 1 .
As an illustration, if is the angle between ∇ 1 and ∇ BS , and these observations have standard deviations of and 2 (because BS is a sum of ranges) respectively, then the eigenvalues of P as defined in Eq. ( 18), for the position vector in the plane containing the stations and target, will be Clearly, when = 0 , i.e. when the transmitter and second receiver are in the same location, one of the eigenvectors goes to infinity, as we have no information along that direction.From there, as cos( ) increases, the largest eigenvalue decreases significantly.It can be shown that if , which shows how an increase in baseline L leads to a decrease of the largest eigenvalue, which again represents the variance along a particular direction in the estimated state. ,

Quality of a Metric
Since a manoeuvre detection metric is essentially a binary classifier, a natural tool to assess their accuracy is a receiver operator characteristic (ROC) curve.To quantify the quality of a metric with a single value, the area under the ROC curve (AUC) can be used.Let G m be a random variable representing the value of a metric when a manoeu- vre occurs, and G 0 when it does not.The AUC equals P(G m > G 0 ) [17].
To approximate the AUC for analytical purposes, we assume that the metric in consideration can be approximated as a quadratic form: where A is the Hessian matrix of G, and x 0 and x f are the deviation between the observed state and a pair of reference states consistent with a no manoeuvre scenario.We will now write the vector of observed states as Further, we assume that x 0 and x f are observed with noise that follows a Gaussian distribution with zero mean and, as before, covariance matrices P 0 and P f respec- tively.From these assumptions, both G 0 and G m follow generalised chi-squared dis- tributions, with known mean and variance [18].
To make the analysis easier, consider the following change of coordinates: where Δx = x f − F( x 0 ) represents the deviation in the final state that is due to the manoeuvre, and is the state transition matrix [16].The measured Δx will have a covariance matrix given by P in Eq. ( 15).The Hessian matrix can now be transformed into these coordinates By taking the bottom six rows and rightmost six columns of matrix Â , we get ÂΔ , the component of the Hessian that depends only on Δx , leading to the following approximation of G: By making a Gaussian approximation to G m − G 0 , results that the AUC grows monotonically with the Mahalanobis distance between the real value of G(Δx) and the Gaussian approximation to the distribution of G m − G 0 , given by: The value of q{G} can then be used as a relative measure of the sensitivity of a met- ric to a particular manoeuvre.Its formula in Eq. ( 27) shows that the quality of a metric depends on the particular manoeuvre being applied (through Δx ), on when it is observed (through x ), on the observation uncertainty (through P ), and on the met- ric itself (through ÂΔ ).The dependency on P in particular contains in it the depend- ency on the radar system, as described in Sect.3. The addition of a bistatic receiver will reduce the eigenvalues of P 0 and P f , thus increasing the value of q{G} .For small Δx , the following approximation is valid, This takes a particularly simple expression for G MD , where it becomes where d is the number of dimensions in the state variable, 6 in our case, and ǦMD is the value of G MD calculated using the exact values of the state vectors.Note that it depends on Δx.
As previously stated, and evidenced by Eq. ( 27), which metric is best may depend on factors such as which manoeuvre actually takes place, represented by Δx .For example, the quantity q can be used to find the manoeuvres for which G 2 is better than G MD , by finding the values of x f for which q{G MD } − q{G 2 } is negative.These are the linear combinations of the eigenvectors associated with the negative eigenvalues of A q {G MD } − A q {G 2 } .In Sect.5.3, a scenario is set up with synthetic data.For an inclination change manoeuvre, the metric G MD performs better, but nonetheless a manoeuvre can be found such that q{G MD } − q{G 2 } is negative.Simulating this manoeuvre and sampling G MD and G 2 results in an ROC curve where G 2 is cleary producing better results than G MD .This illustrates the purpose of introducing this quantity q.It provides a slightly more general description of the quality of a metric at determining whether a manoeuvre occurred or not, by clarifying the dependency of its accuracy on the actual manoeuvre Δx. ( 36 Page 14 of 25

Results
For these tests, we assume that the observations are affected by Gaussian noise.The standard deviations for the angle measurements are given, similarly to [8], by assuming the values of , and are observed with n = 3 independent observations, sepa- rated by an interval of time of Δt = 20s .We assume these observations are affected by independent Gaussian noise with standard deviations σ = σ = 1 arcminute and σ = 100 m .This information is turned into a tracklet , , , α, β, ρ by performing a least squares fit to the data.From this process the standard deviations on the rates are also obtained.For the bistatic range, we have BS = 2 and  ρBS = 2 ρ .Note that with these values, when the range is above approximately 344 km, the variance in the state estimation along the range direction will be lower than in the orthogonal directions.
An orbit without and with a manoeuvre is simulated, and observations y 0 and y f are sampled with N = 100 samples for each case.These observations are simu- lated as having independent Gaussian noise with the standard deviations mentioned before.From these noisy observations, the states x 0 and x f are calculated, using maximum likelihood estimates when this problem is over-determined.

Comparison of Monostatic with Bistatic Radar: Test Case 1
Following the discussion in Sect.3, we present here the results where a spacecraft in equatorial orbit is observed as it crosses a meridian where two ground stations lie, the transmitter at 30° N, and the receiver on the equator.In this case, for simplicity, we assume that the initial state is known perfectly, and there is only uncertainty in the measurement of x f .Figure 2 shows the Receiver Operator Characteristic (ROC) curves for the G 2 metric, for monostatic, bistatic only, and the combination of the two.These curves show that using bistatic radar, in this configuration, leads to a significant improvement compared to monostatic, since the geometry of the bistatic radar in this case is such that it allows better discrimination along the direction in which the state is varying due to the manoeuvre, as explained in Sect.3.For instance, at a 10% false positive rate (FPR), monostatic radar measurements only allow a 60% true positive rate (TPR), while with bistatic radar, a TPR of 85% is achieved at the same FPR, increasing to 95% with multistatic.

Manoeuvre Detection Methods Comparison: Test Case 2
We now present a more realistic test case.The orbit is a circular orbit with 5000 km altitude and 80 degrees inclination.The transmitter is the Millstone Hill Steerable Antenna (MISA), located at the MIT Haystack Observatory, Massachusetts, USA, and the receiver is the Westerbork Synthesis Radio Telescope (WSRT) located in Westerbork, Netherlands.The measurements are separated by approximately 7 h, taking place approximately two orbits apart.The manoeuvre in question is an inclination change of 1°.This scenario is illustrated in Fig. 3.
The ROC curves are shown for the monostatic observations in Fig. 4 and for combined bistatic and monostatic observations in Fig. 5.As expected, G symb is unsuitable for this test case, since this manoeuvre is not an in-plane one.Figure 4 also shows that G Δv , being closer to measuring the actual cost of a manoeuvre, leads to better accuracy than using G 2 , although this advantage is less noticeable in Fig. 5.Even though G 2 and G FABLE are both simulating continuous trajectories, they are optimiz- ing different quantities, so it is interesting that their results are quite close.But in both cases, the best method for higher values of the true positive rate was G 2 , which can be understood as a consequence of it taking into account the uncertainty as well as applying a control distance metric.The addition of uncertainty in x 0 could explain why the G MD metric is no longer as good for these test cases, since the propagation of the state is a nonlinear process which makes the distribution of state deviation no longer Gaussian.In the definition of G 2 it is only assumed that the observations have Gaussian error, which is still valid in this case.
Figures 6 and 7 show histograms with the predicted delta-v cost of the manoeuvre in question, compared to its real value.The presence of noise in the measurements tends to lead to over-estimation of the delta-v.For the no-manoeuvre case, in blue in both Fig. 4 ROC curves for some manoeuvre detection methods, for monostatic observation data, test case 2 Fig. 5 ROC curves for some manoeuvre detection methods, for bistatic+monostatic observation data, test case 2. The curves are all very close to representing a perfect classifier figures, the correct value would be 0 m/s, while for the manoeuvre case in orange, the delta-v was just above 100 m/s.The results in these figures show that the addition of bistatic radar leads to much closer results not just in manoeuvre detection, but also in its reconstruction.Numerically, the root mean squared error for the case with manoeuvre goes from 73.8 to 24.7 m/s by the addition of bistatic observations, which shows the benefit of the additional accuracy provided by bistatic observations.

Quality of a Metric Results
We apply the theory in Sect. 4 to obtain, for test case 2 where typical manoeuvres produce a better result with G MD , a manoeuvre for which G 2 is the better, thus dem- onstrating the accuracy of that theory.For this demonstration, we use the same conditions as in test case 2, except there is only uncertainty in the second observation so as to align with the assumptions made in Sect.4, and the value of the standard deviations for the measurement errors were reduced by 20, so that the small x f approximation can be valid.The eigenval- ues of A q {G MD } − A q {G 2 } are The only negative eigenvalue is orders of magnitude smaller than the larger positive ones.This too supports the idea that G MD is the more robust metric of the two.How- ever, by choosing a x f aligned with the eigenvector corresponding to that negative eigenvalue, chosen such that √ (t f − t 0 )G 2 is 13.9 m/s, it is possible to find a manoeuvre for which G 2 is better, resulting in test case 2a, as shown in Fig. 8.That figure also includes the manoeuvre from test case 2, scaled in the same way, labeled test case 2b.

Cosmos 2542 Satellite Shadowing: Test Case 3
As a realistic test case, we consider the Russian satellites Cosmos 2542 and Cosmos 2543, which in January 2020 manoeuvred so as to shadow the American satellite KH-11, also known as USA 245 [19].By consulting the website "in-the-sky", 2 while the historical ephemerids are not available, it is possible to see graphs of the mean altitude, eccentricity, and inclination.The values around 19th January suggest an inclination change manoeuvre from around 97.9 to 97.6°, and a decrease in the apoapsis of about 58 km, costing respectively around 33 and 15.7m/s of delta-v if they were performed separately and by a high thrust engine.Using this information, we model this spacecraft performing these two manoeuvres, and being observed before and after the manoeuvres, by the TIRA system located at Fraunhofer Institute.
The results are in Fig. 9 for monostatic observations.Clearly, the statistical based methods G MD , G MD and G 2 are far superior, showing perfect results.To better com- pare these methods, and get a better overview of the effect of having additional stations, we repeat the experiments considering two receivers, the Effelsberg 100-m Radio Telescope and the Chilbolton Observatory, see Fig. 10.The baselines of the resulting bistatic pairs are of 20 and 600 km respectively.The noise on the angular measurements is also increased from one arcminute to 0.5°.Figures 11 and 12 show the results for G MD and G MD for different combinations of receivers.The metric G 2 is not included due to the very long convergence times when the noise level is this high.
Two conclusions can be reached from these results.First is that G MD is superior.We note that the MC estimate of the covariance matrix used to estimate G MD was obtained with 1000 sample points, which was a higher number of state propagation functions that the optimization in G MD ever required for this case.Secondly, the longer baseline system improves the results significantly more than the shorter one, as expected from the discussion in Sect.3.
In addition, for this test case, the eigenvalues of q{G MD } were found to always be greater than q{G 2 } , indicating that there is no maneouvre, in this observation conditions, for which G 2 is better at detecting manoeuvres than G MD .Further- more, the hessian of G MD was found to be identical to that of G MD with rela- tive difference in the order of 10 −8 , suggesting that they have the same Hessian, Fig. 10 Ground track and locations of transmitters and receivers for test case 3 Fig.11 ROC curves for test case 3 with angular observation noise with standard deviation of 0.5°, for metric G MD although a theoretical proof of this is outside the scope of this work.The fact that despite this, the metrics have very different performances, shows the limitations of this method of estimating the quality of a metric, limitations which we speculate are related to the small Δx approximations taken when obtaining it.

QoM Sensitivity Analysis
The QoM presented in Sect. 4 was developed by making several assumptions.In this section we test the validity of the Gaussianity approximations.There are two such approximations in use.Firstly, the likelihood functions of the measured xf and pre- dicted xf final states are assumed to follow a Gaussian distribution.After the metric function is approximated as a quadratic form (Eq. 21), statistical analysis reveals that its distribution is a generalized chi-squared distribution [18].Its mean and variance are then obtained following this theory, and the MD between the real value of the metric when the manoeuvre occurs and the Gaussian distribution with this mean and variance corresponds to q{G} .Therefore, the accuracy of this Gaussian distribu- tion at estimating the AUC is directly related to the accuracy of q{G}.
Therefore, in this section we compare estimates of the AUC obtained by three methods: • Numerically-the random variables G 0 and G m are sampled with 1000 sample points each, and the fraction of pairs for which G m > G 0 is the estimate of the AUC; • As a generalised chi-squared distribution-using the methods developed in [20], for which code was made available on the Matlab file exchange, 3 the parameters of the distribution of G m − G 0 are computed, and using the same library the AUC is computed as well; • As a Gaussian distribution with the same mean and variance as the chi-squared estimate.
Another assumption made in developing q was the second order Taylor approximation of G as a function of Δx , in Eq. ( 21).This section is not concerned with evalu- ating the validity of that approximation, so the tests are all conducted on the metric G MD , which can be written exactly in this form.Another detail is that instead of obtaining P as in Eq. ( 15), the unscented trans- form [21] is used to propagate the uncertainty.This removes the need to compute and reduces linearisation errors, but does not completely remove them.Finally, note that for G MD , ÂΔ = 2 P−1 .For test case 3, outlined in Sect.5.4, let Δx ref be the change in final state caused by the manoeuvre described therein.Let also ref be the measurement noise standard deviation defined for the same test case.Results have been obtained with the methods outlined above for a scenario with noise only on the final measurement, and another with noise on both measurements.For all cases, the AUC is estimated for various values of Δx and which are scalar multiples of Δx ref and ref .
The results in Fig. 13 were obtained with noise only in the measurement of x f .The only source of non-linearity is in the measurement function, as there is no noise being propagated.The chi-squared and Gaussian approximations have similar predictions.Therefore, the final Gaussian approximation used to define q is not the main source of error.
The shaded area indicates the 2 sigma (95%) confidence interval on the likelihood of the numerical estimates under the assumption that G m − G 0 follows the chi- squared distribution.For low noise values, it is found that the chi-squared and the numerical results are in agreement, but as the noise increases this stops being the case, visually indicated by the numerical estimates being outside this shaded area.As the noise increases, so does the mismatch between numerical and analytical estimates of the AUC, because the non-linearity effects become more pronounced.There are two factors contributing to this.When the manoeuvres are too small (large), the AUC is close to 0.5 (1), and large changes to Δx result in small changes to the AUC.The numerical estimates are too coarse to reveal these variations, so the only useful data for validating these methods corresponds to the region where the AUC is in between 0.5 and 1.The values of Δx within this region increase as increases.This means both the noise in the measurements and the magnitude of the manoeuvre, two factors that increase the effect of non-linearities, increase as for the area of the AUC plots that allows comparing the various methods.
Adding noise in the initial state measurement means the non-linear propagation has the effect of making the state distribution depart even more from a Gaussian.
Results obtained as before are in Fig. 14.The estimates of the AUC depart from numerical estimates for even lower amounts of noise.The theoretical derivation and the simulated results in this paper show the extent to which the use of bistatic radar receivers improves the accuracy of manoeuvre detection with radar measurements.The relationship between this accuracy and the geometry of the observations, in particular its improvement with increases in baseline (i.e., the distance between the transmitter and receiver(s)), and on the manoeuvre being detected, was also shown with theoretical arguments and synthetic data.
This was complemented with an investigation of some manoeuvre detection methods from the literature, in addition to methods we developed.These were then compared using simulated data.The results found that in general, an accurate modelling of the uncertainty is key to obtaining good results.
We also introduced a method for analytically assessing the quality of a metric for the purposes of manoeuvre detection, which depends on the response of the metric to deviations in the state, on the manoeuvre being performed, and on the measurement uncertainty.This quantity allows quickly evaluating how accurate a metric would be for all possible manoeuvres, avoiding the possibility that results obtained for a specific scenario are only valid for the specific manoeuvre that was tested.There are limitations.The first is that it requires that the metric is a function of class C 2 of the measurements.It also needs to be approximated by a quadratic form, which makes its applicability limited to small manoeuvres.A more conceptual limitation is that it does not provide a simple quantification of the effect of the conditions the observation takes place in, i.e., the time and location of the satellite when it is observed, which in turn depends on its orbit, etc. Overcoming these limitations is a consideration for future work.
Future work will also include the application of state-of-the-art uncertainty quantification methods to manoeuvre detection.In particular, we intend to consider filter based techniques that do not require the assumption of Gaussian distributions, such as Gaussian mixture model filters.

Fig. 1
Fig. 1 Configuration of receiver and transmitter, showing the gradients of the range measurements and the direction of the state deviation

Fig.
Fig. ROC curve for different observation conditions

Fig. 3
Fig. 3 Ground track for the test case 2, showing the locations of the transmitter, receiver, and measurements

Fig. 6 Fig. 7
Fig. 6 Histograms of delta-v control cost as estimated using G Δv , for the monostatic case

Fig. 12
Fig.12ROC curves for test case 3 with angular observation noise with standard deviation of 0.5°, for metric G MD

Fig. 13 Fig. 14 3 36
Fig.13 Comparison of estimates of the AUC for test case 2 with measurement noise only on x f