1. Introduction
Radio astronomy presents an opportunity to study the early Universe, as radio emission generated billions of years ago from neutral hydrogen is redshifted to low radio frequencies (100–200 MHz) (Morales & Wyithe Reference Morales and Wyithe2010). These signals are expected to be weak, meaning that extremely sensitive instruments are needed to detect them and extract astrophysical information from vast data collections. In order to address this challenge, the most sensitive radio telescope of all time is being constructed, known as the Square Kilometre Array (SKA) (Caiazzo Reference Caiazzo2017).
Civilisation also heavily utilises the radio spectrum for the purpose of television broadcast, FM radio, mobile phone communications, radio-navigation, and space-based communications, to name a few. In space-based communications, satellites transmit radio frequency energy to the surface of the Earth, to a ground station or device. If the satellite is transmitting information whilst passing through the receiving beam of a radio telescope, the strength of its transmission can be many orders of magnitude greater than the signals from even the strongest astronomical sources. This poses a significant and well-known risk for radio astronomy research (Gilloire & Sizun Reference Gilloire and Sizun2009; Zheleva et al. Reference Zheleva, Anderson, Aksoy, Johnson, Affinnih and DePree2023; Di Vruno et al. Reference Di Vruno, Winkel, Bassa, Józsa, Brentjens, Jessner and Garrington2023; Grigg et al. Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023).Footnote a
The rate at which human-made objects are launched into space is also rising exponentially.Footnote b All of these need to be tracked to maintain a harmonious space environment, and the more objects there are, the more comprehensive the system tracking them needs to be. At the time of writing this paper, space-track.org tracks approximately 10 200 active satellites, 18 700 debris objects, and 16 700 objects which are variably tracked due to their small size or incomplete data, all in Earth orbit. The intrinsic risk of interference at radio frequencies from satellites operating in Earth orbit is increasing, simply by virtue of the greatly increased numbers of operating objects.
The International Telecommunication Union (ITU) manages satellite spectrum allocations in the radio band under its radiocommunication sector (ITU-R).Footnote c In terms of protections for radio astronomy from Earth, only approximately 3.7% of the spectrum allocations in the frequency range 50–350 MHz receive some protection for radio astronomy. Primary or secondary allocations for radio astronomy exist at 73.00–74.60 MHz, 150.05–153.00 MHz, and 322.00–328.60 MHz across specific geographic regions in the ITU-R radio regulations.Footnote d Footnote 5.149 in this document also states that ‘administrations are urged to take all practicable steps to protect the radio astronomy service from harmful interference’ and applies to the frequency ranges 73.00–74.60 MHz and 322.00–328.60 MHz.
In the 50–350 MHz band, several allocations to the mobile-satellite, meteorological-satellite, and amateur-satellite services exist in the radio regulations, both as outright allocations and under footnote 5.254
$^{\text{d}}$
for defence use. These can significantly impact telescopes such as the SKA, for which the low-frequency half (SKA-Low) is being built in Western Australia. This is due to the historical nature of the regulations and the historical claims to the radio astronomy spectrum being limited and also due to the rapid acceleration of satellite and launch technologies.
Historically, radio telescopes have successfully managed the limited regulatory protection of the spectrum for radio astronomy by operating from remote areas, largely free from Radio Frequency Interference (RFI). Examples include the Murchison Widefield Array (MWA, Tingay et al. Reference Tingay2013) and Australia Square Kilometre Array Pathfinder (ASKAP, Johnston et al. Reference Johnston2008) both in Western Australia, the MeerKAT telescope (Jonas Reference Jonas2009) in South Africa, and the Very Large Array (Condon et al. Reference Condon, Cotton, Greisen, Yin, Perley, Taylor and Broderick1998) in New Mexico, USA. With satellites essentially visible from anywhere on the Earth, geographic isolation does not mitigate against RFI originating from satellites, highlighted by work using the MWA and Engineering Development Array 2 (EDA2, Tingay et al. Reference Tingay2013; Prabu et al. Reference Prabu, Hancock, Zhang and Tingay2020b; Tingay et al. Reference Tingay, Sokolowski, Wayth and Ung2020; Grigg et al. Reference Grigg, Tingay, Sokolowski and Wayth2022).
The SKA-Low telescope has a planned observing frequency of 50–350 MHz, and satellites may impact the pursuit of its headline science goals. This type of scenario is already a reality. An example is the ALMA radio telescope in Chile, which needed to introduce restrictions on observing at the same frequency as CloudSat transmissions, as there was a non-negligible risk of damage to its receivers.Footnote e
The work we describe in this paper is a first step towards characterising the impact of all satellites, specifically across the full SKA-Low frequency range of 50–350 MHz at the SKA-Low site. In this work, we test various methods of acquiring and processing data to understand the optimal way to conduct a more comprehensive survey in the future.
The radio telescopes we use are the Aperture Array Verification System, version 2 (AAVS2, van Es et al. Reference van Es, Marshall, Spyromilio and Usuda2020) and the Engineering Development Array, version 2 (EDA2, Wayth et al. Reference Wayth2021). These are pictured in Fig. 1. The two SKA-Low prototype stations were chosen because they are representative of a single SKA-Low station. The two telescopes are on the same site as the SKA-Low in Western Australia, Inyarrimanha Ilgardi Bundara, the CSIRO Murchison Radio-astronomy Observatory.Footnote f This site has been designated a Radio Quiet Zone by national legislation (Wilson et al. Reference Wilson, Chow, Harvey-Smith, Indermuehle, Sokolowski and Wayth2016). This makes the AAVS2 and EDA2 ideal for understanding the environment in which the SKA-Low will operate, especially as it pertains to sources of RFI.
Figure 1. A comparison of the EDA2 (left) and AAVS2 (right) low-frequency radio telescopes used in this work.
Two methods are explored to detect satellites with a radio telescope. The first is the detection of direct transmissions from the satellite. The mechanisms behind these types of radio emissions are defined as follows: ‘intentional electromagnetic radiation’ (IEMR) is intentional transmission from a satellite’s antenna at its designated downlink frequency; ‘interference’ is any additional energy from the satellite’s transmission antenna that is outside the designated downlink frequency; and ‘unintentional electromagnetic radiation’ (UEMR) is energy emitted from another component of the satellite that is not the antenna (for example the propulsion system). IEMR is regulated by the ITU’s Radio Regulations (ITU-R),Footnote g whereas interference and UEMR are not currently regulated by any international regulatory bodies. As per the ITU-R definition, anything which is IEMR is henceforth referred to as ‘emission’.
The second method is the detection of satellites via the reflection of terrestrial sources of radio emission. This is a form of radar, and in radar parlance is described as passive non-coherent bi-static radar (Griffiths Reference Griffiths, Sidiropoulos, Gini, Chellappa and Theodoridis2014). It is passive, meaning that no signal is deliberately emitted by the observer to make a detection, and is bi-static, meaning that the source (illuminator) and the receiver are spatially separated. It is non-coherent because the received signal is not compared to reference signals (the illumination signal) to achieve phase coherence in order to allow for a range measurement.
One prominent source of illumination in this scenario is commercial or private FM radio transmitters distributed throughout Australia, which can output up to 250 kW of power over 0.2 MHz of bandwidth.Footnote h This transmission is directed over the land, but due to complex transmission beams, significant energy can be directed away from the Earth (Hennessy et al. Reference Hennessy, Rutten, Young, Tingay, Summers, Gustainis, Crosse and Sokolowski2022). The radio waves can reflect off objects in the local airspace or in Earth orbit and subsequently be detected by radio telescopes operating in this frequency range. Other sources of illumination are also possible, for example DTV broadcasting, and many other licenced sources of terrestrial and airborne transmission.
Table 1. Surveyed frequencies
Figure 2. The 15 frequencies of datasets used in this work alongside noteworthy spectrum allocations.
Previous studies have developed processes for the detection of satellites in AAVS2 and EDA2 datasets, based on both these methods of detection. An initial study by Tingay et al. (Reference Tingay, Sokolowski, Wayth and Ung2020) began by acquiring three days worth of
${\approx}$
1 MHz bandwidth data overlapping several known FM transmitters using the EDA2 and making seven ‘by eye’ detections. A follow up study by Grigg et al. (Reference Grigg, Tingay, Sokolowski and Wayth2022) developed and utilised autonomous detection algorithms on the same data set, resulting in a 6 fold increase in the detection rate. Work by Sokolowski et al. (Reference Sokolowski2021) used the AAVS2 and EDA2 in tandem over a small selection of frequencies within the SKA-Low observing bandwidth and detected 30 unique objects (satellites or debris) as part of a larger survey of RFI at the site. Most recently, Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023) detected IEMR and UEMR from Starlink satellites using the EDA2, from some of the same datasets that will be reported on in the current paper. Separately to these, Di Vruno et al. (Reference Di Vruno, Winkel, Bassa, Józsa, Brentjens, Jessner and Garrington2023) and Bassa et al. (Reference Bassa, Di Vruno, Winkel, Józsa, Brentjens and Zhang2024) also detected UEMR from Starlink satellites using the LOFAR radio telescope (van Haarlem et al. Reference van Haarlem2013).
Building on our previous work using AAVS2 and EDA2, our current work uses twelve frequencies (
${\approx}$
1 MHz bandwidth per frequency) which are listed in Table 1 and illustrated in Fig. 2. These frequencies were selected based on a range of factors, in an attempt to span the full SKA-Low frequency range at discrete frequency intervals. Our purpose is to take the first step in systematically surveying the full SKA-Low frequency range, in order to understand the impacts that RFI from satellites (and other airborne objects such as aircraft) may have on the SKA-Low.
Improving autonomy of the detection algorithms to delineate between satellites and RFI is one of the key goals of this work, which is quantified with a simulation to estimate a misidentification rate. We explore data processing and signal identification methods and compare them with other recent studies (Grigg et al. Reference Grigg, Tingay, Sokolowski and Wayth2022; Sokolowski et al. Reference Sokolowski2021). This study will be a key step in informing the design of a larger, more comprehensive survey in the future. This survey will also be a baseline against which to compare future monitoring of RFI from satellites.
Section 2 presents the methods we have used to collect and process data from the AAVS2 and the EDA2. Section 3 presents our results, and discussion of the results is in Section 4. In Section 5, we draw our conclusions from this work.
2. Methods
This study was designed to utilise frequency channels within the SKA-Low frequency range of 50–350 MHz. Frequency channels which were known to be contaminated with RFI (from satellites, aircraft, meteors, and terrestrial transmitters) were intentionally chosen, along with others which were thought to contain little RFI. This Section will examine the choice of frequencies, data processing techniques, and a simulation used to quantify the expected misidentification rate.
2.1. Instrumentation
The two radio telescopes used for this work were the AAVS2 (van Es et al. Reference van Es, Marshall, Spyromilio and Usuda2020) and EDA2 (Wayth et al. Reference Wayth2021). They both feature an arrangement of 256 antennas within a diameter of
${\approx}$
38 m and
${\approx}$
35 m, respectively. The difference between the two telescopes is that the AAVS2 uses SKA prototype antennas (known as SKALA) and the EDA2 uses MWA bowtie antennas. The functional bandwidth of the two telescopes is comparable to the SKA-Low being approximately 50–350 MHz.
In previous internal tests, it was found that the AAVS2 had a higher sensitivity at the horizon. This increased the detrimental effect of transmission directly from terrestrial radio transmitters, resulting in a higher percentage of unusable data. The AAVS2 was therefore only used to acquire higher frequency data, making use of its wider frequency range compared to the EDA2.
2.2. Choice of frequencies
The choice of specific frequencies for this work was driven by their practical significance and scientific relevance. These are listed in Table 1 and illustrated in Fig. 2. In total, almost 20 days of telescope recording time was collected, comprising 1.6 million images over a cumulative total of
${\approx}$
12 MHz of bandwidth.
The ‘notes’ column in Table 1 describes the reasoning behind the choice of each frequency. These will be briefly explained.
We choose frequencies at 73.4, 150.8, 324.2, and 325 MHz to overlap with ITU-R protected frequencies for radio astronomy to search for satellite transmissions or UEMR.
Work by Sokolowski et al. (Reference Sokolowski2021) had previously identified satellites at frequencies of 159.4 and 229.7 MHz. These were again chosen to compare satellite identifications in our study with Sokolowski et al.’s work. Differences in the number and type of detections made in these two studies prompted a reprocessing of two
${\approx}$
130 h datasets acquired by the AAVS2 and EDA2 from Sokolowski et al.’s paper (on 2020-06-26 both at 159.4 MHz) to determine if our new algorithms result in a higher detection rate and lower misidentification rate. Note that a STARLINK train has been discovered transmitting UEMR in the 2021-11-16 159.4 MHz dataset, but is omitted from this analysis as it has already been reported by Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023).
It is interesting to understand whether radio astronomy research can be performed at or near frequencies which are allocated to satellite transmission. The frequency of 137.5 MHz was chosen for this reason, as many satellites have their designated downlink frequency within this band. This same dataset was examined in the recent paper by Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023), only describing the STARLINK satellite transmission. We now describe the other non-STARLINK detections in this dataset.
It is known that FM radio waves reflect off satellites (Prabu et al. Reference Prabu, Hancock, Zhang and Tingay2020a; Hennessy et al. Reference Hennessy, Rutten, Young, Tingay, Summers, Gustainis, Crosse and Sokolowski2022; Tingay et al. Reference Tingay, Sokolowski, Wayth and Ung2020), but it is unknown whether reflections of digital television transmission can also be detected. The frequency of 185.2 MHz was chosen to overlap with four DTV transmitters in Western Australia, which are listed in Table 2. Note that satellite detections in the FM band had been described in Grigg et al. (Reference Grigg, Tingay, Sokolowski and Wayth2022), which was why longer observational periods were not acquired at FM frequencies.
Table 2. Notable terrestrial transmitters. Nominal bandwidth of the transmitters is 0.2 MHz.
The two frequencies 146.1 and 160.2 MHz were pre-existing in the EDA2 data archive and were made available for this work.
Three frequencies were also chosen to test a new method of image processing called ‘frequency differencing’, explained in more detail in Section 2.5. The 96.9 and 98.4 MHz datasets were chosen to coincide with passes of the International Space Station (ISS; NORAD 25544) at frequencies which overlapped the terrestrial FM transmitters listed in Table 1. The 136.7 MHz dataset coincides with a pass of the satellite SOLRAD 7B (NORAD 1291). Instead of recording data at a single ‘coarse’ 0.926 MHz frequency channel, these four datasets were captured using 32 ‘fine’ frequency channels (28.9 kHz bandwidth per channel). As the 32 fine channels could be averaged into a single 0.926 MHz coarse channel, these datasets were also used in the analysis in Section 2.3.
2.3. Data processing
Section 2.3 is focused on the analysis using the time differencing technique, which is described below.
2.3.1. Pre-processing
The antenna voltages acquired from the AAVS2 and EDA2 were cross-correlated for all baselines. Table 1 lists the correlated integration time for each observation, with each acquisition capturing 0.926 MHz of total bandwidth. For the single channel datasets, these visibilities output by the correlator were calibrated using the procedure described by Sokolowski et al. (Reference Sokolowski2021), which is an extension of the procedure used in Benthem et al. (Reference Benthem2021). This calibration corrects for phase variations (due to residual unaccounted for cable delays) and amplitude (to bring to an absolute flux density scale). It also assumes a quiet Sun model, and while this may not give the highest level of accuracy near solar maximum, it is adequate for this analysis where the astrometry of the signals is the most important factor.
Subsequent to calibration, visibilities were then imaged and cleaned with the radio astronomy data processing package MIRIAD (Sault, Teuben, & Wright Reference Sault, Teuben, Wright, Shaw, Payne and Hayes1995), for which the procedure is outlined in the work of Tingay et al. (Reference Tingay, Sokolowski, Wayth and Ung2020). This produced full-sky images in both the XX (east-west-oriented) and YY (north-south-oriented) orthogonal polarisations.
For the four datasets with 32 channels, a different correlation package was used to handle the fine channelisation. This meant that the only difference with the processing to the image stage was that amplitude calibration was not performed. A simplified amplitude calibration was performed instead at the image stage to allow the amplitudes in these datasets to be comparable to the single channel datasets. This will be briefly explained.
The 2022-12-18 96.9 MHz dataset had Hydra A visible in the images. The Hydra A beam corrected flux density was fitted for in each polarisation, and using the polynomial correction function and coefficients measured by Perley & Butler (Reference Perley and Butler2017), the measured image intensities were multiplied by the correction factors for each polarisation. The 2022-12-16 96.9 MHz dataset had no astronomical radio sources visible in the dataset, which meant that the same correction had to be applied from the dataset two days earlier. The 2022-10-18 98.4 MHz dataset had Fornax A visible and was also corrected for with the same method, also using coefficients measured by Perley & Butler (Reference Perley and Butler2017). The last dataset at 136.7 MHz had the Sun visible in the images. The same quiet Sun model from Tingay et al. (Reference Tingay, Sokolowski, Wayth and Ung2020) was used to derive the bulk correction applied to the images for each polarisation.
For the primary satellite detection and identification, a technique called ‘time differencing’ was used to process the images. Satellites move very quickly over the sky compared to the movement of astrophysical sources, from the point of view of an observer on the Earth. Low Earth orbit (LEO) satellites can for example can pass overhead up to
$\approx$
1 degree s
$^{-1}$
, while the sky moves up to 15 degrees h
$^{-1}$
. Therefore images of the sky, taken close to each other in time, are differenced to isolate the signals which traverse the sky quickly. This is shown in the left column of Fig. 3 and is a technique which has been used in previous similar studies (Tingay et al. Reference Tingay, Sokolowski, Wayth and Ung2020; Grigg et al. Reference Grigg, Tingay, Sokolowski and Wayth2022; Grigg et al. Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023). Other variations in amplitudes between the two images will also be visible in the time differenced images. Examples include other transmitters which vary in directionally transmitted power, and scintillation of astrophysical sources.
Figure 3. The time and frequency differencing techniques are demonstrated with these images of the full sky, encompassed by the blue grid. The red arrows show the location of the ISS. For time differencing, image two is subtracted from image one with 40 s separating the two images with a central frequency of 96.9 MHz over a 28.9 kHz bandwidth. This gives image three which is the ‘time differenced’ image and clearly shows the separation of the reflected FM transmission from the ISS. For frequency differencing, image four (96.90 MHz) shows the ISS, while in image five (97.09 MHz) the ISS is not visible (both have a 28.9 kHz bandwidth). Image five is subtracted from image four to give the ‘frequency differenced’ image shown in image six. It is important to note that the time difference method retains the negative subtracted ISS from the previous frame.
All datasets in Table 1 were time differenced with a nominal separation of 40 s between the two differenced images, for every image other than the first 40 s of the dataset. This allowed a satellite in low earth orbit to be sufficiently separated in the two images used for differencing. Image three in Fig. 3 shows this effect with the positive (white) amplitude emission being separated from the subtracted negative (black) amplitude emission. This met a trade off between the satellite’s emission being separated enough to not overlap, whilst not introducing artefacts from imperfect background subtraction due to the movement of the sky between the two images. The satellite detection and identification algorithms were performed on these time differenced images.
In this paper, we refer to ‘candidates’ as the detection of some signal for which the origin is not yet certain. This could be from a satellite, or noise, or any high flux density signal which is detected by the algorithms outlined in Section 2.3.2. The candidates are then subject to a sorting algorithm which takes all of the candidates and algorithmically sorts them to determine if any are likely to be caused by a satellite. If a specific satellite is attributed to a given detection then we call it an identification. The identification algorithms are described in Section 2.3.3.
2.3.2. Candidate detections
A brute force fitting algorithm was applied to the images to make detections. The XX and YY orthogonal polarisations were run separately through this algorithm. The XX polarisation images are more sensitive to the North and the South, and the YY images are more sensitive to the East and West. If these were combined to form a Stokes I image, that extra sensitivity to the horizon due to one polarisation would be averaged out by the lack of sensitivity in the other, which is the reason the two polarisations are kept separated.
The satellite positional information used in this work was in the form of two line element sets (TLEs),Footnote i which describe six instantaneous Keplerian orbital elements at an epoch. Using the SGP4 propagation routine in the Skyfield python library (Rhodes Reference Rhodes2019), these TLEs were used to calculate the positions of satellites for a specific image. All TLEs between two days either side of the start and end times of a given survey were collected. A perigee limit of 2 000 km was placed on the TLEs, as satellites in higher orbits would be moving too slow to be visible in the time range of the difference images.
The brute force fitting then proceeds as follows: for each time differenced image, the pixel locations of all satellites visible above the horizon were calculated. An elevation cutoff of 20
$^\circ$
was imposed to prevent transmitters on the horizon being detected, as well as limiting the effect of photometric centroid shift which is explained further in Appendix A. The two dimensional elliptical Gaussian function shown in Equation (1) is a reasonable approximation of the primary lobe of the point spread function (PSF) of the AAVS2 and EDA2. An elliptical Gaussian function is needed, as a LEO satellite at a higher frequency can move further than the width of the PSF.
For example, a satellite at a zenithal range of 400 km will travel 2.2
$^\circ$
in 2 s (the integration time of the image). At a frequency of 328.1 MHz, the width of the PSF is
$\approx1.5^\circ$
. This motion, convolved with the PSF will create a streak, which when fitted with an elliptical Gaussian has fitting uncertainties two orders of magnitude below the 5
$\%$
threshold.
This Gaussian does not need to be normalised, as the peak flux density of the fitted Gaussian gives the flux density of the received signal. An attempted fit for this function at the location of every satellite on the sky is made, meaning that all sources of radio emission at that frequency will be fit, regardless of their origins.
\begin{equation} I_{\text{model}} = A \cdot \exp\left(-\left[\frac{(x - x_{0})^{2}}{2 {\sigma_{x}}^2} + \frac{(y - y_{0})^{2}}{2 {\sigma_{y}}^2}\right]\right), \end{equation}
The fitting algorithm used was from the Scipy library (Virtanen et al. Reference Virtanen2020), which in this case estimates the Gaussian’s amplitude (A), x coordinate (
$x_{0}$
), y coordinate (
$y_{0}$
), and the standard deviation of the Gaussian in the y direction (
$\sigma_{y}$
), which is related to the full width at half maximum (FWHM) of the Gaussian by FWHM
$= 2\sigma\sqrt{2ln(2)}$
.
$\sigma_{x}$
is known from the width of PSF at the given frequency, and as such is not fitted for in this equation. An estimate of each of these was provided to constrain the fits, with A estimated to be the maximum pixel value in a
$3\times 3$
kernel centred on
$x_{0}$
and
$y_{0}$
, the predicted TLE pixel values of the satellite, and
$\sigma_{y}$
was estimated from the width of the PSF. The standard deviation on each fitted parameter was calculated by taking the square root of the diagonalised covariance matrix of the fit. These standard deviations were treated as a measure of the uncertainty on each of the fitted parameters.
If a fit was successful for a given TLE position of a satellite, it had to pass through several constraints before it was deemed a candidate:
-
1.
$A \gt $
1 Jy/beam; -
2.
$A_{err} \gt 0.05 A$
and
$\sigma_{err} \gt 0.05 \sigma$
; -
3. The separation between the TLE predicted position and the fitted position could not be greater than the larger of two pixels or three degrees. Three degrees can be a fraction of a pixel at lower elevations, which is why a pixel and degree threshold is given; and
-
4. If the RMS of the difference image was less than 10 000 Jy/beam
$^{-1}$
(meaning that radio galaxies and the Sun were likely visible if overhead), the coordinates of the fit could not be inside an exclusion radius around the Sun and other radio galaxies which was defined as a function of the point spread function width at that frequency.
The final part of this first processing step was to perform a beam correction on the detected flux density of the candidate for each polarisation. This was calculated by dividing A by the value of the normalised antenna beam at the azimuth/elevation location of the candidate for each polarisation.
2.3.3. Identifications
In this Section, a qualitative description of how identifications are made from the list of candidates is described, with the reproducible mathematical steps and quantitative analysis available in Appendix B.
For each survey, the list of candidates needed to be sorted to ensure satellites could be correctly identified. In this list, there may be multiple candidates detected at the same location, at the same time. This needs to be reduced to one candidate (for an identification) or zero candidates (for RFI). We now refer to the TLE predicted position of a satellite as the ‘predicted’ position, and the position of the fitted Gaussian (
$x_{0}$
and
$y_{0}$
) as the ‘measured’ position. Each candidate will have this predicted and measured position information.
This proves to be valuable, as when the direction of travel of the predicted and measured positions is similar over time, it increases the probability of a correct identification. If the predicted and detected positions are travelling in different directions over the sky, then these candidates can be discarded. Therefore, knowing the instantaneous direction of travel of multiple candidates potentially associated with the same satellite can be used to discard infeasible candidates for that satellite.
The list of candidates is then sorted, so all candidates for each satellite are put into 15 min bins which are designed to be long enough that any object in the TLE list will pass over the sky completely during this time. The TLE list includes satellites which have a perigee up to 2 000 km. A satellite at a range of 2 000 km at the zenith would take at a maximum nine minutes to travel from 20
$^\circ$
elevation (the elevation cutoff) on one side of the sky to 20
$^\circ$
elevation on the other side of the sky. The 15 min window ensures overlapping bins are not necessary. The total number of candidates in each bin (N) is stored with the individual candidates. If there are candidates for more than one satellite at the same location at the same time, the candidate with the most number of detections across the pass is kept, with the others being discarded.
Fig. 4 shows an example of this. In this figure, satellite ‘y’ transmits 20 times as it passes over the visible horizon, and the uncertainty in measuring those positions in the data is shown by the jitter of the points around the predicted path of the satellite. Satellite ‘x’ has a similar path over the sky, and assuming the two satellites are passing with the same angular speeds at the same time, the predicted path of satellite ‘x’ overlaps the six measured positions in the shaded region at the same time at satellite ‘y’. This means that
$N = 20$
for satellite ‘y’ candidates and
$N = 6$
for satellite ‘x’ candidates. It is important that the satellite ‘x’ candidates are discarded as they are misidentifications, and none of the filters so far will remove them. In practice, if a satellite transmits hundreds of times during a pass, the number of additional candidates due to other satellite passes crossing over these can be quite large.
Figure 4. An illustration of a scenario where the predicted paths of two satellites ‘x’ (green) and ‘y’ (red) are similar over the sky. The blue points are the measured positions of satellite y in 20 separate images. The six measured positions in the shaded region are flagged as candidates for both satellites.
As another constraint, if
$N \lt 4$
, then these candidates would be discarded, as the algorithm performs best when N is high. There was no evidence that certain satellites transmitted less frequently than this cadence over the course of a pass.Footnote
j
The predicted and measured positions of a satellite may follow a similar direction over the sky, but be offset in time, suggesting that these candidates need to be discarded. This was calculated by comparing the rate of change of the angle of travel of the satellite relative to the observer as a function of time, and the candidates were discarded if this exceeded a threshold.
It is therefore assumed that if the measured positions of a set of candidates have
$N \gt 4$
, have a similar path over the sky to its predicted positions, have a similar time incidence to its predicted positions, and have more candidates across the pass than any other coincident satellite, that these candidates are classed as identifications. These identifications are stored in a database with all accompanying metadata and the analysis on these is presented in Section 3.
2.4. Estimating a misidentification rate
Due to the enormous number of images and satellite trajectories in this analysis, it is inevitable that misidentifications are made. A misidentification is defined as either incorrectly labelling RFI as a satellite, or when one satellite is incorrectly identified as another. Estimating this percentage is important, so a simplified simulation was performed to better understand the occurrence of these misidentifications.
For this simulation, synthetic satellite signals were injected into the 229.7 MHz dataset described in Table 1. This dataset was chosen because it already contained known transmission from a number of satellites and was an adequate length of 23 h (
${\approx}$
42 000 images). This simulation ran completely independently from the results reported for the dataset in Section 3.
The parameters of the simulation were as follows:
-
• 100 unique satellites (of the 18 540 with TLEs available under the constraints in Section 2.3.2) were randomly chosen and checked to make sure they were not known transmitters;
-
• Each satellite was assigned a random flux density to be injected between 120 and 500 Jy/beam
$^{-1}$
(this was not attenuated as a function of range or elevation); and -
• An injection pattern was chosen for each satellite (only injected when a satellite’s elevation was
$ \gt 20^\circ$
):
-
- every image;
-
- every three consecutive images then skipping the next three (and repeated); and
-
- for one image then skipping the next four (and repeated).
Once these signals had been injected into the dataset, both the detection and identification algorithms from Sections 2.3.2 and 2.3.3 were run over the data. Only the results for these 100 randomly selected satellites are reported here.
In total, 97 of the 100 randomly selected satellites reached an elevation of
$ \gt 20^\circ$
within the dataset. For the XX polarisation, 219 of the 232 injected passes were recovered with a single misidentification. For the YY polarisation, 213 of the 232 injected passes were recovered, also with a single misidentification. Both misidentifications occurred on the same pass of the same satellite and were due to three STARLINK satellites having extremely similar predicted positions.
This simulation therefore estimates that the percentage of misidentifications is likely to be
$ \lt 1\%$
, with a
${\approx}$
93% recovery rate in the number of passes of satellites. Note that misidentifications are more likely at lower elevation satellite passes due to the slant orthographic projection used in images which has more degrees per pixel at lower elevations.
Table 3. Detections per dataset
2.5. Frequency differencing experiment
Time differencing, although a successful technique for increasing the signal-to-noise ratio (SNR) of detections in images, has fundamental drawbacks. Image three in Fig. 3 shows the ISS in the current image (white) and the subtracted ISS from the previous frame (black). If the time between the subtracted images is reduced, these will begin to overlap at lower elevations where the angular speed of the satellite on the sky is slower. If the time between the subtracted images is increased, the sky moves more in this time and introduces more noise into the difference image. There are also more complex effects from photometric centroid displacement which are explained in Appendix A.
Satellite and FM radio transmissions are narrow band in nature, compared to astrophysical sources of radio emission. If the satellite emission can be isolated in one frequency channel, this can be subtracted from a nearby frequency channel at the same time step which does not contain the satellite’s emission. As the sky and the satellite now do not move, in theory, the SNR should be higher.
Frequency differencing was explored as an alternate method to time differencing. In this experiment, to ensure a fair comparison, both the time and frequency difference were done using fine channel (24.8 kHz) data. The frequency difference was performed by choosing a single fine channel where the satellite was visible, and subtracting from it another where the satellite was not visible from it for the same time step. The time difference was performed by choosing a single fine channel where the satellite is visible, and subtracting an image taken 40 s previously from it. This is illustrated in Fig. 3.
Two scenarios were considered. The first was using the ISS with reflection of terrestrial FM transmission. In theory, the ISS should reflect transmission from the FM transmitter, and be visible at frequencies within its
${\approx}$
0.1 MHz bandwidth. An image at the same time step, but at a frequency outside of the transmitter’s bandwidth, can then be subtracted from this to do the frequency differencing.
The second scenario was using SOLRAD 7B as a known transmitter. From preliminary observations it was known that the emission received from SOLRAD 7B showed periodicity at 136.82 MHz over a narrow bandwidth which would be a suitable case for testing the frequency differencing approach for a transmitter. After quality control was conducted on the dataset, it was found that ORBCOMM FM 110 (NORAD 41182) was transmitting more brightly than SOLRAD 7B. The transmission from these two satellites overlapped, but ORBCOMM FM 110 has the highest received flux density at 136.85 MHz. This satellite was also included in the analysis.
To calculate the SNR of the satellite in each image, the fitted amplitude (A) in Equation (1) is divided by an estimate of the noise of the image. To calculate the noise, the standard deviation (
$\sigma$
) of all pixels in the difference image is calculated after running two passes of removing all pixel intensities which are
$ \gt 3\sigma$
.
Each time step where an identification of the targeted satellite was made in both the time and frequency differenced images, the SNR of the identification was calculated. The ratio of the frequency difference SNR with the time difference SNR was averaged across all time steps and is presented in Section 3.5. This aimed to give an average ratio of the two SNRs across the entire pass. Note that this is tested on the XX polarisation images only.
3. Results
In total, 152 unique satellites were identified across these datasets. These results are displayed per dataset in Table 3 and per NORAD I.D. in Appendix C. This comprised 14,165 and 14,840 individual identifications in the XX and YY polarisations, respectively. The closest identification of a satellite, from the distance predicted by the TLE, was BRICSAT (NORAD 44355) at 316 km, and the furthest was COSMOS 2453 (NORAD 35500) at 3,013 km. There was also the identification of an ex-geosynchronous satellite EUTE 1-F5 (NORAD 19331) which was identified manually and is the longest range identification of a satellite using these techniques.
Fig. 5 shows the spatial spread of identifications over the sky. There is uniform coverage over the full sky, with identifications consistently being made down to the elevation cut-off of 20
$^{\circ}$
, showing that the identification algorithms are not biased to certain areas of the sky.
Figure 5. A plot of the azimuth/elevation locations of all 29 005 individual identifications of satellites across all 18 datasets analysed in this work. The colourbar shows the logarithmic flux density of the identifications in Jy/beam. The black circle surrounding the detections shows the 0
$^{\circ}$
elevation isoline.
The flux density of all identifications as a function of range is shown in Fig. 6.
Figure 6. The flux density and range of all 14 165 XX and 14 840 YY polarisation identifications of satellites across all 18 datasets analysed in this work.
The large majority of identifications were due to the direct detection of transmission from satellites. There was evidence of both IEMR and UEMR. Due to the high volume of identifications, only the most significant findings will be discussed, with the full results available in Appendix C. This Section summarises the results per dataset, and the implications of these results are discussed in Section 4.
3.1. Known downlink frequencies
The known downlink frequencies from this survey include the 146.1, 136.7, and 137.5 MHz frequency bands. All of these displayed a high number of identifications. Due to the widespread occurrence of high intensity signals, conducting radio astronomy at these frequencies would be exceedingly difficult.
3.1.1. 137.5 MHz
The variations in this dataset were extreme compared to other datasets, with some image RMSs as low as 10 Jy/beam
$^{-1}$
and as high as 10
$^6$
Jy/beam
$^{-1}$
. There are very short periods where no satellites are transmitting overhead and the Milky Way is clearly visible, but the majority of the time there could be one to 10 satellites visible. It was quite common to find a satellite transmitting so strongly that it would obscure emission from other transmitting satellites. For example, the identification of the ISS was only possible because nothing else was transmitting overhead at the same time.
At 137.5 MHz, the majority of satellites identified were ORBCOMM and SPACEBEE. There were also many identifications of STARLINK satellites in this dataset but not discussed in this work as they were reported in Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023). Detected intensities were frequently on the order of 10
$^7$
Jy/beam
$^{-1}$
, with most satellites transmitting periodically.
The ex-geosynchronous satellite EUTE 1-F5 was identified serendipitously when manually assessing the quality of this dataset. This satellite was placed in a graveyard orbit in the year 2000. At a range of
${\approx}$
36 000 km, it moves very slowly across the sky, meaning that time differencing will cause the emission to be completely subtracted. This particular satellite was transmitting strongly enough that its small movement across the sky was seen in the difference images and extremely obvious in the undifferenced images. It appeared to be rotating slowly as the flux density of the signal was sinusoidal over a period of multiple minutes (this was difficult to constrain due to constant high flux density signals from other satellites). This is expected to be easily identifiable when frequency differencing is used in the future.
The identification of PEGASUS R/B (NORAD 22491) was interesting, as it was not anticipated that a rocket body would be transmitting, especially one launched in 1993. It is unlikely to be reflecting terrestrial transmission as it was detected at a TLE predicted range of
${\approx}$
1 500 km, slightly further than any previous detection of the ISS (the largest satellite in orbit) using this technique (Grigg et al. Reference Grigg, Tingay, Sokolowski and Wayth2022). It was manually verified for quality control and was the only TLE that matched up well for the period it was visible. space-track.org and n2yo.com both list the name of this satellite as PEGASUS R/B, but an amateur observational database db.satnogs.org lists the name as OXP1. On another amateur satellite website,Footnote
k
there is a note: ‘An unknown emitter on 137.05 MHz, initially reported by Raydel, has been identified by experts on the Hearsat list as OXP-1 that was launched in 1993. The search found that the closest match was for the Pegasus R/B but a comparison of radar cross Section suggests the ‘correct assignments are 22491 for OXP-1 and 22489 for the Pegasus R/B stage’.
Fifty four unique SPACEBEE satellites were identified in this dataset. Their received emissions were sporadic and consistently exceeded
$10^{7}$
Jy/beam
$^{-1}$
. With a longer acquisition time, there would likely be additional SPACEBEE satellites, as there are missing SPACEBEE satellites from the same launch. Identifications of SPACEBEE 63 and SPACEBEE 55 (NORADs 47442 and 47459, respectively) were interesting because they were launched two years earlier than most of the other identified SPACEBEE satellites. The reason why only two were identified when there were still others in orbit from the same launch is unknown.
3.1.2. 136.7 MHz
In the short 12 min dataset at 136.7 MHz, which was targeted on SOLRAD 7B, there were also many identifications of the satellite ORBCOMM 110. This satellite was also identified in the 137.5 MHz dataset. It is listed as transmitting at 137.288 MHz, so the mechanism behind the propagation of this transmission is unknown at this time. If this dataset were longer, there would have been more identifications made, as there were also many short bursts of radio energy which were difficult to identify without seeing the full passes of these satellites.
Fig. 7 shows that the intensities of SOLRAD 7B in the two polarisations varied sinusoidally in time, out-of-phase by 90
$^{\circ}$
. This suggests that the satellite may no longer have stabilisation control and could be tumbling. More follow up observations would be needed to confirm this.
Figure 7. The flux density of detections of SOLRAD 7B in the 136.7 MHz dataset. ORBCOMM FM110 is also transmitting strongly between 0 to 150 s.
3.1.3. 146.1 MHz
The quality of these data were high, with satellites being the main source of RFI. The majority of satellites identified at this frequency are cubesats or for amateur use, with a radar cross Section of
$ \lt $
1 m
$^{2}$
. All of the identified satellites at this frequency have their designated downlink frequency within the 0.926 MHz bandwidth of the acquisition. The cadence and measured flux density of the transmissions varied greatly, proving the robustness of the identification algorithms. Many of the cubesats had their transmission powers publicly available and all were less than 0.5 W in power.
BGUSAT, MAX VALIER SAT, and TRITON 1 are known to transmit broadband UEMR (Grigg et al. Reference Grigg, Tingay, Sokolowski and Wayth2022), and at this frequency they are all observed at their operational downlink frequency. The flux density of the received transmission is much higher at this frequency than any other frequency they have been observed at in this study and in Grigg et al. (Reference Grigg, Tingay, Sokolowski and Wayth2022). This could imply that the transmission detected at other frequencies is from a different mechanism and is likely UEMR. Unlike the UEMR from Starlink satellites (Grigg et al. Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023), the transmission is visible for short multi-second bursts, and then disappears for short periods of time. For example, in general across all frequencies TRITON 1 disappears for longer periods than BGUSAT and MAX VALIER SAT.
OSCAR 11 (NORAD 14781) was launched in 1984 and is also an example of a long decommissioned satellite still emitting at radio frequencies. This satellite was identified over three passes, and these were all times when the satellite was illuminated by the Sun, powering its battery. These were the only passes above 20
$^{\circ}$
elevation during the span of the dataset.
HAMSAT (NORAD 28650) is yet another example of a satellite which has been decommissioned and yet is still transmitting. Its official websiteFootnote
l
states the reason for it being decommissioned as a ‘failure of on-board lithium batteries’. There were two passes of the satellite during the dataset at 46
$^{\circ}$
elevation, where it was illuminated by sunlight, and 34
$^{\circ}$
elevation, where it was eclipsed by the Earth. It was identified in the first pass in the YY polarisation and was not identified in the second pass, implying that the Sun may temporarily power the satellite, causing it to transmit.
3.2. Protected frequencies
The protected frequencies were 73.4, 150.8, 324.2, and 325 MHz. No satellites were identified in any of these datasets.
3.2.1. 324.2 & 325 MHz
The two datasets at 324.2 MHz and the other at 325 MHz were of high quality and had extremely low levels of RFI if any. The only high flux density signals in the datasets came from the Sun and its sidelobes.
3.2.2. 150.8 MHz
There were two datasets acquired at 150.8 MHz. In the 150.8 MHz dataset acquired on 2023-06-13, there were two instances where an object appeared to move over the sky. These both followed the north/south common aircraft flight path over the telescope which is illustrated in Figure 8 in the paper by Grigg et al. (Reference Grigg, Tingay, Sokolowski and Wayth2022). There are some local mobile/defence transmitters across WA which have powers up to 50 W which could be causing these weak reflections. Follow up observations could be made to discern which transmitter the reflected transmission originated from, but is out of scope for this work.
Both datasets had 32, 28.9 kHz fine channels of data available, and although the main analysis processed these combined into a single 0.926 MHz channel, the fine channels were manually surveyed to determine if there was any narrowband emission in the data. It was noticed that there were some occasional broadband (covering all frequency channels) bursts, while others were narrow band (one or two fine frequency channels). There was no apparent temporal association between any of these bursts.
3.2.3. 73.4 MHz
Lastly, the 73.4 MHz dataset exhibited the least amount of RFI out of all the datasets in this work. There was no visible RFI in the 37 min of the dataset, which makes it ideal for other radio astronomy analyses.
3.3. Comparison frequencies
3.3.1. 229.7 MHz
The 229.7 MHz dataset overlaps with the Australian digital television broadcasting band which extends up to 230 MHz. The Australian Radio Frequency Spectrum PlanFootnote m lists the uses of this frequency for ‘broadcasting’ (fixed and mobile) and aeronautical radionavigation. Overall the quality of this dataset was high, and the highest measured flux densities were from the Sun and Milky Way.
All identified satellites at this frequency were launched by the Commonwealth of Independent States (formerly USSR) and were of the COSMOS, GONETS, and STRELA series. All of these detected satellites displayed a very predictable pattern of transmitting within a single image, every 60s. This is why the identification counts are low for each satellite in Table C1. An amateur satellite websiteFootnote n lists the reason for the transmission pattern as: ‘Message Signal every 60s, silent when receiving data message from ground’.
Sokolowski et al. (Reference Sokolowski2021) shows identifications at the same frequency a year and a half earlier. All of these identifications were also made in this work, apart from COSMOS 2438 (NORAD 32955). The EDA2 is most sensitive towards the zenith, and as the satellite was not detected during two high elevation passes, means it may have stopped transmitting during this year and a half period. It would be unlikely to be silent due to receiving data from the ground as it was travelling over Australia at this time. Our work additionally identifies some of the GONETS M series which are not identified in the work of Sokolowski et al. (Reference Sokolowski2021) (GONETS M17 and 22 were launched after the Sokolowski et al. Reference Sokolowski2021 dataset).
Twenty three aircraft were also detected in this dataset at low intensities. 21 of these were on the common north/south flight path shown in Figure 6 in Grigg et al. (Reference Grigg, Tingay, Sokolowski and Wayth2022). They were mostly visible in the southern half of the sky, implying they are reflecting transmission closer to Perth (
${\approx}$
600 km south) where the aircraft would be taking off/landing. There was no evidence of any satellites reflecting this same terrestrial transmission.
3.3.2. 159.4 MHz
There were three datasets used at this frequency. One was the
${\approx}$
23 h dataset acquired for this work on 2021-11-16 and the others were two
${\approx}$
133 h datasets acquired by Sokolowski et al. (Reference Sokolowski2021) on 2020-06-26 with the AAVS2 and EDA2 acquiring data simultaneously.
The 2021-11-16 dataset yielded identifications of the satellites BGUSAT and TRITON 1, outside of these satellites’ designated downlink frequency. This is another dataset showing detections of these two producing UEMR.
It was noticed in the study by Sokolowski et al. (Reference Sokolowski2021) that there were many additional identifications of rocket bodies, debris objects, and geostationary orbit satellites. As none of these were identified in the 2021-11-16 dataset, there was the question of why these were being identified in the 2020-04-10 and 2020-06-26 datasets in the Sokolowski et al. (Reference Sokolowski2021) work. The authors of Sokolowski et al. were kind enough to allow the two datasets from 2021-11-16 with the AAVS2 and EDA2 acquiring simultaneously for
${\approx}$
133 h to be made available for reprocessing with our new algorithms.
In our reprocessing of the data, for the polarisations XX/YY, BGUSAT was identified 316/486 times in the AAVS2 dataset and 469/515 times in the EDA2 dataset. Sokolowski et al. (Reference Sokolowski2021) combines both polarisations into I (XX + YY) and does not differentiate between the datasets, but makes 499 identifications of BGUSAT. The timestamps for the two telescopes are on a very slightly different cadence, which makes combining the identification counts between the two datasets difficult, but there are more identifications of BGUSAT using our new algorithms.
Doing the same for TRITON 1, in our reprocessing of the data, for the polarisations XX/YY, TRITON 1 was identified in 158/180 images in the AAVS2 dataset and 240/237 images in the EDA2 dataset. Sokolowski et al. (Reference Sokolowski2021) makes 191 identifications of TRITON 1. Again there are more identifications of the satellite using our new algorithms.
The only satellite which is identified less with the new algorithm is MAX VALIER SAT. There were six individual identifications only in the YY polarisation in the EDA2 dataset with our algorithm, whereas there were 17 individual identifications in Sokolowski et al. (Reference Sokolowski2021).
In the EDA2 dataset, SWIATOWID (NORADs 117) was additionally identified with our new algorithms.
There were no rocket bodies, debris objects, or geostationary orbit satellites identified with these new algorithms. This will be discussed in Section 4.
3.4. Additional frequencies
3.4.1. 160.2 MHz
The 160.2 MHz dataset was recorded over 55 h. There were 29 unique satellites identified over the two polarisations. Most of the identifications were Starlink satellites. The character of the emission is similar to the character reported at 159.4 MHz in the paper by Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023). This emission is constant from image-to-image, unlike the transmission every 100 s that was exhibited at 137.5 MHz in the same study. The study also suggests that the emission is likely caused by the propulsion or avionics system and is also likely to be broadband, meaning that it was classed as UEMR.
This case is different. In the Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023) study, the identified Starlink satellites at 159.4 MHz are type ‘v1p5’, whereas those identified in this dataset are all ‘v2_mini’. All of these ‘v2_mini’ Starlink satellites were on one of four payloads which launched on 2023-02-27, 2023-04-19, 2023-05-19, and 2023-06-04. Our dataset was acquired on 2023-06-23 for just over two days. This means that the Starlink satellites that were launched on 2023-02-27 and 2023-04-19 were no longer being boosted to their operational altitude. Two examples can be seen for NORAD I.D.s 55695Footnote o and 56301Footnote p showing their historic altitudes. The propulsion system may still be engaged when these satellites are identified, but implies that this UEMR from Starlink satellites is likely to be ongoing when they arrive at their operational altitude.
BLUEWALKER 3 (NORAD 53807) was also identified at this frequency. In its technical specifications document submitted to the Federal Communications Commission, there is nothing mentioned about communications lower than 750 MHz.Footnote q There was no evidence of reflection of terrestrial transmission off aircraft at this frequency, putting in question whether the detection of BLUEWALKER 3 was reflected energy. Its highest recorded flux density was in the YY polarisation at 1,043 Jy/beam. The source of this signal is currently unknown.
Table 4. Frequency differencing experiment results.
The last two objects which were TRITON 1 and OBJECT DFootnote r . TRITON 1 is a known broadband transmitter across many frequencies (Grigg et al. Reference Grigg, Tingay, Sokolowski and Wayth2022), but OBJECT D (either unnamed or unidentified) has not been previously identified, and very little information is available on it other than it was in a launch of other similarly named satellites from China on 2018-12-07.
3.4.2. 185.2 MHz
The quality of this dataset was high in general. This being said, there were two transmitters on the southern horizon which were visible for the second half of the dataset (night time local time). Measured as a clockwise bearing angle from the vertical (north) in the images, the two transmitters are on a bearing of 182
$^{\circ}$
and 195
$^{\circ}$
. Taking the transmitters from Table 2, the bearing angle on the images of the Morowa transmitter is 195
$^{\circ}$
(
${\approx}$
300 km away) and the bearing of the two Perth transmitters is 185
$^{\circ}$
(
${\approx}$
600 km away). These line up well with the transmitters seen on the horizon in the images. The effect of the two Perth transmitters interfering with each other can also be seen in the images.
The only satellite identified in this dataset was BGUSAT. There were also eleven aircraft detected, with most on a north/south or east/west trajectory. This dataset was acquired in February 2020, which was around the time that strict restrictions were put on International flights in Western Australia, which could explain the low detection rate of aircraft considering the dataset was 40 h in length. Although the two Perth transmitters are 50 kW in power, comparable to transmitters in the FM band, the power is spread over a larger bandwidth (100 kHz for FM vs.
${\approx}7$
MHz for DTV). This could explain the lack of satellite identifications from the reflection of this DTV transmission.
3.5. Frequency differencing
Table 4 summarises the results from the four datasets which were used to test the frequency differencing technique.
3.5.1. Reflected FM transmission
The three passes of the ISS at the two different frequencies showed between a 39% increase to an 11% decrease in SNR of the ISS. At these FM frequencies, there were no astronomical sources visible in the images, which meant that the noise was dominated by sidelobes and the FM transmitters visible on the horizon. In some cases in the time difference, the direct transmission from these FM transmitters was partially subtracted as it was visible in both images used in the differencing. For frequency differencing, the direct transmission would never be subtracted as it was only visible in one of the images used in differencing by design. This explains why there was a 11% decrease in the SNR of the ISS in the frequency differenced images.
This being said, frequency differencing increases the average SNR of detections by 39% and 15% in the other two datasets. This was more due to a reduction in the noise than an increase in the signal used to calculate the SNR.
3.5.2. Transmission
For the pass of SOLRAD 7B, frequency differencing increased the average SNR by 123%, and the pass of ORBCOMM FM 110 by 55%. The reason for such a different result for both satellites is that the intensities of SOLRAD 7B oscillate as shown in Fig. 7. This means that SOLRAD 7B’s measured flux density will be low in some images, but the subtracted flux density from the image 40 s prior will be high, resulting in a low SNR in the time difference.
The lack of high flux density transmitters on the horizon, along with astronomical sources of radio emission being slightly visible in the images results in the frequency difference increasing the SNR of the two satellites. This will be the favoured approach for detecting direct transmission from satellites in the future.
4. Discussion
The results presented in this work demonstrate significant improvements in our algorithms for the detection and identification of satellites. Imaging the whole sky presents advantages in that multiple satellites can be characterised simultaneously in the same image, providing rich datasets worthy of many types of astrophysical analysis.
A total of 152 unique satellites were identified across frequencies over a 250 MHz range. The 18 datasets analysed varied considerably in the type of signals received by the telescope. Different sources of radio energy such as terrestrial transmitters, meteors, the Sun, astrophysical sources of broadband synchrotron emission, aircraft, lightning, and other anthropogenic sources were handled to robustly identify satellites. Different satellite transmitting cadences were handled by the algorithm to identify satellites which transmitted sporadically over small portions of their passes.
There is a slight difference in the total number of identifications in each polarisation, with 14 165 individual identifications in the XX polarisation and 14 840 in YY. This is more skewed in some individual datasets shown in Table 3. This can be influenced by the antenna beam response in each polarisation, as well as physical properties of the satellite itself (polarisation characteristic and orientation of the satellite’s antenna). For example, the 229.7 MHz dataset has
${\approx}$
43% more identifications in the XX polarisation. All of the satellites identified at 229.7 MHz followed a north/south path over the sky, meaning the satellite will have slower apparent angular speed in the north and south of the sky where the XX polarisation is more sensitive. This can explain the difference observed between the two polarisations. An analysis of how the antenna’s physical characteristics influences detection rates would be interesting but is out of scope for this work.
The majority of identifications were in LEO as shown in Table C1. For example, identifications at 229.7 MHz show that medium earth orbit (MEO) satellites can also routinely be identified. Since the majority of satellites are in LEO, being able to identify satellites beyond LEO shows good promise for future broader frequency surveys. The manual identification of the (now in graveyard orbit) geosynchronous satellite EUTE 1-F5 shows that satellites with a range of over 36 000 km can be detected. Tests show that it could repeatedly be detected in non-differenced images of the sky, showing that this technique could be applied to geosynchronous/geostationary range searches in the future.
Fig. 6 shows that the measured flux density of satellite transmission can regularly exceed 10
$^6$
Jy/beam
$^{-1}$
, with many higher than 10
$^8$
Jy/beam
$^{-1}$
. The 137.5 MHz frequency band is host to many of these satellites, but is also widespread among other frequency bands. This indicates challenges that will be faced in attempts to make astrophysical observations at these frequencies. Satellite constellations such as ORBCOMM and SPACEBEE are legally allowed to transmit at these frequencies, but their transmission obscures astrophysical signals over the whole sky in most instances.
4.1. Comparison of algorithms from Sokolowski et al. (Reference Sokolowski2021)
Reprocessing of the two five and a half-day datasets from the work of Sokolowski et al. (Reference Sokolowski2021) shows a marked decrease in the number of misidentifications. Our work introduces an additional final step in the identification algorithm that takes candidate detections and ensures they follow the expected trajectory of the satellite over the sky, which eliminates many of these misidentifications. This new algorithm also identifies more of the same objects, showing that they are robust over a long observing period. The simulation run to estimate a misidentification rate also shows an extremely low occurrence of these.
4.2. Incomplete database information
At the time of writing this paper, space-track.org tracks
${\approx}$
45 000 objects in Earth orbit. The maintenance and validation of this database is an enormous effort and is in-part reliant on satellite operators being able to identify their own spacecraft. Our detection of PEGASUS R/B appearing to ‘transmit’, along with the information from the amateur satellite page showing that it was likely OXP-1, is an example of inconsistencies that could appear in this database. There were also at least ten passes of objects that travelled across the sky at angular speeds of LEO satellites, which could not be attributed to anything in the TLE list obtained from space-track.org.
4.3. Satellite behaviour
The lack of identifications in any of the protected frequency bands is encouraging. The degree of protection over the 73.4, 150.8, 324.2, and 325.0 MHz frequencies varies for each, but it is promising that in these datasets there are no identifications made. These results can serve as a baseline for future studies to ensure satellites are transmitting safely and correctly.
4.4. Satellite operational status
A number of satellites were identified which have been officially decommissioned, and yet are still emitting at radio frequencies. There could be a couple of reasons for this. Although decommissioned, these satellites may still be performing a limited subset of functions, or loss of control of the satellite may mean that they will keep transmitting until they lose power or ultimately decay. For the three listed examples of SOLRAD 7B, HAMSAT, and OSCAR 11, they all display the following characteristics: they were all launched some time ago (1965, 2005, and 1984, respectively), are publicly listed as decommissioned, all are illuminated by the Sun when detected in the data (the Sun potentially temporarily charges up the battery), and they all show sinusoidal variation in their lightcurves (potentially due to lack of stabilisation control).
Satellites that remain active after decommissioning need to continue to be tracked and monitored but are shown to sometimes continue to take up an allocation in an already congested radio allocation zone. The legislative landscape and permit allocation procedures have changed significantly since some of these satellites were launched and will continue to change in the future. There needs to be stronger enforcement of procedures for the decommissioning of satellites so they do not continue to be a source of RFI for radio astronomy.
4.5. UEMR
In the work by Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023), the authors quote a discussion with SpaceX (the company operating the STARLINK constellation), in which SpaceX engineers propose that the propulsion or avionics system is the mechanism behind the likely broadband UEMR detected at 159.4 MHz. In this particular example, the STARLINK satellites were being boosted to their operational altitude. Our new results show that the STARLINK satellites are also transmitting at 160.2 MHz when they are already at their operational altitude. Multiple passes of the same STARLINK satellites are also identified, bringing into question whether the propulsion system is the cause of the observed signal.
None of the model ‘v1p5’ Starlink satellites which were detected at 159.4 MHz in Grigg et al. (Reference Grigg, Tingay, Sokolowski, Wayth, Indermuehle and Prabu2023) were detected at 160.2 MHz where we detect the ‘v2_mini’ model Starlink satellites in this work. For the ‘v2_mini’ Starlink satellites which were identified above 85
$^{\circ}$
elevation, the average flux density was was 78 Jy/beam
$^{-1}$
at a range of 532 km (EIRP
$\approx$
2.6 μW). The work by Bassa et al. (Reference Bassa, Di Vruno, Winkel, Józsa, Brentjens and Zhang2024) confirms that the types of UEMR transmitted by the version 1 and version 2 are at different frequencies, and that the version 2 (mini and mini direct-to-cell) satellites were detected on average 32 times brighter than the version 1.0 satellites. Further monitoring of the STARLINK satellites will be crucial for assessing the impact of their UEMR on radio astronomy.
In addition to STARLINK, there were other satellites also transmitting broadband UEMR. BGUSAT, MAX VALIER SAT, and TRITON 1 continue to show evidence of being detected in many frequency channels. At the date of writing this paper, BGUSAT has decayed, but MAX VALIER SAT and TRITON 1 will continue to be sources of interference for radio astronomy research. BLUEWALKER 3 also shows evidence of an unknown source of radio frequency transmission at 160.2 MHz, which is not listed as a downlink frequency for the satellite. An in-depth investigation of whether satellites are transmitting outside of their designated downlink frequency bands is beyond the scope for this work, but Table C1 provides detailed information on the identification of individual satellites that interested readers may use to assess this further.
4.6. Aircraft
Aircraft prove to be an ultra-bright source of RFI at frequencies of high powered terrestrial transmitters, and this has been characterised for the FM band in Grigg et al. (Reference Grigg, Tingay, Sokolowski and Wayth2022). The datasets at frequencies of 185.2 and 229.7 both show evidence of aircraft, although at lower intensities than in the FM band. This again showcases the potential to use bi-static radar in the form of a radio telescope to passively search for and locate aircraft. With larger baselines, it would be possible to perform accurate parallax measurements of aircraft to also determine range.
4.7. Frequency differencing
Testing of the frequency differencing method has shown that in general, the SNR of satellite detections is increased compared to time differencing. The caveat to this is that time differencing may give a better SNR when there are high measured flux densities from terrestrial transmitters visible on the horizon. Frequency differencing additionally minimises the impact from broadband sources of radio emission, and removes the subtracted satellite emission shown in image 3 in Fig. 3. Ultimately, frequency differencing will be used in future all-sky surveys for the characterisation of satellites with the EDA2.
4.8. Future work and improvements
These algorithms are now ready to be deployed in an automated manner over a more targeted survey on key science frequencies of the SKA-Low observing bandwidth of 50–350 MHz. With the addition of 32 channels of 28.9 kHz bandwidth data, such a study would provide very detailed information about the implications of satellites on radio astronomy observations at these frequencies. This will be the next work in this series and will form an enhanced baseline for satellite activity, which can be repeated in the future to assess the change in radio frequency environment.
Although the predicted misidentification rate is
$ \lt 1\%$
, this may not hold for datasets with a high number of visible satellites or high levels of RFI. As frequency differencing will be the favoured algorithm for processing the next larger survey, a more comprehensive simulation will be performed in that work.
Time differencing still holds value, for if a satellite produces broadband UEMR, then it is unlikely to be detected with frequency differencing. A potential improvement could be to decrease the time between the two images used in the time difference for higher frequencies as the size of the PSF in degrees decreases at higher frequencies.
Future work will likely attempt to expand on the polarisation analysis, looking to compare the observed polarisation response to physical antenna characteristics. Using far longer integration times to do time differencing may prove useful for detecting satellites in orbits out to MEO and geostationary. Building a database of unidentifiable detections will also be useful for future orbit determination analysis. TLEs could then be derived for these objects.
Although transmission bandwidth of some satellites is known, currently it is difficult to determine if this leaks into adjacent smaller frequency channels when the current observing frequency bandwidth is 0.926 MHz. There is also time variability of the received transmission which is difficult to model. The future survey will have finer frequency resolution across multiple channels which will enable a more detailed study of transmission bandwidths.
Lastly, the identifications made in this study have been stored in a database for future use. Future identifications will be added to this database, which may reveal temporal trends over the coming years. This can be used in conjunction with conventional approaches to space situational awareness such as radar and optical telescopes.
5. Conclusion
Using all-sky imaging to comprehensively characterise and identify satellites across the SKA-Low science frequencies with two orthogonal polarisations has made critical insights into the behaviour of satellites at these low frequencies. This study offers evidence to the challenges that future low frequency science will face from the exponentially growing population of satellites in Earth orbit. Examples are also shown of satellites transmitting UEMR and so-called ‘zombie satellites’ transmitting after their decommissioning or demise. The identification algorithms put forward in this work are ready to be deployed in a systematic and targeted survey of the SKA-Low observing frequencies to assess the widespread effects of satellite transmission. LEO and MEO automated detection and identification is currently achievable, with a detection in a graveyard geosynchronous orbit paving the way for future algorithmic development. Ultimately, we will develop comprehensive tools for space situational awareness, offering unique scientific insight into the behaviour of satellite operations.
Acknowledgement
We would like to acknowledge the Wajarri Yamatji people as the traditional custodians of the Murchison Radio Observatory site, Inyarrimanha Ilgardi Bundara (Lonsdale et al. Reference Lonsdale2009). The data processing for this work was run at the DUG Technology Perth HPC centre. We also wish to thank the developers of several key packages to this work, namely MIRIAD (Sault et al. Reference Sault, Teuben, Wright, Shaw, Payne and Hayes1995) for the preprocessing, and the Python packages Skyfield (Rhodes Reference Rhodes2019), Astropy (Astropy Collaboration et al. 2018), and Numpy (Walt, Colbert, & Varoquaux Reference van der Walt, Colbert and Varoquaux2011).
Data availability
Not applicable.
Appendix A. Sinusoidal motion in measured positions
This Appendix investigates a phenomenon, where the offset between the predicted and measured locations of a satellite varied sinusoidally in time for some cases. This effect was strongest at low elevations. An example of this can be seen in the left panel in Fig. A1, as it also manifests in the values of
$\theta_{P}$
and
$\theta_{M}$
.
Figure A1. Two scatter plots showing a real example of detections of ORBCOMM FM113 in the 137.5 MHz dataset (left) and a simulated pass of ORBCOMM FM113 (right). The scatter plots show
$\theta_{P}$
(red points) and
$\theta_{M}$
(blue points). The green points are from the linear line of best fit for the
$\theta_{M}$
values. The residuals are between the
$\theta_{M}$
and
$\theta_{P}$
values. The image number is zeroed on the image of the first candidate.
There are two effects here which are important to separate. The first is due to photometric centroid shift. Image three in Fig. 3 shows an ideal separation between the positive and negative intensities of the satellite due to time differencing. When the separation between these two decreases, it can cause an asymmetry in the intensity profile of the positive signal, which results in a shift in its centroid. The closer these two are, the greater this observed shift is.
The second effect is due to sub-pixel sampling of the fitted centroid. At frequencies of 73.4 and 325 MHz, the sampling of pixels is 2.2 and 0.5 degrees/pixel, respectively. At lower frequencies, it can take a satellite several images to fully cross a pixel. This sub-pixel sampling imparts a sinusoidal pattern on the offsets between the predicted and measured positions over time.
To test these effects, a simulation was run. A pass of the satellite ORBCOMM FM113 (NORAD 41185) was chosen, as its plotted values for
$\theta_{M}$
showed this sinusoidal variation as can be seen in the left panel in Fig. A1. To recreate this as simply as possible, the images were replaced with arrays of zeros, and a synthetic signal was introduced at the predicted locations of the satellite. This included a positive signal at the location in the current image, and a negative signal at the location of the satellite in the subtracted image used in the time differencing. The injected signal was a circular Gaussian multiplied with a Bessel function as an approximation of the point spread function. This was parameterised as follows:
\begin{align} I_{\text{model}} = J_0\left(\frac{r}{\sigma_{\text{B}}}\right) \cdot A \exp\left(-\frac{1}{2 \sigma_{\text{G}}^{2}} \left[(x - x_{0})^{2} + (y - y_{0})^{2}\right]\right), \\[-12pt] \nonumber \end{align}
where
$r = \sqrt{(x - x_0)^2 + (y - y_0)^2}$
is the radial distance from the centre
$(x_0, y_0)$
, A is the amplitude,
$\sigma_{\text{B}}$
is the standard deviation of the Bessel function,
$\sigma_{\text{G}}$
is the standard deviation of the Gaussian function and
$J_0(x)$
is a zero-order Bessel function of the first kind:
\begin{align} J_0(a) = \sum_{n=0}^{\infty} \frac{(-1)^n a^{2n}}{2^{2n}(n!)^2}.\end{align}
The following values estimated from the real measurements of the pass of ORBCOMM FM113:
$A = 50,000,000$
,
$(x_0, y_0)$
were the TLE predicted positions of the satellite,
$\sigma_{\text{B}}$
was derived from the size of the PSF in degrees for this frequency, and
$\sigma_{\text{G}} = 4$
(appropriate taper for this frequency).
The right panel in Fig. A1 shows that the simulation displays the same sinusoidal pattern at the beginning and end of the pass of the satellite as the real case on the left. This effect can therefore be explained with geometry and is not related to fluctuations in the position of the received signal.
The photometric centroid shift is minimised by choosing a large enough time offset between the two images being used in the time difference, and imposing an elevation cutoff helps when the satellite is moving slower over the sky at low elevations. The sub-pixel sampling effect is minimised by fitting the straight line (green line in Fig. A1) to the
$\theta_{P}$
and
$\theta_{M}$
values. This performs well in real examples.
In summary, understanding both the photometric centroid shift and sub-pixel sampling effects allows us to manage their effects to ensure that satellite identification remains robust and accurate.
Appendix B. Technical details for making identifications
This Appendix provides the quantitative, technical analysis of how identifications are made to accompany Section 2.3.3.
To repeat for context, we now refer to the TLE predicted position of a satellite as the ‘predicted’ position, and the position of the fitted Gaussian (
$x_{0}$
and
$y_{0}$
) as the ‘measured’ position. Each candidate will have this predicted and measured position information.
The instantaneous direction of travel of a satellite over the sky in an image is calculated as the bearing angle between two adjacent candidates and is denoted by
$\theta$
.
$\theta$
is the bearing angle measured clockwise between
$\hat{y}$
, a unit vector in the y direction, and
$\overrightarrow{C_{1}C_{2}}$
, the vector between the position of the candidate at time
$t_{1}$
(
$C_{1}$
) and the position of the next candidate at time
$t_{2}$
(
$C_{2}$
).
$\theta$
is calculated for both the predicted (
$\theta_{P}$
) and measured (
$\theta_{M}$
) positions for each candidate, resulting in
$N-1$
$\theta$
s for N candidates. This process is illustrated in Fig. B1.
Figure B1. An illustration showing how
$\theta_{M}$
and
$\theta_{P}$
are calculated for three candidate detections
$C_{1}$
,
$C_{2}$
, and
$C_{3}$
from three separate images. The circular outline shows the full visible sky horizon to horizon and the curved black dotted line is the predicted path of a satellite across the sky from TLE information. These candidates are for the same satellite, where the blue points mark the measured locations, and the red points mark the predicted locations. The distance between these is exaggerated in this figure.
Figure B2. A plot showing
$\theta_{P}$
(red points) and
$\theta_{M}$
(blue points) calculated for passes of ORBCOMM FM117 (left) and SPACEBEE-111 (right). The green points are from the linear line of best fit for the
$\theta_{M}$
values. The residuals are between
$\theta_{M}$
and
$\theta_{P}$
values. The image number is zeroed on the first image the candidate is detected.
$\phi$
is measured as the clockwise bearing angle of each of the lines of best fit.
The list of candidates is then sorted so that all candidates for each satellite are put into 15 min bins which are designed to be long enough that any object in the TLE list will pass over the sky completely during this time, and where each bin represents individual passes over the sky for each satellite. The number of candidates for the same satellite in each bin were summed (N) and stored with each candidate. If
$N \lt 4$
, then these candidates would be discarded, as the algorithm performs best when N is high. There was no evidence that certain satellites transmitted less frequently than this cadence over the course of a pass.Footnote
s
$\theta_{P}$
and
$\theta_{M}$
are then calculated for each group of candidates. Fig. B2 shows real examples of what this looks like for the pass of a satellite that transmits every consecutive image across its path (left) and a satellite that transmits periodically across its path (right).
$\theta_{P}$
and
$\theta_{M}$
gradually increase or decrease depending on the direction of travel of the satellite over the sky.
There is uncertainty in the measured positions of the satellite (due to noise and fitting) and uncertainty in the TLE predicted positions (described in detail for LEO satellites in Geul, Mooij, & Noomen Reference Geul, Mooij and Noomen2017). Although these will perturb the
$\theta_{P}$
and
$\theta_{M}$
values in Fig. B2, they will still both follow the same overall direction for a correct identification. Fitting a straight line to
$\theta_{P}$
and
$\theta_{M}$
provides a way to ‘smooth’ out these uncertainties and proves to be robust for low N. This being said, the algorithm performs best when N is high.
The bearing angle measured clockwise from a vertical unit vector in the
$\theta$
- time plane (shown in Fig. B2) with the fitted line across all
$\theta_{P}$
and
$\theta_{M}$
is calculated and referred to as
$\phi_{P}$
and
$\phi_{M}$
, respectively. The difference in these values,
$|\phi_{P} - \phi_{M}|$
, or
$\Delta\phi$
, shows the variation in travel direction over the sky between the measured and predicted positions. This value is stored with every candidate detection.
The next step was to reduce the number of candidates by ensuring that no candidates had the same location and time. The signals detected in the data cover at least the width of the PSF on the sky, or even longer if in LEO. This means that there are likely going to be multiple satellites over this PSF width, for which multiple may have satisfied the candidate criteria so far. Fig. 4 shows a simplified example of this.
Each candidate now has N and
$\Delta\phi$
calculated. If there was more than one candidate for the same time and location, only the candidate with the highest N was kept. If N was equal for two candidates, then the candidate with the lowest value for
$\Delta\phi$
was kept. Therefore, in Fig. 4, the six candidates for satellite ‘x’ would be discarded, because their N (six) is less than the N (twenty) of satellite ‘y’.
Once these discards had been made for the whole dataset, each candidate had N and
$\Delta\phi$
recalculated within each 15 min bin. Candidates with
$N \lt 4$
were again discarded, as well as candidates with
$\Delta\phi \gt 10$
.
The mean of the residuals between the fitted line to
$\theta_{M}$
and the
$\theta_{P}$
values was calculated and denoted as
$\mu_{r}$
. When the predicted trajectories of two satellites are very similar in position, but offset in time, this value of
$\mu_{r}$
will become larger. In Fig. B2, the fitted values for
$\theta_{P}$
and
$\theta_{M}$
will be offset by a fixed distance proportional to this timing offset. Candidates with
$\mu_{r} \gt 3^{\circ}$
are therefore discarded. This value was data driven and derived by observing examples where two satellites had very similar orbits. These identifications are then stored in the database.
An important caveat is that this algorithm was parameterised to prioritise making higher confidence identifications. Some of the parameters could have been relaxed, which would have made more identifications, especially in the noisier datasets, but at the expense of allowing more misidentifications into the results.
Appendix C. Satellite detection list by NORAD
This Appendix shows in detail the list of identifications in full in Table C1.
Table C1. Individual identifications per NORAD I.D. for each dataset
Open access
