© The Authors 2025

1 Introduction

Radio astronomy has profoundly transformed our understanding of the universe by enabling the study of celestial objects through the detection of radio waves. These extremely faint cosmic signals face a growing threat from radio frequency interference (RFI), which is becoming increasingly problematic as the radio spectrum becomes more crowded due to advancements in modern technology. Defined as unwanted radio emissions that disrupt astronomical observations (Kocz et al. 2010), RFI arises from various sources such as terrestrial broadcast systems, cellular networks, satellite communications, and everyday electronic devices. Thus, it presents a formidable challenge in radio astronomy, compromising the quality and quantity of astronomical data, and thereby hindering our exploration of the universe. As the sensitivity of modern radio telescopes continue to improve, these instruments face an increasingly invasive RFI environment. These radio signals, which often overlap with the frequency bands used for science, tend to obscure, distort, or mimic the faint astrophysical emissions astronomers aim to detect.

The proliferation of communication technologies has exacerbated this issue. For instance, satellite constellations designed for global internet coverage, such as SpaceX’s Starlink, have been shown to unintentionally emit RFI at low frequencies including in the protected bands of radio astronomy (Di Vruno et al. 2023; Grigg et al. 2023; Bassa et al. 2024). No doubt, the increasing prevalence of these satellite constellations will lead to a worsening of the RFI environment for telescopes across the globe. Similarly, terrestrial sources such as mobile phone networks and radar systems add to the pervasive RFI environment, particularly in densely populated areas. Although international guidelines, such as those established by the International Telecommunication Union, provide some safeguards for radio astronomy, the dynamic and evolving nature of RFI demands more adaptive and sophisticated mitigation strategies.

This paper is structured as follows. Section 2 discusses the current state of post-correlation RFI mitigation approaches. Section 3 briefly describes the trajectory-based RFI subtraction and calibration (TABASCAL) algorithm and its intrinsic workings. In Section 4.1, we describe the basics of Bayesian inference (i.e. Gaussian processes), how they are used to interpolate, and how they can be used as a fringe rate filter. In Section 4.2, we describe the full Bayesian forward model, how we set our priors based on sound theoretical footing, some of the intricacies when modelling RFI sources accurately, and the computational considerations to allow this method to effectively scale up. In Section 4.3, we discuss parameter transformations used to effectively perform optimisation in such a large parameter space as well as an optional method to obtain posterior uncertainty estimates in a scalable manner. In Sections 5.1 and 5.2, we describe the simulation set we generated to test TABASCAL and how we compare our results. In Section 6.1, we analyse the performance of TABASCAL in the raw, recovered visibilities. In Sections 6.2 and 6.3, we discuss the performance of TABASCAL recovered visibilities in the image domain, both in terms of image noise and artefacts and point source recovery. TABASCAL’s phase calibration capabilities are analysed in Section 6.4. The robustness of the prior distribution and potential improvements of the method are discussed in Section 7. Finally, in Section 8, we summarise the method and our key results.

2 RFI mitigation

There are numerous RFI mitigation strategies and in general it needs to be a multi-faceted approach at any given observatory. For reviews on RFI mitigation strategies see Kesteven (2010), Briggs & Kocz (2005) & Fridman & Baan (2001). Among the various mitigation strategies we are interested in offline post-correlation strategies. Post-correlation approaches to RFI mitigation fall broadly into two categories: (1) flagging and (2) subtraction or removal. Flagging is the process of identifying contaminated data samples and flagging them in a way that they are not included in any further downstream data analysis. The second category is removal or subtraction, which our method TABASCAL falls into.

Currently, flagging methods fall into three categories: (1) traditional, (2) machine learning, and (3) modified likelihood. One of the most prevalent methods for traditional flagging is AOFLAGGER (Offringa 2010), based on the SUMTHRESH-OLD (Offringa et al. 2010) algorithm. In the past several years, there has been a lot of work on using machine learning for RFI flagging, such as the work of Mosiane et al. (2016) and, specifically, deep learning triggered by Akeret et al. (2017). Notable deep learning examples are from Kerrigan et al. (2019) and (Vafaei Sadr et al. 2020) as they apply to radio interferometric data. Finally, standing as potential flagging counterparts to our method are modified likelihood models such as Mhiri et al. (2022); Leeney et al. (2023); Anstey & Leeney (2024); Mhiri et al. (2024). We note that the BEAMS formalism, introduced in Kunz et al. (2007) to deal with supernova contamination, can enable the efficient computation of the full Bayesian posterior for RFI flagging. We refer to these methods as modified likelihood as they do not assume the standard Gaussian likelihood appropriate for thermal noise alone; rather, they modify it to include a contribution from outliers such as RFI contamination.

RFI flagging is a staple at radio observatories across the world, however, flagging comes at a cost. For example, with respect to observations by the MeerKAT telescope (Jonas 2018), Sihlangu et al. (2021) presented a detailed analysis of 200 TB of data, (around 1500 observation hours). This analysis showed that more than 23% of all L-band data is currently lost to RFI, with 37% of the band subject to persistent RFI, mainly from satellites. This particularly affects neutral hydrogen intensity mapping efforts in our local universe (Cunnington et al. 2022; Engelbrecht et al. 2024). From existing RFI mitigation approaches at the VLA telescope (Napier et al. 1983) in the 1–5 GHz band, a maximum loss of the order of 30% is expected (Selina et al. 2020). These losses impact the sensitivity required by astronomers to carry out their sciences and, in turn, require additional observing time to compensate for the contaminated data, thus reducing the overall efficiency of radio telescope usage. Data loss can significantly impact the calibration process in radio astronomy, potentially introducing errors in flux density measurements and source localisation (Rau et al. 2009). Epoch of reionisation (EoR) science is an area where precise reconstruction of the cosmic signal is of paramount importance as the signal is buried in the noise. Offringa et al. (2019) shows that even flags from real data applied to a RFI-free simulated dataset can introduce biases that inhibit our ability to detect the EoR.

The alternative to flagging is subtraction. Over the years, many methods for RFI subtraction have been proposed. Some online methods include spatial nulling (Kocz et al. 2010) and cancellation using a reference horn (Mitchell & Robertson 2005). In the offline setting, there is subspace projection (Shiyu et al. 2016), post-correlation filtering (Offringa et al. 2012; Helmboldt et al. 2019), fringe fitting (Athreya 2009), fringe rate CLEANing (Kogan 2010), and deep learning (Zhang et al. 2024). Spatial nulling and filtering depends on the instrument having a multi-beam receiver. Although effective, this limits its use to such instruments, however, it comes at the cost of creating primary beam irregularities. Subspace projection methods also show great promise; however, they no longer show good separation performance when the RFI signal becomes decorrelated within a single time dump. This limits their applicability as many RFI signals decorrelate on the typical integration time scales used. Reducing these times increases data usage dramatically and can quickly become infeasible. Post-correlation filtering, fringe-fitting, and fringe rate CLEANing all depend on the difference in fringe rate between the signal-of-interest and the RFI signal. These are powerful techniques, however, they operate on a single baseline. Since the differential fringe rate can be low on certain baselines, these methods become ineffective in this case. In many instances, this is also where the RFI strength is highest. A method is thus required that avoids many of these pitfalls.

Our proposed method, TABASCAL, is an extension of the method by the same name introduced in Finlay et al. (2023, hereafter Paper I). TABASCAL performs an antenna decomposition of the signal from RFI sources, thus taking advantage of all baseline information simultaneously to separate the astronomical and RFI signals. TABASCAL is in essence a progression, although extensive, of the method proposed by Perley & Cornwell (2003) with additional features such as the potential to calibrate off the RFI signal. Additionally, TABASCAL avoids the need for high time resolution data and therefore is applicable to more standard post-correlation data, including archived observations.

TABASCAL draws on ideas from the work of Roth et al. (2023) with respect to Gaussian processes, but places them in the context of fringe filtering and fitting. It is conceivable to envision a single method combining TABASCAL with the work of Roth et al. (2023) (imaging) and Leeney et al. (2023) (flagging unmodelled effects) to form a unified Bayesian approach to the radio interferometric data reduction process.

3 Overview of TABASCAL

TABASCAL is a method used to separate RFI from astronomical signals in post-correlation radio interferometry data in an offline setting. Additionally, under the right circumstances, it can also jointly perform phase calibration. It employs a Bayesian framework to effectively cast the use of fringe rate filters as Bayesian priors. In this way, multiple filters can be applied simultaneously in multiple directions.

3.1 TABASCAL I and II

This paper is the second in a series. TABASCAL I (hereafter Paper I) addressed the problem of deriving antenna gains from calibrator observations in the presence of RFI moving on regular trajectories relative to the phase centre (e.g. satellites, towns, etc.). Therefore, TABASCAL I operates in the situation where the astronomical visibilities are known (i.e. a calibrator source) and the unknown quantities to be estimated are the antenna gains and the RFI signal parameters. In TABASCAL II, the astronomical visibilities are also unknown. This is the primary distinction between TABASCAL I and II. Currently, TABASCAL estimates a single polarisation product, for example, the XX polarisation, across all antenna pairs (baselines). In the future, TABASCAL will take advantage of the full polarisation nature of the data to enhance its separation capabilities as RFI is often strongly polarised compared to astronomical signals.

Traditionally, the antenna gains, derived from the calibrator observation, would be applied to a target observation in a transfer calibration sense; namely, first-generation calibration (1GC) as defined in Smirnov (2011b). In contrast, in this paper, titled TABASCAL II, we are analysing the target science observations directly. Calibration using the target observation data directly is referred to as self-calibration (selfcal) or second-generation calibration (2GC) in Smirnov (2011a). In 2GC, the unknown quantities are the antenna gains as well as the astronomical visibilities. In our case, we also have unknown RFI in the target observations. Traditional self-calibration is an iterative, cyclic process. Initially, gains are fixed and the visibilities are imaged. The imaging process acts like a prior in the sense that it searches for solutions that are sparse in the image domain. After this, astronomical visibilities are fixed for each round when calculating gain solutions such that the problem is over-constrained¹. In contrast, TABASCAL II solves for the astronomical visibilities, antenna gains and the RFI signals jointly in one pass. As such, TABASCAL II requires prior information. The prior information provided for the astronomical and RFI signals is in the form of the scale of the signal and its time variability, namely, fringe rate filters.

As described above, TABASCAL I and II solve a similar problem but in two different situations, namely, calibrator and target observations. They employ largely the same data model, however, there are some notable differences in TABASCAL II versus TABASCAL I, respectively:

1. Astronomical visibilities are modelled versus known a priori;

2. Estimated RFI signals are complex-valued versus real-valued;

3. Satellite trajectories are calculated using accurate two-line elements (TLEs) instead of the simpler circular orbit model;

4. RFI trajectory errors are absorbed into the complex-valued RFI signal versus a parametrised trajectory model;

5. RFI signal interpolation is performed using a time-correlated covariance function versus linear interpolation;

6. Antenna gains are interpolated to observation times versus a parameter for each time step;

7. We did not fit for the trajectories of the satellites here, assuming them to be previously known within some degree of error. If necessary, discovery can be performed using the methods of the TABASCAL I paper.

All of the differences noted above (except the first, i) are improvements on the previous work and can be retroactively applied there.

3.2 TABASCAL II’s working principle

To learn the RFI signal, gains and astronomical visibilities, TABASCAL II effectively applies a fringe rate filter in multiple directions simultaneously, jointly fitting for the signals from these different directions. These fringe rate filters are applied through the use of the appropriate, derived priors on the associated parameters. By design, the astronomical visibilities for a tracking, fringe stopped interferometer, with a limited field of view, are expected to have a fringe frequency very close to zero (Offringa et al. 2012). The maximum expected fringe rate is proportional to the projected baseline length and the field of view of the telescope. For RFI sources, we manually fringe stop in the direction of each source and then apply a fringe rate filter in this direction. Fringe stopping in a particular direction causes the fringe rate for that source to be close to zero. To effectively separate RFI sources whose fringe rate is close to that of the astronomical sources or each other, each RFI source has its signal modelled at antenna level. When doing this, we need to account for fringe winding loss as RFI signals will decorrelate when the integration time is sufficiently long compared to the source’s fringe rate. We do this by interpolating the RFI signals to sub-integration time scales and then integrating back to the data rate. Using an antenna based decomposition allows us to account for direction dependent gains in the direction of that specific source that may vary across antennas, such as the primary beam, ionospheric effects, and position errors. In essence, TABASCAL extends the concept of peeling (Intema et al. 2009) to RFI sources. The signal separation capabilities of TABASCAL are hinged on fringe rate filtering and peeling of a time-varying source.

