Treatment of Extended Kalman Filter Implementations for the Gyroless Star Tracker

The literature since Apollo contains exhaustive material on attitude filtering, usually treating the problem of two sensors, a combination of state measuring and inertial devices. More recently, it has become popular for a sole attitude determination device to be considered. This is especially the case for a star tracker given its unbiased stellar measurement and recent improvements in optical sensor performance. The state device indirectly estimates the attitude rate using a known dynamic model. In estimation theory, two main attitude filtering approaches are classified, the additive and the multiplicative. Each refers to the nature of the quaternion update in the filter. In this article, these two techniques are implemented for the case of a sole star tracker, using simulated and real night sky image data. Both sets of results are presented and compared with each other, with a baseline established through a basic linear least square estimate. The state approach is more accurate and precise for measuring angular velocity than using the error-based filter. However, no discernible difference is observed between each technique for determining pointing. These results are important not only for sole device attitude determination systems, but also for space situational awareness object localisation, where attitude and rate estimate accuracy are highly important.


Introduction
To improve attitude determination performance, statistical estimation and filtering techniques are adopted by the space system designer. The most prevalent filter is the Kalman Filter, proposed by Rudolf Kalman during the 1950s [1]. The theory was developed and applied to space flight through a series of NASA reports [2,3] during the early 1960s. The earliest known study to apply the Kalman Filter for the attitude determination problem was in 1970 [4], where Farrell was one of the first to recognise the merit of the technique.
Since the Apollo era, the Kalman Filter and its modifications have been proposed for a variety of attitude determination systems, each adopting a different set of attitude determination sensors. The work by Farrell studied the crude use of sun sensors and magnetometers [4]. One of the first uses of a star sensor was used with a gyroscope by the Aerospace Corporation [5], applying a discrete Kalman Filter.
During the space race, satellite applications in low and medium Earth orbits were discovered. Since many missions required only an Earth-pointing attitude, not experiencing many complex dynamic motion routines, the gyroscope was neglected, and the attitude calculated using orientation only sensors. Gai in the 1980s [6] attempted to use a star sensor that, at the time, could only receive a single star measurement at a fast enough rate for attitude estimation. The attitude accuracy was poor, applying a simple finite difference-based approach and assuming negligible noise.
Spacecraft torques in the control system loop, as well as external torques caused by the space environment such as solar wind or aerodynamic drag [7,8], have also been only take place if a control torque is acted on the satellite. Thus, the expected angular acceleration is already known by the system.
The Calculated Reference of Stellar System (CROSS) star tracker is used as a test platform in this article. CROSS is a Wide Field of View (WFOV) star tracker being developed by the University of Sydney that seeks to be a high performing, versatile and competitive attitude determination unit for spacecraft designers [30][31][32].
In this paper, the additive and multiplicative EKF are applied to evaluate and determine the best approach to the gyroless star tracker attitude estimation problem. The novelty and originality of the work are summarised by: 1.
Novel implementation of the multiplicative EKF for the case of the gyroless star tracker.

2.
An in-depth, independent analysis to the state-of-the-art additive EKF for the gyroless star tracker by [17,18], comparing to the multiplicative EKF approach as well as an unfiltered approach (i.e., Linear Least Squares (LLS) only).

3.
Analysis employs both simulation and real night sky testing. Provides a case example of the novel approach of star tracker testing by [33].
The significance of the research is summarised as: • Enabling more compact and lower power attitude determination systems. • Improving attitude estimation accuracy and precision, a topic of interest for not only satellite attitude control but SSA applications. • Reducing bias to attitude and rate estimates.
The paper is structured by firstly outlining both EKF models. It then implements each model by simulation and real night-sky testing, as well as comparing to a static LLS approach, to discuss and evaluate the best approach. Theoretical arguments are also considered in this discussion. Recommendations of the best approach are made at the conclusion.
Given that gyroscope bias is now removed, and the EKF is estimating the angular velocity state, additive and multiplicative filters are relabelled as the state and the stateerror EKF. This naming change highlights that the multiplicative approach estimates the state error, rather than directly the state.

Model and Methods
This section introduces the conventions adopted to describe the spacecraft attitude. A sensor model is then introduced for the star tracker, which includes star position error considerations. The section closes by describing the two approaches to the EKF implementation, which contains a modified version of implementations in [15,17] to treat the performance of the gyroless star tracker case.

