PRECISE LINE-OF-SIGHT MODELLING FOR ANGLES-ONLY RELATIVE NAVIGATION

This work presents a precise analytical model to reconstruct the line-of-sight vector to a target satellite over time, as required by angles-only relative navigation systems for application to rendezvous missions. The model includes the effects of the geopotential, featuring: the analytical propagation in the mean relative orbital elements (up to second-order expansion), the analytical two-way osculating/mean orbital elements conversion (second-order in J 2 and up to a given degree and order of the geopotential), and a second-order mapping from the perturbed osculating elements’ set to the local orbital frame. Performances are assessed against the line-of-sight reconstructed out of the precise GPS-based positioning products of the PRISMA mission. The line-of-sight modelled over a far-range one day long scenario can be ﬁtted against the true one presenting residuals of the order of ten arc-seconds, which is below the typical sensor noise.


INTRODUCTION
Angles-only navigation plays a relevant role to treat the problem of space debris in the Low Earth Orbit (LEO) region. Regarding active debris removal, as recently shown by the AVANTI (Autonomous Vision Approach Navigation and Target Identification) in-flight demonstration, 1 an onboard autonomous angles-only relative navigation system is a convenient solution to safely approach a noncooperative flying object till the distance range where the full pose becomes strictly necessary. As for space situational awareness, spaceborne relative orbit estimation based on angles-only observations can be exploited by space-based architectures, to complement the existing ground-based services. 2 The relative orbit estimation problem from bearing-only observations is weakly observable. At practical level, during a rendezvous, few sporadic manoeuvres can be performed to disambiguate the possible solutions in range. Several theoretical studies focus on the more convenient manoeuvring strategies to improve angles-only observability. [3][4][5][6] In-flight demonstrations of such approach are provided by ARGON (Advanced Rendezvous Demonstration using Global Positioning System and Optical Navigation) 7 and AVANTI, 1 where the manoeuvres executed to perform the rendezvous also supported the convergence of the relative navigation solution. Despite the feasibility of this method, the unknown manoeuvre execution errors worsen the achievable navigation accuracy. An alternative approach to improve the observability property of the problem, valid also for manoeuvre-free arcs, is to consider the non-linearities introduced by perturbations (e.g., through mean to osculating elements' transformations for the geopotential) and by the orbit curvature in the modelling of More in details, the relative motion is propagated analytically in the doubly-averaged Relative Orbital Elements' (ROEs) space, through closed-form first-order state transition matrix, including the secular effects due to J 2 , J 2 2 , J 4 , and J 6 , 11 and through a second-order state transition tensor, accounting for the J 2 effects. The adopted formulation is valid for whatever eccentricity of the reference orbit, and outperforms the available ROE-based first-order models. 1,12 At the same time it improves the Gim-Alfriend 13 and Yang et Al. 14 relative motion models, either in order a/o in considered perturbations, while preserving the compact formulation deriving from the parametrisation in ROEs. The two-way conversion between osculating and mean orbital elements is carried out through an analytical and compact algorithm that combines a second-order Lie-based approach to cancel the J 2 effect, to the Kaula's linear method for the remaining terms of the geopotential. 15 Indeed, this improves the overall modelling accuracy by drastically reducing the artificial drift introduced by the transformation errors, regardless the moment of the orbit when the conversion takes place. Lastly, the non-linear line-of-sight is recovered through a second-order expansion from the osculating ROE set to the rectilinear Cartesian local orbital frame. This mapping improves other available algorithms, 13, 14 either in order a/o in accuracy performance. Moreover, a compact formulation is proposed exploiting the properties of three-dimensional rotations in the expansion of the anomalies.
Accuracy results are provided by comparing the modelled line-of-sight against the one reconstructed out of the GPS-based relative positioning products of the PRISMA mission, 16 which took place in an orbit environment highly representative for future active debris removal missions. Few data sets have been selected to represent critical conditions for the relative navigation system at far-range. Such data sets have been used to assess the overall modelling accuracy as well as to compare the proposed model against available methods from the literature. Results show that the proposed fully analytical modelling of the line-of-sight presents residuals of few arc-seconds when fitted against true data sets over an entire day.
After this introduction, a first section presents the modelling framework, with special focus on the development of the mapping from osculating ROEs into the relative position in the orbital frame. The following section collects the results of the performed analysis.

LINE-OF-SIGHT MODELLING
Following the notation of Gaias et Al., 17 the angles-only observations at each instant of time consist of two angle measurements, i.e. azimuth η and elevation ψ, which subtend the Line-Of-Sight (LOS) unit-vector to the target satellite expressed in the sensor frame u s (as depicted in Figure 1 of Reference 17): Neglecting the camera offset, the LOS unit-vector is related to the relative position x between target and chief satellites in the rectilinear co-moving orbital radial-tangential-normal (RTN) frame (centred on the chief) by: where R s RTN denotes the rotation from RTN to the sensor frame, in this work assumed as: In the following, α = (a, u, e x , e y , i, Ω) T is the set of osculating Keplerian non-singular elements, with a the semi-major axis, u = ω + M the mean argument of latitude, ω the argument of the perigee, M the mean anomaly, e x = e cos ω and e y = e sin ω the x and y components of the eccentricity vector, and i the inclination. The dimensionless relative state in ROEs δα is defined as: where ∆· denotes the difference between quantities of the target and chief c satellites, δλ is called the relative mean argument of longitude, and the vectors (δe x , δe y ) T and (δi x , δi y ) T are respectively known as the relative eccentricity and inclination vectors.
The analytical modelling of the angles-only observationsz is obtained executing the chain of actions depicted in Figure 1. This algorithm encompasses three functional components, namely: the propagation of the relative orbit in the mean space (core part within vertical lines), the conversion between mean/osculating orbital elements (T 2 and T −1 2 ), and the mapping from osculating ROEs at time to the relative position x in the RTN frame (T 4 •T 3 ).  The absolute orbit of the chief satellite is known in the Earth Mean Equator and Equinox of J2000 (EME) reference system. Thus, to obtain the osculating orbital elements (OEs) α c , first the rotation R from EME to the true of date reference system is performed. Secondly, the Cartesian absolute state is transformed into the corresponding set of osculating elements through the transformation T 1 .
The OEs transformations, T 2 for the direct conversion from mean to osculating elements, and T −1 2 for the inverse one, are performed analytically thought the KA-l×m algorithm of Reference 15. It combines a Hamiltonian approach applied to the J 2 problem to the second-order with Kaula's linear perturbation method for the remaining terms of the geopotential, being l and m respectively order and degree of the geopotential terms accounted in the corrections.
Once in the mean space, the only orbital elements that present a secular variation are are Ω, ω, and M due to spherical Earth (M ) and to even zonal harmonics only (all). Therefore, the analytical model of the mean relative motion in ROEs is obtained by performing a Taylor expansion of the mean chief orbit. By retaining only the first-order term, the state-transition matrix Φ is derived, including the secular effects due to J 2 , J 2 2 , J 4 , and J 6 . By performing the expansion to the second order, the state-transition tensor Ψ is derived. In this work it accounts for the unperturbed and J 2 (to the first-order) terms.

Mapping the ROEs into the RTN frame
The mapping from the osculating ROEs at time to the relative position in the local rectilinear RTN frame is performed into the following steps. First, the transformation T 3 is required to pass from δα to the set of elements ∆oe = (∆a, ∆θ, ∆i, ∆e x , ∆e y , ∆Ω) T , where θ is the true argument of latitude. This is needed since the ROEs are defined using the relative mean argument of longitude (i.e., a function of M ), whereas the observations are taken on the true osculating orbit. To the first-order, the transformation T 3 is given by: with: where η = 1 − e 2 x − e 2 y , R is the distance between chief and Earth, p is the semi-latus rectum, n is the mean motion, and V r and V t are respectively the radial and transverse components of the chief orbital velocity. For small values of the orbit eccentricity, which is the case of active debris removal applications a/o of several formation-flying activities in LEO, the true argument of latitude θ can be computed from M using series expansions in the eccentricity: 18 In this work only the terms up to e 2 are kept.
The transformation T 4 first maps the relative state ∆oe at time into the local curvilinear orbital frame. Afterwards, such quantity is expressed it into the local rectilinear RTN frame. Accordingly, the first phase of T 4 corresponds to the geometric transformation Σ(t) of Gim and Alfriend, 13 when the expansions are carried out to the first-order. Its extension to the second-order, instead, is given by Yang at Al. 14 A compact version of the second-order mapping, moreover generalized to include the effects of zonal terms higher than J 2 , is derived here as follows. Note that for the approach of Figure 1, only the relative position is required, whereas here the complete transformation is discussed.
The mapping from ∆oe to the relative state (position and velocity) in the local curvilinear orbital frame is obtained by equating the osculating inertial position and velocity of the deputy satellite written in the chief-local frame to the expressions derived from a geometric transformation obtained by expanding the chief inertial position and its direction cosine matrix. 13,14 This latter quantity defines the orientation of the local orbital frame of the chief C with respect to the inertial one I, and it is given by the Euler 3-1-3 rotation function of θ, i, and Ω. By defining Λ = R I C , (i.e., the rotation from C to I), Λ(θ, i, Ω) is a 3-directional rotation ∈ SO 3 and from the orthogonality property it derives that: where the notation δ• (1) is used for the first-order expansion in θ, i, and Ω, which is the virtual variation of Λ. The vectorg associated to the skew symmetric matrix [g×] (with g i as i-th column), is given by: The second-order expansion of Λ produces the matrix F: Thus, the second-order mapping, 14 to deliver the curvilinear relative state (x,ẋ) T , can be compactly written as: Here, the terms deriving from the expansion of quantities in the orbital plane are given by: with ∆ oe = (∆a, ∆θ, ∆e x , ∆e y ) T , ∇• the gradient, and H • the Hessian. Whereas, the columns of F required in Eq. (11) are here explicitly written for components: The remaining term is the perturbed angular rate of C written in the local frame. Accordingly, it would be function ofΩ,θ, andi, computable from the Lagrange planetary equations subject to a given disturbing function. Nevertheless, taking into account the osculating orbit constraint (i.e., t = 0), for a conservative system, it simplifies to: where µ ⊕ is the Earth's gravitational parameter. The expression ofΩ J 2 is provided in Eq. (14) of Reference 13. The contributions to the radial component of the angular rate due to J 3 and J 4 are given by: where R ⊕ is the mean Earth's radius.
The use of a curvilinear local frame allows obtaining a precise result in the along-track direction, 10 and thus it is fairly accurate for large bounded relative orbits. 13 Nevertheless, for large along-track separations, the following correction of the radial components becomes necessary to account for the curvature of the orbital path: where ρ =x,ρ =ẋ, ϑ =y/R, andθ =ẏ/R. Equation (16) transforms the curvilinear radial component into a rectilinear frame. The so obtained relative position x = (x,y,z) T in RTN is used in Equations (2) and (1), to produce the modelled observationsz.
Note that, in order to overcome the limitation in accuracy of so far available LOS modelling based on the mapping of an OE-based relative state, the LOS unit-vector was usually computed by retrieving the absolute state of the target satellite (i.e., applying a T −1 1 to the osculating α d ), by subtracting to it the absolute state of the chief, and then by rotating the obtained relative state in the chief-centred RTN frame. This non-linear transformation produces the exact relative state and therefore will be used to compute the true LOS out of the true states of the satellites. The step T −1 1 requires the numerical solution of the Kepler's equation; whereas the RTN frame centred on the chief is given by Eq. (2) of Reference 9. Examples of exploitation of this non-linear transformation are provided by the onboard navigation system of the AVANTI experiment 19 or by the Algorithm 1 of Reference 2.
For the specific application under study, (i.e., LOS modelling for near-circular relative motion), a simplified non-linear transformation to deliver only the relative position in the RTN frame is given by: where Λ • is the Euler 3-1-3 rotation introduced before, R • is the satellite-Earth distance, and θ is computed through Eq. (7). This simplified non-linear transformation exploits the fact that the absolute osculating elements α c and α d are available in the analytical framework of Figure 1 before computing δα. In the near-circular far-range cases discussed in the next section, Eq. (17) introduces an error with respect to the exact non-linear relative state at sub-millimetre level in normal and radial directions and sub-decimetre level in along-track direction. Accordingly, this output can be regarded as accurate as the true LOS, and thus, Eq. (17) represents an alternative to the mapping of the ROEs into x through T 4 •T 3 , when the latter is not accurate enough.

RESULTS
To assess the accuracy of the proposed analytical modelling of the LOS, in this section a comparison is performed against the line-of-sight reconstructed out of the GPS-based relative positioning products of the PRISMA mission. These products are accurate to the sub-centimetre level, 16 and thus represent the true orbit of the satellites of the PRISMA formation. In addition, given its orbital scenario, PRISMA is extremely representative for future missions exploiting LOS navigation in LEO (e.g., active debris removal missions), considering the weak effect of the differential aerodynamic drag. 20 Among the available products, these two data sets have been selected: • 5-Mar-2011, see Figure 2, with 5 hours of drift with large relative semi-major axis (i.e., aδa > 500 m); • 17-Feb-2011, see Figure 3, with 5 hours of bounded relative motion at large relative mean longitude (i.e., aδλ > 30 km).
In the aforementioned plots the ROEs computed from the precise orbit determination (POD) products are plotted over time. The occurrence of manoeuvres is marked through vertical dashed lines.  The first analysis regards the accuracy of the proposed mapping T 4 •T 3 , obtained through Equations (5)-(7) and (11)-(13) and (16). Accordingly the true osculating ROEs are taken as input and the results are assessed in the forms of errors w.r.t. the true LOS. In order to show the improvements of the proposed method, the following available algorithms are employed for comparison: • GSOC: the mapping developed at the German Space Operations Center 21 and used in the first prototype of ROE-based angles-only relative navigation filter 17 employed in the ARGON experiment; 7 • GA-curv: the first-order Σ transformation of Gim and Alfriend, 13 to deliver a relative state in the curvilinear orbital frame; • GA-rect: the GA-curv corrected by Eq. (16); • YLZ: the second-order mapping of Yang et Al., 14 to deliver a relative state in the curvilinear orbital frame.    17) is not plotted, since it is basically coincident with the exact solution, thus it would provide a constant zero error in both figures. Note that, given the sensor orientation of Eq. (3), the y-axis of the sensor frame is directed as -R, whereas the x-axis to the -N. By referring to Figures 4 and 5, the accuracy in the normal direction is comparable among the methods, as shown by e Az and e N . In the along-track T direction, the GSOC mapping achieves a poor result since it works with the mean orbit (using the osculating mean argument of latitude u); whereas the remaining methods are all very accurate, thanks to the use ofy. In the radial direction the proposed method is much more accurate, thanks to the correction of the orbital curvature. Note that, for its derivation, the YLZ and T 4 are the same in this application (i.e., the relative velocity is here not used) to the net of Eq. (16). This explains the difference only in radial direction, decreasing over time for the data set of March 2011 (i.e., large δa case). The GA mappings, instead, suffer from the lack of the second-order correction (in both data sets aδλ is greater than 10 km) and, for the curv case, from the lack of the correction of Eq. (16). As a whole, the proposed mapping (in black in the plots) remarkably outperforms the others relative OE-based mappings, when dealing with large relative orbits. At the same time all these algorithms that map the ∆oe state into the relative position x introduce an oscillating error in the normal direction proportional to the size of the maximum displacement. This can be noted by relating the results of e Az and e N in Figures 4 and 5 to the magnitude of the relative inclination vector, whose components are shown in Figures 2 and 3. This is due to how thez is computed, which realizes an approximation of the rectilinear true z. In far-range rendezvous scenarios, usually the out-of-plane size of the relative orbit is already reduced to ≈1 km (at conclusion of the orbit phasing transfer). In these cases, the maximum azimuth error is less than the typical noise of available sensors. For example, the camera employed for ARGON and AVANTI exhibits a line-of-sight noise of about 40 arc-seconds at far-range (corresponding to less than half-apixel). 10,22 On the other hand, in case of very large out-of-plane motions, the analytical non-linear simplified transformation of Eq. (17) can be used instead of the mapping. This allows basically removing the error source related to the modelling of the measurements, at the cost of working with absolute OEs instead of with ROEs, which reduces the geometrical insight in the relative problem.
The second performed analysis concerns the overall accuracy obtainable through the analytical modelling of the algorithm sketched in Figure 1. In this case, the accuracy results from the performances of the OEs conversions (i.e., T 2 and T −1 2 ), of the mean relative orbit propagation (i.e., Φ and Ψ), and of the just evaluated mapping. In order to complement the overall error budget with the single error contributions, several models are again considered. By referring to Table 1, M1 is the model adopted in ARGON, though here using a slightly improved first-order relative motion model. M2 is the framework exploited by the authors to deal with angles-only initial relative orbit determination in Reference 10, though here applying the correction of Eq. (16) to the GA mapping. The next two cases are a realization of the current framework. In particular, M3 is the computationally lightest version accounting for only J 2 ; whereas M4 includes geopotential effects up to order-6 degree-6. Finally, M5 is a variation of M3, where the transformation of Eq. (17) is used instead of the ROE-based mapping T 4 •T 3 . This case has been introduced to isolate the error contribution due to the analytical propagation from the error introduced by the mapping.  By considering the data set of March 2011 (i.e., large δa case), the overall accuracy measured in errors in azimuth and elevation w.r.t. the true LOS is reported in Figure 6-left. The poor result of M1 motivated the use for AVANTI of the non-linear LOS reconstruction for the on-board relative navigation system, 19 as well as of the numerical integration for the ground-based precise relative orbit determination layer. 9 The error in elevation for M2 is mainly due to the first-order only map-ping. This is more evident by looking at Figure 7, where the osculating ROEs are compared. Since all models M2-M4 exploit OEs conversion algorithms second-order in J 2 (M5 is not mentioned as it is equivalent to M3 for the propagation), the accuracy at this point is almost the same (see Reference 15 for more details about the effects of the T 2 algorithms). It should be emphasized, however, that the method of M2 requires numerical iterations for the inverse transformation, whereas M3-M5 are fully analytical. At ROEs propagation level, the improvement brought by employing the KA-6×6 algorithm instead of KA-2×0 is visible in the accuracy of the aδλ component: a more accurate value of aδa at the initial time reduces the along-track error over time. 15 This can also be noted in Figure 6-right, sub-plot e T . However, since in KA-6×6 the Kaula-based corrections are performed only on the semi-major axis component, some residuals oscillations appear in the error in radial direction. As a matter of fact, the trade-off between M3 and M4 regards the achievable gain in along-track precision against the increase of computations to carry out the periodic corrections (due to geopotential terms higher than order-2 degree-0) in the Kaula phase. The results of M5 follow the trend of the ones of M3, since they employ the same propagation model, with reduced amplitude of the error oscillations due to the use of Eq. (17). This is more evident from Figure 8, where the plots refer to the data set of February 2011, given the larger magnitude of the inclination vector (models M1-M2 are not shown to focus on small error values). Note that being the relative motion almost bounded, the accuracy gain in along-track using M4 becomes negligible. By referring to the e R plots of Figures 6 and 8  So far the comparative analysis has been performed taking a fixed initial state at a randomly chosen initial time. Figure 9, instead, presents the observation residuals obtained when the LOS modelled through M3 is fitted against the true values from the POD products. For both data sets, the residuals amount to few arc-seconds over the 5 hours, when the computationally lightest model is used. Moreover, the residuals in azimuth are almost of the same magnitude, despite the difference in size of the relative inclination vector for the two data sets.
In relative orbit determination problems, an important aspect is represented by the choice of the length of the data set to be processed. This is generally related to the trade-off between requirements from the data editing and errors introduced by the propagation method. Accordingly, a last analysis is here performed regarding the LOS modelling accuracy over extended manoeuvre-free arcs.  this end the additional data set depicted in Figure 10 is considered: on that day no manoeuvres were performed and the orbit drifted to decrease of circa 10 km the along-track separation. By modelling the LOS through M3-M5 over the whole day, the errors in relative position are given in Figure 11left. One can note that, M4 achieves a better result than M3. Though a certain along-track error is accumulated due to the non-modelled effects acting on the relative dynamics. Despite this, the LOS fitting based on the computationally lightest propagation option M3, achieves observation residuals within ±5 arc-seconds. This score has the same order of magnitude of the one achieved by M5, and lies well within the 40 arc-seconds noise threshold of the camera sensor employed at far-range for ARGON and AVANTI. 10 Figure 11. LOS modelling (left) and LOS fitting (right) via M3 and M5 for the data set of Figure 10.

