Simulation of Dense Star Map in Deep Space Based on Gaia Catalogue

TL;DR

This paper introduces a high-precision simulation framework for deep-space star fields using the Gaia catalog to improve space situational awareness and target detection.

Contribution

The novel integration of a full optoelectronic imaging chain with dynamic platform disturbances and realistic noise models using the Gaia catalog.

Findings

The model achieves SNR error of less than 10% when validated against ground-based telescope data.

Sub-pixel centroiding accuracy exceeds 0.01 pixels, demonstrating high simulation fidelity.

The framework supports realistic simulation of dim stars at low magnitudes with unprecedented precision.

Abstract

High-fidelity star field simulation is paramount for target detection and space situational awareness (SSA) in geostationary and deep-space environments. However, accurately modeling the synergistic effects of ultra-dense stellar backgrounds and complex platform perturbations remains a formidable challenge. This paper proposes an integrated simulation framework that leverages the Gaia catalog to generate high-precision stellar environments. The core methodological novelty lies in the end-to-end coupling of a full optoelectronic imaging chain with dynamic platform disturbances, effectively bridging the gap between theoretical orbital dynamics and realistic sensor responses. Distinguishing itself from conventional models, our approach uniquely integrates radiative transfer and high-fidelity noise suites—including photon shot noise and non-uniform stray light—while utilizing the Gaia catalog to achieve unprecedented precision in simulating dim stars at low magnitudes. The fidelity of the proposed model was quantitatively validated against empirical data from a ground-based wide-field telescope (GTC). Experimental results, derived from multiple simulation realizations, demonstrate high consistency with real-world observations, achieving a Signal-to-Noise Ratio (SNR) error of less than 10% and a sub-pixel centroiding accuracy exceeding 0.01 pixels. This work provides a robust, high-fidelity data synthesis tool that significantly advances the development of target detection algorithms and the performance optimization of space-based optical sensors.

Figures

Full text

1. Introduction

With the rapid development of aerospace technology, the scope of human space activities has gradually expanded from Low Earth Orbit (LEO) toward more distant outer space. This expansion brings increasingly complex security challenges. On the one hand, the proliferation of defunct satellites, rocket bodies, and fragmentation debris in high-altitude orbits poses a severe threat to operational spacecraft. Due to their long distances and faint signatures, traditional ground-based monitoring methods struggle to achieve all-time, high-precision cataloging and early warning capabilities. On the other hand, the risk of impact from deep-space threats, namely Near-Earth Objects (NEOs), is ever-present. Consequently, establishing a forward-looking deep-space defense system has become a focus for nations worldwide.

Against this background, enhancing High-Orbit Space Situational Awareness (SSA) capabilities is particularly critical. Compared to ground-based systems, which are affected by atmospheric turbulence and limited observation windows, space-based optical surveillance platforms offer inherent advantages, such as immunity to cloud cover, all-weather operation, and close-proximity observation.

However, as the detection performance of space-based optical systems improves, the influence of the stellar background also gradually increases. In the early optical era, the detection limit for stars was primarily determined by the light-gathering capability of the human eye and small telescopes; Ptolemy’s star catalog from 150 AD only recorded about 1000 stars, with a limiting magnitude of approximately 6. With the application of photographic technology, the stellar detection limit rapidly increased, and late 20th-century catalogs like the Digitized Sky Survey (DSS) [1] and Guide Star Catalog (GSC) [2] could reach a detection limit of magnitude 20. Today, with the widespread use of devices like CCD and CMOS, the Sloan Digital Sky Survey (SDSS) [3] can achieve a limiting magnitude of 22.5 in the r-band. Therefore, stronger detection capability inherently implies more stellar interference.

Existing star field simulations typically rely on catalogs like Tycho-2 or UCAC4, which are sufficient for basic attitude determination but inadequate for high-fidelity situational awareness. Moreover, conventional frameworks often employ simplified point spread functions (PSFs) and decoupled kinematic models, causing poor energy concentration and significant Signal-to-Noise Ratio (SNR) errors when rendering dense, faint star fields. Accurately detecting dim high-orbit targets demands physically rigorous simulations of these faint backgrounds.

To address this, we propose an advanced simulation framework built strictly upon a complete optoelectronic imaging chain. By utilizing the Gaia catalog, our method accurately renders low-magnitude stars, significantly enriching background complexity. Unlike traditional approaches, we tightly couple precise energy transfer and dynamic motion models directly into the imaging link. This integration yields superior physical fidelity—demonstrated by higher energy concentration and substantially lower SNR errors—providing a robust data foundation for evaluating dim-target detection algorithms in deep space.

This paper directly addresses this issue by proposing a star catalog-based simulation technique for dense stellar backgrounds.

The subsequent structure and main content of this paper are organized as follows:

1. Chapter 2 introduces and analyzes the relevant theoretical concepts involved in the simulation technique.
2. Chapter 3 details the modular construction of the actual simulation model.
3. Chapter 4 presents the experimental design of the simulation model, along with a fidelity assessment using real observational data from the Ground-based Wide-field Telescope (GTC). This chapter also provides simulation schemes and results for several specific scenarios.

