A Star Sensor On-Orbit Calibration Method Based on Singular Value Decomposition

The navigation accuracy of a star sensor depends on the estimation accuracy of its optical parameters, and so, the parameters should be updated in real time to obtain the best performance. Current on-orbit calibration methods for star sensors mainly rely on the angular distance between stars, and few studies have been devoted to seeking new calibration references. In this paper, an on-orbit calibration method using singular values as the calibration reference is introduced and studied. Firstly, the camera model of the star sensor is presented. Then, on the basis of the invariance of the singular values under coordinate transformation, an on-orbit calibration method based on the singular-value decomposition (SVD) method is proposed. By means of observability analysis, an optimal model of the star combinations for calibration is explored. According to the physical interpretation of the singular-value decomposition of the star vector matrix, the singular-value selection for calibration is discussed. Finally, to demonstrate the performance of the SVD method, simulation calibrations are conducted by both the SVD method and the conventional angular distance-based method. The results show that the accuracy and convergence speed of both methods are similar; however, the computational cost of the SVD method is heavily reduced. Furthermore, a field experiment is conducted to verify the feasibility of the SVD method. Therefore, the SVD method performs well in the calibration of star sensors, and in particular, it is suitable for star sensors with limited computing resources.


Introduction
A star sensor is a navigation system that obtains the attitude information of the carrier by the observation of stars. They are the most accurate optical attitude sensors, at present [1]. Due to their high navigation accuracy, strong autonomy, and lack of cumulative error, they are favored in the aerospace industry [2,3]. As the "eyes" of the spacecraft, the accuracy of the star sensor determines the performance of the spacecraft directly. The star sensor is an optical device whose accuracy depends on the imaging quality and the accuracy of the optical parameters, including focal length, principal point, and distortion [4,5]. Therefore, calibration is one of the key technologies for a star sensor.
The calibration methods of star sensors can be divided into two categories: ground-based calibration and on-orbit calibration. Generally, multi-frame observations should be accumulated for the ground-based calibration method; hence, this method often relies on a fixed platform and complex experimental installations [6,7], and the cost of ground-based calibration is high. Additionally, as the working environment of a star sensor is different from the calibration environment, the parameters may change in orbit [8,9]. In most cases, on-orbit calibration is carried out with observed data during

Camera Model
It is generally accepted that a conventional camera model includes three parts: extrinsic parameters, intrinsic parameters, and distortion parameters [22,23]. Notably, a star sensor is mainly used to calculate the attitude or the rigid rotation between the camera coordinate and the inertial coordinate, so there is no need to solve for the extrinsic parameters in the calibration. Neglecting the extrinsic parameters, the coupling effect between intrinsic and extrinsic parameters can be avoided to improve the calibration accuracy. Therefore, in this paper, the camera model of a star sensor only consists of intrinsic and distortion parameters.
The reference frames of the camera model involve three coordinate frames: the camera frame, the physical image frame, and the image frame, as shown in Figure 1. Figure 1a shows the relationship between the camera frame and the physical image frame, and Figure 1b shows the relationship between the physical image frame and the image frame. These reference frames are defined based on the pinhole camera model, where O-XYZ is the camera reference frame, the origin O of the camera frame is the projection center, and the Z-axis is the optical axis. The physical image frame o'-th u'v' is related to the detector plane. This frame is parallel to the XOY plane of the camera frame; the intersection of this plane and the optical axis is the principal point o'; and the u'and v'-axes are parallel to the Xand Y-axes of the camera frame, respectively. The distance between these two planes along the optical axis is the focal length f . The uand v-axes of the image reference frame o-uv (in Figure 1b) are parallel to the u'and v'-axes of the physical image frame, respectively, and the origin of the image reference frame is located in the top left of the image; so, the coordinates of the principal point in the image frame are Let w = [X, Y, Z] T be an arbitrary unit star vector [16] with respect to the camera reference frame, and its projection on the image frame is p = [u, v] T . The perspective projection relationship between w and p can be represented as [22,24]: where [u, v, 1] T are the homogeneous coordinates of the point p, and f u and f v are: f where f is the focal length and D u and D v are the horizontal and vertical pixel sizes, respectively. The relationship between the factors D u and D v is: where s is the aspect ratio (the height of a pixel compared to the width of that pixel) [22]. Normally, the reference values of D u and D v are given in the datasheet of the camera. From Equations (2)-(4), f u and f v are linearly dependent on the parameters s and f , and so, u 0 , v 0 , s, and f are the intrinsic parameters to be calibrated. As a result of the imperfection of the camera lenses, distortion should be considered in the camera model. Lens distortion is usually expressed as: where u and v are the distortion-free coordinates in Equation (5), u d and v d are the corresponding coordinates with distortion, and δ u (u, v) and δ v (u, v) are the distortions in the u and v directions, respectively. It has been shown that higher order distortions may cause numerical instability [25], so here, we consider only the first-and second-order radial distortion: where k 1 and k 2 are the distortion coefficients of the radial distortion, and r 2 is defined as [24]: With the model above, given the position of the observed star centroid in the image frame p d = [u d , v d ] T , we can calculate the corresponding unit star vector in the camera frame w = [X, Y, Z] T . We express the relationship as: where F(·) is the back-projection function with distortion.

