© The Authors 2024

1 Introduction

Since the launch of the Iridium satellite constellation comprising 66 satellites in 1997 (James 1998), and up to the current date of 26 October 2023, with Starlink having deployed over 4900 satellites, the impact of these satellite constellations on astronomical observations has gained increasing attention from astronomers and astronomy enthusiasts (Barentine et al. 2023; Falchi et al. 2023). Proposals for additional satellite constellations similar to the Starlink project, such as OneWeb and Amazon’s Project Kuiper, are on the rise. Projections suggest that by 2030, the number of low-Earth-orbit satellites dedicated to communication will surpass 100000 (Lawrence et al. 2022).

Satellites produce bright trails on images captured by ground-based and space telescopes due to sunlight reflection, leading to contamination of imaging and spectral data (Kruk et al. 2023; Mróz et al. 2022; Kovalev et al. 2023); the diffuse reflection of space debris resulting from satellite collisions contributes to an increase in the brightness of the night sky, reducing imaging contrast (Kocifaj et al. 2021; Barentine et al. 2023); communication satellites also introduce significant radio frequency interference for radio telescopes (Di Vruno et al. 2023; Grigg et al. 2023); and short-term occultations of observed targets by satellites can lead to inaccurate photometry (Groot 2022).

The initial Starlink V1.0 satellites had a visual magnitude of 5.1, which decreased to 7.1 after the application of low-reflectance coatings. The final Starlink VisorSats have a visual magnitude of 6.0 (Halferty et al. 2022). Furthermore, the median visual magnitude of OneWeb satellites is recorded at 8.23 (Krantz et al. 2023). Additionally, in their pursuit of enhanced communication signals, SpaceMobile designed the BlueWalker 3 satellite with an expansive size of 64 square meters, resulting in a visual magnitude of 0.2 (Nandakumar et al. 2023). Moreover, high-orbit satellites can reach a visual magnitude of 12 or even brighter (Zhao et al. 2021; Kerr et al. 2021). It should also be noted that because telescopes require compensation for the Earth’s rotation during observation, geostationary satellites will leave trails in telescope imaging. Therefore, most satellites will indeed have a significant impact on astronomical observations.

In the face of the escalating severity of satellite interference, effective measures have been proposed to reduce its impact. Osborn et al. (2022) achieved satellite transit forecasting using two-line element (TLE) data to assist in observational planning and avoidance of satellites. Hu et al. (2022) introduced a method that does not rely on precise satellite orbit prediction. Their method involves constructing a satellite distribution map and devising observational strategies to avoid regions with a high satellite density.

Meanwhile, various models have been proposed for quantifying the impact of satellites on astronomical observations. McDowell (2020) introduced distinct categories of satellites in orbit and assumed an evenly spaced satellite distribution along the orbital paths. These authors used this assumption to develop a transient satellite distribution model and analysed the number of illuminated Starlink satellites at varying latitudes, times of the year, and times of the night. In a related work, Hainaut & Williams (2020) introduced the concept of a uniform satellite distribution model, wherein satellites are uniformly dispersed across the Earth’s surface. These authors investigated how the count of observable satellites is affected by solar latitude and estimated the brightness of satellites. Building upon this work, Bassa et al. (2022) refined the satellite uniform distribution model, considering circular orbits with a certain inclination angle for the satellites, enabling a more precise estimation of both the probability distribution of satellites and their brightness.

To analyse the impact of all satellites on observations, it is necessary to consider not only low-Earth-orbit satellites but also high-Earth-orbit satellites and satellites on highly elliptical orbits. A satellite distribution probability model based on elliptical orbits is required for this purpose. In our research, we optimised the circular orbit model by considering satellite elliptical orbits. We constructed a satellite probability distribution model based on these elliptical orbits. Using real satellite distribution data, we provided the probability of observing satellites from various ground stations at different viewing angles. Additionally, we conducted a reliability check of our approach by comparing it to observation results. We used imaging data obtained from LAMOST’s acquisition cameras (Cui et al. 2012; Hu et al. 2023) to analyse the probability of satellite observations from the Xinglong Observatory, validating the effectiveness of our method.

This paper is organised as follows: The first section provides an overview of the data sources we used in our research, including satellite orbit data and data from the LAMOST guide cameras. The second section delves into our satellite probability distribution model, which is based on elliptical orbits. We explain the methodology and equations used to construct this model. The third section focuses on our analysis of the impact of satellites on observations at the Xinglong Observatory, specifically using data from the LAMOST acquisition cameras. We compare the observed results with the predictions from our satellite probability distribution model. Finally, we summarise our conclusions based on the discussion.

