An Extended Kalman Filter-Based Attitude Tracking Algorithm for Star Sensors

Efficiency and reliability are key issues when a star sensor operates in tracking mode. In the case of high attitude dynamics, the performance of existing attitude tracking algorithms degenerates rapidly. In this paper an extended Kalman filtering-based attitude tracking algorithm is presented. The star sensor is modeled as a nonlinear stochastic system with the state estimate providing the three degree-of-freedom attitude quaternion and angular velocity. The star positions in the star image are predicted and measured to estimate the optimal attitude. Furthermore, all the cataloged stars observed in the sensor field-of-view according the predicted image motion are accessed using a catalog partition table to speed up the tracking, called star mapping. Software simulation and night-sky experiment are performed to validate the efficiency and reliability of the proposed method.


Introduction
Star sensors are widely used for attitude determination in both orbiting and interplanetary spacecraft [1]. A star sensor typically operates in two modes: initial acquisition and tracking mode. The difference between them is whether the approximate attitude is known. In initial acquisition without a-priori attitude input, stars inside the whole field of view (FOV) are acquired and mapped to the matching stars in the sensor's star catalog by a search over the whole celestial sphere. The lost-in-space algorithms are well summarized by Benjamin [2]. In the case sufficient stars have been identified, the star sensor starts the calculation of attitude and goes into tracking mode. In contrast to initial acquisition mode, the sensor performs an attitude tracking and not a star tracking in tracking mode. The current attitude and angular rate is used to determine an expected attitude for the following measurement frame. The sensor uses this expected attitude and onboard star catalog to determine the star track windows to get the stars to be measured and evaluated. This method shortens the time to find the measurement stars and will also to eliminate spurious objects. Tracking mode is more efficient and reliable than initial acquisition mode and is the key operation mode in the lifetime of a star sensor. Attitude Tracking Algorithm directly affects the star sensor's performance.
At present, a few papers are published for tracking. Li [3] and Jiang [4] used the locations of recognized stars (reference stars) in previous star image as the centers of star track windows. If there is only one observed star in the window, the star is matched with the reference star corresponding to the track window. The radius of the window is in proportion to star sensor's angular rate. The larger the radius is, the more stars will be observed in one track window. In that case the angular distance between the observed stars and reference stars are used to find the correct match which is complicated and time consuming. In order to speed up the accessing the cataloged stars to identify the stars which enter or exit the sensor FOV, Jiang assigned the cataloged stars into sub-catalog. A partition table including all sub-catalogs is generated which is memory-consuming. Finally, according to stars transformation between inertial-based coordinate system and image plane, the star sensor's attitude is calculated by QUEST algorithm [5]. Samaan [1] presented two new recursive methods for identifying the stars within the FOV during tracking. The spherical polygon approach is used to sort the star position vector components in x, y and z-axes and then to access the cataloged stars within the FOV. However, it needs angular velocity obtained from the rate gyro and runs slower than jiang's method [4]. The star neighborhood approach is also used as a second method for star identification. In this method Samaan accesses all cataloged stars in the union of the neighborhoods of stars identified in the previous frame. Like Jiang [4], there comes the need for more memory to save the "star neighbor table". The attitude estimator adopted in both the proposed algorithms is the "second Estimator of the Optimal Quaternion" [6].
To improve the tracking performance under high dynamic condition, Sun [7] proposed a star tracking method which combines the spot-target-based optical flow analysis and Kalman filter to determine the angular velocity. Effective positions and size of the region of interest for detecting the star spot are provided. However, the algorithm is not suitable for use in orbit as it requires a large amount of calculation. Yu [8] proposed a multiexposure imaging based star tracking method for intensified star trackers which is not applicable to traditional star trackers without intensified image detector due to the low star sensitivity.
Typically, the key issues in tracking mode are determining the star locations within the FOV and an optimal attitude rapidly. A novel attitude tracking method is presented herein. Two techniques are introduced to improve the efficiency and reliability. Firstly, the star sensor is modeled as a nonlinear stochastic system with the state estimate providing the three degree-of-freedom attitude quaternion and angular velocity. The stochastic model uses the extended Kalman filter (EKF) [9] to estimate the attitude for the current and following frame. Secondly, inspired by [1,4], a partition table is established and accessed by spherical polygon approach involving a tradeoff between time versus memory consuming. The new tracking algorithm, presented here, is supported by both the software simulation and the hardware test results.
The remainder of this paper is organized as follows: Section 2 describes the algorithm, including the tracking model, attitude estimation, star catalog partition and initial angular velocity estimation criteria, in detail. Section 3 presents the simulation and experiment results to demonstrate the performance of the proposed algorithm. Finally, some conclusions are drawn.

