Parallel vision-based pose estimation for non-cooperative spacecraft

Abstract

This article proposes a relative pose estimation method between non-cooperative spacecrafts based on parallel binocular vision. As the information of non-cooperative spacecraft in space is not accessible, the target is considered to be freely tumbling in space. The line feature of non-cooperative target is used to extract the feature points first; then the stereo matching and three-dimensional restructuring are taken for the feature points; finally, an algorithm based on parallel binocular vision algorithm is used to calculate the relative pose between the target coordinate and the world coordinate. The experimental results show that the proposed method has high-accuracy real-time performance.

Introduction

With the development of technology, more than 80–130 satellites are launched into space every year. The fuel running out, system failing, and other problems may lead to satellite failure. ¹ To explore the feasibility of repairing the malfunctioned satellite on orbit, Germany, America, and other countries have researched the capture method for the malfunctioned satellites, such as DEOS project and MiTEx system. The research on the three-dimensional (3D) pose estimation is one of the key technologies to ensure the security and reliability of the malfunctioned satellite rendezvous and docking on orbit. ² The use of vision systems for autonomous landings or pose estimation of aircrafts and spacecrafts has been earlier suggested. ³ An important task of computer vision is to measure 3D pose parameters of one object. When the measurement system comprises two or more cameras, 3D pose parameters can be obtained by intersecting the cameras at three or more feature points on the object. ⁴

A number of pose estimation methods have been reported. Segal et al.^5,6 developed a stereovision-based filtering algorithm for estimating the state between two non-cooperative spacecrafts. Using two cameras mounted on a chaser satellite, the relative state with respect to a target satellite, including the position, attitude, and rotational and translational velocities, is estimated. Du et al. ⁷ proposed a method based on two collaborative cameras to determine the pose of a large non-cooperative target. The corresponding algorithm of image processing and pose measurement is addressed. Chen et al. ⁸ and Yu et al. ⁹ established a binocular stereovision model and proposed a non-cooperative target relative pose measuring method based on stereovision and RANSAC algorithm. Liu et al. ¹⁰ proposed a novel method based on multi-feature information fusion to estimate the non-cooperative spacecraft. Xu et al. ¹¹ achieved the target recognition and pose measurement based on the stereovision; the guidance, navigation and control (GNC) of the chaser; and the coordinated plan and control (CP&C) of space robot. According to the most rigid objects’ feature on the ground with the mutually perpendicular or nearly perpendicular planar structures, Lv et al. ¹² adopted the normals of planar areas on the target surface to estimate the 3D pose of the target with ladar range image. Sultan et al. ¹³ proposed a monocular camera vision system for a 6-degree of freedom (DOF) drawing robotic arm to estimate 3D pose of the end effecter robustly. Kirac et al. ¹⁴ had demonstrated an implementation of regression forests for estimation of the articulated 3D pose of the human hand. Yu et al. ¹⁵ proposed an estimation method of relative pose based on stereovision and Kalman filter, which consists of three parts: body frame reestablishment, translational parameter estimation, and rotational parameter estimation. It is presented for the final phase of the rendezvous and docking of non-cooperative satellites.

Compared with the monocular vision model, the binocular one possesses many merits: the binocular cameras can be used for the data fusion in the stereovision to enhance the robustness of the measuring precision and algorithm; it is free to any designs of the marker, so there is no strict geometric restriction between markers; the binocular cameras can backup information from each other to increase the reliability of measuring results.

In view of the non-cooperative space target, this article proposes a relative pose estimation between spacecrafts. On the basis of the parallel optical axis, this algorithm is relieved from the external calibration of the camera, which is very convenient to the on-orbit service. Taking full advantage of the target line characteristics and extracting the intersection of the associated lines as the feature points, the feature point extraction based on the geometric line feature can reduce the workload of the general epipolar matching. Moreover, the accuracy and stability of the feature points’ extraction are higher in theory.

As shown in Figure 1, we can see that the line feature of the aircraft is very obvious. The tripod to connect the solar array and the aircraft body is a weak part of the whole structure. Once damaged, the tripod must influence the on-orbit operation of the spacecraft. Relative to the other features of the spacecraft, the tripod is a simple structure and an easy identification. So the tripod of spacecraft is selected as the research object. The obtained relative position and attitude of the tripod can be served as the basis for the further relative pose estimation of the whole spacecraft.

Figure 1.

Spacecraft model.