TABASCAL II has phase calibration capabilities that depend on the signal-to-noise ratio (S/N) of the RFI signal. This is due to the sensitivity of the predicted visibilities to small changes in the position of the RFI sources. This is thanks to the antenna decomposition of the signal for these sources. The priors used on the astronomical and RFI signal strengths are non-informative. As such, TABASCAL II does not perform amplitude calibration. However, this is possible if a known source is in the field or more a priori information is known about an RFI source. As for phase calibration, TABASCAL II’s capabilities are dependent on the strength of the RFI sources. If the signal is too weak, we do not have a high-enough S/N. In this case, TABASCAL II should return gain estimates that are consistent with the prior which is set by traditional 1GC or the use of TABASCAL I when RFI sources are present.

4 Method

Our proposed method, TABASCAL, makes use of a Bayesian forward model to predict the RFI contaminated visibilities from a set of parameters (including nuisance parameters) to model the RFI signal. The parameters of this model are then estimated by optimising over the posterior distribution to obtain a maximum a posteriori (MAP) estimate. Thereafter, Gaussian constrained realisations can be used to estimate the posterior covariance. All nuisance parameters can be ignored, leaving us with calibrated astronomical visibilities along with optional error estimates. All time-varying signals, such as antenna gains, astronomical visibilities, and the RFI signal at the antennas, are modelled using Gaussian processes where prior knowledge of their time variability can be encoded. We expect the antenna gains to vary on a longer time scale compared to the astronomical visibilities which in turn are expected to vary more slowly than the RFI signal. Encoding this prior information into the probabilistic model allows us to break the inherent degeneracies when estimating more parameters than data points, as is the case here.

In this section, we review the required concepts to build our Bayesian forward model and how to obtain a MAP and covariance estimate. We start, in Section 4.1, by introducing Bayes theorem, Gaussian processes (GP), and two different implementations of GPs. Next, in Section 4.2, we introduce our full Bayesian forward model with all the separate components including their prior distributions. We also explain how to determine the prior distributions and finish by discussing the computational considerations for the two different GP methods. Finally, in Section 4.3, we explain how optimisation is performed using standardised parameters as well as how to estimate the posterior covariance in a scalable manner.

4.1 Concepts

4.1.1 Bayesian overview

The central object in Bayesian statistical methods is the posterior distribution, 𝒫(θ|𝒟), namely, the probability distribution over model parameters, θ, given some observed data, 𝒟. Bayes theorem gives $$𝒫\left(\mathit{\theta }|𝒟\right)=\frac{ℒ\left(𝒟|\mathit{\theta }\right)\Pi \left(\mathit{\theta }\right)}{𝒵\left(𝒟\right)},$$(1)

which is constructed from the likelihood, ℒ(𝒟|θ), the prior, Π(θ), encoding our beliefs about the parameters before seeing the data, and the evidence or marginal likelihood, 𝒵(𝒟), which is a normalising constant that is not required to find the MAP estimate or to draw samples from the posterior. The likelihood term is derived from the noise distribution that we assume in our data model and is what links information from our data to our model parameters, θ. The prior distribution is where we can include any additional information to drive our solutions in a desired direction from information we have about the problem without the knowledge of the data. The fundamental piece of prior information we use in this analysis is the expected time variability of the different signals present in our observations.

4.1.2 Gaussian process

A Gaussian process is a collection of random variables, any finite number of which have a joint multivariate normal distribution. It defines a distribution over functions, in the sense that for any finite set of input points, the corresponding outputs follow a joint Gaussian distribution. In our work, we consider GPs defined over a one-dimensional temporal domain. Analogous to the normal distribution, a GP is fully specified by its mean function and covariance function. The covariance function encodes the pairwise correlations between function values at different time points.

A common class of covariance functions is the Matérn class of covariance functions. The squared exponential (SE) is the limiting case of these and leads to the smoothest functions. In this work, we use the SE covariance function which is defined as $$\kappa \left(t_1,t_2\right)={\sigma }^2\exp \left[\frac{-{\left|t_1-t_2\right|}^2}{2{\tau }^2}\right].$$(2)

Therefore, prior knowledge of the time variability, between t₁ and t₂, and the overall scale of a function and/or signal can be encoded through the definition of an appropriate covariance function. Then, τ determines the length scale of the GP solutions and is referred to as a hyperparameter in the context of our work.

For example, in this work we wish to estimate the gain phase values over the observation period. If we expect these to vary by less than 10° per hour², then we would choose σ = 10° and ${\tau }_{{\varphi }_G}=1$ hour as a prior covariance between the gain phase values on a single antenna. In this work these parameters are determined by a fit to data from adjacent calibration observations which is covered in Paper I. Alternatively, these parameters could also be varied and have an associated hyperprior on them based on the estimate from the calibration observations.

4.1.3 Interpolation and inducing points

GP regression can be viewed as an interpolation when the covari-ance of the known points is assumed to be noise-free. Given a covariance function, κ(t₁, t₂), locations of known points, t′, and new, interpolated function locations, t, an interpolator, ${\mathit{I}}_{tt\prime }$, can be defined as $$\begin{array}{l}\mathit{y}\left(\mathit{t}\right)={\mathit{C}}_{tt\prime }{\mathit{C}}_{t\prime t\prime }^{-1}\mathit{y}\prime \\ ={\mathit{I}}_{tt\prime }\mathit{y}\prime .\end{array}$$(3)

The evaluation of the covariance function, κ(t₁, t₂), over all combinations of locations in t′ gives the covariance matrix, ${\mathit{C}}_{t\prime t\prime }$. Then, y represents the interpolated function values at the points, t, and y′ represents the function values at the locations, t′. Furthermore, ${\mathit{C}}_{t\prime t\prime }$ is the covariance of the prior distribution of the function values at the locations, t′. In this work, the values y′ and locations t′ are referred to as inducing points (Lawrence et al. 2002). The y′ values will be parameters in our model and subsequently estimated. The t′ locations are fixed in our model and regularly spaced at intervals of approximately τ with the end-points of our observation interval included. The y values are the interpolated points which could be at the data rate as is the case for our gains, or at a higher sampling rate as is the case for our RFI signals.

4.1.4 Fourier-based GP method

A particular class of covariance functions commonly used are homogeneous/stationary covariance functions, namely, they depend only on the difference in time. Matérn covariance functions form part of this class. Evaluating a homogeneous covariance function on a periodic domain, with regular spacing, leads to a circulant covariance matrix. Such a matrix is trivially diagonalisable with a Discrete Fourier Transform which can be implemented using the fast Fourier transform (FFT). The Wiener–Khinchin theorem states that the Fourier transform of a stationary covariance function is the power spectrum of the Gaussian process (GP). The relation between the covariance matrix, C_tt, and the power spectrum, Ρ_ηη, is given by $${\mathit{C}}_{tt}=\mathit{F}{\mathit{P}}_{\eta \eta }{\mathit{F}}^H,$$(4)

where F is the Fourier transform, F^H is its Hermitian transpose (inverse), and Ρ_ηη is a diagonal matrix. Then, Ρ_ηη is therefore the covariance in Fourier space η, namely, the fringe rate space when used for visibilities as is the case for us.

4.2 Bayesian forward model

To compute the likelihood in Equation (1), we need a model that predicts the data given model parameters. In this work we simulate data from a radio interferometer, referred to as visibilities. Visibilities are a measure of the signal coherence from the sky in two locations, namely, the locations of a pair of antennas called a baseline. The visibility for an ideal interferometer is defined as $$V\left(\mathit{u},\lambda \right)={\displaystyle {\iint }_{lm}\mathit{I}\left(\mathit{l}\right)\exp \left[-\frac{2\pi i}{\lambda }\left(\mathit{l}-{\mathit{l}}_0\right)·\mathit{u}\right]\frac{dldm}{n},}$$(5)

where V is the visibility, I is the sky surface brightness distribution, l = (l, m, n) are the sky coordinates, u = (u, υ, w) are the baseline coordinates, and λ is the observation wavelength. Here, l₀ is the direction of phase tracking centre and is used to fringe stop in that direction. The visibilities can be thought of as the Fourier transform of the sky surface brightness. This becomes true in the limiting case when the telescope is co-planar (w = 0) or we only consider a small field of view (n ≈ 1). In this context, we are only dealing with the visibilities in our proposed method; however, the link to the sky should not be forgotten and will be used in our analysis of the method.

Our visibility model for the RFI contaminated visibilities takes the simple functional form of $$V_{pq}=G_p\left(V_{pq}^{\text{AST}}+V_{pq}^{\text{RFI}}\right)G_q^{*},$$(6)

where G_p are the complex-valued antenna gains on the antenna, p, while $V_{pq}^{\text{AST}}$ represents the astronomical visibilities and $V_{pq}^{\text{RFI}}$ are the RFI visibilities on the baseline formed by the antennas p and q. All of these terms vary over time but at different rates. The complex-valued noise in the visibilities is assumed to be independent and identically distributed Gaussian noise with standard deviation of σ_n in both the real and imaginary components (Thompson et al. 2017, Chapter 6). The likelihood is therefore $$p\left({\mathit{V}}^{\text{OBS}}|\mathit{\theta }\right)={\left(2\pi {\sigma }_n^2\right)}^{-N_D}\exp \left[\frac{{\left|{\mathit{V}}^{\text{OBS}}-\mathit{V}\left(\mathit{\theta }\right)\right|}^2}{2{\sigma }_n^2}\right],$$(7)

where V^OBS is the vector of all observed visibility data, V is the vector of model visibilities, θ is the vector of all the model parameters, and N_D is the number of data points (i.e. the dimensionality of V^OBS). Therefore, N_D = N_TN_A(N_A – 1)/2. N_T is the number of time steps in the observation and N_A is the number of antennas in the array. It should be noted that the likelihood in Equation (7) is for complex values and, therefore, N_D is the number of complex-valued data points. Finally, σ_n is the noise on the real and imaginary components individually.

Fig. 1 Power spectrum prior used for our astronomical visibilities (orange). The faint black curves are the true calculated power spectrum for all baselines in one of our simulations. The blue curve shows the median power spectrum and the black dashed line shows the maximum expected fringe rate of the longest baseline in our simulation. The parameters for the example priors shown are P0 = 107 Jy2, η0 = 1 mHz. and γ = 5, 2, and 1.

4.2.1 Astronomical visibility model

The astronomical visibilities are modelled with a GP in Fourier (fringe rate) space, ${\tilde{V}}_{pq}\left(\mathit{\eta }\right)$. Then an inverse Fourier transform is applied to obtain functions of time, V_pq(t), for each individual baseline. The power spectrum therefore defines the prior covariance in Fourier space, and the prior mean for each component is 0, namely, $$p\left({\tilde{V}}_{pq}^{\text{AST}}\right)=𝓝\left(\mathbf{0},{\mathit{P}}_{pq}\right),$$(8)

where 𝓝(μ, Σ) represents the multivariate normal distribution with mean μ and covariance Σ. P_pq is diagonal and is the expected power spectrum of V_pq(t). Obtaining the astronomical visibilities over time is done with $$V_{pq}^{\text{AST}}\left(t\right)={\mathit{F}}^H{\tilde{V}}_{pq}^{\text{AST}}\left(\mathit{\eta }\right).$$(9)

The power spectrum is the squared magnitudes of the visibilities in fringe rate space. The analytical form of the power spectrum used in this work is $$P\left(\eta \right)=\frac{P_0}{2}\left[\exp \left(-\frac{{\eta }^2}{2{\eta }_0^2}\right)+{\left(1+\frac{{\eta }^2}{{\eta }_0^2}\right)}^{-\gamma }\right],$$(10)