Algorithm Description
In this section, we initially describe the tracking model of the star sensor. We then provide some implementation details, particularly on attitude estimation, star catalog partition and initial angular velocity estimation. Firstly, the EKF algorithm for attitude estimation is designed. Secondly, a novel star catalog partition method is presented to speed up accessing the cataloged stars while having a small memory requirement. Thirdly, the initial angular velocity estimation algorithm is described.

Tracking Models
After a full-sky star identification is finished to establish an initial attitude and angular velocity, the star sensor enters the tracking mode and recursively performs the following three steps as shown in Figure 1. 1 From EKF the attitude quaternion and its error covariance is extrapolated to the next frame in EKF. The expected attitude is used to access onboard partition table and star catalog by spherical polygon approach [1] to determine the centers of the star track windows. Meanwhile, the error covariance is used to determine the radius of the track windows. 2 Get the stars to be measured [10] and evaluated in the track windows of star image. If there is only one star in a reference star's neighborhood, the star in star image is matched with the star Sensors 2017, 17, 1921 3 of 18 in reference star image. Further, the match can be checked by the inter-star angles between the reference stars and the observed stars. 3 The coordinates of the matched stars in the observation star image are used as the measurements for EKF to estimate the optimal attitude quaternion of current frame and the angular velocity.
Sensors 2017, 17,1921 3 of 18 ③ The coordinates of the matched stars in the observation star image are used as the measurements for EKF to estimate the optimal attitude quaternion of current frame and the angular velocity.

EKF Based Attitude Estimation
The quaternion q q q q = q and the angular velocity selected as state variables: The quaternion q representing a rotation is given by [11] where u is the axis of rotation and θ is the angle of rotation about u . The quaternion satisfies a single constraint, which is A system model, that is, a function φ and additive noise σ required to calculate the next state in terms of the present state can be defined as [9]: σ is assumed to be additive Gaussian noise with zero mean and covariance matrix Q . The state transition function φ extrapolates from the state at time interval 1 k − to the next state at time interval k . The angular velocities are assumed constant so that The quaternion propagation in time must be derived. The relation between q and the angular velocity w is known to be [11]:

EKF Based Attitude Estimation
The quaternion q = q 0 q 1 q 2 q 3 and the angular velocity w = w x w y w z T are selected as state variables: s = q 0 q 1 q 2 q 3 w x w y w z T The quaternion q representing a rotation is given by [11] where u is the axis of rotation and θ is the angle of rotation about u. The quaternion satisfies a single constraint, which is q = 1.
A system model, that is, a function φ and additive noise σ required to calculate the next state in terms of the present state can be defined as [9]: σ is assumed to be additive Gaussian noise with zero mean and covariance matrix Q. The state transition function φ extrapolates from the state at time interval k − 1 to the next state at time interval k. The angular velocities are assumed constant so that w i (t k ) = w i (t k−1 ).
The quaternion propagation in time must be derived. The relation between q and the angular velocity w is known to be [11]: where the symbol '⊗' represents quaternion multiplication. Solving for q when w is constant gives: Since 2Ω Ω is an orthogonal and skew-symmetric matrix, its eigenvalues are ±i, ±i [12]. After solving for the closed form of the matrix exponential [13], the solution, when simplified, is: The state transitions for each variable have now been defined. The state transition function is: where: Q tran is the 4-by-4 quaternion transition matrix. The linearized state transition matrix Φ k is computed as the partial derivative of the state transition function with respect to each state variable and evaluated at time t = t k . The time dependency requires that it be computed at each state update. Formally, Φ k is a matrix of dimension 7-by-7: where Q tran is defined above in Equation (7). Q σ is defined as: Each column may be computed separately as given below. The partial derivative of Q tran q(t k ) with respect to each w i is: The partial derivatives of Ω with respect to w x , w y and w z are: The measurement update equation is [9]: The measurement noise v is zero mean, Gaussian sequence, with covariance matrix R. Measurement function h k comprises the 3-D rotation and perspective projection functions shown in Figure 2.
The partial derivatives of Ω with respect to x w , y w and z w are: The measurement update equation is [9]: q v (12) The measurement noise v is zero mean, Gaussian sequence, with covariance matrix R .
Measurement function ℎ comprises the 3-D rotation and perspective projection functions shown in Figure 2. The rotation matrix which describes the rotation from celestial reference frame to star sensor reference frame could be defined in terms of a unit quaternion [14]: q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  (14) The star locations ( ,y ) i i x is then calculated by using: (15) where f is the camera focal length.
The linearized measurement matrix is the partial derivative of the measurement function with respect to each state variable and evaluated at time t = tk. Formally, Hk is determined as: Hk is a 2n-by-7 matrix where n is the number of matched stars. For each matched star, the partial derivative of each coordinate needs to be computed: The rotation matrix which describes the rotation from celestial reference frame to star sensor reference frame could be defined in terms of a unit quaternion [14]: Let r = (r x , r y , r z ) be the 3-D coordinate of a unit vector in star sensor reference frame and s = (s x , s y , s z ) in celestial reference frame: The star locations (x i , y i ) is then calculated by using: where f is the camera focal length. The linearized measurement matrix is the partial derivative of the measurement function with respect to each state variable and evaluated at time t = t k . Formally, H k is determined as: H k is a 2n-by-7 matrix where n is the number of matched stars. For each matched star, the partial derivative of each coordinate needs to be computed: The partial derivatives of the observed stars x and y with respect to each element of r are given below: The partial derivatives of r with respect to each state variable can be calculated: The partial derivatives of R with respect to each rotational quaternion variable are:  Table 1. The (−) and (+) notations represent the estimate before and after the measurement update, respectively.
The partial derivatives of r with respect to each state variable can be calculated: The partial derivatives of R with respect to each rotational quaternion variable are:  Table 1. The (−) and (+) notations represent the estimate before and after the measurement update, respectively. k t r a n k 1

