Image Degradation Model for Dynamic Star Maps in Multiple Scenarios

Abstract

To meet the ground test requirements of star sensors, we establish the star map simulation algorithm and the interactive interface in multiple scenarios. The combination of the degradation model of star points, the imaging noise model, and the attitude disturbance model is introduced to solve the problem of different patterns of noise existing in the actual measurement, improving the traditional simulation model. In addition, a user-friendly interface design makes it easier for both scholars and average individuals to understand the parameters and then generate static single-frame star maps—or a series of dynamic sequence star maps—under various conditions. The results of the proposed star map simulation method are highly comparable to the actual captured star images, and this method can be applied for the tests and calibrations of star sensors.

1. Introduction

Star sensors play an essential role in aerospace attitude detection. Star sensors need to undergo rigorous performance tests before the engineering of applications. Generally, there are three ways to calibrate star sensors with high accuracy: numerical simulation, hardware simulation, and outfield observation experiments. Outfield observation experiments are conducted by placing a star sensor in an open area with low atmospheric turbulence and relatively weak stray light. Such an approach is costly and time-consuming and is usually only used as the last step in a star sensor validation experiment. Hardware simulation refers to utilizing a starfield simulator in a laboratory environment to simulate an infinity star point. The star sensor carries out identification and attitude calculation by observing the simulated star images. The disadvantages of this method are that the starfield simulator is expensive, and the parameters of the starfield simulator, such as field of view and aperture—which are difficult to match precisely with the star sensor—are not adjustable. Therefore, the versatility of the hardware simulation is limited. Compared with these two methods, numerical simulation is a highly flexible and the least resource-intensive approach. In routine performance tuning or laboratory experiments, simulated star maps are widely used because of their low cost and ease of use [1,2,3,4].

Many scholars have invested in studying the use of star sensors from digital simulation, i.e., the study of star map simulation algorithms. In a star map simulation, the establishment of the noise model is essential. To solve this problem, some scholars simulated the Gaussian grayscale distribution of star point images by considering the stellar image shift caused by the satellite motion [5]. Some scholars have introduced new theories, such as the effect of rotation around the optical axis on star point imaging, the star color index, and the idea of black-body radiation, to evaluate the quality of star points in a simulation result [6,7]. To enlarge the simulated scenarios in the simulation system, nebulae, moonlight, and Earth-obscured background are added into the field of view in the star sensors [8,9]. The literature [10] proposed a fast pixel-discretization algorithm based on the convolutional surface model to simulate dynamic real-time star maps and established a complete description of the star point motion trajectory model. However, the parameters in this algorithm are difficult to understand, limiting its usage. Researchers [11] designed a star map simulation algorithm for arbitrary exposure time length in a high-fidelity manner. However, this method is performed in an ideal state, without incorporating various error factors into the model.

Here, we design a set of algorithms and interfaces for simulating star maps in multiple scenarios from a numerical simulation approach. The algorithm includes the simulation of static single-frame star maps and dynamic sequence star maps, taking into account various noise sources. The noise may come from the electronic noise caused by the imaging system, the overall noise amplitude variation due to atmospheric condition or stray light, and the attitude disturbance caused by the different motion states between the imaging platform and the tracking star target. With various parameter inputs, noise levels and dynamic trajectories are visualized, and star maps close to the actual images are generated, which can provide a source of data for the parametric testing of the star sensor.

2. Motion Blur Degradation Model

2.1. The Process of Image Degradation

For a ground-based star sensor, during the process of image acquisition, there will be many challenges, such as the shift and superposition of light beams. All these difficulties can result in the degradation of image quality [12]. Figure 1 shows the general model of image degradation.

In Figure 1, if h (x, y) is a linear spatially invariant process, then the degraded image in the spatial domain is represented as:

g (x, y) = h (x, y) × f (x, y) + n (x, y)

(1)

in which f (x, y) represents a static, two-dimensional image, which degenerates to g (x, y) under the interference of additive noise n (x, y) through blur kernel h (x, y). * denotes a general convolution operator.

2.2. Simulations of Defocus Factor