where P₀ controls the overall power of the signal, namely, the variance of the signal in Fourier space. This power spectrum essentially looks like a leg with a bent knee. Here, η₀ controls the position of the knee in η-space and γ controls the angle of the lower leg. Since the angle of the lower leg determines the proportion of lower and higher frequency modes in the signal, γ therefore controls the smoothness of the signal with larger values increasing smoothness. Figure 1 (orange, purple, and green curves) shows examples of the power spectrum. Then, P₀ should be chosen based on our expectation of the sky signal itself. We used P₀ = 10⁷ throughout this work, which corresponds to a prior standard deviation of $\sqrt{P_0}/N_T=7$ Jy. Here, γ is user-tunable and exists mostly for numerical stability, however, it should remain above 2. Due to our knowledge of the instrument and our pointing direction, we can appropriately choose the η₀ parameter as simply the maximum expected fringe rate for a given baseline. Excluding any strong astronomical sources in our sidelobes, the maximum fringe rate for a baseline would be expected when a strong source is on the edge of our field of view (FoV). In this case, the fringe rate is given by $$v_f={\omega }_E\frac{\left|\mathit{u}\right|}{\lambda }\sin \left({\theta }_{\text{FoV}}/2\right),$$(11)

where ω_Ε is the rotation rate of the Earth in rad s⁻¹, |u| is the baseline length, λ is the observation wavelength and θ_FoV = 1.22D / λ is the field of view of the telescope with a dish diameter of D. For our simulations with a maximum baseline of 800 m, 25 cm wavelength, and a 13.5 m dish diameter, this works out to ν_f = 2.6 mHz. We have used η₀ = 1 mHz throughout this work.

In practice, because the FFT forces periodic solutions on a finite interval, we end up with edge effects in the estimated astronomical visibilities in the time domain. Additionally, many baselines have fringe periods much larger than the simulated observations we are using. To remedy this, we make use of padding.

4.2.2 Gain model

The antenna gains are modelled with a GP using inducing points G_p(t′) which are then interpolated to the points G_p(t). This is done using the interpolator defined by their covariance function as is described in Section 4.1.3. A separate interpolator is used for the gain amplitudes compared to the gain phases as these are expected to vary at different rates and have different variances. However, the same interpolator is used for all antennas as these are expected to vary at the same rate. Therefore, $$\left|G\left(\mathit{t}\right)\right|={\mathit{I}}^{\left|G\right|}\left|G\left(\mathit{t}\prime \right)\right|,$$(12)

and $${\varphi }_G\left(\mathit{t}\right)={\mathit{I}}^{\varphi }{\varphi }_G\left(\mathit{t}\prime \right).$$(13)

The gain amplitudes |G| and phases ϕ_G are then combined as $$G=\left|G\right|\exp \left[i{\varphi }_G\right]$$(14)

to form the complex gains at each antenna over time.

The prior distribution for the antenna gains is chosen based on the estimates obtained from the expected calibration observations that would sandwich the target observation in question. A GP can be fitted to the calibration estimates at the times t″ of the calibration observations in the standard way as is well described in Rasmussen & Williams (2005). We briefly summarise the procedure below.

Given gain estimates y″ with error covariance Σ″ at calibration times t″ we can fit the appropriate SE kernel parametrised by l and σ, giving us $C″=C_{t″t″}$. The gain estimates from the calibration observation can be obtained, even in the presence of RFI, using the method described in Paper I. The prior mean μ′ and prior covariance Σ′ for our inducing points y′ are therefore given by $$\mathit{\mu }\prime ={\mathit{C}}_{t\prime t″}{\left(\mathit{C}″+\mathbf{\Sigma }″\right)}^{-1}\mathit{y}″,$$(15)

and $$\mathbf{\Sigma }\prime =\mathit{C}\prime -{\mathit{C}}_{t\prime t″}{\left(\mathit{C}″+\mathbf{\Sigma }″\right)}^{-1}{\mathit{C}}_{t\prime t″}^T.$$(16)

The prior distribution for our inducing points y′ at the times t′ is then $$p\left(\mathit{y}\prime \right)=𝓝\left(\mathit{\mu }\prime ,\mathbf{\Sigma }\prime \right).$$(17)

The priors for the gain amplitudes and gain phases are assumed independent of each other, as is the case for all prior terms that make up the full Bayesian model. When calculating the prior covariances and interpolators I^|G| and I^ϕ via Equation (3), the σ and l values estimated from fitting a GP to the calibration portions are used.

4.2.3 RFI visibility model

The RFI visibility model is constructed from two parts, these are then combined and then averaged to effectively model fringe winding (time-smearing) of the RFI visibilities. The complex-valued RFI signal at each antenna is modelled as a GP using inducing points and the trajectories of the RFI satellite sources are modelled using two-line element sets (TLEs) with a simplified general perturbation (SGP) model (Hoots & Roehrich 1980). From the trajectory, geometric phase delays are calculated. For a real observation, the actual position of the satellite can differ from the predicted position using TLEs. This position error leads to a differential fringe rate and phase offset between the model and the truth. When this error is not too large, this can be accounted for by the fitted complex signal at each antenna. Alternatively, an implementation of SGP in the auto-differentiation framework PyTorch has already been released by Acciarini et al. (2025) and could be included in the future to account for this.

The overall RFI visibility model looks like $$V_{pq}^{\text{RFI}}\left(t\right)=\frac{1}{\Delta t}{\displaystyle {\int }_{t-\Delta t/2}^{t+\Delta t/2}{\displaystyle \sum _s{A}_{ps}\left(t\right)A_{qs}^{*}\left(t\right)K_{ps}\left(t\right)K_{qs}^{*}\left(t\right)dt,}}$$(18)

where Δt is the integration time of each observed data point, A_ps(t) is the RFI signal at antenna p from the source labelled s, and K_ps(t) is the geometric phase delay term induced by the trajectory of the RFI source. Each RFI source is modelled as a point source, this is valid for nearly all RFI sources except those present within the telescope site Paper (I).

The integral in Equation (18) is calculated using a left Riemann sum with the sampling points defined with regular spacing. The sampling frequency is determined by the component with the fastest time variability which is typically the geometric phase delay term on the longest baseline. The sampling rate, v_s, required to maintain closure relations Perley & Cornwell (2003) is $$v_s>\pi \left|v_f\right|\sqrt{\frac{\left|V_{\text{Inst}}^{\text{RFI}}\right|}{6{\sigma }_n}},$$(19)

where $\left|V_{\text{Inst}}^{\text{RFI}}\right|$ is the maximum instantaneous visibility amplitude of a source and ν_f is its fringe frequency. In addition to this, there is an error, ϵ, introduced by the numerical integration. For a Riemann sum, this is $$ϵ\le \frac{\left|V_{\text{Inst}}^{\text{RFI}}\right|\Delta t}{v_s}.$$(20)

Thus, to maintain an integration error below the visibility noise, σ_n, we have an upper bound defined by the numerical integration scheme used. This results in bounds for v_s as $$\pi \left|v_f\right|\sqrt{\frac{\left|V_{\text{Inst}}^{\text{RFI}}\right|}{6{\sigma }_n}}<v_s\le \frac{\left|V_{\text{Inst}}^{\text{RFI}}\right|\Delta t}{{\sigma }_n}.$$(21)

This bound can be improved upon if using a more accurate integration scheme such as Simpson’s rule; thus, leading to a more computationally efficient implementation.

The fringe frequency is determined by the observation wavelength, λ, the baseline length, B, the velocity of the RFI source relative to the baseline, and the distance of the RFI source from the baseline. The fringe frequency of a moving RFI source is given by $$v_f=-\frac{1}{\lambda }\left(\frac{\left({\mathit{r}}_s-{\mathit{r}}_p\right)·{\stackrel{·}{\mathit{r}}}_s}{\left|{\mathit{r}}_s-{\mathit{r}}_p\right|}-\frac{\left({\mathit{r}}_s-{\mathit{r}}_q\right)·{\stackrel{·}{\mathit{r}}}_s}{\left|{\mathit{r}}_s-{\mathit{r}}_q\right|}+{\omega }_E{u}_{pq}\cos {\delta }_0\right),$$(22)

where r_s and ṙ_s are the position and velocity of the RFI source labelled s, respectively, and r_p and r_q are the positions of the antennas labelled p and q respectively. These should be expressed in the International Terrestrial Reference Frame coordinates. Then, ω_Ε is the rotation rate of the Earth, u_pq is the u-component of the baseline, pq, and δ₀ is the declination of the phase centre. The first two terms are due to the motion of the satellite and the last term is due to the natural fringe frequency induced by the rotation of the Earth (Thompson et al. 2017). For an RFI source that is stationary with respect to the telescope, ṙ_s = 0, and therefore only the last term in Equation (22) remains. This is the more common result shown (Thompson et al. 2017, Chapter 4.3) and other works.

4.2.4 RFI geometric phase delays

The RFI geometric phase delay model is a per antenna term that is dependent on the trajectory of the RFI source. In this work, we consider RFI sources to be in the near-field. Given than the Fraunhofer distance for an interferometer is $$d_F=2{\left|\mathit{u}\right|}^2/\lambda ,$$(23)

where |u| is the baseline length and λ is the observational wave-length. For the MeerKAT telescope with baselines up to around 8 km and maximum wavelength of approximately 1 m, the Fraunhofer distance is 128 000 km. This is well beyond the orbital height of even geostationary satellites. As such, to maintain the correct phases for the longest baselines we need to treat all RFI sources as near-field. The computational cost of using a near-field formulation is larger than that of a far-field but overall becomes a negligible cost in comparison to the calculation of per baseline visibilities from antenna-based signals for the RFI sources. This is due to it scaling with the square of the number of antennas. The geometric phase delay on a single antenna is defined as $$K_{ps}\left(t\right)=\exp \left[-\frac{2\pi i}{\lambda }\left(\left|{\mathit{r}}_s\left(t\right)-{\mathit{r}}_p\left(t\right)\right|+w_p\left(t\right)\right)\right],$$(24)

where r_s(t) is the position of the RFI source labelled s, r_p(t) is the position of the antenna labelled p, and w_p is the w-component of the antenna position. The w-component term comes from the delay compensation introduced in a fringe-stopping interferometer (τ₀ = w_p/c), namely, from phase tracking in the pointing direction. This is a near-field version of the traditional phase delay term for astronomical sources that is present in Equation (5). It assumes spherical wavefronts with a direct path of propagation between the source and receiver.

4.2.5 RFI signal model

The complex-valued RFI signal at each antenna is modelled using a Gp over time using inducing points A_p(t′). The modelled RFI signal represents the intrinsic signal emitted by the RFI source in the direction of the antenna, including any signal corruptions along the way, until it is measured in the antenna by the receiver, as well as any position errors of the source or receiver. Therefore, signal corruptions such as Faraday rotation, ionospheric effects, the antenna voltage pattern, and the antenna gains are all included in the RFI signal A_p(t). For a comprehensive overview of possible signal corruptions see Smirnov (2011b).

We use inducing points as described in Section 4.1.3 to reduce the number of parameters in our model based on the expected variability of the RFI signal. Typically, the deciding factor in the RFI signal variability is the movement of the RFI source through the antenna primary beam sidelobes. Analogously to the gains, the RFI signal inducing points are interpolated to the desired sampling points with $$A_p\left(\mathit{t}\right)={\mathit{I}}^A{A}_p\left(\mathit{t}\prime \right),$$(25)

where the locations t are determined by satisfying Equation (21). The prior distribution of the RFI signal at each antenna is $$p\left(A_p\left(\mathit{t}\prime \right)\right)=𝓝\left(0,{\mathit{C}}_A^{\prime }\right),$$(26)

where $C_A^{\prime }$ is the covariance calculated from the SE covariance function evaluated at the inducing points t′. Therefore $C_A^{\prime }\left({\sigma }_{\text{RFI}}^2,{\tau }_{\text{RFI}}\right)$ is an M × M real-valued matrix where M is the number of inducing points. ${\sigma }_{\text{RFI}}^2$ can be set to the level at which the telescope response becomes non-linear. If the RFI signal is strong enough to push into the non-linear regime of the telescope, this method will not work in its current formulation. Currently, TABASCAL assumes linearity, however, the non-linearity could in future be modelled.

Assuming the fastest variation in the RFI signal is due to the movement through the primary beam sidelobes, then τ value can be chosen based on the apparent sidelobe traversal rate. To approximate this, we divide the approximate angular sidelobe spacing, $\frac{\lambda }{D}$, by the apparent angular velocity, $\frac{\left|{\mathit{\upsilon }}_{\text{RFI}}\right|}{R_{\text{RFI}}}$, of the RFI source, and further divided by four, to give four samples per cycle. Therefore, we have $${\tau }_{\text{RFI}}\approx \frac{\lambda R_{\text{RFI}}}{4D\left|{\mathit{\upsilon }}_{\text{RFI}}\right|},$$(27)

where D is the dish diameter.

