Star Tracker Geometric Calibration Through Full-State Estimation Including Attitude

A star tracker calibration method using star images is presented in this paper. Unlike previous works, the proposed method estimates all parameters and the attitudes at once in a single least-squares formulation for the optimal calibration, which can be easily converted to a recursive estimation form. In addition, this paper presents a method to estimate the overall star tracker performance for attitude determination from the calibration results. Since the proposed method uses star images only, it can be applied to both on-orbit and ground star tracker calibration. The simulations show improvements in calibration performance about four times compared to the previous calibration method. The calibration experiments with actual star images are conducted to test its application.


Introduction
A star tracker is one of the most accurate sensors for spacecraft attitude determination [7,8]. Its typical accuracy ranges from 10 arcsec to sub-arcsec, depending on its size and purpose. Since the star tracker is an optical sensor that estimates geometric star vectors from the images captured by a star camera, the geometric defection called camera parame- ters, such as lens distortion and misalignment, can affect the sensor measurements. Thus, the calibration for the camera parameters is one of the most important aspects that dictate the star tracker's performance.
The star tracker calibration is a process of identifying the factors that make measured star positions different from their actual position. The types of calibration are divided into 'on-ground' and 'on-orbit' calibration according to the purpose, and also can be categorized as attitude dependent (ADC) and attitude independent (AIC) calibration according to the principle [4,16]. The ADC methods generally define the star measurements as the pixel points for star centroids on the detector coordinate. ADC technique needs the sensor's attitude information to compare the measurements with the actual star vector on the different coordinate system. On the other hand, The AIC method eliminates the necessity of attitude information by defining the angle between two stars as the measurements [4].
Several star tracker calibration methods have been proposed in recent decades. The on-ground calibrations typically have been performed by using the ADC methods. A fair number of the on-ground calibration uses some devices such as a gimbal or collimated beam in a laboratory-based on the estimated misalignment of the devices [9,14,19]. Wei et al. used a nonlinear least-squares algorithm to obtain optimal solutions for 14 parameters in a star sensor calibration process [18]. However, a coupling issue occurred between the intrinsic and extrinsic calibration parameters, including correlations between calibration parameters and misalignment of the rotary stage [11]. Although some papers proposed decoupling approaches [3,20], these are not free from coupling issues when using a gimbal mount for calibration. Instead of the artificial experimental setup, the natural starlight can also be used for ground star tracker calibration [5], which has advantages in costs as this does not require expensive equipment.
These approaches are also used in on-orbit calibrations with star images; however, obtaining an accurate attitude becomes an issue in this case. This is a "Chicken or the egg" problem: good attitude estimates are needed for accurate star tracker calibration, but accurate calibration parameters are needed to precisely estimate the attitude. As star trackers are the most accurate attitude sensors on most satellites, it is practically difficult to obtain an accurate attitude without them. Several studies address this issue and many use two or more steps to estimate the attitude and parameters separately. Pal and Bhat first estimated the attitude with some of the parameters before estimating the others [13], with a similar approach taken by Geng et al. [22]. Li et al. used the QUEST algorithm to obtain the attitude and calibrate the other parameters separately in multiple steps [6]. These multi-step approaches have been used in practice, but their convergence and optimality have not yet been proven.
The AIC method is another approach to tackle the "Chicken or the egg" problem. It has been widely used in the literature as well as in practice, especially for on-orbit star tracker calibrations [2,10,12,23]. As this is free from the spacecraft attitude, the attitude determination error does not affect the convergence of the calibration parameters, and the algorithm becomes simple and convenient. However, this raises a new question concerning its optimality. As mentioned in [2], the AIC performance depends on the star-pair selection. It may be possible to find an optimal way of selecting pairs, but this has not been reported.
This paper presents a low-cost and straightforward onground star tracker calibration method as the ADC. Unlike other ADC methods, Not using any device to control or give the attitude, a star tracker which takes the night sky alone is assumed. We derive a complete set of derivative equations for the measurements with respect to the calibration parameters and the attitude quaternions. For the full-state estimation, the measurement is defined as the star-centroid coordinates in detector space. Then, the weighted least-squares algorithm is applied to estimate the attitude quaternions of each star image as well as the camera parameters at once. As it uses the full-state estimation, the optimality of the estimates is guaranteed in the linear-estimation point of view. The standard deviations of the calibration parameters are compared with the reference study, which can be told a typical AIC method. The attitude determination performance of the calibrated star tracker is discussed with the results from simulations and outdoor experiments.

