Establishing the Capabilities of the Murchison Widefield Array as a Passive Radar for the Surveillance of Space

This paper describes the use of the Murchison Widefield Array, a low-frequency radio telescope at a radio-quiet Western Australian site, as a radar receiver forming part of a continent-spanning multistatic radar network for the surveillance of space. This paper details the system geometry employed, the orbit-specific radar signal processing, and the orbit determination algorithms necessary to ensure resident space objects are detected, tracked, and propagated. Finally, the paper includes the results processed after a short collection campaign utilising several FM radio transmitters across the country, up to a maximum baseline distance of over 2500 km. The results demonstrate the Murchison Widefield Array is able to provide widefield and persistent coverage of objects in low Earth orbit.


I. INTRODUCTION
R ADIO telescopes, sensitive instruments designed to detect the faintest signals from the most distant cosmic objects, are proving to be capable receivers for the purpose of radar. Although radio telescope arrays have been used as imaging interferometers, recent interest in phenomena such as pulsars, fast radio bursts, and cosmic rays, is motivating the development of high time-resolution capabilities for interferometric arrays [1]- [3]. These high time-resolution observation modes allow the direct implementation of radar functionality. Thus, researchers are investigating the use of arrays of radio telescopes as radar receivers for the surveillance of space [4]- [7]. This is certainly true for the Murchison Widefield Array (MWA), a low-frequency radio telescope in remote Western Australia, a precursor to the forthcoming Square Kilometre Array (SKA) [8]. Even before the MWA commenced operations in 2013, its use as a passive radar for the surveillance of space has been considered and planned [9]. Over the last decade, work has been undertaken to investigate the MWA's ability to detect and track resident space objects (RSOs) in low Earth orbit (LEO) from the reflections of terrestrial transmissions, particularly FM radio [4], [10], [11]. Initial simulations predicted the MWA would be able to detect FM radio reflections from an RSO of a radar cross-section that is as small as 0.5 m 2 at a range of 1000 km [10].
Rather than relying upon pre-existing sources of illumination, many radio telescopes are being utilised as active radar receivers as part of a bistatic system in conjunction with a cooperative radar transmitter [6], [12]. This configuration is beneficial, as the radar transmitter is ideal for the purpose of providing high-powered illumination of the RSOs. However, this is often achieved through the use of narrow-beam transmitter antennas. A direct beam will illuminate only a narrow volume and so the system will only be capable of surveilling a small region, especially when combined with a narrow receiver beam. This is congruous with the operation of the vast majority of space-surveillance radar systems which operate with a single narrow beam, or a fence, configuration [13].
Passive radar systems are being used for space surveillance purposes with wide-area receiver systems, utilising a (noncooperative) broad-beam transmitter at a considerable standoff distance. Such systems include the Low-Frequency Array (LOFAR) radiotelescope, among others [14]- [17]. When combined with a receiver with a wide field of regard, a wide area of broadcast transmission ensures a significant volume above the receiver is illuminated. However, this broad coverage also means the relative power directed toward RSOs can be reduced when compared to the cooperative bistatic configuration. Despite this lower directed power, there are benefits to utilising these types of transmitters. A vast coverage volume ensures RSOs may be detected for a significant amount of time, enough to provide accurate orbit estimates from a single pass without any prior knowledge. Additionally, terrestrial broadcasters, such as TV or radio, consist of a distributed network of many transmitters, allowing multiple transmitters to be used as simultaneous sources of radar illumination of the volume above the receiver. This paper details the use of the MWA as a passive radar for the surveillance of space, illustrated in Figure 1, and is structured as follows: Section II describes the MWA configuration and its unique aspects as a passive radar receiver. Section III covers the signal processing pipeline and includes the techniques to integrate range-compressed pulses to efficiently match objects in orbit, as well as detailing the complete orbit determination (OD) step in Section IV. This OD work builds on earlier work investigating uncued detection and initial orbit determination (IOD) [18]. The steps taken to overcome the challenges of using non-cooperative FM radio transmitters (an analogue broadcasting method) amidst the network's unique configuration in Australia are covered in Section V. The paper includes some illustrative results in Section VI from a short collection campaign utilising four transmitters around the country at various baseline distances of over 2500 km. The paper concludes with Section VII.

