On-Orbit Geometric Distortion Correction on Star Images through 2D Legendre Neural Network

Abstract

The star tracker is a prerequisite device to realize high-precision attitude determination for a spacecraft. However, due to the errors in optical lens machining, optical path assembly, and temperature alternation, optical instruments suffer from some amount of optical geometric distortion, resulting in declining star tracker accuracy. The on-orbit distortion correction of star images is indispensable for precise performance. In this paper, a novel single-layer 2D Legendre neural network (2DLNN) to automatically correct the geometric distortion of the star tracker is proposed. An offline training method grounded on batch star images and an online training algorithm based on sequential star images are designed, respectively. The 2DLNN realizes the ground-based and on-orbit online correction of optical geometric distortion for the star tracker. The 2DLNN features self-learning, lifelong learning, and good adaptability. The single-layer neural network is simple, quick convergence, which is suitable for on-orbit implementation. The simulations demonstrate that the average distortion error can be reduced to less than 0.04 px after ground-based training. In the earth-orientation mode of the LEO satellite, the on-orbit sequential training algorithm can converge in 2500 star images under 1 frame/s. The proposed 2DLNN can achieve high-precision correction at the sub-pixel level, effectively improving the star tracker’s attitude determination accuracy.

1. Introduction

Star trackers are the most accurate attitude measurement instruments for spacecraft. They are widely applied in multiple fields and missions, such as self-contained navigation, environmental monitoring, global communication, national security, and deep-space exploration [1,2,3,4,5]. In recent years, with the release and implementation of satellite constellations [6,7], the developing trend of star trackers is miniaturization and low cost. However, the requirements for star trackers are surging rather than falling, which brings greater challenges to star trackers.

Regardless of the the craftsmanship of the optical system of the star tracker, all-optical instruments, whether ground-based or space-based, suffer from some amount of optical geometric distortion. Optical lenses with a large field of view (FOV) will generate a large amount of nonlinear distortion, and the distortion can reach more than two pixels at the image plane edge. These distortions have multiple origins that induce different behavior. They may be static, resulting from unavoidable errors in the shape and placement of the optics or from the imperfect fabrication of the detectors. In general, laboratory calibration and testing are conducted to obtain true values of all parameters, and these errors contribute little to the distortion of the star images [8,9,10,11,12,13].

Distortion may be dynamic; many factors, such as intense vibration during launching, component aging, temperature effects [14,15,16] and a variable environment during operation, alter the systematic errors and change the camera parameters on orbit [17,18], resulting in the star spot imaging being uneven in energy distribution and asymmetrical shape. Star trackers use a detector with a rolling shutter, usually active pixel sensors (APS). Under dynamic conditions, it is well-known that the detected star image is deformed. The centroid of the starspot is shifted [19,20,21]. Consequently, the accuracy of extracting the centroid is directly impacted. To guarantee and ensure the attitude accuracy and the safety of the spacecraft, it is essential to establish and optimize the distorted models to complete the correction of distorted star images.

The on-orbit calibration or distortion correction for spacecraft attitude star trackers was first proposed by Samaan [22], in which the method makes use of residuals between measured inter star cosine angles and the inter star cosine angles between the corresponding cataloged stars to learn the calibration corrections on orbit. The approach makes use of the truth that inter star angles are an invariant of rotational transformations. Therefore the correction process is naturally divided into two major parts, calibration of the principal point offset and focal length and then calibration of the focal plane distortions.

Based on Samaan’s method, various on-orbit calibration algorithms have been developed. The attitude-dependent methods [23,24,25,26], gyroscopes, infrared horizon sensors, sun sensors and other attitude devices are required to determine the initial attitude information first. However, the error of the star tracker frame attitude affects the correction inevitably. Thus, the attitude-independent approach has been proposed. Many attitude-independent algorithms estimate the star tracker optical parameters by directly utilizing the star coordinate in the star image plane and the corresponding star vector in the inertial coordinate frame. The most commonly used of these methods are the least-squares method and the Kalman Filter algorithm [17,18,22,23,27,28,29,30,31,32,33].

In these methods, the simple second-order radial distortion model is considered. To correct the starspot position error caused by the rolling shutter exposure, the authors of [20] used a first-order approximation of the real motion of the stars across the detector. The proposed strategy set the centroids to the same time instant. In Ref [21], the authors corrected the starspot position at an asynchronous time to the same time based on time-domain constraints. The focal length errors and the principal point will affect the accuracy of star identification, and the above methods are based on the assumption that the principal point and focal length change greatly. In this paper, the errors caused by the change of principal point are considered systematic errors, and the random error is not considered.

Methods based on polynomial fitting have been developed to correct other complex distortions. Ordinary polynomials fit the mapping relationship between distorted stars and reference stars. Wei [34] used the improved genetic algorithm to obtain the optimal coefficient, while Liu [35] utilized machine learning to solve the polynomial fitting. However, according to the author of the paper [36], it has been proven that the solution of ordinary polynomials is the system of linear equations, known as the normal equations.

As the order or the numbers of polynomials increases, however, solving a large system of equations becomes unstable and inaccurate. Compared with ground-based correction, the biggest problem of on-orbit correction is that the types of distortion are often unknown. Some distortions are stable, such as radial distortion and tangential distortion caused by the change of the optical system parameters. There are also some image distortions caused by other factors, such as perspective distortion.

