A Systematic Error Compensation Method Based on an Optimized Extreme Learning Machine for Star Sensor Image Centroid Estimation

: As an important error in star centroid location estimation, the systematic error greatly restricts the accuracy of the three-axis attitude supplied by a star sensor. In this paper, an analytical study about the behavior of the systematic error in the center of mass (CoM) centroid estimation method under di ﬀ erent Gaussian widths of starlight energy distribution is presented by means of frequency ﬁeld analysis and numerical simulations. Subsequently, an optimized extreme learning machine (ELM) based on the bat algorithm (BA) is adopted to predict the systematic error of the actual star centroid position and then compensate the systematic error from the CoM method. In the BA-ELM model, the input weights matrix and hidden layer biases parameters are encoded as microbat’s locations and optimized by utilizing the strong global search capacity of BA, which signiﬁcantly improves the performance of ELM in terms of prediction accuracy. The simulation result indicates that our method can reduce the systematic error to less than 3.0 × 10 − 7 pixels, and its compensation accuracy is two or three orders of magnitude higher than that of other methods for estimating a star centroid location under a 3 × 3 pixel sampling window.


Introduction
Reliable and accurate attitude determination plays a significant role in aerospace missions. Star sensors provide the most accurate three-axis attitude information when compared with other attitude measurement devices such as the sun sensor, magnetometer, and gyroscope [1,2]. Therefore, star sensors are widely equipped on orbiting satellites and interplanetary spacecraft [3][4][5]. Attitude estimation by star sensors proceeds by comparing star locations in the image taken by star sensors with those in the predefined on-board catalogue. The accuracy of estimating star centroid locations is one critical factor which directly affects the performance of the star sensor [6]. In this paper, we apply the extreme learning machine optimized by the bat algorithm (BA-ELM) to improve the star centroiding accuracy.
To solve the revolution limitation problem of the image plane, defocusing technology, expanding the star image spot to cover neighboring detector pixels, is adopted to obtain centroid precision at the subpixel level [7]. In the past 40 years, several subpixel centroid approaches [8][9][10] have been

Error Analysis of the Center of Mass Method
The CoM method is adopted by most star sensors to estimate the star centroid position with subpixel accuracy since it is easy to implement and exhibits better real-time performance. Therefore, the error of the CoM method is utilized for all discussions in this paper. When the CoM method is applied to digital images captured by the star sensor, the model of the method is described aŝ where (x g ,ŷ g ) is the actual star centroid position, n w is the number of validated discrete pixels in the sampling window, and (x i , y i ) and I i are the geometric center coordinates and the detected signal intensity, respectively, of the ith pixel.
From Equation (1), we can see that the centroiding accuracy of the CoM method is associated with three factors: the size of the sampling window, the coordinates of the validated pixel, and the corresponding detected signal intensity. The use of both the sampling window and the geometric center coordinates is the nature of the CoM method, while the uncertainty in the detected signal intensity is related to random noise caused by dark current noise, radiation noise, readout noise, etc. As a result, the error of the CoM method includes two parts, namely, the systematic error and the random error. The systematic error is derived from discrete approximation and sampling window truncation in the centroiding process, while the random error is introduced as image sensor noise.
Firstly, the random error resulting from uncertainty in the detected I i is analyzed in this section. Since there are various noises in the sensor image, the measured signal intensity I i of the ith pixel includes two components: the starlight intensity E i and the noise intensity σ I . According to the statement in Reference [15], the random error σ rand in the x direction of the CoM method can be given as where Through Equation (2), it can be seen that the σ rand value is determined by the signal-to-noise ratio (SNR), so improving the SNR, such as by optimizing the circuits or choosing low-noise image sensors, is an effective method to reduce the random error [20]. Scholars have proposed some methods to eliminate the random error, and further details concerning these are described in [12,21].
In this paper, we focus on the systematic error in the CoM method. The frequency field method and simulations are utilized to present an explicit analysis of the two types of error sources that bring about the systematic error, namely, the discrete approximation error and the sampling window truncation error.

