Orbit determination from one position vector and a very short arc of optical observations

In this paper we address the problem of computing a preliminary orbit of a celestial body from one topocentric position vector and a very short arc (VSA) of optical observations. Using the conservation laws of the two-body dynamics, we write the problem as a system of 8 polynomial equations in 6 unknowns. We prove that this system is generically consistent, namely it admits solutions at least in the complex field. From this system we derive a univariate polynomial $\mathfrak{v}$ of degree 8 in the unknown topocentric distance at the mean epoch of the VSA. Through Gr\"obner bases theory, we show that the degree of $\mathfrak{v}$ is minimum among the degrees of all the univariate polynomials solving this problem. The proposed method is relevant for different purposes, e.g. the computation of a preliminary orbit of an Earth satellite with radar and optical observations, the detection of maneuvres of an Earth satellite, and the recovery of asteroids which are lost due to a planetary close encounter. We also show some numerical tests in the case of asteroids undergoing a close encounter with the Earth.


Introduction
The problem of computing the orbit of celestial bodies has attracted the interest of scientists since a long time, see Laplace (1780), Lagrange (1783), Gauss (1809).The recent improvements in the observation technology have posed new interesting mathematical problems in this field.The number of asteroids observations performed by modern telescopes is very large.Usually they can be grouped into very short arcs (VSAs) of optical observations, however it is not easy to determine whether VSAs collected in different nights belong to the same observed objects.In general, from a VSA we can not compute a reliable preliminary orbit, but we can try to put together different VSAs to perform this task, see e.g.Milani et al. (2005).Assuming that two VSAs belong to the same asteroid, we can write polynomial equations to compute a preliminary orbit using the conservation laws of the two-body dynamics.Recently, different ways to combine these integrals of motion have been developed and tested, see Taff and Hall (1977); Gronchi et al. (2010Gronchi et al. ( , 2011Gronchi et al. ( , 2015Gronchi et al. ( , 2017)).
In this paper we address the problem of computing the orbit of a celestial body (hereafter, the OD problem) from one topocentric position P 1 = (α 1 , δ 1 , ρ 1 ) at epoch t 1 and a VSA of optical observations, from which we can derive an attributable A 2 = (α 2 , δ 2 , α2 , δ2 ) at the mean epoch t 2 of the arc.Here, α, δ, ρ denote respectively right ascension, declination and topocentric distance of the body, while α, δ stand for the angular rates.
Using the conservation of angular momentum, energy and Laplace-Lenz vector we write the OD problem as a system of polynomial equations in the unknowns ρ1 , α1 , δ1 , ρ 2 , ρ2 , z 2 , where z 2 is an auxiliary variable, and the other variables allow us to obtain an orbit in spherical coordinates at the two epochs.
We prove that this polynomial system is consistent, namely it generically admits solutions, at least in the complex field, and we obtain a univariate polynomial v of degree 8 in the unknown range ρ 2 to solve the OD problem.Through a computer algebra software we are also able to show that the degree of v is minimum among the degrees of all the univariate polynomials in ρ 2 solving this problem.
The proposed method is relevant for different purposes, such as the computation of a preliminary orbit of an Earth satellite with radar and optical observations, the detection of maneuvres of an Earth satellite, the recovery of asteroids which are lost due to a planetary encounter, if the latter can be modeled as an instantaneous change of direction of the velocity vector like in Öpik theory ( Öpik, 1976).Here we show the results of the application of our algorithm to the orbits of some near-Earth asteroids, whose positions have been changed to increase the effect of the close encounter with the Earth.
The paper is organized as follows: after recalling the expressions of the Keplerian integrals in Section 2, we introduce the OD problem in Section 3 as an overdetermined polynomial system, whose consistency is shown in Section 4, where the univariate polynomial v is derived.The minimality of the degree of v is proved in Section 5.In Section 6 we discuss the selection of the solutions.Finally, in Section 7 we present the numerical tests.

Keplerian integrals
Let us consider a celestial body whose dynamics can be modeled by Kepler's problem where r and µ denote respectively its position and the gravitational parameter.Equation (1) admits the first integrals corresponding to the angular momentum, the energy, and the Laplace-Lenz vector of the body, respectively.We call Keplerian integrals these constants of motion.Note that we can write

The OD problem
Let us introduce a reference frame whose origin is the center of force of the Keplerian dynamics.We consider the Keplerian motion of a celestial body around the origin observed from a moving point of view.Assume the position q and the velocity q of the observer are known functions of time.The position of the observed body is given by where ρ represents the topocentric distance, and êρ = (cos δ cos α, cos δ sin α, sin δ) is the line of sight.
Let us assume that we know the position vector , and a very short arc of optical observations of the same object, from which we compute the attributable at the mean epoch t2 .We can combine these data to compute the quantities ρ1 , ξ 1 , ζ 1 , ρ 2 , ρ2 , where which are missing to have a complete set of orbital elements at both epochs.For this purpose, we use the conservation of the integrals listed in (2).Hereafter, we shall use subscripts 1, 2 for all the quantities relative to epochs t 1 , t2 .The dependence of position and velocity on the unknowns is given by where ê⊥ 2 = α2 cos δ 2 êα 2 + δ2 êδ 2 is a known vector.Note that r 1 = q 1 + ρ 1 êρ 1 is known.To write a polynomial system, we replace the term µ |r 2 | with an auxiliary variable z 2 , independent from ρ 2 , and set Lemma 1.The following relations hold: Proof.From the definitions of L 1 and L 2 we have Relations (5) immediately follow from the definitions of E 1 and E 2 .

Consistency of the equations
We will show the following property: Theorem 1.The overdetermined polynomial system (6) is generically consistent, i.e. it always has solutions in the complex field for a generic choice of the data P 1 , A 2 , q 1 , q1 , q 2 , q2 .
The proof makes use of the results presented in the following subsections.We start by noting that, because of relations (5), equation of 7 equations in 6 unknowns.Therefore, to prove Theorem 1 we show that system (7) is generically consistent.Moreover, for a generical choice of the data P 1 , A 2 , q 1 , q1 , q 2 , q2 , system (7) is equivalent to q 1 = q 2 = q 3 = q 4 = q 5 = q 6 = q 7 = 0, where and

The angular momentum
The angular momenta c 1 and c 2 , written in terms of the unknowns, become where D 1 , D 2 are defined as in (9), and Therefore, we can write with

Elimination of variables
Using relations (11), ( 12), the polynomials q 1 , q 2 , q 3 defined in (8) can be written as where In particular, the occurence of the variables ρ1 , ξ 1 , ζ 1 , ρ2 in q 1 , q 2 ,q 3 is at most linear.From the definitions above, we obtain Therefore, the coefficients Q (1) 100 and Q (1) 010 are both vanishing iff In fact, the second condition in the alternative above is equivalent to Thus, for a generic choice of the data, the angular momentum equations allow us to eliminate ρ1 , ρ2 and one variable between ξ 1 and ζ 1 .Assuming Q (1) 100 ̸ = 0, we choose to eliminate ξ 1 , which can be written as Substituting ( 13) in equations q 3 = q 2 = 0 and assuming W 12 ̸ = 0 we find ρ1 = − 1 Hence, by using relations ( 13), ( 14), and (15) we can eliminate the variables ξ 1 , ρ1 , ρ2 in the generators q 4 , q 5 , q 6 , q 7 .The resulting polynomials, named q4 , q5 , q6 , q7 , can be written as

A redundant relation
Here, we prove that equation q 4 = 0 can be dropped from system (7) without losing any solution.In particular, we show: the polynomial q 4 is generated by q 1 , q 2 , q 3 , q 6 .
Proof.Let us consider the polynomials q4 , q6 obtained from q 4 , q 6 by eliminating ρ1 , ρ2 , ξ 1 through relations q 1 = q 2 = q 3 = 0. We note that q4 = for some polynomials a j , b j .We prove that q6 divides q4 ; in particular, For this purpose, we first note that system q 1 = q 2 = q 3 = 0 is equivalent to c 1 = c 2 and, inserting the last relation in the expressions of q 4 , q 6 , after eliminating ρ1 , ρ2 , ξ 1 , we obtain where c 1 , c 2 , r 2 , ṙ1 , ṙ2 are meant as functions of ζ 1 , ρ 2 only.More details about these computations can be found in Appendix A.
We have

The latter relation immediately yields (16).
Setting we can write which concludes the proof.

The univariate polynomial
We compute a univariate polynomial u(ρ 2 ) of degree 8, which is a consequence of the equations in (7).Let us consider the polynomial system q 1 = q 2 = q 3 = q 5 = q 6 = q 7 = 0 (17) including all the polynomials in (8) except q 4 .To solve the OD problem, we solve the polynomial system defined by ( 17).Now, we compute the resultant of q5 and p 6 with respect to ζ 1 , which is denoted by v(ρ 2 ).For this purpose, the terms in q5 and p 6 are grouped in the following way: where .
Therefore, we have which has degree 8.

An optimal property
We show that the univariate polynomial v(ρ 2 ) has the minimal degree among all the univariate polynomials in the variable ρ 2 that are algebraic consequences of More precisely, let us introduce the polynomial ideal We will show the following result: Theorem 2. For a generic choice of the data we can find a Gröbner basis of the ideal I for the lexicographic order Proof.As a consequence of Proposition 1 we obtain I = ⟨q 1 , q 2 , q 3 , q 5 , q 6 , q 7 ⟩, where we dropped q 4 from the set of generators.Let us consider the elimination ideal We compute a Gröbner basis G ′ of the ideal I ′ for the lexicographic order ζ 1 ≻ ρ 2 using the software Mathematica1 : in this basis we have a univariate polynomial for some constant coefficients u i .The roots of u are the only values of the ρ 2 components of the solutions of system (17).We note that, since the resultant v is a univariate polynomial of degree 8, v is necessarily proportional to u, i.e. v = κu, κ ∈ R\{0}.
In place of the generators q 1 , q 2 , q 3 , can consider where the polynomials h 1 , h 2 , h 3 are defined by relations ( 13), ( 14), ( 15).The set is a Gröbner basis of the ideal I.In fact, the leading term of every polynomial in I is divisible by the leading term of one polynomial in G.
As a consequence of the previous theorem we obtain the following property: Corollary 1.The polynomial v has the minimum degree among the univariate polynomials in ρ 2 belonging to I.

Selecting the solutions
From the solutions of ( 7), we first discard the ones with non-real components and the ones with ρ 2 , z 2 ≤ 0. From the remaining solutions we can compute Keplerian orbital elements at epochs where c is the speed of light and ρ 2 is the value of ρ 2 for the j-th solution.We emphasize that the accepted solutions of the system are such that the Keplerian elements at the two epochs are the same.In fact, from c 1 − c 2 = 0 we get where i, Ω, a and e denote the inclination, ascending node, semi-major axis and eccentricity.In particular, from z 2 2 |r 2 | 2 − µ 2 = 0, which is a consequence of system (7), we get z 2 = ± µ |r 2 | .Having discarded the solutions with z 2 ≤ 0, we obtain which implies e 1 = e 2 , ω 1 = ω 2 .
Since the eccentrities e 1 , e 2 are equal, then also We observe that for the other orbit determination methods using the Keplerian integrals (Gronchi et al., 2010(Gronchi et al., , 2011(Gronchi et al., , 2015) ) the trajectories of the solutions at the two epochs are not necessarily the same.
In the following subsections we explain how to choose or discard the remaining solutions.

Selection without covariance
In case of multiple solutions, labeled with (j), we make our choice according to the following procedure.We propagate each of the computed orbits, referring to epoch t (j) 2 , backward to epoch t1 in the framework of Kepler's dynamics, and consider the norm of the differences ρ (j) 1 − P 1 , where is the topocentric position of the j-th solution propagated to epoch t1 .Note that in general we have |ρ We select the solution attaining the minimum value min j |ρ (j) 1 − P 1 |.

Selection using covariance
Assume that the data D = (P 1 , A 2 ), where be a solution of Note that, generically, the solutions S of system (18) correspond to the S components of the solutions of system (17). Let car , E (2) car ) and att , E (2) att ) be the vectors of the Cartesian coordinates and the attributable elements at epochs t1 and t2 .2Let us introduce the transformation T car att : E att → E car by ( 3), ( 4) for both epochs, and consider the map Ψ defined by Φ = Ψ • T car att .Equation Φ = 0 is equivalent to the system We introduce the vector V 1 = ( α1 , δ1 , ρ1 ).
The covariance matrix of the Cartesian coordinates at epoch t1 is with ∂E (1) car
For a given vector u ∈ R 3 we define the map Then, using ûT = −û, we have and Equations ( 7) do not set any constraints on the mean anomalies ℓ 1 , ℓ 2 of the pairs of Keplerian orbits that we compute.To properly select a solution we integrate back each orbit E (2) car computed at epoch t2 , to the epoch t1 of the orbit E (1) car , together with its covariance matrix Γ (2) .Then, we compute a predicted position vector P 1,p and its 3 × 3 covariance matrix Γ P 1,p .To define a metric to choose/discard solutions we use a three-dimensional version of the attribution algorithm (Milani and Gronchi, 2010), and introduce the identification penalty where We keep only the solutions with the lowest values of χ 3 or, if we wish to admit multiple solutions, we keep the ones whose value of χ 3 is below a certain threshold.

Numerical tests
In this section we show an application of our algorithm to asteroids undergoing a close approch with the Earth.We consider all the near-Earth asteroids present in the NEODyS database3 to the date of December 6, 2023, and select the 1305 asteroids with Earth MOID4 d min smaller than 10 −3 au.For each selected asteroid we compute the points P A and P ⊕ , lying on the osculating orbits of the asteroid and the Earth at epoch t 0 , where the minimum distance d min is attained.Then, we find the closest epoch t 1 > t 0 when the Earth arrives at P ⊕ and change the asteroid phase so that it would arrive at P A at the same epoch t 1 in the framework of the Sun-asteroid two-body dynamics.In this way, we try to enhance the effect of the close approach with the Earth at t 1 .A full n-body propagation, from t 0 to t 1 , is used to compute the topocentric position vector P 1 from the point of view of the Pan-STARRS 1 telescope (et al., 2019).Continuing the propagation from t 1 to a successive epoch t 2 , we derive an attributable A 2 by coordinate change.As a result, the components α 2 , δ 2 , α2 , δ2 of A 2 are not affected by interpolation errors, which are usually introduced when an attributable is obtained from the astrometric observations.
Assume we can model a close approach at epoch t 1 by an instantaneous velocity change, like in Öpik's theory ( Öpik, 1976).Then, our algorithm can be applied using P 1 , A 2 as input data, because the values of the Keplerian integrals are conserved in (t 1 , t 2 ].At t 1 , when the velocity instantaneously changes, the position vector P 1 remains unchanged. In practice, the close encounter is not instantaneous, but takes a finite time: we select the epoch t 2 when the close approach is over, as described below.In a typical close encounter, the value of Tisserand's parameter T changes quickly, going back approximately to its pre-encounter value when the encounter is over.In Cavallari et al. (2023, Sec.3), for a given value C of the Jacobi constant, the authors provide a value d C of the geocentric distance attaining | dT dt | ≤ ε, for a fixed small parameter ε, so that we can think that the encounter is over when the asteroid reaches that distance from the Earth.Following this approach, we choose t 2 as the time when the value of the distance becomes d C .
In Fig. 1 we show the effect of the close approach on these orbits.On the left, we plot the distribution of the relative differences between the semimajor axes a 1 , a 2 obtained by the known orbits propagated at epochs t 1 and t 2 using the software OrbFit 5 .On the right, we plot the distribution of the distance δ(T 1 , T 2 ) between the propagated trajectories T 1 , T 2 , where In ( 19) the subscripts 1, 2 of the Keplerian elements a, e, i, Ω, ω and the trajectories T refer to the epochs t 1 , t 2 .In most cases the relative change in semimajor axis is within 20%.However, there are a few cases where the change is larger, with one extreme case (Fig 1,left) where it passes from a 1 ∼ 2.84 to a 2 ∼ 0.69 au.The approximate values of the mean and the standard deviation of (a 1 − a 2 )/a 1 are about 0.012 and 0.0478, respectively.Next, we show the performace of our algorithm using P 1 , A 2 as input data, which have been computed for each of these orbits by a full n-body propagation with OrbFit.In case of multiple solutions, we select the best one according to the procedure explained in Section 6.1.
In 8 cases, out of 1305, we could not obtain an orbital solution.However, by increasing the time span [t 1 , t 2 ] by 10%, while keeping the same value of t 1 , thus changing only A 2 , we recovered an orbital solution in all these cases.In one case the computation of the second epoch t 2 with Newton's method failed.However, employing the starting guess value for t 2 , computed with a geocentric two-body dynamics, we obtained an orbital solution also in this case.
In Fig. 2 we show the distribution of the differences between the Keplerian elements computed with our algorithm (with subscript c) and the same elements at t 2 obtained by propagation.In Fig. 2 (bottom right) we also show the distribution of the distances δ(T 2 , T c ) between propagated (T 2 ) and computed (T c ) trajectories, where the function δ is the same as in ( 19).
The results of this preliminary statistical test are satisfactory: the mean and the standard deviation of the distributions displayed in Fig. 2  Finally, we show the dependence of the results of our algorithm on the type of close encounter.Following Cavallari et al. (2023, Sect. 3), close encounters can be classified as deep/shallow and fast/slow, according to the values of the minimum distance q from the Earth center and the corresponding geocentric velocity υ.Here, we compute an approximation of q and υ, still denoted by these symbols, using the following procedure: q is the pericenter distance of the geocentric two-body orbit computed from the Keplerian elements at t 1 , and υ is the corresponding two-body velocity.In Fig. 3 we plot the values of q and υ using a color scale representing the values of the trajectory distance δ(T c , T 2 ).
For values of q greater than 0.001 au, we obtain better results (e.g.lower values of the trajectory distance) for higher values of the velocity υ.For q smaller than 0.001 au, the results appear worse, and there is no apparent correlation with velocity in the figure.

Conclusions
We have introduced a method to compute preliminary orbits with one topocentric position vector P 1 and a very short arc of optical observations, from which we can derive an attributable A 2 .This method is based on polynomial equations coming from the first integrals of Kepler's problem, i.e. angular momentum, energy, and Laplace-Lenz vector.
Using the conservation laws of these integrals, after introducing the auxiliary variable z 2 , we obtain a polynomial system that always has solutions, at least in the complex field, even if P 1 , A 2 do not correspond to the same celestial object.There are some checks that can be performed to accept or reject solutions.
We applied this algorithm to the computation of 1305 NEA orbits, whose phase was changed in order to enhance the close encounter effect with the Earth: these preliminary results are satisfactory, see Section 7. The ideal situation for the application of the proposed algorithm in this context would be given by an instantaneous effect of the close encounter, which is not the case.The selected epoch t 1 does not exactly correspond to the time of passage at the MOID because we use a full n-body propagation starting from t 0 , so that at t 1 both the osculating trajectory and the time law along it may have changed.Moreover, from this preliminary test, we saw that also varying t 2 may affect the results.Therefore, a more detailed study is necessary to better understand the applicability in case of close encounters and the reliability of the computed orbits.On top of that, the sensitivity of the algorithm to astrometric errors has still to be investigated.
Applications of the method introduced in this work to the orbit computation of Earth satellites undergoing a maneuvre are also possible, and are worth to be investigated in a future work.

Figure 1 :
Figure 1: Effects of the close approach on the semimajor axes (left) and on the whole trajectories (right).

Figure 2 :
Figure 2: Distribution of the differences in orbital elements a (top left), e (top right), i (center left), Ω (center right), ω (bottom left) between the propagated trajectory and the one computed with our algorithm at epoch t 2 .Distribution of the distances δ(T 2 , T c ) between propagated and computed trajectories (bottom right).

Figure 3 :
Figure 3: Values of q and υ in a log-log plot.The colors represent the values of the distance δ(T 2 , T c ) between known and computed trajectories.
for some bivariate polynomials h 4 , h 6 , h 7 .For later use we observe that ij .

Table 1 :
are shown in Table1.Mean and standard deviation of the distributions represented in Fig.2.