II. THE MURCHISON WIDEFIELD ARRAY USED AS A RADAR RECEIVER
The MWA is the low-frequency precursor to the SKA, operating in the frequency range of 70-300 MHz. The main scientific goals of the MWA are to detect radio emissions from neutral hydrogen during the Epoch of Reionisation, to study the Earth's Sun and its heliosphere, Earth's ionosphere, and radio transient phenomena, as well as to map the galactic and extragalactic radio sky [19].
The telescope is strategically located in a legislated radioquiet zone, the Murchison Radio-astronomy Observatory, far from cities and other sources of electromagnetic interference.
The MWA consists of 4096 dual-polarised wideband dipoles configured into subarrays of 4 × 4 square grids, referred to as tiles. The grid spacing for each tile is 1.1 m (corresponding to a half-wavelength separation for 136 MHz), and the 16 dipoles are attached to a 5 m × 5 m steel mesh ground plane. These 256 tiles are distributed over ∼15 square kilometres. The current phase of the MWA, Phase II, allows for the use of half of the tiles at any one time, in compact or extended configurations [20]. Figure 2 shows the layout of the compact configuration of the Phase II MWA.
Each tile's 16 antennas are combined with an analogue beamformer to form a tile beam in a particular direction in the sky. The beamwidth is approximately 40°at the zenith for FM frequencies.
The full frequency of each tile is directly sampled at 327.78 MHz covering the array's frequency range of interest, and a polyphase filter bank (PFB) channelises this data to 256 × 1.28 MHz-wide critically sampled coarse channels. The MWA is able to transfer a user-defined subset of 24 of these 256 channels in real time through to the high timeresolution voltage capture system (HTR-VCS), allowing for an instantaneous bandwidth of 30.72 MHz.
The coarse channels, consisting of 1.28 MHz HTR-VCS data, are then further downsampled to select the various FM frequencies of interest for radar processing.
For typical space surveillance applications, the standard utilisation is a (wide) beam-stare mode, with the tile's analogue beam directed near the zenith. The net output is 128 channels with a bandwith of (typically) 100 kHz. These data are then calibrated to remove any impact of the distributed array and phase-align each channel to the antenna. Once calibrated, the data are able to be used for precise electronic beamforming with the distributed array forming a large aperture providing very narrow surveillance beams. The transmitter reference signals are collected directly at sites close to the transmitter with a small software-defined radio (SDR) system. These SDR reference systems are GPS-disciplined to allow synchronisation with the MWA's surveillance data. Additionally, despite its location in a radio-quiet zone, the MWA is able to detect transmissions from the nearest radio transmitters 300 km away, allowing for the direct observation of nearby reference signals.