We aim to design an algorithm that can automatically correct the star coordinates of the on-orbit star tracker and is insensitive to the type of image distortion. In this paper, a novel on-orbit geometric distortion correction method that combines 2D Legendre polynomials and a neural network is proposed. To avoid the problem of instability and inaccuracy, 2D normalized Legendre polynomials that are orthogonal polynomials are used to replace the ordinary polynomials. A 2D Legendre neural network (2DLNN) is constructed to obtain the optimal coefficient of the Legendre polynomial.

The neural network has strong adaptability and fast convergence speed, and can also correct any distorted star images. For the different distortions that may occur in the lifelong period of on-orbit star trackers, we can adaptively learn the distortion model by adjusting the number of neurons. Simulations and experiments are conducted to show effectiveness and superiority. The proposed correction method has the potential to become a universal geometric distortion correction framework for star trackers.

Our main contribution is to propose learning-based methods to correct a wide range of geometric distortions of star trackers blindly. More specifically, we propose:

• A single-layer 2D Legendre neural network, which learns the coefficients of 2D Legendre polynomial.
• An offline training method grounded on the batch star images, which realizes the initial correction of the star tracker.
• An on-orbit online learning algorithm based on the sequential star images, which realizes lifelong autonomous distortion correction of the star tracker.

The paper is organized as follows. In Section 2, measurement models related to the star tracker are described. In Section 3, the architectures of 2DLNN are described. Two training methods of offline and online learning are given in detail. In Section 4, the simulation results are shown, and in Section 5, our conclusions and final remarks are reported.

2. Measurement Model of the Star Tracker

2.1. Ideal Model

The star tracker is designed to measure the stars’ direction vector. The ideal measurement model of the star tracker can be considered as a pinhole imaging system [1], as shown in Figure 1.

The i-$\mathrm{th}$ measurement star vector ${\omega }_i\in R^{3\times 1}$ in imaging frame $(X,Y,Z)$ can be expressed by Equation (1),

$${\omega }_i=\frac{1}{\sqrt{x_i^2+y_i^2+f^2}}\left[\begin{array}{c}-x_i\\ -y_i\\ f\end{array}\right]$$

(1)

where $(x_i,y_i)$ is the i-$\mathrm{th}$ star centroid, and f is the focal length of the optical system.

The relationship between ${\omega }_i$ and $v_i$ can be expressed as follows,

$${\omega }_i=A{v}_i$$

(2)

where $v_i$ is the direction vector of the navigation star, and A is the attitude matrix of the star tracker. When the number of navigation stars exceeds two, the attitude matrix can be estimated by the QUEST [37,38] algorithm.

2.2. Typical Star Image Distortion Model

The existence of errors and noise in the optical imaging system is inevitable. Noise in the CMOS APS sensor is a combination of dark current, single point noise, Gaussian readout noise, Poisson noise, fixed pattern noise, etc. [39]. The aberrations in the optical system usually consist of radial distortion [40], tangential distortion [41], thin prism distortion [42] and image transformation [43]. Due to the various influencing factors in space, the optical system parameters may change, and various types of geometric distortion will inevitably occur in star images. Hence, the pinhole model cannot fit with real star imaging exactly.

The relationship between the distorted star spot coordinates $(x_i^{\prime },y_i^{\prime })$ and the real star spot coordinates $(x_i,y_i)$ can be expressed by Equation (3),

$$\begin{array}{cc}\hfill x_i^{\prime }& =x_i+\Delta x_i\hfill \\ \hfill y_i^{\prime }& =y_i+\Delta y_i\hfill \end{array}$$

(3)

where $(\Delta x_i,\Delta y_i)$ is the distorted magnitude of star coordinates caused by the change in the parameters of the star tracker. In this paper, five distortions and their causes are briefly described. The manifold distorted star images, the distorted flow of pixels, and the changing of starspot coordinates are shown in Figure 2. Several common distortions of the on-orbit star tracker are shown below.

• Radial Distortion
  The vertical magnification of the optical system in different fields of view is different. The light is more crooked from the center to the edge of the lens, which results in radial distortion and is the main source of distortion. The expression of radial distortion is shown in Equation (4),

  $$\begin{array}{cc}\hfill \Delta x_i^r& =\sum _{k=1}^n{q}_k{x}_i{r}^{2k}\hfill \\ \hfill \Delta y_i^r& =\sum _{k=1}^n{q}_k{y}_i{r}^{2k}\hfill \end{array}$$

  (4)

  where n is the order of Taylor expansion, $r=\sqrt{{(x_i-x_0)}^2+{(y_i-y_0)}^2}$ is the distance of the distortion star and the center of the image plane, and $q_i$ is the i-$\mathrm{th}$ radial distortion coefficient. Barrel distortion typically has a negative term for $q_i$, whereas pincushion distortion will have a positive value.
• Tangential Distortion
  The coaxial error of the optical components in the optical system leads to tangential distortion, and the expression of this type of distortion is shown in Equation (5),

  $$\begin{array}{cc}\hfill \Delta x_i^t& =p_1(3x_i^2+y_i^2)+2p_2{x}_i{y}_i+O\left[{(x_i,y_i)}^4\right]\hfill \\ \hfill \Delta y_i^t& =2p_1{x}_i{y}_i+p_2(x_i^2+3y_i^2)+O\left[{(x_i,y_i)}^4\right]\hfill \end{array}$$

  (5)

  where $p_1$ and $p_2$ are the tangential distortion coefficients.