Star Catalog Partition
In order to speed up accessing the cataloged stars to identify the stars which enter or exit the sensor FOV, the celestial sky is divided nearly evenly spaced and completely covering of the sphere. Ideally, the vertices of a regular polyhedron represent evenly spaced points on a sphere. Unfortunately there are a limited number of regular polyhedrons, none of which provide enough points. However, the faces of the polyhedron can be divided into equal polygons with the vertices projected on the sphere to get enough points. Clearly, the minimal nonuniformity will result from the polyhedron with the most and consequently smallest faces. Also, the triangular or square faces would be the simplest to subdivide. The icosahedron, a regular 20-sided polygon whose faces are equilateral triangles, fits both requirements. The partitioning steps are as follows: 1 The celestial sphere can be divided into twenty equivalent zones by an inscribed icosahedron which is shown in Figure 4a. A pyramid is formed by connecting the center of the celestial sphere with the three vertices on one face of the icosahedron. The pyramids divide the celestial sphere totally into twenty equivalent zones. 2 Each of twenty regions is divided into several congruent equilateral triangles. As shown in the Figure 4b, each edge is divided into N + 1 segments, and each region is divided into (N + 1) 2 triangles. 3 When projected on the sphere, the vertices of these smaller triangles become the n = 10N 2 + 20N + 12 centers of sub-catalogs that completely cover the sphere, called sub-catalog centers. Each star in the catalog is then assigned to the sub-catalog of closest vertex. The sub-catalog centers number is computed as follows:

Star Catalog Partition
In order to speed up accessing the cataloged stars to identify the stars which enter or exit the sensor FOV, the celestial sky is divided nearly evenly spaced and completely covering of the sphere. Ideally, the vertices of a regular polyhedron represent evenly spaced points on a sphere. Unfortunately there are a limited number of regular polyhedrons, none of which provide enough points. However, the faces of the polyhedron can be divided into equal polygons with the vertices projected on the sphere to get enough points. Clearly, the minimal nonuniformity will result from the polyhedron with the most and consequently smallest faces. Also, the triangular or square faces would be the simplest to subdivide. The icosahedron, a regular 20-sided polygon whose faces are equilateral triangles, fits both requirements. The partitioning steps are as follows: ① The celestial sphere can be divided into twenty equivalent zones by an inscribed icosahedron which is shown in Figure 4a. A pyramid is formed by connecting the center of the celestial sphere with the three vertices on one face of the icosahedron. The pyramids divide the celestial sphere totally into twenty equivalent zones. ② Each of twenty regions is divided into several congruent equilateral triangles. As shown in the Figure 4b, each edge is divided into N + 1 segments, and each region is divided into ( ) 2 N 1 + triangles. ③ When projected on the sphere, the vertices of these smaller triangles become the 2 n 10N 20N 12 = + + centers of sub-catalogs that completely cover the sphere, called sub-catalog centers. Each star in the catalog is then assigned to the sub-catalog of closest vertex. The subcatalog centers number is computed as follows: The element of the partition table is structured as Table 2, which is arranged in descending order of the x component of the center vector.  table is structured as Table 2, which is arranged in descending order of the x component of the center vector. If the boresight of a star sensor is roughly known, the subcatalogs around the boresight are to be found quickly, so the guide stars within the FOV can be accessed immediately. The term θ is used to represent the maximum angle between the guide stars and the center of subcatalog. As shown in Figure 5, the necessary condition for stars to enter the FOV is that the angle between the center of the subcatalog and the boresight of the star sensor denoted by d is less than θ + FOV/2.
If the boresight of a star sensor is roughly known, the subcatalogs around the boresight are to be found quickly, so the guide stars within the FOV can be accessed immediately. The term θ is used to represent the maximum angle between the guide stars and the center of subcatalog. As shown in Figure 5, the necessary condition for stars to enter the FOV is that the angle between the center of the subcatalog and the boresight of the star sensor denoted by d is less than θ + FOV/2.
Similarly, let the angle between the boresight and the OYZ plane be The necessary condition for stars to enter the FOV can be expressed as Equation (23): Let the angle between the boresight and the OXY plane of celestial reference frame be α z , then the lower and upper bound of z component of the center vector are respectively: Similarly, let the angle between the boresight and the OYZ plane be α x , and the angle between the boresight and the OZX plane be α y , then the lower and upper bound of x and y component of the center vector are respectively: The necessary condition for stars to enter the FOV can be expressed as Equation (23): Hence, the search of subcatalogs around the boresight is accomplished using spherical polygon approach [1] and standard Binary Search technique. The above star catalog partition method greatly reduces searching time, while needs less memory space comparing the methods in ref. [1,4], as shown in Section 3.

Initial Angular Velocity Estimation
Initial angular velocity is derived from the quaternion difference. Taking account of Equation (5), the quaternion difference is written as: Taylor series expansion of ∆q 1 (t) is: The linear term is used to approximate w x : Similarly, we have: Typically, a star sensor slews within 8 • /s [15][16][17][18] and operates with an update rate of 10 Hz inducing some error less than 1 arcsec/s, so this approximation is reasonable.
Equation (25) can be expanded as:

Track Window Radius Determination
In the tracking mode, the stars are measured and evaluated in the star image tracking windows. If there is only one observed star in a window, the star is matched with the reference star corresponding to the track window. The larger the radius, the greater the probability that the window contains false stars due to hot spots, space radiation, stray light, etc. In order to improve tracking reliability, the tracking window radius should be as small as possible. The radius of the window is determined by the accuracy of the predicted attitude. It is possible to estimate the average star prediction accuracy, E boresight , by Equation (27): where E boresight is the predicted boresight accuracy obtained from the EKF error covariance extrapolation, N pixel is the number of pixels across the focal plane, and N star is the number of stars detected on the camera image. In order to improve the location accuracy of the star, the spot size of the star image is enlarged to a few pixels ranging from 3 to 5. Hence, the tracking window radius could be determined by Equation (28):

Simulation and Experiment
Several simulations and night sky tests were performed to evaluate the tracking algorithm presented in this paper. The simulations were implemented in MATLAB in the Microsoft Windows environment on a Core II 2.8 GHz PC.