The flow chart of pose estimation for non-cooperative spacecraft is as shown in Figure 2.

Figure 2.

The flow chart of pose estimation for non-cooperative spacecraft.

This article is organized as follows. How to extract the feature point is introduced in section “Feature extraction.” Section “3D relative pose estimation” represents the 3D pose estimation algorithm for space non-cooperative target. In section “Experimental validations,” the semi-physical simulation experiment is conducted to verify the proposed algorithm. The last section summarizes the whole article and gives the conclusion.

Feature extraction

The line feature is the basis of the proposed pose estimation algorithm. We use the line feature extraction to obtain the feature points that are very important to computing the pose parameters. First of all, the high-precision boundary points are extracted by the boundary tracking method. Second, the boundary lines are obtained by the least squares fitting of the boundary points. Finally, the associated line intersections are extracted as the feature points.

Edge feature extraction

Hypothesis that the image area element consists of four adjacent pixels, f(m, m + 1), f(m + 1, n + 1), f(m, n), f(m + 1, n). Use plane ax + by + c to fit four adjacent pixels of the image area element ΔS and use ax + by + c to fit f(x, y); the mean square deviation is

$$\begin{array}{cc}\multicolumn{1}{c}{\epsilon }& ={\displaystyle \sum _{x,y\in \Delta s}}{[g(x,y)-f(x,y)]}^2\hfill \\ \multicolumn{1}{c}{}& ={[\mathrm{am}+\mathrm{bn}+c-f(m,n)]}^2\hfill \\ \multicolumn{1}{c}{}& +{[a(m+1)+\mathrm{bn}+c-f(m+1,n)]}^2\hfill \\ \multicolumn{1}{c}{}& +{[\mathrm{am}+b(n+1)+c-f(m,n+1)]}^2\hfill \\ \multicolumn{1}{c}{}& +{[a(m+1)+b(n+1)+c-f(m+1,n+1)]}^2\hfill \end{array}$$

(1)

In order to achieve the most optimum anastomosis, the mean square deviation needs to be minimized. The a, b, and c in equation (1) are, respectively, derivations. There are

$$\begin{array}{c}\hfill a=\frac{f(m+1,n)+f(m+1,n+1)}{2}-\frac{f(m,n)+f(m,n+1)}{2}\hfill \\ \multicolumn{1}{c}{b=\frac{f(m,n+1)+f(m+1,n+1)}{2}-\frac{f(m,n)+f(m+1,n)}{2}}\\ \multicolumn{1}{c}{c=\frac{1}{4}[3f(m,n)+f(m+1,n)+f(m,n+1)-f(m+1,n+1)-\mathrm{ma}-\mathrm{nb}]}\end{array}$$

(2)

Based on the definition of gradient, the gradient magnitude of the plane ax + by + c is

$$G[g(x,y)]={[{\left(\frac{{\partial }_g}{{\partial }_x}\right)}^2+{\left(\frac{{\partial }_g}{{\partial }_y}\right)}^2]}^{1/2}={(a^2+b^2)}^{1/2}$$

(3)

It follows that a is the difference of the average value between two columns and b is the difference of the average value between two rows. Since the difference is established on the basis of smooth, this method is lower sensitive to the noise compared with the differential operator.

Because the plane is the best approximation of image gray level in the neighborhood of 2 × 2, a plane gradient can be reasonably regarded as the approximation of image gradient values at the neighborhood center point (m + 1/2, n + 1/2). The idea of edge detection based on the surface fitting is to use a plane or curved surface to approximate an image area element. Then, this point gradient is replaced by the plane gradient to implement the edge detection.

According to the principle of the edge feature extraction, the results are as shown in Figure 3.

Figure 3.

Edge feature extraction: (a) left image and (b) right image.

Discrete line feature extraction

On the basis of the edge feature’s extraction in “Edge feature extraction,” the discrete line feature is extracted. The starting point for each line can be served as the basis of the line feature extraction.

The criterion for the selection of edge feature can be shown as in Figure 4. The boundary starting points need to be in the eight connected regions. Meanwhile, when turning around the point in clockwise or anticlockwise rotation, there must be three adjoining pixels whose gray value is 0. The pixels marked by “1” are available and the pixels marked by “2” are not available. Following the edge line to track the edge points, the coordinate sequence of the edge points should be given, which is the basic requirement of edge tracking.

Figure 4.

Starting point.

The searching of the edge points’ step is as follows:

