Relative Pose and Inertia Determination of Unknown Satellite Using Monocular Vision

. The paper proposes a two-stage algorithm for autonomous relative motion determination of noncooperative and unknown object ﬂ ying in space. The algorithm is based on image processing and can be applied to motion determination of space debris with unknown geometry and dynamic characteristics. The ﬁ rst stage of the algorithm is aimed at forming a database of possible reference points of the object during continuous observation. Tensor of inertia, initial velocity, and angular velocity of the object are also estimated. Then, these parameters are used in the second stage of the algorithm to determine the relative motion in real time. The algorithm is studied numerically and tested using the video of the Chibis-M microsatellite separation.


Introduction
A large amount of space debris flying around the Earth prevents intensive exploration of near-Earth orbits.Space community aims to remove the most dangerous fragments and inactive satellites.Many missions were proposed in the last decade.One of the most perspective concepts is to launch the special "chaser" satellite that will collect debris and move them to the low Earth orbit, so they will be rapidly deorbited by atmospheric drag.The most crucial and complicated part of this algorithm is the debris collection [1], but for the task solving it is necessary to determine the relative motion of the debris of unknown shape and dynamic characteristics.
In recent years, relative position and attitude estimation of spacecraft have been extensively studied, and many methods have been proposed.Chaser satellite must be equipped with the onboard sensor for relative pose determination, i.e., relative attitude and relative center of mass position.The most reliable sensor type is an optical one that usually comprises a CCD camera.In the paper [2], an overview of the existing visual-based systems is presented, and in the paper [3], the state of the art of the cooperative and uncooperative satellite pose determination techniques are presented.A set of papers are devoted to the developing and flight testing of the algorithms for relative center of mass position determination of the unknown target using anglesonly visual navigation [4][5][6].The relative pose navigation system based on image processing was implemented in several projects; for example, PRISMA [7,8], it is used during docking stage of the "Progress" to ISS [9] and also of the cargo spacecraft ATV [10][11][12].In the paper [13], a set of the relative initial pose estimation techniques are compared.A number of papers are devoted to the developing of image-based navigation for satellites using laboratory facility [14][15][16][17][18][19].In the papers mentioned above, the motion determination algorithms are usually based on recognition of specially designed satellites' marks which positions are known in the target-related reference frame.Moreover, the dynamic characteristics such as tensor of inertia or position of the center of mass are used for target satellite motion prediction, i.e., the 3D model of the target is known.In the case of the unknown target before the pose determination, a prior monitoring stage should proceed and the model of the target is estimated.The well-known concept of the simultaneous localization and mapping (SLAM) can be adopted for these tasks [20].A set of papers are devoted to the problem using images obtained by multiple cameras installed onboard of the chaser satellite [21][22][23].For example, a SLAM approach based on the extended Kalman filter and inertia tensor determination using stereovision measurements is proposed in the paper [22].A feature-based SLAM algorithm using the stereovision system is developed in [23].There are a very few papers on the pose determination of unknown object using monocular vision measurements.In the paper [24], a purely monocular-based vision SLAM algorithm for pose tracking and shape reconstruction of an unknown target is presented; however, it is assumed that the target rotates with constant angular velocity and the target inertia tensor is not estimated.The fusion of the data obtained from the monocular camera and a range sensor is utilized in the algorithm for the pose, and shape reconstruction of a noncooperative spacecraft is considered in the paper [25]; the unscented Kalman filter is used for the estimation of 6n + 17 variables, where n is the number of the feature points.The current paper is quite similar to [25], but we divided the pose determination algorithm into two stages-the most computational power consumed determination of the reference point position and inertia tensor using least square method and the second, determination of the relative pose only based on the data obtained by the first stage.So, it is not necessary to estimate the reference points and the inertia tensor all the time since it is assumed that they are constant since the unknown satellite is considered as a rigid body.
In this paper, a noncooperative and unknown object relative motion determination algorithm based on the image processing of satellite illuminated by the Sun is considered.Initially, the target satellite has unknown geometry and mass distribution.It is assumed that the object has some features (surfaces, lines, or points) which can be recognized in the picture.The algorithm consists of two stages.First, the reference features' positions in target satellite reference frame are estimated during continuous observation.Also, the tensor of inertia is determined using video information.During this stage, the database of observable reference features is formed.These features will be used later for relative motion estimation in real time.The pose vector of a relative state that consists of the center of mass position, linear velocity, attitude quaternion, and angular velocity vectors is estimated.The main idea of the first stage of the algorithm is to find such dynamical parameters and positions of reference points that minimize the difference between the points' pixel coordinates on frames and the estimated by model ones.
The second stage of the algorithm is briefly described as follows.For a given initial relative state vector, the satellite reference point positions are projected to the video frame.Reference points' projections are searched then in the neighborhood of these positions.Found reference points' positions are inaccurate due to a noise and discrete sensing element.Least-square method (LSM) is used to minimize the sum of squared distances between the predicted and obtained reference points' positions in the frame.In such a way, the relative position and attitude are roughly estimated.This estimation is used as the Kalman filter input to compute the state vector with higher precision.In the paper, the relative motion determination accuracy is investigated depending on the parameters of the objects.The algorithm is tested using the video of the Chibis-M microsatellite separation.