2 Data

In this section, we provide the data pertinent to the development and validation of the probability model for satellite distribution. The data encompass information about satellite orbital and the LAMOST acquisition camera.

2.1 Satellite orbital data

A satellite’s orbit can be determined by six parameters, commonly including orbital inclination, right ascension of the ascending node (RAAN), orbital eccentricity, argument of perigee, perigee altitude, and true anomaly, as illustrated in Fig. 1. Satellite orbit data are stored as TLE files with three lines, each containing 69 characters. From these files, we extract the necessary data, which include the orbital inclination (i), argument of perigee (aop), and orbital eccentricity (e). Additionally, the number of orbits a satellite completes in a day (n) is provided, enabling the calculation of the satellite’s orbital period (T), as well as the semi-major axis (a) and semi-minor axis (b) of the orbit, where $$T=\frac{24\times 60\times 60}{n}\text{, }$$(1) $$a={\left(\frac{GM{T}^2}{4{\pi }^2}\right)}^{1/3},$$(2)

and $$b=a\times \sqrt{1-e^2}\text{. }$$(3)

We downloaded orbit data for all in-orbit satellites from the space-track website¹ on 1 September 2023, amounting to a total of 25 606 objects. This dataset includes geostationary satellites, medium-Earth-orbit satellites, low-Earth-orbit satellites, and high-elliptical-orbit satellites, among others. After removing satellite debris, often designated with the abbreviation ‘deb’ in the name from the TLE data, there were 14214 remaining objects. Among these, 9381 were identified as low-Earth-orbit satellites, 201 as high-elliptical-orbit satellites, and 1023 as geostationary satellites. Table 1 provides the number of satellites for prominent satellite mega-constellations.

Fig. 1 Parameters required to determine the orbit of a satellite. These include the orbital inclination, which is the angle between the satellite’s orbital plane and the Earth’s equatorial plane. Additionally, the satellite’s RAAN denotes the point where the satellite crosses from the southern hemisphere to the northern hemisphere while rotating anticlockwise with respect to right ascension. The argument of perigee signifies the angle between the line connecting the Earth’s centre to the nearest point on the orbit (perigee) and the line connecting the Earth’s centre to the ascending node. Other essential parameters include the eccentricity and semi-major axis of the orbit.

2.2 LAMOST acquisition camera data

LAMOST is a fibre-fed spectroscopic telescope consisting of optical paths Ma and Mb, along with a focal surface equipped with 4000 optical fibres. It boasts a five-degree field of view for wide-field sky surveys, making it particularly susceptible to the impact of satellites. The focal surface accommodates eight acquisition cameras, positioned at 1.5 degree and 2.5 degree field positions, recording the tracks of satellites. Using satellite trails in the imaging data, we can estimate their influence on the sky survey telescope. Each camera comprises 1024 × 1024 pixels, with each pixel measuring 24 micrometres, corresponding to a field of view of 0.025 arcseconds. Figure 2 presents a sample image influenced by satellites as captured by the acquisition cameras.

Fig. 2 Satellite tracks recorded on the focal surface of the acquisition cameras of LAMOST. The images in the bottom-left and bottom-right corners document the trajectories of two satellites obtained during observations in the same celestial region.

3 Method

In this section, we describe the extension of the method proposed by Bassa et al. (2022) to elliptical orbits, enabling ELT data, including data for satellites on highly eccentric orbits. We constructed a geographic coordinate system with the Earth’s centre as the origin, and the satellite’s coordinates can be expressed as (lat,lon,R). In this representation, (lat, lon) denote the latitude and longitude of the intersection point between the satellite-to-Earth centre line and the Earth’s surface, while R represents the distance from the satellite to the Earth’s centre. The orbital trajectory of a satellite can be expressed by the elliptical equation: $$R=\frac{a(1-e^2)}{1+e\cos (\theta -aop)},$$(4)

where θ represents the angle between the line connecting the satellite to the centre of the Earth and the line connecting the ascending node to the Earth’s centre. The variable aop stands for the argument of perigee. We can express θ using the satellite’s latitude ϕ in the celestial coordinate system, as illustrated in Fig. 3. In the spherical triangle O – ABC, the angle θ of the ellipse is related to the latitude in the geographic coordinate system, as $$\frac{\sin \theta }{\sin {90}^{\circ }}=\frac{\sin \varphi }{\sin i},$$(5)