• Thin Prism Distortion
  The thin prism distortion is caused by the tilt of the optical component, which is equivalent to inserting a thin prism into a non tilted optical component. The thin prism distortion can be expressed by Equation (6),

  $$\begin{array}{cc}\hfill \Delta x_i^{tp}& =s_1(x_i^2+y_i^2)+O\left[{(x_i,y_i)}^4\right]\hfill \\ \hfill \Delta y_i^{tp}& =s_2(x_i^2+y_i^2)+O\left[{(x_i,y_i)}^4\right]\hfill \end{array}$$

  (6)

  where $s_1$ and $s_2$ are the thin prism distortion coefficients.

Generally, typical distortions of the star trackers are often a superposition of the three distortions described above. Image transformation is also considered in this paper to test the extendability and adaptability of the proposed method. Different image transformations can be expressed as follows,

• Shear Distortion

  $$\left[\begin{array}{c}{x}_i^{s_h}\\ y_i^{s_h}\\ 1\end{array}\right]=\left[\begin{array}{ccc}tan{\varphi }_x& 0& 0\\ 0& tan{\varphi }_y& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\\ 1\end{array}\right]$$

  (7)

  where $tan{\varphi }_x$ and $tan{\varphi }_y$ are the coefficients of shear distortion.
• Perspective Distortion
  Perspective distortion can be regarded as an imaging plane transformed into a new imaging plane [8], and the relationship between two imaging planes is transformed by a $3\times 3$ transformation matrix.

  $$\left[\begin{array}{c}{x}_i^{\prime }\\ y_i^{\prime }\\ w\end{array}\right]=\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& 1\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\\ 1\end{array}\right]$$

  (8)

  $$\begin{array}{c}\hfill x_i^{p_e}=\frac{x_i^{\prime }}{w}=\frac{h_{11}{x}_i+h_{12}{y}_i+h_{13}}{h_{31}{x}_i+h_{32}{y}_i+1}\\ \hfill y_i^{p_e}=\frac{y_i^{\prime }}{w}=\frac{h_{21}{x}_i+h_{22}{y}_i+h_{23}}{h_{31}{x}_i+h_{32}{y}_i+1}\end{array}$$

  (9)

  where $h_{11}, h_{12}, h_{13}, h_{21}, h_{22}, h_{23}, h_{31}, h_{32}$ are the elements of the transformation matrix.

2.3. Legendre Polynomial Distorted Modeling

The coordinate of measured star i is noted as $(X_i^{mea},Y_i^{mea})$, the corresponding reference coordinate of star i is $(X_i^{ref},Y_i^{ref})$, by using the distortion function as $D_j$,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{X}_i^{mea}& =D_{j,X}(X_i^{ref},Y_i^{ref})\\ Y_i^{mea}& =D_{j,Y}(X_i^{ref},Y_i^{ref})\end{array}\right.\end{array}$$

(10)

where $D_{j,X}$ and $D_{j,Y}$ are the j-$\mathrm{th}$ distortion model of X direction and Y direction, respectively.

Here, we can use different functions to describe the distortion of the star images. Most commonly, 2D polynomials up to third or fifth order, depending on the authors, are used as image distortion models. For example, Wei [34] used a third-order polynomial fitting to correct distortion. The improved genetic algorithm was adopted to obtain the optimal solutions. Liu [35] used machine learning to solve the coefficients of ordinary polynomials. However, polynomials are unstable and inaccurate for solving a large system of equations.

To avoid this, orthogonal polynomials should be used. The authors of [44] compared the efficiency of different models: 2D Cartesian polynomials, bivariate B-spline, and 2D Legendre polynomial. They concluded that the 2D Legendre polynomial basis provides faster convergence and lower residuals. Simultaneously, the authors of [45] defined an orthogonal polynomial, which is well adapted to describe the distortion on square images. Finally, the ascending polynomial degree organization of the basis is convenient for characterizing the distortion using a limited number of modes.

Each Legendre polynomial $P_n\left(x\right)$ is an n-$\mathrm{th}$ degree polynomial and may be expressed in one dimension by Equation (11),

$$P_n\left(x\right)=\frac{1}{2^n}\sum _{k=0}^n{\left(\begin{array}{c}n\\ k\end{array}\right)}^2{(x-1)}^{n-k}{(x+1)}^k,x\in [-1,1]$$

(11)

As each mode contains the same amount of distortion in both directions (x and y), the polynomials need to be normalized. The final basis of 2D Legendre polynomials can be expressed by Equation (12),

$$b_{Nmodes}=\left\{\begin{array}{cc}{L}_0\left(x\right)\times L_0\left(y\right)& \\ L_1\left(x\right)\times L_0\left(y\right)& \\ L_0\left(x\right)\times L_1\left(y\right)& \\ L_1\left(x\right)\times L_1\left(y\right)& \\ L_2\left(x\right)\times L_0\left(y\right)& \\ L_0\left(x\right)\times L_2\left(y\right)& \\ \dots \end{array}\right\}=\left\{\begin{array}{cc}\frac{1}{2}& \\ \frac{\sqrt{3}}{2}x& \\ \frac{\sqrt{3}}{2}y& \\ \frac{\sqrt{3}}{2}xy& \\ \frac{\sqrt{5}}{4}(x^2-1)& \\ \frac{\sqrt{5}}{4}(y^2-1)& \\ \dots \end{array}\right\},(x,y)\in {[-1,1]}^2$$