Calibration Reference: Singular Values
This calibration method is based on the invariant characteristic of the singular-value decomposition of the star sensor. The principle is presented as follows [19].

Singular-Value Decomposition of the Star Sensor
Define w i as an observed star vector in the camera frame and v i as a guide star vector in the inertial frame. The transformation between these two frames is: where W and V are the column vector matrices: and C is the direct cosine matrix (DCM), which indicates the transformation from the inertial frame to the camera frame. Hence, C is an orthogonal matrix. The SVD is a general decomposition method for matrices. With the SVD, the matrices W and V can be decomposed into: where P w and P v are 3 × 3 orthogonal matrices of left singular vectors p wi and p vi (i = 1,2,3), Q w and Q v are N × N orthogonal matrices of right singular vectors q wi and q vi (i=1,2,. . . , N), and Σ w and Σ v are 3 × N diagonal matrices where the diagonal elements are the singular values σ wi and σ vi (i=1,2,3) of W and V, respectively. For no less than three stars in the field of view (FOV), there are three non-zero singular values, and the SVD is unique. The singular values have the following property.

Invariant Singular Values
Post-multiplying Equation (9) by W T , we obtain: Substituting Equations (12) and (13) into Equation (14) yields: or: where S w and S v are diagonal matrices with eigenvalues σ 2 wi and σ 2 vi (i=1,2,3) of WW T and V V T , respectively. As C is an orthogonal matrix, Equation (14) is a similarity transformation. Thus, the eigenvalues of the WW T and V V T are equal; namely, Meanwhile, the positive singular values of W and V are equal: Therefore, we can conclude that the singular values of the star vectors remain constant under coordinate transformation.

On-Orbit Calibration Method Based on Singular-Value Decomposition
The singular values can be calculated using the data obtained by the star sensor itself, and the invariance of the singular values provides a new way to calibrate the optical parameters of star sensors in orbit.
Supposing that SV(·) indicates the singular value-solving operator, Equation (18) can be presented as: Substituting Equation (8) into Equation (19) yields: where P d = [P d1 , P d2 . . . P dN ] is the collection of observed star coordinates in the image frame. After star identification, the observed star coordinates P d and the corresponding star vectors V in the star catalog are matched with each other. Thus, according to Equation (20), the singular values can be obtained by the star vectors V, and they, also, can be calculated by the camera parameters and the observed star coordinates P d . The accuracy of the star catalog is very high, and the singular values calculated by V have good precision. Obviously, Equation (20) is suitable for the measurement equation.
The measurement equation is non-linear, and thus, we estimate the camera parameters based on the extended Kalman filter (EKF) [26]. The state transition and measurement models are: where , and x k and x k−1 are the state parameters of the star images marked k and k-1, respectively. Suppose that the camera parameters are constant, and let I 6×6 be an identity matrix. Then, z k = [σ 1 , σ 2 , σ 3 ] T can be calculated by the guide star vectors V; h(x k ) is a simplified representation of SV(F(s, f , u 0 , v 0 , k 1 , k 2 , P d )); and n c represents the measurement error caused by the noise. Starting with initial estimates of the noise covariance P 0 and the parameters x 0 , which may be obtained from the ground calibration, we process the calibration frame-by-frame. For the k th star image, the EKF prediction equations are: where Q is the covariance matrix of the prior estimation error. In this model, the parameters are constant, and so, the theoretical Q is a null matrix. However, if this is so, the estimator is just a sequential least squaresestimator, and the artificial process noise is a common technique for forgetting old measurements. The EKF update equations are: where R is the covariance matrix of measurement noise and H k is the Jacobian matrix of the measurement The model of the h is complex, so we adopt numerical differentiation to calculate the Jacobian matrix.

