Vision-based localization of the center of mass of large space debris via statistical shape analysis

The current overpopulation of artificial objects orbiting the Earth has increased the interest of the space agencies on planning missions for de-orbiting the largest inoperative satellites. Since this kind of operations involves the capture of the debris, the accurate knowledge of the position of their center of mass is a fundamental safety requirement. As ground observations are not sufficient to reach the required accuracy level, this information should be acquired in situ just before any contact between the chaser and the target. Some estimation methods in the literature rely on the usage of stereo cameras for tracking several features of the target surface. The actual positions of these features are estimated together with the location of the center of mass by state observers. The principal drawback of these methods is related to possible sudden disappearances of one or more features from the field of view of the cameras. An alternative method based on 3D Kinematic registration is presented in this paper. The method, which does not suffer of the mentioned drawback, considers a preliminary reduction of the inaccuracies in detecting features by the usage of statistical shape analysis.


Introduction
Since the launching of Sputnik, the first artificial satellite, by the USSR in 1957, a very large number of satellites have been launched into space by several countries and with several purposes. In particular, from the second half of the twentieth century up to the present, the number of artificial satellite in Low-Earth Orbit (LEO) has dramatically increased. Some of them are now inactive and uncontrolled. As the number of failed spacecraft increases significantly every year, also the risk of very dangerous collisions between relatively large space objects increases. A very famous collision that happened in 2009 between Cosmos 2251 and Iridium 33 caused an environmental disaster due to the fragmentation of the two bodies [1].
In this context major space agencies such as NASA and ESA are interested in the possibility of removing relatively large space debris and they planning the first missions. Under this spur, several concepts of debris removal systems were developed [2]. It is common in particular the idea of docking passive spacecraft with active chasers spacecraft.
For a good planning of the docking between spacecraft, it is important to estimate the complete relative dynamic state between the chaser and target spacecraft. In particular, the precise estimation of the position of the center of mass of the target leads to the possibility of predicting torques coming from the interaction forces between spacecraft.
This estimation is particularly challenging because of the absence of helpful sensors that allow the retrieval of information about the kinematics and dynamics of the target spacecraft. Astrometric measurements from ground-stations are essential to determinate approximately the orbit of the targets, but the accuracy of those measures is only sufficient for the purpose of planning the rendezvous with the chaser spacecraft. Some examples of methods to make an initial orbit determination of a space object are presented in [3] and in [4]. Once the chaser reaches an orbit that is close to the one of the target, a proximity tracking of its natural features become possible if the chaser equips visionbased sensors.
There are in literature some works regarding the complete state estimation of passive targets that rely on 3D sensors. Tracking is pursued in real-time with different techniques but the possession of a 3D CAD model of the target is required. For instance, in [5] an active laser ranging device is adopted to pursue the estimation of the complete state of a known but non-cooperative spacecraft. The ICP algorithm and a Kalman filter are combined together in a closed-loop fashion to make the method tolerant of potential failed detections of the objects. Active sensors represent a fairly reliable solution for tracking the motion of space debris but, as a drawback, they slightly add complexity in the chaser design due to increasing energy consumption and, secondarily, costs.
Against, looking to the state of the art, it is possible to find examples that demonstrates that passive sensors maintain comparable performances to active ones. One of the most significant works is illustrated in [6]: the dynamic relative state of two non-cooperative satellites is estimated through the application of an iterated extended Kalman filter, which is fed with measurements from a stereo-rig system on the chaser spacecraft. Moreover, no prior assumptions about the shape of the target are required by this estimation technique However, the last mentioned solution, as any other in the state of the art, do not consider the possibility of losing temporary the availability of measurements or having modifications of the measured quantities. To better explain this last concept, imagine that one or more features disappear suddenly from the field-of-view of the stereo-rig system; at the same time, the tracking sensors could detect completely new features. Besides, also occlusions can happen, compromising totally the acquisition of the measurements. In these conditions, no known methods are able to recover the complete dynamic state of the non-cooperative target.
Some previous works of the authors of this paper ( [7], [8], and [9]) contemplates the mentioned tracking conditions and provides methods for recovering the attitude dynamics of space debris, but not the orbital one. Therefore, the literature lacks of methods that consider difficult tracking conditions for estimating the location of the center of mass (CoM) of space debris.
The aim of the present paper is to introduce an approach that does not rely on typical recursive Bayesian observers, but considers only the registration of the target kinematics to locate the CoM, assumed coincident with the center of gravity of the passive object of interest. In particular, this approach relies on the identification of the axes and poles of the successive finite rotations of the target in several different time intervals.
Considering no relative translation between the chaser and the target, all the axes intersect in a point of the body that coincides with the CoM. When a relative translation exists, the axes do not pass any more for the CoM but its distance from the axes is positively related to the relative translational velocity. Nevertheless, if the entity of the translation is within certain bounds, it is possible to find an efficient geometric estimator of the CoM location.
Because of the evident unavailability of real measured data, the proposed method is validated on simulated data. These data include the presence of a noise that is coherent to the uncertainty associated to plausible stereo-vision systems for space applications.
The problem of finding one screw axis between two poses has not an exact solution when some noise corrupts the coordinates of the points representing the poses. This issue is often faced by estimating an optimal screw axis in the least-squares sense. On the contrary, in this paper the authors propose to filter the data in a pre-processing stage to obtain series of rigid clusters of homologous points belonging to the target. This approach corresponds to infer from raw data the original shape of the examined body. In this way, a unique screw axis underlies the finite movement between two filtered poses. While the common least-squares estimation techniques lead to a total failure in localizing the CoM, the proposed method is useful for a reliable recovery of the position of the searched point. The remainder of this paper is organized as follow: section 2 describes how to simulate the measurements for the method validation. Section 3 illustrates a simple technique to identify screw axes from the coordinates of three couples of homologous points. Then, it shows a possible estimator of the CoM location. The proposed estimator is tested on ideal data (not corrupted by noise or missing samples due to occlusions) and its accuracy is briefly discussed. Section 4 presents one of the most famous techniques for the statistical shape analysis of homologous point clusters: the general Procrustes analysis (GPA). Section 5 shows the results of the proposed methodology applied on simulated data. Finally, section 6 provides the conclusions of this work.

