Nomenclature

$a({M_V})$

star magnitude-dependent constant

$AMP(m,n,{x_0},{y_0})$

total energy (amplitude) received at the star tracker pixel location (m,n)

$C_B^{STi}$

direction cosine matrix defining the body frame to the orientation of the ith star tracker frame

$C_{STi}^{ECI}$

direction cosine matrix defining orientation of the ith star tracker frame to the Earth-centered inertial (ECI) coordinate frame

d

detector length

$\Delta {q_{xi}},\Delta {q_{yi}},\Delta {q_{zi}}$

deviation of the ith star tracker along x, y and z axes

$\Delta t$

integration step

$erf(x)$

Gaussian error function $ \left(= \dfrac{2}{{\sqrt \pi }}\int\limits_0^x {{e^{ - {p^2}}}dp}\right) $

(m,n)

pixel location at the mth row and the nth column of the captured image frame

$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{ECI \to B}^{STi}$

star motion rates from the ECI frame to the body frame

$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{B \to STi}^{STi}$

star motion rates from the body frame to the star tracker frame

$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{ECI \to B}^B$

spacecraft body rates from the ECI frame to the body frame

${\omega _{xi}},{\omega _{yi}}$

star motion rates across the star tracker focal plane

$PSF(x,y,{x_0},{y_0})$

point spread function at the star tracker’s pixel location (m,n)

$ratio({\omega _{xyi}})$

amplification factor from the star tracker cross-axis error to the boresight axis error

${R_{STi}}$

noise covariance matrix of the ith star tracker

${\sigma _x},{\sigma _y}$

Gaussian radiuses with respect to the star tracker’s pixel

T

total integration time period

$({x_0},{y_0})$

star location in the x-y plane during integration time

${\hat x_0},{\hat y_0}$

estimated star locations in the x-y plane

${x_{st}}$

true star position

2. Spacecraft attitude determination

Consider a spacecraft attitude determination subsystem which consists of a set of four single-axis Litef FOG gyros (4-for-3 redundancy) and two star trackers. Figure 1 shows the functional block diagram of our attitude determination design which involves two gyro stellar inertial attitude estimates (GS IAEs). Each GS IAE uses quaternion provided by the individual star tracker as shown and the sensed spacecraft body rates provided by 3-axis gyros to produce the best estimate of spacecraft attitude and body rates. GS IAE jointly employ gyros and star sensors to secure precise attitude information since star trackers will degrade their attitude accuracies when the spacecraft is subject to high-rate motion. The obtained estimates of spacecraft attitude and body rates are optimally combined in the attitude and rate fusion algorithm. The details of GS IAE implementation and the attitude and rate fusion algorithm can be found in our previous work [Reference Wei, Tsai, Yeh, Jan and Wu2]. Other related IAE designs using star trackers, gyros, and EKF can be found in Refs Reference Lefferts, Markley and Shuster1, Reference Wu and Hein4.

Figure 1. Attitude determination functional block diagram.

3. Star motion induced error and compensation

3.1 Star motion induced centroiding error

A typical star tracker, as shown in Fig. 2, focuses starlight onto an active pixel array converts the received photons to electrical signals, processes the signals, and provides time-tagged star position and star magnitude via I/O bus to the spacecraft control processor (SCP). For most star tracker designs, star image is purposely defocused so that its image will occupy several pixels across the detector arrays. As a consequence, sub-pixel resolution can be obtained in determining the star centre position from the defocused image using the so-called centroiding algorithm.

Figure 2. Star tracker with the active-pixel array.

An intuitive point spread function (PSF) to describe the Gaussian star image can be described as

(1) \begin{align}PSF(x(t),y(t),{x_0},{y_0}) = PSF(x(t),{x_0})PSF(y(t),{y_0})\end{align}

where $PSF(x(t),{x_0}){\rm{ = }}\dfrac{1}{{{\sigma _x}}}\sqrt {\frac{{a({M_V})}}{{2\pi }}} {e^{ - \frac{1}{2}{{\left( {\frac{{x(t) - {x_0}}}{{{\sigma _x}}}} \right)}^2}}}$ and $PSF(y(t),{y_0}) = \dfrac{1}{{{\sigma _y}}}\sqrt {\frac{{a({M_V})}}{{2\pi }}} {e^{ - \frac{1}{2}{{\left( {\frac{{y(t) - {y_0}}}{{{\sigma _y}}}} \right)}^2}}}$ . The total energy or amplitude received at the pixel location (m,n) is then given by