Further Study on the SVD Calibration Model
Concerning Equation (28), if all the stars in the FOV are used to calculate the singular values, it is easy to find that there are only three rows in the Jacobian matrix H k . The order of H k is less than the number of the state parameters, so the observability of the estimator is bad. In order to solve this problem, the stars should be put into several groups, such that the order of H k can satisfy the estimator. There are many different combinations of the stars, however; how to choose the proper combination is discussed below.

Observability Analysis
Observability is an indicator to evaluate the feasibility of the system; namely, with different models, the same input derivation may cause different output derivations. If the magnitude of the output derivation is bigger, then the observability is better and the system is more feasible, and vice versa. In this section, we define the observability as the infimum of the output derivation [27].
According to the definition of the observability, we can obtain: where δx is the input derivation vector with the same derivation for all the parameters, δz k is the output derivation, and H k is the Jacobian matrix. As the precision of different parameters varies considerably in practice, this means δx cannot represent the actual derivation of different parameters. Thus, δx should be weighted according to the magnitudes of different precision, where W is the diagonal weighting matrix whose elements are the typical precision of the parameters (see Table 1), according to the star sensor in the field experiment. Hence, the observability matrix is: Here, the SVD is used again. With the SVD of the observability matrix, Equation (29) can be represented as: where P k and Q k are orthogonal matrices of left singular vectors and right singular vectors, respectively. To make sure that H k is observable, Σ k is defined as a 6 × N diagonal matrix where the diagonal elements are the non-zero singular values σ i (i=1∼6). As P k and Q k are orthogonal matrices, we can calculate: where δx 2 and δz k 2 are the two-norms of δx and δz k , respectively. The infimum of the output derivation is: where σ min is the minimum singular value (MSV) of the observability matrix. It is obvious that, with larger σ min , the infimum of the output derivation is bigger and the observability is better. Therefore, we adopt σ min as the indicator to evaluate the performance of the system and find a suitable calibration model.

Star Combination Models
With regards to no less than three stars, there are three non-zero singular values, and so, the combination can be formed by three stars, four stars, and so on. If all the combinations are used to constitute the measurement models, the computational cost is huge. In addition, the combinations are not independent, and the information for calibration does not increase linearly with the number of combinations. Here, we discuss several models based on the different number of stars in the combination.
The observability analysis of these models is shown in Figure 3, where the MSVs of the Jacobian matrix of each model are calculated frame-by-frame. The MSVs of Model-1, Model-2, and Model-3 interweave with each other, and thus, we suspect that, for different star images, combinations of different star numbers have different observabilities. The MSVs of Model-4 are bigger, in most cases, which means that the observability of Model-4 is better than the others (this may be because Model-4 combines the advantages of various star-number combinations), and so, we take Model-4 as the optimal combination model. This is only a rough discussion of the combination models; further studies on this problem are still underway.

Singular Value Selection
During the research, one unanticipated finding was that the sensitivities of the three singular values were quite different, as shown in Figure 4: σ 1 , σ 2 , and σ 3 are three singular values with descending order, and the observability of σ 1 is worse than σ 2 and σ 3 . Thus, we should consider whether all three singular values are suitable for the estimation of the camera parameters, and we discuss this question with respect to the physical interpretation of the SVD of the star vector matrix W (see Figure 5).  The black vector depicts the star vector w i . According to the definition of the SVD, the left singular vectors p wi (i = 1, 2, 3) are unit vectors, so the square of the maximum singular value of W is: where p max represents the left singular vector associated with the maximum singular value σ 1 . Equation (34) illustrates that p max is the vector that maximizes the sum of projections of each star vector w i onto p max . In the same way, p min , associated with the minimum singular value σ 3 , is the vector that minimizes the sum of projections of each star vector w i onto p min . Then, p middle , which is associated with the intermediate singular value, is perpendicular to the plane generated by p max and p min . An illustration of p max , p middle , and p min is shown in Figure 5.
As p max and w i are unit vectors, Equation (34) can be rewritten as: where α i is the angle between w i and p max , which represents the projection of the star vector w i onto p max . Therefore, the derivation of σ 1 is the magnitude of sin α i . As the plane generated by p middle and p min is perpendicular to p max , σ 2 and σ 3 are proportional to sin α i , and the derivations of σ 2 and σ 3 are on the magnitude of cos α i . With regards to a small FOV, sin α i is smaller than cos α i , and so, the derivation of σ 1 is smaller than σ 2 and σ 3 , which is in line with the observability analysis. As a result, in order to reduce the calculation, we can employ σ 2 and σ 3 for the measurements. As shown in Figure 5, the status of σ 2 and σ 3 is similar. The values of σ 2 and σ 3 are the sum of the projections of sin α i onto the p middle and p min axes, respectively. Depending on the position of the two axes, the observabilities of σ 2 or σ 3 may be large or small, resulting in an unstable effect when using σ 2 or σ 3 alone.

