Blurred Star Image Processing for Star Sensors under Dynamic Conditions

The precision of star point location is significant to identify the star map and to acquire the aircraft attitude for star sensors. Under dynamic conditions, star images are not only corrupted by various noises, but also blurred due to the angular rate of the star sensor. According to different angular rates under dynamic conditions, a novel method is proposed in this article, which includes a denoising method based on adaptive wavelet threshold and a restoration method based on the large angular rate. The adaptive threshold is adopted for denoising the star image when the angular rate is in the dynamic range. Then, the mathematical model of motion blur is deduced so as to restore the blurred star map due to large angular rate. Simulation results validate the effectiveness of the proposed method, which is suitable for blurred star image processing and practical for attitude determination of satellites under dynamic conditions.


Introduction
As high accuracy and high reliability devices, star sensors play an important role in attitude determination for satellites in celestial navigation system (CNS). The main steps for attitude determination include star point location, star identification and attitude tracking [1]. Based on the captured star images, stars can be located and identified. Whether the attitude determination of satellite OPEN ACCESS is successful or not, the pattern recognition is very important for star images in the field of view (FOV) [2]. It indicates that only the available and recognizable star images can ensure star sensor can give an accurate satellite attitude [3], so it is critical to improve the accuracy of star point location.
In the past many years, some algorithms have been developed to extract star centroids utilizing initial star images. Reference [4] shows a new sub-pixel interpolation technique to process image centroids. Reference [5] gives a method of enhancement of the centroiding algorithm for star tracker measure refinement. An analytical and experimental study for autonomous star sensing, including the star centroid process, is presented in [6]. However, these studies are generally used under static conditions or balanced processes. For many agile maneuver satellites, star sensors work under dynamic conditions as a result of the rotation of the star sensor along with the satellite. Therefore, various noises in star field caused by dynamic factors may affect the quality of imaging. Moreover, due to the large angular rate, the star point moves on the focal plane during the exposure time, which may lead to two changes for the star point: the position shifting on the focal plane and the limited starlight energy dispersing into more pixels. As a result, the SNR (Signal to Noise Ratio) of blurred star images may decrease and the measured star centroid would be inaccurate.
Dynamic conditions stress the need for a very accurate and robust processing method for blurred star images. At present, many investigations tend to concentrate on the exploration and analysis of star location under dynamic condition [3]. For example, reference [7] shows the simulation analysis of dynamic working performance for star trackers. Reference [8] provides an analysis of the star image centroid accuracy of an APS star sensor under rotation conditions. Blind iterative restoration of images with spatially-varying blur is the research topic in reference [9]. However, most of them are limited to the useless locating capability when the angular rate is larger than 2°/s. In [3], the method is effective to estimate attitude, but it has the contrary effect when the angular rate is low.
The main theme of this paper is to overcome the difficulties arising from dynamic imaging blur of star sensors, including denoising and restoration of blurred star images by estimating the angular rate. On the one hand, if the angular rate is in the dynamic range of the star sensor, a proposed adaptive wavelet threshold is used for denoising according to the characteristics of the blurred star image, which can guve the accurate centriod within sub-pixels. On the other hand, if the angular rate is larger than the dynamic range, the restoration algorithm based on the angular rate is used to process the -tail‖ star image. As will be seen later in this paper, the adaptive method outperforms other denoising methods in terms of Power Signal-to-Noise Ratio (PSNR) and visual qualities, and the large variation of the angular rate has little effect on the star centroid determination based on the restoration method. This paper is divided into six sections. Following this Introduction, the imaging theory of star sensors is outlined in Section 2, as well as the characteristic of blurred star images under dynamic conditions. Then the method of denoising based on adaptive wavelet threshold is described in detail in Section 3. The restoration method is developed in Section 4 by analyzing the Point Spread Function (PSF) of motion blurred star images, where a Wiener filter with an optimal window is used to overcome the edge error. In Section 5, simulation results are shown to demonstrate the proposed methods. At the end, conclusions are drawn in Section 6.