and therefore the orbit’s altitude can be represented by the satellite’s latitude. The probability dp of a satellite being within a range of dθ angles can be calculated as the ratio of the area S _dθ subtended by dθ degrees between the elliptical orbit and the total area S ellipse of the ellipse, as shown in Eq. (6), assuming the satellite fully complies with Kepler’s Second Law. $$\mathrm{d}p=\frac{S_{\mathrm{d}\theta }}{S_{\text{ellipse }}}=\frac{R^2\mathrm{d}\theta }{2\pi ab}.$$(6)

Taking the differential of Eq. (5), we obtain $\mathrm{d}\theta =\frac{\cos \varphi }{\cos \theta \sin i}d\varphi$. Substituting this into Eq. (6) results in $$\frac{\mathrm{d}p}{\mathrm{d}\varphi }=\frac{R^2}{2\pi ab}\frac{\cos \varphi }{\cos \theta \sin i}.$$(7)

Taking into account that the same ϕ corresponds to two different θ, denoted as θ₁ and θ₂, with corresponding satellite altitudes R₁ and R₂, the above equation should be modified to: $$\frac{\mathrm{d}p}{\mathrm{d}\varphi }=\frac{R_1{}^2}{2\pi ab}\frac{\cos \varphi }{\cos {\theta }_1\sin i}+\frac{R_2{}^2}{2\pi ab}\frac{\cos \varphi }{\cos {\theta }_2\sin i}.$$(8)

Considering that, in general, within a given orbital plane, there will be multiple satellites with the same orbital inclination, and accounting for the Earth’s rotation, the right ascension of the ascending node for satellites will continuously change. Consequently, the relative positions of satellites will also vary. As a simplifying assumption, we consider that the probability of a satellite being at a given latitude is the same for different longitudes. This means that the probability of a satellite appearing at a certain position (lat, lon, R) on its orbital surface is independent of longitude (lon) and depends only on latitude (lat) and the distance from the Earth’s centre (R). Thus, we can obtain the probability of a satellite in a unit solid angle as: $$\begin{array}{rl}{p}_{\text{sat }}& =\frac{\mathrm{d}p}{\mathrm{d}\mathrm{\Omega }}=\frac{\mathrm{d}p}{\mathrm{d}\varphi \times 2\pi \times \cos \varphi }\\ & =\frac{R_1{}^2}{4{\pi }^{2}ab}\frac{1}{\cos {\theta }_1\sin i}+\frac{R_2{}^2}{4{\pi }^{2}ab}\frac{1}{\cos {\theta }_2\sin i}.\end{array}$$(9)

According to the formulas presented above, we can calculate the probability distribution of satellites corresponding to various longitudes and latitudes.

Fig. 3 Projection of a satellite orbit on the celestial sphere. In the geographic coordinate system, the projection of the ascending node on the celestial sphere is denoted by point A, and the projection of the satellite’s position is point C. In this system, i denotes the orbital inclination, and ϕ represents the satellite’s latitude.

4 Result and discussion

In this section, we present the key conclusions of our method and validate the model using observational data.

Fig. 4 Observable satellite count at different latitudes, considering all satellites, with the equatorial region being most severely affected. The upper figure depicts the count of observable satellites within a 30 degree field of view centred on the zenith at the observatory, while the lower figure represents the count of observable satellites within a 60 degree field of view centred on the zenith at the observatory.

4.1 Satellite probability distribution

Equation (9) addresses the probability distribution for an individual satellite, taking into account a total of N satellites in orbit. The overall observation probability signifies the count of observable satellites, necessitating the summation of the observation probabilities for each of the N satellites. Taking into account the existence of N satellites, for a telescope with a field of view A and an exposure time t_exp, L represents the width of the telescope’s field of view. Here, ω_sat,j denotes the angular velocity of the satellite and (A + Lω_sat,j × t_exp) represents the total field of view traversed by the telescope during the exposure time. Therefore, the number of satellites observable in a single exposure is given by the following formula: $$N_{\text{sat}}={\displaystyle \sum _{\text{sat}=1}^N}{\displaystyle \sum _{\mathrm{j}=1}^2}{\rho }_{\text{sat},\mathrm{j}}\times (A+L{\omega }_{\text{sat},\mathrm{j}}\times t_{\text{exp}}).$$(10)