III. SIGNAL PROCESSING
An issue faced when using a sensor such as the MWA as a space-surveillance radar is the tradeoff between needing to integrate for a longer amount of time for increased sensitivity, and the changing geometry that orbital motion imparts. The precise and narrow beamwidth, resulting from such a large aperture, means that RSOs in LEO will occupy a beam for only a brief instant. For even a short coherent processing interval (CPI), the object in orbit will need to be tracked spatially throughout the CPI. Similarly, the object will need to be tracked in delay and Doppler parameter space as well. This is the fundamental challenge with the detection of RSOs with the MWA; the need to extend the CPI to detect smaller and more distant RSOs is balanced against the difficulty of coherently forming detection signals. This section covers efficient and scalable methods for dealing with this problem by matching the radar-receiver parameters to the orbital motion. This type of track-processing is achievable given the a priori information on orbits, whereas for uncued detection and searches, earlier studies have detailed practical methods for forming hypothesised orbits [18].
Given a large number of potential orbits, either from known tracks, searching a volume of orbits around known tracks, or uncued search hypotheses, the surveillance data detailed in Section II, and the collected reference signals, are able to be coherently matched to detect an RSO in that orbit.
Rather than directly processing across the received signals, range-compressed pulses are formed for each antenna. These pulses are then matched to the RSO motion [21], [22]. The signals are split into M pulses, each with a duration of τ , such that the CPI length is given by M τ . The pulse length, τ , needs to be sufficiently short to ensure that any change in the Doppler frequency across the pulse is insignificant [22]. Decreasing the pulse length is equivalent to increasing the maximum unambiguous velocity coverage. For a CPI length of 3 s, M will need to be in the order of 40,000 pulses in order to unambiguously span potential orbital velocities.
Given a reference signal s r and N tiles that each have a received signal s n (n ranging from 0 to N − 1), the range-compressed pulses are formed by correlation. Given the sample rate B, each pulse consists of Bτ samples, and the pulse compression forming the range-compressed pulse stack is obtained with the following formula: where t is the fast-time (or delay) sample index, m is the slowtime pulse index, and t is the fast-time correlation index for the two pulses. Note that the sample rate, B, is treated here as equal to the signal bandwidth, although in practice the true signal bandwidth will vary with the analogue content. This pulse stack is formed for each tile, n, creating a compressed pulse cuboid. Typically these pulses are coherently integrated simply with a Fourier transform (FT), resolving  Doppler. However, for rapidly changing geometries, this will not be sufficient, and instead the phase resulting from the target's motion needs to be matched from pulse to pulse.
Constant radial motion will result in a linear phase rate across subsequent pulses, which will be coherently matched by the FT. However, with orbital motion, the target's radial slant range changes rapidly, as will its Doppler signal. This results in a complicated signal with the Doppler changing rapidly, which is typically treated as a polynomial phase signal. For orbital motion, matching higher-order terms in this polynomial phase signal results in significant increases in the signal-tonoise ratio (SNR) [23], [24]. The moving object's phase will be different for each antenna as well, resulting in differing polynomial phase signals from each tile which need to be matched in order to be combined.
Although it is possible to form a single beam (and even a moving beam), and then search in Doppler-rate terms (or indeed vice versa), the parameter space is far too large and this approach is intractable. Instead, each tile's matching polynomial phase signal is determined from a hypothesised orbit in order to best match the orbital motion. This essentially matches the Doppler phase signal and the spatial phase signal in one process.
The bistatic radar configuration is illustrated in Figure 3, with a target at position r, relative to the centre of the Earth, with the velocityṙ. The receiver's Cartesian location is given by r rx and the transmitter's by r tx . The position vectors from the receiver and the transmitter to the RSO are ρ rx , and ρ tx , respectively. The polynomial phase coefficients are derived from this geometry with an orbital motion model. The bistatic radar configuration is illustrated in Fig. 3, with a target at position r, relative to the centre of the Earth, with velocityṙ. The receiver's cartesian location is given by r rx and the transmitter's by r tx . The position vectors from the receiver and the transmitter to the RSO are ρ rx , and ρ tx , respectively. The polynomial phase coefficients are derived from this geometry with an orbital motion model.
For the relatively short duration of a single CPI, the motion model assumes that the only force acting on the object in orbit is Earth's gravity. The acceleration due to the Earth's gravitÿ r, is given by the inverse square law: where µ is the standard gravitational parameter for Earth. This type of orbit, a Keplerian orbit, describes two-body motion and is the simplest orbital model. An orbit is fully determined by the state vector x = rṙ , and so every future or past motion parameter, and then the phase parameter can be determined from the state vector x at a given time.
and this is similarly the case for the transmitter's terms ρ tx , ρ tx ,ρ tx and ... ρ tx . Additionally, ... r is determined from (2) and is given by: Note thatṙ rx ,r rx and ... r rx are the (known) motion terms for the receiver on the Earth's surface, and this is similarly the case forṙ tx ,r tx and ... r tx for the transmitter. In this reference frame with a rotating Earth, the MWA is travelling (instantaneously) 47 m/s faster than the southern-most transmitter used in Section VI.
The expressions for the instantaneous bistatic delay and Doppler are The spatial parameters, azimuth and elevation, and their rates, are determined directly from the orbit as well, ensuring that any RSO is spatially tracked throughout the CPI. These parameters are determined from the receiver's slant range vector rotated from an Earth-centred inertial (ECI) geocentric equatorial reference frame to a south-east zenith (SEZ) topocentric-horizon frame. The rotated vector q and its subsequent ratesq andq are given by the following formula: where D is the SEZ to ECI rotation matrix [25]. Note in some publications q and ρ are instead written as ρ sez and ρ eci respectively.
From these topocentric pointings, the azimuth and elevation, θ and φ, respectively, and their rates, are determined. Given q = [q S , q E , q Z ] T , the expressions for the spatial parameters are: Note that the slant ranges (and their rates) will be unchanged by the rotation, q = |q| = ρ rx , and again, this is also the case for the subsequent rates.
A previous study assumed that two terms for each of the spatial parameters were sufficient [24]. However, some particularly fast moving objects, such as rocket bodies in geosynchronous transfer orbits, require additional parameters. The full angular accelerations were required to detect an SL-12 rocket body (NORAD 20082) travelling at 9.6 km/s (in the ECI reference frame).
An equivalent approach is to rotate the tile locations to the ECI frame, in which case the spatial parameters and subsequent beamforming will have a topocentric right ascension and topocentric declination (rather than the azimuth and elevation) [18].
Given these spatial parameters, the polynomial phase coefficients can be determined for both Doppler and spatial aspects. The Doppler phase coefficients for the first four terms are given by their respective-order Taylor series terms: Similarly, the expression for the spatial coefficients is given by: where u n is the location of the nth tile, and k,k andk are the wavevector and its rates determined from (13) . (26) This full set of matching phase signals ensures that a potential orbit determined by the state vector [r,ṙ] will be completely tracked both spatially and in Doppler across every pulse and every tile. This phase-matching matrix can then be applied to the data, by applying the polynomial phase signal correction to each tile by forming the Hadamard product between the two. These signals are then coherently combined by summing across each tile to form a single, fully matched, slow time-series for a single range bin, χ[m]. The range bin delay sample is determined by the time delay t D , from (8), as well as the sample rate B.
Finally, the Fourier transform will integrate the fully orbitalmatched slow-time pulses resolving the slow frequency, or Doppler, signal. This full process is illustrated in Figure 4. This final signal is passed through a constant false alarm rate (CFAR) detector to produce detections of the matched orbit.
Note that the d 0 term is not included, since the final phase is not of immediate interest; however, the inclusion of d 1 ensures any matched target returns will be close to zero Doppler. This allows for a pruned FT implementation. That is, only those frequency bins sufficiently close to zero Doppler need to be determined to cover any potential orbital velocity offset, and also large enough to encompass sufficient bins allowing for accurate threshold estimation for the CFAR detector.
The process in (27) and Figure 4 only samples a single range bin; creating an entire delay-Doppler map would serve little purpose, since the orbit-derived parameters used to generate that map would only be relevant to a single range. A common pitfall with passive radar and analogue signals is that the ambiguity function is content-dependant; depending on the specific audio signal, the range resolution can be quite poor [26]. By processing a single range bin, this issue can be avoided, or at least moderated. Although signal content might result in a poor ambiguity function, by matching orbits directly, the orbitderived parameters will vary across range bins and only one orbit will have the best-matched parameters.
By forming detections with orbit-derived parameters, every detection will be associated with an orbital track with some confidence, since the beamforming has followed the orbit through the CPI. This associated trajectory greatly assists ongoing tracking and detection-track association. Additionally, having a known trajectory estimate is required for the OD step outlined in Section IV.
This process is entirely flexible and the motion model can be extended to more complicated orbital models, such as incorporating an oblate Earth or other perturbing forces. The measurement model can also be tailored, rather than applying far-field beamforming which is suitable for the MWA's compact configuration. The matched signals are able to be readily extended to near-field beamforming. Instead of calculating beamforming coefficients as well as Doppler coefficients, Doppler coefficients ( (19)-(22)) can be determined for each tile's location and the resulting matched signal can be determined as before.