Attitude Definitions
Attitude is typically expressed as a matrix that transforms a vector from one coordinate frame to another. In the case of spacecraft attitude, the transformation typically relates the sensor or spacecraft frame to a celestial body frame, such as for the Earth or Solar System centre. The celestial frame will be denoted by r, and the sensor or spacecraft frame will be denoted by b. The relation between each frame and the attitude may be written as, where A is the attitude matrix.
To aid in the representation of attitude, consider an axis e for a rotation to act on. The angle of rotation may then be denoted by θ. The rotation of an arbitrary vector x about the rotation vector is illustrated in Figure 1. The attitude matrix may be expressed as a function of the rotation vector and angle, A = A(e, θ). It is more convenient, however, to express the rotation as a quaternion. The quaternion q may be expressed in terms of the rotation vector as, The quaternion must satisfy the unity constraint |q| = 1. Maintaining this constraint in estimation and filtering can prove problematic, as introduced in Section 1.
The attitude matrix and quaternion are related by, Derivations of attitude kinematics are not included in this description, as they are extensively written in the literature [15,34], but are instead stated. The time derivative of the attitude quaternion may be expressed in terms of the angular velocity, ω = ω x ω y ω z T , by,q where Ω is a matrix representation of the tensor product operation between the angular velocity and quaternion, A linearised form of the state transition over the duration of a time step ∆t can be derived. Assuming a constant angular velocity over the time step and that the angular rotation is small, the closed form may be expressed as, where Φ qq is the quaternion transition matrix and I 4 is the 4 × 4 identity matrix.
The angular velocity may also be estimated assuming it is constant over the time step. The quaternion rate is expressed as, The quaternion rate may then be estimated using this equation if the time step is justifiably small. Equation (4) may then be rearranged in terms of the angular velocity, using a similar approach to Equation (6), where,

Star Tracker Model
The star tracker captures images of the stars using a standard camera and lens assembly. The image is then analysed to identify stars against a known catalogue. Using the measured unit vectors of stars in a sensor frame and the known unit vectors from a catalogue in a global frame, an attitude may be calculated.
The measured star image is represented by the sub-pixel coordinates of the stellar source (α, γ). The sensor frame unit vector is related to the sub-pixel coordinates by, where f is the focal length. The body and celestial frame unit vectors are then related by Equation (1). The reverse relation to Equation (10) uses the unit vector normalisation condition, |b| = 1. So, the body unit vector may be calculated by, A common measurement error model for the measured image coordinates α and γ is [11,14], This model is appropriate for sensor measurements that are no greater than 15 • from the boresight. The star tracker measurements considered in this paper are 10 • from boresight. The eigenvalues and eigenvectors can be determined from the matrix to know the maximum error of the image.

Extended Kalman Filter Model
Two approaches to applying EKF to attitude determination systems are considered in the literature, the state-based and state error-based, as are first introduced in Section 1.
A stand-alone star tracker will be initially introduced, using the state-based approach. The popular error state variant will then be discussed. The angular velocity cannot be directly measured by the star tracker. An estimate by finite difference, applying Equation (8), is adopted. The EKF is then compared to static estimation by LLS, using the popular Davenport q method described in [15].

State-Based Estimation
Using a stand-alone star tracker in the attitude determination system may estimate both the attitude as a quaternion and the angular velocity. This section adopts approaches initially proposed in [17].
The state variables may be expressed as, where q = q 0 q 1 q 2 q 3 T and ω = ω x ω y ω z T . The state transition, or propagation update, may be modelled using, where φ is the state transition function and σ is the additive Gaussian noise. Using Equation (6), the state transition function may be expressed as, The EKF is an approach to using non-linear system models in the linear Kalman Filter technique. It does this by a process of linearisation, as has already been applied to Equation (6). The overall state transition matrix is represented by, where Φ qq is already known from Equation (6). The partial derivative of q(t + ∆t) with respect to ω(t) is, and thus for i = x, y, z, where, Since ω is assumed to be constant over the time step, Φ ωq = 0 and Φ ωω = I 3 .
In a similar form to the state or propagation update, as in Equation (14), the measurement update is given by where h is the measurement function and v is the measurement noise, modelled as a zero-mean Gaussian noise. The measurements of the star tracker are given in sub-pixel coordinates calculated from the star source brightness distribution. The expected measurement h(x) is derived from the identified stars. The star catalogue coordinates are given in the inertial frame r. These are transformed to the star tracker sensor frame by the attitude matrix A using Equation (1). The ith star locations known in the derived expected image with coordinates α i and γ i are obtained using Equation (10).
A process of linearisation is also required for the measurement update, and so the partial derivative of the measurement function with respect to the state variable produces the measurement matrix, where n is the total number of stars in the field of view. Using the chain rule, the partial derivatives with each state variable may be calculated using, So, the partial derivative of the measured sub-pixel star locations α and β are, The partial derivative of b with respect to the angular velocity ω is, and with respect to the quaternion q is, where for each quaternion the partial derivative of the attitude matrix is, The complete methodology of the EKF model is summarised by Table 1. The + and − superscripts indicate the pre-and post-measurement update, respectively, used by the state vector x and covariance matrix P. The precidicted state vector is expressed by z.
Given the risk to ill-conditioning of the covariance matrix P owing to non-linearities and the re-normalisation of the quaternion, P is checked for positive definiteness by a Cholesky Factorisation attempt. If the attempt failed, the covariance matrix is reset by negating all diagonal negative terms and setting all non-diagonal terms to zero. In the reported results of this work, this reset was never necessary.