Simulation
The star sensor configuration used for the simulations reported in this paper uses a 14.5 degrees circular FOV and 44.27 mm focal distance with a focal plane consisting of 2048 × 2048 pixels. A guide star catalog is made of 4338 stars below magnitude 6, based on the Smithsonian Astrophysical Observatory (SAO) J2000 star catalogue. From the discussion below parameter N in the star catalog partition method is 4, so the number of sub-catalogs divided is 252.
Simulation has been performed to evaluate the performance of the attitude tracking method under a variety of conditions. The results are obtained by varying initial conditions to characterize errors, convergence, and stability in the attitude quaternion estimate. Root mean square errors are presented.
Initial conditions requiring specification include the initial attitude and the error covariance matrix which are obtained from the iterative identification algorithm [19] and the QUEST algorithm [5]. Process noise given by the covariance matrix Q and measurement noise given by the covariance matrix R is also specified initially, remaining constant throughout. Table 3 gives the initial states for the simulation, as well as the assumed initial error covariance diagonal terms and the process noise covariance diagonal terms. (1 × 10 −9 , 1 × 10 −9 , 1 × 10 −9 , 1 × 10 −9 ) (1 × 10 −6 , 1 × 10 −6 , 1 × 10 −6 ) (rad/s) 2 Actual initial attitude quaternion is (0.54435, 0.15051, −0.07010, 0.82226), and true angular rate is (−0.03, 0.04, −0.02) rad/s. The angular velocity is assumed to be constant in the tracking models. Taking into account the actual conditions, the fluctuation in angular velocity is assumed to be 0.001 rad per second. Hence, the process noise covariance diagonal terms are {1 × 10 −9 , 1 × 10 −9 , 1 × 10 −9 , 1 × 10 −9 , 1 × 10 −6 , 1 × 10 −6 , 1 × 10 −6 }, while the off-diagonal terms are zero. The star measurement accuracy is determined by the star signal-to-noise ratio (SNR). The higher the SNR is, the more accuracy a star sensor can achieve. The measurement noise is a Gaussian with a standard deviation ranging from 0.04 to 0.18 pixel [18]. Since the noise for the points is independent with respect to x and y as well as each other, the measurement error covariance matrix is a straightforward diagonal matrix whose diagonal elements are the noise variances while the off-diagonal terms are zero. The simulated sample interval is 0.1 s. Figures 6 and 7 show the results of a sample run where the true states are shown along with the estimated values. The simulation was run for a time of 250 s with a measurement interval of 0.1 s. The tracking is always successful. Figure 8 gives the attitude estimation errors over the same period of time from the same sample data represented in Euler angles around the x-, yand z-axis, whose standard deviations are 0.5, 0.5, and 5.7 arc seconds respectively. Figure 9 gives the estimation errors for angular velocity about the x-, yand z-axis whose standard deviations are 2.6 × 10 −5 , 2.3 × 10 −5 , 1.3 × 10 −4 rad/s respectively. The accuracy along the boresight of the star sensor (z-axis) is less than the accuracy cross the boresight (x-and y-axis), which is an inherent characteristic of a star sensor. The ratio of the uncertainty along the boresight to the uncertainty cross the boresight spans from 7 to 90 [20]. As can be seen, the estimates track the state variables reasonably well. Figure 10 shows the maximum prediction error of the star locations in each tracking cycle. It is noted that the initial angular velocity estimations are inaccurate so that the prediction errors in the first frame are larger which are 2.27, 1.45 pixels respectively about x-and y-coordinate. The prediction errors in the next frame are less than 0.2 pixel.
boresight to the uncertainty cross the boresight spans from 7 to 90 [20]. As can be seen, the estimates track the state variables reasonably well. Figure 10 shows the maximum prediction error of the star locations in each tracking cycle. It is noted that the initial angular velocity estimations are inaccurate so that the prediction errors in the first frame are larger which are 2.27, 1.45 pixels respectively about x-and y-coordinate. The prediction errors in the next frame are less than 0.2 pixel.  track the state variables reasonably well. Figure 10 shows the maximum prediction error of the star locations in each tracking cycle. It is noted that the initial angular velocity estimations are inaccurate so that the prediction errors in the first frame are larger which are 2.27, 1.45 pixels respectively about x-and y-coordinate. The prediction errors in the next frame are less than 0.2 pixel.        The smaller the error is, the smaller the size of track windows is, and the more robust the tracking is. As the star location prediction error is relatively small, as shown in Figure 10, the radius of the track windows is set to 16 pixels ensuring a unique match and not updated according to the EKF error covariance. The above results are for a single representative run and are not necessarily typical. A further measure of performance is the error over a large number of trials. Using the same initial conditions given above, a series of tests were performed to experimentally measure the error. Figures 11  and 12 show the calculated root mean square (RMS) error over 100 sample runs. The square root of the corresponding mean diagonal element from the calculated covariance error matrix is also shown in each figure. In the nonlinear EKF, the covariance matrix depends on the measurements and does not necessarily represent the actual error covariance [21]. In these tests, the covariance matrix values correspond closely to the actual error values for the attitude estimation while minor differences are present for the angular velocity estimation.
conditions given above, a series of tests were performed to experimentally measure the error. Figures 11  and 12 show the calculated root mean square (RMS) error over 100 sample runs. The square root of the corresponding mean diagonal element from the calculated covariance error matrix is also shown in each figure. In the nonlinear EKF, the covariance matrix depends on the measurements and does not necessarily represent the actual error covariance [21]. In these tests, the covariance matrix values correspond closely to the actual error values for the attitude estimation while minor differences are present for the angular velocity estimation.  and 12 show the calculated root mean square (RMS) error over 100 sample runs. The square root of the corresponding mean diagonal element from the calculated covariance error matrix is also shown in each figure. In the nonlinear EKF, the covariance matrix depends on the measurements and does not necessarily represent the actual error covariance [21]. In these tests, the covariance matrix values correspond closely to the actual error values for the attitude estimation while minor differences are present for the angular velocity estimation.