Tracking simulation
The methods that are presented in this work consider as input the positions of several features of the target in many time samples. These data are assumed to be measurable by using a stereo-rig system that is mounted on a chaser spacecraft. Thus, to create a realistic set of input data it is necessary to simulate the relative orbital motion between two close satellites. Moreover, also the attitude dynamics of the target have to be modeled. Indeed, considering the ideal case in which no relative orbital dynamics exists between the satellites, the chaser would see the features of the target rotating around the searched CoM. Otherwise, the trajectories of the features would be the result of the combination of both the target rotation and translation from the chaser point of view.
The attitude dynamics of a satellite is well-approximated by the Euler's equation for torque-free rigid bodies: where is the principal inertia matrix of the target and is the angular rate of the target expressed in a principal body-fixed coordinate system . In the entire paper, the upper prescript on a vector indicates the coordinate system in which the vector is expressed.
The resulting angular rate from equation (1) determines the variation of the attitude of the target according to the following Darboux's equation: where is the orthonormal direction cosines matrix that maps a vector expressed in an inertial frame into the same vector expressed in the coordinate system . The notation × indicates the following skew-symmetric matrix: Equations from (1) to (3) are useful for the simulation of the attitude dynamics of the target. For simulating the formation of chaser and target, it is convenient to consider the motion of the chaser, assumed as a material point, about the CoM of the target. This mentioned relative dynamics is given by the following equation: where ℓ is the position of the chaser in the local-vertical-local-horizontal (LVLH) ℓ coordinate system of the target, is the true anomaly of the target, is the planetary constant, is the distance between the CoM of the target and the Earth's center. Finally, has the following expression: The position of the target CoM can conveniently be expressed in a coordinate system that is placed on the chaser and whose axes are aligned to the ones of the inertial Earth centered system : Once the initial conditions for equations (1), (2), and (4) are set, it is possible to evaluate at any sample time the position of any feature whose known constant position in the body-fixed reference frame is . Hence: For a better comprehension of the mentioned quantities and coordinate systems refer to figure 1. The computation of the second member of the equation (7) is the last step of the data generation procedure, which includes the solution of all the mentioned differential and algebraic equations. To guarantee the robustness of the presented approach, we split the simulation time in several intervals in which different cluster of features are considered. In addition, for some time intervals, we considered a complete unavailability of data simulating the occurrence of occlusions. Figure 2 illustrates this concept. Besides, figure 3 shows an example of the computed trajectories for five particular features during the period of one orbit of the target. From that figure it is possible to evidence the joint effect of the attitude and orbital dynamics.