IV. ORBIT DETERMINATION
Given an orbital track, either from an RSO catalogue or an initial orbit hypothesis, the six dimensional positions and velocities are determined and the processing steps described in Section III are applied. The position and velocity state vectors, or indeed the orbital elements, can be adjusted to form search volumes for RSOs, either to detect manoeuvred targets or to update an old track. If an RSO is detected, a series of associated measurements will be produced, although the process utilises the six-dimensional state vectors. The measurements are produced in the standard radar measurements of azimuth, elevation, bistatic-range, and Doppler.
The number of measurement parameters is extendable in many measurement dimensions. As covered in previous sections, there is a need to account for higher-order motion parameters such as Doppler rates as well as spatial rates. These dimensions are not searched in the processing steps, as that would result in an intractable search space. Instead, only azimuth, elevation, bistatic range, and Doppler are the adjusted measurement parameters (either directly or via the orbital elements); the first three are searched over as part of the Cartesian location and the latter, Doppler, is searched over via the FT as part of the final Doppler-resolving step. If sufficient computational resources existed, it may be possible to independently search through higher-order parameters, which would allow them to be included as part of the OD step.
Orbit determination is achieved using the batch leastsquares method outlined in [27]. This method fits an orbit to a track (or collection) of measurements z. The measurements vector z consists of k measurements such that z = z 0 z 1 . . . z k−1 T , each observed at times t 0 , t 1 ... t k−1 , with each measurement consisting of the detected delay, Doppler, azimuth and elevation such that If the function f maps a state vector x to its respective measurement parameters at times t 0 , t 1 . . . t k−1 (for a single pass, two-body orbit propagation is used such that f consists of Equations (8), (9), (13) and (16); for longer-term orbit determination, more complicated models need to be used), the best orbital fit is the state vector which, when propagated, minimises the residuals between the measurements and the predicted measurements: As f is highly non-linear, finding a general minima is not trivial and instead a solution is found by linearising all quantities around an initial state vector x 0 . This initial solution may be provided a priori from a source such as a previous pass or a space catalogue, or instead from the detections directly using an IOD method. The residuals, , can then be approximated: with ∆x = x − x 0 being the difference between x and the reference state vector, and ∆z = z − f (x 0 ) being the difference between the actual measurements and the predicted measurements for the reference orbit. Additionally, the Jacobian F = ∂f (x) ∂x | x=x0 consists of the partial derivatives of the modelled observations with respect to the state vector. Now, the orbit determination step is achieved by solving a linear least-square problem, with Equation (29) simplified as: x 0 Reference Orbit ∆x ls With each detection and track associated with a state vector, x 0 , a solution to Equation (33) is readily determined [28]. This process is illustrated in Figure 5.
In order to compare different measurement types equally, the residuals are normalised by scaling the measurements (and thus, the Jacobians) by the mean measurement error σ i . In Equation (33), F and ∆z are replaced by F = ΣF and ∆z = Σ∆z, with Σ being the diagonal matrix Σ = diag(σ −1 0 , σ −1 1 , . . . , σ −1 k−1 ).The mean measurement errors in Σ are determined experimentally and also from the measurement resolutions, with the range measurement accuracy determined by the signal bandwidth, the Doppler resolution determined by the CPI length, and the azimuth and elevation resolutions determined by the size of the array aperture. In terms of signal processing, the only aspect that is within the system's control is the Doppler resolution through adjusting the CPI length. It is also possible to further scale the errors using each detection's SNR, such that stronger detections contribute to the orbital fit to a greater extent than the weaker detections [29].
With this process, the linearised error covariance matrix is now obtained with the following formula: The diagonal elements of this covariance matrix yield the standard deviation of the estimate of each element of the state vector.