State Error Estimation
For the stand-alone star tracker, the state error estimation approach uses the formalisation from [15]. However, the approach is made novel by replacing the gyroscope bias state by a direct estimate of the angular velocity, making various modifications to the model to reflect this state estimate change. The state variables may be expressed as in the state-based approach, However, the filter will not act directly on the state, but the state errors. The state error vector is expressed as, where δθ = δθ x δθ y δθ z T and δω = δω x δω y δω z T . States The propagation update is rewritten as, The linearised state transition matrix is then, An alternative approach is adopted by [15] to derive the state transition matrix. Using the linearised state space estimate equation for the rate of state change, where F(t) is the Fisher information matrix, which provides a linearised matrix operator on the propagated state via, It is known from [15] that the rate of attitude error change is related to the attitude error and angular velocity terms by, where [ω×] is the matrix operator of the cross product, written as, So the Fisher information matrix is expressed as, where I 3 is the 3 × 3 identity matrix. The Fisher information matrix and the transition matrix is related by the expression, Considering each side and maintaining equivalence, the transition matrix constituents may be expressed as, The state update is still performed using the approach in Equation (6). The angular velocity is assumed to be constant over the time-step. The measurement update uses a similar form to the state-based approach, expressed as, However, instead of considering the measured pixel coordinates, it considers the measured star unit vectors directly. The unit vectors are calculated using the image measurements α and γ using Equation (10). The measurement function is defined as, The linearised measurement matrix is solved by considering the partial derivatives of the measurement function with respect to the state error vector, This is solved directly in [15], using the chain rule and so the state error vector is expressed as, The attitude error represents the difference between the true and estimate quaternion by the relation, assuming the attitude error is suitably small, where Ξ was defined in Equation (9). Similarly, the angular velocity is related by, The methodology of this approach is summarised in Table 2. States

Results and Discussion
This section presents the results of analysing each estimation model. The results include simulations of the star tracker, as well as real night sky testing using the CROSS star tracker [30]. In the real night sky test, the Earth's rotation is used for dynamic testing under a constant rotation. This unfortunately constrains performance analysis to precision only and not accuracy, but still allows meaningful conclusions to be derived.

Simulations
Both models presented in Section 2.3 are simulated and compared. The error models described in Section 2.2 are used to produce representative measurement errors. The model errors are the same for each filter approach considered. As mentioned, results are also compared to a static least square estimator, to establish a baseline for performance.
Each model starts at an initial attitude of q 0 = √ 2 2 0 0 1 1 T . The star tracker is set initially at a constant rotation in the y-axis, or pitch direction, with a speed of 36 arc-sec/s. A noise of 1 × 10 −4 pixels is set at the centroid of each star. The simulations are run over a period of 90 min, with measurements captured once every second, i.e., 1 Hz. The star tracker is set to a square field of view of 20 • and a stellar magnitude maximum of 5.0. The attitude does not considerably improve with more than 8-10 stars [15,31], so a higher magnitude limit is not necessary.
The results are presented in Figure 2 and Table 3. The best performance for angular velocity is found for the state-based EKF compared to static LLS and state error. A slight improvement from state-based and a large improvement from static when comparing the speed error and standard deviation of Table 3. However, in terms of the attitude, the results are not so distinct. The yaw error is clearly the worst, being the direction perpendicular to the boresight. In all directions, a slight improvement might be noted in the yaw error by adopting a state approach. The state error approach corrects by small increments at each time step. With significant noise in each stellar measurement, the restricted error term may not be enough to correct the noise.