3D kinematic registration of rigid poses
The kinematic registration is a procedure that consists in finding one synthetic representation of a possible motion that happened between the current pose of a rigid body a past one of the same body. Thus, it is a procedure that produces an alternative representation of the motion of a rigid body. Actually, performing the kinematic registration of a body does not mean to estimate the real evolution of the posture of the body. However, it will be shown that performing such registration leads to information that is useful to estimate the unknown location of the CoM of the target. The very well-known Mozzi-Chasles' theorem states that two poses of a rigid body can be mapped one with respect to the other through one translation and one rotation about the direction of the translation. Thus, it is possible to evidence the existence of one axis, whose direction and location characterize the displacement of a rigid body. If this displacement happens during an interval of time whose duration tends to zero, then this axis is identified as the instantaneous screw axis. Otherwise, if the duration is finite, it is named as the finite screw axis.
To characterize the pose of one object it is sufficient to know exactly the Euclidean coordinates of three points that belong to it. Indeed, it is possible to build a reference triad of orthogonal axes by making appropriate cross products of vectors connecting the points. Thus, the knowledge of the positions of three couples of homologous points at two separate instants is sufficient for the computation of the finite screw axis. There are several methods in the literature that address this computation (see [10] for instance). This paper considers a previously developed method consisting of the 3D generalization of the planar Reauleaux's rule [11]. This rule individuates the pole of planar rotations by intersecting the bisecting lines connecting two couples of homologous points. Therefore, according to this rule, two poses of a rigid body overlaps through a single planar rotation. Actually, if the real motion between the poses was a pure rotation, the found pole would be exactly the CoM of the body. On the contrary, the existence of a translational motion leads the found pole distancing from the actual CoM. In particular, in the instantaneous case, the distance is proportional to the ratio between the linear and the angular speed of the body. In the 3D case, two poses overlaps through one translation and one rotation. However, the rotation develops on a plane, which is perpendicular to the direction of the translational displacement. This means that the Reuleaux's rule remains applicable once the effect of the translation is removed in some way. For instance, given three vectors connecting a couple of homologous point ( and ′ ), the vectors 1 = 2 − 1 and 2 = 3 − 1 , span a plane which is perpendicular to the screw axis . Thus, the direction of the screw axis is given by the following expression: Equation (8) provides the direction along with projecting two of the three couples of homologous points on the plane spanned by 1 and 2 . This projection eliminates the translation, so the mentioned Reulaeaux's rule holds on the projected points for the computation of the pole of the remaining rotation. This pole is one of the infinite points belonging to the screw axis.
Differently from the pure Reauleaux's rule, its 3D generalization admits a translational displacement to overlap two different poses of a rigid body. However, this translation must be along . Actually, if the real motion between the poses involved a misaligned translational component, the screw axis would not pass for the CoM of the body. In particular, in the instantaneous case, the distance from the axis to the CoM is proportional to the component of the linear speed that is orthogonal to the axis.
Furthermore, even if no translations occur, the generalization does not individuate the CoM but only a screw axis. Thus, even in the case of pure rotations, it is necessary to compute two screw axes considering three different poses. The axes will intersect in the searched CoM. However, in the general case, two successive screw axes are non-coplanar; so, the CoM cannot be individuated.
Assuming that the translational motion is relatively slow, as it should be the relative orbital motion of two satellites after rendezvous, the authors of this paper proposed to perform a pseudo-intersection between successive screw axes [12]: the CoM is estimated as the middle point of the shortest line segment connecting the non-coplanar axes (see figure 4). The work in [12] showed that if two satellites remain close, the line segments computed from ideal data are extremely short (less than 2 mm). Moreover, their length do not depend by the sampling time of the input data but depend only on the simulated dynamics of the target and chaser. On the other hand, the main source of error in the estimation of the location of the CoM is the value of the linear velocity of the target relative to the chaser. In particular the most influencing value is the magnitude of the component that is orthogonal to the instantaneous screw axis, whose direction is the same of the angular velocity vector. The last statement is confirmed by application examples. For instance figure 5 illustrates the comparison between the mentioned component of the velocity and the error in localizing the CoM by performing the pseudo-intersection on two different sets of ideal trajectories (no noise and no missing samples).