Equations of Motion
Consider two satellites, a chaser and a target, flying along the circular LEO.Let the relative translational motion be described by well-known Clohessy-Wiltshire equations [26]: where x, y, and z are coordinates of the target satellite with respect to the orbital reference frame O 1 xyz, origin O 1 is placed in the chaser satellite center of mass which moves along a circular orbit of radius r 0 and with orbital angular velocity ω, axis O 1 z is directed along chaser satellite radiusvector from the center of the Earth, O 1 y is a normal vector to the orbit plane, and O 1 x is the vector product of O 1 y and O 1 z (see Figure 1).
Consider the attitude motion of both satellites.If there is no control torque, then after neglecting perturbations except gravitational torque, the attitude motion could be described by where ω is an absolute angular velocity of the target spacecraft, J is its inertia tensor, A is a direction cosine matrix that is used for the target spacecraft frame rotation relative to the orbital frame, e 3 is a unit vector directed from the Earth to the satellite center of mass written in the frame О 2 у 1 у 2 у 3 , and О 2 is a target satellite center of mass.

Measurement Model
Consider two satellites move so close to each other that details of the sunlit spacecraft structures can be recognized.
The coordinate system О 1 х 1 х 2 х 3 is the chaser satellite body-fixed frame, and the system О 2 у 1 у 2 у 3 is the target body-fixed frame (Figure 1).Let the CCD camera's optical axis be parallel to the Ох 3 axis.If we shift the coordinate frame along the О 1 х 3 axis and take the center О р of the CCD matrix as the center of the coordinate frame, we get a coordinate frame of the camera (Figure 2).Suppose the distance О 1 О р is known. 2 International Journal of Aerospace Engineering Consider a feature point Р 1 of the target satellite which has a radius-vector T in the frame О 2 у 1 у 2 у 3 .On the other hand, the point has coordinates T in the frame О р х 1 х 2 х 3 .Figure 2 illustrates these two vector relations and equality This vector expression can be written in the frame О р х 1 х 2 х 3 as follows: where A is a transformation matrix from the frame О 2 у 1 у 2 у 3 to the frame О р х 1 х 2 х 3 .Vector О р О 2 and matrix A are relative position and relative attitude of the two satellites correspondingly.
Consider an image of the reference feature on the CCD matrix.According to projective geometry, the point Р 1 has the coordinates in the picture defined by where f is the camera focus distance with negative sign.The coordinates in image (in pixels) can be measured.The measurements of the reference features can be processed to estimate relative position and relative orientation.
As a reference, features of the target object can be considered special surfaces of unique shape and illumination (for example, solar panels), lines, or special points.Consider for concretization the tips of the antennas as reference features (Figure 3).