The symbol ρ_sat denotes the satellite number density at the telescope’s location, which can be calculated through p_sat, where $$\begin{array}{rl}{\rho }_{\text{sat}}& =\sum _{\mathrm{j}=1}^2{\rho }_{\text{sat},\mathrm{j}}=\sum _{\mathrm{j}=1}^2{p}_{\text{sat},\mathrm{j}}\times \frac{d_j^2\times \cos {\alpha }_{2,j}}{R_j^2\times \cos {\alpha }_{1,j}},\\ & =\sum _{\mathrm{j}=1}^2\frac{R_j{}^2}{4{\pi }^{2}ab}\times \frac{1}{\cos {\theta }_j\sin i}\times \frac{d_j^2\times \cos {\alpha }_{2,j}}{R_j^2\times \cos {\alpha }_{1,j}},\end{array}$$(11) $${\omega }_{\mathrm{s}\mathrm{a}\mathrm{t},\mathrm{j}}=\frac{2\pi ab}{T{R}_j^2},$$(12)

where α₂ represents the angle between the satellite-to-Earth centre vector and the normal vector to the satellite orbital surface, the term α₁ represents the angle between the satellite-to-observation point vector and the normal vector to the satellite orbital surface, and d is the distance between the satellite and the observation point. We present a detailed calculation of cos α₁, cos α₂, and d in Appendix A, and α₁, α₂, and d are shown in Fig A.1.

Fig. 5 Observable satellite count considering only low-Earth-orbit satellites. The image reflects the results considering only low-Earth-orbit satellites, with the impact being most severe around latitudes of ±50 degrees and ±80 degrees, consistent with the orbital inclinations of Starlink and OneWeb satellites. The upper figure depicts the count of observable satellites within a 30 degree field of view centred on the zenith at the observatory, while the lower figure represents the count of observable satellites within a 60 degree field of view centred on the zenith at the observatory.

4.2 Satellite possibility map

The ρ_sat can serve as an indicator of the impact of satellites on the observation site. Considering the intricate nature of a satellite’s orbital shell, which typically consists of two rotationally symmetric surfaces, we calculate the heights R₁ and R₂ of these two orbital planes at the zenith position of the observation site. In two observation modes, where the observed elevation angle is greater than 30 degrees or greater than 60 degrees, the orbital surface can be approximated as two spherical surfaces with radii R₁ and R₂. Under this approximation, cos α₂ = 1, and $${\rho }_{\text{sat }}={\displaystyle \sum _{j=1}^2}{p}_{\text{sat,}\mathrm{j}}\times \frac{d_j^2}{R_j^2\times \cos {\alpha }_{1,j}}.$$(13)

We integrate ρ_sat within 30 degrees and 60 degrees near the zenith of the observation site. We performed integration using TLE data for 14 214 satellites, including all satellites. The results are summed up, as shown in Fig. 4.

The results indicate that the maximum number of observable satellites is near the equator, potentially influenced by the presence of numerous geostationary satellites. Moreover, the number of observable satellites at elevation angles surpassing 30 degrees surpasses those at angles exceeding 60 degrees.

Additionally, we also computed the results considering only low-orbit satellites. We integrated the data for 9314 low-orbit satellites, reflecting the future impact of satellites on telescope observations. The results are shown in Fig. 5. The results indicate that future satellite impact will be particularly severe in regions around ±50 degrees and ±80 degrees. These correspond to the typical orbital inclinations of Starlink and OneWeb satellites. Due to the higher orbit of OneWeb satellites, more observable satellites are present in ±80 regions. Comparing the integration results for low-orbit satellites with those including all satellites, it is evident that the influence of high-orbit satellites still dominates the overall impact.

4.3 Observation results

We conducted our research using data obtained from the acquisition camera of the LAMOST telescope from January 2023 to June 2023. Initially, we excluded images that could not be used to determine the presence of satellite trails, such as those taken during the transition between sky regions and images with camera malfunctions. The remaining dataset consists of 226 630 images from eight acquisition cameras, with a total exposure time of 7 295 597 s. This corresponds to a telescope observation time of 911 949.625 s. This non-integer value results from the exclusion of certain camera images that were unsuitable for determining the presence of satellite trails.