(2) \begin{align}AMP(m,n,{x_0},{y_0}) = \frac{{a({M_V})}}{{2\pi {\sigma _x}{\sigma _y}{T^2}}}\int\limits_0^T {\int\limits_{m - 0.5}^{m + 0.5} {{e^{ - \frac{1}{2}{{\left( {\frac{{x(t) - {x_0}}}{{{\sigma _x}}}} \right)}^2}}}dx} dt\int\limits_0^T {\int\limits_{n - 0.5}^{n + 0.5} {{e^{ - \frac{1}{2}{{\left( {\frac{{y(t) - {y_0}}}{{{\sigma _y}}}} \right)}^2}}}dy} dt} } \end{align}

Based on the received total amplitudes from each pixel, the estimated star locations can be computed using the centre of moment formula as follows

(3) \begin{align}{\hat x_0} = \frac{{\sum\nolimits_m {\sum\nolimits_n {mAMP(m,n,{x_0},{y_0})} } }}{{\sum\nolimits_m {\sum\nolimits_n {AMP(m,n,{x_0},{y_0})} } }},\,{\hat y_0} = \frac{{\sum\nolimits_m {\sum\nolimits_n {nAMP(m,n,{x_0},{y_0})} } }}{{\sum\nolimits_m {\sum\nolimits_n {AMP(m,n,{x_0},{y_0})} } }}\end{align}

By omitting the time index, as the pixel number m×n approaches to infinity and the detector length d approaches to zero, they become

(4) \begin{align}{\hat x_0} = \frac{1}{T}\int\limits_0^T {\int\limits_{ - \infty }^\infty {x \cdot AMP(x,{x_0}){\rm{ }}dx} }\ {\rm{ }}dt,\,\,{\hat y_0} = \frac{1}{T}\int\limits_0^T {\int\limits_{ - \infty }^\infty {y \cdot AMP(y,{y_0}){\rm{ }}dy} }\ {\rm{ }}dt\end{align}

In practice, only the finite number of detectors can be used in (3) to determine the star locations. One of the key design parameters in the current state-of-the-art star tracker is the selection of sigma values (or Gaussian radiuses) which spread star images over several pixels. Figure 3 illustrates four distributions of the received amplitudes over the one-dimensional detector array versus sigma values for a stationary star with ${x_0} = - 0.25$ pixel.

Figure 3. Amplitude distributions as function of sigma values.

As shown with the sigma value equal of 0.7 (in terms of the detector width), 99.74% of the star amplitudes (or energy) will be captured within five pixels. In this case, star locations with better than one tenth of sub-pixel resolution can be obtained just using 5×5 pixel arrays. However, when a star is moving across the detector during the star integration time period, the number of pixels needed to contain more than 99.74% of its received energy will increase. Figure 4 shows the resultant amplitude distribution over one-dimensional detector array for a stationary star with ${x_{st}} = 0.25d\;{\rm{and}}\;\sigma = 0.65$ (Fig. 14(a)) and a moving star shifting at the constant rate of 80 pixels/sec during 0.1 second integration time (Fig. 14(b)). The result indicates that more than 11 pixels are needed for this moving star to accurately determine its location. Although the centroiding accuracy can be retained by using a larger number of pixel arrays for a moving star, however, the position error due to detector noises will be increasing as well.

Figure 4. Effect of star motion on pixel amplitude distribution for the (a) stationary star, (b) moving star.

3.2 One dimensional pixel amplitude computation

3.2.1 Pixel amplitude calculation for stationary stars

Consider a stationary star located at the centre of detector array shown as in Fig. 5 and let ${x_{st}}$ denote its true location. According to (2), the amplitude of the centre pixel or pixel #0 can be written as

(5) \begin{align}AMP(0,{x_{st}}) = \frac{{a( {{M_V}})}}{2}\left\{ {erf({p_2}) - erf({p_1})} \right\}\end{align}

where ${p_1} = \dfrac{{ - 0.5 - {x_{st}}}}{{\sqrt 2 {\sigma _x}}}\ {\rm{ and }}\ {p_2} = \dfrac{{0.5 - {x_{st}}}}{{\sqrt 2 {\sigma _x}}}$ . For the ith pixel, its amplitude $AMP(i,{x_{st}})$ can be computed as