Given the combination of a geometric delay term and use of a GP to model the RFI signal at an antenna, this can be thought of as a fringe rate filter, relative to the array centre, in the direction determined by the geometric delay term. Since we use the SE covariance function, this corresponds to an SE power spectrum with η₀ = 1/2πτ. Effectively, this fringe stops in the expected direction of the RFI source and applies a fringe filter determined by the GP. Following this logic we can see that anything that causes a differential fringe rate, from this direction, can be effectively modelled by the GP, for that particular source. A differential fringe rate can be caused by many things in the signal chain including the variability of the RFI signal at emission (duty cycle), estimated position errors, ionospheric fluctuations in the direction of the source, beam attenuation from the source’s antenna, as well as, the telescope primary beam as is accounted for in Equation (27).

4.2.6 Computational considerations

In the previous sections, we describe the GP models used for the astronomical visibilities, the antenna gains, and the RFI signals. Two particular types of GP models were used, namely the Fourier-based method and the inducing points method. These particular GP formulations are not requirements of our method, however, they were chosen with computational considerations in mind.

Since the astronomical visibilities vary on different time scales, ideally, a separate prior covariance is used for each baseline. Unfortunately, this would lead to unfavourable scaling when storing and/or calculating the astronomical visibility covariances. The benefit of restricting ourselves to GPs with covariance functions that are diagonal in Fourier space, is specifically in terms of computational and memory requirements. Given N_T time steps at which we want to infer the signal, a general GP requires us to calculate $𝒪\left(N_T^2\right)$ terms forming the covariance matrix. Inverting this matrix would have $𝒪\left(N_T^3\right)$ computational cost. Since F can be applied via the Fast Fourier Transform (FFT) and Ρ_ηη is diagonal, the inverse of C_tt in this case can be computed in 𝒪(N_T log N_T) with only 𝒪(N_T) memory requirements.

For the RFI signal and antenna gains, the inducing points GP method was used. Since the expected behaviour of these components is roughly the same on each antenna, we are able to use the same prior covariance across the antennas and subsequently the same interpolator as defined in Equation (3). There is little additional value in changing to the Fourier-based GP method as the number of covariance matrices to compute (and potentially store) is independent of the number of antennas/baselines. This could of course change if we are dealing with, for example, a particularly unstable antenna, or a more turbulent ionosphere above a certain set of antennas.

As can be seen from Equation (3), the calculation of an interpolation matrix involves the inversion of a real-valued matrix of size N × N, where N is the number of known locations, t′, one wants to interpolate from. The real-valued interpolator matrix ${\mathit{I}}_{tt\prime }$ has a size of M × N, where M is the number of points to interpolate to, t. It is therefore computationally favourable to reduce N as far as possible without negatively impacting the final solution. A reasonable choice is N = ΔΤ/τ where ΔΤ is the interval over which to interpolate and τ is the correlation length used in the SE covariance function.

For the antenna gains, the number of interpolation locations (inducing points) can be kept very low and is typically kept at N = 3 for a given target observation block of around 15 minutes. However, for the RFI signal, the number of inducing points N and interpolated values M could become very large. If this is the case, it becomes favourable to use a Fourier-based GP model with some adaptations to reduce edge effects and parameter count. We do not investigate this alternative here, however, we do note its potential to improve the scaling of this method in unfavourable scenarios.

4.3 Posterior approximation

In Section 4.2, the different components that make up our Bayesian model are defined. This includes the likelihood and all of the prior terms which are assumed to be independent of one another. The full prior distribution is therefore simply the product of the individual prior distributions. The posterior distribution, calculated with Equation (1), gives us the updated distribution over our model parameters, after the inclusion of information supplied by our data. The central challenge after defining our probabilistic model becomes estimating the posterior distribution.

There are a number of ways to estimate the posterior distribution. The most rigorous is to use a Markov chain Monte Carlo (MCMC) scheme such as Hamiltonian/hybrid Monte Carlo (Brooks et al. 2011), where samples can be drawn directly from the posterior. However, MCMC techniques can be slow to converge and often infeasible in high-dimensional settings such as for our problem. Alternatives include: variational inference (Blei et al. 2017), the Laplace approximation (Tierney & Kadane 1986), and the simplest, approximating with a delta distribution, namely, a maximum a posteriori (MAP) estimate only. In this work, we stick to MAP estimation for computational reasons. The Laplace approximation approximates our posterior with a Gaussian centred on the MAP point. We have implemented a scalable method to estimate the posterior using the Laplace approximation, which can be optionally run after the optimisation step.

Our posterior approximation method is therefore summarised as using a non-linear optimisation routine to estimate the MAP point and thereafter, optionally, estimating the posterior covariance from some notion of the posterior information. In this section, we describe standardised coordinates which help to decorrelate our parameter space and serve to precondition our optimisation routine, and then we describe the method used to estimate the posterior covariance in a scalable way.

4.3.1 Standardised coordinates and optimisation

Standardised coordinates (Knollmüller & Enßlin 2019) are a coordinate system in which the prior distribution follows a simple and uncorrelated model. Most commonly, this is a standard normal distribution. This is also referred to as a non-centred parametrisation (Betancourt & Girolami 2015). In a hierarchical model, the prior distribution itself is also parametrised, in this case, it is referred to as the reparametrisation trick (Kingma & Welling 2013). In this work, we do not employ a hierarchical model, however, our method still benefits from a standardised parametrisation given the correlations in our prior. Given a prior distribution, 𝓝(μ, Σ), over a set of parameters Ψ we can define new standardised parameters $\tilde{\psi }$ as $$\tilde{\psi }={\mathit{L}}^{-1}\left(\psi -\mathit{\mu }\right),$$(28)

where L is a matrix square root of Σ = LL^T. Therefore, the prior distribution over the standardised parameters is $$p\left(\tilde{\psi }\right)=𝓝\left(0,1\right),$$(29)

where 0 is the vector with all zeros and 1 is the identity matrix, i.e. the prior on $\tilde{\psi }$ is the standard normal distribution. $\tilde{\psi }$ can now be used as the base parameters of our model and all parameters described in Section 4.2 can be calculated using $$\psi =\mathit{\mu }+\mathit{L}\tilde{\psi }.$$(30)

Since our model makes heavy use of GP priors, our parameter space is strongly correlated in the space where the prior dominates. This makes both sampling and optimisation in such a space very inefficient due to slow convergence. This is often aided by preconditioning the optimisation routine. By using standardised coordinates, we achieve the same objective. Both methods are equivalent when the pre-conditioner is L⁻¹, however, by explicitly doing this ourselves we have the option to use many out-of-the-box optimisation routines that are available. In our work, we have used the AdaBelief optimiser (Zhuang et al. 2020).

4.3.2 Covariance estimation

The standard way of performing a Laplace approximation is to calculate the posterior Fisher information, which is defined as $$\mathit{I}\left(\mathit{\theta }\right)=\text{𝔼}\left[{\left(\frac{\partial }{\partial \mathit{\theta }}\log p\left(\mathit{\theta }|𝒟\right)\right)}^2\right],$$(31)

and then inverting the resulting matrix to get the covariance matrix. In Equation (31) above, p (θ|𝒟) is the posterior distribution over model parameters θ. In a model like ours, where the number of parameters Ν_P is on the order of 10⁵–10⁶, this is infeasible due to the computational cost of inversion being $𝒪\left(N_P^3\right)$ and memory requirements of $𝒪\left(N_P^2\right)\approx 1$ TB. An alternative method is to draw samples from the Gaussian distribution with the desired covariance. To do this we employ the method of Gaussian constrained realisations (Hoffman & Ribak 1991). We derive this method for the non-linear case in Appendix A. The main result is $$\Delta \mathit{\theta }=\widehat{\mathbf{\Sigma }}\left({\mathit{J}}^T{\mathbf{\Sigma }}_N^{-1}\Delta \mathit{\psi }+{\mathbf{\Sigma }}_{\Pi }^{-1}\Delta \mathit{\varphi }\right),$$(32)

where Δθ are the perturbations about the MAP point obtained from optimisation. $\widehat{\mathbf{\Sigma }}$ in Equation (32) is defined as I⁻¹(θ) and applied in the equation in an implicit manner so that is is never directly evaluated in its entirety for the reasons given above. We also have $$\begin{array}{ccc}\Delta \mathit{\psi }\sim 𝒩\left(0,{\mathbf{\Sigma }}_N\right),& \text{and}& \Delta \mathit{\varphi }\sim 𝒩\left(0,{\mathbf{\Sigma }}_{\Pi }\right),\end{array}$$(33)

where Σ_N is the visibility noise covariance and Σ_Π is the prior covariance. Fortunately, both of these distributions can be trivially and scalably sampled due to their diagonal covariances in standardised coordinates. Obtaining the posterior covariance estimate is then done as, $$\widehat{\mathbf{\Sigma }}=\text{𝔼}\left[\Delta \mathit{\theta }\Delta {\mathit{\theta }}^T\right].$$(34)

Since $\widehat{\mathbf{\Sigma }}$ has such a high dimensionality, this method allows us to efficiently evaluate marginalised covariances without evaluating the entire posterior covariance matrix.

5 Data simulations and evaluation set

In this section we briefly describe the simulation set used to analyse the performance of TABASCAL across a wide range of RFI strengths. We also describe our comparison methods used to compare against TABASCAL’s performance. Most notably, the uncontaminated data which serves as a benchmark for statistically perfect signal separation.

5.1 Simulation set description

To evaluate TABASCAL rigorously, we tested on 92 simulated observations with widely varying RFI signal amplitudes. We have simulated observations with a 32 antenna (from the core) MeerKAT (Jonas 2018) array and a single frequency channel. Figure 2 shows the antenna layout used in our simulations. Each observation is 15 minutes long with two second integration times. A summary of the telescope parameters used is shown in Table B.1 in Appendix B. The gains were set to vary linearly in time with starting values for each antenna drawn from 𝓝(1, 0.05²) and 𝓤(−90°, 90°), for the amplitudes and phases respectively. 𝓤(−90°, 90°) represents the uniform distribution between the bounds −90° and 90°.

A detailed description of the simulations used in this work is provided in Paper I. The only differences in this work are:

1. the use of 32 antennas instead of 64 for computational and storage reasons;

2. the change of astronomical flux distribution from an exponential to a power law;

3. the use of TLE-based satellite trajectories as opposed to a simple circular orbit.

To accurately simulate the contribution from RFI sources we use a geometric delay term, as given in Equation (24), which considers RFI sources to be in the near-field for the interferometric array. Additionally, we include the effect of fringe winding due to the non-sidereal movement of the RFI source through the integration used in Equation (18) where the required sampling is determined by Equation (19). Finally, in Paper I, it was determined all RFI sources outside of the telescope site itself would be in the far-field of a single antenna where the primary beam term for astronomical sources is valid. However, due to RFI sources being in the near-field on the array level, the angle of signal arrival is different for each antenna and, as such, included in the simulations. The intrinsic RFI signal is assumed to be constant in time in the simulations. This will not be true in real data as the RFI source has its own emission pattern and duty cycle. However, in our Bayesian forward model, we allow for this effect by modelling the RFI signal (A_ps(t)), per source (s), and per antenna (p) to be time varying with a user tunable parameter for its variation.

The included satellite-based RFI sources are selected from the GPS satellites that passed the target direction within 45 degrees. This corresponded to between two and nine satellites for each 15 minute simulation. Each observation has 100 point sources uniformly distributed within the field of view (1.42 degrees) with minimum angular separations of 180″. The theoretical synthesised beam size is around 70″ for an uncontaminated dataset. The source fluxes are drawn from a power-law distribution given by $$p\left(S\right)\propto S^{-\beta },$$(35)

where β = 1.6 (Intema et al. 2011) and, minimum and maximum fluxes set to $10{\sigma }_I=14$ mJy and 1 Jy respectively. σ_I is the theoretical image noise for an RFI-free dataset with visibility noise of σ_n = 0.65 Jy per 2 second integration time, as is used in our simulations. The mean astronomical visibility amplitude is 1.1 Jy and the mean RFI visibility amplitude varies from 10⁻³ Jy to 10³ Jy.

Each observation uses an independent sky, gain, and noise realisation. The RFI power was artificially varied to give a large range in RFI visibility amplitudes for testing purposes. For the simulations used in this paper, the maximum fringe frequency, ν_f, of all satellites was approximately 0.6 Hz. From Equation (21), we calculate the maximum S/N for our sampling rate to be 4.4 × 10⁵ and a minimum S/N of 260.

