https://orcid.org/0000-0002-8824-9956
Abstract
The determination of the full population of resident space objects (RSOs) in low Earth orbit (LEO) is a key issue in the field of space situational awareness that will only increase in importance in the coming years. We endeavour to describe a novel method of inferring the population of RSOs as a function of orbital height and inclination for a range of magnitudes. The method described uses observations of an orbit of known height and inclination to detect RSOs on neighbouring orbits. These neighbouring orbit targets move slowly relative to our tracked orbit, and are thus detectable down to faint magnitudes. We conduct simulations to show that, by observing multiple passes of a known orbit, we can infer the population of RSOs within a defined region of orbital parameter space. Observing a range of orbits from different orbital sites will allow for the inference of a population of LEO RSOs as a function of their orbital parameters and object magnitude.
1 INTRODUCTION
The number of resident space objects (RSOs) in Earth orbit is increasing year on year (Mark & Kamath 2019; Blake 2022; Bernhard, Deschamps & Zaccour 2023). There are data bases which attempt to track the population of these objects (e.g. space-track.org) but the fraction which are fully characterized falls off dramatically with size. In low Earth orbit (LEO) particularly, targets smaller than |${\sim }\, 10$| cm have low completeness (Pardini & Anselmo 2021). These targets are hard to track since most are not placed in orbit by design, but are debris associated with specific launches or satellites. Optical observations of these objects rely on reflected sunlight, thus their magnitude is proportional to their reflecting surface area and albedo. This means the smallest objects can be below the threshold of traditional observational methods.
The increase of objects in LEO is leading towards an increasing risk of in-orbit collisions (Pardini & Anselmo 2020; Acciarini et al. 2023). Already, in recent years there have been multiple accidental collisions (as well as a few anti-satellite demonstrations), each of which has increased the budget of LEO debris by many thousands of pieces (Springer et al. 2010; Pardini & Anselmo 2011). Due to their high relative velocities, LEO object collisions can produce catastrophic events (Olivieri & Francesconi 2020), both for the colliding objects, and by potentially increasing the density of debris to such an extent that entire sections of Earth orbit become unusable, known as the Kessler Syndrome (Kessler & Cour-Palais 1978). Avoiding future collisions requires a well populated data base of RSOs, detailing their physical parameters and, importantly, their orbital elements (Blake et al. 2023). A knowledge of the debris population will afford satellite operators the chance to perform orbital manoeuvres to avoid collisions (Slater, Byram & Williams 2006; Uriot et al. 2022). Accurate RSO orbital elements will improve the effectiveness and efficiency of these manoeuvres (Dutta & Misra 2023).
Because of the small size and high orbital velocity of many LEO objects, they can be challenging to observe optically. The most straightforward way to observe a population of objects is to use a stationary telescope and detect targets as they cross the field of view (e.g. Chote, Blake & Pollacco 2019; Blake et al. 2021). However, even for targets of sufficient size that the flux would be detectable, the relative velocity between the target and a stationary telescope may be sufficiently high that the flux is spread over so many pixels that each one falls below the detectability limit of the telescope. There have been multiple attempts to resolve these issues (Suthakar et al. 2023), including use of machine learning for detection (e.g. Wright, Rocchetto & Waldmann 2023), observing from space (e.g. Clark et al. 2022), and exposure stacking; both targeted (Zhu et al. 2023) and blind (Cooke et al. 2023). This paper describes a new approach, utilizing neighbouring orbits. To avoid the issue of high relative velocities we employ a tracking telescope set-up, targeting a well defined orbit. We then attempt to recover RSOs on similar orbits, which will result in much reduced relative velocities. The method is thus capable of fainter (and thus smaller) target detection and is more focused on particular orbital regimes.
The neighbouring orbits detection technique as discussed herein is broadly simulation-based. We focus on a simulation-based approach to reduce the dependence on significant amounts of observational data which can be both expensive and limited. As this paper is a discussion of the technique and its applications, and not a full survey, this is reasonable. When discussing target detection in detail, we use example data from a specific telescope set-up to demonstrate the capabilities and limitations of the technique. The use of observational data and an expanded survey will be the basis of future work.
This paper begins with a description of neighbouring orbits and orbit generation (Sections 2 and 3), before detailing the simulation in full (Sections 4 and 5). It then continues into a discussion of the method used to detect targets on these neighbouring orbits (Sections 6 and 7) and how to convert detections into an inferred population (Section 8). Finally, we present our simulation results and conclusions, highlighting some future work (Sections 9 and 10).
To simplify, Sections 2 and 3 outline our assumptions and definitions. The simulation is described in Section 4 and the detection process described in Section 6. Section 8 then extrapolates detections into results which are described in Sections 9 and 10. Section 5 builds upon the base simulation, expanding the approach to different orbits and observing locations. Section 7 further details our detection method, justifying our methodology and providing avenues for adjusting thresholds.
2 ORBITAL PARAMETERS
To simplify our approach to this problem, we make the assumptions that all orbits are circular (i.e. eccentricity, |$e=0$|). This means that we can define an orbit, and an object’s position on it, using four key parameters. These four parameters are height h, inclination i, a modified right ascension of the ascending node |$\Omega$|, and a modified true anomaly |$\nu$| (modifications are based on the chosen observing site and are discussed in Section 3 and Appendix A). For computational ease, as well as simplifying the definition of other orbital parameters, we define a height h and then set the apoapsis and periapsis equal to |$h+R_\oplus \pm \epsilon$| where |$\epsilon$| is a negligibly small value so as to simulate a circular orbit without forcing |$e=0$|. Since the orbit is, in practice, circular, |$\nu$| does not actually affect the orbit and is instead used to assign a given object to a particular part of its orbit, allowing us to compare targets on identical orbits but with different epochs. Fig. 1 displays the effect of changing the parameters h, i, |$\Omega$|, and |$\nu$| on a satellite’s orbit. We define an orbit O as a unique combination of the four orbital parameters which is written as |$O = (h, i, \Omega , \nu )$| with h given in units of kilometres and |$i, \Omega , \nu$| in units of degrees.
Impact of orbital parameter variations. Lines show orbits offset in the relevant parameter. For (d), points show the difference in orbit start point, since the orbit itself is unchanged by the |$\nu$| parameter. Lines progress from red (inner) to blue (outer); h ranges from |$0.2-0.4R_\oplus$|, i ranges from 80° to 60°, |$\Omega$| ranges from 10° to 30°, |$\nu$| ranges from 10° to 30°. Figures are generated using the Mayavi software package (Ramachandran & Varoquaux 2011).
3 TLE GENERATION
To use orbital parameters to predict the position of an RSO at any given time, we first generate a two line element set (TLE). A TLE is simply a data format which encodes the specific parameters of an orbit. Specialized codes can then be used to predict the coordinates of an object on said orbit at any point in time. To generate the required TLEs we use the faketle tool from SATTOOLS (Bassa 2023) which requires the four orbital parameters described above, as well as the epoch. We wish the generated orbit to have height h and inclination i and to pass directly overhead (i.e. across the zenith) of our observing site (with latitude = |$\delta$| and longitude = |$\alpha$|) at a given epoch. We can use these parameters (h, i, |$\alpha$|, |$\delta$|, epoch) as well at the zenith crossing requirement, to determine the final two parameters |$\Omega$| and |$\nu$|. |$\Omega$| and |$\nu$| are given by the following equations:
where
and |$A = \sin {\alpha }\cos {\delta }\cos {\theta }$|, |$B = \cos {\alpha }\cos {\delta }\cos {\theta }$|, and |$C = \sin {\theta }\sin {\delta }$|. For a full definition and derivation, see Appendix A.
The five orbital parameters are then used to generate the corresponding TLE which, in turn, is used to determine the coordinates of an object on said orbit at any point in time. The conversion of TLE to coordinates is done using Skyfield (Rhodes 2019). The orbit is designed to cross the zenith of our observing site at the chosen epoch, but the amount of time for which it is observable before and after depends on both the observing site, and the specific orbital parameters. For a target to be observable it must satisfy two conditions. First, it must be at least 20° above the horizon of the observing site. This is due to a combination of instrumental and airmass restrictions. Secondly, it must be sunlit. For a target to be observable it must reflect sunlight which can then be seen by the telescope, therefore the vector from the Sun to the target must not be intercepted by the Earth, while the Sun must also be at least 6° below the horizon as seen from the observing site. Fig. 2 shows these criteria for an orbit with |$h=850$| km and |$i=99^{\circ }$|, observed from our site on La Palma.
Observing constraints for a target on an orbit with |$h=850$| km and |$i=99^{\circ }$| observed from a site on La Palma.
From Fig. 2(a), we see that our target is above the required altitude limits for 460 s per pass, symmetrically around the chosen epoch. Fig. 2(b) shows how the fraction of this period which is also sunlit changes throughout a night, beginning at fully observable after sunset, falling to fully unobservable during the middle of the night, and then rising again before sunrise. Successive passes are separated by 460 s, the length of time for which a target is above |$20^{\circ }$| altitude. The number of fully observable passes per night is 25. Successive passes refer to successive observations of theoretical targets with the same orbital elements (h and i) but different ascending nodes and true anomalies (|$\Omega$| and |$\nu$|), not passes of the same target (which would be separated by |$\sim 90$| min at LEO). Each pass tracks the path of a different theoretical target, which cross the meridian in 460 s intervals (|$\Omega$| and |$\nu$| are adjusted to maintain a zenith-crossing object at an offset epoch).
4 OFFSET TLE SIMULATION
To determine whether an object in a neighbouring orbit is detectable, we need to generate a second TLE and a corresponding set of RA/Dec coordinates, at the same time stamps as used for the principal orbit. A neighbouring orbit refers to an orbit whose four orbital parameters are similar to those used in the principal orbit and are defined in terms of offsets from the principal orbit values. For each of the four orbital parameters we define a maximum offset value (|$h_{\rm off}$|, |$i_{\rm off}$|, |$\Omega _{\rm off}$|, and |$\nu _{\rm off}$|) and an offset step size (|$h_{\rm step}$|, |$i_{\rm step}$|, |$\Omega _{\rm step}$|, and |$\nu _{\rm step}$|). The maximum offset value is found using an iterative method and is such that increasing the offset of that orbital parameter any further (positively or negatively) will not result in any further detectable targets. It is dependent on the height and inclination of the orbit, and the latitude of the site, as discussed in Section 5. The offset steps are constant and chosen to be 2 km in h and |$0.1^{\circ }$| in i, |$\Omega$|, and |$\nu$|. We test both positive and negative offsets, thus the number of tested offsets for each parameter is given by |$p_n = 2(p_{\rm off}/p_{\rm step}) + 1$| (where p is a placeholder for any one of the four orbital parameters and the |$+1$| term accounts for an offset of zero). We test all combinations of the four parameter offsets, resulting in |$N_{\rm neighbours}$| neighbouring orbits being tested, where |$N_{\rm neighbours} = h_ni_n\Omega _n\nu _n$|.
For each combination of offset orbital parameters we then produce a TLE as described in Section 3 and generate RA/Dec coordinates for each of the time steps used for the principal orbit. We then convert these RA/Dec coordinates into frame positions, assuming that the frame is centred on the RA/Dec coordinates of the principal orbit. The equations are reproduced here, and described in full in Govind et al. (2018), |$(\alpha _1, \delta _1)$| are the RA/Dec coordinates of the principal orbit, and |$(\alpha _2, \delta _2)$| are the RA/Dec coordinates of the neighbouring orbit (all coordinates are given in radians).
where |$x_{\rm sf}$| and |$y_{\rm sf}$| are scale factors to convert radians to pixels given by
and |$(x_c, y_c) = \left(\frac{9600}{2}, \frac{6422}{2}\right)$| are the coordinates of the centre of the frame (frame size used is for the CLASP instrument, specifications detailed in Table 1). From the x/y coordinates of the neighbouring orbit we can determine its relative separation from the principal tracked orbit and, assuming linear motion between successive time stamps, its relative velocity. Fig. 3 shows a plot of the relative motion of two targets on neighbouring orbits. The principal orbit is |$O_{\rm principal} = (850, 99, 29.8, 28.9)$| and the two targets have orbits of |$O_1 = O_{\rm principal} + (+2, +0.1, +0.1, -0.1)$| (right hand side, green/yellow) and |$O_2 = O_{\rm principal} + (-2, +0.1, -0.1, +0.1)$| (left hand side, red/yellow). We see here how only relatively minor changes in the orbital parameters can drastically affect the resulting motion across the FoV.
Location of targets on neighbouring orbits across the FoV of an instrument tracking the principal orbit for a full pass. The central, red point denotes the centre of the FoV (tracking |$O_{\rm principal} = (850, 99, 29.8, 28.9)$|) and the red box, its edges. The blue axes lines describe the central axes of the frame. The tracks of the two neighbouring targets are shown coloured by their velocities relative to the centre of the frame. Target 1, |$O_1 = O_{\rm principal} + (+2, +0.1, +0.1, -0.1)$|, is on the right, coloured green/yellow and target 2, |$O_2 = O_{\rm principal} + (-2, +0.1, -0.1, +0.1)$| is on the left, coloured red/yellow.
Instrument specifications.
| Parameter | CLASP |
|---|---|
| Diameter | 36 cm |
| FoV | |$2.63^{\circ } \times 1.76^{\circ }$| |
| Frame size | |$9600 \times 6422$| pixels |
| Focal ratio | F/2.2 |
| Detector | QHY600M sCMOS |
| Pixel size | 3.76 μm |
| Gain | 0.42 e|$^-$|/ADU |
| FWHM | 3.6 pixels |
| Zero point | 23.01 mag |
| Parameter | CLASP |
|---|---|
| Diameter | 36 cm |
| FoV | |$2.63^{\circ } \times 1.76^{\circ }$| |
| Frame size | |$9600 \times 6422$| pixels |
| Focal ratio | F/2.2 |
| Detector | QHY600M sCMOS |
| Pixel size | 3.76 μm |
| Gain | 0.42 e|$^-$|/ADU |
| FWHM | 3.6 pixels |
| Zero point | 23.01 mag |
Instrument specifications.
| Parameter | CLASP |
|---|---|
| Diameter | 36 cm |
| FoV | |$2.63^{\circ } \times 1.76^{\circ }$| |
| Frame size | |$9600 \times 6422$| pixels |
| Focal ratio | F/2.2 |
| Detector | QHY600M sCMOS |
| Pixel size | 3.76 μm |
| Gain | 0.42 e|$^-$|/ADU |
| FWHM | 3.6 pixels |
| Zero point | 23.01 mag |
| Parameter | CLASP |
|---|---|
| Diameter | 36 cm |
| FoV | |$2.63^{\circ } \times 1.76^{\circ }$| |
| Frame size | |$9600 \times 6422$| pixels |
| Focal ratio | F/2.2 |
| Detector | QHY600M sCMOS |
| Pixel size | 3.76 μm |
| Gain | 0.42 e|$^-$|/ADU |
| FWHM | 3.6 pixels |
| Zero point | 23.01 mag |
For a target on a neighbouring orbit to be detectable it must satisfy two criteria. First, its x/y location must place it within the bounds of a frame centred on a target on the principal orbit. That is, we require |$0 \le x \lt 9600$| and |$0 \le y \lt 6422$|. Additionally, the relative velocity of the target must be below a certain threshold. Targets moving too fast, relative to the tracked target, will have their flux spread over too many pixels to be detectable using the detection method described in Section 6. The exact velocity threshold depends on the instrumental specifications and the detection goals and is discussed further in Sections 6 and 7, but the default value used here is equal to 10 pix s−1. We also require that these thresholds are met for a minimum number of consecutive frames. Based on the detection method described in Section 6, we require a target to be within the FoV and below the limiting relative velocity value for at least 20 consecutive frames. Additionally, we make the assumption that targets on neighbouring orbits are sunlit at the same time as a target on the principal orbit. This is a result of the relative speed of the TLE being significantly larger than the rotation speed of the Earth on the small FoV scale being considered. Simulations have confirmed that this is a reasonable assumption, which can significantly reduce computational runtime for testing orbital offset combinations.
Each combination of offset parameters therefore produces a binary value of detectability. We do not account for targets which are detectable for a larger fraction of a pass than others, and simplify our results to detectable versus non-detectable offset parameter combinations. Fig. 4 shows the resulting detectability map for a principal orbit with |$O_{\rm principal} = (850, 99, 29.8, 28.9)$| observed from La Palma using CLASP. The results are presented in a 2D histogram in |$h-i$| space with the colour representing the number of |$\Omega -\nu$| combinations that produce a detectable target for each combination of |$h-i$|. For these parameters we find that detectable orbital parameters are as follows: |$h = 850^{+156}_{-126}$| km, |$i = 99.0^{+3.9}_{-4.0}\, ^{\circ }$|, |$\Omega = 29.8^{+3.1}_{-3.2}\, ^{\circ }$|, |$\nu = 28.9^{+1.7}_{-1.4}\, ^{\circ }$|. Note, however, that not every combination of offset parameters within these ranges will be detectable, just that these are the extreme values for each parameter independently.
Target detectability as a function of |$h_{\rm offset}$| and |$i_{\rm offset}$| for |$O_{\rm principal} = (850, 99, 29.8, 28.9)$| observed from La Palma (|$\delta =29^{\circ }$|) using CLASP. Colour designates relative probability, i.e. number of |$\Omega -\nu$| combinations resulting in a detection for each |$h-i$| pair. Values assume a limiting relative velocity of 10 pix s−1.
The bimodality of target detectability seen in Fig. 4 is due to the interplay between offsets in the four orbital elements. Near the very centre of the figure, the offsets in h and i are small, giving rise to an orbit very similar to the central one. Therefore, it takes only a small offset in |$\Omega$| or |$\nu$| to generate an orbit which is very similar in motion but sufficiently offset on the sky as to be outside the FoV and thus, undetectable. A slightly larger offset in either h or i gives a more significantly different orbit, creating a track which changes in position and velocity relative to the central one. Since the relative position and motion changes during the visible window, there are more combinations of |$\Omega$| and |$\nu$| offsets which produce a target orbit that is both sufficiently close and moving sufficiently slowly relative to the central one as to be detectable for at least a portion of this observing window. Increasing the h and i offsets further creates an orbit so different from the central one that no combinations of |$\Omega$| and |$\nu$| offsets can result in a target that is detectable.
The approximate symmetry is due to the fact that, on small scales, it matters not whether the compared orbit is offset positively or negatively to the central one. The relative difference between the central and offset orbit is broadly unchanged. The result is a broadly symmetrical bimodal distribution of detectability, peaking at low, but non-zero, h and i offsets, separated by a depressed region in the centre of the |$h-i$| offset parameter space and surrounded by a region of low detectability at the extreme edges.
5 EFFECT OF HEIGHT, INCLINATION, AND LATITUDE
The detectable ranges for |$h_{\rm offset}$|, |$i_{\rm offset}$|, |$\Omega _{\rm offset}$|, and |$\nu _{\rm offset}$| are dependent on both the chosen principal values of h and i, and the latitude of the observing site. The first impact is that the chosen values of h and i affect the amount of time for which a target is visible per orbit (we do not consider the effect of changing |$\Omega$| or |$\nu$| as they are not chosen directly, but calculated from h and i). Larger values of h mean higher orbits, and thus, slower moving targets, since the period of an orbit, |$P \propto a^{3/2} = (R_\oplus + h)^{3/2}$|. This means that a target remains above 20° for a greater amount of time per orbit. The number of fully observable passes per night remains roughly constant, since higher orbits mean that more time is required per pass, but also that a target remains sunlit for a greater amount of time after sunset and before sunrise. The combination of these two effects effectively cancel out. For example, observing an orbit with |$h=750$| km allows for 410 s per pass and 26 passes per night. An orbit with |$h=950$| km allows 500 s per pass and 25 passes per night. As described in Section 3, each pass tracks the path of different theoretical target (with the same h and i values), which cross the meridian in 410 s (at |$h=750$| km) or 500 s (at |$h=950$| km) intervals.
An increase in latitude of the observing site has no impact on the maximum observable time per pass (under the assumption that a pass crosses through the zenith), but it does affect how many fully visible passes are possible per night. The length of twilight depends on Sun elevation and is longer closer to the poles, allowing for more passes per night. At a latitude of 29° we can observe 25 fully visible passes per night. At a latitude of 50° this increases to 33 passes per night and at a latitude of 75° we can observe 76 passes per night. It should be noted that these simulations were carried out during the Northern hemisphere winter months at a solar declination angle of −21° (close to mid-winter in the Northern hemisphere, when the length of twilight is at its shortest). During the summer months twilight lasts longer and thus more observable passes per night would be expected.
Alongside the effects on the available observing time per pass, the chosen values of h, i, and |$\delta$| can change the range of detectable orbital parameter offsets. Some of this is coupled to the change in observable time per pass. For example, consider a combination of orbital parameter offsets that produces a target that is only detectable for 10 s, from |$Z-230$| to |$Z-220$| (where Z is the time at which the tracked orbit reaches zenith). If the chosen h, i, |$\delta$| values allow for only 410 s of observation per pass (i.e. 205 s before and after zenith) this offset orbit will not be considered detectable since the observable period begins 15 s after the detectable period has ended. If, however, the chosen values allow for 500 s of observation time per pass (i.e. 250 s before and after zenith) the offset orbit will be considered detectable since the detectable portion occurs within the observable period. There are also visible time independent effects however. As discussed above, latitude has no effect on observable time per pass, but does affect the range of offset parameter combinations which produce detectable targets. Fig. 5 shows the same sort of plot as displayed in Fig. 4 but with different latitudes. The chosen orbit is |$O_{\rm principal} = (850, 99, 29.8, 28.9)$| with |$\delta = 50^{\circ }$| (Fig. 5a), and |$\delta = 75^{\circ }$| (Fig. 5b) as opposed to |$\delta = 29^{\circ }$| used in Fig. 4. For ease of comparison, Appendix B reproduces Figs 4, 5(a) and (b) on the same axis scale.
Target detectability as a function of |$h_{\rm offset}$| and |$i_{\rm offset}$| for |$O_{\rm principal} = (850, 99, 29.8, 28.9)$|. Colour designates relative probability, i.e. number of |$\Omega -\nu$| combinations resulting in a detection for each |$h-i$| pair. Values assume a limiting relative velocity of 10 pix s−1. For ease of comparison, Appendix B reproduces Figs 4, 5(a) and (b) on the same axis scale (Fig. 12).
As we can see from this figure, the chosen value of |$\delta$| can significantly affect the resulting detectable ranges for |$h_{\rm off}$|, |$i_{\rm off}$|, |$\Omega _{\rm off}$|, and |$\nu _{\rm off}$|. We find the detectable ranges to be
As |$\delta$| increases, the detectable ranges of the four orbital parameters decreases. However, we see that although there are fewer allowed |$h-i$| combinations, the density of detectable combinations actually increases. For example, we see that for |$\delta =29^{\circ }$| the |$h-i$| combination with the highest density of |$\Omega -\nu$| combinations has 24 valid |$\Omega -\nu$| combinations. For |$\delta =50^{\circ }$|, this increases to a maximum of 39 valid |$\Omega -\nu$| combinations and for |$\delta =75^{\circ }$| we find a maximum of 97 valid |$\Omega -\nu$| combinations. Fig. 6 shows the relative motion of a target as observed from three different latitudes. The principal orbit is |$O_{\rm principal} = (850, 99, 29.8, 28.9)$| and the neighbouring targets have |$O = O_{\rm principal} + (+2, +0.1, +0.1, -0.1)$|. The relative motions are calculated as observed from |$\delta =29^{\circ }$|, |$50^{\circ }$|, and |$75^{\circ }$|.
Relative motion of targets on neighbouring orbits as observed from different latitudes. General form is similar to Fig. 3. The principal orbit is |$O_{\rm principal} = (850, 99, 29.8, 28.9)$| and offset orbit is |$O = O_{\rm principal} + (+2, +0.1, +0.1, -0.1)$|. The three tested latitudes are |$\delta =29^{\circ }$| (inner data points, shown in green/yellow), |$\delta =50^{\circ }$| (middle data points, shown in red/yellow), and |$\delta =75^{\circ }$| (outer data points, shown in blue/green).
The cause of the difference in detectable offset ranges is a combination of two major effects. The first is that, as |$\delta$| increases, the relative offset between two orbits with different inclinations increases. Two orbits, differing in only inclination, will overlap where the orbit crosses the equator and have maximum separation furthest from the equator [see Fig. 1(b) for a representation of this]. Therefore, when observed from higher latitudes, an offset in inclination will have an increased effect of relative separation and motion between two targets, leading to a reduced inclination offset range for detectable targets. The second effect is that the relative offset between orbits with different |$\Omega$| values decreases as the observing latitude increases. Orbits that differ only in |$\Omega$| will have their maximum separation at the equator and minimum separation closest to the poles (see Fig. 1c for a representation of this). Therefore, observing from a higher latitude site will result in a reduced effect from |$\Omega$| offsets leading to a higher density of |$\Omega -\nu$| combinations being detectable. These effects are complicated due to the fact that an offset in any single parameter affects the detectable offsets in every other parameter, but these factors are the driving forces behind the latitude-dependent effects seen in Figs 5 and 6. The optimal choice of observing site depends on the exact goals of the survey and is discussed further in Section 10.
6 DETECTION
Targets are identified in our observational data through the following method. Individual frames are first reduced using a combination of bias, dark, and flat fields. Master bias and dark frames are obtained once, near the start of the observing period, and are compiled from an average of multiple individual frames. A master flat field frame is obtained every evening that observations occur, again, being an average of multiple individual frames. Observational data are gathered using the CLASP instrument on La Palma (see Table 1 for instrument specifications).
To identify slow moving targets, i.e. those objects on similar orbits to our tracked orbit, and for the purposes of this project, we employ a frame stacking method. Frames are arranged sequentially and, after data reduction of each frame, we median stack rolling batches of n frames where n is an independent variable, with the default value being 10 (see Section 7 for an in depth discussion of the optimal choice of variables). Each subsequent batch contains the last |$n-1$| frames from the previous batch along with the next frame in the observation sequence. These stacked frames are then searched using the Python library for Source Extraction and Photometry (sep; Bertin & Arnouts 1996; Barbary 2016), and we limit ourselves to smaller signals that have limited extension. Contamination from stellar signals are less of an issue, since they are mainly removed due to the stacking procedure and the limit of signal extension.
For each detection, we record the detection parameters, as well as a list of the pixel coordinates of each pixel determined to be part of the detections. The results of SEP can be used to infer the magnitude and relative velocity of each detection. We then move into an additional step, required to limit false positives, and increase confidence in our detections. True detections are likely to appear in more than one frame and move across these frames in a predictable pattern, whereas false positives are more likely to appear in few, or single, frames and may appear to move in a more random pattern. To this end, we implement a bespoke clustering algorithm. For each detection, we define a physical and temporal buffer zone. We then search this buffer zone for other, similar, detections. If one is found, we assume it belongs to the same object as the original detection, and add it to the cluster of detections. The process is then repeated using this new detection as the origin of the ‘buffer’ zone. This continues until no new detections are found, and we are left with a cluster of detections, similar in both physical and temporal space. Clusters which contain a suitable number of individual detections are determined to be real, and clusters which do not reach this threshold, are rejected as false positives. Fig. 7 shows an example injected signal through each step in the process. Panel (a) shows a single image with the injected signal, while panel (b) shows a median stack of |$n=10$| consecutive frames. Panel (c) then shows the sum of all stacked frames in which the signal is detected and panel (d) shows the calculated cluster, with each cluster pixel coloured by the frame number in which it was identified.
Detection of an injected signal. Panel (a) single 0.5 s exposure with injected signal. Panel (b) median stack of |$n=10$| frames. Panel (c) sum of all stacked frames in which signal is detected. Panel (d) calculated cluster, with each cluster pixel coloured by the frame number in which it was identified.
7 INJECTION/RECOVERY TESTS
It is important to know the theoretical faintest object that can be detected using the frame stacking method discussed above. The limiting magnitude is a function of the object’s relative orbit (i.e. its relative speed with respect to the tracked orbit) and a range of observational and processing parameters including the number of threshold crossing pixels required for a detection (|$n_{\rm pix}$|), the number of frames per stack (n), and the exposure time of the frames (t). An additional approach that could affect detectability is binning pixels before attempting detection. This technique is not explored here but a discussion of pixel binning and its effects on the recovery of streaked targets is given in Cooke et al. (2023). Also important are a range of telescope-specific parameters such as the noise properties of individual frames (|$\sigma$|), as well as the Point Spread Function/Full-width Half Maximum (PSF/FWHM) and Zero Point (ZP). We use the CLASP telescope on La Palma for this work. The relevant Instrument specifications are given in Table 1.
We have bench-marked this detection method through a range of injection/recovery tests. We simulate a range of signals, each corresponding to a target of a given magnitude, moving at a given relative velocity across the frame [streaks are simulated following the procedure laid out in Cooke et al. (2023)]. We then employ the above detection method to determine the reliability with which it can be recovered. Targets are injected multiple times, on different parts of the frame. The probability for recovery is simply the fraction of injections for a given magnitude/velocity combination that are recovered. Fig. 8 shows the recovery probability for a tested set of injected signals spanning a range of magnitudes and velocities. The red line corresponds to a 50 per cent recovery fraction contour, assuming a linear interpolation between distinct points and the blue line shows a polynomial fit. These results employ |$n = 10$|, |$n_{\rm pix} = 50$|, and |$t=0.5$|s.
Recovery probability for injected signals with a range of magnitude and velocity values. The red, piecewise line indicates a contour of 50 per cent recovery probability, assuming a linear interpolation between distinct points and the smooth, blue line shows a polynomial fit to the same data.
From Fig. 8 we see that the limiting magnitude for 50 per cent recoverability falls with increasing target velocity. This is to be expected, since the faster a target is moving, the more pixels its flux is spread over, and thus the fainter the signal it produces. Additionally, faster moving objects move across more pixels between subsequent exposures, making it less likely that stacking frames will be able to align the brightest parts of a signal, reducing the brightness of the signal in the stacked image. We also see that, somewhat counterintuitively, the limiting magnitude also falls off for the slowest moving targets. This is due to the fact that the signals in subsequent images are so well aligned that the resulting signal in the stacked image is too small (in terms of number of pixels) to be identified by our detection method (to minimize the number of false positives we define a minimum number of pixels required for a detection, see Fig. 9b).
Limiting magnitude as a function of relative velocity of an object. Lines are polynomial fits to data. Sub-figures show how the relation is affected by different values of n, |$n_{\rm pix}$|, and t. The highlighted tracks corresponds to |$n = 10,\ n_{\rm pix} = 50,\ t = 0.5$|s.
The best-fitting line for 50 per cent recovery, calculating magnitude M as a function of velocity v, is a cubic, and has the form
with |$a=0.006515$|, |$b=-0.1445$|, |$c=0.5864$|, and |$d=15.57$|. Inverting the above equation to allow the determination of the maximum velocity for a given magnitude (and selecting for the correct cubic root) gives
with |$\displaystyle p=\frac{3ac-b^2}{3a^2}$| and |$\displaystyle q=\frac{2b^3-9abc+27a^2(d-M)}{27a^3}$|. For derivation of cubic roots, see Cox (2012).
To determine how this limiting magnitude varies with the observational and processing parameters described above, we repeat the injection and recovery tests with a range of parameter values. Fig. 9 shows how the limiting magnitude is affected by each of these variables in turn. For each sub-figure, two variables are held constant at their default values while the third is varied. Default values are defined as |$n_{\rm pix} = 50,\ n = 10,\ t = 0.5$|s. In each sub-figure the highlighted line (thicker line, outlined in pink) shows the limiting magnitude for the default values.
From Fig. 9(a), we see that lowering the value of n leads to the detection of fainter targets at the high velocity end of the distribution. This is due to the fact that targets need to remain on the same pixels for fewer frames, thus they can be correspondingly fainter and still be detected. However, at the slower velocity end, the detection limit actually decreases. This is due to the fact that the flux from these targets is now spread over so few pixels, that their extended size on a frame is below the threshold size for detection. Based on simulations (see Section 4), we see that the relative velocity of an object changes during the observing window. Therefore, targets which are, at one point, moving too slowly to be detected, are likely to increase their relative speed during observations, moving them into a region of parameter space with a higher recovery probability, resulting in a successful recovery. Additionally, lowering the number of frames leads to a less successful removal of background stellar signals and thus an increase in the number of false positives. Increasing the value of n has the opposite effect, with a brighter detection limit for faster moving targets and a slightly fainter limit for the slowest moving targets. Using |$n=10$| gives the best trade off between magnitude limit and false positives.
From Fig. 9(b), we see that lowering the value of |$n_{\rm pix}$| leads to a fainter magnitude limit at all velocities. This is as expected, since allowing smaller signals to be picked up means that fainter targets, which have a smaller size on the frame, can be detected. The consequence is that reducing this threshold leads to a drastic increase in the number of false positive detections. Conversely, increasing the value of |$n_{\rm pix}$| means that targets must be brighter to be successfully detected, but reduces the number of false positives. We find that |$n_{\rm pix}=50$| gives a good compromise between increased detectability and a reduced number of false positives.
From Fig. 9(c), we see that lowering the value of t results in improved detections for the faster moving targets. This is due to shorter exposures meaning that targets move across fewer pixels in individual frames (trails are shorter). This increases the chance that they remain on similar parts of the image in successive frames, thus they are more likely to be detectable in stacked images. At the slower end of the velocity range the detectability falls off with shorter exposures, due to the reduced size of the signal on an image, moving it below our detection threshold. Reducing the exposure time also leads to a higher data rate. In a realistic scenario, data storage and processing time are important considerations, both of which increase as exposure time decreases. For longer exposure times, the data requirements are reduced but the detectability suffers due to the inverse of the effects mentioned above. Using |$t=0.5$|s gives a compromise between these considerations.
8 DETECTIONS TO NUMBERS
To convert detections into a prediction for the number of objects within a certain region of orbital parameter space, we need a reasonable underlying distribution. One approach is to use a uniform distribution in the four orbital parameters described above, h, i, |$\Omega$|, and |$\nu$|. This however, is not realistic for all parameters, since some orbits are more favoured than others, for example, polar orbits (with |$i\sim 90^{\circ }$|) are significantly more heavily populated than others. Additionally, depending on the purpose of particular space missions, some orbital heights are likewise more heavily utilized. To determine a reasonable underlying distribution, we utilized the ESA Meteoroid and Space Debris Terrestrial Environment Reference model v8.0.3 (MASTER; Horstmann et al. 2016). MASTER uses a combination of confirmed objects with TLEs and simulations to produce a theoretical population of orbiting debris (MASTER also includes a population of meteoroids). For |$\Omega$| and |$\nu$| we assume a uniform distribution but for h and i we use the MASTER data base to infer the underlying distribution of satellites and debris in Earth orbit. We use the MASTER population for all objects larger than 1 cm in size, which is a combination of known objects and simulation. We assume that the underlying distribution of MASTER objects is a reasonable approximation to the true underlying distribution, even if the actual number of objects is less certain. Utilizing the MASTER distribution will allow us to make our own inferences about the number of targets in a particular region of parameter space.
On large scales we find significant structure in the distribution of h and i of RSOs. Fig. 10 shows a portion of a 2D histogram of the MASTER population |$h-i$| parameter space, coloured by the relative number of objects in each bin.
Relative number of objects as function of orbital height and inclination (full height and inclination ranges are truncated for visibility) from the MASTER population. Only objects |$\ge 1$| cm in size are shown.
In this figure we see that some |$h-i$| bins are much more heavily populated than others. On smaller scales however, the distribution is much less structured. Within each bin shown in Fig. 10 the distribution of specific objects is much more random. Therefore, we assume that, though the number of objects per bin can be inferred from the MASTER population, the distribution of objects within each bin is uniform. The bin sizes we use as a separation between these two regimes are 25 km in h and 0.5° in i.
To use this distribution to convert detections to targets we proceed in the following way. We assume N objects total, within the tested region of |$h-i$| parameter space. The extent of the |$h-i$| region about which inferences can be made is a result of the minimum and maximum detectable values of |$h_{\rm offset}$| and |$i_{\rm offset}$| found by the simulations discussed in Section 4.
This region of parameter space is then split into A offset parameter bins, each corresponding to a different tested combination of |$h_{\rm offset}$| and |$i_{\rm offset}$|. For each of these offset parameter bins, we determine which MASTER population bin it falls into (the MASTER population bins are shown in Fig. 10), noting that since the offset bins are smaller there will be |$b = \frac{25}{2}\times \frac{0.5}{0.1} = 62.5$| offset bins within each MASTER population bin.
The fraction of the N total objects that fall into each offset bin is given by |$M_a/bM$| where |$M_a$| is the number of MASTER population objects in the MASTER bin that contains the current offset bin a, and M is the total number of MASTER objects in all A offset bins within the tested |$h-i$| parameter space. Thus each offset bin contains |$N_a = NM_a/bM$| objects.
Within an individual parameter offset bin, the probability of any individual object being detected is proportional to the number of |$\Omega -\nu$| combinations that result in a detectable target for the relevant |$h-i$| pair. In fact, the exact probability is equal to the fraction of all possible |$\Omega -\nu$| combinations that give a detectable result, i.e. |$P_a = C_{a,\rm detectable}/C_{\rm total}$| where |$C_{a,\rm detectable}$| is the number of detectable |$\Omega -\nu$| combinations for |$h-i$| bin a and |$C_{\rm total} = (360/0.1)^2$|, is the total number of |$\Omega -\nu$| combinations, since the |$\Omega$| and |$\nu$| offset bins are 0.1|$^{\circ }$| wide and can have any value from |$0^{\circ } \le \Omega ,\nu \lt 360^{\circ }$|.
The expected number of detections within a single offset bin, |$d_a$|, is thus equal to the number of objects in the relevant bin, multiplied by their probability for detection, |$d_a = N_aP_a$|. The total number of detections D, is thus the sum of |$d_a$| over all A bins, |$D = \sum {N_aP_a}$|. Expanding the terms and solving for N, thus gives a prediction of the number of objects N as a function of the number of detections D;
For a single pass the above equation holds true, i.e. tracking a TLE through one pass (as described in Sections 3 and 4) and making D detections will result in a predicted N objects in the relevant region of |$h-i$| parameter space. For multiple passes however, the average number of detections per orbit should instead be used. Since the above detectable offset combinations refer to a detection during a single pass, detections on subsequent passes give further information about the underlying population. Therefore, across p passes,
where D is the total number of detections across all p passes. Additionally, since the number of detections found can only be an integer number, the recorded value D should actually be seen as a lower bound on the true number of detections |$D_T$|, where |$D \le D_T \lt D+1$|. As such, the true number of objects, |$N_T$|, will lie in the range |$N(D) \le N_T \lt N(D+1)$|.
9 RESULTS
Based on the arguments above, it is possible to predict the number of RSOs in a particular region of orbital |$h-i$| parameter space by observing a number of passes of a theoretical orbit. Depending on the latitude of the observing site and the tracked |$h-i$| pair, the investigated region differs and the inferred population changes. To build up a full population of RSOs on all orbits would require significant investment of observing resources, tracking multiple passes of multiple theoretical orbits. To increase the resolution of the investigated |$h-i$| regions, the observations are better taken from higher latitude observing sites. Table 2 gives the results of the above simulation for a range of |$h/i/\delta$| combinations, as well as for a number of different maximum relative velocity thresholds and their corresponding magnitudes (using equation 6). The number of targets given assumes 1 detection per 20 passes, but can be simply scaled by the number of detections found (the chosen number of detections is simply an example, 20 passes was chosen as at least this number of fully observable passes per night is possible from most latitudes). The h, i, |$\Omega$|, |$\nu$| region columns give the maximum positive and negative offsets from the tracked orbital parameters that produce a detectable target.
Table of predicted number of RSOs inferred from observations of a theoretical tracked orbit. First four columns give the tracked orbit parameters, observation latitude, and limiting relative velocity used. Column five gives the corresponding limiting magnitude (using equation 6). Column six gives the predicted number of RSOs assuming 1 detection in 20 passes. Remaining columns give maximum positive and negative orbital element offsets that produce detectable targets.
| h (km) | i (|$^{\circ }$|) | |$\delta$| (|$^{\circ }$|) | max vel (pix s−1) | magnitude | N targets | h region (km) | i region (|$^{\circ }$|) | |$\Omega$| region (|$^{\circ }$|) | |$\nu$| region (|$^{\circ }$|) |
|---|---|---|---|---|---|---|---|---|---|
| 550 | 99 | 75 | 10.0 | 13.50 | 52510–105020 | |$h^{+48}_{-46}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+1.5}_{-1.5}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 2.5 | 16.24 | 405597–811193 | |$h^{+58}_{-54}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.2}_{-1.3}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 5.0 | 15.71 | 148832–297665 | |$h^{+58}_{-62}$| | |$i^{+1.2}_{-1.2}$| | |$\Omega ^{+1.8}_{-1.8}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 750 | 99 | 50 | 7.5 | 14.59 | 82380–164760 | |$h^{+66}_{-64}$| | |$i^{+1.2}_{-1.4}$| | |$\Omega ^{+1.9}_{-2.1}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 750 | 99 | 50 | 10.0 | 13.50 | 54370–108740 | |$h^{+72}_{-68}$| | |$i^{+1.5}_{-1.6}$| | |$\Omega ^{+2.3}_{-2.4}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 99 | 75 | 10.0 | 13.50 | 20767–41535 | |$h^{+68}_{-66}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 115 | 50 | 10.0 | 13.50 | 82984–165969 | |$h^{+82}_{-78}$| | |$i^{+1.4}_{-1.6}$| | |$\Omega ^{+2.7}_{-2.9}$| | |$\nu ^{+1.0}_{-1.0}$| |
| 850 | 95 | 75 | 10.0 | 13.50 | 31325–62650 | |$h^{+70}_{-68}$| | |$i^{+0.4}_{-0.5}$| | |$\Omega ^{+1.5}_{-1.7}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 29 | 2.5 | 16.24 | 695484–1390968 | |$h^{+106}_{-98}$| | |$i^{+2.2}_{-2.2}$| | |$\Omega ^{+2.1}_{-1.9}$| | |$\nu ^{+1.3}_{-1.1}$| |
| 850 | 99 | 29 | 5.0 | 15.71 | 254061–508121 | |$h^{+138}_{-118}$| | |$i^{+3.2}_{-3.0}$| | |$\Omega ^{+2.9}_{-2.5}$| | |$\nu ^{+1.6}_{-1.3}$| |
| 850 | 99 | 29 | 7.5 | 14.59 | 131839–263679 | |$h^{+148}_{-126}$| | |$i^{+3.6}_{-3.6}$| | |$\Omega ^{+3.0}_{-2.9}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 29 | 10.0 | 13.50 | 80392–160784 | |$h^{+156}_{-126}$| | |$i^{+3.9}_{-4.0}$| | |$\Omega ^{+3.1}_{-3.2}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 50 | 2.5 | 16.24 | 287108–574217 | |$h^{+58}_{-58}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.3}_{-1.3}$| | |$\nu ^{+0.5}_{-0.4}$| |
| 850 | 99 | 50 | 5.0 | 15.71 | 111222–222444 | |$h^{+62}_{-66}$| | |$i^{+1.3}_{-1.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 850 | 99 | 50 | 7.5 | 14.59 | 60585–121171 | |$h^{+70}_{-70}$| | |$i^{+1.5}_{-1.5}$| | |$\Omega ^{+2.3}_{-2.3}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 50 | 10.0 | 13.50 | 40188–80377 | |$h^{+78}_{-74}$| | |$i^{+1.9}_{-1.8}$| | |$\Omega ^{+3.0}_{-2.8}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 99 | 75 | 10.0 | 13.50 | 16218–32437 | |$h^{+76}_{-72}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 103 | 75 | 10.0 | 13.50 | 27063–54127 | |$h^{+78}_{-76}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.9}_{-0.8}$| |
| 950 | 99 | 75 | 10.0 | 13.50 | 13884–27769 | |$h^{+88}_{-86}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+2.0}_{-2.1}$| | |$\nu ^{+0.9}_{-0.8}$| |
| h (km) | i (|$^{\circ }$|) | |$\delta$| (|$^{\circ }$|) | max vel (pix s−1) | magnitude | N targets | h region (km) | i region (|$^{\circ }$|) | |$\Omega$| region (|$^{\circ }$|) | |$\nu$| region (|$^{\circ }$|) |
|---|---|---|---|---|---|---|---|---|---|
| 550 | 99 | 75 | 10.0 | 13.50 | 52510–105020 | |$h^{+48}_{-46}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+1.5}_{-1.5}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 2.5 | 16.24 | 405597–811193 | |$h^{+58}_{-54}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.2}_{-1.3}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 5.0 | 15.71 | 148832–297665 | |$h^{+58}_{-62}$| | |$i^{+1.2}_{-1.2}$| | |$\Omega ^{+1.8}_{-1.8}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 750 | 99 | 50 | 7.5 | 14.59 | 82380–164760 | |$h^{+66}_{-64}$| | |$i^{+1.2}_{-1.4}$| | |$\Omega ^{+1.9}_{-2.1}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 750 | 99 | 50 | 10.0 | 13.50 | 54370–108740 | |$h^{+72}_{-68}$| | |$i^{+1.5}_{-1.6}$| | |$\Omega ^{+2.3}_{-2.4}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 99 | 75 | 10.0 | 13.50 | 20767–41535 | |$h^{+68}_{-66}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 115 | 50 | 10.0 | 13.50 | 82984–165969 | |$h^{+82}_{-78}$| | |$i^{+1.4}_{-1.6}$| | |$\Omega ^{+2.7}_{-2.9}$| | |$\nu ^{+1.0}_{-1.0}$| |
| 850 | 95 | 75 | 10.0 | 13.50 | 31325–62650 | |$h^{+70}_{-68}$| | |$i^{+0.4}_{-0.5}$| | |$\Omega ^{+1.5}_{-1.7}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 29 | 2.5 | 16.24 | 695484–1390968 | |$h^{+106}_{-98}$| | |$i^{+2.2}_{-2.2}$| | |$\Omega ^{+2.1}_{-1.9}$| | |$\nu ^{+1.3}_{-1.1}$| |
| 850 | 99 | 29 | 5.0 | 15.71 | 254061–508121 | |$h^{+138}_{-118}$| | |$i^{+3.2}_{-3.0}$| | |$\Omega ^{+2.9}_{-2.5}$| | |$\nu ^{+1.6}_{-1.3}$| |
| 850 | 99 | 29 | 7.5 | 14.59 | 131839–263679 | |$h^{+148}_{-126}$| | |$i^{+3.6}_{-3.6}$| | |$\Omega ^{+3.0}_{-2.9}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 29 | 10.0 | 13.50 | 80392–160784 | |$h^{+156}_{-126}$| | |$i^{+3.9}_{-4.0}$| | |$\Omega ^{+3.1}_{-3.2}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 50 | 2.5 | 16.24 | 287108–574217 | |$h^{+58}_{-58}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.3}_{-1.3}$| | |$\nu ^{+0.5}_{-0.4}$| |
| 850 | 99 | 50 | 5.0 | 15.71 | 111222–222444 | |$h^{+62}_{-66}$| | |$i^{+1.3}_{-1.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 850 | 99 | 50 | 7.5 | 14.59 | 60585–121171 | |$h^{+70}_{-70}$| | |$i^{+1.5}_{-1.5}$| | |$\Omega ^{+2.3}_{-2.3}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 50 | 10.0 | 13.50 | 40188–80377 | |$h^{+78}_{-74}$| | |$i^{+1.9}_{-1.8}$| | |$\Omega ^{+3.0}_{-2.8}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 99 | 75 | 10.0 | 13.50 | 16218–32437 | |$h^{+76}_{-72}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 103 | 75 | 10.0 | 13.50 | 27063–54127 | |$h^{+78}_{-76}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.9}_{-0.8}$| |
| 950 | 99 | 75 | 10.0 | 13.50 | 13884–27769 | |$h^{+88}_{-86}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+2.0}_{-2.1}$| | |$\nu ^{+0.9}_{-0.8}$| |
Table of predicted number of RSOs inferred from observations of a theoretical tracked orbit. First four columns give the tracked orbit parameters, observation latitude, and limiting relative velocity used. Column five gives the corresponding limiting magnitude (using equation 6). Column six gives the predicted number of RSOs assuming 1 detection in 20 passes. Remaining columns give maximum positive and negative orbital element offsets that produce detectable targets.
| h (km) | i (|$^{\circ }$|) | |$\delta$| (|$^{\circ }$|) | max vel (pix s−1) | magnitude | N targets | h region (km) | i region (|$^{\circ }$|) | |$\Omega$| region (|$^{\circ }$|) | |$\nu$| region (|$^{\circ }$|) |
|---|---|---|---|---|---|---|---|---|---|
| 550 | 99 | 75 | 10.0 | 13.50 | 52510–105020 | |$h^{+48}_{-46}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+1.5}_{-1.5}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 2.5 | 16.24 | 405597–811193 | |$h^{+58}_{-54}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.2}_{-1.3}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 5.0 | 15.71 | 148832–297665 | |$h^{+58}_{-62}$| | |$i^{+1.2}_{-1.2}$| | |$\Omega ^{+1.8}_{-1.8}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 750 | 99 | 50 | 7.5 | 14.59 | 82380–164760 | |$h^{+66}_{-64}$| | |$i^{+1.2}_{-1.4}$| | |$\Omega ^{+1.9}_{-2.1}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 750 | 99 | 50 | 10.0 | 13.50 | 54370–108740 | |$h^{+72}_{-68}$| | |$i^{+1.5}_{-1.6}$| | |$\Omega ^{+2.3}_{-2.4}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 99 | 75 | 10.0 | 13.50 | 20767–41535 | |$h^{+68}_{-66}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 115 | 50 | 10.0 | 13.50 | 82984–165969 | |$h^{+82}_{-78}$| | |$i^{+1.4}_{-1.6}$| | |$\Omega ^{+2.7}_{-2.9}$| | |$\nu ^{+1.0}_{-1.0}$| |
| 850 | 95 | 75 | 10.0 | 13.50 | 31325–62650 | |$h^{+70}_{-68}$| | |$i^{+0.4}_{-0.5}$| | |$\Omega ^{+1.5}_{-1.7}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 29 | 2.5 | 16.24 | 695484–1390968 | |$h^{+106}_{-98}$| | |$i^{+2.2}_{-2.2}$| | |$\Omega ^{+2.1}_{-1.9}$| | |$\nu ^{+1.3}_{-1.1}$| |
| 850 | 99 | 29 | 5.0 | 15.71 | 254061–508121 | |$h^{+138}_{-118}$| | |$i^{+3.2}_{-3.0}$| | |$\Omega ^{+2.9}_{-2.5}$| | |$\nu ^{+1.6}_{-1.3}$| |
| 850 | 99 | 29 | 7.5 | 14.59 | 131839–263679 | |$h^{+148}_{-126}$| | |$i^{+3.6}_{-3.6}$| | |$\Omega ^{+3.0}_{-2.9}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 29 | 10.0 | 13.50 | 80392–160784 | |$h^{+156}_{-126}$| | |$i^{+3.9}_{-4.0}$| | |$\Omega ^{+3.1}_{-3.2}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 50 | 2.5 | 16.24 | 287108–574217 | |$h^{+58}_{-58}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.3}_{-1.3}$| | |$\nu ^{+0.5}_{-0.4}$| |
| 850 | 99 | 50 | 5.0 | 15.71 | 111222–222444 | |$h^{+62}_{-66}$| | |$i^{+1.3}_{-1.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 850 | 99 | 50 | 7.5 | 14.59 | 60585–121171 | |$h^{+70}_{-70}$| | |$i^{+1.5}_{-1.5}$| | |$\Omega ^{+2.3}_{-2.3}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 50 | 10.0 | 13.50 | 40188–80377 | |$h^{+78}_{-74}$| | |$i^{+1.9}_{-1.8}$| | |$\Omega ^{+3.0}_{-2.8}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 99 | 75 | 10.0 | 13.50 | 16218–32437 | |$h^{+76}_{-72}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 103 | 75 | 10.0 | 13.50 | 27063–54127 | |$h^{+78}_{-76}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.9}_{-0.8}$| |
| 950 | 99 | 75 | 10.0 | 13.50 | 13884–27769 | |$h^{+88}_{-86}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+2.0}_{-2.1}$| | |$\nu ^{+0.9}_{-0.8}$| |
| h (km) | i (|$^{\circ }$|) | |$\delta$| (|$^{\circ }$|) | max vel (pix s−1) | magnitude | N targets | h region (km) | i region (|$^{\circ }$|) | |$\Omega$| region (|$^{\circ }$|) | |$\nu$| region (|$^{\circ }$|) |
|---|---|---|---|---|---|---|---|---|---|
| 550 | 99 | 75 | 10.0 | 13.50 | 52510–105020 | |$h^{+48}_{-46}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+1.5}_{-1.5}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 2.5 | 16.24 | 405597–811193 | |$h^{+58}_{-54}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.2}_{-1.3}$| | |$\nu ^{+0.5}_{-0.5}$| |
| 750 | 99 | 50 | 5.0 | 15.71 | 148832–297665 | |$h^{+58}_{-62}$| | |$i^{+1.2}_{-1.2}$| | |$\Omega ^{+1.8}_{-1.8}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 750 | 99 | 50 | 7.5 | 14.59 | 82380–164760 | |$h^{+66}_{-64}$| | |$i^{+1.2}_{-1.4}$| | |$\Omega ^{+1.9}_{-2.1}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 750 | 99 | 50 | 10.0 | 13.50 | 54370–108740 | |$h^{+72}_{-68}$| | |$i^{+1.5}_{-1.6}$| | |$\Omega ^{+2.3}_{-2.4}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 99 | 75 | 10.0 | 13.50 | 20767–41535 | |$h^{+68}_{-66}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 750 | 115 | 50 | 10.0 | 13.50 | 82984–165969 | |$h^{+82}_{-78}$| | |$i^{+1.4}_{-1.6}$| | |$\Omega ^{+2.7}_{-2.9}$| | |$\nu ^{+1.0}_{-1.0}$| |
| 850 | 95 | 75 | 10.0 | 13.50 | 31325–62650 | |$h^{+70}_{-68}$| | |$i^{+0.4}_{-0.5}$| | |$\Omega ^{+1.5}_{-1.7}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 29 | 2.5 | 16.24 | 695484–1390968 | |$h^{+106}_{-98}$| | |$i^{+2.2}_{-2.2}$| | |$\Omega ^{+2.1}_{-1.9}$| | |$\nu ^{+1.3}_{-1.1}$| |
| 850 | 99 | 29 | 5.0 | 15.71 | 254061–508121 | |$h^{+138}_{-118}$| | |$i^{+3.2}_{-3.0}$| | |$\Omega ^{+2.9}_{-2.5}$| | |$\nu ^{+1.6}_{-1.3}$| |
| 850 | 99 | 29 | 7.5 | 14.59 | 131839–263679 | |$h^{+148}_{-126}$| | |$i^{+3.6}_{-3.6}$| | |$\Omega ^{+3.0}_{-2.9}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 29 | 10.0 | 13.50 | 80392–160784 | |$h^{+156}_{-126}$| | |$i^{+3.9}_{-4.0}$| | |$\Omega ^{+3.1}_{-3.2}$| | |$\nu ^{+1.7}_{-1.4}$| |
| 850 | 99 | 50 | 2.5 | 16.24 | 287108–574217 | |$h^{+58}_{-58}$| | |$i^{+0.9}_{-0.9}$| | |$\Omega ^{+1.3}_{-1.3}$| | |$\nu ^{+0.5}_{-0.4}$| |
| 850 | 99 | 50 | 5.0 | 15.71 | 111222–222444 | |$h^{+62}_{-66}$| | |$i^{+1.3}_{-1.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.6}_{-0.5}$| |
| 850 | 99 | 50 | 7.5 | 14.59 | 60585–121171 | |$h^{+70}_{-70}$| | |$i^{+1.5}_{-1.5}$| | |$\Omega ^{+2.3}_{-2.3}$| | |$\nu ^{+0.6}_{-0.6}$| |
| 850 | 99 | 50 | 10.0 | 13.50 | 40188–80377 | |$h^{+78}_{-74}$| | |$i^{+1.9}_{-1.8}$| | |$\Omega ^{+3.0}_{-2.8}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 99 | 75 | 10.0 | 13.50 | 16218–32437 | |$h^{+76}_{-72}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+1.7}_{-1.7}$| | |$\nu ^{+0.7}_{-0.7}$| |
| 850 | 103 | 75 | 10.0 | 13.50 | 27063–54127 | |$h^{+78}_{-76}$| | |$i^{+0.3}_{-0.3}$| | |$\Omega ^{+2.0}_{-2.0}$| | |$\nu ^{+0.9}_{-0.8}$| |
| 950 | 99 | 75 | 10.0 | 13.50 | 13884–27769 | |$h^{+88}_{-86}$| | |$i^{+0.4}_{-0.4}$| | |$\Omega ^{+2.0}_{-2.1}$| | |$\nu ^{+0.9}_{-0.8}$| |
To correctly make inferences about the total number of targets in a particular orbital region we must first limit the relative velocity of any detections. The predicted number of targets is inferred from the number of detections but requires an assumption of an upper limit to the relative velocity of a target. From Figs 8 and 9 however, we see that it is possible for a target moving faster than the imposed velocity threshold to be detected, if the target is sufficiently bright. Therefore, we must first filter our detections, rejecting any detection moving above the chosen relative velocity threshold. Since the relative velocity of a target changes over a single pass, this means we reject only those targets whose slowest detected relative velocity is above the threshold. This ensures the detections we count fall into the simulated parameter space, from which we can make inferences about the total number of targets.
Conversely, for targets moving below the relative velocity threshold, the limiting magnitude is not constant, and thus the magnitude to which the detection sample will be complete depends on the exact velocity for each target, not just the defined threshold. Thus, a target moving below the velocity threshold may or may not be detectable depending on its magnitude. To account for this, we must again filter our detection sample, rejecting detections for which |$M\gt M(v)$|, i.e. targets fainter than the magnitude corresponding to the defined relative velocity threshold. This ensures that our sample of detections is complete down to a magnitude of |$M(v)$|. The numbers given in column six of Table 2 are therefore interpreted as being the predicted number of targets within the given |$h-i$| parameter space brighter than the given magnitude.
The result of this filtering approach means that we can filter our detections multiple times, each time producing a slightly different inferred population of targets. From a full set of detections from p passes, we can select multiple different relative velocity thresholds, filter the detections accordingly, and thus infer multiple different RSO populations without requiring additional data. As such, Table 2 shows the results of using different velocity thresholds. As the chosen velocity threshold decreases, the magnitude to which the sample is complete becomes fainter, and the resolution on the |$h-i$| region of parameter space becomes higher.
10 DISCUSSION AND CONCLUSIONS
The work presented here is a description and simulation of a novel method for inferring RSO populations by observing neighbouring orbits. By making certain assumptions, this technique allows one to infer the number of RSOs in particular |$h-i$| orbital parameter bins. The method is statistical, extrapolating detections based on an underlying distribution, to determine the full population. A complete survey would require observing multiple passes while tracking orbits with a range of |$h-i$| combinations. Ideally, this survey would also be carried out across multiple observing sites situated across a range of latitudes so as to analyse |$h-i$| bins with a range of resolutions (as seen above, observing the same |$h-i$| region at different latitudes allows inferences about different sized regions around the chosen |$h-i$| pair).
Since the proposed method relies on detecting objects with low relative velocity, it can detect fainter than other methods at LEO. Objects in LEO have high orbital velocities and as such, spread their flux over a large number of pixels when observed with a stationary telescope. This can be improved upon using bespoke techniques [i.e. blind stacking, as discussed in Cooke et al. (2023)] but the target streaks are still hundreds of pixels long in a single exposure. The above method detects neighbouring orbit targets with streaks less than 10 pixels per exposure. Depending on the chosen threshold for relative velocity, this method can achieve 50 per cent recovery of targets down to |${\sim }\, 16^{\rm th}$| magnitude, a significant improvement over other methods tested on the same equipment and at the same observing site (Cooke et al. 2023). The downside, of course, is that the parameter space of detections is much more constrained, only able to detect targets moving on similar orbits to the targeted one. Therefore, there is a trade-off between different methods, depending on the goals of the survey; a more complete population approach, or a targeted analysis of RSOs on specific orbits.
This paper lays out the neighbouring orbits technique in terms of simulations. The next step would be to begin a campaign of observations of a range of |$h-i$| orbital pairs and make predictions based on detections. Once a set of detections has been established, the inferred number of detections can be calculated and compared to other estimates, generated through theoretical and/or simulated populations, resulting in an independent validation of these methods. A full survey of all potential orbital heights and inclinations would be a significant undertaking, but focusing on a smaller subset of |$h-i$| pairs of particular interest would be more feasible. A future work for this project will be to begin a campaign of this type, most likely situated at the observing site on La Palma (|$\delta =29^{\circ }$|). Collaboration with, or data from, other observing sites would allow for a more complete picture of the orbital distribution of RSOs to be obtained.
Acknowledgement
BFC and PC acknowledge support from a Science and Technology Facilities Council CLASP award (grant ST/V002279/1) and from the Defence Science and Technology Laboratory (UK). JAB gratefully acknowledges support from the Defence Science and Technology Laboratory (UK) and the Science and Technology Facilities Council (grant ST/Y50998X/1).
This work has made use of data obtained using the Warwick CLASP test telescope operated on the island of La Palma by the University of Warwick in the Spanish Observatory del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.
The authors would like to thank an anonymous referee for their useful comments.
For the purpose of open access, the author has applied a Creative Commons Attribution (CC-BY) licence to any Author Accepted Manuscript version arising from this submission.
DATA AVAILABILITY
Relevant data are freely available under a CC-BY-4.0 license at https://github.com/BenCooke95/RSO_populations. Further details may be granted upon reasonable request to the corresponding author.
References
APPENDIX A: ORBITAL PARAMETER DERIVATIONS
A1 Definitions
Let the relevant orbit be circular, with radius r and inclination |$\theta$|. The orbit should pass through the zenith of the observing site at the chosen epoch.
Let the observing site, have |$\rm {latitude}=\delta$|, |$\rm {longitude}=\alpha$|.
Let a defined orbit cross the equator in the positive z-direction at |$(r, \frac{\pi }{2}, \phi _e)$|.
Let |$\phi _{RAAN}$| be the Right Ascension of the Ascending Node (RAAN) of a polar orbit that passes directly over observing site.
Orbital parameter |$\Omega$| is equal to the RAAN of the observing site (equal to the sidereal time) plus the difference between longitude of the observing site and longitude of the point at which the orbit crosses the equator in the positive z-direction. Orbital parameter |$\nu$| is equal to the great circle angle between the observing site and the point at which the orbit crosses the equator in the positive z-direction.
A2 Derivations
Define a circular orbit as the intersection of a sphere of radius r at the origin, and a plane. To define a plane, we need three distinct points on plane:
Origin: |$P_1 = (x_1, y_1, z_1) = (0, 0, 0)$|
Observing site: |$P_2 = (x_2, y_2, z_2) = (r, \theta _2, \phi _2)$|
Point in orbit at maximum z: |$P_3 = (x_3, y_3, z_3) = (r, \theta _3, \phi _3)$|
In Cartesian coordinates:
Equation of plane, P, from |$P_1, P_2, P_3$|:
For any point |$Q(x, y, z) = Q(r, \theta , \phi )$| on plane P:
A2.1 Calculation of |$\phi _e$| and n
Let Q be the point at which orbit crosses equator in positive z direction. Thus |$Q(x, y, z) = Q(r, \frac{\pi }{2}, \phi _e)$|.
Thus, substituting equations (A1) and (A3) into equation (A2):
Note that |$\theta _3 = \theta$|. From symmetry, |$\phi _3$| is exactly |$\frac{\pi }{2}$| around equator from |$\phi _e$|:
Thus:
Note that |$P_2$| represents the observing site, thus |$\phi _2=\alpha$|, |$\theta _2=\frac{\pi }{2}-\delta$|. Therefore |$\sin {\theta _2} = \cos {\delta }$|, |$\cos {\theta _2} = \sin {\delta }$|. Substituting into equation (A6) gives:
Where |$A = \sin {\alpha }\cos {\delta }\cos {\theta }$|, |$B = \cos {\alpha }\cos {\delta }\cos {\theta }$|, |$C = \sin {\theta }\sin {\delta }$|, defined in terms of |$\alpha$|, |$\delta$|, and |$\theta$|.
Using the double angle formula:
Let |$\phi _{e2} = \frac{\phi _e}{2}$|, therefore:
Thus, taking the positive root:
Let |$\phi _d = \alpha -\phi _e$|. Thus, finally:
A2.2 Calculation of |$\sigma$| and |$\nu$|
Utilize the great circle angle formula [for example derivation see Rajpoot (2016)]. |$\theta$| measured from equator to pole, |$\phi$| measured around equator. |$\sigma$| is angle between |$P(r, \theta _1, \phi _1)$| and |$Q(r, \theta _2, \phi _2)$|.
P is point at which orbit crosses equator, thus |$\theta _1 = 0$|, |$\phi _1 = \phi _e$|. Q is observing site, thus |$\theta _2 = \delta$|, |$\phi _2 = \alpha$|. Therefore:
Orbital projections. Green circle represents the Earth, black, the Cartesian coordinate axes and construction lines, and red, the orbit, with an arrow showing direction of motion (dotted portion shows orbit behind the earth). Three individual points are also shown; blue, for the observing site (middle point in (a) and (b)), magenta, for the orbital point with maximum z-coordinate (upper point in (a), lower point in (b)), and cyan, for the point at which the orbit cross the equator in the positive z-direction (lower point in (a), upper point in (b)). The observing site in Fig. A1(b) appears to lie off the line of the orbit due to angular projection effects. Projections are constructed using |$\theta =20^{\circ }$|, |$\alpha =\delta =40^{\circ }$|, |$r=1.1R_\oplus$| for illustration purposes.
APPENDIX B: TARGET DETECTABILITY COMPARISON PLOTS
Reproduction of Figs 4, 5(a) and (b) on the same axis scale. Target detectability as a function of |$h_{\rm offset}$| and |$i_{\rm offset}$| for |$O_{\rm principal} = (850, 99, 29.8, 28.9)$|. Colour designates relative probability, i.e. number of |$\Omega -\nu$| combinations resulting in a detection for each |$h-i$| pair. Values assume a limiting relative velocity of 10 pix s−1.