2. Theoretical Analysis

This chapter establishes a comprehensive theoretical framework for a full-link physical simulation model tailored for space-based observation systems. This model serves as the foundation for subsequent implementation and performance characterization. To meet the rigorous detection requirements for high-orbit space debris and deep-space small celestial bodies, the analysis encompasses the entire physical chain—extending from orbital dynamics resolution to optoelectronic imaging conversion.

2.1. Analysis of Target Characteristics

A systematic characterization of the targets, stellar backgrounds, and imaging mechanisms is provided. As illustrated in Figure 1, the simulation explicitly accounts for the relative motion between the observation platform and space target. The arrows in the figure indicate the direction of this motion. Due to the immense distances involved, the stellar background is treated as a stationary celestial reference frame, while the dynamic interaction is primarily dictated by the satellite’s orbital maneuvers and attitude perturbations.

2.1.1. Target Motion State

To support the high-fidelity simulation, the relative motion between the observation platform and the space target is first modeled. In an ideal scenario, the interaction between the Earth and the space target is simplified as a two-body problem, where the target is assumed to move in a circular orbit. Under these assumptions, the orbital velocity $[eqn]$ and the angular velocity $[eqn]$ of the target are defined as [4]:

$$[eqn]$$

where $[eqn]$ is the Geocentric Gravitational Constant, $[eqn]$ denotes the mean radius of the Earth, and $[eqn]$ represents the orbital altitude of the target. The relative distance $[eqn]$ between the observation platform and the target can be derived using the Law of Cosines:

$$[eqn]$$

where $[eqn]$ is the distance from the Earth’s center to the observation platform (e.g., the Earth-Moon distance for specific lunar-based scenarios), and $[eqn]$ is the angular separation between the platform and the target relative to the geocenter.

In the triangle defined by the observation platform, the space target, and the Earth, the Law of Sines is applied as follows:

$$[eqn]$$

$[eqn]$ is the angle subtended by the Earth and the space target at the observation platform. Therefore, the angle $[eqn]$ , which corresponds to the tangential component of the space target relative to the observation platform, is expressed as:

$$[eqn]$$

Consequently, the angular velocity of the space target relative to the observation platform can be expressed as:

$$[eqn]$$

In the actual simulation process, we employ two modes to simulate the target’s motion trajectory: directly inserting virtual trajectories onto the image plane, and calculating the target trajectory based on TLE data and then mapping it onto the image plane. The former method involves directly generating several sets of trajectories for the target moving on the image plane (with controllable trajectory types). The latter requires calculation combining the camera’s internal and external parameters.

2.1.2. Target Energy Distribution Analysis

Due to the target’s significant distance from the detector, it forms a point source on the image plane, exhibiting an imaging distribution essentially identical to that of a star. The energy distribution of a single stellar image is described by the response function of a point source, which is the Point Spread Function (PSF). The Optical Transfer Function (OTF) is its corresponding Fourier Transform. The stellar light intensity distribution is simulated using the circular aperture diffraction model, yielding the following equation:

$$[eqn]$$

where $[eqn]$ is the peak intensity, and $[eqn]$ is the first-order Bessel function

In ideal star maps, stellar targets and faint small targets are represented as point sources. When a survey camera utilizes long integration times—typically when staring at stars to observe faint objects—targets with relative motion develop motion smear (or streaking). The motion (smear) of the streaked target, measured in pixels per frame interval, is calculated as:

$$[eqn]$$

where $[eqn]$ is the angular velocity of the target relative to the image plane; and $[eqn]$ is the angular resolution of a single detector pixel. The simulation results for the streaked target are shown in Figure 2.

2.2. Analysis of Stellar Intensity

2.2.1. Stellar Intensity Analysis

Stars are approximated as point sources located at an infinite distance from the detector. In the dark background of deep space, their radiation propagates through the imaging chain, influenced by the optical system’s Optical Transfer Function (OTF) and the resulting Point Spread Function (PSF). Consequently, the energy is finally dispersed into a Gaussian spot (or blob) on the image plane. However, the actual size of the resulting image may be affected by multiple factors, including the optical system characteristics, the detector, and the jitter (vibration) of the target and the observation platform. Its two-dimensional Gaussian distribution on the image is expressed as [5]:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the distribution means (i.e., the pixel center); and $[eqn]$ and $[eqn]$ are the distribution standard deviations (i.e., the standard deviation of the PSF).

2.2.2. Band Impact Analysis

The simulation magnitude is based on the G-band magnitude data provided in the Gaia catalog. This band is the core photometric system of the Gaia catalog, covering 3301050 nm, with the effective response mainly concentrated between 400900 nm, and peak sensitivity at approximately 673 nm. The specific value of the magnitude is calculated by integrating the star’s photon flux in the G-band [6], using the formula:

$$[eqn]$$

$[eqn]$ : The measured flux in the G-band (instrument-calibrated) $[eqn]$ : The zero-point of the G-band, calibrated using standard stars

To simulate the star’s appearance on the detector, the catalog magnitude $[eqn]$ is converted into the corresponding focal plane flux (pixel intensity) using the following relationship:

$$[eqn]$$

where $[eqn]$ is the reference flux for a zero-magnitude star in the G-band, $[eqn]$ is the effective aperture area, $[eqn]$ is the optical transmittance, $[eqn]$ is the quantum efficiency, and $[eqn]$ is the integration time. This ensures a physically consistent mapping from catalog data to simulated digital numbers (DN).

Within the defined spectral range, the diffraction limit is evaluated to justify the PSF modeling. The maximum radius of the star’s Airy disk is calculated at the upper-bound wavelength of 900 nm as follows:

$$[eqn]$$

where $[eqn]$ is the focal length and $[eqn]$ is the aperture diameter. In this simulation, these diffraction parameters are utilized to determine the kernel size of the Gaussian PSF, ensuring that the energy concentration of the simulated stars accurately reflects the theoretical resolution of the optical system.

Among the wavelength-related factors, the diffraction radius remains smaller than the pixel size even in the maximum diffraction state, satisfying the envelope requirement that the energy concentration is greater than 90%. Therefore, a central wavelength of 673 nm will be adopted for the subsequent simulation.

2.3. Imaging Model Analysis of Deep Space Exploration System

2.3.1. Analysis of Factors Affecting Geometric Performance

The satellite orbital dynamics model is utilized to generate high-precision attitude and orbital parameters, which are essential for accurate platform positioning. To ensure simulation fidelity, the model accounts for various orbital perturbations and the precise timing offsets between different astronomical time systems (e.g., UTC, TAI, and TDT). Furthermore, system compensations are integrated to account for the sensor’s relative imaging velocity and viewing angle variations.

The modeling of satellite attitude characteristics is tightly coupled with the closed-loop attitude control system (ACS), as illustrated in Figure 3. Rather than being treated as static parameters, the satellite attitude is dynamically simulated within the control loop. Importantly, the attitude control errors and high-frequency disturbances (e.g., micro-vibrations) generated by this model are directly propagated into the imaging chain. These disturbances are translated into time-varying line-of-sight (LOS) jitter and integrated over the exposure time, thereby simulating realistic motion blur and geometric distortion in the synthetic star maps. This coupling ensures that the simulated images accurately reflect the dynamic instability of the observation platform.

2.3.2. Sampling Model Analysis

Discrete sampling is an inherent characteristic of all optoelectronic imaging systems [7]. The discrete distribution of the detector elements dictates that the scene is sampled in two spatial directions. The sampling model is given by:

$$[eqn]$$

where $[eqn]$ is the spacing between adjacent detector elements in the $[eqn]$ direction; and $[eqn]$ is the spacing between adjacent detector elements in the $[eqn]$ direction.

2.3.3. Analysis of Noise Models

Noise sources include: dark current noise (detector), shot noise (detector, circuit), readout noise (detector), amplifier noise (circuit), quantization noise (circuit), and pattern noise (detector) [8]. The first five are additive noise, and the last is multiplicative noise. The models for additive and multiplicative noise are as follows:

$$[eqn]$$ $$[eqn]$$

Dark current typically originates from surface defects at the SiO_2_/Si interface, surface generation, thermal generation, and defects in the semiconductor manufacturing process. The dark current generated in a single pixel during the integration time, expressed in electrons, is:

$$[eqn]$$

where $[eqn]$ is the detector integration time; $[eqn]$ is the dark current density; $[eqn]$ is the detector pixel area; and $[eqn]$ is the charge carried by one electron.

In addition to dark current density, the pixel area and integration time also influence the magnitude of dark current noise. Since the detector randomly generates electrons, the dark signal produces additional noise known as dark current shot noise, which is therefore described by the Poisson process:

$$[eqn]$$

where $[eqn]$ , shot noise is time-related, making it temporal noise.

Dark current and shot noise are generated in the integration process of the detector, and the magnitude is related to the integration time and the size of the detector.

Readout noise is generated during the charge transfer process, is independent of the signal, and is analyzed using a white noise model. Its value is usually provided in the detector manual (readout noise, noise equivalent electrons, or noise floor).

Amplifier noise includes on-chip amplifier noise and off-chip amplifier noise. On-chip amplifier noise mainly includes $[eqn]$ noise and white noise; off-chip amplifier noise is similar to on-chip amplifier noise and is treated as white noise here.

Quantization noise is generated during AD conversion. Its noise standard deviation is inversely proportional to the number of quantization bits; the more quantization bits, the lower the noise standard deviation, as shown in the formula below. This type of noise is treated as white noise.

$$[eqn]$$

where $[eqn]$ is the saturation electron count, and $[eqn]$ is the number of quantization bits.

Pattern noise primarily includes Fixed Pattern Noise (FPN) and Photo Response Non-Uniformity (PRNU).

Using $[eqn]$ to represent dark current non-uniformity, the fixed pattern noise is:

$$[eqn]$$

Using $[eqn]$ to represent photo response non-uniformity, the PRNU noise is:

$$[eqn]$$