• Searching the untracked edge points by image-scanning in clockwise direction. The point, if detected, can be defined as the starting point of tracking, recorded as P _0.

• Searching the eight adjacent pixel points in a counterclockwise direction, it will turn up to step (c) if the change from the white pixel point to the red pixel point occurs in the position of the next edge point. Otherwise, the starting point tracked would be an isolated point if there is no red pixel in the eight adjacent pixel points, thus the tracking should be over.

• Moving forward to the next edge point P_n . If P_n−1  = P_n , P_n  = P ₁, then the tracking should be over.

The results of the discrete line feature extraction are as shown in Figure 5.

Figure 5.

Discrete line feature extraction.

Line feature serialization

A function relationship between x and y by the linear equations y = a ₀ + a ₁ x. a ₀ represents the intercept and a ₁ represents the slope. For group N, data precision obtained by measurement of (x_i , y_i ), i = 1, 2, ..., N, x_i value is believed to be accurate, and the error can be linked with y_i .

When we use the method of the least squares to estimate the parameters, the y_i ’s deviation and sum of squares are to be required for the minimum. The beeline fitting precision of observation is given by

$${\displaystyle \sum _{i=1}^N}{[y_1-(a_0+a_1{x}_1)]}^2|{}_{a=\hat{a}}$$

(4)

Based on the requirements, the formula should be

$$\begin{array}{cc}\multicolumn{1}{c}{}& \frac{\partial }{\partial a_0}{\displaystyle \sum _{i=1}^N}{[y_i-(a_0+a_1{x}_i)]}^2|{}_{a=\hat{a}}\hfill \\ \multicolumn{1}{c}{}& =-2{\displaystyle \sum _{i=1}^N}(y_i-{\hat{a}}_0-{\hat{a}}_1{x}_i)=0\hfill \\ \multicolumn{1}{c}{}& \frac{\partial }{\partial a_1}{\displaystyle \sum _{i=1}^N}{[y_i-(a_0+a_1{x}_i)]}^2|{}_{a=\hat{a}}\hfill \\ \multicolumn{1}{c}{}& =-2{\displaystyle \sum _{i=1}^N}(y_i-{\hat{a}}_0-{\hat{a}}_1{x}_i)=0\hfill \end{array}$$

(5)

After finishing getting normal equations