Coordiante Frames
In developing a set of equations to be mechanized by a celestial navigation system and star sensor or in studying the behavior of a given system, it is necessary to introduce several sets of orthogonal coordinates: Inertial coordinate system (X L -Y L -Z L ) has its origin at the center of the Earth and is non-rotating with respect to the fixed stars. Its x-axis is in the equatorial plane and the z-axis is normal to that plane; and the y-axis complements the right-handed system. Star sensor coordinate system (X s -Y s -Z s ) has its origin at the center of mass of the star sensor. Its x-axis points along longitudinal axis of the star sensor; the z-axis is perpendicular to the longitudinal plane of symmetry and is along the boresight of the star sensor; and the y-axis completes a right-handed system. Focal plane coordinate system (X p -Y p ) has its origin at the center of the focal plane. Its x-axis points along longitudinal axis of the focal plane; the y-axis is perpendicular to the longitudinal plane. Figure 1 illustrates the general large FOV star sensor for attitude determination. After capturing stars in the real sky and imaging by the star sensor, the attitude information is completed by an autonomous procedure (including image pretreatment, star centroiding, star map matching, attitude determination, etc.). The schematic of imaging is also shown in Figure 1, where f is the focal length of lens and α is the angle of FOV. According to the coordinates of star points in the focal plane coordinate system X p -Y p , it is easy to obtain the coordinates matrix S of star points in the star sensor coordinate system X s -Y s -Z s . Combined with stars coordinates G in the inertial coordinate system X L -Y L -Z L , the attitude A can be determined as the form of 3 × 3 direct cosine matrix by:

Characteristics of Star Image under Dynamic Conditions
Under static conditions, the distribution of star points is generally represented as a two-dimensional Gaussian with a 3 × 3 or 5 × 5 dispersion circle by defocusing technology [10], so that the accuracy of star centroid can be kept within a sub-pixel level. However, under dynamic conditions, the original star image is perturbed and blurred by various additive noises, which mainly include photon response uniform noise, photon shot noise, dark current noise, readout noise, etc. [11].
At present, the dynamic range of a large FOV star sensor is about 3-5°/s [8]. Suppose the angular rate is w. If w ranges by 1-1.5°/s under dynamic condition, the rotation of the star sensor has little effect on the star images. However, a star sensor may lose tracking using the technique under static conditions due to the various noises caused by the dynamic conditions. On the other hand, if w is larger than the dynamic range, the star point constantly shifts in the focal plane and appears to trail badly during exposure time, which may affect the star centroid accuracy and even result in the failure of attitude determination.
Based on the foregoing discussions, denoising and deblurring are two crucial parts for the pretreatment of blurred star images. Performance parameters of the star sensor used in this paper are shown in Table 1. The CMOS image sensor chip is the STAR 1000 from the Cypress Semiconductor Co. [12].

Denoising Modeling Based on Wavelet Transform
Supposing the size of the clear image f(i, j) is N × N, a common model of the corrupted image g(i, j) is mathematically defined as: (2) where 0  i, j  N-1, and n(i, j) is additive random noise and independent of f(i, j). The goal is to remove n(i, j) and estimate f(i, j) which minimizes the Mean Squared Error (MSE) [13].
In general, the important information of f(i, j) is mostly distributed as a smooth signal at low frequency, while n(i, j) is distributed at high frequency. Based on this, a two-dimensional (2-D) discrete wavelet transform (DWT) can be implemented to transform g(i, j) into the wavelet domain.
Then, wavelet coefficients denoting different scales and orientations can be obtained with the use of the Mallet algorithm [14]. Figure 2 shows the subbands of the orthogonal DWT of three levels. LL 3 is an approximation subband (or the resolution residual) which contains the low frequency portion of g(i, j).The subbands HL k , LH k , HH k (k = 1,2,3) respectively denote the details of vertical, horizontal and diagonal orientations, where k is the scale and size of a subband at scale k is N/2 i × N/2 i . There is no space here to go into detail on this method, and for this level of detail the reader can refer to [13,15] for more information. It is important to point out the small coefficients mostly due to noise and large coefficients due to important signals. Accordingly, denoising can be accomplished by thresholding these coefficients.

Threshold Selection
Thresholding is simple because it operates on one wavelet coefficient at a time. The method of using an adaptive threshold to implement denoising described by Lakhwinder Kaur et al. [16] appears more suitable, in which threshold choice is: where  y is the standard deviation of each subband, and  is the scale parameter for each scale computed by: (4) where J is the largest scale, and L k is the length of subband at the scale of k (k = 1,2…J). The noise variance  2 is obtained by: Studies in [17] indicate that the square error relating to HH 1 of g(i, j) almost equals the noise variance  2 . On the other hand, the larger the decomposition level is, the smaller the weight of noise in the coefficient variance will be. For this reason, it is more convenient to complete denoising for the star images than general images because the contrast between star signal on a black background and noise is more distinct, even when the star point is blurred under dynamic condition. We pay attention to T N in scale J, where has the approximation of star points. This proposed method processes each coefficient in scale J using a different threshold. It can be executed mainly by the following steps: 1. Apply an M × M local window to compute σ 2 lJ , which denotes the coefficients variance of window l in scale J. M is determined by the square root of the number of pixels occupied by the star point, and generally is not more than seven. 2. Compute noise variance  2 according to Equation (5). 3. Obtain the threshold by: (6) where Th l is the threshold in window l of scale J.