General Procrustes analysis to regularize noise on trajectories
The results that are shown in figure 5 refers to data that do not suffer the corruption provided by uncertainties due to the imperfections of the sensors used for the data acquisition. The same procedure applied to corrupted data produces unbounded error. The main reason of this sensibility of the procedure to noise is mainly attributed to an apparent loss of stiffness of the acquired shape of the target. Indeed, the distance between two generic points of a rigid body have to remain constant during the motion of the object. The randomness of the noise on the acquired trajectories causes instead a violation of this important property. Moreover, the presented solution to find screw axes is radically based on the nondeformability of the moving body.
Therefore, to maintain bounded the error it is absolutely needed to remove some of the noise while the distance between the acquired features is constrained to remain constant. However, to perform this filtering, the various poses should be comparable such to determinate one particular set of features, which is optimally mappable to all the poses through simple isometries, i.e. combinations of rotations and translations.
Given two different poses of an object, each characterized by a corrupted set of points outlining the surface of the object, a first operation that is useful to make comparisons is the overlap of the centroids of the two sets. The overlapping is obtained by simply subtracting the vector 2 that connects the centroids to all the points characterizing the second pose. If the positions of those points are the columns of a configuration matrix 2 , this operation corresponds to subtracting to 2 the matrix , whose columns are all equal to the vector 2 . Hence, after filtering-out the translational information, one has two configurations matrices 1 and 2 ′ , which should be separated by one rotation about the coincident centroid of the two sets of points. Actually, due to the corruption, there are no rotations that exactly maps the two matrices. Thus, the best possible operation to compare the poses consists in finding the orthonormal rotation matrix 2 that most closely maps 2 ′ to 1 This last statement correspond to the so-called Orthogonal Procrustes problem, whose formulation is the following: This problem is a convex constrained quadratic optimization whose closed-form solution is shown in [13]. The solution comes from the singular value decomposition of the matrix 1 2 ′ to find its nearest orthogonal matrix. In particular it holds: where Σ and are two orthogonal matrices that represent two following rotations, while is a diagonal matrix that represent an intermediate scaling. Thus, it holds: The same applies to the configuration matrices 3 ′ , 4 ′ …, ′ with respect to 1 for the entire duration of the time interval in which a certain set of features is detected. Once a new set is detected, the comparison procedure restart using a new reference configuration matrix (see figure 2).
Obviously, the transformed sets of homologous points do not matches exactly but they form a group of samples of a random variable, which is called shape. The mean of these samples can be considered as the optimal shape, which is represented by the resulting mean configuration matrix ̅ 1 .
The estimation of the original shape of a moving body from the detection of its natural features or landmarks takes the name of statistical shape analysis. In particular the described procedure is better known as general Procrustes analysis (GPA) [14]. Figure 6 illustrates this technique. Note that after the computation of the mean configuration matrix ̅ 1 , it is possible to perform the following operation: where = 2, 3, … , .
The new obtained positions of the features are such to preserve the values of the distances between the features itself. This regularization of the raw acquired data allows the efficient determination of the screw axes through the mentioned 3D generalization of the Reuleax's rule.

Estimation of the trajectory of the CoM of the target
Although the GPA has a significant positive effect on the capacities of the method exposed in section 2, the uncertainties on the determination of the screw axes remain non-negligible. An example of the results of such determination after the GPA of raw simulated data is in figure 7: the directions and the positions of the poles are plotted as signals in the time domain; moreover, they are compared to their reference values.
The typical regularity of the attitude dynamics of space debris leads one to understand that the irregular behaviour shown by the screw axes in figure 7 is due to residual uncertainties. A low-pass filter is then exploited to remove these residuals. Phase delays are not acceptable because they would provide bias to the final estimated trajectory of the CoM. Therefore, zero-phase filtering is considered.
The pseudo-intersection of the filtered screw axes provides, as shown in section 2, the estimated trajectory of the CoM of the target. However, as visible from the example in figure 8, the resulting signal contains a significant high-frequency noise. Moreover, the signal presents missing samples due to the introduced occlusions.
The final trajectory is recovered by polynomial curve-fitting. The error related to the example in figure 8 is shown in figure 9. This error has a non-random behaviour but it remains quite bounded. Considering the previous methods listed in the introduction, the one proposed in this paper showed a competitive accuracy. Notwithstanding, the level of required assumptions to let the method working is drastically lowered. In particular, there is no need to detect always the same features. Moreover, the proposed method is robust to occlusions or any kind of sudden loss of data.

Conclusions
The goal of this work was the investigation of methods and algorithms for localizing the CoM of a space debris during an active removal mission. Space debris mitigation is one of the most important topics that are currently discussed within the space community. An accurate knowledge of the CoM location helps the capture and the removal of the object significantly. Indeed, the contact between two noncooperative systems severely undermines the stability of the joint system. The underlying assumption under which the methods were developed regards the possibility of tracking the positions of several features of the passive object. However, part of the research work was spent to overcome one of the most common but restrictive assumptions in the literature, that is to consider the possibility of tracking one fixed set of features of the objects.
In this work, a pure kinematic estimation approach was developed. The approach is based on the computation of the finite screw axes that represent the motion between the various detected poses of a space debris. In particular, the combined usage of statistical shape analysis and low-pass filters was decisive for guaranteeing high accuracy in case of corrupted input data. The errors obtained, although distributed in a non-catalogued way, resulted bounded. The found accuracy is similar to one of the best methods known; however, the level of the necessary hypothesis was significantly lowered.