(6) \begin{align}AMP(i,{x_{st}}) = \frac{{a( {{M_V}})}}{2}\left\{ {erf({p_{2i}}) - erf({p_{1i}})} \right\}\end{align}

where ${p_{1i}} = \dfrac{{i - 0.5 - {x_{st}}}}{{\sqrt 2 {\sigma _x}}}\ {\rm{ and }}\ {p_{2i}} = \dfrac{{i + 0.5 - {x_{st}}}}{{\sqrt 2 {\sigma _x}}}$ .

Figure 5. One-dimensional pixel array.

3.2.2 Pixel amplitude calculation for moving stars

For a nonstationary star, ${x_{st}}$ will be time-varying during the integration time period on the ergotic hypothesis of the stationary random process, the statistic is calculated by averaging the amplitude of each pixel received over the integration time as

(7) \begin{align}AMP(i,{x_{st}}) = \frac{{a( {{M_V}})}}{T}\int\limits_0^T {\left\{ {erf({{\bar p}_{2i}}) - erf({{\bar p}_{1i}})} \right\}dp} \end{align}

where ${\bar p_{1i}} = \dfrac{{i - 0.5 - {x_{st}}(t)}}{{\sqrt 2 {\sigma _x}}}\ {\rm{ and }}\ {\bar p_{2i}} = \dfrac{{i + 0.5 - {x_{st}}(t)}}{{\sqrt 2 {\sigma _x}}}$ . There is no closed-form solution to solve for this integration; it can be realised by taking the following approximation with a sufficiently large N:

(8) \begin{align}AMP(i,{x_{st}}) \approx \frac{{a( {{M_V}})}}{N}\sum\limits_{n = 0}^{N - 1} {\left\{ {erf({{\bar p}_{2i,n}}) - erf({{\bar p}_{1i,n}})} \right\}} \end{align}

where ${\bar p_{1i,n}} = \dfrac{{i - 0.5 - {x_{st}}({t_n})}}{{\sqrt 2 {\sigma _x}}},{\rm{ }}{\bar p_{2i,n}} = \dfrac{{i + 0.5 - {x_{st}}({t_n})}}{{\sqrt 2 {\sigma _x}}},{\rm{ }}{t_n} = n\Delta t\ {\rm{ with }}\ \Delta t = \frac{T}{N}$ .

3.3 One-dimensional model

A simulation model was generated to (i) compute the amplitudes each pixel received according to (8), (ii) compute the star centroid position using (3), and (iii) compute the difference (the centrioding error) between the obtained star centroiding position and the star position at the half of integration time for a given star at various rates. Figure 6 shows the resultant centroiding error as function of the number of pixels used in the centre of moment calculation and non-uniformity noises with a stationary star, where $a( {{M_V}}){\rm{ = }}1$ and d = 27.78arc sec. For stationary stars, seven-pixel arrays are sufficed for accurately (less than 0.4arc sec) determining their positions with low noises (e.g. signal-to-noise ratio less than 0.003, 3 sigma) or without noise. For a star moving at the rate of 80pixels/sec, its centroiding error as functions of number of pixels used and non-uniformity noises is shown in Fig. 7. As shown, more than 11 pixels are needed to accurately determine its position in this case. Figure 8 shows the centroiding error as function of star motion rates using 13-pixel arrays for centroiding calculation. It shows that the centroiding error increases significantly when the star moving rate is greater than 0.8deg/sec with noise-to-signal ratio equal to 0.003 (3 sigma). The star motion rate having significant centroiding errors will be further reduced with the reduced pixel arrays for centroiding calculation as shown in Fig. 9 with the star moving rate 0.6deg/sec versus 0.8deg/sec when only 11 pixel arrays were used in the centroiding calculation.

Figure 6. Centroiding error as functions of number of pixels used and non-uniformity noises for a stationary star.

Figure 7. Centroiding error as function of number of pixels used and non-uniformity noises with a star with the moving rate 0.61728deg/sec.

Figure 8. Centroiding error as function of star motion rates using 13-pixel arrays for centroiding calculation.

Figure 9. Centroiding error as function of star motion rates using 11-pixel arrays versus 13-pixel arrays for centroiding calculation.