Night Sky Experiments
The night sky tests have been performed at the Xinglong Astronomical Observatory of the Chinese Academy of Sciences. The star sensor was placed on a rotating table with vertical axis as shown in Figure 13. The parameters of the star sensor are listed in Table 4.  Limited by the sensitivity of the image sensor, the angular rate is set to 2 degrees per second around y-axis of the star sensor. The estimated quaternion is plotted in Figure 14. In order to assess the 3-axis attitude measurement accuracy we create a reference quaternion (qr) which represents as good as possible the trajectory track under investigation. For that purpose the four elements of the measured quaternions were fitted over their time stamp to a 5th degree polynomial. The difference between the reference quaternion which is smooth due to 5th order Limited by the sensitivity of the image sensor, the angular rate is set to 2 degrees per second around y-axis of the star sensor. The estimated quaternion is plotted in Figure 14.  Limited by the sensitivity of the image sensor, the angular rate is set to 2 degrees per second around y-axis of the star sensor. The estimated quaternion is plotted in Figure 14. In order to assess the 3-axis attitude measurement accuracy we create a reference quaternion (qr) which represents as good as possible the trajectory track under investigation. For that purpose the four elements of the measured quaternions were fitted over their time stamp to a 5th degree polynomial. The difference between the reference quaternion which is smooth due to 5th order In order to assess the 3-axis attitude measurement accuracy we create a reference quaternion (qr) which represents as good as possible the trajectory track under investigation. For that purpose the four elements of the measured quaternions were fitted over their time stamp to a 5th degree polynomial. The difference between the reference quaternion which is smooth due to 5th order polynomial fitting and the measured quaternion both at the same time stamp represents the attitude error which can be expressed as Euler angle rotation around x-, y-, z-axis as shown in Figure 15. The standard deviations are 5.3, 9.8, and 50.8 arc seconds respectively, which are larger than simulation errors due to the larger star position error caused by low SNR value, the less star number detected in the field of view caused by low star sensitivity and the more star vector uncertainty caused by atmosphere turbulence. It can be seen that the y-axis error is larger than the x-axis error due to introduction of system error of rotating table. It is shown in Figure 16 that the prediction error of the star locations is less than 3 pixels.
The tracking time in every cycle is shown in Figure 17. The tracking time in every cycle keeps less than 40 milliseconds, so this algorithm is suitable for real-time tracking. the field of view caused by low star sensitivity and the more star vector uncertainty caused by atmosphere turbulence. It can be seen that the y-axis error is larger than the x-axis error due to introduction of system error of rotating table. It is shown in Figure 16 that the prediction error of the star locations is less than 3 pixels.
The tracking time in every cycle is shown in Figure 17. The tracking time in every cycle keeps less than 40 milliseconds, so this algorithm is suitable for real-time tracking.  atmosphere turbulence. It can be seen that the y-axis error is larger than the x-axis error due to introduction of system error of rotating table. It is shown in Figure 16 that the prediction error of the star locations is less than 3 pixels.
The tracking time in every cycle is shown in Figure 17. The tracking time in every cycle keeps less than 40 milliseconds, so this algorithm is suitable for real-time tracking.

Comparison and Analysis
The time to predict the frame of the stars at each time step constituted with the time to search the partition table and the time to access the catalog is related with the number of subcatalogs. The larger the subcatalog number is, the longer the time searching the partition table is. However, the average number of the stars in one subcatloag decreases as the subcatalog number increases, and it will take shorter time to access the catalog. To determine the optimum number of subcatalogs, the distribution of the stars on the celestial sphere was assumed to be uniform, and the time to simulate