Notation
This section briefly defines the mathematical notation for quaternion and vector operations used in the subsequent sec- and a quaternionq is defined as where q is called the vector part, and q 4 is called the scalar term. Here the attitude is represented by the quaternion, which is defined as (2) where The unit vectorê is the axis of rotation, and θ is the angle of rotation. So, the attitude quaternion always satisfies the unity norm constraint An overbar is used to denote quaternions in this paper. Letp, q, andr be quaternions and let T p , T q , and T r be the corresponding direction cosine matrices (DCM). The quaternion product operator ⊗ is defined such thatr =p ⊗q corresponds to the product T r = T p T q . The skew-symmetric matrix that satisfies the cross product v × b = [v×]b is given by for a given 3-dimensional vector v. Using the notation defined above, a DCM corresponding to an attitude quaternionq can be calculated as The A(q) denotes the DCM corresponding toq. From (6), it is inferred that bothq and −q represent the same DCM. This meansq = −q in the sense of rotations. The inverse quaternionq −1 ofq is defined as a quaternion that satisfies q −1 ⊗q =q ⊗q −1 = 0 0 0 1 T .
From (2), (6), and (7), the inverse quaternionq −1 is obtained as In addition, the estimate of a state x is denoted by the hat symbol asx.

Camera and Star Measurement Model
The camera model is an eight-parameter inverted version of the model proposed by Wang et al. [17] and used by Enright et al. [2]. Unlike [22], this model does not include tangential distortions, but its effectiveness was empirically verified in [1,2].
The eight-parameter model is well-described in [2,17], so this paper only defines the notation and presents the equations, which are nearly identical to [2] with minimal comments for the readers' convenience.
A star measurement from a star tracker is the coordinate of a star centroid on the detector which is denoted as where m c is the column coordinate and n c is the row coordinate. The units of m c and n c are in pixels. The camera model has eight parameters to correct for lens distortion and alignments, which are m 0 , n 0 , f , b 2 , b 4 , a 1 , a 2 , and g y . The m 0 and n 0 define the optical center in vector form as The f is the focal length of the lens illustrated in Fig. 1a.
The b 2 and b 4 are the radial distortion coefficients in vector form as These correspond to the radial distortion effect as in Fig. 1c. The linearized rotation of the optical axis, shown in Fig. 1b, is represented by a 1 and a 2 as Fig. 1d shows g y , which is a relative skew factor between the pixel spacing in the column and row directions.
The inverted camera model reconstructs the 3-dimensional unit star vector in the detector frame s D from the measure- where is the pixel pitch of the detector. Applying an optical-axis rotation correction with a 1 and a 2 gives the displacement from the optical center as where Then, the radial correction factor B is calculated as where With the radial correction factor and the rotation-corrected coordinate, we construct the star vector in a pin-hole camera model with a focal length f as The unit star vector is obtained by normalizing r D as The star measurement model, which is the inverse of (13)- (20), is a star centroid coordinate on a detector as calculated from a unit star vector. For a given unit star vector in a reference frame s r , the star vector in the detector frame is calculated as whereq is the attitude quaternion from the reference frame to the detector frame. Then, r D is calculated with a given focal length f as With r 2 defined by (17) becomes We determine B by solving the equation where B is numerically calculated using the simple Newton-Raphson method. With B, the parameters of U and V are obtained as From (15) and (16), Finally, the star measurement m is written as where ν i j is the additive measurement of the j-th star in the i-th image noise given as The partial derivatives of m with respect to the eight parameters and attitude is derived to apply a nonlinear leastsquares method. Though complicated, it is possible to derive the derivatives in the form of