In general, h (x, y) in the degradation model includes defocus blur and motion blur. To realize the out-of-focus impact more accurately, we use a point spread function model to simulate the image degradation with the following formula [13,14]:

$$g_{ij}=\frac{\mathrm{A}}{{2\pi \sigma }^2}\underset{\mathrm{i}-0.5}{\overset{\mathrm{i}+0.5}{{\displaystyle \int }}}\underset{\mathrm{j}-0.5}{\overset{\mathrm{j}+0.5}{{\displaystyle \int }}}exp[-\frac{{{(\mathrm{x}}_{\mathrm{i},\mathrm{j}}{ - \mathrm{x}}_{\mathrm{m}})}^2{+(\mathrm{y}}_{\mathrm{i},\mathrm{j}}{ - \mathrm{y}}_{\mathrm{m}}{)}^2}{{2\mathsf{\sigma }}^2}]\mathrm{dxdy}$$

(2)

where A denotes the energy grayscale coefficient, which is related to the total illumination, i.e., magnitude, of the imaging point on the photosensitive surface of the star sensor. σ denotes the Gaussian dispersion radius of the energy distribution in the star point region, indicating the degree of defocusing. (x_i,j, y_i,j) is any pixel within the range of the diffuse pixel point. (x_m, y_m) denotes the projected position of the star on the imaging plane of the star sensor. We can evaluate a static star point in a simulation in terms of two key parameters: the gray energy factor A and the Gaussian dispersion radius σ [15].

Table 1. Parameters of the simulated image.

Position	A	σx	σy
(15, 15)	613	1.0	1.0
(30, 30)	2113	2.0	1.5
(40, 10)	4113	2.5	3.0

lists the parameters of the simulated image in Figure 2, which shows the combination of different values of A and σ, representing various results of image degradation that may occur under real circumstances.

2.3. Integrated Gaussian Noise Model

In addition to the simulation for the background star point and target star point, the star sensor simulation also needs to simulate the effect of noise to improve the simulation’s authenticity. The random noise contained in the star map is mainly ambient noise and detector noise, such as dark current noise and output noise. These types of noise can be defined as Gaussian noise, due to their randomness [16]. The literature [2,7] states that the simulation of ground-based imaging systems must consider the effects of atmospheric influences and uneven background illumination, such stray light, moonlight, and sunlight, on the captured images. These effects will lead to a change in the overall brightness value of the image from weak to strong. Here, we use Equation (3) to represent the integrated Gaussian noise model:

N_img = f (x, y) × ψ× η (0, std), η~Gaussian (0, std)

(3)

where f (x, y) is the original image. ψ is designed to control the shape of the noise area, and it can be expressed by functions. Here is an example of a linear expression. If we set ψ_y (w, k) = w × y + k, a mathematical model of monotonic increase or monotonic decrease, the image result will have a gradient effect on the y-axis. That is, the magnitude of the image will change from dark to light or from light to dark. ψ can also be expressed in terms of other functions, such as polar coordinates. With the use of polar coordinates, a circular, radial noise pattern can be clearly expressed. η (0, std) is the amount that controls the magnitude of the noise, which follows a Gaussian distribution. The default mean of η is 0, and std represents the standard deviation. The larger the value of std, the larger the range of noise fluctuation.

Figure 3 shows the simulations of the integrated Gaussian noise model with different inputs. The simulated image is 8 bits, and the image array is 512 × 512, where the Gaussian dispersion radius of the star point is between 1 and 4. The star point circled in the figures is the same star point (with the same A, σ, and location) in the three images.

3. Simulations for Dynamic Sequence Star Maps

3.1. The Energy Distribution Model of Trajectory

Because of the relative motion between the ground-based imaging system and the celestial body, this study will complete a star map simulation in a dynamic environment. In the dynamic environment, the target star point will be dragged to form a trajectory on the imaging plane during the camera exposure time, which is called a “smearing star point.” The length and direction of this trailing path can reflect the motion information of the target star point [17]. To summarize, if there is a relative motion between the star point and the ground-based imaging system, the target is in the form of a trajectory line segment. Particularly if the imaging system has tracked and located a target object, the target star point appears relatively stationary in the image.