A. Multistatic Orbit Determination
The OD approach is readily extendable to incorporate multistatic returns, with the detections from additional transmitters simply providing extra measurement parameters to fit the orbit. Having each detection associated with a state vector allows for the easy association of multistatic measurements.
A given detection's position is defined by the narrow beamwidth of the electronically steered beam and its intersection with the isorange ellipsoid defined by the time delay from Equation (8). However, the range resolution achievable using FM radio signals is quite poor, and although the large aperture allows spatially accurate beams, the volume of the intersection will extend radially. This segment, when intersected with subsequent ellipsoids, will not dramatically improve the estimation, as each subsequent coarse ellipsoid will still be intersected with the identical narrow beam.
Conversely, the target's velocity estimate can be dramatically improved. By expressing (9) in terms of the orbital velocity,ṙ, it is clear that every Doppler measurement f D defines a plane of potential velocities: Additional Doppler returns from different sensors will drastically constrain the extent of possible velocities. Two measurements will define a line, and three or more detections will completely determine (or even overdetermine) the velocity. An accurate velocity estimate is the most important aspect of the orbital state vector estimate, as errors in the velocity estimation will produce increasingly erroneous position estimates when propagated forward. Even if multistatic detections are not coincident, the resulting orbit will be improved for having multiple Doppler measurements to constrain the region of possible velocities.
Other benefits of multistatic observations include the diversity of coverage, both spatially and in terms of signal content, as well as resilience by making the most of the vast amount of energy being radiated outward. If the bistatic configuration is relying on narrow elevation sidelobes, there will be gaps in illumination. Using multiple transmitters will help ensure any gaps are covered by at least one other sensor. Additionally, analogue FM radio is not necessarily ideal for radar due to the content-dependant nature of the ambiguity function, and diverse options (even multiple stations from the one tower) will provide resilience in the event one station's content is not suitable. FM radio signals have been shown to be well suited for distributed bistatic radar systems [30].

B. Initial Orbit Determination
The process outlined in Section III requires an orbit to match the parameters in order to form detections. Although these orbits can come from a catalogue of tracks, there will still need to be additional work for uncued detection. The sixdimensional search space of all potential orbit state vectors is currently an unreasonably large search space. However, earlier work has shown that the application of radar constraints can greatly simplify the process. By looking for RSOs at their point of minimum bistatic range, and constraining ranges in orbital eccentricity, this search space can be reduced further [31]. Akin to creating a spatial fence with receiver beams, this approach only searches a narrow region in orbital parameter space. RSOs that do not pass through the right parameter space on one pass will pass through eventually. Tracked RSOs can then be treated as normal, per the general steps in this section.
There may be thousands of objects from a typical space catalogue observable from the MWA at any one point in time, and matching these tracks (and also searching a region around these tracks) would potentially require matching a hundred thousand orbits. An uncued search of a region, even including a constrained orbital search space, will potentially require matching one million orbits.
Depending on the available computer memory, a region's hypothesised orbits and resulting parameters can be stored. The phase matrix coefficients, or indeed the full phase matrices, are applicable for all subsequent CPIs and thus, storing them saves on ongoing computation. The Australian continent provides a unique context for large-scale FM radio-based passive radar. The majority of the population is concentrated in cities near the coast, particularly in the south-east of the country. The FM radio transmitters are similarly located in the population centres. This is illustrated by a map of the most powerful FM transmitters in Figure 6. The MWA is naturally situated far away from these powerful transmitters. This isolation benefits the MWA's astrophysics goals by reducing incident radio frequency interference (RFI). Such a location would ordinarily be less than ideal for terrestrial passive radar purposes, with most passive radar systems designed to be located within the footprint of the illuminating transmitter. However, for detecting satellites at LEO altitudes, such a separation becomes a significant advantage.