3.4 Adaptive EKF implementation

The following is a follow-up to our previously proposed GS IAE implementation [Reference Wei, Tsai, Yeh, Jan and Wu2], the noise covariance matrices ${R_i}$ of the i-th star tracker used for computing the Kalman gain matrices upon the star tracker updates are given by

(9) \begin{align}{R_i} = \left[ {C_{STi}^{ECI}} \right]\left[ {{R_{STi}}} \right]\left[ {C_{ECI}^{STi}} \right]\end{align}

where the $3 \times 3$ matrix ${R_{STi}} = diag(\Delta q_{xi}^2,\Delta q_{yi}^2,\Delta q_{zi}^2)$ . For the definitions of direction cosine matrices $C_{STi}^{ECI}$ and $C_{ECI}^{STi}$ see Ref. Reference Dong and Colibrys14. Design of EKF is well known we thus refer the interested readers to Ref. Reference Wei, Tsai, Yeh, Jan and Wu2 for the details of implementation.

Taking into account the star motion-induced errors appeared in the star trackers, we propose here an adaptive EKF (AEKF) implementation by specifying the following three time-varying functions when computing the noise covariance matrices to remedy potential attitude performance degradation using the traditional EKF design:

(10) \begin{align}\Delta {q_{xi}} = {f_x}({\omega _{xi}}),\;\;\Delta {q_{yi}} = {f_y}({\omega _{yi}}),\;\;\Delta {q_{zi}} = {f_z}({\omega _{xi}},{\omega _{yi}})\end{align}

where ${\omega _{xi}},{\omega _{yi}}$ are star motion rates $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{ECI \to STi}^{STi}$ across the x-y focal plane, which can be estimated using the measured spacecraft body rates ${\hat{\omega}^B}$ :

(11) \begin{align}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{ECI \to STi}^{STi} &= \left[ {\begin{array}{*{20}{c}}{{\omega _{xi}}}\\ \\[-7pt] {{\omega _{yi}}}\\ \\[-7pt] {{\omega _{zi}}}\end{array}} \right] = \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{ECI \to B}^{STi} + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{B \to STi}^{STi} \nonumber\\[3pt] &= \left[ {C_B^{STi}} \right]\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{ECI \to B}^B \nonumber\\[3pt] &= \left[ {C_B^{STi}} \right]\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over \omega } _{}^B \end{align}

Based on the star tracker specifications [15], the profiles of ${f_x}({\omega _{xi}})$ and ${f_y}({\omega _{yi}})$ are depicted in Fig. 10. The function ${f_z}({\omega _{xi}},{\omega _{yi}})$ is given by

(12) \begin{align}{f_z}({\omega _{xi}},{\omega _{yi}}) = ratio({\omega _{xyi}})\max \left\{ {{f_x}({\omega _{xi}}),{f_y}({\omega _{yi}})} \right\}\end{align}

where ${\omega _{xyi}} = \max \left( {{\omega _{xi}},{\omega _{yi}}} \right)$ and $ratio({\omega _{xyi}})$ is an amplification factor from the star tracker cross-axis error to the boresight error, which is described by Fig. 11. The ratio increases as the star motion rate increases, since the effective star separation angle will be decreasing accordingly. The formulation above resulting in adaptive feature of the EKF related to star motion rates.

Figure 10. Derived functions for $\;{f_x}({\omega _{xi}})$ and $\;{f_y}({\omega _{yi}})$ .

Figure 11. Boresight-axis error to cross-axis error ratio as function of star motion rates.

In the simulation model described below, we use the 11^th order polynomial fitted to (10) to model the star tracker motion-induced errors. However, for practical concern of the limited computational resource, a lower-order (i.e. 4^th order) polynomial is suggested when computing the noise covariance matrices for real-world implementation. The computation time of the 4^th order Runge-Kutta with the 4^th-order polynomial is faster than the 11^th-order polynomial; the simplified one is about two times faster.

4. Simulation results