Real Night Sky
The EKF approaches may also be compared and further analysed through samples of real night sky images. The EKF is also compared to a LLS based approach that does not consider any prior information.
The test is performed by fixing the CROSS star tracker to a tripod and placing it with the camera facing approximately at the zenith. The camera's horizontal field of view is approximately 19.9 • and the lens focus is located a distance to the sensor that intentionally blurs the stars in the image, aiding sub-pixel centroid analysis. After each captured image, each star is processed by measuring the unit pointing vector and identifying it to an equivalent catalogue of vectors. These vectors are then employed in the LLS or EKF. An illustration of the assembly is presented in Figure 3. As the Earth rotates, the stars move across the camera's field of view at the known speed of Earth's rotation. To compute a constant attitude throughout the test duration, the Earth's rotation, as well as empirically known nutation and precession angular fluctuations, may be corrected through a series of matrix operations to the measured attitude matrix. The residuals and precision of the measured attitude may then be calculated, providing an evaluation of the star tracker's performance. The corrections to the attitude matrix A may be given by, where R, N and P are the Earth rotation, nutation and precession matrices, respectively. Further description of each correction is described in [33]. The spread is calculated by fitting a linear trend line to the data set, measuring the goodness-of-fit, and transforming this to an error variance. Each image was captured on a hill in the Inner West of the City of Sydney, Australia, on 10 May 2021. The sky contained some light pollution, but the astronomy weather forecast website, Meteo Blue, reported the presence of no planetary bodies in the night sky, with good visibility conditions and little to no cloud cover. It is considered that atmospheric noise may lead to worse performance, but should not significantly affect comparisons of each approach. The camera board is a FLIR Blackfly S containing a Sony IMX264 image sensor [35], 2448 × 2048 px, and is accompanied by a Scorpion Imaging lens [36], with an approximate diameter of 17 mm and 2/3 sensor size. Each image is captured with a 15 s interval and an exposure time of 0.15 s.
Each approach is compared, also considering a shorter time difference between each image, and so a faster rotation speed. This is accounted for by reducing the time factor between each correction. Tthe measured attitude in each direction, yaw, pitch and roll, and the angular velocity is considered. The angular velocity of the LLS approach is calculated by taking the simple difference between the measured Euler angles and dividing by the time step.
It should be noted that the dependence on update rate has minimal influence on performance. Farrenkopf [15,37] calculated the steady-state error in the state error EKF for a gyroscope and star tracker Kalman filter, where, where σ v is the gyroscope noise. For the gyroless case, σ v may be treated as a process noise in the propagation step. The dependence on the time update ∆t is the same though, where it is the quartic root.
The results for the attitude and angular velocity are presented in Figure 4. The results are also tabulated in Table 4. The error is calculated for the attitude by fitting a trend line and considering the residual. The error in the angular velocity measurement is calculated by the residual to the known Earth rotation. The standard deviation represents the precision of the data set series to the trend line. The attitude in Figure 4 shows a slow drift from the initial attitude. This may be due to external factors such as the wind or a general slip on the tripod, and is highlighted by the data spread being greater than the slip or drift speed. The speed results are calculated after a convergence time of 60 s.
Comparing the attitude results of the EKF to the LLS static approach, no distinct improvement is derived from using a filter. Differences between the state error and statebased approaches are also not so clear. This conclusion was also demonstrated in Section 3.1.
The convergence time of the angular velocity is much longer for the EKF, but especially for the state-error approach. This is reasonable, as the state error correction is restricted to only small increments. However, once converged, as seen in Table 4, the EKF precisions are much higher than the static LLS. It may be concluded that for high spin rates, a Kalman filter is more suitable to estimate angular velocity.
Overall, the state approach is better in terms of angular velocity, as in Table 4. This result is suspected to be due to the restricted magnitude of the state error, which is insufficient to correct for propagation approximation error. The difference in attitude ac-curacy and precision between state-based and state error techniques after convergence in Figure 4 and Table 4 is not strong enough to significantly prefer a state-based over a state error-based approach. It could be argued that the star tracker could also be tested against more erratic behaviour with a changing angular velocity. However, such phenomena are almost always caused by a deliberate torque in the spacecraft attitude controller, which can be directly measured and fed into the attitude determination system. This would improve results. Higher speeds also are not appropriate, as significantly smeared images at high spin rates cannot successfully identify an attitude [38].

Conclusions
This work compares state and state error-based Kalman filter formulations to determine the best approach for the gyroless star tracker. Each approach is traditionally known as additive and multiplicative EKF, but is relabelled for the gyroless case where the spin rate is a directly estimated state. This treatment is important for new applications of star trackers in SSA, where spacecraft spin speeds are an important parameter, as well as traditional attitude determination.
Approaches consider both simulation and real night sky analysis. The real night sky case uses a novel technique for testing in the absence of a known truth orientation. The state-based Kalman approach was determined to perform the best in each case. It is not restricted by the error term and reset step. This approach may encounter ill conditioning of the covariance matrix and requires renormalisation, increasing the error margin. However, the magnitude of these errors are considered minor and were not noticeable in the results reported. Moreover, resetting the covariance matrix proved unnecessary for the analyses reported here.
Further work will see implementation of this Kalman filter on an upcoming satellite mission. It will also see implementation of the state estimate in direct SSA applications.