(12)

where $N_{modes}$ is the number of modes considered, it determines the dimension of the vector $b_{Nmodes}$ and can be calculated by Equation (13),

$$N_{modes}=\frac{(n+1)(n+2)}{2}$$

(13)

where n is the order of the Legendre polynomial, and $L_n$ is the n-$\mathrm{th}$ degree normalized Legendre polynomial,

$$L_n=\frac{P_n}{\parallel P_n\parallel }$$

(14)

$${\parallel P_n\parallel }^2=\frac{2}{2n+1}.$$

(15)

The distortion modes are referred to using their index mode (m) defined as Equation (16),

$$m(k,l)=\frac{l(l+1)}{2}+k,0\le k\le l\in \mathbb{N},0\le m\le N_{modes}$$

(16)

Hereafter, the mode of every Legendre polynomial can be expressed by Equation (17),

$$b_m(x,y)=L_{l-k}\left(x\right)\times L_k\left(y\right)$$

(17)

finally, the coordinates of the i-$\mathrm{th}$ star in the j-$\mathrm{th}$ frame, on which a distortion function $D_{j,X}$ and $D_{j,Y}$ can be written as Equation (18),

$$\left\{\begin{array}{cc}{X}_i^{mea}=D_{j,X}(X_i^{ref},Y_i^{ref})=X_i^{ref}+{\displaystyle \sum _{l=0}^n}{\displaystyle \sum _{k=0}^l}{\omega }_{i,m,X}{L}_{l-k}\left(X_i^{ref}\right)L_k\left(Y_i^{ref}\right)& \\ Y_i^{mea}=D_{j,X}(X_i^{ref},Y_i^{ref})=Y_i^{ref}+{\displaystyle \sum _{l=0}^n}{\displaystyle \sum _{k=0}^l}{\omega }_{i,m,Y}{L}_{l-k}\left(X_i^{ref}\right)L_k\left(Y_i^{ref}\right)\end{array}\right.$$

(18)

where ${\omega }_{i,m,X}$ and ${\omega }_{i,m,Y}$ are referred to as the distortion coefficients of the X and Y directions, respectively.

The statistical properties of the noise that we should consider to build the distortion model must be consistent with the noise present in the measurement of star positions. We assume here a noise that follows a Gaussian white noise. The final distortion model is

$$\left\{\begin{array}{cc}{X}_i^{mea}=D_{j,X}(X_i^{ref},Y_i^{ref})=X_i^{ref}+{\displaystyle \sum _{l=0}^n}{\displaystyle \sum _{k=0}^l}{\omega }_{i,m,X}{b}_m(X_i^{ref},Y_i^{ref})+n_{i,X}& \\ Y_i^{mea}=D_{j,X}(X_i^{ref},Y_i^{ref})=Y_i^{ref}+{\displaystyle \sum _{l=0}^n}{\displaystyle \sum _{k=0}^l}{\omega }_{i,m,Y}{b}_m(X_i^{ref},Y_i^{ref})+n_{i,Y}\end{array}\right.$$

(19)

where $(n_{i,X},n_{i,Y})$ is the observation error caused by imaging noise and other factors, which is regarded as the independent Gaussian noise.

3. Legendre Neural Network

The fully-connected neural network implicitly learns the relationship between the camera and the global coordinate system [33], which is proven effective; however, the structure is complicated. Traditional neural networks have drawbacks, such as difficulty determining initial weights and the number of layers and neurons. In general, the convergence time is also long when trained by the backpropagation algorithm.

In the previous section, we defined the direct model. To solve the distortion problem, we now aim to inverse this model by minimizing a criterion. In this paper, a single-layer neural network based on the Legendre polynomial is developed to avoid the problems of traditional neural networks and is named 2DLNN. The single-layer neural network is simple enough to reduce the calculation time to implement the distortion correction algorithm in orbit.

No complicated activation function is included, making it possible to have a rapid response to obtain the results from the input data. The hyper-parameters of the network do not require much tuning. In addition, the authors of [46,47] had proven that a function could be approximated by an orthogonal function set (such as the Legendre function), and the coefficients are unique and bounded. Orthogonal polynomials functions have recursive properties to determine their expansion terms. The flow chart of the geometric distortion correction method based on 2DLNN is shown in Figure 3.

3.1. Network Architectures

In this paper, the optimal coefficients of the Legendre polynomial are solved to fit the distortion of star images. The distortion modeling of the star images can be described as Equation (19), ${\left({\omega }_{i,m,X/Y}\right)}_{0\le m\le N_{modes}}$ are the parameters that need to be estimated. The input of the 2DLNN is the coordinates $(X_{mea},Y_{mea})$ of measuring stars, the output is the offset $(\widehat{\Delta x},\widehat{\Delta y})$ of coordinates between the measured stars centroid and the centroid of reference stars. Consequently, the 2D Legendre neural network is constructed to learn the mapping of measuring stars’ coordinates and coordinates offset. The estimation offset is the 2D vector field that specifies where pixels in the star image with distortion should move to find the corrected coordinates.

According to Equation (19), a single-layer Legendre neural network as shown in Figure 4 can be built.

The 2DLNN is trained to minimize the error between the estimated displacement of coordinates and the position offset. The offset is the distance between measuring stars and reference stars. Equation (19) can be abbreviated as Equation (20),

$$F_X=X^{mea}-X^{ref}=\sum _{l=0}^n\sum _{k=0}^l{\omega }_m{b}_m(x,y)+n_x=W_X^T{b}_m(x,y)+n_x$$

(20)

The 2DLNN can be expressed by Equation (21), taking the distortion in the x direction as an example, and the distortion in the y direction can be obtained in the same way.

$${\widehat{F}}_{X/Y}=\sum _{l=0}^n\sum _{k=0}^l{\widehat{\omega }}_m{b}_m(x,y)={\widehat{W}}_{X/Y}^T{b}_m(x,y)$$

(21)

where ${\widehat{F}}_{X/Y}$ is the estimation error of starspot position, ${\widehat{W}}_{X/Y}$ is the estimated coefficients of 2D Legendre polynomials.

According to Figure 4, the input of the network is the measuring starspot coordinates $(x,y)$, and the output is the coordinate offset $(\widehat{\Delta x},\widehat{\Delta y})$. The processing elements are duplicated from the expansion terms of the normalized Legendre polynomials. The n-$\mathrm{th}$ processing element is identical to the n-$\mathrm{th}$ expansion term of the 2D normalized Legendre polynomials.

The weights concatenating the mid-layer and the out-layer are needed to update. The single-layer neural network has two sets of adjustable weights which are $({\omega }_{m,X},{\omega }_{m,Y})$, and these weights will be trained to approach the desired weights or coefficients ${\omega }_m$ in the Equation (19) by the training algorithm. Since each output has unique and independent weights, the two sets of weights can be trained separately.

3.2. Ground-Based Batch Learning

After having built the 2D normalized Legendre neural network, it is necessary to find a method to adjust the weights ${\widehat{W}}_{X/Y}$ to approach the desired weights W defined in Equation (20). The training method of 2DLNN in the ground-based calibration is introduced first. The star images can be captured from random areas of the sky, and the distribution of stars is easier to cover the whole field of view of the star tracker, which is conducive to learning the global distortion model of the star image. Therefore, it is a point-to-point training method in the ground-based experiment, which learns the mapping from the whole star image datasets.

The estimation error E can be expressed by Equation (22),

$$E=F_X-{\widehat{F}}_X={(W_X-{\widehat{W}}_X)}^T{b}_m(x,y)+n_x$$

(22)

The loss $J_p$ is defined as Equation (23),

$$J_p=\frac{1}{2}{E}^2$$

(23)

The gradient descent algorithm is adopted for the weights update. According to the chain rule, the extreme value of $J_p$ can be obtained by calculating the first-order partial differential of $J_p$, as shown in Equation (24),

$$\nabla J=\frac{\partial J\left(\widehat{W}\right)}{\partial \widehat{W}}=\frac{\partial \left(\frac{1}{2}{E}^2\right)}{\partial E}\frac{\partial E}{\partial \widehat{W}}=E{b}_m(x,y)$$

(24)

The Adam algorithm [48] is adopted to update the weights of 2DLNN. Adam has the attractive benefits of being straightforward to implement with low memory requirements, and the hyper-parameters have intuitive interpretation and typically require little tuning. It is suitable for on-orbit applications. Its convergence has been proven theoretically. Therefore, the weights updating of 2DLNN can be expressed by Equation (25),

$${\widehat{W}}_t=\tau \times {\widehat{W}}_{t-1}-\eta \times \frac{{\widehat{m}}_t}{(\sqrt{{\widehat{v}}_t}+\epsilon )}$$

(25)

where $\tau =0.999$ is the forgetting factor, ${\widehat{m}}_t$ and ${\widehat{v}}_t$ are the parameters of Adam. $\eta =0.0001$ is a positive gain, referred to as the learning rate or step size. Equation (25) is the training algorithm used to train the weights. The author of Adam proved that $\widehat{W}$ would decrease gradually and reach a stable status. The Adam is described more fully in the Appendix A Equation (A1); however, a brief summary is given here.

3.3. On-Orbit Online Learning

Due to the limitation of satellite storage, the star tracker cannot store a large number of star images and starspot coordinates. It is unrealistic to use the point-to-point training method for the on-orbit star tracker. We first calculate the starspot coverage rate at a specific time for the star tracker orbiting the earth. Then, the online training method is proposed.

3.3.1. Coverage of Star Spot

It is impossible to verify the algorithm directly on orbit, we simulate the state of the star tracker orbiting the earth. To ensure that the proposed method can learn the distortion model of the whole focal plane, stars in a period need to cover most of the star image. The coverage rate of starspot is expressed by a ratio of the area $Are{a}_{sp}^t$ covered by stars in FOV to the area of the entire focal plane in time t, as shown in Equation (26),

$$C_{star}=\frac{Are{a}_{sp}^t}{Are{a}_{fp}}$$

(26)

where $Are{a}_{sp}^t$ is the area of the contour profile (not convex hull) of all stars within the time interval t, $Are{a}_{fp}$ is the area of the focal plane, and the unit of them is $\left(p{x}^2\right)$.

The star images are photoed along the orbit. The Monte Carlo algorithm is used to calculate the average star coverage rate in the time t, which is shown in

Table 1. Average coverage rate of star spot for different times.

Time (s)	1	50	500	2000
Coverage rate	$42.267\%$	$71.232\%$	$98.262\%$	$99.309\%$

. The sampling time is 1 s, and the orbit is the Starlink (altitude: 550 km, inclination: 53°), and the right ascension and declination of the orbit are generated from the System Tool Kit (STK) software.

The superposition of starspots in FOV at a different time is shown in Figure 5. They are collected along the Starlink orbit, and the sampling start is determined by the Monte Carlo method. Scatters with the same colors represent stars from the same star image. The coverage area is the envelope of these superimposed stars, as shown in Figure 5a the pale blue background. The coverage rate of stars increases with the increase of star images. Although the coverage rate of 500 frames is high, there are obvious holes in the envelope, and the superposition of 2000 frames covers the entire FOV.

The stars can be regarded as passing through the FOV of the star tracker at a certain angular velocity. In the engineering application, the probability of starspot being correctly identified increases when the stars are in a position with less distortion, and the correctly identified stars can be tracked in real-time. When the tracked star moves to a position with large distortion, it can continue to be used for distortion correction. In Figure 5, it can be seen that the stars enter the FOV from multiple directions so that 2DLNN can learn and correct the global distortion.

3.3.2. On-Orbit Learning Algorithm

The geometric distortion correction of the on-orbit star tracker is implemented frame-to-frame along the orbit. The loss should be determined by the root mean square error (RMSE) of all correctly identified stars in each frame of the star image. The loss is expressed by Equation (27),

$$J_f=\sqrt{\frac{1}{N}\sum _{i=0}^N{\left(E^i\right)}^2}=\sqrt{\frac{1}{N}\sum _{i=0}^N{(F^i-\widehat{F^i})}^2}=\sqrt{\frac{1}{N}\sum _{i=0}^N{\left\{{(W-\widehat{W})}^T{b}_m(x^i,y^i)\right\}}^2}$$

(27)

where N is the number of correctly identified stars in a single star image, i represents the centroid coordinates of i-$\mathrm{th}$ measuring star, and N is different in each star image.

The gradient of online training can be also obtained according to the chain rule, which is expressed by Equation (28),

$$\nabla J_f=\frac{\partial J_f\left(\widehat{W}\right)}{\partial \widehat{W}}=\frac{\partial \left(\sqrt{\frac{1}{N}{\displaystyle \sum _{i=0}^N}{\left(E^i\right)}^2}\right)}{\partial E^i}\frac{\partial E^i}{\partial \widehat{W}}=E^i{b}_m(x^i,y^i){\left(N\sum _{i=0}^N{\left(E^i\right)}^2\right)}^{\frac{1}{2}}$$

(28)

The optimal coefficient $W^{*}$ of the 2D Legendre polynomial is obtained by calculating the minimum gradient of the loss.

$$W^{*}=\underset{W}{argmin}\left(J_{p/f}\left(W\right)\right)$$

(29)

Similarly, the Adam algorithm (Equations (25) and (A1)) is also adopted to update the weight of the network to obtain the optimal solution $W^{*}$.

3.4. Optimum Order of the Legendre Polynomial

The orders of the Legendre polynomial will have an impact on the results of 2DLNN. The order determines the number of neurons in the neural network, which directly impacts the ability of distortion correction. If the order is too small, the star tracker distortion model cannot completely be expressed by the Legendre polynomials. If there are too many neurons, the network structure is large, and the dataset would not be enough to fit the network. Therefore, it is necessary to find a suitable order for the better ability of 2DLNN.

We analyze the influence of polynomial order on training results. Taking radial distortion as an example, 2DLNN with orders increasing from 1 to 10 are trained, respectively. The corresponding network is applied to the testing dataset. The training results of different orders are shown in Figure 6.

Figure 6a is the result of the training dataset, and Figure 6b is the performance of the corresponding network on the testing dataset. When the order is small, the rate of convergence of 2DLNN is slow. The low-order 2DLNN is not convergent on the testing data, indicating that the low-order 2DLNN is underfitting. When the order is 5 or larger, the error remains the minimum on training data (black curve in Figure 6a) and testing data (red circle in Figure 6b). In practical application, the common orders in the optical distortion model are 4 to 6. Considering the time and memory complexity, the optimum order of the 2DLNN should be 5.

3.5. Pseudocode

We summarize our methods for geometric distortion correction of star images in Appendix B Algorithm A1.

4. Experiment and Analysis

Several sets of experiments were designed to verify the validity of the proposed method. They are ground-based and on-orbit experiments. The star images can be captured from random or continuous sky regions in the real night sky observation experiment. The stars in the star images from random sky areas are uniformly distributed, and the numbers are enough for batch learning. The star images obtained from the continuous sky areas were generated along the orbit, which is suitable for frame-to-frame online learning. In the earth-orientation mode of the LEO satellite, the star images generated along the orbit were used to analyze the implementation of the on-orbit training method.

The measurement model was adopted from Figure 1 in Section 2.1. The parameters of the star tracker utilized in the experiments are shown in

Table 2. Star tracker parameters.

Parameters	Value	Unit
Detector area	2048 × 2048	pixels
Pixel size	7.45	$\mathsf{\mu }$m
Field-of-View (FOV)	20	degree
Focal length	43.20	mm
Max magnitude	6.0	Mv

. The optical system distortion (radial distortion, tangential distortion and thin prism distortion) and image transformation (shear transformation and perspective transformation) were tested, respectively. The parameters of various distortions are shown in

Table 3. Parameters of optical system distortion.

Distortion	Value
pincushion	$q_1=\phantom{-}4.638\times {10}^{-6}$,    $q_2=\phantom{-}4.988\times {10}^{-9}$
barrel	$q_1=-4.638\times {10}^{-6}$,    $q_2=-4.988\times {10}^{-9}$
tangential	$p_1=\phantom{-}1.210\times {10}^{-6}$,    $p_2=-2.460\times {10}^{-6}$
thin prism	$s_1=\phantom{-}1.290\times {10}^{-6}$,    $s_2=\phantom{-}1.670\times {10}^{-6}$

and

Table 4. Parameters of star image transformation.

Distortion	Value
shear	$s_1=0.36051$,         $s_2=0.00789$
perspective	$h_{11}=0.999$,    $h_{12}=0.000625$,    $h_{13}=0.0005$
	$h_{21}=0$,          $h_{22}=0.999248$,    $h_{23}=0.0010$
	$h_{31}=0$,          $h_{32}=0.001250$,    $h_{33}=1$

, the distortion model is adopted in Section 2.2.

All experiments and algorithms were simulated in PyCharm CE (V2021.3) with Python 3.8 & Pytorch on a PC with a 2.30 GHz Intel Xeon Gold 6139 CPU of 32 GB RAM, NVIDIA Quadro P4000 GPU of 8 GB RAM and Windows 10 operational system.

4.1. Training Dataset of Star Images

Since real star images are rarely downloaded from on-orbit star trackers, it is unrealistic to use real star images as datasets. In this paper, the star image datasets were collected and generated based on the real on-orbit star tracker parameters. At the same time, the star images datasets were mixed with star images obtained in other ways. The parameters setting of the star tracker in the datasets were the same. The star image datasets are mainly composed as follows.

• Inspired by the literature [18], star images can be captured randomly from Starry Night, a professional and powerful astronomical software, and the digital platform is utilized to obtain the distorted data.
• According to the parameters of the real on-orbit star tracker, the method [49] for generating simulation star images was adopted. Random background noise and random position noise were added to these star images.
• Star images were taken from the real night sky observation experiments, using the star trackers’ parameters in Table 2.

Star images were distorted according to the real distortion model of the star tracker in orbit, which is mainly caused by the change in optical system parameters. The star spot extraction algorithm [39,50,51,52] and the star identification algorithm [53,54,55] were adopted to obtain navigation stars. Then, the transfer matrix from celestial sphere frame to star tracker frame can be calculated using QUEST [37,38] or TRIAS [56] method. The reference starspot coordinates can be obtained through the coordinate transforming equations. The offset distance between the measured starspots and the reference starspots was used as the training label of 2DLNN. Meanwhile, the stars’ centroid coordinates were used as the training data.

The datasets were divided into three parts: ground-based star images from random sky areas (named after GRD) and continuous sky areas (named after GCD). The on-orbit star images along the Starlink orbit (named after OD).

4.2. Simulation of Ground-Based Correction

There is no need to select specific starspots; all available stars in the star images can participate in network model fitting. In the ground-based experiment, we combine the available star coordinates into a complete dataset (namely GRD and GCD). The centroid coordinates of the stars obtained from the distorted star images were used as the input of 2DLNN. Meanwhile, the coordinate offsets were used as the training label. The error will be trained to the minimum. The dataset is randomly divided into training and testing datasets according to the ratio of 8:2.

The hyper-parameters of the 2DLNN were: $batch\_size=128$, $learning\_rate=1\times {10}^{-4}$. The order of the Legendre polynomial is $n=5$; hence, the number of neurons of the Legendre neural network is $N_{modes}=21$. Six distortion models were simulated, and the network structure was the same. The only difference is the input and output data and the weight of the trained model (i.e., the coefficient of Legendre polynomial).

The distortion correction is first implemented for the dataset GRD. The stars in the FOV of the star tracker are evenly distributed, which is conducive to learning the global distortion of the star image. The training process is shown in Figure 7.

Similarly, batch learning was implemented to train the distortion model with the dataset GCD. The training process is shown in Figure 8.

The initial weights of 2DLNN were the same for the different distortions. It can be seen from Figure 7 that the initial estimation errors are different due to different distortion models. The initial weights can be randomly generated or obtained through a swarm intelligence optimization approach, such as the sparrow search algorithm (SAA [57]). After training, the estimation error of the model will converge to less than $3\times {10}^{-5}$ (normalized pixel).

Figure 8 considers the possibility of on-orbit implementation. We will discuss on-orbit distortion correction in Section 4.3. The ground-based experiment demonstrates that the method can fit different distortions, and the star position error can reduce to a small value after correction. The weight, trained by the ground-based experiment, can be used as the initial values of on-orbit correction. However, with different environments between ground and on-orbit, it is not easy to obtain on-orbit parameters beforehand in the laboratory. Therefore, on-orbit distortion correction is online continuation learning.

4.3. Simulation of On-Orbit Correction

Based on the conclusion in Section 3.3.1, the dataset OD is adopted to train the on-orbit distortion model. Unlike the ground-based experiment, the on-orbit geometric distortion correction is frame-to-frame training. The on-orbit training method of 2DLNN is Equation (27), and the training process is shown in Figure 9.

From Figure 9, considering the resource needed for correction, about 2000 to 2200 images are needed in the rough correction. Less than 2500 star images are needed totally to achieve satisfactory correction results. Also, taking radial distortion as an example, the average running time required to train a single star image only using CPU is $0.6$ ms, which is suitable for on-orbit execution.

The results of the geometric distortion correction for a single star spot are shown in Figure 10, the star spot is intercepted from the radial distortion star image. For the original distortion, the distance between the observed star and the reference star is $2.619\times {10}^{-3}\times 1024\approx$ 2.683 pixels. After correction, the star position error is $4.751\times {10}^{-5}\times 1024\approx$ 0.048 pixels, and the relative error of star position can be reduced by (2.683 − 0.048)/2.683 × $100\%\approx 98.18\%$. It shows that the proposed geometric distortion correction algorithm is effective and accurate.

The predicted distortion flow of the image plane is shown in Figure 11. The first line is the distortion flow of the focal plane before correction, and the second line is the distortion flow learned by the proposed method. The average coordinate error of starspot is shown in

Table 5. Average position error of different models.

	B	P	${\mathit{T}}_{\mathit{A}}$	${\mathit{T}}_{\mathit{P}}$	${\mathit{S}}_{\mathit{H}}$	${\mathit{P}}_{\mathit{E}}$
x	0.02131	0.03172	0.04677	0.05288	0.03489	0.03171
y	0.03502	0.03051	0.04931	0.02662	0.03882	0.04632

. At the edge of the star image plane, the distortion caused by the optical system is about 2 pixels, while the maximum distortion of image transformation is about 3 pixels. After correction using the proposed algorithm, the star position error can be less than $0.04$ pixels. It indicates that the proposed method successfully predicts and corrects the star position error. In addition, for other unknown distortion models, such as distortion composed of a set of basal errors or other irregular distortions, the 2D Legendre neural network can also learn the distortion mapping.

4.4. Comparison with Previous Techniques

Piecewise linear mapping algorithm based on Delaunay triangulation [58], ordinary polynomial mapping method based on improved genetic algorithm [34] were applied to the geometric distortion correction of star images. The previous techniques were implemented on the radial distortion of the star tracker optical system. We only compared the distortion caused by the change of the star tracker optical system parameters, that is, radial distortion (B and P), tangential distortion ($T_A$), and thin prism distortion ($T_P$).

The imaging data and the star spot coordinates in the comparison experiments were the same. The average position error and execution time of three methods after correction for different distortion models are shown in

Table 6. Average position error of different methods and average time consumed.

Algorithm	Direction	Average Position Error (Pixels)				Average Time (s) 1
		$\mathit{B}$	$\mathit{P}$	${\mathit{T}}_{\mathit{A}}$	${\mathit{T}}_{\mathit{P}}$
PWL	x	0.85242	0.92072	0.13455	0.13696	0.0203
	y	0.83611	0.85956	0.25855	0.17416	0.0167
OP-IGA	x	0.67175	0.70221	0.10717	0.10839	0.0606
	y	0.65959	0.67365	0.20652	0.13932	0.0482
Ours	x	0.02131	0.03172	0.04677	0.05288	0.0061
	y	0.03502	0.03051	0.04931	0.02662	0.0059

¹ The average time consumed of four distortion on CPU for a single star image.

, which are calculated in the direction of x and y, respectively. It can be seen that the proposed method has significant advantages in correcting the position of starspots and star images. The average execution time is quick, which is suitable for onboard implementation.

5. Conclusions

The star tracker is a prerequisite device to realize high-precision attitude determination for a spacecraft. However, star tracker accuracy is often influenced by optical lens machining, optical path assembly, and temperature alternation errors. Optical lens parameter changes and image distortion will have an impact on the accuracy of the star tracker. To improve the accuracy of the star tracker and the availability of star images, a novel general framework based on 2DLNN for the geometric distortion correction of on-orbit star images was proposed. An offline training method grounded on batch star images and an online training algorithm based on sequential star images were designed, respectively.

2DLNN is a single-layer neural network with fast convergence, self-learning, lifelong learning, and good adaptability, and it is suitable for on-orbit implementation. The simulations demonstrate that the average position error of distortion can be reduced to less than 0.04 px after correction. In the earth-orientation mode of the LEO satellite, the on-orbit sequential training algorithm can converge in 2500 star images under the condition of $1\mathrm{frame}/\mathrm{s}$. The proposed method has the potential to become a general framework for the geometric distortion correction for star images. The method proposed in this paper can be extended to the image correction of satellite-carried optical sensors and real night observation experiments.