We manually labelled these images, identifying 423 images containing satellite trails, for a total of 427 trails. Satellite trails are straight lines with light intensity evenly distributed along their trajectory, making them easy to distinguish from phenomena such as stellar diffraction spikes, CCD smear, malfunctions in the guiding-tracking system, and asteroid trails. Stellar diffraction spikes show uneven light intensity, becoming lighter and wider closer to the star. Malfunctions in the guiding-tracking system will cause all stellar images on the CCD to elongate. Asteroid trails are shorter and thinner. CCD smear, always connected to an overexposed star, aligns with the readout direction of the CCD. The distribution of observation time for telescopes and the distribution of the 427 satellite trails are shown in the top-left plot of Fig. 6.

The number of satellite trails recorded by the acquisition cameras is only a subset of all the satellites passing through the focal surface. Using the count of satellite tracks captured by the acquisition cameras, we can estimate the total number of satellite trails acquired across the entire focal surface. This issue can be reformulated as the probability that any straight line traversing the focal surface does not intersect with the eight fixed-position small circles on the focal surface. This probability can be expressed as $$P=\frac{{\int }_0^{2\pi }{\int }_0^{\pi }p(\theta ,\beta )\mathrm{d}\theta \mathrm{d}\beta }{2{\pi }^2},$$(14)

where θ and β are as shown in Fig. 7, p(θ,β) = 1 when the line does not intersect the small circle, and p(θ,β) = 0 when the line crosses the small circle. We conducted a differential simulation of θ and β, and we calculated the distance between the straight line and the midpoint of each camera; we then compared the distance with the size of the camera to determine whether the line intersects with the camera. The actual imaging units of the camera are square, and we created inscribed and circumscribed circles for the squares, taking the radius of the circle as the size of the camera. The results indicate a probability range of 0.8333 to 0.8797 for the likelihood of not being imaged by any of the eight cameras and passing through the focal surface.

According to Eq. (10), the total number of observed satellites during n observations can be expressed as: $$\begin{array}{rl}& N_{n-\mathrm{o}\mathrm{b}\mathrm{s}}=\sum _{\mathrm{o}\mathrm{b}\mathrm{s}=1}^n{\left(N_{\mathrm{s}\mathrm{a}\mathrm{t}}\right)}_{\mathrm{o}\mathrm{b}\mathrm{s}}\\ & =\sum _{\text{obs }=1}^n{[\sum _{\text{sat }=1}^N\sum _{j=1}^2{\rho }_{\text{sat,}j}\times (A+L{\omega }_{\text{sat,}j}{t}_{\text{exp }})]}_{\text{obs }}\\ & =n\times (\overline{\sum _{\text{sat }=1}^N\sum _{j=1}^2{\rho }_{\text{sat},j}}\times A+\overline{\sum _{\text{sat }=1}^N\sum _{j=1}^2{\omega }_{\text{sat},j}}\times L{t}_{\text{exp }})\\ & =n\times (\overline{\rho }\times A+\overline{\rho \omega }\times L{t}_{\exp })\text{. }\\ & \end{array}$$(15)

The averages $\overline{\rho }$ and $\overline{\rho \omega }$ represent the mean values over n observations and can be expressed as the average of the integrated results over the observed sky regions. Considering that the latitude of the LAMOST telescope is 40.3959 degrees north, the average elevation angle during observations from January to June 2023 is 51.11 degrees, with a median of 50.01 degrees. We can then calculate the $\overline{\rho }$ and $\overline{\rho \omega }$ of four observation scales and satellite types, and the results are shown in Table 2.

We also considered whether satellites could be illuminated during the observation process. It is necessary to consider whether the line connecting the satellite and the sun is blocked by the Earth, and whether the angle formed by the Sun–satellite–Earth is greater than 90 degrees. If the distance from the Earth’s centre to the satellite along the line to the Sun is greater than the Earth’s radius, or if the angle formed by the Sun–satellite–Earth is greater than 90 degrees, the satellite is illuminated by the Sun. We considered scenarios with satellite altitudes of 1000 km and 10000 km in different observation modes. Through simulation, we calculated the proportion of satellites illuminated at different times within various observation ranges at the Xinglong Observatory (altitude 900 m, longitude 117.58 degrees, latitude 40.40 degrees), where the LAMOST telescope is located; the results are shown in the top-right plot of Fig. 6.

