Smearing model and restoration of star image under conditions of variable angular velocity and long exposure time

Ting Sun; Fei Xing; Zheng You; Xiaochu Wang; Bin Li

doi:10.1364/OE.22.006009

1. Introduction

With the development of earth-observation, deep-space exploration and celestial navigation, the requirements for attitude measurement are rapidly increasing. The star tracker is supposed to have the highest accuracy under stable conditions among different types of attitude measurement devices [1,2]. Dynamic performance is becoming its major limitation for further applications. In the dynamic case, the shape of the star spot on the image sensor will be elongated over several or scores of pixels at high angular velocity, the energy will disperse and the SNR of the star spot will decrease, resulting in difficult extraction, and inaccuracy or even failed star spot position value. To detect faint stars and to reduce the size of the star tracker, the long exposure time method is usually employed. However, this will make the situation worse. Several solutions can be employed to improve the dynamic performance of the star tracker. These solutions include adopting a new type of image sensor with higher sensitivity, and enhancing [3] or recovering the blurred star spot energy with image processing algorithm. Among these solutions, star image restoration is the most convenient and efficient, for which dynamic Point Spread Function (PSF) estimation and restoration algorithm are crucial. Shen [4] and Xing [5] described predictive mathematical models of the star spot under dynamic conditions. Samaan [6] discussed a predictive centroiding method with the image smear effect. Wu [7] and Zhang [8] discussed the degradation models and corresponding PSFs, and they respectively adopted the constrained least square filter and Wiener filter for star image restoration. Whereas, degradation model under motion conditions of variable angular velocity is still urgently needed to be discussed, which can avoid bringing errors in the position accuracy of the restored star spot.

A smearing model estimation and restoration method of the star spot under conditions of variable angular velocity is proposed in this paper. Simulated dynamic star images and laboratory experiments are performed to validate this approach. The simulation and experiment results demonstrate that through presented method, the star spot can be recovered from a long asymmetrical trail to the shape of Gaussian distribution under the conditions of variable angular velocity and long exposure time. The energy of the star spot can be centralized, resulting in high SNR and high position accuracy. These factors are especially important to the subsequent processing and the whole performance of the star tracker.

2. Position model of star spot

2.1 Star tracker measurement model

The star tracker is a high-accuracy attitude measurement device, which considers stars as the measuring targets. It obtains the direction vector from the celestial inertial coordinate system by detecting different stars on the celestial sphere. After many years of astronomical observations, star direction vectors on the celestial sphere are very accurate. As shown in Fig. 1, navigation star S_i with direction vector v_i in the celestial coordinate system can be detected by the star tracker, whereas the vector of its image can be expressed as w_i in the star tracker coordinate system.

Fig. 1 Ideal imaging model of the star tracker.

Download Full Size | PDF

The position of the principal point of the star tracker on the image plane is defined as (x₀, y₀). The position of the image point of navigation star $S_{i}$ on the image plane is (x_i, y_i). The focal length of the star tracker is L_f. Vector w_i can be expressed using Eq. (1) [9]:

w_{i} = \frac{1}{\sqrt{{(x_{i} - x_{0})}^{2} + {(y_{i} - y_{0})}^{2} + L_{f}^{2}}} [\begin{matrix} (x_{i} - x_{0}) \\ (y_{i} - y_{0}) \\ - L_{f} \end{matrix}] .

Equation (2) elucidates the relationship between w_i and v_i under ideal conditions, where A is the attitude matrix of the star tracker.

w_{i} = A v_{i} .

When the number of navigation stars is more than two, the attitude matrix can be solved using the QUEST algorithm [10]. Thus, the optimal attitude matrix in the inertial space of the star tracker can be calculated.

2.2 Position Analysis and motion model of star spot