First Stage: Reference Features and Motion Parameter Determination Algorithm
Let the target satellite be uncooperative with unknown reference point positions (R p1 in (3)).Assume that the chaser satellite observes the target during continuous time.In that time, the reference features or its candidates can be tracked on the CCD images.The new feature candidates can arise in the field of view or disappear in the shadow during the observation.As a result, a sequence of the measured coordinates of the features is available for processing.The problem is to estimate the positions of the reference features using these measurements.It can be done by LSM using the equations of motion if initial conditions and tensor of inertia of the target satellite are known.If the parameters are unknown, it is possible to estimate them jointly with the reference point positions.First, consider how to track the reference feature position from image to image.When the first image of the target is obtained, it is reasonable to filter it for noise elimination.Then, the image can be converted to the binary one by using predetermined threshold.On the binary image, it is easy to find boundaries of the detectable reference feature positions.If the feature is a point, then the center of the mark can be calculated as a center of illumination or center of the boundary.Assume all positions of the reference points x i , y i on the first image are measured.On the next image, we will search the corresponding feature points in the certain ε -neighborhood of the x i , y i .If the feature point is found in ε-neighborhood of the previous point, then assume that it is the same features and collect it in special time column with coordinates of the ith feature points on the sequence of frames.If there is no feature point in ε-neighborhood of the x i , y i , assume that the point is shadowed or it is a point with not stable luminosity properties and it cannot be a reference point.Such an algorithm of feature point association fails if the target object moves or rotates too fast and it is applicable only in the close range relative motion when all the details of the target surface can be recognized in the image.
As a result of the feature point association algorithm, the database of the possible reference points is obtained.It is reasonable to process the feature points with the lengthiest records.It means that some of the point was observed for a rather long time and it is the points with stable luminosity.
Consider a vector of estimated parameters: where r 0 is a vector to the target satellite center of mass at the moment of the first image acquisition, v 0 is an initial vector of relative velocity, q 0 is a vector part of mutual attitude 3 International Journal of Aerospace Engineering quaternion at the first moment, ω 0 is an initial vector of angular velocity, J is a tensor of inertia of the target object, and R k is the radius-vector of the kth reference point in the body-fixed frame.If the vector f is defined, the measured coordinates of the reference points can be calculated by integrating equations of motion ( 1) and ( 2) and using measurement model ( 3) and (4).So, the propagated measurement vectors ρ t j i xt j p,i , ŷt j p,i of ith reference point at time step t j are functions of the estimated vector f.
Consider an error function: where x i t j , y i t j are measured and collected in database reference points coordinates and xt j p,i , ŷt j p,i are propagated values of corresponding measurements.The problem of minimization of the function Φ f can be solved only numerically.When one solves the task of minimization of ( 6), the problem of the observability of the vector f must be discussed.In the paper [27], the covariance matrix nonsingularity method was applied to the observability problem of pose estimation using line-of-sight measurements.The main idea of the method is to construct the Fisher information matrix and to define its rank.The Fisher matrix is defined as follows: If the rank of the F is equal to the length of vector f, then f is observable, i.e., it is possible to achieve vector f from measurements.Numerically, it has been studied that if there is more than three observable reference points, it is possible to estimate vector f.
The scheme of proposed algorithm is presented in Figure 4.A video with continuous target satellite observation is required for the algorithm.Then, the separate frames are extracted from the video.The frames are processed to have a binary contour image, where the centers of the feature points are calculated like described above.Its pixel coordinates are collected in the feature point database.When all International Journal of Aerospace Engineering the frames are processed, the feature point database is filtered to eliminate all the noise data like a shortly observed feature points.Finally, the resulting database is used for solving the minimization problem (6).So, the estimated feature point positions, tensor of inertia, and initial conditions are the result of the algorithm.
The proposed algorithm has a set of restrictions.The satellite on the image can be recognized only if the background is dark and direct sunlight does not saturate the image.The temporary occlusion of the reference points can occur during the motion.However, the first stage of the algorithm uses a continues sequence of the images for reconstructing the 3D model of the reference point position, so occlusion of the reference points might affect the provided accuracy but this will not lead to the confusion of the position of the object.

Numerical Simulation
Consider a numerical example of the algorithm.Let's simulate the relative motion of the microsatellite Chibis-M during its separation from the cargo-vehicle Progress [28].Let the tips of the antennas be the reference points.
As initial values of the estimated vector f 0 were chosen, next parameters are as follows: The initial value of the reference point coordinates is random.The algorithm processed a sequence of the video frames.After applying the reference feature association algorithm, the database of the measurements is constructed.Table 1 demonstrates a part of the database, where coordinates (in pixels) are collected for each reference point.As one can see at frames 18 and 19, two points disappeared since it has been shadowed.The database measurements were used for estimation of the parameter vector f.Figures 5-7 demonstrate the convergence process of the solution of the  5 International Journal of Aerospace Engineering minimization problem.One can see that the errors of the components of the vector f steadily decrease as the number of iterations increases.