V. CONTINENTAL RADAR
Typical FM radio transmitting antennas consist of six to eight antennas combined to form a beam. These antennas are typically angled (as well as electronically beam-steered) towards the ground in order to direct the maximum amount of Fig. 7. Typical FM array's in situ measured antenna pattern via the SixArms airborne radio measurement system [32]. energy to the population [33]. This can pose a challenge for satellite illumination, as every effort is made to minimise the amount of wasted energy being radiated outward from Earth, with the main elevation sidelobes being as low as −15 dB compared to the main lobe [33]. However, these sidelobes will still provide sufficient illumination, and perhaps some Australian transmitters are not so directive, especially given the widespread and low population density found in some rural areas. Figure 7 details a typical FM transmitter pattern. The main beam is directed to the ground; however, sidelobe levels are not insignificant at elevations of up to 15°. At higher elevations, considerable energy is still being radiated in some sidelobes. The patterns of other transmitters may not be so directional. Fig. 8. Incident power on an object above a receiver, for an equivalent isotropically radiated power (EIRP) of 100 kW, with a beam pattern as shown in Figure 7. Note that this figure is based on a spherical Earth model. Instead of relying on sufficient sidelobes for target illumination, it is possible to make use of the main lobe for illumination. Given the number of transmitters in the southeast of Australia, and the lack of large interfering transmitters in between that region and the MWA, a target above the MWA will be illuminated by the main beam of potentially dozens of significant sources. This is illustrated in Figure 8, showing the incident power on an RSO at various altitudes above a receiver with a wide range of potential transmitter-receiver separation distances. It shows that the loss incurred by the increased transmitter to target distance, ρ tx , is more than offset by the transmitter gain of the main lobe, from Figure 7. Indeed, the MWA has been used to detect FM radio returns from the moon, undoubtedly from transmitters half-way around the world [34]. The use of the main lobe is illustrated in Figure 9, showing the SNR of the International Space Station (ISS), utilising a transmitter in Mount Gambier, over 2500 km away from the MWA. This shows that the strongest detections occur at the lowest transmitter elevation angles (despite the larger signal path loss); it even highlights the diffraction of signals along the Earth's surface with detections occurring at transmitter elevations of almost 2°below the horizon.
Finally, the majority of FM transmitters transmit a vertically polarised signal, so the direction relative to the MWA will determine the best receiver polarisation. The vertical polarisation to the south will be coplanar with the MWA's north-south polarisation, and there is a similar situation for transmitters to the east and the MWA's east-west polarisation. However, impacts from effects such as Faraday rotation due to the ionosphere and other factors mean that the presence of matching transmitter and receiver polarisations is not necessarily a decisive factor in the detected SNR.