Discrete Approximation Error
The star image sampling process includes three steps for star sensors, as illustrated in Figure 1. Firstly, the ideal star signal intensity function I(x) projected on the surface of the Charge-coupled Device (CCD) is convoluted with the pixel response function p(x) to generate the continuing pixel signal function f (x). Then, after multiplying the pixel sampling function s(x), f (x) is discretized and the corresponding sampled pixel signal function f s (x) is obtained. Finally, f s (x) is multiplied by the window response function r(x) to give the windowed pixel signal function g(x) utilized in star centroid estimation by the CoM method. The above process can be written as Appl. Sci. 2019, 9,4751 4 of 21 where f (x) = I(x, x 0 ) ⊗ p(x) (5) where I(x, x 0 ) denotes the ideal star signal intensity function with ideal centroid position x 0 in the x direction of the image plane.
Appl. Sci. 2019, 9,  where I(x, ) denotes the ideal star signal intensity function with ideal centroid position in the x direction of the image plane. We assume that the fill factor is 100%, pixel sensitivity is uniform, and each pixel in the image plane has the same photon response. Then, p(x) and r(x) are equal to the rectangle function rect(x), s(x) is equal to a comb function comb(x) with 1⁄ sampling frequency (T is the length of a pixel), and Equation (4) can be written as If we represent all distances in unit T, Equation (6) can be rewritten as In this section, we discuss the systematic error under the case where almost the complete sampled pixel signal (x) is included in the window, so (x) can be considered equivalent to g(x); then According to the properties of Fourier transform, the estimated centroid is represented as follows [15]: where G(s) is the Fourier transform of g(x).
We set a pixel center as the coordinate system origin O; then the function f(x) of the star spot located at d can be expressed as the function (x) shifted by offset d from the origin O, i.e., From Equations (5) and (10), we notice that d = . Thus, the Fourier transform of the function f(x) is where (s) is the Fourier transform of (x). From Equations (8) and (11), G(s) is given as We assume that the fill factor is 100%, pixel sensitivity is uniform, and each pixel in the image plane has the same photon response. Then, p(x) and r(x) are equal to the rectangle function rect(x), s(x) is equal to a comb function comb(x) with 1/T sampling frequency (T is the length of a pixel), and Equation (4) can be written as If we represent all distances in unit T, Equation (6) can be rewritten as In this section, we discuss the systematic error under the case where almost the complete sampled pixel signal f s (x) is included in the window, so f s (x) can be considered equivalent to g(x); then According to the properties of Fourier transform, the estimated centroidx g is represented as follows [15] where G(s) is the Fourier transform of g(x).
We set a pixel center as the coordinate system origin O; then the function f (x) of the star spot located at d can be expressed as the function f e (x) shifted by offset d from the origin O, i.e., From Equations (5) and (10), we notice that d = x 0 . Thus, the Fourier transform of the function where F e (s) is the Fourier transform of f e (x). From Equations (8) and (11), G(s) is given as Substituting Equations (12) and (13) into Equation (9) yieldŝ Therefore, the systematic error is obtained as Equation (15) is the general expression of the systematic error in the CoM method considering only the sampling frequency limitation. Then, we can derive a more explicit relationship by substituting the specific function f (x) of the star spot in the image sensor into the above equation.
The detected star is viewed as a point source for star sensors, and it usually covers a certain pixel region in the actual star images due to the defocus measure. Therefore, a Gaussian function is a reasonable approximation to describe the signal blurring effect [7,22], and the 2-D signal intensity distribution of the star spot can be written as where I 0 is the total energy of the starlight, (x 0 , y 0 ) is the ideal star centroid position, and σ PSF is the Gaussian width. For the 1-D case, Equation (16) can be reduced to Since I(x, 0) is an even function, the function f e (x), which is the convolution result of I(x, 0) and p(x), is also even. Thus, we have From Equation (5), it follows that where F {} denotes the Fourier transform operation. Then Therefore, the theoretical expression of the systematic error with a Gaussian distribution of starlight energy can be written as Actually, when the condition σ PSF > 0.2 is satisfied, exp[−2(πσ PSF ) 2 ] exp[−2(πnσ PSF ) 2 ] (n ≥ 2); then Equation (21) is simplified to Equation (22) indicates that the relationship between σx g and x 0 can be described as a sinusoid function under sampling frequency limitation. Moreover, the amplitude of σx g decreases exponentially as the Gaussian width σ PSF increases. To obtain the relationship in an actual situation, numerical simulations were performed in the following part.
During the experiment, the sampling window was fixed at 3 × 3 pixels, and the ideal point spread function (IPSF) model [23] was utilized to simulate a space discretized digital gray star image spot, that is, the signal intensity detected by the ith pixel was After drawing a group of relationship curves between σx g and x 0 under the different Gaussian widths, the corresponding 3-D numerical simulation results are illustrated in Figure 2.
Equation (22) indicates that the relationship between σ and can be described as a sinusoid function under sampling frequency limitation. Moreover, the amplitude of decreases exponentially as the Gaussian width increases. To obtain the relationship in an actual situation, numerical simulations were performed in the following part.
During the experiment, the sampling window was fixed at 3 × 3 pixels, and the ideal point spread function (IPSF) model [23] was utilized to simulate a space discretized digital gray star image spot, that is, the signal intensity detected by the ith pixel was After drawing a group of relationship curves between and under the different Gaussian widths, the corresponding 3-D numerical simulation results are illustrated in Figure 2. In Figure 2, there is an approximately sinusoid relationship when the Gaussian width is at a small value ( < 0.3 pixels), and the numerical simulation results are also consistent with the expression in Equation (22). The relationship curve tends to be linear when the Gaussian width is relatively large ( > 0.6 pixels), which results from the appearance of truncation error caused by the limitation of the sampling window size; this will be discussed in detail in the next section.