$$\{\begin{array}{cc}\multicolumn{1}{c}{}& {\hat{a}}_{0}N+{\hat{a}}_1\sum x_i=\sum y_i\hfill \\ \multicolumn{1}{c}{}& {\hat{a}}_0\sum x_i+{\hat{a}}_1\sum x_i^2=\sum x_i{y}_i\hfill \end{array}$$

(6)

The equations can be obtained from the best estimation of linear parameters a ₀ and a ₁

$$\begin{array}{cc}\multicolumn{1}{c}{{\hat{a}}_0}& =\frac{(\sum x_i^2)(\sum y_i)-(\sum x_i)(\sum x_i{y}_i)}{N(\sum x_i^2)-{(\sum x_i)}^2}\hfill \\ \multicolumn{1}{c}{{\hat{a}}_1}& =\frac{N(\sum x_i{y}_i)-(\sum x_i)(\sum y_i)}{N(\sum x_i^2)-{(\sum x_i)}^2}\hfill \end{array}$$

(7)

The edge features are the discrete feature points on the boundary. In order to get the edge feature lines, we need the discrete feature points to fit the least squares.

Based on the least squares fitting principle, the discrete feature points are selected to do least squares fitting. The result is shown in Figure 6.

Figure 6.

Edge feature.

Extraction and matching of feature points

Based on the least squares fitting, we obtain three edge feature lines on the same plane. By solving three linear equations, we can calculate three intersections of simultaneous pixel coordinates. The results are shown in Figure 7.

Figure 7.

Feature points’ extraction.

The epipolar constraint can greatly limit the search range of the corresponding point, especially when the binocular stereovision system is in a parallel alignment configuration. In the system, the x axes of two cameras are the same, so the polar line of image plane becomes parallel lines. Each epipolar equation does not need to be calculated when matching; we only need to scan the image matching the feature points. The flow chart of the extraction and matching of the feature points is as shown in Figure 8.

Figure 8.

The flow chart of the extraction and matching.

3D relative pose estimation

As shown in Figure 9, O_t-X_tY_tZ_t is the target coordinate and O_W-X_WY_WZ_W is the world coordinate. The relative pose refers to the position and attitude compared to the reference object. The transformation between the two coordinates can be represented mathematically. Therefore, if we know the coordinate and three non-collinear feature points in two coordinates, we can obtain the relative position and orientation. According to the relationship among various points between object coordinate and world coordinate, we can calculate the relative pose between world coordinate and target coordinate.

Figure 9.

Relative coordinate.

3D position in world coordinate

After the completion of the feature matching, the positions of the feature points are rebuilt in the observation coordinates according to the measuring principle based on the stereovision. Thus, 3D structure of the non-cooperation target in space can be reconstructed.

The transformation from world coordinate system to the camera coordinate system: since the world coordinate system is a coordinate system used to describe the illusion of camera position, we can take the binocular vision of the left image to the camera coordinate system as the world coordinate system. In this case, the transformation from the world coordinate in the parallel binocular vision system to the left camera coordinate can be expressed as follows

$$\left[\begin{array}{c}\hfill x_c\hfill \\ \multicolumn{1}{c}{y_c}\\ \multicolumn{1}{c}{z_c}\\ \multicolumn{1}{c}{1}\end{array}\right]=\left[\begin{array}{cc}\hfill R\hfill & \hfill T\hfill \\ \multicolumn{1}{c}{0^T}& \hfill 1\hfill \end{array}\right]\left[\begin{array}{c}\hfill x_w\hfill \\ \multicolumn{1}{c}{y_w}\\ \multicolumn{1}{c}{z_w}\\ \multicolumn{1}{c}{1}\end{array}\right]$$

(8)

where R is a rotation matrix, and there is no rotation relationship because of the parallel binocular vision of coordinates system. R is an identity matrix, T is a translation vector, and the mold is the base distance b. Because it is in right translation, the calculation should take a negative sign. 0 ^T is the transverse vector of 0.

The transformation from the camera coordinate to the imaging plane coordinate system: because the model is the pinhole imaging model, the transformation from the camera coordinate to imaging plane coordinate system agrees to the perspective projection, which can be expressed in homogeneous coordinates as follows

$$Z_c\left[\begin{array}{c}\hfill x\hfill \\ \multicolumn{1}{c}{y}\\ \multicolumn{1}{c}{1}\end{array}\right]=\left[\begin{array}{c}\hfill -f\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{-f}\\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]\left[\begin{array}{c}\hfill x_c\hfill \\ \multicolumn{1}{c}{y_c}\\ \multicolumn{1}{c}{z_c}\\ \multicolumn{1}{c}{1}\end{array}\right]$$

(9)

In which, f represents the focal length of the camera, (x_c , y_c , z_c ) is the arbitrary P 3D coordinate in camera coordinate system space, and (x, y) is the two-dimensional (2D) coordinate of P in imaging plane coordinate system.

The convention from imaging plane coordinate system to computer image

$$\left[\begin{array}{c}\hfill u\hfill \\ \multicolumn{1}{c}{v}\\ \multicolumn{1}{c}{1}\end{array}\right]=\left[\begin{array}{ccc}\hfill 1/\mathrm{dx}\hfill & \hfill 0\hfill & \hfill u_0\hfill \\ \multicolumn{1}{c}{0}& \hfill 1/\mathrm{dv}\hfill & \hfill v_0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill 1\hfill \end{array}\right]\left[\begin{array}{c}\hfill x\hfill \\ \multicolumn{1}{c}{y}\\ \multicolumn{1}{c}{1}\end{array}\right]$$

(10)

where (u ₀, v ₀) is an imaging plane coordinate system origin in the computer image coordinate system coordinates, often taken as the center point of a computer image coordinates. (d_x , d_y ) are the distance between the pixels adjacent to the imaging plane coordinate system in the x- and y-directions.

Substituting formulas (8) and (9) into formula (10), we can get equation (11)

$$z_c\left[\begin{array}{c}\hfill u\hfill \\ \multicolumn{1}{c}{v}\\ \multicolumn{1}{c}{1}\end{array}\right]=\left[\begin{array}{cccc}\hfill -f/\mathrm{dx}\hfill & \hfill 0\hfill & \hfill u_o\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill -f/\mathrm{dy}\hfill & \hfill v_o\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{cc}\hfill R\hfill & \hfill T\hfill \\ \multicolumn{1}{c}{0^T}& \hfill 1\hfill \end{array}\right]\left[\begin{array}{c}\hfill x_w\hfill \\ \multicolumn{1}{c}{y_w}\\ \multicolumn{1}{c}{z_w}\\ \multicolumn{1}{c}{1}\end{array}\right]$$

(11)

Setting equivalent transformation as follows

$$M_1=\left[\begin{array}{c}\hfill -f/\mathrm{dx}\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{-f/\mathrm{dx}}\\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\hfill u_0\hfill \\ \multicolumn{1}{c}{v_0}\\ \multicolumn{1}{c}{1}\end{array}\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right]$$

(12)

$$M_2\left[\begin{array}{cc}\hfill R\hfill & \hfill T\hfill \\ \multicolumn{1}{c}{0^T}& \hfill 1\hfill \end{array}\right]$$

(13)

Equation (11) can be abbreviated as follows

$$Z_c\left[\begin{array}{c}\hfill u\hfill \\ \multicolumn{1}{c}{v}\\ \multicolumn{1}{c}{1}\end{array}\right]=M_1{M}_2\left[\begin{array}{c}\hfill X_w\hfill \\ \multicolumn{1}{c}{Y_w}\\ \multicolumn{1}{c}{Z_w}\\ \multicolumn{1}{c}{1}\end{array}\right]=M\left[\begin{array}{c}\hfill X_w\hfill \\ \multicolumn{1}{c}{Y_w}\\ \multicolumn{1}{c}{Z_w}\\ \multicolumn{1}{c}{1}\end{array}\right]$$

(14)

wherein R is a rotation matrix, T is a translation vector, M ₁ is the camera internal matrix, M ₂ is the camera external matrix, and M = M ₁ × M ₂ is the camera projection matrix. In the parallel binocular vision system, both cameras are parallel movement and have identified the left camera coordinate system for the world coordinate system, thus the rotation matrix R is a unit matrix and the mold of translation vector is the distance b between the two cameras. Equation (11) can be further simplified into formula (15)

$$\{\begin{array}{c}\hfill Z_w·u_1=\frac{-f·x_w}{d_x}+\frac{f·b}{d_x}+u_0·Z_w\hfill \\ \multicolumn{1}{c}{Z_w·u_1=\frac{-f·y_w}{d_y}+v_0·Z_w}\end{array}$$

(15)

Finishing the above formulas’ calculation between the world coordinate system and the image coordinate system, we can obtain equation (16)

$$\{\begin{array}{c}\hfill x_w=\frac{-b·(u_1-u_0)}{u_2-u_1}\hfill \\ \multicolumn{1}{c}{y_w=\frac{-b·d_y(v_1-v_0)}{d_x·(u_2-u_1)}}\\ \multicolumn{1}{c}{z_w=\frac{b·f}{d_x·(u_2-u_1)}}\end{array}$$

(16)

According to equation (16), all the point’s 3D coordinates in world coordinate can be obtained.

Non-cooperative target pose solver

Setting P ₂ as the original point, the plane where P ₁, P ₂, P ₃ lies is the x–y plane, X ₂ X ₃ is the x-axis, y is perpendicular to P ₂ P ₃, the y-axis positive direction is the same as the world coordinate system, and z-axis is perpendicular to the x–y plane establishing the standard coordinate system, as shown in Figure 10.

Figure 10.

The relative coordinate system.

The relationship between the target coordinate and the world coordinate system is

$$\left[\begin{array}{c}\hfill x_w\hfill \\ \multicolumn{1}{c}{y_w}\\ \multicolumn{1}{c}{z_w}\end{array}\right]=R\left[\begin{array}{c}\hfill x_t\hfill \\ \multicolumn{1}{c}{y_t}\\ \multicolumn{1}{c}{z_t}\end{array}\right]+T$$

(17)

[x_w , y_w , z_w ] ^T is the coordinate in the world coordinate, [x_t , y_t , z_t ] ^T is the coordinate in the target coordinate, R is the orthogonal rotation transformation matrix, and T is the translation vector.

Taking the coordinate of P ₁, P ₂, P ₃ in world coordinate system and target coordinate system into formula (17), and combining the feature of orthogonal matrix, then we can get the value of R.

The rotation order is z, x, y. The corresponding rotation angles are γ, α, β. According to the Euler angle method, the parameters can be obtained from the following formula

$$R=\left[\begin{array}{ccc}\hfill \cos \beta \cos \gamma -\sin \alpha \sin \beta \sin \gamma \hfill & \hfill \cos \beta \sin \gamma +\sin \beta \sin \alpha \cos \gamma \hfill & \hfill -\sin \beta \cos \alpha \hfill \\ \multicolumn{1}{c}{-\cos \alpha \sin \gamma }& \hfill \cos \gamma \cos \alpha \hfill & \hfill \sin \alpha \hfill \\ \multicolumn{1}{c}{\sin \beta \cos \gamma +\cos \beta \sin \alpha \sin \gamma }& \hfill \sin \beta \sin \gamma -\cos \beta \sin \alpha \cos \gamma \hfill & \hfill \cos \alpha \cos \beta \hfill \end{array}\right]$$

(18)

From formula (18), we can get

$$\alpha =a\sin (r_{23})$$

(19)

$$\beta =-a\tan \left(\frac{r_{13}}{r_{33}}\right)$$

(20)

$$\gamma =-a\tan \left(\frac{r_{21}}{r_{22}}\right)$$

(21)

Experimental validations

The experiment platform of space control technology mainly contains a 6-DOF platform, a stereo camera, and a motion controller. The spacecraft model on the 6-DOF platform can realize all kinds of accurate poses (Figure 11).

Figure 11.

Experiment platform.

The relationship between the vision system and the 6-DOF platform is shown in Figure 12. The pose of the platform is adjusted to the initial status by the motion controller. At this time, all the pose parameters in the coordinate system of the 6-DOF platform (X_p-Y_p-Z_p ) are 0. By camera calibration, ¹⁶ we can obtain the relationship between the initial status of the 6-DOF platform and the vision system. The matrix transform is given in formula (22). R ₁ is a 3 × 3 orthogonal rotation transformation matrix, and T ₁ is the translation vector

$$\left[\begin{array}{c}\hfill x_w\hfill \\ \multicolumn{1}{c}{y_w}\\ \multicolumn{1}{c}{z_w}\end{array}\right]=R_1\left[\begin{array}{c}\hfill x_p\hfill \\ \multicolumn{1}{c}{y_p}\\ \multicolumn{1}{c}{z_p}\end{array}\right]+T_1$$

(22)

Figure 12.

Relationship between stereovision system and 6-DOF platform.

Experiment result

The 3D pose parameters include the position and attitude parameters. The x, y, and z are used to characterize the position information of the target spacecraft in the world coordinate. The yaw, pitch, and roll are used to characterize the attitude angles of the target aircraft relative to the world coordinate.

The experimental results show that the error of relative pose and the distance are linear relationship. The analysis of Figure 13(a)–(c) shows that in the distance from 2 to 4 m, the position errors are close to 0. As shown in Figure 13(d)–(f), in the distance from 1 to 3.5 m, the error of pitch angle is about 1°, the error of yaw angle is less than 1°, and the error of roll angle is less than 0.5°.

Figure 13.

Data analysis: (a) analysis of x data, (b) analysis of y data, (c) analysis of z data, (d) analysis of pitch data, (e) analysis of yaw data, and (f) analysis of roll data.

Comparison with existing method

The experiment results in Segal et al. ⁵ show that the accuracy of the relative position estimation can be up to 0.02 m, and the accuracy of the relative attitude estimation can be within 2°. The accuracy of the proposed algorithm is significantly higher than that in Segal et al. ⁵ Especially in the distance, there are more obvious advantages in the proposed algorithm. The attitude accuracy comparison in the same distance is shown in Table 1.

Table 1.

Comparison of algorithm accuracy.

Distance	2 m	3 m	4 m
Accuracy in this article	0.08°	0.55°	0.82°
Accuracy in Segal et al. 5	0.25°	0.75°	1°

Conclusion

This article proposes a stereovision-based method to estimate the relative pose between non-cooperative spacecrafts. The surface fitting method is applied to detect the edge of the object. The lines are fit by the points on the edge. The intersections of the lines are taken as the feature points of the object. The matching of the feature points between the two images is based on the epipolar geometry relationship. Since the parallel vision-based system is adopted, the external camera parameters need not be calibrated in the process of the 3D reconstruction. The 3D pose parameters of the spacecraft are estimated by the relative coordinate values in world coordinate and object coordinate. The accuracy of the method is validated on the experiment platform of space control technology and the results show that the higher accuracy is the advantage comparing with the existing algorithm.

Inheriting this article’s work, the future work will concentrate on the speed of the algorithm to achieve a real-time system.