The two parameters, FPN and PRNU, typically refer to dark current non-uniformity and photo response non-uniformity, respectively [9]. The standard deviation used for calculation is usually obtained by dividing by 5~6 times.

The Poisson probability distribution model is:

$$[eqn]$$

The Gaussian probability distribution model is:

$$[eqn]$$

2.3.4. Quantitative Model Analysis

The Digital Number (DN) value of the quantized digital image is given by:

$$[eqn]$$

where $[eqn]$ is the electrical signal generated by the detector; LSB is the Least Significant Bit or quantization interval; and int is the integer truncation operator (or rounding-down operator). For uniform quantization, the LSB is shown below:

$$[eqn]$$

where $[eqn]$ is the saturation electron count (or full-well capacity).

The aforementioned noise sources are integrated into the simulation following the physical image-formation sequence. First, photon shot noise is modeled as a Poisson process during signal accumulation. Next, dark current and FPN are added to simulate thermal and non-uniformity effects. Finally, readout noise and quantization errors are introduced during the signal conversion and digitization stages. This step-by-step approach ensures the synthetic images accurately reflect the stochastic and systematic uncertainties of the detector.

3. The Establishment of Simulation Model

The complete image simulation process begins by inputting the orbital data of the observation satellite and targets, alongside star catalog data. Subsequently, the orbital position data for the satellite, targets, and stars are resolved. Within the deep-space background simulation module, simulation is performed based on the characteristics of the targets and the stellar background to obtain object-space simulation data. As the optical telescope in the detection system collects optical signals from space targets, non-target optical signals from the space environment (such as solar radiation) are captured simultaneously in the form of stray light; therefore, it is necessary to model the stray light entering the system [10]. Finally, the detector noise generated during the optoelectronic signal conversion process is modeled to complete the imaging link modeling of the deep-space detection system.

Regarding the simulation content, the settings are divided into objective and subjective factors. Objective factors refer to results that are not influenced by parameter inputs, primarily involving the information input and orbital position simulation modules, including platform coordinates, star catalogs, and target orbital data. Subjective factors are simulation components influenced by subjective inputs, mainly including the deep-space background simulation module, stray light modeling, and the signal electron count model, which incorporate camera parameters, noise parameters, and other variables, as shown in Figure 4.

To ensure the authenticity of the relative positions between the stars, the observation platform, and the observation targets, the inputs for the simulation model consist of public TLE data and star catalogs. For the simulation of stellar positions and energy within dense star field images, a more extensive star catalog is utilized; specifically, the Gaia catalog is adopted as the primary input source for stellar position and magnitude information.

3.1. Orbit Position Simulation Module

3.1.1. Resolution of Star Positions on the Image Plane

Stars can be considered to be at an infinite distance from both the detector and the Earth’s center. Thus, the phase angle of a star relative to the satellite in the Inertial Coordinate System is equivalent to the star’s Right Ascension (Ra) and Declination (Dec) data [11].

$$[eqn]$$

where RA is the Right Ascension in star catalog coordinates, Dec is the Declination in star catalog coordinates, and $[eqn]$ is the corresponding unit direction vector of the star relative to the detector in the Earth-Centered Inertial (ECI) frame.

After inputting Ra, Dec, and magnitude data from the Gaia catalog, the first step is to resolve the star’s position on the image plane. The calculation steps are as follows [12]:

1. RA and Declination (Dec) are the fundamental parameters of the celestial coordinate system used to describe the positions of celestial bodies on the celestial sphere.

Right Ascension (RA) is analogous to longitude on Earth; it uses the vernal equinox as the zero point and is measured eastward along the celestial equator. It is typically expressed in time units (hours) or degrees, ranging from 0° to 360°, as shown in Figure 5.

Declination (Dec) is analogous to latitude, representing the angular distance of a celestial body relative to the celestial equator, ranging from −90° to +90°.

The celestial coordinate system is a system used to describe the positions of objects on the celestial sphere. It is centered on the Earth, projecting the positions of celestial bodies onto an imaginary sphere. This system is closely related to the Earth’s axis of rotation and the equatorial plane. Conversion from RA and Dec to Celestial Rectangular Coordinates [13]:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ , $[eqn]$ : are the input Right Ascension and Declination in degrees. $[eqn]$ represents the output unit direction vector in the celestial coordinate system.

2.Transformation from Celestial Coordinates to the Satellite Inertial Coordinate System

As the satellite moves in space, its inertial reference frame is typically aligned with the celestial coordinate system (the J2000 coordinate system is adopted in this paper). In this step, the celestial rectangular coordinate system obtained in the first step can be directly utilized. Specifically, the unit vector of the star’s coordinates in the inertial reference frame is defined as: $[eqn]$

3.Transformation from the Inertial Coordinate System to the Satellite Body Coordinate System

The satellite body coordinate system is typically defined such that the optical axis is aligned with the Z-axis, while the X and Y axes lie within the image plane. The satellite’s attitude is defined by a rotation matrix and a displacement vector. The formula for mapping coordinates from the J2000 coordinate system to the satellite body coordinate system is as follows:

$$[eqn]$$