The Second Stage: Relative Pose Determination Algorithm
When the reference point positions are determined, the problem of real time motion estimation can be solved [28].Due to the estimation errors and noise in the CCD matrix, obtained reference point positions have error.So, a rootmean-square method should be applied.After determination, at least of three reference points on image, evaluate the simulated coordinate corrections using least squares procedure.Define an error functional as where xi and ŷi are the simulated coordinates of the reference marks and x c,i and y c,i are the centers of the objects in the picture in the ε-neighborhood of the reference marks.
According to (3), point's position in the frame depends on the position of the spacecraft center of mass.
xi = xi x c , y c , z c , q 1 , q 2 , q 3 , ŷi = ŷi x c , y c , z c , q 1 , q 2 , q 3 , 10 where x c , y c , z c are the coordinates of the spacecraft center of mass and q 1 , q 2 , q 3 are the vector part of the quaternion q = q 0 q 1 q 2 q 3 that describes its attitude.Let f = x c , y c , z c , q 1 , q 2 , q 3 denote the vector of parameters.
The problem of minimization of the function Φ f can be solved only numerically.
Because of the noise in the reference point measurements, the resulting vector of parameters f i should be filtered.At each step, one has propagated value f − i obtained by integration of the relative motion equations.It is reasonable to use a constant gain Kalman filter to estimate the state vector f + i Constant gain Kalman filter does not require too much computational power and can be easily implemented on-board.So, the final estimated value can be calculated as where K is a constant gain matrix.Note that the developed approach of the relative position determination has some restrictions.For example, it is necessary to find, at least, four reference marks that are not uniplanar to determine spacecraft position explicitly.In addition, the algorithm cannot operate if the Sun lights the camera or the second spacecraft is in the eclipse.If the second spacecraft has a symmetric form, there is an ambiguity when solving (3).The shooting frequency should be higher than the rotational velocity of the spacecraft.6 International Journal of Aerospace Engineering

Algorithm Validation Using the "Chibis-M" Launch Video
The developed algorithm was applied for video processing of the "Chibis-M" separation from the cargo vehicle "Progress-13M".Figure 8 illustrates the satellite "Chibis-M" in the docking hatch of cargo vehicle "Progress-13M" during undocking from ISS.After undocking, the cargo vehicle was in a circular orbit at the altitude about 600 km, and the "Chibis-M" was separated on the 25th of January, 2012 at 2:18:00 UTC. Figure 9 shows the direction of the "Progress-13M" frame axes.
Figure 10 illustrates the frame from the video of the "Chibis-M" separation.The video was ruined by sunlight, but the long antennas are well visible.So, it is possible to recognize the frame, and Figure 11 shows a binary image of the frame.
Next, the binary image is used as input for the algorithm described above, and the current relative phase state is defined using the least squares procedure (Figure 12).
Unfortunately, due to the overexposure of the video, the reference point on the images could be recognized only during three seconds (it is about 100 frames).The obtained frames were processed by the first stage of the proposed algorithm for the reference point database creation, and for its position estimation, along with initial condition of the relative motion equations and inertia tensor determination.As a result, the reference point position is estimated with the 3 cm accuracy.An error in the estimated principal moments of the inertia tensor is about 0.3 kg•m 2 .Figures 13 and 14 show the results of relative pose determination of the Chibis-M satellite at the second stage of the algorithm.
The root-mean-square error of determination is about 1 cm/s for the linear velocity and is about 0.5 deg/s for Y component of the angular velocity and 4 deg/s for its X component and Z component (not a high determination accuracy  7 International Journal of Aerospace Engineering of the angular velocity components is specified by the noise on the reference mark images).The value of the angular velocity is equal to 4.24 deg/s.In 30 seconds after the separation, the navigation system was actuated; the magnitude of the angular velocity was measured by a rate-of-turn gyroscope and equaled 4.5 deg/s.Thus, the angular velocity evaluated by the algorithm proved to be close to the angular velocity measured by the rate-of-turn gyroscope.It indicates a satisfactory precision of the algorithm.

Conclusion
The proposed monocular vision-based two-stage relative pose determination algorithm is a convenient instrument for relative navigation in case of operations with noncooperative and unknown target satellite.It can be also useful for space removal missions when the space debris has an unknown shape.The algorithm has a number of limitation: it can be applied only for short relative distances, the target satellite must have some well-recognizing features, and the satellites' relative motion could not be too fast.Nevertheless, the numerical simulation shows that the algorithm is stable, after the feature point position determination, it estimates the relative motion in real time.

Figure 4 :
Figure 4: Scheme of an algorithm.

Figure 5 :
Figure 5: The difference between the estimated initial (a) center of mass position and (b) velocity vectors.

Figure 6 :
Figure 6: Difference between diagonal elements of tensor of inertia.

Figure 7 :
Figure 7: The difference between the estimated position of (a) R 1 and (b) R 2 reference points.

Figure 9 :Figure 10 :Figure 11 :
Figure 9: "Progress-13M" position relative to the Sun and the orbital frame during the undocking of the "Chibis-M."

Figure 12 :
Figure 12: The process of relative position determination.The asterisks are the pictures of the initial approximation of the reference marks; daggers are the reference mark estimation.

Figure 13 :
Figure 13: The relative position of the "Chibis-M" center of mass.

Table 1 :
Database of the measurements of the reference points (in pixels).