Sampling Window Truncation Error
From Equations (21) and (22), we can see that the systematic error in the CoM method will be reduced by increasing the Gaussian width of the star spot when only considering the sampling frequency limitation. Nevertheless, the received starlight energy becomes weaker for a pixel as the Gaussian width increases, which results in deterioration of the SNR of the star image, resulting in In Figure 2, there is an approximately sinusoid relationship when the Gaussian width σ PSF is at a small value (σ PSF < 0.3 pixels), and the numerical simulation results are also consistent with the expression in Equation (22). The relationship curve tends to be linear when the Gaussian width σ PSF is relatively large (σ PSF > 0.6 pixels), which results from the appearance of truncation error caused by the limitation of the sampling window size; this will be discussed in detail in the next section.

Sampling Window Truncation Error
From Equations (21) and (22), we can see that the systematic error in the CoM method will be reduced by increasing the Gaussian width of the star spot when only considering the sampling frequency limitation. Nevertheless, the received starlight energy becomes weaker for a pixel as the Gaussian width increases, which results in deterioration of the SNR of the star image, resulting in greater random error. Therefore, a sampling window is utilized to limit the number of pixels around the actual star spot. When the sampling window size is relatively small, starlight energy falls out the window, and the signal of the detected star may be truncated asymmetrically; this causes a truncation error. Figure 3 illustrates how window truncation leads to error in the centroid calculation process.
ci. 2019, 9, x FOR PEER REVIEW 7 ow, and the signal of the detected star may be truncated asymmetrically; this causes a trunc . Figure 3 illustrates how window truncation leads to error in the centroid calculation proc n Figure 3a, it can be seen that the discrete star spot signal function e(k) detected by CCD p ies about four pixels in the image plane, while the width of the sampling window is only s (Figure 3b); as a result, the discrete signal function e(k) is truncated asymmetrically whe detected signal is multiplied by the window response function. As shown in Figure 3c l pixel detected signal in e(k) located outside the sampling window r(x) will be exclud lating the star centroid position. The nonsymmetrical truncation causes the estim oid to be smaller than the ideal centroid , that is, a truncation error is introduced int estimated results of the CoM method.
ere, we analyze the numerical simulation results under Gaussian widths of 0.3 and s to present the influence of the truncation error on the relationship curve of and .
igure 4 shows the respective pixel detected signals in the three-pixel sampling window w eal star centroid positions are = 2, 2.25, 2.5 under Gaussian widths = 0.3, 0.671 p an see that the complete pixel detected signal is included in the sampling window with xels during the process of moving the star centroid from the pixel center to the edge (Figur ce the width of the effective pixel detected signal is less than three pixels when σ is value. In this case, the pixel detected signal will not be truncated asymmetrically by ling window; thus, the systematic error is associated with the discrete approximation error lationship curve between and is sinusoid ( Figure 5a). However, for the pixel det l with = 0.671 pixels, the width is relatively lager than three pixels. As a result, the det In Figure 3a, it can be seen that the discrete star spot signal function e(k) detected by CCD pixels occupies about four pixels in the image plane, while the width of the sampling window is only three pixels ( Figure 3b); as a result, the discrete signal function e(k) is truncated asymmetrically when the pixel detected signal is multiplied by the window response function. As shown in Figure 3c, the partial pixel detected signal in e(k) located outside the sampling window r(x) will be excluded in calculating the star centroid position. The nonsymmetrical truncation causes the estimated centroidx g to be smaller than the ideal centroid x 0 , that is, a truncation error is introduced into the final estimated results of the CoM method.
Here, we analyze the numerical simulation results under Gaussian widths σ PSF of 0.3 and 0.671 pixels to present the influence of the truncation error on the relationship curve of σx g and x 0 . Figure 4 shows the respective pixel detected signals in the three-pixel sampling window when the ideal star centroid positions are x 0 = 2, 2.25, 2.5 under Gaussian widths σ PSF = 0.3, 0.671 pixels. We can see that the complete pixel detected signal is included in the sampling window with σ PSF = 0.3 pixels during the process of moving the star centroid from the pixel center to the edge (Figure 4a-c), since the width of the effective pixel detected signal is less than three pixels when σ PSF is at a small value. In this case, the pixel detected signal will not be truncated asymmetrically by the sampling window; thus, the systematic error is associated with the discrete approximation error, and the relationship curve between σx g and x 0 is sinusoid ( Figure 5a). However, for the pixel detected signal with σ PSF = 0.671 pixels, the width is relatively lager than three pixels. As a result, the detected signal is truncated asymmetrically when it locates at a non-central position of the pixel, as shown in Figure 4e,f. Under this condition, the systematic error is mainly determined by window truncation error, the relationship curve becomes linear (Figure 5b), and the maximum error value appears when the star centroid is located at the pixel edge. the relationship curve between and is sinusoid ( Figure 5a). However, for the pixel detected signal with = 0.671 pixels, the width is relatively lager than three pixels. As a result, the detected signal is truncated asymmetrically when it locates at a non-central position of the pixel, as shown in Figure 4e,f. Under this condition, the systematic error is mainly determined by window truncation error, the relationship curve becomes linear (Figure 5b), and the maximum error value appears when the star centroid is located at the pixel edge.

Materials and Methods
The key to compensating systematic error in the CoM method lies in predicting the corresponding systematic error of the actual star centroid position accurately and rapidly.
Researchers have proposed several compensation methods in recent years, such as the BP method [16], analytical compensation (AC) method [17], bivariate polynomial method [18], and LSSVR compensation method [19], and they can reduce the systematic error to some extent, but there are some shortcomings in the compensation model of these methods. For example, the poor performance indices and low learning rate of the BP algorithm, along with how easily it becomes trapped in a local optimum, limit the compensation accuracy of this method; the simplified approximation and iteration estimation in the analytical compensation method lead to a reduction in the prediction accuracy and high on-line computational complexity, which is not suitable for the on-orbit embedded system; the bivariate polynomial compensation template is valid only for some specific Gaussian width cases, and its application range is thus limited; and the problem of scientifically setting the penalty factor and kernel parameter in LSSVR still remains unsettled, meaning that model training is more difficult because of a time-consuming parameter selection process. ELM [24,25] has attracted attention in robot control [26], human face recognition [27], medical diagnosis [28], sales forecasting [29], and protein structure prediction [30] fields, among others, due to its simple training process and excellent generalization ability when compared to other traditional algorithms [31]. However, the randomness of input weights and hidden layer biases become a bottleneck restricting the stability and prediction accuracy of the ELM network. In order to overcome

Materials and Methods
The key to compensating systematic error in the CoM method lies in predicting the corresponding systematic error σx g of the actual star centroid positionx g accurately and rapidly. Researchers have proposed several compensation methods in recent years, such as the BP method [16], analytical compensation (AC) method [17], bivariate polynomial method [18], and LSSVR compensation method [19], and they can reduce the systematic error to some extent, but there are some shortcomings in the compensation model of these methods. For example, the poor performance indices and low learning rate of the BP algorithm, along with how easily it becomes trapped in a local optimum, limit the compensation accuracy of this method; the simplified approximation and iteration estimation in the analytical compensation method lead to a reduction in the prediction accuracy and high on-line computational complexity, which is not suitable for the on-orbit embedded system; the bivariate polynomial compensation template is valid only for some specific Gaussian width cases, and its application range is thus limited; and the problem of scientifically setting the penalty factor and kernel parameter in LSSVR still remains unsettled, meaning that model training is more difficult because of a time-consuming parameter selection process. ELM [24,25] has attracted attention in robot control [26], human face recognition [27], medical diagnosis [28], sales forecasting [29], and protein structure prediction [30] fields, among others, due to its simple training process and excellent generalization ability when compared to other traditional algorithms [31]. However, the randomness of input weights and hidden layer biases become a bottleneck restricting the stability and prediction accuracy of the ELM network. In order to overcome the adverse effects of extreme learning parameters and improve the performance of compensation methods, a systematic error compensation method based on BA-ELM is given herein to eliminate the systematic error of the star centroiding method considering both compensation accuracy and on-line computational load; the BA-ELM model is utilized to predict the systematic error and then improve the accuracy of the CoM method. Figure 6 shows the ELM network structure, which is formed by three layers: the input layer with n nodes, the hidden layer with L nodes, and the output layer with m nodes. The neurons in adjacent layers are fully linked. The input weight values and hidden layer biases in ELM are set randomly, and the activation function in the hidden layer maps the input space onto a new space, then a linear combination is performed on the new space. Thus, the parameter which has to be solved is the weight matrix that links the hidden and output neurons, which can be determined through the least square approach.

Extreme Learning Machine
Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 20 layers are fully linked. The input weight values and hidden layer biases in ELM are set randomly, and the activation function in the hidden layer maps the input space onto a new space, then a linear combination is performed on the new space. Thus, the parameter which has to be solved is the weight matrix that links the hidden and output neurons, which can be determined through the least square approach. According to the theory proposed by Huang [21], the ELM can approximate all given training samples without error, which implies that there are , , and such that Equation (25) can be written in a linear system as where H is the output matrix of the hidden layer and is defined as where is the m-dimensional output vector of the network with input vector is the weight vector between the jth hidden neuron and the output layer, w j = w j1 , w j2 , . . . , w jn T is the weight vector between the jth hidden neuron and the input layer, b j is the bias parameter of the jth hidden neuron, w j · x i denotes the inner product of w j and x i , and g h (x) is the activation function in the hidden layer.
According to the theory proposed by Huang [21], the ELM can approximate all given training samples without error, which implies that there are w j , b j , and β j such that Equation (25) can be written in a linear system as where H is the output matrix of the hidden layer and is defined as where β is the weight matrix connecting the hidden layer and the output layer, and T is the target output matrix. After assigning the coefficient matrices for random parameters, the matrix β can be estimated analytically by finding the minimum norm least squares solution: where H + is the Moore-Penrose generalized inverse of H.

Bat Algorithm
The BA was proposed by Cambridge University scholar Yang [32] and has been approved as a useful method to solve the optimization problem. The principle of BA lies in the echolocation behavior of the microbats. The microbats emit sound pulse, and the sound pulse will be bounced back when it hits objects; then, bats listen for the echoes to estimate the surrounding scenario and search for their prey.
We assume that the microbats pursue the optimal solution in a d-dimensional space, and the velocity and location of the ith bat at step t are v t i and x t i , respectively. Then the updated velocity v t+1 i and location x t+1 i at step t + 1 are given by where x * is the current global best solution in the bat group. F i is the frequency of the sound pulse emitted by the ith microbat, and it can be obtained as shown in Equation (31): where F max , F min are the maximum and the minimum sound frequency of the microbat.
After selecting a solution from the current best solutions, each bat performs a local search and produces a related new solution x new utilizing the random walk according to Equation (32): where A t is the average loudness of the sound pulse of the microbats at step t.
If the new solution is better than the current best solution, it will be accepted as an optimal solution with a certain probability. Meanwhile, the emission rate r t i and loudness A t i of the ith microbat at step t + 1 are updated by where α and γ are constants and 0 < α < 1, γ > 0. Based on the above description of BA, its pseudo code is shown in Algorithm 1. In fact, BA is much superior to other popular metaheuristic approaches such as particle swarm optimization (PSO), genetic algorithm (GA), and Harmony search [32]. Therefore, due to its strong ability in searching for optimal solutions, BA is utilized in engineering and industrial applications. Rank the bats and find the current best x* end while

Optimized Extreme Learning Machine
The random selected parameters in the basic ELM may cause unacceptable prediction results and poor stability of the network [33][34][35]. To promote prediction accuracy, the parameters were coded as the microbat locations, with the corresponding fitness values being the root-mean-square error (RMSE) between the prediction value and actual value. BA was utilized to search for the optimal parameter combination during off-line training. The fundamental steps of BA-ELM are summarized as follows: Step 1: Define the computational parameters relating to BA in the BA-ELM model, such as the maximum iteration number N_max and population size P_num.
Step 2: Define the ELM neutral network structure and the activation function (x).
Step 3: Set k = 1. The population set is randomly generated as the initial locations for the microbat group, in which each location of a bat possesses all random parameters. The position of the ith bat in the kth generation is denoted by ( ) and defined as where ,( , ) and ,( , ) are the biases and input weights, respectively, of the jth hidden node for ( ).
Step 4: Determine the output weights through Equations (36) and (37), and then evaluate the fitness value of the position of ( ) by Equation (38):

Optimized Extreme Learning Machine
The random selected parameters in the basic ELM may cause unacceptable prediction results and poor stability of the network [33][34][35]. To promote prediction accuracy, the parameters were coded as the microbat locations, with the corresponding fitness values being the root-mean-square error (RMSE) between the prediction value and actual value. BA was utilized to search for the optimal parameter combination during off-line training. The fundamental steps of BA-ELM are summarized as follows: Step 1: Define the computational parameters relating to BA in the BA-ELM model, such as the maximum iteration number N_max and population size P_num.
Step 2: Define the ELM neutral network structure and the activation function g h (x).
Step 3: Set k = 1. The population set is randomly generated as the initial locations for the microbat group, in which each location of a bat possesses all random parameters. The position of the ith bat in the kth generation is denoted by θ i (k) and defined as θ i (k) = w 1,(i,k) , . . . , w j,(i,k) , . . . , w L,(i,k) , b 1,(i,k) , . . . , b j,(i,k) , . . . , b L,(i,k) . (35) where b j,(i,k) and w j,(i,k) are the biases and input weights, respectively, of the jth hidden node for θ i (k).
Step 4: Determine the output weights through Equations (36) and (37), and then evaluate the fitness value of the position of θ i (k) by Equation (38): where t s and o s are the target value and output value, respectively, of the sth training sample, H (i,k) is the output matrix for θ i (k), and H + (i,k) is the Moore-Penrose generalized inverse of H (i,k) .
Step 5: Initialize the location θ i and fitness value Fit[θ i ] of the ith bat as θ i (k) and F[θ i (k)], respectively, in which k = 1.
Step 6: Set k = k + 1, i = 1. If the maximum iteration N_max is reached, go to Step 12; otherwise, go to Step 7.
Step 9: Update the location θ i and fitness value Fit[θ i ] of the ith microbat as θ i (k) and Here, rand is a random number in [0, 1] and A i is the loudness of the sound pulse of the ith microbat.
Step 10: Find the current best location among the microbats.
Step 11: Set i = i + 1. If the population size P_num is reached, go to Step 6; otherwise, go to Step 7.
Step 12: Export the optimized ELM neutral network.

BA-ELM for Systematic Error Compensation
The BA-ELM model was utilized as the estimating function in our systematic error compensation method. In practice, the ideal star centroid position x 0 in the captured image cannot be obtained, but we can calculate the actual star centroid positionx g by the CoM method. According to Equation (24), the functional relationship between σx g andx g estimated by BA-ELM can be written as where w * j and b * j are the optimal parameters in the ELM obtained by BA. In practical operations, if we inputx g into the BA-ELM model, it will output the corresponding predicted systematic error δ BA−ELM (x g ), and then we have the compensated positionx c g by Equation (41). Through this way, the systematic error of the CoM method can be eliminated to the maximum extent.

Simulations and Discussion
In this section, experiments were conducted to evaluate the capacity of the proposed BA-ELM compensation method. The experiments included three steps. Firstly, we collected training samples before employing the BA-ELM to estimate the function. Secondly, some parameters which have an impact on estimating function performance were determined. The cross-validation method was utilized to optimize these parameters and obtain the BA-ELM model whose fitting accuracy and prediction accuracy met the requirements. Thirdly, we used the simulated star images to validate the BA-ELM systematic error compensation method in improving the accuracy of the CoM method. Finally, the performance of our method and that of other compensation methods were compared and discussed. All these simulations were written in Matlab R2015a and performed on an Intel Core 3.4 GHz computer.

Preprocessing the Sample Data
To estimate the relationship between the ideal star centroid position x 0 , the actual star centroid positionx g , and the systematic error σx g under the different Gaussian widths through the BA-ELM method, we devised a numerical simulation method to obtain the relationship function among them. In this paper, the sampling window was set to be 3 × 3 pixels, and two situations were simulated. The first case was when σ PSF is smaller than the size of the sampling window, in which the systematic error is only caused by approximation error. The second case was when σ PSF becomes larger than the size of the sampling window, and the systematic error includes both approximation error and truncation error. Thus, σ PSF was set to be 0.3 and 0.671 for these cases, respectively. The same method can be employed to obtain the compensation template and eliminate the systematic error under different Gaussian widths.
In our experiment, the star spot was projected onto the position (5, 7) of the image plane, and the star centroid position in the x direction ranged from 5 to 6. The pixel was subdivided evenly into 1000 parts, so the ideal star centroid position x 0 was 5.0001, 5.0002, 5.0003, ···, 6, that is, the simulation step was 0.001 pixels. In each simulation step, the actual star positionx g and its systematic error σx g defined in Equation (15) were recorded. The 1000 samples were acquired to estimate the relationship function betweenx g and σx g when the Gaussian width was 0.3 or 0.671, as shown in Figure 7a,b, respectively, and they were used as input training samples for the BA-ELM model. 1000 parts, so the ideal star centroid position was 5.0001, 5.0002, 5.0003, ···, 6, that is, the simulation step was 0.001 pixels. In each simulation step, the actual star position and its systematic error defined in Equation (15) were recorded. The 1000 samples were acquired to estimate the relationship function between and when the Gaussian width was 0.3 or 0.671, as shown in Figure 7a,b, respectively, and they were used as input training samples for the BA-ELM model. In Figure 7, the maximum systematic error was 0.054 pixels when was 0.3, while under the condition = 0.671, the maximum systematic error increased to more than 0.1 pixels, which is large enough to degrade the star sensor performance. We proposed the BA-ELM compensation method to decrease the systematic error. The 1000 samples were utilized to train the BA-ELM to estimate the relationship curves, and then the fitting accuracy and prediction accuracy were calculated in order to evaluate the BA-ELM compensation method. In Figure 7, the maximum systematic error was 0.054 pixels when σ PSF was 0.3, while under the condition σ PSF = 0.671, the maximum systematic error increased to more than 0.1 pixels, which is large enough to degrade the star sensor performance. We proposed the BA-ELM compensation method to decrease the systematic error. The 1000 samples were utilized to train the BA-ELM to estimate the relationship curves, and then the fitting accuracy and prediction accuracy were calculated in order to evaluate the BA-ELM compensation method.

The Fitting Accuracy of the BA-ELM
There are some key parameters influencing the BA-ELM performance, such as the hidden neuron number (HN_num), activation function g h (x), population size (P_num), and iteration number (N_max). During the process of optimizing the parameters, the RMSE in N training samples was adopted as the evaluation criterion: where δ BA−ELM (x g ) denotes the prediction output of the BA-ELM with inputx g . In consideration of the accuracy and computation time of the BA-ELM for estimating the function, we set HN_num = 25, P_num = 25, and N_max = 50 in the model. The prediction accuracies of the sig activation function and the sin activation function were compared by calculating their RMSE values, and the RMSE of the sin function was larger than that of the sig function by two orders of magnitude, so the sig function was selected as the activation function g h (x). Figure 8 shows the performance of the BA-ELM in the function estimation.  Figure 9. Here, we defined the fitting error as the difference between the actual systematic error and the predicted systematic error ( ) at the sampling positions.
From Figure 9, it can be seen that when the values were 0.3 and 0.671, the maximum fitting errors were 2.75 × 10 pixels and 2.88 × 10 pixels-both smaller than 3 × 10 pixels. The relationship functions in Figure 8 almost overlap with the fitting curves estimated by the BA-ELM model; this illustrates that the BA-ELM can achieve excellent fitting accuracy under Gaussian widths σ PSF of 0.3 and 0.671. Detailed results on the fitting errors of the BA-ELM model are shown in Figure 9. Here, we defined the fitting error as the difference between the actual systematic error σx g and the predicted systematic error δ BA−ELM (x g ) at the sampling positions. From Figure 9, it can be seen that when the σ PSF values were 0.3 and 0.671, the maximum fitting errors were 2.75 ×10 −7 pixels and 2.88 ×10 −7 pixels-both smaller than 3 × 10 −7 pixels. Gaussian widths of 0.3 and 0.671. Detailed results on the fitting errors of the BA-ELM model are shown in Figure 9. Here, we defined the fitting error as the difference between the actual systematic error and the predicted systematic error ( ) at the sampling positions.
From Figure 9, it can be seen that when the values were 0.3 and 0.671, the maximum fitting errors were 2.75 × 10 pixels and 2.88 × 10 pixels-both smaller than 3 × 10 pixels.

The Prediction Accuracy of the BA-ELM
The high fitting accuracy in the training samples cannot illustrate the performance of the BA-ELM model comprehensively. In practical operations, the BA-ELM is to be utilized to predict the systematic error of the star centroid at random positions. Therefore, the prediction accuracy is a critical parameter to evaluate the BA-ELM model. Firstly, the prediction accuracy is defined in this section.

The Prediction Accuracy of the BA-ELM
The high fitting accuracy in the training samples cannot illustrate the performance of the BA-ELM model comprehensively. In practical operations, the BA-ELM is to be utilized to predict the systematic error of the star centroid at random positions. Therefore, the prediction accuracy is a critical parameter to evaluate the BA-ELM model. Firstly, the prediction accuracy is defined in this section.
Letx g be the actual star centroid position obtained by the CoM method, and let δ BA−ELM (x g ) be the predicted systematic error of the BA-ELM model with inputx g ; we have the compensated star centroidx c g asx c g =x g − δ BA−ELM (x g ).
Then, the prediction error ξ of the BA-ELM model can be expressed by the following equation: where x 0 is the ideal star centroid position. After obtaining a trained BA-ELM model with optimal parameters, its prediction accuracy was calculated. In the simulations, 600 positions on the image plane were selected randomly, and then we made a star spot centroid traverse these positions sequentially. The coordinates in the x direction of these selected positions ranged from 1 to 512. Figure 10 shows the respective results of the prediction accuracy in the 600 random positions under σ PSF = 0.3 and σ PSF = 0.671. The red lines were formed by connecting all ideal coordinates (x 0 , σx g ) successively, while the blue lines were formed by connecting the compensated coordinates (x c g , δ BA−ELM x g ). From Figure 10a2,b2, which are enlarged pictures relative to Figure 10a1,b1, we can see that the blue lines nearly the overlap red lines; this means that each compensated star positionx c g is extremely close to its corresponding ideal position x 0 in all test positions. The prediction error is illustrated in Figure 11.
( , ) successively, while the blue lines were formed by connecting the compensated coordinates ( , ( )). From Figure 10a2,b2, which are enlarged pictures relative to Figure 10a1,b1, we can see that the blue lines nearly the overlap red lines; this means that each compensated star position is extremely close to its corresponding ideal position in all test positions. The prediction error is illustrated in Figure 11.  As shown in Figure 11, the maximum prediction errors of the BA-ELM were 2.5 × 10 −7 pixels and 3.0 × 10 −7 pixels when the values were 0.3 and 0.671, respectively. The test results indicate that the proposed compensation method based on the BA-ELM model can achieve excellent prediction accuracy under different Gaussian widths.

The Compensation Performance in Star Images
After evaluating the compensation performance on single star spots, we tested the BA-ELM compensation method on a simulated star image. For the simulations, the configuration of the star sensor is shown in Table 1, the field of view (FOV) direction was selected randomly, and the As shown in Figure 11, the maximum prediction errors of the BA-ELM were 2.5 × 10 −7 pixels and 3.0 × 10 −7 pixels when the σ PSF values were 0.3 and 0.671, respectively. The test results indicate that the proposed compensation method based on the BA-ELM model can achieve excellent prediction accuracy under different Gaussian widths.

The Compensation Performance in Star Images
After evaluating the compensation performance on single star spots, we tested the BA-ELM compensation method on a simulated star image. For the simulations, the configuration of the star sensor is shown in Table 1, the field of view (FOV) direction was selected randomly, and the Smithsonian Astrophysical Observatory (SAO) J2000 star catalog was used to synthesize star images. To completely verify the systematic error compensation method under different Gaussian widths, two simulated star images were generated, as shown in Figures 12 and 13, whose Gaussian widths σ PSF were 0.3 and 0.671 and FOV directions were (269.85, −55.64) and (315.84, −26.360), respectively.  As shown in Figure 11, the maximum prediction errors of the BA-ELM were 2.5 × 10 −7 pixels an × 10 −7 pixels when the values were 0.3 and 0.671, respectively. The test results indicate th e proposed compensation method based on the BA-ELM model can achieve excellent predictio curacy under different Gaussian widths.

. The Compensation Performance in Star Images
After evaluating the compensation performance on single star spots, we tested the BA-ELM mpensation method on a simulated star image. For the simulations, the configuration of the sta nsor is shown in Table 1, the field of view (FOV) direction was selected randomly, and th ithsonian Astrophysical Observatory (SAO) J2000 star catalog was used to synthesize star image completely verify the systematic error compensation method under different Gaussian width o simulated star images were generated, as shown in Figures 12 and 13, whose Gaussian width were 0.  Principle point (256 pixels, 256 pixels) Pixel size 15 mm × 15 mm Magnitude threshold 6.0 Mv In Figure 12, the 10 brightest stars were selected, and their ideal position in the x direction (IX sition in the x direction before compensation (XBC), position in the x direction after compensatio AC), and errors in the x direction before compensation (EXBC) and after compensation (EXAC) a ted in Table 2. In Figure 13, the 10 brightest stars were selected, and their ideal position in the x direction (IX sition in the x direction before compensation (XBC), position in the x direction after compensatio AC), and errors in the x direction before compensation (EXBC) and after compensation (EXAC) a ted in Table 3. Through comparing the EXBC and EXAC values in Tables 2 and 3, it was noticed that th oposed BA-ELM compensation method can drastically reduce the systematic error and achiev gh accuracy in estimating star centroid locations. In Figure 12, the 10 brightest stars were selected, and their ideal position in the x direction (IX), position in the x direction before compensation (XBC), position in the x direction after compensation (XAC), and errors in the x direction before compensation (EXBC) and after compensation (EXAC) are listed in Table 2. In Figure 13, the 10 brightest stars were selected, and their ideal position in the x direction (IX), position in the x direction before compensation (XBC), position in the x direction after compensation (XAC), and errors in the x direction before compensation (EXBC) and after compensation (EXAC) are listed in Table 3.  Tables 2 and 3, it was noticed that the proposed BA-ELM compensation method can drastically reduce the systematic error and achieve high accuracy in estimating star centroid locations.

Comparison and Discussion
To verify the superiority of the BA-ELM compensation method, comparison studies were conducted with other methods, namely, the BP method, AC method, and LSSVR compensation method. In addition, to further illustrate the better prediction accuracy and generalization abilities of the BA-ELM, the performance of the ELM method in eliminating the systematic error was also evaluated in the comparison experiment. The results on the compensation accuracy of these methods are shown in Table 4. From Table 4, we can see that the prediction error of our method was much smaller than those of the other compensation methods-the BP method (2.0 × 10 −3 pixels), the AC method (5.0 ×10 −4 pixels), the LSSVR compensation method (6.0 ×10 −5 pixels), and the ELM method (4.0 ×10 −5 pixels). Our method needs longer training time because of the parameter optimization, but the compensation model training is an off-line process.

Conclusions
In this paper, a comprehensive study on the systematic error in the CoM method considering both the detector sampling frequency limitation and the sampling window size limitation was presented. The BA-ELM was applied to predict the corresponding systematic error of the actual star centroid position and then improve the star centroiding accuracy of the CoM method. The input weights matrix and hidden layer biases were optimized with BA in an off-line training process, which promoted the prediction accuracy and stability of the basic ELM model. Several simulations were implemented to test the BA-ELM compensation method, and the results indicate that the proposed method can effectively eliminate the systematic error under different Gaussian widths of starlight energy distribution, and it achieves higher accuracy in estimating the star centroid location after compensation than do other methods. Future work will focus on optimizing the search ability of BA to obtain better parameters in ELM and further validate the compensation method with a real star sensor. In summary: (1) Through a frequency field and numerical simulations approach, the characteristics and causes of two types of systematic error in the CoM method, namely, discrete approximation error and truncation error, were illustrated, and the relationships between the systematic error and actual star centroid position under different Gaussian widths were obtained. This helped us to design the compensation method. (2) BA was adopted to improve the prediction accuracy and stability of the ELM. The random parameters in the BA-ELM model were optimized during the off-line training process. Through the verification of the results, the proposed BA-ELM compensation method was shown to effectively eliminate the systematic error and enhance the accuracy of the CoM method in estimating star centroid locations. (3) Comparing the simulation results of the BP neural network method, analysis compensation method, LSSVR compensation method, and basic ELM showed that the BA-ELM compensation method performed better in terms of prediction accuracy than did the other methods during the systematic error compensation process.