We found TABASCAL is capable of successfully removing RFI up to an S/N of 9.4 × 10³ equating to amean RFI S/N of 2.2 × 10³. However, this is on simulations using a limited sampling rate due to compute constraints.

Fig. 2 Antenna layout of the 32 MeerKAT antennas used in our simulations.

5.2 Evaluation set description

To analyse the performance of TABASCAL over a range of RFI strengths, we have binned the simulations over RFI S/N into five bins. The RFI S/N is calculated as the average RFI visibility amplitude across all data points divided by the standard deviation of the visibility noise. The average RFI visibility amplitude is defined as $$\left|V^{\text{RFI}}\right|=\frac{1}{N_D}{\displaystyle \sum _{i=1}^{N_D}\left|V_i^{\text{RFI}}\right|,}$$(36)

and |V^AST|, the average astronomical visibility amplitude, is similarly defined. Therefore the RFI visibility S/N is defined as $$\text{S}/\text{N}\left(\left|V^{\text{RFI}}\right|\right)=\frac{\left|V^{\text{RFI}}\right|}{{\sigma }_n}.$$(37)

We compare TABASCAL to three other cases: (i) ‘uncontaminated’, (ii) ‘perfect 3σ flagging’, and (iii) AOFLAGGER. The perfect 3σ flagging and AOFLAGGER cases use visibilities that have been correctly calibrated; namely, the true calibration solutions have been applied to the uncalibrated, simulated visibilities. The resulting data is then flagged for RFI using either a perfect 3σ flagging or AOFLAGGER. Flagged data is not used in any downstream processing. This is where TABASCAL is able to gain an advantage as it does not flag the data but rather corrects it so that more data is available for downstream processing. Perfect 3σ flagging, is referred to as such because the true calibration solutions have been applied and the flagging has had access to the true (noise-free) astronomical visibilities. This is obviously not a realistic scenario but is included to show the limits of flagging as a method to address the removal of RFI as compared to using TABASCAL, where the visibilities are corrected and therefore no flagging is applied.

Perfect 3σ flagging is where we have flagged based on the difference between the true astronomical visibilities, V^AST, and the correctly calibrated visibilities, V^CAL. We have flagged where the amplitude of this is greater than three times the true noise amplitude, σ_n. Mathematically, this is represented as $${\text{Flag}}_i=\left|V_i^{\text{CAL}}-V_i^{\text{AST}}\right|>3{\sigma }_n.$$(38)

The AOFLAGGER runs show a slightly more realistic situation. As with the perfect 3σ flagging, we use the correctly calibrated visibilities but then the AOFLAGGER algorithm has been applied. We make use of three passes on the data, each with a dedicated strategy taken directly from CARACAL (Jόzsa et al. 2020), a radio interferometry data-reduction pipeline software³. This shows a more realistic scenario, however, the strategies used have not been optimised for our dataset. Although they are for generic MeerKAT data, which is what we have simulated, except for one channel only.

The uncontaminated case refers to the sum of the true astronomical visibilities and the same noise realisation used in the other RFI contaminated cases. This stands as our reference point and should be the limit that can be reached by TABASCAL when it is working perfectly and no significant amount of information about the astronomical visibilities has been given.

From the 92 simulations, 86% (79) resulted in a ${\overline{\chi }}^2<\mathrm{1.1.}{\overline{\chi }}^2$ is the χ² per data point and is defined as $${\overline{\chi }}^2=\frac{1}{N_D}{\displaystyle \sum _i^{N_D}{\left(\frac{V_i^{\text{OBS}}-{\widehat{V}}_i}{{\sigma }_n}\right)}^2,}$$(39)

where N_D is the number of real-valued data points, $V_i^{\text{OBS}}$ are the data, ${\widehat{V}}_i$ are the estimated visibilities, and σ_n is the noise standard deviation. The visibilities are split up into the real and imaginary components here.

The following results only consider these 79 simulations where TABASCAL successfully converged according to this criteria. In the next section, we show that TABASCAL recovers astronomical visibilities comparable in accuracy to the situation where an idealised telescope has observed the same sky with no RFI contamination, namely, the uncontaminated case, effectively allowing us to ‘see through’ the satellites in the contaminated observation.

Table B.3 shows the prior parameters that have been used throughout the results presented in Section 6. These priors are chosen either to be non-informative, based on expectations of the signal due to telescope and observation considerations in the case of the RFI and astronomical signals, or based on calibration observation information for the gains. The astronomical signal prior is described in Section 4.2.1, the RFI signal prior is described in Section 4.2.5, and the prior on the gains is described in Section 4.2.2.

6 Results

In this section, we discuss the performance of TABASCAL in terms of the astronomical visibility recovery (Section 6.1), the resulting image quality (Section 6.2) and subsequent point source recovery statistics (Section 6.3). We compare the performance of TABASCAL to the uncontaminated ideal case, a perfect 3σ flagging and calibration situation, and finally to correctly calibrated data that was then flagged using AOFLAGGER (Offringa 2010). We finish this section with an analysis of the gain phase calibration capabilities of TABASCAL by leveraging the RFI S/N in Section 6.4.

Fig. 3 Example predictions from TABASCAL for two different baselines (left and right) on an observation containing six satellites. The orange dotted curve shows the TABASCAL prediction and the blue curve shows the true value. The top panel shows this for the RFI visibility magnitude and the lower panels show the real part of the astronomical visibility. The top panel corresponds to a TABASCAL run where γ = 5, as defined in Equation (10). The three lower panels show the results when varying the prior parameter γ for the astronomical visibilities. Finally, γ controls the smoothness of the solutions.

6.1 Visibilities

Due to the presence of multiple RFI sources with differing fringe rates, the resulting RFI visibility signal is non-trivial. In Fig. 3, we show an example of the TABASCAL prediction for the RFI visibility magnitude and real part of the astronomical visibility for two distinct (uυ) baselines, a short (29 m) and a long baseline (715 m). For the astronomical visibility prediction, in the lower three panels, we have changed the prior hyperparameter γ which controls the smoothness of the estimated astronomical visibilities. Equation (10) gives the functional form of the prior covariance in fringe rate space, namely, the power spectrum. Here, P₀ has been fixed to 10⁷ Jy², corresponding to a prior standard deviation of σ^AST ≈ 7 Jy in the time domain (real space). The mean astronomical visibility magnitude is approximately 1 Jy across all observations and baselines. We found no significant effect on the solutions when varying this hyperparameter unless it was chosen to be too small. The main consideration when choosing this hyperparameter should be to encompass the expected signal, namely, σ^AST > 1 Jy for our simulations. We expect a reasonable estimate can be found from a neighbouring uncontaminated channel.

To be noted is the difference in time variations of the true signals between the left and right hand side panels (short and long baselines respectively). From the lower three/six panels, we observe the effect of varying γ on the TABASCAL predicted astronomical visibilities. As γ is increased, the predicted visibilities become smoother and have lower variance. The trade-off here becomes increased correlation between samples. We found that the subsequent imaging results, in Section 6.2, showed equal image noise for all values of γ tested which varied between 1 and 5. As can be seen from the lowest panel in Fig. 3, γ = 1 leads to the astronomical visibilities starting to fit the noise, this can also lead to residual RFI signal leaking into the astronomical visibilities. As γ → ∞, the power spectrum (fringe filter), as defined in Equation (10), becomes Gaussian. This is the desired result, however, this leads to numerical instabilities due to the variance becoming 0 for large η values. As such, the recommendation is to make γ as large while avoiding these numerical instabilities. Throughout the rest of these results we use γ = 5.

For all of these solutions η₀ is kept fixed at 1 mHz, and is chosen according to Equation (11) which is dependent on the telescope configuration and the expectation of no strong off-axis sources. η₀ can have a unique value for each baseline, as calculated using Equation (11). However, due to the fringe period being much longer than our 15 minute observations for many baselines, this would lead to mostly constant solutions on the shorter baselines. Therefore, we have chosen η₀ to correspond to one of the longer baselines. We note that using a unique η₀ as calculated using Equation (11) leads to comparable results and is expected to be the better choice when using TABASCAL on longer observations.

In the top panel of Fig. 3, we can clearly see the accuracy of the RFI prediction on a rather complicated signal from the combination of six satellites. We have chosen a simulation where the RFI S/N is around 2 to show its performance in the lower-mid S/N range. When looking at solutions for an S/N of >10, the errors are no longer visible. This is because the prediction error is constant and does not scale with the RFI strength.

In the top panel of Fig. 4, we show the error distribution in the astronomical visibility prediction from TABASCAL. The individual errors, ϵ_i, are defined as $${ϵ}_i={\widehat{V}}_i^{\text{AST}}-V_i^{\text{AST}},$$(40)

where ${\widehat{V}}^{\text{AST}}$ is the TABASCAL prediction. The real and imaginary parts are concatenated. The black curve represents the distribution for the uncontaminated case where there is only thermal Gaussian noise added. We see that the error distribution from TABASCAL predicted astronomical visibilities has Gaussian errors down to 3σ relative to their fitted distributions. The 3σ error limit is shown as dashed vertical blue bars. We also see that the error distribution is almost identical at all RFI scales except for the largest bin. This is likely because of the limited sampling rate used in our RFI simulations. The accuracy of our RFI simulations diminishes at higher S/N due to this limited sampling rate in our simulations when averaging the sub-integration samples. Section 4.2.3 describes these bounds in detail. We only have theoretical guarantees of simulation accuracy up to RFI S/N of 260. Unfortunately, data volumes above the sampling frequency used would have led to a reduced number of simulations being possible. These deviations from Gaussianity appear to have little effect on the subsequent imaging and point source recovery results.

In the middle panel of Fig. 4, we see that AOFLAGGER does not manage to flag much of the low level RFI contaminated visibilities. This is likely due to the lack of an uncontaminated reference for the algorithm as we only include a single contaminated frequency channel in our simulations. For the perfect 3σ flagging case in the bottom panel, we see the hard 3σ flagging threshold (with respect to the noise distribution) used and note the increased number of very low level visibility errors below the threshold. For both the flagging cases, we see that the error distribution follows very closely to that of the uncontaminated case for RFI strengths below an S/N of 1 as would be expected.

Of note is the reduced error variance for the TABASCAL predicted visibilities vs the uncontaminated case. As can be seen already from Fig. 3, we are able to trade off error variance with increased correlation over time in our solutions through the variation in the γ parameter in the prior. Due to this trade-off, we find that subsequent imaging and point source recovery results remain stable. This is due to the gridding step in imaging which uses a convolution kernel to resample the visibilities onto a grid. We are effectively applying a convolution in the time axis through our astronomical visibility prior.

Fig. 4 Errors in the astronomical visibility predictions from TABAS-CAL and the other cases for comparison as the S/N of the RFI is varied. The black curve shows the distribution of the visibility noise used in the observations and therefore corresponds to the errors in the uncontaminated case. In each panel five bins of RFI strength are used where multiple observations are bundled together. The coloured curves show the Gaussian fit to the error distributions, (i) The top panel shows the errors in the TABASCAL predicted visibilities, (ii) The middle and bottom panels show the errors from the AOFLAGGER and perfect 3σ flagging cases where flagged data is not included in the histograms.

6.2 Imaging

TABASCAL is used to estimate the uncontaminated astronomical visibilities, however, one of the main data products from radio interferometry are images. To evaluate the performance of TABASCAL in this regard, the visibilities from all simulated observations for all four cases have been imaged using WSCLEAN (Offringa et al. 2014). The same imaging parameters⁴ were used for all cases and observations. In Figures 5 and 6, we show an example set of images for all four cases being compared. We have used an observation where the mean RFI visibility amplitude is 1.2 and 15 Jy corresponding to S/N s of 1.8 and 23. These images were chosen based purely on S/N and therefore should be a representative sample of other results at similar S/N values.

TABASCAL shows very comparable image quality to the uncontaminated case at both S/N levels. In contrast, the flagging cases are showing higher image noise and image artefacts. At RFI S/N of 1.8, visible image artefacts of diagonal banding are present for both cases, as well as, more than 50% higher image noise. These become especially apparent at a higher S/N of 23 where the AOFLAGGER image becomes unusable with strong diagonal striping and an image noise eight times higher than TABASCAL. At this S/N level, even the perfect 3σ flagging case has six times higher image noise. This really shows the limits that flagging as a mitigation strategy has, and in all likelihood, this data would not even be included in any downstream data reduction. A multi-frequency flagging strategy would most likely flag all of this data resulting in no possible image in this particular frequency channel. In contrast, TABASCAL shows effectively equivalent image quality to the uncontaminated case and even outperforms it in the specific example shown for the low S/N case.