VI. RESULTS
The results in this section are obtained from 20 min of data consisting of a series of five minute dwells collected in December 2019. As detailed in Section II, these dwells were recorded and channelised in real-time and transferred to storage in Perth. Subchannels were selected such that the full national FM band of approximately 20 MHz was collected. The MWA's analogue beamformer was directed to point at the zenith. In addition to MWA observations, several transmitter reference signals were collected from around the country. These are outlined in Table I and shown in Figure 10. Despite being located in a radio-quiet observatory, a sufficient signal is able to be observed from the nearby transmitters in Geraldton, some 300 km to the south-west of the MWA. Reference observations were collected in Perth utilising receiver hardware identical to that used in the MWA. Located at the Curtin Institute of Radio Astronomy (CIRA), this MWA-like reference receiver is synchronised with the MWA via GPS.
Transmitter reference signals were also collected at locations near Albany and Mount Gambier with a SDR setup. These remote SDR nodes are all GPS-disciplined in order to maintain synchronisation, the collections were manually triggered, and each node was able to record a reference with a bandwidth of 10 MHz. Although 10 MHz is insufficient to collect the full FM band, it is generally sufficient to collect every high-powered FM station from a single site. The transmitter near Mount Gambier, over 2500 km away from the MWA and situated in the south-east of South Australia, is one of the highest power radio transmitters in the country. It should be noted that for each transmitter there are many different FM radio stations, all potentially being broadcast at different power levels. The figures in Table I are simply the maximum licensed power level from that tower, and the true levels may, in fact, be lower.
Data from these remote SDR devices were then transferred to servers at CIRA in Perth, alongside the MWA-collected data, allowing offline space surveillance processing. This is achieved by downsampling the MWA data, as well as the SDR transmitter reference data, to narrowband signals (typically 100 kHz) matching known FM stations, and then undertaking radar processing as detailed in Section III, utilising a 3 s CPI.  Table I. Figures 11 and 12 illustrate some of the aggregate detections from these observations. These detections are formed by parsing the Doppler signal data (from Equation (28)) through a cell-averaging CFAR detector. In passive radar processing (and indeed, noise radar), target signals exist against a noise/clutter pedestal floor formed by the cross correlation of the reference signal against other unwanted/mismatched signals. The CFAR detector estimates this floor from a local threshold region around the cell that is being evaluated, and the SNR is determined by the peak signal against this floor. For these results, we have used a very conservative threshold of 16 dB, greatly minimising the presence of any false detections. System performance would be improved with a more realistic threshold; however, care would need to be taken to ensure false detections are not incorporated into orbital estimates. The MWA is able to maintain tracks of many targets at various ranges. During these short and targeted dwells, the MWA was able to detect every large RSO that was within the MWA's main beam at a range of less than 1000 km. The USSPACECOM catalogue defines a large object as having a median radar cross section (RCS) of 1 m 2 or greater and a medium object having a median RCS between 0.1 m 2 and 1 m 2 ; however, these values are for microwave frequencies which differ from those used in this paper, and should only be taken as a general indicator of size [35]. Additionally, many detections are found outside these limits, including the detection of medium RSOs, RSOs at longer ranges, and indeed RSOs well outside of the main receiver beam. This is consistent with the earlier theorised performance [10]. Predictions of the large RSOs with a closest approach of less than 1000 km indicate the MWA would detect over 1800 RSO passes per day, when used in a beam-stare mode pointing at the zenith. However, detections outside these conservative limits, as well as the ability to rapidly adjust the analogue beamforming, suggest the true number will be larger. Figure 13 shows the detections of an outbound pass of COSMOS 1707 (NORAD 16326), a large (now defunct) satellite. The detections were formed utilising the transmitter near Albany and show the bistatic range, bistatic Doppler, azimuth and elevation. Tracked for almost 90 s, the RSO passes the closest bistatic approach (at zero Doppler) and moves north to the closest approach to the receiver (at its maximum elevation). Figure 14 shows the accuracy of the orbit generated from the COSMOS 1707 measurements. The top row shows the accuracy of the positional and the velocity covariance, from (34). The bottom row shows the accuracy of the position and velocity estimate in comparison to the twoline element (TLE) ephemeris. The two rows are in general agreement as to the resulting accuracy and the results are significantly improved when compared to the initial study [4]. These results are typical of most of the objects the MWA detects with a bistatic configuration.
With many more transmitters at the radar's disposal, there is a scope for increased coverage. Figure 15 shows the SNR of a pass of the International Space Station for every transmitter collected in this campaign. It shows three minutes of detections with almost 50 s of complete overlap for each bistatic pair. The SNR fluctuations shown (for all transmitters) highlight the variable nature of the illuminator coverage due to changes in the transmitter beampattern as well as variations in bistatic RCS. There may be additional contributing factors such as Faraday rotation. The ISS is detected well outside of the MWA's receiver beam, to an elevation of as low as 5°above the horizon.
Although only the ISS and the Hubble Space Telescope were large enough to be detected simultaneously using all transmitters, approximately three quarters of all the detected targets had associated detections from another transmitter. Additionally, every transmitter was able to detect objects that were not detected by any other transmitter, including the comparatively weaker Geraldton site. This highlights that FM broadcast transmissions do not uniformly cover the volume above the MWA, and results will improve with more transmitters being utilised. Despite only being licensed to transmit up to a maximum of 50 kW, the particular configuration of the Albany transmitter, and its elevation sidelobes, produced the largest number of detections of all the transmitters listed in Table I. Figure 16 shows the detection's measurements for the same RSO pass as in Figure 13; however, this time, the detections from the Perth transmitter are included. The spatial parameters are near-identical as expected, with the differing geometry resulting in differing delay and Doppler tracks. When these are combined together in the OD stage, the results are significantly improved. Figure 17 shows the accuracy of the combined orbit, equivalent to Figure 14 showing a single bistatic case. The combined orbit is significantly more accurate than either individual bistatic pairs, particularly the determined velocity. A single detection from Perth (at the 17 s mark) reduces the velocity covariance by an order of magnitude, matching the expectations outlined in Section IV-A.
As mentioned earlier, multistatic detections do not need to be coincident to improve the overall orbit. Figures 18 and  19 show the results from the detections of NADEZHDA 5 (NORAD 25567), a far smaller (albeit still classified as large) RSO at a range of 1000 km. The figures again show detections from both the Albany and Perth bistatic pairs, but instead of being coincident, the set of the detections are separated by over a minute. However, just as before, the multistatic detections greatly improve the accuracy of the orbit, confirmed both by  the reduced covariance as well as compared to the TLE.
An example of the three transmitters is simultaneously shown in Figures 20 and 21. In this example, the satellite OPS 5721 (NORAD 9415) is detected for approximately 20 s with the Albany illuminator; however, these detections are supplemented by a small number of detections achieved utilising the Perth and Mount Gambier illuminators. Despite the short period of detections, the resulting orbit is very accurate when compared against the TLE. Indeed, after only five seconds, the resulting orbit utilising detections across all three transmitters is very accurate.   There are complications when comparing and assessing determined orbits, especially when comparing them to the TLEs. Looking at the covariance of the multistatic results in Figure 17, the increasing number of detections improves the estimate, especially for velocity. However, when compared to the TLE, the error does not improve; rather, it plateaus. This could be due to many factors; however, these results  are within the accuracy of the TLEs themselves, as positional errors generally vary from a minimum error of approximately 1 km at the TLE's epoch up to 5 km, depending on the age of the TLE [36], [37]. These uncertainties could potentially mask any systematic biases or offsets, either from the system itself or from the ionosphere [38], [39]. Longer surveillance campaigns are needed to properly assess any potential systemic issues and to fully evaluate the accuracy of short-arc orbit determination.
The true RCS sizes of objects are challenging to estimate, particularly for a passive radar. Without knowing the precise details of transmitter characteristics, the amount of incident power is not known. Additionally, bistatic RCS is typically a complicated function and without accurate knowledge of the precise size and attitude of RSOs, the bistatic RCS is difficult to determine.
The RCS values for known RSOs can be coarsely estimated with simple shapes, such as cylinders. Comparing the estimates of the detected objects provides additional data points to match against the earlier performance predictions. For example, the medium RCS satellite OV1-5 (NORAD 2122) was detected at a range of 1150 km from the MWA (the corresponding bistatic range was 1600 km) with an SNR  is approximately 0.7 m 2 . For OV1-5, the RSO's length is less than half a wavelength, meaning that for smaller RSOs, the scattering will cause the RCS to decrease rapidly [40]. Conversely, for RSOs that possess trailing antennas that have low RCS at high frequencies, these structures can produce a large RCS at MWA frequencies. Examples such as OV1-5 agree strongly with the initial predictions that the MWA, used as a passive radar, is able to detect objects with an RCS of 0.5 m 2 to a range of 1000 km [10].