To assess the GS IAE performance subject to the above derived star motion-induced errors, we create a simulation model, which consists of three major modules: (i) spacecraft attitude and body rate motion generation; (ii) sensor models for two STR’s and one 3-axis gyro; and (iii) GS IAE algorithms. Two representative bus body inertial rate profiles for a common micro-satellite orbiting in a circular sun-synchronous orbit at the altitude of 570km: one for nadir pointing mode with low body inertial rates as shown in Fig. 12 (roll in blue, pitch in green, and yaw in red) and the other for off-nadir video imaging mode with high body inertial rates shown in Fig. 13. Assume that the three-axis gyros are perfectly aligned with the spacecraft body frame. Two direction cosine matrices transforming STR#1 and STR#2 coordinate frames to the spacecraft body frame are given by

\begin{align*} {M_{STR\# 1toBody}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{ - 0.5407} & {0.4557}& {0.7071}\\ \\[-7pt] {0.8412} &{0.2929} &{0.4545}\\ \\[-7pt] { - 0.0000} &{0.8406}& { - 0.5417}\end{array}} \right],\,{M_{STR\# 2toBody}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{ - 0.5407} &{ - 0.4557} &{ - 0.7071}\\ \\[-7pt] { - 0.8412} &{0.2929} &{0.4545}\\ \\[-7pt] {0.0000} &{0.8406} &{ - 0.5417}\end{array}} \right]\end{align*}

Figure 12. Low body rate profiles used in simulation.

Figure 13. High body rate profiles used in simulation.

The corresponding star motion rate profiles appeared at STR#1 and STR#2 are displayed in Fig. 14. Table 1 lists the major gyro noise parameters as well as star tracker noises considered in the simulation study.

Table 1. Parameter settings for simulation study

Figure 14. (a) Star motion rates on STR#1 and (b) Star motion rates on STR#2 for the satellite mission under consideration.

Three cases were set up for IAE performance evaluation and comparison: Case 1 (low-rate case)-nadir pointing for simultaneous panchromatic (PAN) and multispectral (MS) imaging; Case 2 (high-rate case with EKF) off-nadir pointing for video imaging with the traditional EKF design; Case 3 (high-rate case with AEKF) off-nadir pointing for video imaging with the proposed AEKF design. Table 2 summaries the simulation results based on 1,000 Monte Carlo experiments in each case where the geolocation errors were calculated by using spacecraft altitude as well as the resulting body roll, pitch and yaw errors. All errors are expressed in a 68.27% confidence level.

Table 2. IAE performance with different GS IAE configurations

Case 2 shows the effect of star motion-induced errors on the GS IAE performance – the geolocation error increases from 32.64 to 33.32m when the spacecraft operates at the high-rate regions. Case 3 shows performance improvement using the AEKF design – the geolocation error reduces from 33.32 to 33.05m at high body rates. For this particular spacecraft body rate profile, one notices that the star motion rates in each axis are less than 0.82deg/sec (see Fig. 13) which shows limited benefit of using the AEKF.

To assess advantages of using the adaptive EKF in the GS IAE design, we generate a higher body rate profile as shown in Fig. 15 for a CubeSat mission which is going to perform video imaging with off-nadir pointing. Figure 16 exhibits the corresponding star-motion-rate profiles appeared at STR1 and STR2, in which the maximum star-motion-rate reaches 1.2deg/sec. Table 3 summarises simulation results of the above three cases using the CubeSat body rate profiles. Significant geolocation (or IAE) performance degradation (from 18.27 to 58.30m) due to star motion-induced errors observed from Cases 1 and 2. The advantage of using the proposed AEKF in the GS IAE design, when the spacecraft is subject to a much higher motion rate (greater than 1 deg/sec), is well demonstrated by the result of Case 3, which shows more than 40% performance improvement over that of Case 2.

Figure 15. High body rate profiles used in the simulations for the CubeSat.

Figure 16. (a) Star motion rates on STR#1 and (b) Star motion rates on STR#2 for the CubeSat program.

Table 3. IAE Performance with different GS IAE configurations for the CubeSat

5. Conclusion

In this paper, a general framework of mathematical model describing the star motion-induced error is analysed. Compensation technique as well as implementation of an adaptive EKF in the GS IAE design to mitigate the IAE performance degradation due to star motion-induced errors are presented. The simulation study concludes that GS IAE performance degradation during off-nadir pointing for video imaging is insignificant and the AEKF design presents slight improvement for the described representative space program. However, the results show that it does provide significant IAE performance improvement over 40% for spacecraft, such as CubeSats, which exhibit higher motion rates during off-nadir video imaging operations.