A key metric in evaluating image quality in radio astronomy is the image noise. In this work, the image noise is calculated by taking the standard deviation of the residual image after running WSCLEAN. We found very similar estimates from PYBDSF (Mohan & Rafferty 2015). Figure 7 shows the image noise for all four cases where observations have been binned according to the mean RFI S/N. The uncontaminated case, of course, does not include any RFI. Therefore, when we show the uncontaminated as a function of the RFI S/N, we are displaying the results from the same observations as the other cases (i.e. we have the same pointing and visibility noise realisation). The dotted lines give the median across the images in a bin and the shaded regions give the 68% uncertainty interval. We see that TABASCAL (orange) is statistically consistent with the uncontaminated case for all RFI S/N bins. In contrast, the flagging cases are only consistent with the uncontaminated case when the RFI S/N is below 1. Above an S/N of 1, we start to see the flagging cases deviate significantly from the uncontaminated case. The perfect 3σ flagging case performs significantly better than AOFLAGGER on this particular metric which is more than likely due to AOFLAGGER not flagging all of the RFI in the data as was shown in the middle panel of Fig. 4. Nonetheless, for RFI S/N values greater than approximately 100, the image noise for perfect 3σ flagging is at least an order of magnitude greater than for TABASCAL.

Fig. 5 Images constructed from the same observation with our four different cases. Top left: No RFI contamination. Top right: TABASCAL fully removes the RFI and recovers the astronomical signal with comparable image noise to the uncontaminated data. Bottom left and right: the images from perfect 3σ flagging and AOFLAGGER respectively, showing significant striping due to residual RFI contamination. The mean RFI S/N in this data was 1.8, showing that significant issues occur for traditional flagging methods even at weak RFI. Here, a perfect 3σ flagging means that any RFI with true amplitude greater than 3× the noise is perfectly removed.

Fig. 6 Images constructed from the same observation with our four different cases. Top left: no RFI contamination. Top right: TABASCAL fully removes the RFI and recovers the astronomical signal with comparable image noise to the uncontaminated data. Bottom-left and right: the images from perfect 3σ flagging and AOFLAGGER respectively showing significantly higher image noise of 6× and 8× respectively. The AOFLAGGER image shows significant striped image artefacts which largely invalidates its use for science, with purity and completeness of around 20% and 50% respectively (see Figures 9 and 8). The mean RFI S/N in this data was 23. For large RFI amplitudes the quality differences between TABASCAL and the flagging methods increases further.

6.3 Point source recovery

Although the TABASCAL algorithm we have presented is fully general, to demonstrate its performance we want a specific test case with clear metrics. For this purpose we choose point-source imaging which also has many science applications such as the radio luminosity function (RLF) Sadler et al. (2013); Simpson (2017); Mauch et al. (2020); Tucci & Toffolatti (2021).

In this section, we evaluate the performance of our four cases with respect to point source recovery. The images generated for all 79 simulations from the previous section were used here. For this analysis, we consider three metrics to be of primary importance: completeness, purity, and flux estimation error.

To extract a source catalogue from our images to compare against the input (true) source catalogue we used PYBDSF (Mohan & Rafferty 2015) with the thresh_isl and thresh_pix set to 1.5. To ensure suitable statistics, we generated point sources with a minimum source separation of approximately 3 beam widths and fluxes with an expected S/N of greater than 10. Both of these values are estimated based on an uncontaminated observation. However, both the beam width and image noise are expected to be larger in images with significant amounts of flagged data leading to potentially overlapping sources. We used a 5σ cut (based on the theoretical image noise) on the PYBDSF catalogues to ensure valid source detections. Afterwards, these were matched with the true source catalogues using match_to_catalog_sky from ASTROPY (Astropy Collaboration 2013, 2018, 2022). Source pairings were only considered a match if the angular separation was less than the beam width. Finally, any true source that was matched to more than one detected source was cut to only include the closest match, in angular separation, to the true source.

Catalogues for all observations have been binned into five RFI S/N bins. In the following Figures 8 to 10, the solid lines with dots represent the median within a bin and the shaded region represents the 68% uncertainty interval within a bin. Each bin consists of approximately 16 observations. The same imaging and source extraction parameters were used for all datasets across all four cases.

Completeness, often referred to as recall, is defined as $$\text{Completeness}=\frac{\text{TP}}{\text{TP}+\text{FN}},$$(41)

where TP stands for true positive (the detected sources that are in the true catalogue) and FN stands for false negative (the undetected sources that are in the true catalogue). Therefore, completeness refers to how complete our detected source catalogue is. Figure 8 shows the completeness statistics for our four cases binned with respect to RFI S/N. TABASCAL shows comparable performance to the uncontaminated case for all RFI S/N bins. Both flagging cases have the same performance as TABASCAL for RFI S/N below 1 and start to deviate above this. Significant performance degradation is seen at RFI S/N above 10 with AOFLAGGER consistently performing worse than perfect 3σ flagging as expected.

Purity, often referred to as precision, is defined as $$\text{Purity}=\frac{\text{TP}}{\text{TP}+\text{FP}},$$(42)

where TP is as before and FP stands for false positive (the detected sources that are not part of the true catalogue). Therefore, purity refers to how pure our detected source catalogue is. Figure 9 shows the purity statistics for our four cases binned with respect to RFI S/N. TABASCAL shows statistically comparable performance to the uncontaminated case across all RFI S/N bins. Again, the flagging cases perform comparably to TABASCAL for RFI S/N below 1. However, on this metric the flagging cases so show very steep performance degradation with increasing RFI S/N. This is likely due to the significant image artefacts that can be seen in Figures 5 and 6 for both flagging cases. The images shown are not hand picked and therefore show a representative sample of the broader dataset. At an RFI S/N level of about 30, more than 80% of the sources detected in the flagging cases are fake sources. These are fake sources that exist at a greater than 5σ level relative to their image noise, not the theoretical image noise.

We note here that all images appeared to show some low level of ghost sources. This was true for all cases. This is likely due to the particular set of imaging parameters coupled with the sparse uυ sampling in these observations. This manifests in slightly lower than perfect purity scores considering the 5σ flux cuts that were used along with only including 10σ_I flux sources. In spite of this, we believe these results still give us valuable insight into the performance of TABASCAL and other methods as all cases were affected in the same way. This is further boosted by the fact that the theoretical image noise was attained for all cases with RFI S/N below 1 and for TABASCAL and the uncontaminated across all RFI S/N levels.

Figure 10 shows the mean absolute flux error of the matched sources from the detected catalogue. The solid lines with dots represent the median value within an S/N bin, and the shaded region shows the 68% uncertainty interval. Each bin consists of results from approximately 16 images for each case considered. The mean absolute flux error (MAFE) is calculated as $$\text{MAFE}=\frac{1}{N_s}{\displaystyle \sum _{i=1}^{N_s}\left|{\widehat{S}}_i-S_i\right|,}$$(43)

where Ŝ_i is the predicted total flux for matched source i, S_i is the true flux, and N_s is the number of matched sources from a single image. TABASCAL shows comparable MAFE performance to the uncontaminated case across all RFI S/N levels. Much like all other metrics shown, the flagging cases show comparable performance to TABASCAL at S/N levels below 1 and then quickly deteriorate for RFI S/N levels above this. As with all the other metrics shown AOFLAGGER performs worse than the perfect 3σ flagging case as expected. On the MAFE metric, AOFLAG-GER shows a faster degradation in performance, with increasing RFI S/N, compared to the other metrics relative to TABASCAL. This is likely due to the image artefacts that arise. Examples of these image artefacts are shown in Figures 5 and 6. When considering the purity performance, shown in Fig. 9, in conjunction with the MAFE results this conclusion seems likely. At RFI S/N levels of 40, perfect 3σ flagging performs around five times worse than TABASCAL with AOFLAGGER performing roughly 30 times worse.

Fig. 7 Image noise calculated from the residuals after using WSCLEAN for imaging vs RFI S/N (lower x-axis) and ratio of RFI-to-Astronomical visibility strength (upper x-axis). The solid lines with dots represent the median within an RFI bin containing roughly 16 images each. The shaded region represents the 68% uncertainty interval over the observations in the bin. TABASCAL (orange) is statistically consistent with the uncontaminated case showing successful signal recovery despite the RFI contamination across the entire range of RFI strengths. In this simulation the astronomical visibilities are chosen to be close in amplitude to the visibility noise so the upper and lower x-axis scales are similar but not identical.

Fig. 8 Completeness statistics of source catalogues across RFI S/N bins. Solid lines with dots indicate the median within a bin. The shaded region indicates the 68% uncertainty interval. Each bin consists of statistics from approximately 16 images for each case.

Fig. 9 Purity statistics of source catalogues across RFI S/N bins. Solid lines with dots indicate the median within a bin. The shaded region indicates the 68% uncertainty interval. Each bin consists of statistics from approximately 16 images for each case.

Fig. 10 Mean absolute flux error with respect to all matched sources across RFI S/N bins. Solid lines with dots indicate the median within a bin. The shaded region indicates the 68% uncertainty interval. Each bin consists of statistics from approximately 16 images for each case. TABASCAL (orange) is statistically consistent with the uncontaminated case.

Fig. 11 Evidence that TABASCAL achieves phase calibration off the strong RFI sources: A comparison of the root mean squared error (RMSE) in gain phase estimates and the posterior standard deviations from TABASCAL. TABASCAL shows statistically rigorous phase calibration capabilities that depend on the RFI S/N. Above an S/N of 1. phase constraints are inversely proportional to the RFI S/N. Below an S/N of 1, phase solutions are found within the prior information given.

6.4 Gain phase calibration

Figure 11 shows TABASCAL’s phase calibration capabilities are shown with increasing RFI S/N. The root mean squared error (RMSE) in the gain phase, for a single observation, is calculated as follows: $$\text{RMSE}=\sqrt{\frac{1}{N_A{N}_T}{\displaystyle \sum _{p,i}{\left({\widehat{\varphi }}_{G_p}\left(t_i\right)-{\varphi }_{G_p}\left(t_i\right)\right)}^2}},$$(44)

where ${\widehat{\varphi }}_{G_p}\left(t_i\right)$ is the TABASCAL estimate of the gain phase on antenna p and time t_i, and ${\varphi }_{G_p}\left(t_i\right)$ is the true gain phase value. The posterior standard deviations for a single observation are also combined in a root mean squared sense. Results from all observations are binned according to the mean RFI S/N values in each. The median and 68% uncertainty interval, across observations, are represented with a solid dotted line and shaded region respectively. The RMSE shows good agreement with the posterior errors indicating that TABASCAL’s MAP estimates and errors are statistically consistent. The proportional decrease between the posterior uncertainties and the RFI S/N indicates that the phase calibration capabilities of TABASCAL are a legitimate feature of the model. The prior distribution on the gains was chosen to have standard deviations of 1% and 1 degree in the amplitude and phase respectively. The mean of the prior distribution is set by taking a sample from a distribution with equivalent standard deviation centred on the true value. Therefore, the prior distribution is not centred on the true value but is statistically consistent with the true value.

7 Discussion

7.1 Considerations for real data

In this work, we tested TABASCAL on simulations that include the dominant factors required to mimic real data. However, there are some scenarios and effects that have not been included that may be present in real data that need to be considered. The main considerations are: (1) multipath propagation effects, (2) non-linear telescope response, and (3) unmodelled sources of RFI. We discuss these below.

In our simulations, we only considered the direct line of sight between the RFI source and the antennas. Although our model for the RFI signal for each source at each antenna allows for some path difference, it is not flexible enough to account for multipath effects such as diffractions and reflections from environmental features. The contribution of these effects is expected to be negligible to non-existent for satellite-based RFI sources. Although TABASCAL is general and can be applied to nearly all off-site sources of RFI, satellite-based sources are the main focus of this work and of the largest growing concern given the advent of satellite mega-constellations (Witze 2025).