This section discusses the star spot smearing model when the star tracker is with motion of three-axis angular velocity. As shown in Fig. 2, the coordinate system OXYZ represents the coordinate of the star tracker; O denotes the principal point position of the image plane; OZ is along the optical principal axis of the lens of the star tracker; OX and OY are along the two orthogonal directions of the image plane; OO_C has the same length as the focal length; and $ω_{x}^{t}$ , $ω_{y}^{t}$ and $ω_{z}^{t}$ represent the three axes angular velocities at time t, and are positive in the counter-clockwise direction. We use P to represent the star spot position at time t, and OP is its vector under the star tracker coordinate system. $P^{'}$ represents the image position of the same star after movement of the star tracker at time $t + Δ t$ .

Fig. 2 Position analysis and motion model of star spot.

Download Full Size | PDF

According to Eq. (2), the relationship between vector $w_{i}$ and vector $v_{i}$ can be expressed in Eq. (3). $A^{t}$ denotes the attitude matrix at time t, while $A_{t}^{t + Δ t}$ is the attitude transformation matrix from t to $t + Δ t$ .

{\begin{cases} \vec{O_{C} P} = w_{i}^{t} = A^{t} v_{i} \\ \vec{O_{C} p^{'}} = w_{i}^{t + Δ t} = A_{t}^{t + Δ t} A^{t} v_{t} \end{cases} .

Then, the relationship between $w_{i}^{t}$ and $w_{i}^{t + Δ t}$ in Eq. (4) can be obtained in terms of Eq. (3).

w_{i}^{t + Δ t} = A_{t}^{t + Δ t} w_{i}^{t} .

The attitude transformation matrix $A_{t}^{t + Δ t}$ can be calculated from Eq. (5), and ${\tilde{ω}}^{t}$ is the cross product matrix populated with the components of the angular velocity vector $ω^{t}$ .

A_{t}^{t + Δ t} = I - {\tilde{ω}}^{t} Δ t = I - [\begin{matrix} 0 & - ω_{z}^{t} & ω_{y}^{t} \\ ω_{z}^{t} & 0 & - ω_{x}^{t} \\ - ω_{y}^{t} & ω_{x}^{t} & 0 \end{matrix}] Δ t = [\begin{matrix} 1 & ω_{z}^{t} Δ t & - ω_{y}^{t} Δ t \\ - ω_{z}^{t} Δ t & 1 & ω_{x}^{t} Δ t \\ ω_{y}^{t} Δ t & - ω_{x}^{t} Δ t & 1 \end{matrix}] .

Considering that $w_{i}^{t}$ has related to the position of star spot, a motion position model can be established as shown in Eq. (6).

{\begin{cases} x_{i}^{t + Δ t} = \frac{x_{i}^{t} + y_{i}^{t} ω_{z}^{t} Δ t + L_{f} ω_{y}^{t} Δ t}{(- x_{i}^{t} ω_{y}^{t} Δ t + y_{i}^{t} ω_{x}^{t} Δ t) / L_{f} + 1} \\ y_{i}^{t + Δ t} = \frac{y_{i}^{t} - x_{i}^{t} ω_{z}^{t} Δ t - L_{f} ω_{x}^{t} Δ t}{(- x_{i}^{t} ω_{y}^{t} Δ t + y_{i}^{t} ω_{x}^{t} Δ t) / L_{f} + 1} \end{cases} .

According to the parameters of the star tracker in this paper, $(- x_{i}^{t} ω_{y}^{t} Δ t + y_{i}^{t} ω_{x}^{t} Δ t) / L_{f} \leq 0.001$ . In this case, the value of one-step error is close to negligible, and Eq. (6) can be simplified to Eq. (7).

{\begin{cases} x_{i}^{t + Δ t} = x_{i}^{t} + (y_{i}^{t} ω_{z}^{t} + L_{f} ω_{y}^{t}) Δ t \\ y_{i}^{t + Δ t} = y_{i}^{t} - (x_{i}^{t} ω_{z}^{t} + L_{f} ω_{x}^{t}) Δ t \end{cases} .

3. Restoration method of motion-blurred star image

A crucial problem of image restoration lies in establishing the degradation model. Gonzalez et al [11–13] propose that the process of degradation of the image can be considered as the original image $c_{(x, y)}$ is degraded to $g_{(x, y)}$ through the system H and noise $η_{(x, y)}$ , as shown in Fig. 3.

Fig. 3 Degradation and restoration process model.

Download Full Size | PDF

The value of noise $η_{(x, y)}$ is random. The process can be written using common convolution symbols as follows:

g_{(x, y)} = h_{(x, y)} * c_{(x, y)} + η_{(x, y)} .

The process in the frequency domain is expressed as:

G_{(u, v)} = H_{(u, v)} C_{(u, v)} + N_{(u, v)} .

The matrix form is:

g = H c + η .

The restoration method proposed in this paper focuses on the estimation of dynamic degradation model of the star spot and restoration filter as shown in the block with dash line in Fig. 3. In the process of estimating the dynamic degradation model, motion mode in Section 2.2 is utilized to combine with the static energy diffusion of the star spot. The detailed descriptions are in the following section.

3.1 Dynamic degradation model of the star spot

Stars are far from the earth, thus, starlight rays are considered as parallel light rays and can converge to a point on the focal plane. However, with the existence of the PSF, the energy of the star spot on the image plane can be considered as Gaussian diffusion [14]. The star spot energy Gaussian distribution in the static case is shown in Eq. (11):

\begin{array}{l} g_{(x, y)} = \int_{0}^{Δ t} f (x - x_{C}, y - y_{C}) d t \\ = \int_{0}^{Δ t} (\frac{E_{s u m}^{M_{v}}}{2 π σ_{P S F}^{2}} \exp [- \frac{{(x - x_{C})}^{2}}{2 σ_{P S F}^{2}}] \exp [- \frac{{(y - y_{C})}^{2}}{2 σ_{P S F}^{2}}]) d t . \end{array}

Where $(x_{C}, y_{C})$ represents the real center position of the star spot, $σ_{P S F}$ is the Gaussian radius which represents the PSF of the energy concentration. $E_{s u m}^{M_{v}}$ is the energy-gray coefficient, related to the apparent magnitude of the corresponding star, the quantum efficiency, the integral time, the lens aperture and the optical transmittance.

The star spot energy model in dynamic case [15,16] can be expressed as:

g_{(x, y)} = \int_{0}^{Δ t} f (x - x_{C} (t), y - y_{C} (t)) d t .

Fourier transform of Eq. (11) is obtained as:

\begin{array}{l} G_{(u, v)} = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} [\int_{0}^{Δ t} f (x - x_{C} (t), y - y_{C} (t)) d t] e^{- j 2 π (u x + v y)} d x d y \\ = F_{(u, v)} \int_{0}^{Δ t} e^{- j 2 π [u x_{C} (t) + v y_{C} (t)]} d t . \end{array}

The degradation function $H (u, v)$ can be expressed as:

H_{(u, v)} = \int_{0}^{Δ t} e^{- j 2 π [u x_{C} (t) + v y_{C} (t)]} d t .

Where $(x_{C} (t), y_{C} (t))$ represent the real center position of the star spot at time t, and is in accordance with the movement in Eq. (7), the angular velocity information can be obtained using the gyroscope.

3.2 Restoration method under conditions of variable angular velocity and long exposure time

Image restoration methods include linear and nonlinear restoration methods. Inverse filter, Wiener filter and least-squares filter are common linear restoration methods. Compared with linear restoration methods, nonlinear restoration methods have advantages in suppressing noise amplification and preserving image edge information. In this work, we use a nonlinear restoration filter reference to the Lucy-Richardson filter [17,18]. Lucy-Richardson filter is an effective and applicable algorithm for processing Poisson noise in a star image, and can avoid ringing artifacts that usually exist in linear restoration methods.

In many cases, the image needs to be modeled with the random Poisson. If a random variable X has a Poisson distribution, that is to say its value of the probability can be expressed using Eq. (15):

P (X = k) = \frac{λ^{k} e^{- λ}}{k!}, 0 \leq k < \infty .

Here, we use one-dimensional expression of an image in order to simplify the description. Degradation discrete model can be expressed as Eq. (16). n represents the sequence number of a single pixel. $g_{(n)}$ is the gray value of the nth pixel of blurred image g. $\sum_{i} h_{(n - i)} c_{(i)}$ is the convolution operation result of the degradation function h and the gray value of original image c at the nth pixel.

g_{(n)} = \sum_{i} h_{(n - i)} c_{(i)} + ξ_{(n)} .

When the original image c is given, the distribution function P(g|c) of g is in Eqs. (17) and (18).

a_{(n)} = \sum_{i} h_{(n - i)} c_{(i)} .

If pixels are independent, that is:

P (g | c) = \prod_{n} \frac{a_{(n)}^{g_{(n)}} e^{- a_{(n)}}}{g_{(n)}!} .

$Π$ is the operator of continued product.

According to the joint distribution, maximum-likelihood method can be used to estimate g. Taking the logarithm of above Eq. (18) is as followed:

L = \ln P (g | c) .

L is differential with respect to $c_{(k)}$ , and $\frac{\partial L}{\partial c_{(k)}}$ is set to zero:

\frac{\partial}{\partial c_{(k)}} \ln P (g | c) = \sum_{n} (h_{(n - k)} \frac{g_{(n)}}{a_{(n)}} - h_{(n - k)}) = 0,

or

\sum_{n} h_{(n - k)} (\frac{g_{(n)}}{a_{(n)}} - 1) = 0, k = 0, 1, ..., N - 1.

To facilitate seeking c, Meinel [19]- [20] recommended the use of a multiplication iterative algorithm, the iterative solution in jth iteration is given by:

{\overset{⌢}{c}}_{(k)}^{j + 1} = {\overset{⌢}{c}}_{(k)}^{j} {(\sum_{n} h_{(n - k)} \frac{g_{(n)}}{\sum_{i} h_{(n - i)} {\overset{⌢}{c}}_{(i)}^{j}})}^{m}, k = 0, 1, ..., N - 1.

When m = 1, it is Lucy-Richardson algorithm which is used in this paper. The degenerate function h has been obtained in Section 3.1. The calculation process of restoration method can be summarized in Fig. 4.

Fig. 4 Flow chart of restoration method.

Download Full Size | PDF

Through the above Lucy-Richardson filter method, the energy of a star spot can be restored and concentrated. Thus, the star extraction method can be performed based on the restored star image. In this manner the star spot position can be acquired and can provide more accurate measurement values for further star identification and attitude determination.

4. Simulations and analysis

In order to verify the smearing model and restoration method under conditions of variable angular velocity and long exposure time, simulations and experiments are implemented. The energy and position accuracy of the motion-blurred star spot are analyzed and compared before and after restoration.

4.1 Star image simulation

Static and dynamic simulated star images are shown in Figs. 5(a) and 5(b). Partial detailed views of one star spot under the two conditions are shown in Figs. 5(c) and 5(d). The longitude and latitude of the bore-sight direction are set as (235.279°, 33.160°). The guide star catalog is selected from Tycho2n. The exposure time is 0.192s. The simulation takes into account the FOV, focal length, pixel size, the magnitude of the star, the noise of the image, and other parameters of the star tracker [21–23] to simulate the star image as realistic as possible. The stellar magnitude model refers to Liebe’s conclusion [1]. The model is reliable and frequently-used in star image simulation. The noise model of the simulated star image refers to Hancock’s analysis [24], and is combined with the parameters of the image sensor. Further, we add motion process with variable angular velocity utilizing the content in Sections 2 and 3. The simulation is implemented in MATLAB. While acquiring the star image, 1000Hz angular velocity information can be obtained by MEMS-gyro.

Fig. 5 Static and dynamic simulated star images and their partial detailed views. (a) Static simulated star image. (b) Dynamic simulated star image with rotation of two axes ( $ω_{x} = 2 ° / s, ω_{y} = 1 ° / s, ω_{z} = 0 ° / s$ ). (c) Partial detailed view of one star spot of (a). (d) Partial detailed view of one spot of (b).

Download Full Size | PDF

4.2 Dynamic motion with variable angular velocity

Three groups of dynamic simulations are conducted at different angular velocities and different levels of noise. In each group, the variation curve of angular velocity along with time is recorded, and the original and restored star images are compared. In order to display the results clearly from different dimensions, we adopt comparisons of gray images, contour lines, histograms and the sectional shape the star spot.

Simulation 1

This simulation intends to show the energy distribution and restoration situations of the star spot under the typical condition. The initial angular velocity is set as 1°/s and the angular acceleration is set as 3°/s² as shown in Fig. 6.

Fig. 6 Curve of angular velocity and time.

Download Full Size | PDF

The following results display the intuitive comparison of a star image before and after restoration. It can be seen from Figs. 7(a) to 7(c) that, when under conditions of variable angular velocity, the star energy is not merely dispersed, but also presents non-uniform distributed trail. This effect is more demanding than that of uniform angular velocity and has a more negative impact on accuracy. When adopting the restoration method proposed in this work, as shown in Figs. 7(d) to 7(f), the gray value of the star spot increases from less that 30 to more than 150. The smearing trail area of the star spot concentrates from more than 10 pixels × 5pixels to approximately 5 pixels × 5 pixels, and the restoration effect of the star energy is significant.

Fig. 7 Comparison between original star image and restored star image. (a), (b) and (c) are gray images, contour lines, and histograms of the original star image. (d), (e) and (f) are gray images, contour lines, and histograms of the restored star image.

Download Full Size | PDF

Figure 8 illustrates that asymmetric star energy distribution can be restored to ideal Gaussian distribution under conditions of variable angular velocity with the proposed approach. Thus this method can avoid the effects of motion as far as possible and maintain the high-accuracy performance of the star tracker at the level of static conditions.

Fig. 8 Energy distribution of the star spot. (a) is sectional shape of original star image. (b) is sectional shape of restored star image.

Download Full Size | PDF

Simulation 2

This simulation intends to display the energy distribution and restoration situations under conditions of large angular velocity and large angular acceleration. The initial angular velocity is set as 2°/s and the angular acceleration is set as 4°/s² as shown in Fig. 9.

Fig. 9 Curve of angular velocity and time.

Download Full Size | PDF

From Fig. 10, we can see that even under conditions of large angular and large angular acceleration, the restoration is sufficiently effective to ensure the star energy.

Fig. 10 Comparison between original star image and restored star image. (a), (b) and (c) are gray images, contour lines, and histograms of the original star image. (d), (e) and (f) are gray images, contour lines, and histograms of the restored star image.

Download Full Size | PDF

Figure 11 shows that, the star can be recovered from non-uniform distribution to Gaussian distribution. This has great significance for star extraction and position accuracy.

Fig. 11 Energy distribution of the star spot. (a) is sectional shape of original star image. (b) is sectional shape of restored star image.

Download Full Size | PDF

Simulation 3

This simulation intends to compare the results of different restoration methods when Gaussian noise exists. The initial angular velocity is set as 1°/s and the angular acceleration is set as 3°/s². Figures 12(b) and 12(e) are results of Wiener filter, and Figs. 12(c) and 12(f) are results of our restoration method.

Fig. 12 Comparison between original star image and restored star image. (a), (b) and (c) are comparison of gray star images. (d), (e) and (f) are comparison of histograms.

Download Full Size | PDF

When noise exists in the system, a considerable part of the original image may be submerged, bringing difficulty and error in star centroid determination, whereas the star spot has dramatically increased SNR after restoration. However, frequently-used restoration methods such as Wiener filter and least-squares filter will bring the ringing phenomenon while recovering star energy. This will bring the interference on the subsequent processing as star extraction. The proposed restoration method has good inhibitory effect on the noise, and is very suitable for the restoration of the star spot of the star tracker.

On account of above simulation results, we can draw conclusions that the energy and distribution of the star spot can be guaranteed under the conditions of large angular velocity, large angular acceleration and Gaussian noise, based on the proposed degradation model and restoration method.

4.3 Restoration accuracy discussion

In the restoration accuracy discussion, we conduct three groups of simulations. The simulations show the continuous imaging process of the star tracker. The simulation starts from the angular velocity of 0°/s, until approximately 4°/s. In this work, our accuracy is defined as the difference value between the motion distance of the star spot of two consecutive simulated images and the motion distance calculated from angular information. This difference value can represent the position accuracy of the star spot of the star tracker well.