Results and Discussion
The simulation experiment and field experiments were conducted to verify and analyze the SVD on-orbit calibration method.

Simulation Experiments
Star data were simulated, with a 19.14 • × 11.18 • FOV, for a 1920 × 1080 pixel-array star sensor at a 2-Hz update rate. The other simulation parameters are listed in Table 1. The simulation data consisted of three datasets: a set of 3D star vectors in the inertial frame S 3D , a set of corresponding 2D star coordinates in the image frame S 2D , and the set of 2D star coordinates with normally-distributed noise (with a standard deviation of 0.5 pixel), which is referred to as S 2D . For each experiment, S 3D , S 2D , and S 2D of 2500 images were simulated. These images were divided into two groups, where the first 2400 images were used for calibration (according to the experiment, we found that 2400 images could ensure the convergence of calibration), and the last 100 images were used to evaluate the performance of the calibration. The calibration programs were written in MATLAB and run on an Intel Core i5-8300H CPU.
In order to evaluate the calibration method fully, the evaluation criteria are presented as follows. As the residual of a singular value cannot directly reflect the influence of calibration on the star sensor, the residual of angular distance between stars, which is normally used to evaluate the performance of on-orbit calibration, is adopted in this paper. Thus, two criteria to evaluate the calibration methods were employed: Criterion A: This criterion is based on the simulation data without errors. With the estimation of camera parameters, S 2D can be back-projected to the camera frame, and the angular distances between them are calculated. The residual errors are obtained by these angular distances and the corresponding angular distances calculated by S 3D . Then, the root-mean-squared error (RMSE) of the angular distances (A|δ RMSE ) in each star image is calculated. The mean value (A|µ δ RMSE ) and the standard deviation (A|σ δ RMSE ) of A|δ RMSE in the last 100 images are used as the evaluation indices of Criterion A.
Criterion B: The difference between this criterion and Criterion A is that this model uses the noisy 2D simulation data S 2D to calculate the RMSE of the angular distances (B|δ RMSE ) in each star image. The evaluation indexes of Criterion B are presented as the mean value (B|µ δ RMSE ) and the standard deviation (B|σ δ RMSE ) of B|δ RMSE in the last 100 images. Therefore, A|µ δ RMSE and A|σ δ RMSE are only related to the estimation errors of camera parameters, and these evaluation indexes can only be used in the simulation. Furthermore, B|µ δ RMSE and B|σ δ RMSE are related to both the estimation errors of camera parameters and the centroid noise. Test 1. Performance of different models: The aim of this simulation test was to evaluate the performance of the star combination models presented in Section 3.3.2. A star catalog was formed by taking stars, up to a visual magnitude of 5.5, from the Tycho-2 catalog. The initial values were set as s= 1, f = 15.5 mm, u 0 = 960, v 0 = 540, k 1 = 0, and k 2 = 0. The calibration results are listed in Table 2. The derivation of the parameters shown in Table 2 is expressed as fractional changes, due to the different scales of the parameters. The results showed that the parameters s, u 0 , v 0 , and k 2 obtained by Model-2 and Model-3 had better accuracy than Model-1; f obtained by Model-1 and Model-2 had better accuracy than Model-3; and k 1 obtained by Model-1 had better accuracy than Model-2 and Model-3. On the basis of these results, we conclude that, for different parameters, each of the first three models had advantages and disadvantages. However, the precision of the parameters calibrated by Model-4 was always at an intermediate level, which means that Model-4 combined the properties of the other combination models. With regard to the error evaluation Criterion A, the performance of Model-4 was better than the others. This is consistent with the analysis in Section 3.3.2, so we used Model-4 for the other experiments.