where $[eqn]$ is a 3 × 3 rotation matrix, and $[eqn]$ is a 3 × 1 translation vector describing the displacement from the origin of the reference coordinate system to the origin of the satellite coordinate system. This step yields the final output: the unit direction vector of the star expressed in the satellite body coordinate system $[eqn]$ .

4.Perspective Projection from the Satellite Coordinate System to the Image Plane Coordinate System

A perspective projection is employed to transform the three-dimensional satellite body coordinates into two-dimensional image plane coordinates, specifically onto the normalized unit plane:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ are the unit plane coordinates, representing the normalized image plane coordinates.

5.Transformation from Normalized Image Plane Coordinates to Pixel Coordinates

The normalized image plane coordinates are mapped to the actual pixel coordinates using the intrinsic parameters (or internal parameters) provided by the camera manufacturer:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ denote the ratio of the focal length to the pixel size in the x and y directions, respectively. $[eqn]$ and $[eqn]$ represent the pixel coordinates of the principal point (the point where the optical axis intersects the image plane). The final output is the pixel coordinates $[eqn]$ .

3.1.2. Target Image Plane Position Resolution

The resolution of the target position on the image plane requires mapping the target coordinates in the Earth-Centered Inertial (ECI) frame, derived from TLE data, onto the image plane coordinates [14].

Transformation from ECI Coordinate System to Orbital Coordinate System

In the ECI frame, let the target position vector be $[eqn]$ and the satellite position vector be $[eqn]$ . The target position vector transformed from ECI to ORB is given by:

$$[eqn]$$

where the coordinate rotation matrix is:

$$[eqn]$$

In this expression, $[eqn]$ denotes the orbital inclination, $[eqn]$ represents the Right Ascension of the Ascending Node (RAAN), and $[eqn]$ is the argument of latitude.

The axis swap matrix is defined as:

$$[eqn]$$

2.Transformation from Orbital Coordinate System to Satellite Body Coordinate System

The satellite body coordinate system is obtained by rotating the orbital coordinate system through three directions: roll, pitch, and yaw. Let the three Euler angles be $[eqn]$ , $[eqn]$ , $[eqn]$ , and the rotation matrix be $[eqn]$ . Then:

$$[eqn]$$

If the coordinates of the target in the orbital and body frames are $[eqn]$ and $[eqn]$ respectively, then:

$$[eqn]$$

3.Transformation from Satellite Body Coordinate System to Sensor Coordinate System

The satellite body coordinate system aligns with the sensor coordinate system after the following rotations: A counter-clockwise rotation of angle $[eqn]$ around the Z-axis; A counter-clockwise rotation of angle $[eqn]$ around the Y-axis. The corresponding rotation matrix is:

$$[eqn]$$

If the target coordinates in the body and sensor frames are $[eqn]$ and $[eqn]$ , then:

$$[eqn]$$

4.Transformation from Sensor Coordinate System to Image Plane Coordinate System

The position in the sensor coordinate system undergoes projection and discretization to yield the image plane position. Let the target coordinates be $[eqn]$ . The projected position is:

$$[eqn]$$

Discretization is performed using the following method:

$$[eqn]$$

where $[eqn]$ is the camera resolution, and $[eqn]$ denotes the floor operation (taking the lower limit).

The complete calculation pipeline for stellar image plane coordinates and the corresponding coordinate system definitions are illustrated in Figure 6.

3.1.3. High-Precision Centroid Simulation Based on Bicubic Interpolation

Compared with bilinear interpolation, which considers only a 2 × 2 neighborhood and lacks derivative continuity, bicubic interpolation employs a 4 × 4 neighborhood integrated with cubic polynomials. By fitting the intensity variation trend among 16 adjacent pixels, this approach ensures the continuity of both the reconstructed surface and its first-order derivatives. Consequently, it effectively preserves high-frequency stellar profiles and minimizes sub-pixel centroiding errors. For high-fidelity simulations leveraging the Gaia catalog, such an approach is essential to maintaining a centroid accuracy within 0.01 pixels [15]. For a target point $[eqn]$ , the interpolated pixel value $[eqn]$ is calculated as follows:

$$[eqn]$$

where $[eqn]$ represents the grayscale values of the 16 neighboring pixels, and W is the bicubic weighting function:

$$[eqn]$$

In this formula, $[eqn]$ is an adjustment parameter taken as −0.5 during the simulation, and $[eqn]$ is the horizontal or vertical distance between the target point and the neighboring pixel. Through bicubic interpolation, the centroid simulation error can be controlled within 0.01 pixels.

For a real image signal $[eqn]$ , performing a Taylor expansion at $[eqn]$ yields:

$$[eqn]$$

The approximation error of bicubic interpolation originates from [16]:

$$[eqn]$$

when the third-order derivative of $[eqn]$ is non-zero, the error is primarily determined by the fourth-order term $[eqn]$ :

$$[eqn]$$

For $[eqn]$ , the maximum error is approximately $[eqn]$ . For sub-pixel displacements (i.e., 0.01 pixels), the error term $[eqn]$ .

3.2. Energy Simulation Module

3.2.1. Stellar Energy Simulation