Our simulations (and the forward model) assume the telescope’s response to the RFI signal to be linear. This is known to be an invalid assumption when the RFI signal is strong enough. This is especially true when the RFI signal is so strong that the receiver is saturated and any signal recovery is impossible. In our Bayesian forward model, the potential for non-linear effects is not taken into consideration. In principle, this can be modelled and included into the simulations and forward model. Currently, however, this stands as a limitation for TABASCAL, although there is interest in researching this direction. For now, there is still a large range of RFI signal amplitudes, where the receiver remains linear, that TABASCAL is applicable to and will greatly benefit the sensitivity of observations.

A valid concern with a model-based approach like TABAS-CAL is when unmodelled sources of RFI, such as other satellites, planes, broadcast towers, and so on are present in the data. An impact analysis of such a case has not been performed in this work. Given that TABASCAL is a model-based approach, any unmodelled signal from other sources will lead to a degradation of the results. This will be apparent from the increased χ² per data point when the signal contribution is comparable to or greater than the data noise. A residual fast varying, ‘model-free’ signal component could be included in the future, however, this is likely to degrade the astronomical signal recovery. The best solution is to determine the sources of unmodelled signals and include them into TABASCAL’s model.

7.2 Effect of the prior and reduced information

To establish the limits of the model used within TABASCAL a number of tests were run with respect to reducing prior information. In this section, we also discuss the limitations of the priors and what assumptions are present.

7.2.1 Astronomical visibility prior

Variation of the astronomical visibility prior can be done through the three parameters defined in Equation (10), namely: P₀, η₀, and γ. First, P₀ defines the prior variance of the astronomical visibility signal. The only significant effect found on the solutions when varying P₀ was when it was made too small. This effectively excludes higher signal amplitudes leading to an under-estimation of the visibilities. Then, η₀ can be calculated based on the geometry of the telescope, as per Equation (11). The main assumption in this is that no strong sources are present in the sidelobes. In the presence of such sources, the same equation can be used but the angular separation between the furthest off-axis source and pointing direction should be used instead of the telescope FoV radius. Alternatively, such sources could be modelled separately within the TABASCAL framework. The traditional method for dealing with this is source peeling. We did not test this scenario. For mid-frequency observations where there is a limited FoV, adapting Equation (11) is expected to work well. For wide FoV observations the increased η, namely, widening of the fringe filter, could lead to RFI signal leakage into the predicted astronomical signal.

As stated previously, γ controls the smoothness of the solutions. During our investigations, we found that values of γ ≤ 3 led to a bias in the gain amplitudes at RFI S/N below 1. The gain amplitudes would consistently be overestimated in this range, which led to an underestimation of the astronomical visibilities and subsequently the source flux recovery. For RFI S/N values above 1, this bias went away, indicating that RFI signal plays a role in constraining gain amplitudes to within the prior distribution. γ plays the role of the steepness on the edges of the effective fringe filter induced by the astronomical visibility prior. A steeper drop-off (larger γ) led to more consistent and less biased results. We note that we did not perform any tests with γ > 5. However, a previous (less efficient) implementation where an SE kernel in the time domain was used, did not exhibit this bias. The SE kernel is the limiting case of our power spectrum model where γ → ∞. γ is required in TABASCAL’s current implementation for numerical stability. However, this will be improved in the future. Our recommendation is to use as large a γ as possible.

Tests of limiting cases were performed with respect to the presence or absence of astronomical and RFI signals. We found that excluding astronomical signal in the data led to the same results as have been presented in this paper. This was the case both when including or excluding RFI parameters. Excluding RFI signal led to equivalent results to those of the RFI S/N < 10⁻¹. The same results were found both when including or excluding the astronomical visibility parameters.

7.2.2 RFI signal prior and old TLEs

The prior on the RFI signal consists of two parameters, namely, ${\sigma }_{A_{\text{RFI}}}^2$ and τ_RFI. ${\sigma }_{A_{\text{RFI}}}^2$ defines the prior variance of the RFI signal at an antenna. This roughly corresponds to the standard deviation in RFI visibility amplitude for a single source. Throughout this work, we have used ${\sigma }_{A_{\text{RFI}}}^2={10}^4$ Jy, which is larger than any RFI signal used in this work. Reducing this below the RFI strength in the data leads to an underestimation of the signal and, subsequently, a χ² per data point above 1.1⁵. τ_RFI defines the correlation time of the RFI signal. Alternatively, 1/2πτ_RFI can be thought of as the width of a fringe filter (relative to the array centre) about the estimated direction of the RFI source. This allows for the RFI signal parameter A_RFI to fit for a range of residual signal effects as described in Section 4.2.5. Equation (27) shows the calculation of τ_RFI to account for, at minimum, the primary beam modulation. For our simulations this was calculated to be τ_RFI = 24 s corresponding to 1/2πτ_RFI = 6.6 mHz or 13 mHz on a baseline. When testing TABASCAL with one day old TLEs, corresponding to an average position error of 190 m, we found that 1/2πτ_RFI = 6.6 mHz was enough to account for the position errors in 71% of observations compared to 86% in the base case. We did not test the robustness of solutions to ionospheric fluctuations or beam irregularities. This is left for future work. When reducing τ_RFI to 5 seconds (32 mHz) we found statistically comparable results on all metrics, across all RFI S/N levels relative to the base case. We did not test the lowering of τ_RFI in conjunction with old TLE data to see if the convergence rate was increased.

7.2.3 Gain prior

The prior on the gains consists of two parameters each for the amplitudes and phases, the standard deviation $\left({\sigma }_{\left|G\right|},{\sigma }_{{\varphi }_G}\right)$ and correlation time $\left({\tau }_{\left|G\right|},{\tau }_{{\varphi }_G}\right)$. When increasing σ_|G| from 1% to 5%, we found image noise and completeness to be statistically comparable to the base case. We found flux recovery and purity to be slightly diminished across the RFI S/N range, as expected given the apparent lack of amplitude calibration capabilities. When increasing ${\sigma }_{{\varphi }_G}$ from 1° to 5°, we found decreased performance in the astronomical source flux recovery at RFI S/N below 1 and equal performance to the base case above 1 thanks to the phase calibration capabilities leveraging the RFI signal. The recovered astronomical source fluxes were underestimated (below RFI S/N 1) likely due to the signal decoherence smearing the astronomical sources in the image. When decreasing l to 30 seconds compared to 3 hours, for both the amplitudes and phases, we found no effect on the results. This is an indication that TABASCAL would be able to handle more rapid gain fluctuations; however, this was not tested.

7.3 Computational costs and potential improvements

The computational costs of this method are dominated by the calculation of the RFI visibilities from the interpolated RFI signal, and phases, at each antenna. The number of floating-point operations (FLOPS) can be calculated as $$N_{\text{FLOPS}}\approx 2N_A^2{N}_f{N}_{\text{RFI}}\Delta T{v}_s,$$(45)

where ΔΤ is the total observation time considered and v_s is the sampling frequency determined by the strongest RFI source. The memory requirements, at single precision, are approximately 4N_FLOPS Bytes. The optimisation time, namely, the cost of a MAP estimate, assuming 2000 gradient steps, is approximately a microsecond per FLOP. Here we are using FLOPS to define the problem size and not the actual number of FLOPS performed by the processor. An accurate TLE estimate and good initialisation of the RFI signal parameters requires only 1000–2000 gradient steps. An outdated TLE estimate from the day before typically requires 6000 gradient steps.

In our work we have used a single Nvidia A100 GPU. The problem size used in this work consists of N_A = 32, N_f = 1, N_RFI = 2–9, and ν_sΔΤ = 1300–6200. This resulted in runtimes between 14 and 140 seconds for a 15 minute observation, namely, 6×–60× faster than the observation time. Extrapolating this to an SKA Mid size telescope of N_A = 192, an average of N_RFI = 10 at mid-earth orbit, and mean RFI S/N of 100 (v_s = 10 Hz), we expect approximately 7 seconds of compute per second of observation for each contaminated channel.

The above calculation has ignored the increased baseline length present in a larger array. Longer baselines lead to an increased fringe frequency as shown in Equation (22). Increased fringe rates also lead to more signal decoherence, so these baselines would not be as greatly affected by RFI. Baselines on which the RFI fringe rate is so large that the RFI signal is below the noise already could be excluded. Traditional fringe filtering methods such as that proposed by Offringa et al. (2012) will be very effective on these longer baselines. TABASCAL best serves the shorter baselines where fringe rates are lower, making it harder for traditional methods to separate the astronomical and RFI signals.

The computational efficiency of TABASCAL can be improved in a number of ways, paving the way for its widespread use. Currently, the sampling frequency v_s is chosen based on the fringe rate, v_f, of the fastest fringing source, additionally, the expensive portion (complex multiplication between antenna signals) does not take into account the fringe rate on a particular baseline. Ideally, the sampling frequency for each RFI source and baseline would be used to calculate their associated visibility value, after which, per source visibilities can be summed to give the total RFI visibility. This is the single largest computational improvement that can be made and would open up TABASCAL to be used across an entire array with very long baselines efficiently. The next largest computational improvement would be in the interpolation of the RFI signal. This can be done through the use of Fourier space GPs, as is used for the astronomical visibilities. Currently, interpolation is performed in the time domain where a full matrix multiplication is required at a cost of 𝒪 (ν_sΔΤ²/τ_RFI). A Fourier domain interpolation would only cost 𝒪 (N log(N)) with N = v_sΔT. This is only beneficial when ΔΤ/τ_RFI > log(ν_sΔΤ) which is almost always the case.

7.4 Improved prior information

Increased prior knowledge can improve the computational efficiency of TABASCAL, the stability of the algorithm, as well as potentially its capabilities for amplitude calibration. As stated in the previous section, better TLE estimates lead to faster convergence. This is a statement about the accuracy in positional information of the RFI sources. In some sense, this can be thought of as one component of phase calibration in the direction of the RFI source. If the primary beam is well modelled and included as an independent factor in TABASCAL’s model, τ_RFI can be increased as the RFI signal variation due to the source’s movement through the sidelobes does not need to be accounted for in the GP signal model. Increasing our prior knowledge about the intrinsic RFI signal of a given source could even lead to amplitude calibration capabilities using the RFI signal. For example, knowledge of the beam pattern and pointing direction of an RFI signal can help us further constrain τ_RFI. Improvements like these can help TABASCAL to narrow the effective phased-up fringe filter and improve its convergence and success rate beyond what has been shown in this work.

Currently TABASCAL uses an exceptionally wide prior on the astronomical visibilities about 0. In the future, an informative prior could be used based on imaging of a neighbouring uncontaminated channel. This could further improve the convergence and success rate of TABASCAL.

8 Conclusions

8.1 Method summary

TABASCAL solves the problem of estimating phase calibrated astronomical visibilities in the presence of RFI sources that follow predictable trajectories. To achieve this, TABASCAL requires knowledge of the RFI sources present in the data as well as a reasonable prior estimate of the RFI trajectories. Furthermore, the priors used on the RFI and astronomical signals in TABASCAL can be thought of as independent statistical fringe filters, each phased-up into the direction of the corresponding sources. This gives TABASCAL the power to effectively separate the RFI and astronomical signals while accounting for differential fringe rates in all considered directions.

8.2 Key results summary

Overall summary: TABASCAL is able to ‘see through’ predictable RFI (satellites, ground stations, etc.) as if it were not there. This may be a practical solution to the challenge of ever more satellites. The key results are summarised below.

Astronomical visibility recovery

• Accuracy: TABASCAL achieves accurate recovery of astronomical visibilities, even under strong RFI contamination. Predictions are highly Gaussian out to 3σ.

• Comparison to flagging: Both AOFLAGGER and perfect 3σ flagging struggle at RFI S/N strengths above 1, failing to match TABASCAL’s accuracy due to incomplete RFI mitigation.

Imaging performance

• Image quality: TABASCAL produces images with noise and artefacts comparable to the uncontaminated case across all RFI levels. In contrast, perfect 3σ flagging and AOFLAGGER images exhibit significantly higher noise and artefacts, especially at high RFI S/N.

• Noise metrics: Image noise for TABASCAL remains statistically consistent with the uncontaminated case, unlike flagging approaches, which show significant noise increases at RFI S/N greater than 1.

Point source recovery

• Completeness: TABASCAL maintains high completeness across all RFI levels, comparable to the uncontaminated case and significantly better than flagging methods.

• Purity: TABASCAL consistently outperforms flagging methods with performance consistent with uncontaminated.

• Flux estimation: Flux estimation errors are low for TABAS-CAL, with performance comparable to the uncontaminated benchmark.

Calibration