Comparison and Analysis
The time to predict the frame of the stars at each time step constituted with the time to search the partition table and the time to access the catalog is related with the number of subcatalogs. The larger the subcatalog number is, the longer the time searching the partition table is. However, the average number of the stars in one subcatloag decreases as the subcatalog number increases, and it will take shorter time to access the catalog. To determine the optimum number of subcatalogs, the distribution of the stars on the celestial sphere was assumed to be uniform, and the time to simulate the reference star image was computed against N using a Monte Carlo simulation method. The time is normalization to that accessing all stars in the catalog. The result is shown in Figure 18. For an FOV of Φ20 • , the time consumed is 0.08 when N is 11. According to the results, to achieve the shortest consumed time, N is decided to be in the range of 3 to 5.
To compare the speed of Samaan [1] and Jiang [4], the time to simulate the reference star image was computed against FOV by Monte Carlo simulation. The result is show in Figure 19. The time taken by our method is shorter than ref. [1], but longer than ref. [4]. As FOV increases, the gap between ours and [4] is reduced, and that between ours and [1] is enlarged. For an FOV of Φ20 • , the time consumed is 0.05, decreased by 70 percent compared to [1] and increased by 50 percent more than [4].

Comparison and Analysis
The time to predict the frame of the stars at each time step constituted with the time to search the partition table and the time to access the catalog is related with the number of subcatalogs. The larger the subcatalog number is, the longer the time searching the partition table is. However, the average number of the stars in one subcatloag decreases as the subcatalog number increases, and it will take shorter time to access the catalog. To determine the optimum number of subcatalogs, the distribution of the stars on the celestial sphere was assumed to be uniform, and the time to simulate the reference star image was computed against N using a Monte Carlo simulation method. The time is normalization to that accessing all stars in the catalog. The result is shown in Figure 18. For an FOV of Φ20°, the time consumed is 0.08 when N is 11. According to the results, to achieve the shortest consumed time, N is decided to be in the range of 3 to 5.
To compare the speed of Samaan [1] and Jiang [4], the time to simulate the reference star image was computed against FOV by Monte Carlo simulation. The result is show in Figure 19. The time taken by our method is shorter than ref. [1], but longer than ref. [4]. As FOV increases, the gap between ours and [4] is reduced, and that between ours and [1] is enlarged. For an FOV of Φ20°, the time consumed is 0.05, decreased by 70 percent compared to [1] and increased by 50 percent more than [4]. A comparison in consumed memory among these methods is listed below. For an FOV of Φ20°and catalog star number of 5000, our methods needs 5 k bytes, which is a decrease by 83 percent compared to [4] and by 91 percent compared to [1].

Conclusions
To improve the performance of star sensor in case of high attitude dynamics, a novel attitude A comparison in consumed memory among these methods is listed below. For an FOV of Φ20 • and catalog star number of 5000, our methods needs 5 k bytes, which is a decrease by 83 percent compared to [4] and by 91 percent compared to [1].

Conclusions
To improve the performance of star sensor in case of high attitude dynamics, a novel attitude tracking algorithm based on extended Kalman filter is presented. The approach model the star sensor as a nonlinear stochastic system with the state estimate providing the three degree-of-freedom attitude quaternion and angular velocity based on the extended Kalman filter (EKF). From EKF the attitude quaternion and its error covariance are extrapolated to the next frame used to predict the centers and radius of the star track windows. In order to speed up the prediction a star catalogue partition method by an inscribed regular icosahedron is introduced for a fast search of the star catalog while having a small memory requirement. Then the matched stars in the track windows of star image are measured and used as the measurements for EKF to estimate the optimal attitude quaternion of current frame and the angular velocity.
The simulation and night sky experiment results demonstrate that this algorithm provides effective and reliable attitude tracking. In the simulation the EKF always converge, and the covariance matrix values accord with the actual error values for the attitude estimation while minor differences are present for the angular velocity estimation. In addition, it is noted that the accuracy along the boresight of the star sensor (z-axis) is less than the accuracy cross the boresight (x-and y-axis), which is an inherent characteristic of a star sensor. Moreover, it can be seen that the prediction error of the star locations is less than 3 pixels and the tracking time in every cycle keeps less than 40 milliseconds in the night sky experiment. Hence, the algorithm is particularly suitable for use in orbit as it is computationally efficient and has a small memory requirement. Future work can involve on-orbit validation of the attitude tracking algorithm.