CONCLUSION
This work presented an analytical model of the line-of-sight between two neighbouring satellites as required for angles-only relative navigation systems. The methodology is valid for large far-range relative orbits in the low Earth orbit region. Achievable accuracy and light computational load make the proposed methodology very convenient for future spaceborne applications.
To the overall error budget contribute: the error in the propagation of the relative motion and the error in computing the relative state from the orbital elements' based relative state. The first contribution is minimized considering at least the J 2 effect to the second-order when moving into the doubly-averaged orbital elements' space. Afterwards, the propagation of large relative orbits remains accurate over long time, since the relative motion model includes the first-order state transition matrix and the second-order state transition tensor. At this stage the error amounts to about 5 arc-seconds in radial and normal direction. The error in along-track direction depends on the size of the propagation horizon, due to the effects of the so far not modelled non-conservative perturbations on the relative dynamics. The error introduced by the mapping of the osculating relative orbital elements into the relative position expressed in the local orbital frame depends on the relative orbit geometry. In-plane components are modelled very precisely thanks to the the inclusion of the second-order expansion in the orbital elements and to the modelling of the curvature of the orbital path. The out-of-plane component is affected by an error function of the maximum size of the normal displacement. For far-range scenarios with relative inclination vector up to 1 km, the observation residuals, fitted with respect to the true line-of-sight reconstructed from flight data, lie within a 10 arc-seconds threshold. This figure is way less than the typical noise of camera sensors employed for relative navigation purposes.

ACKNOWLEDGMENT
The work performed at Politecnico di Milano has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 679086 -COMPASS).