Star Image Denoising
Based on the foregoing analysis, the proposed method for star image denoising is summarized as follows: 1. Execute decomposition of the initial blurred star image using a wavelet transform at level K. 2. Compute the noise variance  2 according to Equation (5). 3. Compute the scale parameter β of level K using Equation (4). 4. Use a 4 × 4 square window in LLK to obtain Thl by Equation (6). 5. Process coefficients in scale K using the following threshold function: (7) which keeps the coefficients information if it is larger than threshold; otherwise, it is set to zero. 6. Invert the multiscale decomposition to reconstruct the denoised star image.

Blurred Star Image Restoration
The current angular rate w of satellite can be obtained by the attitude update. As mentioned in Section 2, if w is larger than the dynamic range, the attitude cannot be correctly computed because of the -trailed‖ image. This section mainly focuses on the restoration of motion-blurred star image as a result of large w.

Mathematical Model of Blurred Star Image
Due to the motion of the star sensor during exposure time, what a star sensor captures is a motionblurred image g(x, y). Suppose, a clear image f(x, y) moves on the focal plane, its displacement components of direction x and y can be respectively termed as x(t) and y(t), where t is the movement time during exposure time T. Then, the expression of g(x, y) can be obtained from Equation (8) (9), which is similar as a previous reference [18]: T g x y f x x t y y t dt     (8) If the satellite rotates clockwise about the boresight Z p with an angular rate w z during exposure time T, it seems that the star point displaces anticlockwise on focal plane, which results in the trail effect. The motion-blurred procedure can be illustrated in Figure 3. where X p -Y p is the coordinate system of focal plane, X p ′-Y p ′ is the corresponding coordinate of X p -Y p after rotation, and θ = w z t, where t is the rotation time in T. The star point P shifts to P′ along the circular arc s with radius r in the focal plane. It is reasonable to assume that the rotation axis and the motion parameters are constants, and l is approximate to s. As a result l can be expressed as: (10) where , and (x p , y p ) is the coordinates of P in coordinate system X p -Y p .
Suppose P moves along track l with angle γ to the horizontal axis X p by velocity v, where v = w z r and γ can be obtained by θ, x p and y p . Combined with Equation (8), the PSF of motion blurred star image can be obtained by:

Blurred Star Image Restoration
To accomplish restoration of the original image, the traditional method is to employ Wiener Filtering in the frequency domain [15]. The Wiener filter is intended to be an optimal filter in the sense that it delivers the best estimate of the object in a least-squares sense for additive Gaussian noise. However, the noise n(x, y) is typically unknown in practice and the classical Wiener filter is problematic [19]. Therefore, we use the modified Wiener filter which is given by: (13) where H * (u, v) denotes the complex conjugate of H(u, v) and a can be considered as an adjustable empirical parameter chosen to balance sharpness against noise.
In order to overcome the edge error, a major factor affects the quality in Wiener filter restoration, the optimal window method is used for star image [20]. Then the steps of restoration based on Wiener filtering are detailed as follows: 1. Introduce h(x,y) according to analysis in Section 4.1. 2. Apply the optimal window w(x, y) as a weight factor to g(x, y), then execute the Discrete Fourier Transform (DFT) of g(x, y) and h(x, y). 3. Use the Wiener filter for deconvolution filtering in the frequency domain, and obtain the estimate of F(u, v) by Equation (13). 4. Compute the Inverse DFT (IDFT) of (x,y) to generate f(x, y) by: (14)

Results and Analysis
In order to verify the proposed method when a star sensor works under dynamic conditions, simulations and experiments are implemented to accomplish denoising and restoration according to blurred star images caused by different w. Comparison of PSNR and the star centroid are also analyzed to estimate the effect of algorithm in this section.