To achieve a star streak quickly, the exposure time T is divided into N segments at equal intervals. Accordingly, the trailing path of the star point is also divided into N segments. When interval δ = T/N is short enough, each part of the trailing path can be approximated as a point. Through such an approximation process, the dynamic trajectory of the star point can be approximated as a superposition of N static points. We set the total energy during the exposure time of the star point to ${\mathrm{E}}_0={\displaystyle \sum }_{\mathrm{n}=1}^{\mathrm{n}=\mathrm{N}}{\mathrm{E}}_{\mathrm{i},\mathrm{i}+1}$, and the energy distribution of each static star point simulating the dynamic trajectory is shown in Equation (4). In the interval from t = i to t = i + 1, the simulation formula for the energy of a star point is:

$${\mathrm{E}}_{\mathrm{i}}\left(\mathrm{x}, \mathrm{y}\right)=\frac{{\mathrm{E}}_{\mathrm{i},\mathrm{i}+1}}{{2\mathsf{\pi }\mathsf{\sigma }}^2}{\int }^{ }{\int }^{ }\exp [-\frac{{{(\mathrm{x} - \mathrm{x}}_{\mathrm{i}})}^2{+(\mathrm{y} - \mathrm{y}}_{\mathrm{i}}{)}^2}{{2\mathsf{\sigma }}^2}]\mathrm{dxdy}$$

(4)

Assuming that the dispersed spots of the above Gaussian distribution are equal in radius, the energy distribution function of the star image in each small exposure time is superimposed, and the total energy distribution function of the star image is as follows:

$$\mathrm{E}(\mathrm{x}, \mathrm{y})={\sum }_{\mathrm{n}=1}^{\mathrm{n}=\mathrm{N}}{\mathrm{E}}_{\mathrm{i}}(\mathrm{x}, \mathrm{y})$$

(5)

3.2. Attitude Disturbance Model

When the ground-based imaging system tracks a moving target star point, there are attitude disturbances between the measurement. The attitude disturbance mainly refers to the jitter of the ground-based platform brought about by its working condition and the motion lag of the camera system when it does not accurately track the motion of the target object.

Figure 4 shows a schematic diagram of the projection trajectory of the projection of the star point on the imaging plane under the influence of attitude disturbance. When rotating in any direction within O_sx_sy_sz_s, the coordinate system measured around the star sensor, the angular velocity can be decomposed onto the 3 axes of O_sx_sy_sz_s. By setting the star image on the imaging plane Oy on the fixed axis of the uniform angular velocity motion along the clockwise direction, the angular velocity is recorded as ω_z. By setting the ith instant in time, the star point position vector R_i = (x_i, y_i), $\left|{\mathrm{R}}_{\mathrm{i}}\right|=\sqrt{{\mathrm{x}}_{\mathrm{i}}{}^2{+\mathrm{y}}_{\mathrm{i}}{}^2}$. The angle between the R_i and the Oy axis is denoted θ_i, so the velocity of the star at the ith instant in time can be expressed as:

$$\{\begin{array}{l}{\mathrm{v}}_{\mathrm{ix}}{ = \mathsf{\omega }}_{\mathrm{z}}\left|{\mathrm{R}}_{\mathrm{i}}\right|{\cos \mathsf{\theta }}_{\mathrm{i}}{ + \mathrm{v}}_{0\mathrm{x}}{ = \mathsf{\omega }}_{\mathrm{z}}{\mathrm{y}}_{\mathrm{i}}{ + \mathrm{v}}_{0\mathrm{x}}\\ {\mathrm{v}}_{\mathrm{iy}}{ = -\mathsf{\omega }}_{\mathrm{z}}\left|{\mathrm{R}}_{\mathrm{i}}\right|{\sin \mathsf{\theta }}_{\mathrm{i}}{ + \mathrm{v}}_{0\mathrm{y}}{ = -\mathsf{\omega }}_{\mathrm{z}}{\mathrm{x}}_{\mathrm{i}}{ + \mathrm{v}}_{0\mathrm{y}}\end{array}$$

(6)

Given the time elapsed from the ith instant in time to the i + 1st instant in time, the amount of star point displacement [18] can be expressed as:

$$\{\begin{array}{c}{\Delta \mathrm{x}}_{\mathrm{i}}{=\mathrm{x}}_{\mathrm{i}+1}{ - \mathrm{x}}_{\mathrm{i}}{=\mathrm{v}}_{\mathrm{ix}} ·{ \Delta \mathrm{t}=(\mathsf{\omega }}_{\mathrm{z}}{\mathrm{y}}_{\mathrm{i}}{+\mathrm{v}}_{0\mathrm{x}}) · \Delta \mathrm{t}\\ {\Delta \mathrm{y}}_{\mathrm{i}}{\mathrm{i}=\mathrm{y}}_{\mathrm{i}+1}{ - \mathrm{y}}_{\mathrm{i}}{=\mathrm{v}}_{\mathrm{iy}} ·{ \Delta \mathrm{t}=(-\mathsf{\omega }}_{\mathrm{z}}{\mathrm{x}}_{\mathrm{i}}{+\mathrm{v}}_{0\mathrm{y}}) · \Delta \mathrm{t}\end{array}$$

(7)

To achieve the effect of motion blur, we propose that ${\mathrm{N}}_{\Delta \mathrm{x}}$ be added as a blur factor to the expression of the appropriate amount of the star position:

$$\{\begin{array}{c}{\mathrm{x}}_{\mathrm{i} + 1}{ = \mathrm{x}}_{\mathrm{i}}{ + \Delta \mathrm{x}}_{\mathrm{i}}{ + \mathrm{N}}_{\Delta \mathrm{x}}{,\mathrm{N}}_{\Delta \mathrm{x}}~\mathrm{Gaussian}(0,\mathrm{n})\\ {\mathrm{y}}_{\mathrm{i} + 1}{ = \mathrm{y}}_{\mathrm{i}}{ + \Delta \mathrm{y}}_{\mathrm{i}}{ + \mathrm{N}}_{\Delta \mathrm{y}}{,\mathrm{N}}_{\Delta \mathrm{y}}~\mathrm{Gaussian}(0,\mathrm{m})\end{array}$$

(8)

where the default mean of the Gaussian distributed random number is 0. n and m represent the standard deviation of the set Gaussian random number. The values of n and m can be same. The larger the n and m, the greater the deviation of the simulated star point from the preset position.

In conclusion, the description of the star trajectory can be expressed by Equation (9):

$$\{\begin{array}{c}{x}_i\left(t\right)=x_{i0}+{\int }_{t0}^{t0+T}{v}_{ix}\left(\tau \right)d\tau \\ y_i\left(t\right)=y_{i0}+{\int }_{t0}^{t0+T}{v}_{iy}\left(\tau \right)d\tau \end{array}$$

(9)

Therefore, we use the proposed methods to simulate the dynamic trajectory of a star point. Figure 5, Figure 6 and Figure 7 show the superposition of the different states when the simulated star trajectory is a linear equation according to $\{\begin{array}{c}{\mathrm{x}}_{\mathrm{i}+1}{ = \mathrm{x}}_{\mathrm{i}}{ + 1 + \mathrm{N}}_{\Delta \mathrm{x}}\\ {\mathrm{y}}_{\mathrm{i}+1}{ = \mathrm{y}}_{\mathrm{i}}{ + 1 + \mathrm{N}}_{\Delta \mathrm{y}}\end{array}$. Figure 8 shows the superposition of the different states when the simulated star trajectory is a curvilinear equation according to $\{\begin{array}{c}{\mathrm{x}}_{\mathrm{i}+1}{ = \mathrm{x}}_{\mathrm{i}}{ + 1 + \mathrm{N}}_{\Delta \mathrm{x}}\\ {\mathrm{y}}_{\mathrm{i}+1} = -0{.006\mathrm{x}}_{\mathrm{i}}{}^2 + 2{.7\mathrm{x}}_{\mathrm{i}}{ - 27 + \mathrm{N}}_{\Delta \mathrm{y}}\end{array}$.

4. Interface Function Description

The field of view of the simulated detector in the proposed method is $12° \times 12°$, and the size of its plane array is $512 \mathrm{pixels} \times 512 \mathrm{pixels}$. The magnitude sensitivity is 8.0 apparent magnitude, the simulated imaging system has a frame rate of 50 Hz, and the exposure time for each frame is 20 ms. The interface uses the Tkinter module based on Python 3.9.

There are two functions in the designed interface. First, it is the static single-frame star map simulation. Figure 9 shows the initial interface of this function. In this state, the default number of the target star point is 1. The adjustable parameters are the Gaussian noise level of the image, the number of background stars, the maximum value of target brightness, and the minimum value of target brightness. Users can generate single-frame simulated star maps, with specified parameters, in this mode. Additionally, a function to display images is added, which can also be used separately as a picture viewer. After generating the simulated star maps, the primary information of the target star point, such as brightness information, position coordinates, and Gaussian dispersion radius, can be viewed in “Target Star Information.”

The second function is the dynamic sequence star map simulation, and its user interface is shown in Figure 10. In this mode, the simulation can provide two motion states for a target star point. One of the states is that the ground-based camera system has followed the target point to be measured. Under this state, by default, the target in the first frame is in the center of the image. In addition to the size and brightness of a star point, it is also possible to simulate errors in the actual situation by inputting the magnitude of the attitude disturbance of a target point and the average and standard deviation of Gaussian noise. The default exposure time of each frame of the system is 20 ms. By adjusting the size of the “Dynamic Frame Count,” users can input the imaging time of the simulated ground-based camera system. The other state is that the ground-based camera system has not followed the target point to be measured. Under this state, both the target and background star points have different movement direction speed sizes, compared with the ground-based camera system. The user must input the x-direction and y-direction offset to complete the exact movement simulation of the target star point. Since it would take time to generate a large number of background star points in real-time to meet the requirements, the way to handle the background star points is by loading a pre-generated binary file. This file is the sum of a specified number of star trajectories laid out according to a preset angle. The user can also change this, according to specific requirements.

Figure 11 shows the first, the 90th, and the 150th frames of the 200-frame sequence of star maps simulated under the conditions of the selected states and parameters depicted in Figure 10. The target to be measured is marked with a red circle, and the coordinates of the target star point in different frames are listed in the white textboxes.

5. Results

5.1. Image Quality Assessment

Figure 12 shows simulation results using the method proposed in this paper. Figure 12a,b shows the static single-frame star map simulation results, and Figure 12c shows the dynamic sequence star map. Figure 13 displays actual photos from the Shanghai Sheshan Observatory and from astronomer Michael A. Earl. From the gray histograms of the two types of images, we can determine that the actual photos have a broader gray distribution than that of the simulation results. Furthermore, there are some stripes in the photos. These stripes may come from the telescope lens.

Here, we use two methods to evaluate the image similarity between the simulation results and actual photos. Figure 14 demonstrates the gray histograms of the simulated images in Figure 12. And Figure 15 demonstrates the gray histograms of the real images in Figure 13. First, the Bhattacharyya coefficient which is based on the gray histogram is used to indicate the color distribution of an image. The coefficient ranges from 0 to 1. The closer to 1, the more similar the two images are proved to be in terms of color distribution. The second standard is the cosine similarity, which changes an image into vectors by adding the divided gray level area numbers [19]. The calculation results are shown in

Table 2. Calculation results of image similarity.

Number	Bhattacharyya Coefficient in Gray Histogram	Cosine Similarity
Figure 12a vs. Figure 13a	0.778	0.977
Figure 12b vs. Figure 13b	0.639	0.992
Figure 12c vs. Figure 13c	0.639	0.893

, which point out that there is still a gap between the grayscale distribution of the simulated results and the actual photos, but the overall character of both is very close.

Figure 12c shows the effect of a 2 s imaging time in a dynamic environment simulated using the proposed method’s noise parameter ${\mathrm{N}}_{\mathrm{img}}=\mathrm{f}(\mathrm{x}, \mathrm{y})·1·{\mathsf{\eta }(0,15),\mathrm{N}}_{\Delta \mathrm{x},}{\mathrm{N}}_{\Delta \mathrm{y},}~\mathrm{Gaussian}(0,3)$. Figure 13c shows an actual image taken by Canadian astronomer Michael A. Earl with a ground-based telescope. These are some of the differences between Figure 12c and Figure 13c and the main reasons behind them. Besides the direction of the background trajectories, which depends on the movement states of the simulated conditions, the trajectory style of the target is more noticeable in the actual photos. Because the satellite tumbles at a regular interval within an extended period, the satellite’s orbit produces a periodic flickering phenomenon, which is not addressed by the method in this paper. Next, the noise distribution in Figure 13c is more complicated than in Figure 12c. Figure 13c becomes brighter due to the height of the image, which means the average of Gaussian noise becomes bigger. which means the average of Gaussian noise becomes bigger. In addition, some black dots and black areas exist in the figure in a more obvious way. It is necessary to determine different characteristics and modes of noise to improve the quality of the simulation results.

5.2. Star Centroid Extraction

The centroid of stars is the most important factor in attitude determination. We adopt the common centroid method to evaluate the accuracy of the simulation result of the proposed method. Figure 16a shows an actual image taken by Michael A. Earl. Four distinct targets are circled in the image, and the simulation result is shown in Figure 16b.

Next, before we perform the star centroid extraction process, we need to implement the Gaussian filter to reduce the interference from background noise in Figure 16a. The Gaussian kernel size is $5 \times 5$ and the Gaussian kernel standard deviation is computed with the formula provided by OpenCV (getGaussianKernel(), https://docs.opencv.org/4.0.0/d4/d86/group__imgproc__filter.html#gac05a120c1ae92a6060dd0db190a61afa accessed on 15 November 2021). Finally, we adopt the common centroid method to extract the “center of mass” as follows:

$$\{\begin{array}{c}{\mathrm{x}}_0=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{m}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{ij}}}{{\sum }_{\mathrm{i}=\mathrm{1}}^{\mathrm{m}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{f}}_{\mathrm{ij}}}\\ {\mathrm{y}}_0=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{m}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{y}}_{\mathrm{j}}{\mathrm{f}}_{\mathrm{ij}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{m}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{f}}_{\mathrm{ij}}}\end{array}$$

(10)

where f_ij denotes the value of the image after background noise reduction.

The same steps will apply to the detection of Figure 16b. The results for the two images are shown in

Table 3. Results of star centroid extraction in Figure 16a,b.

	Figure 16a		Figure 16b		Error
No.	x	y	x	y	Δx	Δy
1.	7.329	121.973	7.498	121.983	0.169	0.01
2.	52.113	267.078	52.037	266.969	−0.076	−0.109
3.	163.037	187.762	162.989	187.964	−0.048	0.202
4.	248.171	220.978	247.978	220.970	−0.193	−0.008

.

Where the errors are calculated as the follow equation:

$$\{\begin{array}{c}\Delta \mathrm{x}={\mathrm{x}}_{\mathsf{\beta }} - {\mathrm{x}}_{\mathsf{\alpha }}\\ \Delta \mathrm{y}={\mathrm{y}}_{\mathsf{\beta }} - {\mathrm{y}}_{\mathsf{\alpha }}\end{array}$$

(11)

where x_β and y_β denote the x and y coordinate values of the stars in Figure 16b. x_α and y_α denote the x and y coordinate values of the stars in Figure 16a.

As we can see, because of the noise of figures, the accuracy is not perfect, but it is acceptable for common use. If we perform the image processing in detail and apply a more specific centroid extraction method, the error will be smaller.

6. Conclusions

In summary, this paper carries out a simulation method for dynamic sequential star maps of space targets in multiple scenarios with complex star backgrounds. Numerical simulation provides a more convenient way to generate the abundant data of star maps to promote the study of star sensors. A static star point model and a dynamic sequence star map model are constructed, respectively. We take the random noises, such as detector noise and dark current noise, and uneven environment illumination into account. Meanwhile, the impact of the noise model analysis and platform motion disturbances on imaging will also provide a reference for performance testing and verification of star sensors. In the next stage of the work, we may consider the condition in which multiple interference target star points exist. In addition, it is essential to continue to study the noise sources and patterns in authentic images so that the parameters can be further optimized, making them more similar to those in the actual scenario.