# Relative Pose and Inertia Determination of Unknown Satellite Using Monocular Vision

Authors: Danil Ivanov; Mikhail Ovchinnikov; Marianna Sakovich
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

The paper proposes a two-stage algorithm for autonomous relative motion determination of noncooperative and unknown object flying 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 first 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.


---

## Body

## 1. 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 angles-only visual navigation [4–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–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–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–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.


## 2. 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]:
(1)x¨+2ωż=0,y¨+ω2y=0,z¨−2ωẋ−3ω2z=0,where x, y, and z are coordinates of the target satellite with respect to the orbital reference frame O1xyz, origin O1 is placed in the chaser satellite center of mass which moves along a circular orbit of radius r0 and with orbital angular velocity ω, axis O1z is directed along chaser satellite radius-vector from the center of the Earth, O1y is a normal vector to the orbit plane, and O1x is the vector product of O1y and O1z (see Figure 1).Figure 1
Reference frames.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(2)Jω̇+ω×Jω=3ω2⋅Ae3×JAe3,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, e3 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.


## 3. 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. Consider a feature point Р1 of the target satellite which has a radius-vector R~p1=у1р1,у2р1,у3р1T in the frame О2у1у2у3. On the other hand, the point has coordinates Rp1=x1р1,x2р1,x3р1T in the frame Орх1х2х3. Figure 2 illustrates these two vector relations and equality ОрР1→=ОрО2→+О2Р1→.Figure 2
Camera reference frame.This vector expression can be written in the frameОрх1х2х3 as follows:
(3)Rp1=ОрО2→+A⋅R~p1,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
(4)х~1p=x1p1−x3p/f,х~2p=x2p1−x3p/f,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 (Figure3).Figure 3
Microsatellite “Chibis-M” and the reference features.


## 4. 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 pointsxi,yi on the first image are measured. On the next image, we will search the corresponding feature points in the certain ε-neighborhood of the xi,yi. 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 xi,yi, 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:(5)fr0,q0,v0,ω0,J,R1,R2,…,Rk,where r0 is a vector to the target satellite center of mass at the moment of the first image acquisition, v0 is an initial vector of relative velocity, q0 is a vector part of mutual attitude quaternion at the first moment, ω0 is an initial vector of angular velocity, J is a tensor of inertia of the target object, and Rk 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 ρitjx^p,itj,y^p,itj of ith reference point at time step tj are functions of the estimated vector f.Consider an error function:(6)Φf=∑tj∑ixitj−x^tjp,i2+yitj−y^tjp,i2,where xitj,yitj are measured and collected in database reference points coordinates and x^tjp,i,y^tjp,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:
(7)F=∂2Φ∂fj∂fk.If the rank of theF 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 Figure4. 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 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.Figure 4
Scheme of an 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.


## 5. 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 vectorf0 were chosen, next parameters are as follows:
(8)r0=000m,v0=0.000.000.01ms,q0=0.00.20.0,ω0=0.00.10.0s−1,J=diag1.61.161.86kg⋅m2.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. Table1 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 minimization problem. One can see that the errors of the components of the vector f steadily decrease as the number of iterations increases.Table 1
Database of the measurements of the reference points (in pixels).
Number of frames
16
17
18
19
20
21
22
23
x1
156,36
157,16
157,02
158,63
159,35
160,33
161,25
161,79
x2
188,20
186,49
186,35
184,46
182,41
182,05
181,65
179,91
x3
242,96
243,80
—
—
—
—
—
—
x4
327,80
327,86
328,55
—
—
—
—
—
x5
411,85
411,58
411,67
411,60
411,63
411,63
411,69
411,77
x6
419,66
420,05
420,11
420,05
419,93
420,22
419,98
420,21
y1
217,04
220,09
220,49
225,49
229,19
230,13
231,63
235,26
y2
80,32
80,26
80,57
80,51
81,79
82,13
82,38
82,95
y3
241,28
243,88
—
—
—
—
—
—
y4
271,32
271,31
271,59
—
—
—
—
—
y5
235,25
238,64
239,00
244,15
247,69
249,41
250,43
254,30
y6
78,71
78,63
78,66
78,93
79,07
79,41
79,47
80,00Figure 5
The difference between the estimated initial (a) center of mass position and (b) velocity vectors.
(a)
(b)Figure 6
Difference between diagonal elements of tensor of inertia.Figure 7
The difference between the estimated position of (a)R1 and (b) R2 reference points.
(a)
(b)


## 6. 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 root-mean-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
(9)Φ=∑ix~c,i−x^i2+y~c,i−y^i2,where x^i and y^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.
(10)x^i=x^ixc,yc,zc,q1,q2,q3,y^i=y^ixc,yc,zc,q1,q2,q3,where xc,yc,zc are the coordinates of the spacecraft center of mass and q1,q2,q3 are the vector part of the quaternion q=q0q1q2q3 that describes its attitude. Let f=xc,yc,zc,q1,q2,q3 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 parametersfi should be filtered. At each step, one has propagated value fi− obtained by integration of the relative motion equations. It is reasonable to use a constant gain Kalman filter to estimate the state vector fi+. 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
(11)fi+=fi−+Kfi−fi−,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.


## 7. 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”. Figure8 illustrates the satellite “Chibis-M” in the docking hatch of cargo vehicle “Progress-13M” during undocking from ISS.Figure 8
Microsatellite “Chibis-M” in the docking hatch of cargo vehicle “Progress-13M” during undocking from ISS (credit: SRI of RAS).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. Figure9 shows the direction of the “Progress-13M” frame axes.Figure 9
“Progress-13M” position relative to the Sun and the orbital frame during the undocking of the “Chibis-M.”Figure10 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.Figure 10
The frame from the video of “Chibis-M” separation from the cargo spaceship “Progress-13M.”Figure 11
Binary image from the video shot during the separation.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 (Figure12).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.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·m2. Figures 13 and 14 show the results of relative pose determination of the Chibis-M satellite at the second stage of the algorithm.Figure 13
The relative position of the “Chibis-M” center of mass.Figure 14
The attitude of the “Chibis-M.”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 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.


## 8. 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.



---

*Source: 9731512-2018-10-29.xml*