Test 2. Performance of different singular values:
The simulation conditions were the same as in Test 1. In this test, all combinations of the singular values were used to construct the Jacobian matrix, and the calibration results are as follows.
According to the results in Table 3, the three best combinations were {σ 1 , σ 2 , σ 3 }, {σ 2 , σ 3 }, and {σ 2 }. This result indicates a good agreement with the analysis in Section 3.3.3; namely, the observabilities of σ 2 and σ 3 were better than that of σ 1 . Particularly, in this test, the observability of the Jacobian matrix related to σ 2 was better than σ 3 , so the calibration results achieved only by σ 2 were good enough. Through several experiments, however, we found that the performance of {σ 2 } was unstable, and in this experiment, the relatively high standard deviation A|σ δ RMSE gave some indication of the stability problem. This may be due to some information in σ 1 and σ 3 that could help with calibration, and so, the combinations {σ 1 , σ 2 , σ 3 } and {σ 2 , σ 3 } are recommended. The calibration precisions of {σ 1 , σ 2 , σ 3 } and {σ 2 , σ 3 } were similar, yet the average calibration time per frame of {σ 2 , σ 3 } was shorter than that of {σ 1 , σ 2 , σ 3 }, which means that the computation of {σ 2 , σ 3 } was smaller than {σ 1 , σ 2 , σ 3 }. Thus, the combination {σ 2 , σ 3 } was adopted for the following tests. Test 3. Comparison with the calibration method based on the angular distance: The calibration method based on the angular distance (AD) method is widely used for the on-orbit calibration of star sensors, so we further evaluated the performance of the proposed method by comparing it with the AD method.
We set up three experiment groups, where the difference of the groups was the limiting visual magnitude of the star sensor, which led to different numbers of stars in the FOV. The average numbers of stars in the groups were 31.4, 19.1, and 7.7 with the limiting visual magnitudes of 6, 5.5, and 4.6, respectively. The results are listed in Table 4. It is interesting that, under Criterion A, the performance of the SVD method was close to, or even better than, the AD method, and the SVD method was more stable. However, under Criterion B, the performance of the SVD method was always worse than the AD method. This phenomenon occurred because the AD method uses the angular distance as the object of optimization estimation, so the comprehensive calibration results, considering both the camera parameters and the centroid noise, were better. It must also be mentioned that the estimation of the attitude of the star sensor was based on the star vector itself, rather than the angular distance, such that the accuracy of the camera parameters was more important; namely, Criterion A could better reflect the calibration effect than Criterion B. Therefore, the accuracy of the SVD method was as good as the AD method. The convergence rates of the two methods were also similar; as shown in Figure 6, the RMSE was calculated between the nominal angular distances and the angular distances that were obtained by the back-projection with calibrated parameters, which was the result of Group 2. Nevertheless, the elapsed time of the SVD method was significantly shorter than the AD method, especially for the case with a large number of stars in the FOV. Concerning Group 1, the calculation speed of the SVD method was improved by 94.48%, compared to the AD method.

Field Experiment
The SVD method is a novel on-orbit calibration method for star sensors, and therefore, it is necessary to verify the feasibility of this method by a field experiment. We tested the proposed on-orbit calibration method on the ground; the experiment platform of the field experiment is shown in Figure 7. A motorized zoom lens was adopted; the focal length was adjusted to about 16 mm; the FOV was about 19.14 • × 11.18 • . The resolution of the CCD is 1920 × 1080 and the pixel size was 2.9 µm × 2.9 µm. The experiment was carried out at the update rate of 2 Hz. To avoid the influence of the atmospheric refraction, the star sensor was placed perpendicular to the ground. The average number of observed stars in the FOV was 14.7. As the nominal centroid coordinates of observed stars were unknown, we could only use Criterion B to evaluate the performance of the field calibration experiment. The initial values of the parameters were the same as in the simulation. The estimation results of the parameters are shown in Figure 8, where the arrows in the images represent the points of convergence, and the order of the convergence of the parameters was f , k 2 , k 1 , u 0 , u 0 , and then s. With about 2200 images, the whole estimation was convergent; B|µ δ RMSE was 8.33", and B|σ δ RMSE was 0.803". We found that the performance in the field experiment was better than the simulation. This was because Criterion B employed the noisy simulation data, and this noise might affect the results; namely, it means the centroid errors in the field experiment were smaller than the assumption in the simulation. Therefore, the SVD method can perform well in the on-orbit calibration of a star sensor.

Conclusions
In this study, we proposed an on-orbit calibration method for star sensors by considering the fact that the singular values of the star vectors remain constant under coordinate transformation. The camera model and the principle of the SVD calibration method were presented. Using an observability analysis, the optimal calibration models for the star combinations and the singular value selection were discussed. The results of simulation with different calibration models were in accordance with our analysis, showing that the calibration model we recommended is an optimal model. Compared with the conventional AD method, the SVD method had similar accuracy; the convergence speed (the average calibration time per frame) of the SVD method, however, was significantly shorter than that of the AD method. This means the computational cost of the SVD method was smaller than the conventional methods, and so, the proposed method was more appropriate for star sensors with limited resources. Additionally, a field experiment showed that the SVD method can meet the requirements of on-orbit calibration in star sensors. Although the SVD method performed well, more tests in different situations should be conducted and further studies should be carried out, in order to learn more about this calibration method, as this is the first time it has been used for star sensor calibration.