VII. CONCLUSION
This paper has described the use of the MWA as a passive radar for the surveillance of space with FM radio illumination. The MWA's high time-resolution and receiving capabilities have been described, and the orbital-specific signal processing methods to form radar products have been detailed, from pulse compression through to forming detections. These orbitalspecific methods are required to track an RSO's motion throughout long CPIs to increase SNR. Following detection, the paper details the orbit determination methods, and how multistatic detections can greatly improve the orbital estimate. To demonstrate and verify these methods, this paper includes the results of a short collection campaign, utilising four transmitters across the country. With the data collected during this campaign the MWA was able to detect and accurately track every large object that passed through its main beam at a range of 1000 km or less. It additionally tracked many other objects outside these limits. These results are in agreement with earlier predictions made of the space surveillance capabilities of the MWA when used as a passive radar receiver.
Multiple transmitters were used to form a multistatic radar network, with multistatic detections allowing for rapid and accurate orbit determination, with additional transmitters also providing greater coverage and resilience. By utilising these many transmitters, the MWA is able to provide persistent and widefield coverage of satellites in low Earth orbit. The MWA is able to achieve this coverage using scalable and efficient signal processing matching the radar processing to the RSO's orbit. The widefield coverage ensures that RSOs are tracked for sufficient time to accurately determine an orbit from a single pass.
Although it is not the SKA's intended purpose, using the techniques described here, the SKA will be a highly capable space surveillance sensor when used as a radar [41]. The SKA will be significantly more sensitive than the MWA, and even a modest increase in sensitivity would enable the detection of signals from a higher altitude (for example, from 1000 km to 2000 km). As shown in Figure 8, higher altitudes would intrinsically allow the utilisation of transmitters at even greater distances, including along Australia's eastern seaboard. Conversely, since transmissions reflected from RSOs are a source of interference for astrophysical observations, knowledge of which RSOs are above the horizon will be important information for the removal of RSO effects from SKA observations. Incorporating a system to remove interference reflected from RSOs would be a natural extension for any radio telescope space surveillance radar.