Using the integrated results of the satellite observation probability in Eq. (15), we calculated the number of observable satellites for four observation types based on the distribution of observation time of LAMOST and the number of times the telescope changed observation targets (the number of times the observation sky area was changed), as shown in the bottom-left plot of Fig. 6. We then considered whether the Sun can illuminate the satellites to determine the observable number of satellites. For the low-Earth-orbit satellites, we used illuminated probability data calculated for 1000 km, while for the other satellites we used illuminated probability data calculated for 10 000 km. The results are compared with the actual observational statistics, as shown in the bottom-right of Fig. 6.

Taking into account that the mean zenith angle for LAMOST observations exceeds 50 degrees, the observed results should fall between observations with zenith angles of greater than 30 degrees and those with zenith angles of greater than 60 degrees. The results indicate good consistency between the model and observational outcomes near midnight. However, the observations during dusk and dawn show a higher discrepancy. This discrepancy may be attributed to the possibility that fragments of low-Earth-orbit satellites are sufficiently bright during these times, forming satellite tracks.

5 Conclusions

We propose an evaluation method for telescopes based on TLE data, considering the impact of not only low-Earth-orbit satellites but also high-Earth-orbit and high-eccentricity satellites. The assessment, conducted using satellite data from September 2023, reveals that, when considering all satellites, the regions near the equator are most affected due to the influence of geostationary satellites. In scenarios considering only low-Earth-orbit satellites, the impact is most severe around latitudes of ±50 degrees and ±80 degrees. This is primarily attributed to the influence of Starlink and OneWeb satellites.

To validate the effectiveness of our method, we used observational data from the acquisition cameras of the LAMOST telescope for the period from January to June 2023. We conducted a statistical analysis of satellite trails captured by the guiding cameras during this period. Comparing the observational results with those considering all satellites, our findings reveal a better fit around midnight. However, there is a slight overestimation in the actual observations during the early morning and late evening. This discrepancy may be attributed to the possibility that fragments of low-Earth-orbit satellites are sufficiently bright during these times, forming satellite tracks. This comparison supports the validity of our method in predicting satellite interference in telescope observations.

Our method does not account for the impact of the telescope’s observation limit, assuming that all satellites at different altitudes are visible. The effective aperture of LAMOST used for validation is 4.4 m, and for telescopes with smaller apertures, the impact of high-orbit satellites would be less pronounced. Additionally, our method allows integration over smaller regions, enabling the analysis of satellite distribution probabilities in specific observational areas, which can be useful for satellite observation scheduling.

Fig. 6 Intermediate results of data processing and the comparison of our model’s results with observed results. The top-left subplot illustrates the distribution of observation time of LAMOST, which amounts to 991 949.625 seconds, as well as the distribution of the 427 satellite trails. The bottom-left subplot illustrates the time distribution of the number of observable satellites under different observation modes. This analysis, based on the results from Table 2, incorporates the distribution of observation time of the telescope and the number of times the telescope changed observation targets (the number of times the observation sky area was changed). However, it does not consider whether satellites can be illuminated by the Sun. The top-right plot depicts the time distribution of satellite visibility on March 15, 2023 for satellites at altitudes of 1000 km and 10 000 km. It considers satellites within the range of elevation angles greater than 30 degrees and those within the range of elevation angles greater than 60 degrees. The bottom-right plot compares the predictions of the model, which takes into account whether satellites can be illuminated by the Sun, with the observed results.

Fig. 7 Probability of satellites passing through the focal surface without being captured by the cameras. The left figure illustrates the configuration of LAMOST’s focal surface acquisition cameras. In the right figure, the probability that the lines corresponding to different numbers do not intersect with any of the eight acquisition cameras is depicted.

Acknowledgements