• Amplitude: TABASCAL includes gain amplitude parameters, however, the errors in the estimated gain amplitudes remain at the same level as is provided in the prior distribution. Therefore, TABASCAL shows no amplitude calibration capability in the absence of no known calibrator sources. The direct inclusion of an astronomical source model is expected to lead to amplitude calibration capabilities, namely, the amplitude estimates would have lower error than given in the prior distribution.

• Phase: TABASCAL is able to leverage the RFI signal to constrain phase calibration solutions in a statistically consistent manner. Phase calibration constraints are shown to be directly proportional to the RFI S/N.

Stress testing

• Robustness: TABASCAL shows resilience under stress tests involving old TLE data and variations in prior hyperparameters, maintaining effective RFI mitigation and astronomical visibility recovery.

Overall, TABASCAL demonstrates remarkable, robust performance, achieving a level of visibility and image quality close to the uncontaminated ideal and outperforming traditional flagging methods in all scenarios for moderate and strong RFI contamination (S/N ≥ 1). For low levels of RFI contamination, all methods demonstrate an approximately equal level of performance. In future works, TABASCAL will be tested on real astronomical data and further computational improvements will be made. We wish to see TABASCAL become part of the standard data reduction pipeline in the future, helping to tackle the ever-growing impact of RFI, especially that of satellite mega-constellations.

Acknowledgements

We thank members of the SARAO Data Science team, Radio Astronomy Research Group at SARAO, Niruj Mohan and Robert Mueller for useful discussions. CF and MK acknowledge funding by the Swiss National Science Foundation. We also acknowledge the support of the South African Radio Astronomy Observatory. This research has been conducted using resources provided by the Science and Technology Facilities Council (STFC) through the Newton Fund and SARAO.

Data availability

The simulation and visibility recovery codes are available on their respective GitHub repositories tab-sim⁶ and tabascal⁷.

Appendix A Covariance estimation

We derive the method for scalable covariance estimation here. We first start by linearising our forward model of the visibilities V as a function of the parameters θ about the MAP point $\widehat{\mathit{\theta }}$. $$\mathit{V}\left(\mathit{\theta }\right)=\mathit{V}\left(\widehat{\mathit{\theta }}+\delta \mathit{\theta }\right)$$(A.1) $$=\mathit{V}\left(\widehat{\mathit{\theta }}\right)+\mathit{J}\delta \mathit{\theta }+𝓞\left(\delta {\mathit{\theta }}^2\right),$$(A.2)

and therefore by defining $\delta \mathit{V}=\mathit{V}\left(\mathit{\theta }\right)-\mathit{V}\left(\widehat{\mathit{\theta }}\right)+\mathit{J}\widehat{\mathit{\theta }}$, we get the linearised model as $$\delta \mathit{V}=\mathit{J}\mathit{\theta },$$(A.3)

where $\mathit{J}=\frac{\partial V}{\partial \theta }$ is the Jacobian. The prior distribution on θ is $$p\left(\mathit{\theta }\right)=𝓝\left({\mathit{\mu }}_{\Pi },{\mathbf{\Sigma }}_{\Pi }\right),$$(A.4)

where μ_Π and Σ_Π are the prior mean and covariance, respectively. The likelihood, assuming Gaussian noise with covariance Σ_Ν, is then $$p\left(\delta \mathit{V}|\mathit{\theta }\right)=𝓝\left(\mathit{J}\mathit{\theta },{\mathbf{\Sigma }}_N\right).$$(A.5)

Given that we are working with a linearised model, the posterior distribution is therefore a Gaussian distribution and is defined as follows $$p\left(\mathit{\theta }|\delta \mathit{V}\right)=𝓝\left(\widehat{\mathit{\theta }},\widehat{\mathbf{\Sigma }}\right),$$(A.6)

where $$\widehat{\mathbf{\Sigma }}={\left({\mathit{J}}^T{\mathbf{\Sigma }}_N^{-1}\mathit{J}+{\mathbf{\Sigma }}_{\Pi }^{-1}\right)}^{-1}$$(A.7)

and $$\widehat{\mathit{\theta }}=\widehat{\mathbf{\Sigma }}\left({\mathit{J}}^T{\mathbf{\Sigma }}_N^{-1}\delta \mathit{V}+{\mathbf{\Sigma }}_{\Pi }^{-1}{\mathit{\mu }}_{\Pi }\right).$$(A.8)

We never actually evaluate Equation (A.8) as we already have this point, the MAP $\widehat{\theta }$ from our optimisation. However, if we add Gaussian perturbations Δψ ~ 𝓝 (0, Σ_Ν) and Δϕ ~ 𝓝 (0, Σ_Π) to δV and μ_Π respectively, we obtain the following, $$\widehat{\mathit{\theta }}+\Delta \mathit{\theta }=\widehat{\mathbf{\Sigma }}\left({\mathit{J}}^T{\mathbf{\Sigma }}_N^{-1}\left(\delta \mathit{V}+\Delta \mathit{\psi }\right)+{\mathbf{\Sigma }}_{\Pi }^{-1}\left({\mathit{\mu }}_{\Pi }+\Delta \mathit{\varphi }\right)\right).$$(A.9) $$\Rightarrow \Delta \mathit{\theta }=\widehat{\mathbf{\Sigma }}\left({\mathit{J}}^T{\mathbf{\Sigma }}_N^{-1}\Delta \mathit{\psi }+{\mathbf{\Sigma }}_{\Pi }^{-1}\Delta \mathit{\varphi }\right).$$(A.10)

It can be easily shown that $\text{𝔼}\left(\widehat{\mathit{\theta }}+\Delta \mathit{\theta }\right)=\widehat{\mathit{\theta }}$ and $\text{𝔼}\left(\Delta \mathit{\theta }\Delta {\mathit{\theta }}^T\right)=\widehat{\mathbf{\Sigma }}$ as desired. Therefore, we now have a way to draw samples from our approximated posterior with which we can estimate covari-ances and marginal uncertainties easily. This process of drawing samples by Gaussian perturbations is taken from Papandreou & Yuille (2010).

Given that we have standardised our parameter space and the measurements have independent noise, as is done in Knollmüller & Enßlin (2019), we therefore have $${\mathbf{\Sigma }}_{\Pi }=1,$$(A.11) $${\mathbf{\Sigma }}_N={\sigma }_n^{2}1.$$(A.12)

We can see that obtaining the samples from ${\mathbf{\Sigma }}_N^{-1}\Delta \mathit{\psi }$ and ${\mathbf{\Sigma }}_{\Pi }^{-1}\Delta \Phi$ is both easy and scales linearly with the number of parameters N_P. The expensive part in terms of memory and computation comes from applying the Jacobian J and the inverse posterior information $\widehat{\mathbf{\Sigma }}$ at the MAP location. These are a Jacobian-vector product (JVP) and matrix-vector product (MVP) respectively and can be defined implicitly, namely, without the need to evaluate the Jacobian of size N_P × N_D and the posterior information of size N_p × N_P explicitly. Additionally, evaluating Equation (A.9) requires the inversion of the posterior information matrix, ${\widehat{\mathbf{\Sigma }}}^{-1}$, which would require $𝒪\left(N_P^3\right)$ computation typically. However, since it is by definition a symmetric operator, we can use the conjugate gradient method (Shewchuk et al. 1994) to apply its inverse to a vector. The benefit of conjugate gradient is we do not need to know the posterior information explicitly, and the inversion can be done in $𝒪\left(N_P^2\right)$ computation or less if an approximate solution is acceptable. Since the model is implemented in JAX (Bradbury et al. 2018), where the Jacobian-vector product is an integral part of the framework and implicit operators are easily defined, the implementation of the above sampling technique is therefore trivial. The error in the covariance estimate scales inversely to the number of samples taken, as expected. This type of implementation allows us to balance the required accuracy in the covariance estimate with the available computation while remaining scalable to millions of parameters.

Appendix B Simulation and prior parameters

References

All Tables

All Figures

Fig. 1 Power spectrum prior used for our astronomical visibilities (orange). The faint black curves are the true calculated power spectrum for all baselines in one of our simulations. The blue curve shows the median power spectrum and the black dashed line shows the maximum expected fringe rate of the longest baseline in our simulation. The parameters for the example priors shown are P0 = 107 Jy2, η0 = 1 mHz. and γ = 5, 2, and 1.

In the text

	Fig. 2 Antenna layout of the 32 MeerKAT antennas used in our simulations.
In the text

Fig. 3 Example predictions from TABASCAL for two different baselines (left and right) on an observation containing six satellites. The orange dotted curve shows the TABASCAL prediction and the blue curve shows the true value. The top panel shows this for the RFI visibility magnitude and the lower panels show the real part of the astronomical visibility. The top panel corresponds to a TABASCAL run where γ = 5, as defined in Equation (10). The three lower panels show the results when varying the prior parameter γ for the astronomical visibilities. Finally, γ controls the smoothness of the solutions.

In the text

Fig. 4 Errors in the astronomical visibility predictions from TABAS-CAL and the other cases for comparison as the S/N of the RFI is varied. The black curve shows the distribution of the visibility noise used in the observations and therefore corresponds to the errors in the uncontaminated case. In each panel five bins of RFI strength are used where multiple observations are bundled together. The coloured curves show the Gaussian fit to the error distributions, (i) The top panel shows the errors in the TABASCAL predicted visibilities, (ii) The middle and bottom panels show the errors from the AOFLAGGER and perfect 3σ flagging cases where flagged data is not included in the histograms.

In the text

Fig. 5 Images constructed from the same observation with our four different cases. Top left: No RFI contamination. Top right: TABASCAL fully removes the RFI and recovers the astronomical signal with comparable image noise to the uncontaminated data. Bottom left and right: the images from perfect 3σ flagging and AOFLAGGER respectively, showing significant striping due to residual RFI contamination. The mean RFI S/N in this data was 1.8, showing that significant issues occur for traditional flagging methods even at weak RFI. Here, a perfect 3σ flagging means that any RFI with true amplitude greater than 3× the noise is perfectly removed.

In the text

Fig. 6 Images constructed from the same observation with our four different cases. Top left: no RFI contamination. Top right: TABASCAL fully removes the RFI and recovers the astronomical signal with comparable image noise to the uncontaminated data. Bottom-left and right: the images from perfect 3σ flagging and AOFLAGGER respectively showing significantly higher image noise of 6× and 8× respectively. The AOFLAGGER image shows significant striped image artefacts which largely invalidates its use for science, with purity and completeness of around 20% and 50% respectively (see Figures 9 and 8). The mean RFI S/N in this data was 23. For large RFI amplitudes the quality differences between TABASCAL and the flagging methods increases further.

In the text

Fig. 7 Image noise calculated from the residuals after using WSCLEAN for imaging vs RFI S/N (lower x-axis) and ratio of RFI-to-Astronomical visibility strength (upper x-axis). The solid lines with dots represent the median within an RFI bin containing roughly 16 images each. The shaded region represents the 68% uncertainty interval over the observations in the bin. TABASCAL (orange) is statistically consistent with the uncontaminated case showing successful signal recovery despite the RFI contamination across the entire range of RFI strengths. In this simulation the astronomical visibilities are chosen to be close in amplitude to the visibility noise so the upper and lower x-axis scales are similar but not identical.

In the text

	Fig. 8 Completeness statistics of source catalogues across RFI S/N bins. Solid lines with dots indicate the median within a bin. The shaded region indicates the 68% uncertainty interval. Each bin consists of statistics from approximately 16 images for each case.
In the text

	Fig. 9 Purity statistics of source catalogues across RFI S/N bins. Solid lines with dots indicate the median within a bin. The shaded region indicates the 68% uncertainty interval. Each bin consists of statistics from approximately 16 images for each case.
In the text

	Fig. 10 Mean absolute flux error with respect to all matched sources across RFI S/N bins. Solid lines with dots indicate the median within a bin. The shaded region indicates the 68% uncertainty interval. Each bin consists of statistics from approximately 16 images for each case. TABASCAL (orange) is statistically consistent with the uncontaminated case.
In the text

Fig. 11 Evidence that TABASCAL achieves phase calibration off the strong RFI sources: A comparison of the root mean squared error (RMSE) in gain phase estimates and the posterior standard deviations from TABASCAL. TABASCAL shows statistically rigorous phase calibration capabilities that depend on the RFI S/N. Above an S/N of 1. phase constraints are inversely proportional to the RFI S/N. Below an S/N of 1, phase solutions are found within the prior information given.

In the text