Camera Parameter Calibration with Attitude Estimation
Camera parameter calibration is accomplished using a nonlinear weighted least-squares approach with iterations. We assume that we have m star images and each image has n i stars where i = 1, 2, . . . , m. For each iteration, we calculate the measurement error δm i j for the j-th star in the i-th image as where m i j is the measured star centroid andm i j is the expected value on the detector, which is calculated from (21) to (29) using the current state estimates of (q i ,ô,f ,b,â,ĝ y ).
The derivatives of the measurement error of (31) allows formulating a linearized measurement error system as Equation (33) is a typical linear system of equations δy = H δx + N; thus, δx can be estimated optimally using a weighted least-squares as where With δx obtained from (34), the attitude estimates are updated aŝ for i = 1, 2, . . . , M. The camera parameters are updated aŝ o + = δo +ô, The superscript + means the iterative estimate on the nonlinear least squares algorithm.

Initial Guess of The States
The nonlinear least-squares presented in the previous section runs with the determined state estimates, which means that the first iteration requires an initial guess. To prevent the singular matrix inverse, a non-zerob 0 is required for the first iteration. This can be a problem if the lens is so well-made that the radial distortion is nearly zero, but this does not often occur in real applications. For the initial attitude guess, the q-method or the QUEST algorithm can be used with the initial camera parameters to determine the initial states for the first iteration [15].

Parameter Normalization
As the orders of magnitude for the parameter values can be different from each other, the (H T R −1 H ) in (34) can be numerically ill-conditioned, and the matrix inverse becomes unstable. Enright et al. proposed an effective method to overcome this issue using parameter normalization [2], which is used here.

Attitude Estimation Performance
The goal of a star tracker is to estimate the attitude, where the estimation accuracy is of interest to the users. With the proposed calibration method using a weighted least-squares, the covariance of the state estimates is given by P q 1 q 1 · · · P q 1 q m P q 1 o · · · P q 1 g y . . . . . . P q m q 1 P q m q m P oq 1 P oo . . . . . .
where P q 1 q 1 , . . . , P q m q m are the covariances of the attitude estimates from the calibration. However, these do not truly represent the attitude determination error of a star tracker in practice. Once a star tracker is installed on a spacecraft, the attitude determination performance is measured by comparing the attitude outputs to other sensors such as another star tracker, payload camera, or gyroscopes. Thus, the common, bias-like, slowly varying errors are blended into its misalignment, and only the high-frequency errors are measured on-orbit. Apart from P q 1 q 1 , . . . , P q m q m in (38), a new metric that represents the measurable performance is derived here. We consider an approximated bias attitude quaternion vector of the attitude estimates in the calibration as with an assumption of a large m. Then, the covariance of the attitude bias can be approximated as The practically measurable attitude error without bias, δp i , is calculated as and its covariance is given by where Finally, the expected measurable error covariance is given as This new metric that estimates the practical attitude performance of a star tracker, P pp , is numerically verified in the following section.

Simulation Study
The  [2]. For the additive measurement noise of˚i j in (30), σ i j x = σ i j y = 0.5 pixels is assumed as modern star trackers are capable of star centroid determination with sub-pixel accuracy.

Camera Parameter Calibration
This section presents the calibration performance of the proposed method. The camera parameters are estimated from simulated star images, whose star centroids are generated randomly instead of from a star catalog. Though the simulated approach may be different from actual star patterns, this gives a brief insight into the performance of the proposed calibration method with respect to the number of captured stars. One thousand images are used for each calibration, and additional stars are added gradually as the number of stars increases for consistency. For example, the image of 6 stars includes all the stars from the 5-star image plus one additional star. Figure 2 shows the standard deviation of the camera parameter estimates to the number of captured stars in a given image, which is calculated from (38). For the reader's convenience, the norms of [σ m o , σ n o ] and [σ a 1 , σ a 2 ] are plotted in Fig. 2a, e, respectively. To show the effectiveness of the Fig. 2 Standard deviation for the camera parameter calibration process using 1000 star images full-state estimation approach, the results are compared to those from the AIC method proposed by Enright et al. [2], which uses the identical camera model. As seen in Fig. 2, the standard deviations from the proposed methods are approximately four times smaller than those from the AIC method. This shows that the proposed method outperforms AIC due to the optimality of the full-state estimations. Figure 3 shows the average standard deviation of the quaternion estimates calculated from Some readers may think this is against common belief as the standard deviation is too large compared to actual star trackers. Generally, a modern star tracker gives about 10 arcsec and 100 arcsec accuracies in the cross-boresight and boresight directions, respectively. This issue is discussed in Sect. 6. The practically measurable attitude determination performance is different from the raw attitude covariance P q i q i . Thus, the value of σ qi,avr does not represent actual