When simulating the starry sky in a computer, the simulated star map is displayed in the form of a two-dimensional digital image. In addition to obtaining the position of the star point on the detector image plane, it is necessary to determine the grayscale value corresponding to the pixels. This involves converting the apparent magnitude of the observed star in the star catalog into the grayscale value on the detector image plane.

Stellar magnitudes are categorized into apparent magnitude and absolute magnitude. Apparent magnitude refers to the brightness of a celestial body as seen by an observer, while absolute magnitude is the brightness as seen from a distance of 32.616 light-years. The simulation data source utilizes the apparent magnitude. Based on the relationship where the irradiance decreases by 2.51 times for every one-level increase in magnitude [17], the expression relating magnitude to irradiance is:

$$[eqn]$$

That is:

$$[eqn]$$

where $[eqn]$ is the magnitude of the $[eqn]$ -th observed star, $[eqn]$ is the irradiance of the $[eqn]$ -th observed star, and $[eqn]$ is the irradiance of a zeroth-magnitude star. Using a linear calculation method, the relationship between star point grayscale and irradiance is:

$$[eqn]$$

where $[eqn]$ is the grayscale value of the $[eqn]$ -th star, $[eqn]$ is the maximum irradiance, $[eqn]$ is the magnitude of the $[eqn]$ -th star, and $[eqn]$ is the grayscale value corresponding to the minimum irradiance.

Stars are at vast distances from Earth and can be regarded as point sources for the observer. Due to inherent aberrations in the optical system, the grayscale of the star point diffuses, and the energy distribution of the image spot follows a Point Spread Function (PSF) model. By employing defocusing techniques to spread the star image over multiple pixels, the positioning accuracy of star extraction algorithms can reach sub-pixel levels (e.g., 1/10th of a pixel). To accurately simulate the star map, the grayscale dispersion caused by aberrations, defocusing, and glare must be considered. The energy distribution of the star image point is approximated by a two-dimensional Gaussian function model:

$$[eqn]$$

where $[eqn]$ is the center of the Gaussian surface, $[eqn]$ is the grayscale value at $[eqn]$ , $[eqn]$ is the total energy (grayscale sum), and $[eqn]$ represents the radius of the diffusion spot.

3.2.2. Target Energy Simulation

Target energy is determined via the Signal-to-Noise Ratio (SNR) and the noise background energy. The SNR of the imaging system is defined as the integrated SNR, which is the ratio of the total target energy to the mask energy. Thus, the total energy of the target image spot is:

$$[eqn]$$

Based on the analysis of the optical system’s Modulation Transfer Function (MTF)—assuming an off-axis three-mirror reflective configuration—the Point Spread Function (PSF) for target imaging can be derived as follows:

$$[eqn]$$

Consequently, the energy distribution of the target is expressed as:

$$[eqn]$$

3.3. Establishment of MTF Models

The MTF chain encompasses various stages, including the optical system, detector, and satellite platform motion. The image degradation caused by satellite platform jitter cannot be modeled in isolation; it must be integrated with the detector’s temporal integration effect [18,19,20]. Treating the system as a linear system, the total system MTF is expressed as:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ represent the MTFs of the optical system, detector, circuitry, and the integrated satellite-detector platform, respectively.

3.3.1. Optical System MTF

The optical MTF can be simulated using prior static transfer function data or by modeling effects such as diffraction, aberration, defocus, and fabrication/assembly errors:

1. Diffraction-limited MTF with Central Obscuration

For a circular aperture optical system with a central obscuration, the diffraction-limited MTF in the x-direction is given by

$$[eqn]$$

where $[eqn]$ is the obscuration ratio, defined as $[eqn]$ ; $[eqn]$ is the diameter of the optical system, and $[eqn]$ is the diameter of the obscured portion. Let:

$$[eqn]$$

where $[eqn]$ is the image spatial frequency, and $[eqn]$ is the optical cutoff frequency of the lens, defined as $[eqn]$ , in which $[eqn]$ represents the F-number of the system. The variables $[eqn]$ , $[eqn]$ , and $[eqn]$ are respectively expressed by the following equations:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

2.Diffraction-limited MTF without Obscuration

For an unobstructed circular aperture, the MTF is a circularly symmetric function:

$$[eqn]$$

3.Aberration MTF

For the aberration-related MTF, when $[eqn]$ :

$$[eqn]$$

where $[eqn]$ denotes the root-mean-square (RMS) wavefront error. The relationship between the RMS wavefront error and the peak-to-valley (P-V) wavefront error ( $[eqn]$ ) is given by $[eqn]$ , with the constant $[eqn]$ . This empirical MTF formula is applicable to cases with small wavefront errors, specifically where $[eqn]$ .

4.Defocus MTF

The modulation MTF caused by defocus is expressed as:

$$[eqn]$$

where $[eqn]$ is the first-order Bessel function; $[eqn]$ is the F-number of the optical system; $[eqn]$ is the defocus amount (unit: mm); $[eqn]$ is the Nyquist frequency. After precise focusing, the defocus MTF can be considered as 1.

5.Other MTF Factors