Denoising of Blurred Star Image
Based on the performance of the star sensor shown in Table 1, the SkyMap star map simulation software is used to generate the original star image, as shown in Figure 4. The boresight direction is set as (150°, 15°) and the 14,581 stars brighter than 6.95 m are selected in Tycho2n star catalog.   To assess the performance of the denoising method proposed in this paper, it is compared with several common denoising techniques like BayesShrink [13], SureShrink [21] and Lowpass filter. The fixed threshold Th is used first to segment the background and the star object. Based on Th, different denoising methods are employed to estimate the original clear star image. Figure 5 shows the noisy image and resulting images at σ = 90 and w = 0.6°/s. We can see that the image processed by the proposed method outperforms the others in terms of visual quality. Then, the PSNRs from various methods are compared in Table 2, and the data are collected from an average of four runs. The AdaptThr method, namely, is the proposed adaptive thresholding method.
Results in Table 2 show that the lower w is, the better AdaptThr performs than other methods, especially when σ is large. AdaptThr approximately has the same poor effect of denoising along with the increase of w. Actually, in dynamic condition with high w, star image is not only perturbed by various noises, but also is blurred by the motion of star sensor. This also means that by only using the proposed denoising method under dynamic conditions with high w, one cannot obtain the star centroid accurately, and one also needs to restore the motion-blurred image. In order to further verify the proposed denoising algorithm, a real star image is adopted in this section. Figure 6 shows the original star image obtained by a star sensor and its gray distribution, from which we can see that the background value in the star image is large. What's more, there is a big ‗singularity spot' which is larger and lighter than other star points. After discarding the singularity spot, a clear star image can be obtained as shown in Figure 7. It can be seen that the dim star object is extracted perfectly from the background noise. This confirms the notable effect of the proposed method which can adapt to the complex dynamic conditions.   Table 2 shows that if w is larger than the dynamic range, star image needs to be restored by deblurring rather than by denoising directly. Figure 8(a) is a real star image slected from the original images obtained by the CMOS star sensor in this paper. Based on the supposition w = 10°/s, the blurred star image can be generated according to the degradation model, as shown in Figure 8

Restoration of Blurred Star Image
Gray distributions of the same star point in two different images of Figure 8 are respectively shown in Figure 9, which show that due to the motion blur, the star point smears out intensely, as well as the gray value decreases. In this section, we implement star centroiding [2] to assess the performance of the proposed deblurring method. The comparison results are shown in Table 3, where the angular rate is 10°/s.  Table 3 shows that the extraction errors of δx and δy are mostly larger than a pixel for each star centroid without deblurring. This is because the SNR of star image decreases as a result of the star points smearing significantly. Moreover, six star points fail to be extracted due to the low gray value of dim blurred star points (as shown in Figure 9(b)), which may affect the star recognition and attitude before(w=5) before(w=10) before(w=20) after(w=5) after(w=10) after(w=20) determination. However, after restoration in advance, the star centroid can be obtained accurately for in that the extraction errors of δx and δy are within subpixel range, as well as the dim star points with low gray value are extracted. As can be seen in Figure 10, the extraction error of δx + δy error is larger than three pixels. With the restoration method, all lost star points can be extracted as well as the extraction error of δx + δy is restricted within one pixel. This is because the proposed deblurring method can keep the accuracy of δx and δy within subpixel levels, even when w is larger than the dynamic range, and the variation of angular rate w has little effect on the star centroid.

Conclusions
This article researches how to process blurred star images according to different angular rates of star sensors under dynamic conditions. A new denoising method based on adaptive wavelet threshold is proposed, as well as a restoration method according to large angular rate out of the dynamic range. Experiments on different types of star images have been conducted with the proposed algorithm. The PSNRs of images with different types of angular velocity show the proposed denoising method, in comparison with the normal denoising methods, has good performance, namely, better than PSNRs of other methods under the same conditions when the angular velocity is in the dynamic range, and also in terms of visual quality. Star centroiding against blurred star images have been analyzed to assess the effectiveness of restoration. It is confirmed that the restoration maintains the extraction error within subpixel levels, and the variation of angular velocity has little effect on the accuracy of star centroid, which shows that the proposed method is both effective and feasible. Experimental results show that the processing method according to angular velocity in before/after using the restoration method with different angular velocity are analyzed, and star points which can be extracted in each method are also shown. Without restoration, the larger the angular velocity is, the more star points cannot be extracted while the extraction under dynamic conditions reported in this paper could keep star sensors stable within a certain range and meet the requirements of attitude determination, which needs uninterrupted output data and attitude accuracy of arcsecond level.