Attitude Determination Performance
As discussed in Sect. 6, the practical attitude covariance can be calculated from (44). Monte-Carlo simulations are conducted to verify the new approach, and the results are compared to the standard deviation from (44). In each case,   100,000-star images are generated with perturbed parameters, and the QUEST algorithm estimates the attitude with the calibrated parameters as estimated from 1000 simulated star images. Then, the attitude bias is calculated by averaging the quaternion estimates, and the standard deviation of the 100,000 sample attitude is calculated with respect to the bias. Figure 4 shows the standard deviation from (44) as σ pp , and the Monte-Carlo simulations give σ from the simulation. In this result, the differences between σ from the simulation and σ pp are less than 0.5%, which is negligible considering the nature of Monte-Carlo simulations. The authors believe that this bias-free error is what was observed with real hardware experiments.

Outdoor Experiment
Outdoor experiments are conducted to test the calibration performance of the proposed algorithm. The ArduCam, which is a low-cost, open-source camera module, and the Raspberry Pi 4 are used for the experiments. The main detector of ArduCam, AR0135AT, is CMOS image sensor, 1.2MP,  Table 1 gives the key specifications of AR0135AT, and the initial values of the calibration parameters. The star images are taken every 3 s during the night, and 10 to 15 stars are captured in each image. The full-state parameter calibration is performed with 1000 images and 11,985 star measurements to estimate the eight calibration parameters and attitude. The centroid coordinates of the stars m are calculated using the star identification algorithm. Figure 5 shows all centroid coordinates for the 11,985 stars as captured on the detector. As the star camera is fixed on the ground, the stars flow in one direction over time due to the Earth rotation. 3000 s was sufficient to obtain distributed measurements on most of the detector area.
To show the calibration performance, the angle between the restored star vectors in the detector frame, s D , and the expected star vectors, A(q)s r , is calculated as for the j-th star in the i-th image where s r is the unit reference star vector in the inertial frame from the star catalog. Figure 6 shows the standard deviations σ θ i calculated with where i = 1, 2, . . . , 1000. The maximum and minimum values are 8.3 and 2.7 arcsec, respectively. Despite using commercial off-the-shelf (COTS) star camera modules, the experimental results demonstrate the nominal performance shown in Sect. 7.2 and illustrate that the proposed method can be effectively applied to actual star tracker calibration procedures.

Conclusions
A new star tracker parameter calibration method with star images is presented. The proposed approach is considered an attitude-dependent calibration (ADC) approach since it requires attitude estimates for the calibration. A complete set of partial derivative equations for the measurement errors is derived, which makes it possible to estimate full-states in a single step. The calibration performance of the proposed method is demonstrated using numerical simulations and compared with the result of the typical AIC method. Figure 2 indicates that the standard deviation of each parameter with the proposed method is under a fourth of that with the AIC method. Also, a practical approach for the evaluation of the calibrated star tracker's attitude determination performance is proposed using the covariance matrix of parameter estimates. When the evaluation results are compared with real attitude errors from the Monte-Carlo simulations, the difference appears to be less than 0.5%. The outdoor experimental results demonstrate the feasibility of the proposed calibration method even with low-cost COTS hardware. Though not presented here, the derivative equations can be used for on-board star tracker calibration with recursive estimation approaches or even for spacecraft attitude determination with a Kalman filtering method directly.

A.1 Derivation of Partial Derivatives of A Star Measurement Equation
This section derives the derivative equations for the measurement m with respect to all states to be estimated in the form of in a step-by-step approach. The first step is to take the derivative of (29) without additive noise ν i j .
From (21), it is known that the derivative of s D to the attitude quaternion is given as where dq is the vector part of a small delta-quaternion δq according to [21] as δq =q ⊗q −1 dq 1 .