Thanks are given to the reviewer for the constructive comments and helpful suggestions. Dr. Sun Tianrui from Purple Mountain Observatory has provided valuable suggestions and discussions regarding the visibility of high-orbit satellites. This work is supported by the National Nature Science Foundation of China (Grant No 12203079, 12103072 and 12073047), the Natural Science Foundation of Jiangsu Province (Grants No BK20221156 and BK20210988), and the Jiangsu Funding Program for Excellent Postdoctoral Talent. The Guo Shou Jing Telescope (or LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project was provided by the National Development and Reform Commission. LAMOST is operated and managed by National Astronomical Observatories, Chinese Academy of Sciences.

Appendix A Impact angle of observation

In the geographic coordinate system, the satellite orbital surface is formed by the rotation of an elliptical satellite orbit around the axis of the Earth’s rotation. The orbital surface equation can be expressed in terms of longitude θ and latitude ϕ as follows: $$r(\theta ,\varphi )=\frac{a\ast (1-e^2)}{1+e\ast cos(arcsin(sin\varphi /sini)-aop}.$$(A.1)

Equation 9 provides the angular density distribution of satellites based on the focal point (Earth’s centre) of the satellite orbit. To calculate the observable number of satellites during telescope observations, it needs to be transformed into a satellite probability density distribution centred around the observation point.

According to Formula A.1, the normal vector to the orbital plane can be obtained as $$\nabla r={\mathit{e}}_{\mathit{r}}+\frac{1}{r}\frac{\partial r}{\partial \phi }{\mathit{e}}_{\mathit{\phi },}$$(A.2)

where $$\frac{\mathrm{\partial }r}{\mathrm{\partial }\varphi }=\frac{-a(1-e^2)\ast e\ast \frac{-sin\varphi }{sini}\ast \frac{1}{\sqrt{1-(sin\varphi /sini{)}^2}}\ast \frac{cos\varphi }{sini}}{1+ecos(acrsin(sin\varphi )/sini-aop{)}^2}.$$(A.3)

Therefore, the normal vector to the orbital plane represented in Euclidean space is $$\mathit{s}\mathit{h}\mathit{e}\mathit{l}\mathit{l}=(cos\frac{1}{r}\frac{\mathrm{\partial }r}{\mathrm{\partial }\varphi },0,Rsin\frac{1}{r}\frac{\mathrm{\partial }r}{\mathrm{\partial }\varphi }).$$(A.4)

The coordinates of the observation point and the satellite are given by (0, ϕ_obs, R_obs) and (θ_sat, ϕ_sat, R_sat).

The vector from the observation point to the satellite is represented as $$\begin{array}{rl}& \mathit{o}\mathit{b}\mathit{s}=(Rcos{\varphi }_{\text{s}at}cos{\theta }_{\text{s}at}-rcos{\varphi }_{\text{o}bs},\\ & Rcos{\varphi }_{\text{s}at}sin{\theta }_{\text{s}at}-rcos{\varphi }_{\text{o}bs},Rsin{\varphi }_{\text{s}at}-rsin{\varphi }_{\text{s}at}).\end{array}$$(A.5)

The vector from the Earth’s centre to the satellite is $$\mathit{e}\mathit{a}\mathit{r}\mathit{t}\mathit{h}=(Rcos{\varphi }_{\text{s}at}cos{\theta }_{\text{s}at},Rcos{\varphi }_{\text{s}at}sin{\theta }_{\text{s}at},Rsin{\varphi }_{\text{s}at}).$$(A.6)

According to formulae A.4 and A.5, we can obtain $$\cos {\alpha }_1=\frac{\mathit{o}\mathit{b}\mathit{s}\bullet \mathit{s}\mathit{h}\mathit{e}\mathit{l}\mathit{l}}{\parallel \mathit{o}\mathit{b}\mathit{s}\parallel \ast \parallel \mathit{s}\mathit{h}\mathit{e}\mathit{l}\mathit{l}\parallel },$$(A.7)

and $$\begin{array}{rl}{d}^2& ={(R\times sin{\varphi }_{\text{s}at}-r\times sin{\varphi }_0)}^2\\ & +{(R\times cos{\varphi }_{\text{s}at}\times sin{\theta }_{\text{s}at}-r\times cos{\varphi }_0\times sin{\theta }_0)}^2.\\ & +{(R\times cos{\varphi }_{\text{s}at}\times cos{\theta }_{\text{s}at}-r\times cos{\varphi }_0\times cos{\theta }_0)}^2\end{array}$$(A.8)

According to formulae A.4 and A.6, we can obtain $$\cos {\alpha }_2=\frac{\mathit{e}\mathit{a}\mathit{r}\mathit{t}\mathit{h}\bullet \mathit{s}\mathit{h}\mathit{e}\mathit{l}\mathit{l}}{\parallel \mathit{e}\mathit{a}\mathit{r}\mathit{t}\mathit{h}\parallel \ast \parallel \mathit{s}\mathit{h}\mathit{e}\mathit{l}\mathit{l}\parallel }\mathbf{,}$$(A.9)

and $$\begin{array}{rl}{R}_{\text{s}at}^2& ={(R\times sin{\varphi }_{\text{s}at}-r\times sin{\varphi }_0)}^2\\ & +{(R\times cos{\varphi }_{\text{s}at}\times sin{\theta }_{\text{s}at}-r\times cos{\varphi }_0\times sin{\theta }_0)}^2.\\ & +{(R\times cos{\varphi }_{\text{s}at}\times cos{\theta }_{\text{s}at}-r\times cos{\varphi }_0\times cos{\theta }_0)}^2\end{array}$$(A.10)

Fig. A.1 The positional relationship between the satellite, the observation point, and the Earth’s centre. The α2 represents the angle between the satellite-to-Earth centre vector and the normal vector to the satellite orbital surface. The term α1 represents the angle between the satellite-to-observation point vector and the normal vector to the satellite orbital surface, and d is the distance between the satellite and the observation point.

References

All Tables

All Figures

Fig. 1 Parameters required to determine the orbit of a satellite. These include the orbital inclination, which is the angle between the satellite’s orbital plane and the Earth’s equatorial plane. Additionally, the satellite’s RAAN denotes the point where the satellite crosses from the southern hemisphere to the northern hemisphere while rotating anticlockwise with respect to right ascension. The argument of perigee signifies the angle between the line connecting the Earth’s centre to the nearest point on the orbit (perigee) and the line connecting the Earth’s centre to the ascending node. Other essential parameters include the eccentricity and semi-major axis of the orbit.

In the text

	Fig. 2 Satellite tracks recorded on the focal surface of the acquisition cameras of LAMOST. The images in the bottom-left and bottom-right corners document the trajectories of two satellites obtained during observations in the same celestial region.
In the text

	Fig. 3 Projection of a satellite orbit on the celestial sphere. In the geographic coordinate system, the projection of the ascending node on the celestial sphere is denoted by point A, and the projection of the satellite’s position is point C. In this system, i denotes the orbital inclination, and ϕ represents the satellite’s latitude.
In the text

Fig. 4 Observable satellite count at different latitudes, considering all satellites, with the equatorial region being most severely affected. The upper figure depicts the count of observable satellites within a 30 degree field of view centred on the zenith at the observatory, while the lower figure represents the count of observable satellites within a 60 degree field of view centred on the zenith at the observatory.

In the text

Fig. 5 Observable satellite count considering only low-Earth-orbit satellites. The image reflects the results considering only low-Earth-orbit satellites, with the impact being most severe around latitudes of ±50 degrees and ±80 degrees, consistent with the orbital inclinations of Starlink and OneWeb satellites. The upper figure depicts the count of observable satellites within a 30 degree field of view centred on the zenith at the observatory, while the lower figure represents the count of observable satellites within a 60 degree field of view centred on the zenith at the observatory.

In the text

Fig. 6 Intermediate results of data processing and the comparison of our model’s results with observed results. The top-left subplot illustrates the distribution of observation time of LAMOST, which amounts to 991 949.625 seconds, as well as the distribution of the 427 satellite trails. The bottom-left subplot illustrates the time distribution of the number of observable satellites under different observation modes. This analysis, based on the results from Table 2, incorporates the distribution of observation time of the telescope and the number of times the telescope changed observation targets (the number of times the observation sky area was changed). However, it does not consider whether satellites can be illuminated by the Sun. The top-right plot depicts the time distribution of satellite visibility on March 15, 2023 for satellites at altitudes of 1000 km and 10 000 km. It considers satellites within the range of elevation angles greater than 30 degrees and those within the range of elevation angles greater than 60 degrees. The bottom-right plot compares the predictions of the model, which takes into account whether satellites can be illuminated by the Sun, with the observed results.

In the text

	Fig. 7 Probability of satellites passing through the focal surface without being captured by the cameras. The left figure illustrates the configuration of LAMOST’s focal surface acquisition cameras. In the right figure, the probability that the lines corresponding to different numbers do not intersect with any of the eight acquisition cameras is depicted.
In the text

Fig. A.1 The positional relationship between the satellite, the observation point, and the Earth’s centre. The α2 represents the angle between the satellite-to-Earth centre vector and the normal vector to the satellite orbital surface. The term α1 represents the angle between the satellite-to-observation point vector and the normal vector to the satellite orbital surface, and d is the distance between the satellite and the observation point.

In the text