Assembly and fabrication also cause MTF degradation, which follows a Gaussian function:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ are the assembly and fabrication factors, respectively. The spatial effect caused by the optical system can be obtained by multiplying the various MTF factors:

$$[eqn]$$

3.3.2. Detector MTF

The effects influencing the detector MTF primarily include the spatial integration of pixel size and the photoelectron diffusion effect.

$$[eqn]$$

6.Space integral

During the detector imaging process, the image function is continuously intercepted by an aperture function. This local averaging process can be modeled as a convolution of the input signal with the sampling aperture function. The sampling aperture is typically regarded as a rectangular function, as detailed below:

$$[eqn]$$

The integral is solved using the rectangular rule as follows:

$$[eqn]$$

Taking stray light and irradiance non-uniformity into account, the corrected signal model is expressed as follows:

$$[eqn]$$

where $[eqn]$ is the radiance; $[eqn]$ is the average transmittance; $[eqn]$ is the F-number; $[eqn]$ is the area obscuration factor; $[eqn]$ is the vignetting coefficient; $[eqn]$ is the field of view (FOV) angle; $[eqn]$ is the image plane irradiance of stray light (related to the stray light coefficient $[eqn]$ , where $[eqn]$ ); $[eqn]$ is the detector pixel area; $[eqn]$ is the integration time; $[eqn]$ is the average quantum efficiency; $[eqn]$ is the number of integration stages (or TDI stages); $[eqn]$ is the speed of light; $[eqn]$ is Planck’s constant; $[eqn]$ is the center wavelength; $[eqn]$ is the offset; $[eqn]$ is the photoelectric conversion efficiency.

The expression for the vignetting coefficient t is given as follows:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the radii of the projected circles formed by the oblique beam passing through the entrance window and the entrance pupil, respectively, and $[eqn]$ is the distance between the two circle centers. To account for the non-linear response, the signal model is further refined as follows:

$$[eqn]$$

where the selection of the threshold values $[eqn]$ , $[eqn]$ , and $[eqn]$ depends on the type of detector material, the operating spectral band, and the manufacturing process level.

7.Simulation Methodology

This simulation model is based on stellar magnitude. It requires converting the apparent magnitude of stars into physical luminous flux, and subsequently calculating the signal electron count through imaging mechanism modeling.

First, the stellar brightness—representing the grayscale distribution of the star on the image plane—is derived based on the apparent magnitude formula (the specific transformation relationship between apparent magnitude and image grayscale is detailed in the subsequent section on stellar characteristics). The energy distribution on the image plane is then simulated by incorporating the optical aperture, transmittance, and Point Spread Function (PSF) provided by the camera manufacturer.

Secondly, parameters such as the quantum efficiency, pixel size, and exposure time of the detector are utilized to convert the number of photons reaching the detector into the corresponding number of electrons.

Finally, by processing the data with circuit gain, dynamic range, and noise characteristic models, the signal intensity for each pixel is output. This completes the simulation transformation from stellar apparent magnitude to a two-dimensional image.

4. Experiment

The simulation algorithms were implemented in Python 3.11, leveraging its open-source scientific computing libraries. All experiments were conducted on a workstation equipped with a 13th Gen Intel Core i7-13650HX processor (2.60 GHz), 16 GB of RAM, and an NVIDIA GeForce RTX 4060 GPU. Given the high computational complexity associated with processing the massive Gaia catalog, GPU acceleration (parallel computing) was specifically employed to handle the coordinate transformations of stars from the celestial sphere to the image plane. Additionally, a block-wise computation strategy was implemented to partition the data. This effectively prevents GPU out-of-memory (OOM) errors when calculating the image plane positions for exceptionally dense regions, such as the Galactic Center. This combination of hardware-accelerated parallelization and memory management significantly optimizes the runtime, enabling the efficient computation of stellar coordinates in highly dense star fields.

4.1. Fidelity Testing

Based on the theoretical foundation and modeling established above, dense stellar simulation experiments were conducted to test the imaging of stars in deep-space scenarios. To verify the accuracy and effectiveness of the simulation results, the measured data and corresponding optical axis pointing from the GTC camera were utilized, information related to the GTC camera is shown in Table 1 [21,22]. The imaging results were simulated under identical imaging conditions and simulation models—using stars up to magnitude 15 from the star catalog based on the detection capability of the measured data—and then compared with the actual captured images.

As shown in Figure 7, both the visual effect of individual stellar targets and the relative positions between stars in the simulated image are highly consistent with the GTC measured image. To further validate the model’s fidelity through a statistically unbiased evaluation, we implemented a stratified random sampling strategy for both stellar targets and background regions. For stellar reference selection, detected stars were categorized into several magnitude bins to ensure the evaluation covers the full dynamic range of the sensor. Within each bin, 100 stars were randomly sampled across different quadrants to account for potential field-of-view non-uniformity, with each star cross-referenced to its simulated counterpart via Gaia catalog coordinates. Simultaneously, 20 × 20 pixel background regions were localized using a Monte Carlo-based random sampling approach, while strictly avoiding high-brightness stars and target signals to prevent energy leakage from biasing the noise statistics. The mean grayscale values and standard deviations within these regions were then analyzed, and the results are illustrated in Figure 8 and Figure 9:

The fidelity of the simulation is quantitatively validated through four critical dimensions. To ensure statistical robustness, the reported metrics represent average values derived from multiple independent simulation runs. These runs encompass various sky regions—strategically selected based on different stellar densities—and diverse detection conditions (which are further detailed in Section 4.2.4). First, the average SNR error (<10%) confirms the accuracy of our radiative transfer model. Second, the centroid precision (<0.01 pixels), achieved via bicubic interpolation, ensures sub-pixel geometric consistency. Third, the background statistical errors (<5% and <10%) prove that the noise floor and stochastic fluctuations mimic real sensor behavior. Finally, the energy concentration (>90%) validates the high-fidelity reconstruction of the PSF. Collectively, these metrics demonstrate that the simulated images are not merely visually similar, but physically equivalent to measured data for algorithm performance evaluation.

4.2. Simulation Experiments for Special Scenarios

4.2.1. Simulation of Dim and Weak Targets

For the simulation of dim and weak targets, a directional Gaussian filter is employed to simulate the smearing (trailing) effect of moving targets within the image. The brightness of the target region is then enhanced based on the Signal-to-Noise Ratio (SNR) and local statistical characteristics. Finally, an integrated image containing the moving targets is synthesized. Let $[eqn]$ and $[eqn]$ be the velocity components. The formulas for calculating the smearing length and the smearing angle are as follows:

$$[eqn]$$ $$[eqn]$$

Let $[eqn]$ be the standard deviation of the target region, $[eqn]$ be the average background grayscale, and $[eqn]$ be the binary mask. The target brightness enhancement formula is given as follows and the simulation results are shown in Figure 10:

$$[eqn]$$

4.2.2. Platform Motion Simulation

The impact of platform motion is primarily manifested in the image motion (smearing) of stars and targets. During simulation, it is necessary to evaluate the magnitude of background stellar movement, considering both stellar position displacement and image motion blur. Consequently, the coordinates of the stars and targets on the image plane must be updated frame by frame. The position update formulas are as follows:

$$[eqn]$$

where $[eqn]$ represents the image plane position of the star or target at the $[eqn]$ frame, $[eqn]$ is the direction of image motion, and $[eqn]$ is the image motion distance per frame. The image motion effect is illustrated in the Figure 11.

4.2.3. Star–Target Overlap (Sticking) Simulation

When the stellar density is high, overlapping (sticking) between the target and background stars is inevitable, necessitating the use of mask processing. The masking strategy is as follows: compare the grayscale values of the star and the target within the overlapping region; if the target’s energy in the overlap is higher than that of the star, the original region is replaced by the target’s pixels; otherwise, the original stellar pixels are retained. The logic for updating the image region is expressed as follows and the simulation results are shown in Figure 12 and Figure 13:

$$[eqn]$$

4.2.4. Simulation of Different Detection Limits

Corresponding simulations were conducted to compare the performance under various detection limits of the detector. We performed comparative simulations for two scenarios: magnitude 15 at SNR 3 and magnitude 22 at SNR 3. It can be observed that under the 15-magnitude @ SNR 3 condition, the image spatial duty cycle is lower, and the overall image brightness is relatively dim. In contrast, under the 22-magnitude @ SNR 3 condition, the image spatial duty cycle is significantly higher, resulting in a higher overall brightness; furthermore, multiple stars may overlap within a single pixel, the simulation results are shown in Figure 14.

5. Conclusions

To address the critical requirement for high-fidelity data in high-orbit Space Situational Awareness (SSA) and deep-space target recognition, this paper presents a comprehensive simulation methodology for dense stellar backgrounds. By integrating the Gaia catalog with target Two-Line Element (TLE) data, we developed an end-to-end optoelectronic imaging model capable of accurately simulating faint signals, such as space debris and asteroids, against complex deep-space environments. Quantitative evaluations demonstrate that the proposed method achieves high fidelity across multiple metrics, including Signal-to-Noise Ratio (SNR), geometric precision, and background statistical characteristics, with all relative errors maintained below 10%.

Recognizing the computational challenges posed by the vast volume of the Gaia dataset—particularly when simulating high limiting magnitudes or large image scales—a strategic data optimization framework was implemented. The Gaia catalog was spatially partitioned into manageable sub-sets, with data storage streamlined to retain only essential simulation attributes: Right Ascension (RA), Declination (Dec), and apparent magnitude. Furthermore, the catalog was categorized into discrete subsets based on limiting magnitudes (e.g., 10, 14, 18, and 22), facilitating on-demand data loading for varied mission requirements. This hierarchical optimization significantly enhances data retrieval efficiency and ensures the scalability of the framework for generating high-density star maps.

Simulation results highlight the significant impact of signal aliasing and background interference on dim target detection within high-orbit environments. By covering a broad spectrum of scenarios—ranging from magnitude 15 to 22 with varying target dynamics—this system provides a robust platform for the design optimization of future space-based detection payloads, algorithm training, and limit-capability assessments.