Issue 
A&A
Volume 615, July 2018



Article Number  A66  
Number of page(s)  12  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/201732190  
Published online  17 July 2018 
The variance of radio interferometric calibration solutions
Qualitybased weighting schemes
^{1}
LESIA, Observatoire de Paris, PSL, CNRS, Sorbonne Université, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité,
5 place Jules Janssen,
92195
Meudon,
France
email: etienne.bonnassieux@obspm.fr
^{2}
Department of Physics & Electronics, Rhodes University,
PO Box 94,
Grahamstown
6140,
South Africa
^{3}
GEPI, Observatoire de Paris,
5 place Jules Janssen,
92190
Meudon,
France
^{4}
Station de Radioastronomie de Nançay (SRN, Observatoire de Paris, CNRS USN, PSL Research University, Université d’Orléans,
18330
Nançay,
France
^{5}
South Africa Radio Astronomy Observatory,
3rd floor, The Park, Park Roads,
7405
Pinelands,
South Africa
Received:
27
October
2017
Accepted:
6
February
2018
This paper investigates the possibility of improving radio interferometric images using an algorithm inspired by an optical method known as “lucky imaging”, which would give more weight to the bestcalibrated visibilities used to make a given image. A fundamental relationship between the statistics of interferometric calibration solution residuals and those of the imageplane pixels is derived in this paper. This relationship allows us to understand and describe the statistical properties of the residual image. In this framework, the noise map can be described as the Fourier transform of the covariance between residual visibilities in a new differential Fourier plane. Imageplane artefacts can be seen as one realisation of the pixel covariance distribution, which can be estimated from the antenna gain statistics. Based on this relationship, we propose a means of improving images made with calibrated visibilities using weighting schemes. This improvement would occur after calibration, but before imaging; it is thus ideally used between major iterations of selfcalibration loops. Applying the weighting scheme to simulated data improves the noise level in the final image at negligible computational cost.
Key words: instrumentation: interferometers / instrumentation: adaptive optics / methods: analytical / methods: statistical / techniques: interferometric / radio continuum: general
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0;), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Interferometers sample Fourier modes of the sky brightness distribution corrupted by instrumental and atmospheric effects rather than measuring the sky brightness directly. This introduces two problems for astronomers to invert: calibration and imaging. Both of these problems are illconditioned.
The problem of imaging consists of correcting for the incomplete uvcoverage of any given interferometer by deconvolving the instrument’s pointspread function (PSF) from images. Its poor conditioning comes from our limited a priori knowledge of the sky brightness distribution, combined with large gaps in our uvcoverage, which prevents us from placing strong constraints on image deconvolution. It can be betterconditioned in different ways, including through the use of weighting schemes (see Briggs 1995; Yatawatta 2014, and references therein) to improve image fidelity at the start of deconvolution. When inverting the imaging problem, we often assume that the sky is stable within the domain (i.e. is constant in time and frequency). There are exceptions, such as wideband deconvolution algorithms (e.g. Rau & Cornwell 2011) that explicitly take into account the sky’s frequency dependence, but still assume that the sky brightness distribution does not vary with time.
The problem of calibration is what concerns us in this paper. It consists of estimating and correcting for instrumental errors (which include effects such as antenna pointing errors, but also the phase delays caused by ionospheric activity, troposphere, and other factors). Calibration consists of solving for gain estimates, where a gain models the relationship between the electromagnetic field of an astrophysical source and the voltage that an antenna measures for this source. Because measurements are noisy, calibration often involves some finetuning of solution intervals to ensure that the solutions are wellconstrained while the solution intervals stay as small as signaltonoise allows. The calibration inverse problem involves three competing statistical effects: thermal noise in the measurements, true gain variability, and sky model incompleteness. If gain solutions are sought individually for each measurement, then calibration estimates will be dominated by thermal noise and will not adequately describe the actual gains. Similarly, if a single gain estimate is fitted to too many measurements, the intrinsic gain variability will be averaged out; for example, a choice of time and frequency interval that is too large will cause the solver to estimate a constant gain while the underlying function varies quickly, thereby missing much of the gain structure. This will introduce error, which will be correlated in time and frequency. This occurs, for example, when solving for ionospheric phase delays: in the most extreme case, where the solution interval is significantly larger than the scale of ionospheric fluctuations, its varying phase can average out to zero over the interval in time and frequency. Finally, if the model being fitted is incomplete, unmodelled physical flux will likely be absorbed unpredictably into both the gain solutions and the residual visibilities: this absorption of physical flux into gain solutions is known as source suppression (see Grobler et al. 2014; Kazemi & Yatawatta 2013, and references therein).
In practice, it is reasonable to assume that gain variation is generally slower than some given scale: we can then reduce the noise of our gain estimates by finding a single gain solution for a small number of measurements, assuming that the underlying gain variation is very small and stable over short intervals. This is generally a valid hypothesis, but the specific value for the variation scale can be contentious. Indeed, while the noise level can often be treated as constant throughout an observation, the gain variability itself is generally not constant: there will be timeperiods where the gains will tend to remain constant for longer, and others where variability will be very quick. This means that, for any choice of calibration interval, some gain estimates will be better than others, and almost all could have been improved (at a cost to others) by a different choice of time (and frequency) intervals.
Since we have measurements that are bettercalibrated than others (in that better estimates for their gains were obtained through chance alone), we could, in principle, take inspiration from “lucky imaging” (an opticaldomain method for making good images: for more details, see Fried 1978, and references therein) to weigh our visibilities according to their calibration quality. Those weights would in effect be an improvement of currently existing methods such as clipping noisy residual visibilities: in the extreme case where all visibilities are equallywell calibrated except a few which are extremely noisy, it should be equivalent to clipping. Otherwise, the weights should show at least a slight improvement over clipping.
The key finding of the present paper is a fundamental relationship between the covariance matrices of residual visibilities and the map of the covariance in the image plane: the CovCov relationship between visibility covariance to imageplane covariance. We show that the pixel statistics in the image plane are determined by a noise PSF, convolved with each source in the sky (modelled or not). This noise PSF is the product of the Fourier transforms of the gain covariance matrix with each cell mapped not from uv space to lm coordinates but rather between their respective covariance spaces: from a new differential Fourier plane (henceforth plane) to the imageplane covariance space δlδm. This imageplane covariance space describes the variance in each pixel and the covariance between pixels^{1}. It describes the expected calibration artefacts and thermal noise around each source, does not vary as a function of direction, and is convolved with each source in the field to yield the final error map. Because all unwanted (in our case, unphysical) signal can be thought of as noise, we will refer to the pixel variance map as the noise map.
The notion of a plane arises organically from the framework of radio interferometry: we are associating a correlation between visibilities to coordinates in covariance space, just as we associate the visibilities themselves to the uvdomain. The plane is the natural domain of these correlations. As previously stated, even if all sources in the field are perfectly known and modelled, a poor choice of calibration interval can introduce correlated noise in the residuals, which would then introduce larger variance near sources in the noise map. Conversely, if calibration is perfect, the noise map should be completely flat (i.e. same variance for all pixels), as there would be no noise correlation between pixels.
The main result of this paper consists in describing a new adaptive, qualitybased weighting scheme based on this insight. Using the CovCov relationship, we can create a new weighting scheme by estimating the residual visibility covariance matrix in a given observation. By weighting visibilities so as to change their covariance matrix, one can change the shape of the noise PSF and thus improve the final image: this manifests as either decreased noise or decreased calibration artefacts. This weighting is applied after calibration, but before image deconvolution; applying it will therefore not only improve the residual noise in the image, and thus the sensitivity achievable with a given pipeline, but will also improve deconvolution by minimising calibration artefacts in the field. It should thus effectively remove spurious, unphysical emission from final data products. Estimating the covariance matrix is the main difficulty of our framework; we do not know the underlying covariance matrix, and the conditioning of our estimation thereof is limited by the number of measurements within each solution interval. As such, we have no guarantee that our estimate of the corrected visibility covariance matrix is accurate. This problem can be alleviated, for example, by estimating the covariance matrix for the antenna gains themselves, and using this estimate to build the visibility covariance matrix: this effectively improves conditioning (cf. Sect. 4).
This paper is split into four main sections. In Sect. 2, we derive the CovCov relationship. With its newfound insights, we propose qualitybased weighting schemes with which to improve radio interferometric images in Sect. 3. We follow in Sect. 4 by showing how to estimate, from real data, the covariance matrix from which the qualitybased weights are derived. Our approach seems to give good results. Finally, we close the paper on a discussion of the applicability and limitations of the qualitybased weighting scheme.
2 Building the noise map
In this section, we derive our first fundamental result: the CovCov relationship, Eq. (25), which describes how the statistics of residual visibilities (and thus the antenna calibration solutions, henceforth gains) relate to the statistics of the image plane, that is, of images made using the associated visibilities. The dimensions of the matrices (denoted by boldface capital letters) and vectors (denoted by boldface lowercase letters) used in this paper are given in Table 1, along with the scalar numbers used to denote specific dimensions. All other variables are scalars.
2.1 The CovCov relationship in the δuδv plane
Let us begin by defining visibility gains. Using the Radio Interferometry Measurement Equation formalism for a sky consisting of a single point source (Hamaker et al. 1996; Smirnov 2011; and companion papers), we can write the following relation between the sky and the signal as measured by a single baseline at time t and frequency ν: (1)
All the quantities above are 2 × 2 matrices. Equation (1) implies a linear relationship between the coherency matrix and the visibilities recorded by a given baseline (), with the addition of a thermal noise matrix N, which is also of shape 2 × 2 and contains different complexvalued realisations of the noise in each cell. Since electric fields are additive, the sky coherency matrix can be described as the sum of the contributions from individual sources in directions d in the sky. We also assume that the sky does not vary over time, that is, that is not a function of time. The Jones matrices () contain the antenna gain information in matrix form, while is the Fourier kernel. Let us limit ourselves to the scalar case, which corresponds to assuming that emission is unpolarised. We assume that B^{d} = s_{d}I, where s is the flux of our single point source and I is the 2 × 2 identity matrix. We also assume that J_{p,tν} = g_{p,tν}I, where is the complexvalued gains of antenna p at time t and frequency ν. This means that we assume that the gains are directionindependent, and so becomes J_{...,tν}. Similarly, , the Fourier kernel in the direction of the source, d. The noise matrix N has one realisation of ϵ in each cell, where (2)
where σ is the variance of the thermal noise. Let us denote each pair by τ, and ignore the sky’s frequency dependence. The following scalar formulation is then equivalent to Eq. (1):
Calibration is the process of finding an accurate estimate of for all antennas p, at all times t and frequencies ν. Since we are in a directionindependent regime, the quality of our calibration then determines the statistical properties of the residual visibilities (and the imageplane equivalent, the residual image). The residual visibilities associated with calibration solutions are defined as our measured visibilities minus the gaincorrupted model visibilities. The variable then denotes our calibration estimate for . We now begin to limit the generality of our framework by assuming that all sufficiently bright sources have been modelled and subtracted: unmodelled flux is then negligible. We can then write the residual visibilities as (5)
The flux values in the imageplane pixels^{2} are the Fourier transform of the visibility values mapped onto each pixel. This can be written as
where lm are the directional cosine positions of a given pixel, and the Fourier coefficient mapping a point in Fourier space to a point on the imageplane. The scalar ω_{pq,τ} is the weight associated to a given visibility. Let us now write this using a matrix formalism. The contribution of a single visibility to the imageplane residuals can be written as
where ϵ is a vector of the ϵ of Eq. (2) and is a vector of size n_{pix}, with
and w_{b} is the scalar weight associated with each visibility. By default, : all visibilities then have the same weight, and then becomes the average of all . Here, is a vector of all , and thus of size n_{b}. Similarly, is the Fourier kernel, of size n_{pix} × n_{pix}. Finally, S_{b} is a matrix of size n_{pix}× n_{b}: its purpose is to encode the uvcoverage. Each S_{b} contains only a single nonzero cell, different for different S_{b}. The height (number of rows) of S_{b} is determined by the size of the uvgrid, and its length (number of columns) by the number of visibilities.
The order of operations is thus the following: each residual visibility is assigned some weight w_{b} and its uvcoordinates are set by . The inverse Fourier transform () is then applied to this grid, and so we recover its imageplane fringe. By averaging over all fringes, we recover the dirty image. The residual image will thus depend on three quantities: the residual gains, the flux in the image, and the weighting scheme. Let us consider the relationship between the statistics of residual visibilities and the variance at a given point in the corresponding residual image.
2.2 Statistical analysis
In the following analysis, we treat our gain solutions and thermal noise as random variables in order to compute the covariance matrix of our residual image, . The diagonal of this matrix gives the variance for each pixel, while the wings give the covariance between pixels. Using the property that Cov {Ax} = ACov{x}A^{H}, we can apply the Cov{} operator to Eq. (9) to write
So far, we have only applied definitions. The net effect of (dimensions of n_{pix}× n_{pix}) is to encode where a given baseline samples the uvplane, and map one cell at matrix coordinates (b, b′) from the correlation matrix onto the visibilitygrid. The matrix is not the gridding kernel, but rather the sampling matrix, which determines where we have measurements and where we do not. We can thus write that , where is the value from the appropriate cell and 1 is the vectorofones of appropriate length (here, n_{b}). This allows us to write
Here, is a Toeplitz matrix, that is, a convolution matrix, associated with baseline b. The set of all defines the convolution kernel that characterises the pointspread function (henceforth PSF) associated with a given uvcoverage, of size n_{pix} × n_{pix}. The matrix , meanwhile, is not generally Toeplitz. Its cells can be written as (19)
Let us investigate how the sky brightness distribution (i.e. ddependence) affects the noise map. We can write the sum over bb′ as two sums:one over b, b′ = b and one over b, b′≠b. Thus (20)
The only directiondependent terms in the above are s_{d} and , which are both inside (for both b = b′ and b≠b′). By making the approximation that the Fourier kernels of different sources are orthogonal (i.e. that )^{3} we can write
We note that , since those are the coordinates along the diagonal; for these values, the matrixofones at the centre of becomes the identity matrix. We note also that , since . We can then write Eq. (20) as (25)
where we have now limited our formalism to the case where the sky is dominated by distant pointlike sources. This is our fundamental result: assuming unpolarised emission coming from distant point sources and normallydistributed thermal noise, it gives a direct relationship between the covariance of the residual visibilities and the covariance of the residual imagepixel values. We thus call it the CovCov relationship. It describes the statistical properties of the image plane as the result of a convolution process changing an average noise level at different points in the image plane, allowing us to describe the behaviour of variance and covariance in the image. Let us focus on the first. By applying the operator (which returns the diagonal of an input matrix as a vector) to both sides of Eq. (25), we can find an expression for the variance map in the image plane:
In Eq. (27), we have (31)
For b≠b′, the diagonals of are the Fourier kernels mapping δuδv space to δlδm. Applying can then be thought of as performing a Fourier transform. It is not a diagonal matrix. It behaves as a covariance fringe, allowing us to extend standard interferometric ideas to covariance space. Each fringe can be thought of as a single “spatial filter” applied to the pixel covariance matrix. Just as a given baseline has coordinates in uvspace, a given correlation between baseline residual errors has coordinates in uv correlation space, which we will henceforth refer to as δuδvspace.
This δuδv space warrants further discussion. Figure 1 shows, for a given uvtrack (Fig. 1b), both the corresponding δuδv domain (Fig. 1c) and pointspread function (Fig. 1a). The symmetric, negative uvtrack is treated as a separate track, and thus ignored. This means that we do not fully constrain the noise PSF (since the covariance matrix of the symmetric track is simply the Hermitian of the first), but we do not seek to constrain it in this section, but rather to show that our results hold. We can see that the δuδvtracks are symmetrical about the origin. The δuδv space corresponding to a given uvtrack can thus be most concisely described as a “filled uvtrack”, with its boundaries defined by the ends of the uvtrack. The set of , each of which maps one value of the covariance matrix to a fringe in the image plane, would then characterise a PSF equivalent for the noise distribution, which we refer to as the noise PSF. In our formalism, the only source of statistical effects in the field are calibration errors and thermal noise. The average variance in all pixels will be given by the diagonal of the covariance matrix and the thermal noise, provided that it truly follows a normal distribution. The only effects which will cause the variance in the image plane to vary from one pixel to the next are those mapped onto the covariance fringes; such positiondependent variance fluctuations will be caused by correlated gain errors, which are spurious signal introduced by erroneous gain estimates. Assuming all sources in the field are pointlike and distant, then these variance fluctuations will follow a specific distribution, convolved to every source in the field. Since the variance fluctuations act as tracers for calibration artefacts, artefacts in the image can be understood as one realisation of the variance map, which is characterised by an average level determined by the variance in the gains and thermal noise, and a noise PSF convolved with the sky brightness distribution. The actual artefacts in the image will still be noisy, as a realisation of the true variance map. For the same reason, in the absence of correlated gain errors, are all zero and is a diagonal matrix. We then recover a “flat” noise map: the variance will be the same in all pixels, as the noise PSF is absent. In the ideal case, were we to recover the true value of the gains for all times and frequencies, this would become pure thermal noise.
Fig. 1 Panel a: PSF image of a simulated 1Jy source at phase centre. Colourbar units are in Jansky. Panel b: associated UV track, panel c: the corresponding tracks. The δuδv plane does not have a homogeneous point density, but is denser near its origin: here, this is shown by plotting only one random point in 10 000. 
2.3 Noisemap simulations
We have shown in Eq. (25) that there exists an analytical relationship between residual visibility statistics and imageplane residual statistics. This section gives details of simulations we have performed to support our claims about this “CovCov” relationship. Specifically, we simulate residual visibilities for a single baseline by generating a set of correlated random numbers with zero mean and a distribution following a specified covariance matrix C. It contains a periodic function of period T along the diagonal, which is then convolved with a Gaussian of width σ_{τ} corresponding to the characteristic scale of correlation. The values of these parameters are chosen arbitrarily. A small constant term is added on the diagonal, the net value of which is strictly positive. This simulates a low thermal noise. Finally, singular value decomposition is used to ensure that this matrix is Hermitian positive semidefinite. The net effect is a nonstationary correlation: some residuals are correlated with their neighbours, and uncorrelated with others. An example of this covariance matrix for arbitrary parameter values is shown in Fig. 2. We see that, for any given point, correlation is stronger with some neighbours than others (as determined by σ_{τ}). Samples are drawn as follows: C is built following user specifications as described above. We then find its matrix square root C_{0} so as to apply it to random numbers generated from a normal distribution. We generate 2000 realisations i of our random variables r_{i}:
where is a matrix of dimensions n_{realis} × n_{b}. Since x follows a normal distribution, Cov{x} = I and the covariance matrix of each is, by construction, . The covariance matrix of is thus also C.
As for the uvtrack, our simulations read a single one from a specified dataset. In this case, we read an eighthour uvtrack for a baseline between two arbitrary LOFAR (LOw Frequency Array) stations (specifically, CS001HBA0 and RS310HBA) in an observation of the Bootes extragalactic field. The effective baseline length varies between 37.9 km and 51.8 km. The dataset includes 20 channels, each with a spectral width of 97.7 kHz; the central observing frequency is 139 MHz. The temporal resolution is one measurement per second.
We compare the measured variance map , built by measuring the variance across realisations at each pixel in the imageplane, with the predicted variance map , built using the CovCov relationship (Eq. (25)). Since we are only interested in the variance map, rather than the covariance between pixels, we compute only the diagonal terms,
where the thermal noise is already incorporated into the diagonal of C and Eq. (36) is merely the diagonal operator applied to the CovCov relationship.
Fig. 2 Example of a nonstationary covariance matrix, which can be used to simulate . The colour bar units are Jy^{2}. The correlation scale σ_{τ} is 40 cells, and the variability period is 500 cells. The matrix is made positive semidefinite (and therefore a covariance matrix) through SVD decomposition. The maximum size of the “bubbles” is determined by σ_{τ}. 
2.3.1 Simulation with a single point source
We model our sky as containing a single 1 Jy point source at phase centre: we thus have . The source as seen through the set of uvtracks used in our simulation, along with their corresponding space, are shown in Fig. 1. The simulated noise map is calculated by drawing a large sample (n_{realis} = 2000) of random numbers from the correlated distribution, thereby creating 2000 sets of residual visibilities. By Fouriertransforming the visibilities to the image plane and taking the variance of the values for each image pixel (i.e. each l, m pair) as per Eq. (35), we can estimate . The predicted noise map, meanwhile, was found by assigning each cell of C to the appropriate point in the plane and Fourier transforming from this plane into the image plane, as per Eq. (36). We compare the outcome of simulating a large number n_{realis} of realisations and taking the variance across these realisations for each pixel with mapping the covariance matrix onto the δuδvplane and using the CovCov relationship.The results of our simulations are shown sidebyside in Fig. 4: the predicted and simulated noise PSFs match. The peaknormalised predicted noise PSF is less noisy, as shown in Fig. 3 for different correlation scales. This is expected, since it is calculated directly from the underlying distribution, rather than an estimate thereof. As n_{realis} →∞, we expect the two methods to fully converge. As the maximum characteristic correlation length σ_{τ} increases, the variance becomes ever more sharply peaked.
Since our simulated sky consists of a 1Jy source at phase centre, there is only one noise PSF to modulate the average noise map, and it lies at phase centre. Let us test our formalism further by considering a model with multiple point sources.
2.3.2 Simulation with three point sources
We wish to test our prediction that the noise map can be described as a convolutional process modulating an average noise level. We thus perform another simulation, this time with three 1 Jy point sources. The associated dirty image is shown in Fig. 5d.
This dirty image simply consists of performing a direct Fourier transform (i.e. without using a Fast Fourier Transform algorithm) on simulated visibilities corresponding to these three point sources. We now perform a similar test as above on this field. Firstly, we “apply” gain errors to these visibilities by multiplying our model with our residual gain errors. This allows us to find 2000 realisations of residual visibilities, and find the variance for each pixel across these realisations. This gives us the simulated noise map, shown in Fig. 5a. Secondly, we perform a Direct Fourier Transform (DFT) from the differential Fourier plane to the (l, m) plane as before, assigning one cell of to each point of the differential Fourier plane. This time, however, is not simply unity for all points in the differential Fourier plane. Instead, it is calculated for the threepointsource model, and applied for each point. This gives us the predicted noise map in Fig. 5b. Finally, Fig. 5c shows the absolute value of the difference between the two images. We see that there is some structure present in these residuals: this is expected, as the PSF of the sources in the dirty images clearly overlap. We are thus not quite in the regime where emission is fully spatially incoherent. Nevertheless, our predictions hold to better than 5% accuracy.
It bears repeating that, for correlated noise, this map can be understood as a distribution map for calibration artefacts: the amount of spurious correlated emission seen by each baseline will determine the noise map, and the true imageplane artefacts will then be one set of realisations of this underlying distribution.
3 Adaptive qualitybased weighting schemes
As discussed in Sect. 1, some intervals of an observation will have lower gain variability. These will show up in the gain covariance matrix as intervals with lower variance. Similarly, those with larger intrinsic gain variability will have greater error in their gainestimate. By giving greater weights to the former, and lower weights to the latter, we expect to be able to improve image reconstruction. We thus talk of adaptive qualitybased weighting, as the weights will adapt based on the calibration quality.
The pixel variance is determined by the visibility covariance matrix, as shown in Eq. (27). The diagonal of the visibility covariance matrix will add a flat noise to all pixels, while its wings will determine the calibration artefact distribution,which will be convolved to the sky brightness distribution. We thus have two sources of variance in the image plane. Minimising the farfield noise (i.e. the variance far from sources) in an image would involve downweighting noisier calibration intervals while upweighting the more quiescent ones, without taking noise correlation between visibilities into account. This is because the farfield noise will be dominated by the diagonal component of the covariance matrix (cf. Eq. (27)). By the same token, minimising calibration artefacts would involve downweighting measurements with stronglycorrelated noise, and upweighting the lesscorrelated. This would not, however, minimise the diagonal component: in fact, it will likely exaggerate its upweighting and downweighting. As such, it will increase the constant level of the noise map, but flatten the noisePSF’s contribution. There are thus two competing types of noise that we seek to minimise: uncorrelated noise, which corresponds to δuδv = 0 (i.e. the diagonal components of the gain covariance matrix), and correlated noise, which corresponds to δuδv≠0 (i.e. its wings). Minimising the first will minimise farfield noise without optimally reducing artefacts, while minimising the last will minimise noise near sources at a cost to farfield noise. In the following subsections, we will discuss weighting schemes used to accomplish this.
Fig. 3 Three lines in each figure corresponding to three horizontal crosssections from images in Fig. 4. The units on the yaxis are dimensionless [Jy^{2}∕Jy^{2}]. σ_{τ} is the maximum characteristic error correlation length. In decreasing intensity, they correspond to m = 0″, m = 4″, and m = 8″. The dashed lines correspond to the variance measured with 2000 realisations for each pixel, while the solid line correspondsto the predicted value at that pixel. There are 31 pixels. We do not show crosssections for negative m due to image symmetry about the origin. 
Fig. 4 Simulated noise maps, compared with theoretical prediction. The pixel values are normalised by the average value of the covariance matrix: the units of the colour bar are thus dimensionless (Jy^{2} ∕Jy^{2}). These are on the same angular scale as the source shown in Fig. 1. 
Fig. 5 Noise map of sky with correlated gain errors and three point sources. The colour bars of (a), (b), and (c) have dimensionless units, while that of (d) is in Jansky. We note the presence of structure in the residuals (c): these show the limits of ourhypothesis that sources are spatially incoherent. 
3.1 Optimising sensitivity
The CovCov relationship (Eq. (25)) tells us that, far from any sources, the variance map (Eq. (27)) is dominated by a constant term: the contribution from thermal noise and the diagonal of the residual visibility covariance matrix. Maximising sensitivity far from sources therefore implies minimising . This is equivalent to assigning visibilities weights inversely proportional to their variance: (37)
For each baseline, those times with larger variance in the residuals will be downweighted, and those with smaller variance will be upweighted; this scheme does not require information about the underlying gains, only the error on our solutions. Since we are treating as a constant for all antennas and all times, those times where our gains estimate is closer to the true gains will be upweighted, and those moments where they are farther from the actual gains will be downweighted: hence the term “adaptive qualitybased weighting”. We note that the diagonal of the weighted residuals’ covariance matrix should therefore become constant: this weighting scheme explicitly brings the residuals closer to what is expected in the case of perfect calibration, assuming uncorrelated noise. For the remainder of this paper, we will refer to these weights as sensitivityoptimal weighting.
3.2 Minimising calibration artefacts
Minimising calibration artefacts  optimising the sensitivity near bright sources  means flattening the noise map. Since the noise map can be understood as a noise PSF convolved with all the modelled sources in the sky modulating the background variance level, it will be flattest when its peak is minimised. From the CovCov relationship (Eq. (25)), we can see that, at the peak of the noise PSF (which would be the variance at the pixel where l = m = 0), the Fourier kernel is unity: the variance for that pixel is thus the sum of all the cells in the covariance matrix. By accounting for normalisation, we can write the variance at the centre of the noisePSF as (38)
Our optimality condition is then, after some algebra,
We find that one w which satisfies the above is (41)
where 1 is a vector of ones. These weights depend only on calibration quality: badlycalibrated cells will include spurious timecorrelated signal introduced by trying to fit the noise n on visibilities. Downweighting these cells helps suppress artefacts in the field, at the cost of farfield sensitivity. This weighting scheme is thus only a function of the relative quality of calibration at different times, boosting bettercalibrated visibilities and suppressing poorlycalibrated visibilities. For the remainder of this paper, we will refer to these weights as artefactoptimal weighting.
4 Estimating the covariance matrix
In our simulations, we have worked from a known covariance matrix and shown that our predictions for the residual image’s behaviour hold. With real data, however, we do not have access to this underlying covariance matrix. Since our weights are extracted from said matrix, estimating it as accurately as possible remains a challenge; this is in turn limited by the number of samples that can be used for each cell.
Each cell in the covariance matrix is built by averaging a number of measurements or samples. The more samples available, the better our estimate becomes. Once we have more samples than degrees of freedom, we say that our estimation is wellconditioned. Otherwise, it is poorlyconditioned. In this section, we will discuss ways in which we can improve the conditioning of the covariance matrix estimation.
4.1 Baselinebased estimation
One way to improve the conditioning of our covariance matrix estimation is to make the same hypothesis as the calibration algorithm: we can treat the underlying gains as constant within each calibration interval. Provided this interval is known, this allows us to find a single estimate for each interval block of the covariance matrix, turning a n_{b} × n_{b} matrix into a smaller n_{intervals} × n_{intervals} equivalent, where n_{intervals} is the number of solution intervals used to find the gain solutions. We then improve our conditioning by a factor of n_{int}, which is the number of samples in a calibration interval. The estimate of the covariance matrix is built by applying the covariance operator (42)
where the ⟨…⟩ operator denotes taking the average over the full vector. If the calibration solver’s gain estimates are unbiased (i.e. ) and the model of the sky is sufficiently complete, this quantity should be zero. Having created , which will be of size n_{b} × n_{b}, its cells can now be averaged over blocks of n_{int} × n_{int}. This allows us to estimate the weights for each baseline and each time.
Mathematically, we retrace the steps of Sect. 2. In the absence of directiondependent effects, we define the residual visibilities as before, and use them to define the normalised residual visibilities ρ_{b}:
We then organise the residuals in cells:
The matrix R contains all the residual visibilities within one calibration cell , that is, for where It is therefore of size , where n_{intervals} is the number of calibration intervals in the observation. Normalising the residual visibilities by k_{b} allows us to recover the underlying covariance matrix by multiplying the residual visibility matrix R with its Hermitian conjugate, (47)
We have divided the noise term by the flux model S_{b}, which can be very small in some cells. As such, care must be taken not to cause the relative thermal noise contribution to explode. Those cells where this would occur are dominated by thermal noise and information on the covariance matrix cannot be recovered from them.
In this framework, we simply treat the index as containing all the times and frequencies, for individual baselines, corresponding to a single calibration interval. The matrix is then an estimate of the residual visibility covariance matrix.
4.2 Antennabased estimation
In the subsection above, we assumed that finding one solution per interval would give us strong enough constraints to make the problem of estimating the covariance matrix wellconditioned: this may not be true in all cases. Conditioning may then need to be improved further: in this subsection, we show one way in which this can be done. There are others, for example using the rank of the matrix itself to find betterconditioned estimates of the covariance matrix at a lower resolution (i.e. a single estimate for a greater number of cells). They will not be presented in this paper, but are a possible avenue for future work on this topic.
In estimating the covariance matrix for each baseline and each calibration cell, we are severely limited by the small number of samples in each cell. One way to overcome this problem is to find estimates for the variance of antenna gains and use these to return to the baseline variances. In this formalism, we extend to include all visibilities pointing at a single antenna at a given time. Let us begin by writing an expression for the gain vector, which contains the gains for all antennas and all calibration cells,
and the variance on each antenna in each calibration cell is then (50)
As we can see, Eq. (50) is simply a vector form of Eq. (12). The residual gains of Eq. (12) can now be understood as random samples of the covariance between the gains for antennas p and q at a given time, assuming complete skymodel subtraction. We can thus define the variance sample matrix as an estimate of the variance matrix:
We define the residual matrix as (53)
where we explicitly place ourselves in the limits of our formalism, that is, we do not have directiondependent gains. We now see that at the core of Eq. (53) lies , where . The Kmatrix is defined as (54)
Since the residual matrix depends on the gains, we define the residual visibility vectors as
The matrix contains all the residual visibilities within one calibration cell , that is, for where Let us define as the number of elements in each calibration cell. The residual variance sample matrix can now be built by multiplying the residual visibility matrix with its Hermitian conjugate: (57)
We do this because it allows us to turn a single noise realisation ϵ into a statistical quantity σ. We can relate to the variance of individual antenna gains: (58)
To reach this point, in Eq. (23), we made the hypothesis that the sky brightness distribution is dominated by spatially incoherent emission. Applying this hypothesis again here, we can make the approximation that the cross terms in the sum over d, d′ average to zero: . We then have:
where ° denotes the Hadamard or entrywise product and . Thus, allows us to estimate the variance of each antenna and for each calibration cell by using all the visibilities pointing to that antenna within that calibration cell. With this information, we can then rebuild the baselinedependent matrix, having improved our sampling by a factor of n_{ant}.
5 Applying the correction to simulated data
In this section, we show the impact of our weighting schemes on a noise map made from arbitrarily stronglycorrelated residuals. Here, we assume that our sky contains only a single point source at phase centre: there is thus only a single instance of the noise PSF, placedat phase centre, to modulate the average variance level. We sample this instance by taking a crosssection from to a large l, keeping m constant. The only difference between these crosssections is the weighting scheme applied: unit weights for all visibilities (“Uncorrected”, blue), sensitivityoptimal weights (green), and artefactoptimal weights (red). We plot both the result predicted by the CovCov relationship (solid line) and the variance estimated across 2000 realisations (dashed line): the result is shown in Fig. 6. The two remain in such agreement throughout the crosssection as to be nearly indistinguishable.
There are a few significant points to note on this figure. Firstly, most of the power in the matrix lies along the diagonal; both weighting schemes thus give good improvements in variance across the map. The artefactoptimal weights, while decreasing the peak further, as expected, also increase the noise far from sources. This is due to the fact that the artefactoptimal weights are in a sense more “selective” than the sensitivityoptimal weights; they upweigh and downweigh more severely, and will only result in a constant covariance matrix if this matrix is zero everywhere outside of the diagonal. In effect, the noise map becomes flatter, but much broader.
Fig. 6 Sensitivityoptimal (green) and artefactoptimal (red) weights both give improvements over the unweighted noise map (blue). 
6 Applying the correction to real data
In this section, we show the effect of adaptive qualitybased weighting on real data. The dataset used in this section is a single subband from an eighthour LOFAR observation centred on the Extended Groth Strip (α = 14:19:17.84, δ = 52:49:26.49). The observation was performed on August 28, 2014. The subband includes eight channels of width 24.414 kHz each, for a total bandwidth ranging from 150.2 to 150.5 MHz. The data have been averaged in time to one data point per second. The data were calibrated using Wirtinger calibration (see Tasse 2014; Smirnov & Tasse 2015, and references therein) and a sky model consisting only of a nearby calibrator source, 3C295. A reference image (a cutout of which is shown in Fig. 7a) was made by calibratingthe data according to best practice for LOFAR survey data: one calibration solution per eight seconds and per four channels. The residual data was then corrected by the gain solutions and imaged using Briggs weighting (robust = 0), pixel size of 1.5″, and deconvolved using the default deconvolution algorithm in DDFacet (see Tasse et al. 2018).
The data was then timeaveraged to create a new, 2.4 GB dataset with one data point per eight seconds. Deliberatelypoor calibration was then performed on this dataset, solving for one calibration solution per two minutes (caeteris paribus). The resulting corrected residual data was imaged using the same imaging parameters as the reference image, and a cutout of the result is shown in Fig. 7b. As expected, the very long calibration intervals introduce calibration artefacts into the image. The brightest sources are still visible, but much of the fainter emission is buried under these artefacts. We then have the case where our residual visibilities are dominated by calibration error rather than sky model incompleteness.
Weightswere then calculated based on the badlycalibrated residual visibilities. Figure 7c was made using the same visibilities as Fig. 7b and applying baselinebased, sensitivityoptimal weight. Similarly, Fig. 7d used the poorlycalibrated residual visibilities with the application of baselinebased, artefactoptimal weighting. These weights are likely to be poorlyconditioned. In both cases, all other parameters were conserved.
Applying antennabased sensitivityoptimal weighting to the badlycalibrated data (not shown here) allows us to recover the reference image with only a very small increase in rms (increased by a factor of 1.14). Further testing on complex field simulations will be required to ascertain the usefulness of artefactoptimal weighting; it is likely that it fails to correct the image fully due to the poor conditioning of the covariance matrix used here.
The pixel histograms show us that the weights do not completely mitigate the poor calibration interval choice, but certainly give a dramatic improvement over the unweighted, poorlycalibrated residuals. This is compatible with our statement that the weights give similar residuals in the image with a dramatic improvement in time at some cost to sensitivity. It is interesting to note that while Fig. 7d looks noisier than Fig. 7c, its residual flux histogram is actually closer to that of Fig. 7a.
As for performance, the weights used for Fig. 7d took eight hours of computing time on a single core^{4}, working on a 29 GB dataset, which is not particularly large for LOFAR data. Since the problem is massively parallel, this cost can be alleviated. The main bottleneck is likely due to very poor code optimization. As for the weights used for Fig. 8b, they are computed in 1min6s on the same single core.
7 Discussion
This paper began by investigating the use of an algorithm inspired by “lucky imaging” to improve images made using radio interferometric data. By investigating the statistics of residual visibilities, we have made the following findings:
A relationship between the statistics of residual visibilities and residual pixel values (the “CovCov relationship”).
A description of the noise map in the image plane as a constant variance level modulated by a noise PSF convolved with the sources in the field. This gives the variance in the flux of the image as a function of distance from the sources in the sky for a given calibration.
An adaptive qualitybased weighting scheme, which reduces the noise in the image (and the presence of calibration artefacts) by minimising either the constant noise term or the noise PSF.
While our results are not a panacea for poor calibration, they show that we can not only improve images made with wellcalibrated data, but also mitigate the worst effects of poorlycalibrated visibilities in otherwise wellcalibrated datasets. Provided that the gain variability timescale is long enough at certain points of the observation, we can effectively get images of similar quality using both the “standard” bestpractice calibration interval for LOFAR survey data (calibration solution interval of eight seconds) and a significantly larger solution interval of two minutes (frequency interval unchanged). Of course, if no such stable interval exists, there will be no good intervals to upweigh, and we will be left only with equallypoor data chunks. This means that, in the right conditions, net pipeline time can be sped up by a factor of nearly three, at a slight cost to sensitivity. This increase will be greater than what could be achieved with existing comparable methods such as “clipping”.
We emphasise that the adaptive qualitybased weighting schemes work because the noise map describes the underlying noise distribution, of which calibration artefacts are one single realisation. To fully characterise the artefacts, the correlation between different pixels (i.e. offdiagonal elements of ) must be computed; this has not been done in this paper. Nevertheless, lesser constraints on the spatial distribution of artefacts can be found using only the diagonal elements of . The weighting schemes merely seek to minimise this spatial distribution as much as possible; the end result is fewer artefacts, which can be distributed across a much larger area. This is the source of the dramatic improvement from Figs. 7b and c. We have simply downweighted those visibilities where spurious signal was introduced by the calibration solutions, and upweighted those visibilities where such signal was lesser. Since this spurious signal is the source of calibration artefacts, downweighting the associated visibilities reduces it dramatically. The poor improvement from Figs. 7b and d is likely due to limits in the conditioning of our estimation of the covariance matrix.
The work presented here can be improved upon, notably by working on improving the conditioning of our covariance matrix estimate; for real observations, it is impossible to have more than one realization of each gain value for all antennas. By treating each visibility within a calibration interval as a realization of the true distribution, we can better estimate the covariance matrix per baseline, and thus reach a better estimate of the variance in the image plane. Of course, in practice, we can never access the true, underlying timecovariance matrix for each baseline. Significant hurdles remain:
The impact of sky model incompleteness (since calibration requires a sky model) is ignored in this paper; we start by assuming that we have a complete sky model. In practice, of course, acquiring a complete sky model is often a key science goal in and of itself. The impact of this hypothesis therefore ought to be investigated in future work.
The conditioning of our covariance matrix estimation remains a concern. By using an antennabased approach, we can improve conditioning by a factor of n_{ant}, but this is only one approach among many. Further work is needed to investigate which method, if any, proves optimal.
Fig. 7 Restored images of the centre of the Extended Groth Strip, as seen with an eighthour observation using the full LOFAR array. Panel a: shows an image of the field made with good calibration intervals; wellcalibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Used for comparison with the other images. Units of colour bar are Janskys. This image was made with data calibrated following best practice (solution intervals of 8 s, half the bandwidth). rms = 5.87 mJy beam^{−1}. Panel b:shows an image of the field made with poor calibration intervals; poorlycalibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Units of colour bar are Janskys. This image was made with the same data as for panel a, but averaged in time and calibrated using larger gain solution intervals: two minutes and half the bandwidth. rms = 86.4 mJy beam^{−1}. Panel c: shows image made with the same visibilities and imaging parameters, but with the application of the sensitivityoptimal weighting scheme; image made using the same imaging parameters and corrected visibilities as panel b, using sensitivityoptimal weighting. Units of colour bar are Janskys. rms = 9.69 mJy beam^{−1}. Panel d: similarly, differs from panel c only in that artefactoptimal weights, rather than sensitivityoptimal weights, were used; image made using the same imaging parameters and corrected visibilities as panel b, using artefactoptimal weighting. Units of colour bar are Janskys. rms = 15.8 mJy beam^{−1}. The histograms of pixel values in each image have 1000 flux bins ranging from –0.16 Jy to 0.16 Jy. Their ordinates are in log scale. Pixel size is 1.5′′ in all images. 
Fig. 8 Comparison between the wellcalibrated image (i.e. the same image as Fig 7a) and antennabased sensitivityoptimal weights. Units of both colour bars are Jansky. Panel a: wellcalibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Used for comparison with the other images. Calibration solution intervals used were eight seconds, half the bandwidth. rms = 5.87 mJy beam^{−1}. Panel b: image made using the same imaging parameters and corrected visibilities as Fig. 7b, with the application of antennabased, sensitivityoptimal weighting. Solution interval of two minutes, half the bandwidth. rms = 6.69 mJy beam^{−1}. We see that we recover a very similar image, despite the fact that the data used for the weighted image are averaged by a factor of 8compared to those used for the unweighted image. 
Acknowledgements
This work is based upon research supported by the South African Research Chairs Initiative of the Department of Science and Technology and National Research Foundation. The authors thank the LOFAR Surveys KSP Primary Investigators for allowing us to use their data for our tests and images in this paper. They would also like to thank M. Atemkeng and T. Gobler for fruitful discussions over the course of this work. They helped shape this work through their insightful comments on both style and substance while proofreading.
References
 Briggs, D. S. 1995, in Amer. Astron. Soc. Meet. Abstr., BAAS, 27, 1444 [Google Scholar]
 Fried, D. L. 1978, J. Opt. Soc. Am. (19171983), 68, 1651 [NASA ADS] [CrossRef] [Google Scholar]
 Grobler, T. L., Nunhokee, C. D., Smirnov, O. M., van Zyl, A. J., & de Bruyn A. G. 2014, MNRAS, 439, 4030 [NASA ADS] [CrossRef] [Google Scholar]
 Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kazemi, S., & Yatawatta, S. 2013, MNRAS, 435, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Rau, U., & Cornwell, T. J. 2011, A&A, 532, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smirnov, O. M. 2011, A&A, 527, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smirnov, O. M., & Tasse, C. 2015, MNRAS, 449, 2668 [NASA ADS] [CrossRef] [Google Scholar]
 Tasse, C. 2014, ArXiv eprints [arXiv: 1410.8706] [Google Scholar]
 Tasse, C., Hugo, B., Mirmont, M., et al. 2018, A&A, 611, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yatawatta, S. 2014, MNRAS, 444, 790 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Panel a: PSF image of a simulated 1Jy source at phase centre. Colourbar units are in Jansky. Panel b: associated UV track, panel c: the corresponding tracks. The δuδv plane does not have a homogeneous point density, but is denser near its origin: here, this is shown by plotting only one random point in 10 000. 

In the text 
Fig. 2 Example of a nonstationary covariance matrix, which can be used to simulate . The colour bar units are Jy^{2}. The correlation scale σ_{τ} is 40 cells, and the variability period is 500 cells. The matrix is made positive semidefinite (and therefore a covariance matrix) through SVD decomposition. The maximum size of the “bubbles” is determined by σ_{τ}. 

In the text 
Fig. 3 Three lines in each figure corresponding to three horizontal crosssections from images in Fig. 4. The units on the yaxis are dimensionless [Jy^{2}∕Jy^{2}]. σ_{τ} is the maximum characteristic error correlation length. In decreasing intensity, they correspond to m = 0″, m = 4″, and m = 8″. The dashed lines correspond to the variance measured with 2000 realisations for each pixel, while the solid line correspondsto the predicted value at that pixel. There are 31 pixels. We do not show crosssections for negative m due to image symmetry about the origin. 

In the text 
Fig. 4 Simulated noise maps, compared with theoretical prediction. The pixel values are normalised by the average value of the covariance matrix: the units of the colour bar are thus dimensionless (Jy^{2} ∕Jy^{2}). These are on the same angular scale as the source shown in Fig. 1. 

In the text 
Fig. 5 Noise map of sky with correlated gain errors and three point sources. The colour bars of (a), (b), and (c) have dimensionless units, while that of (d) is in Jansky. We note the presence of structure in the residuals (c): these show the limits of ourhypothesis that sources are spatially incoherent. 

In the text 
Fig. 6 Sensitivityoptimal (green) and artefactoptimal (red) weights both give improvements over the unweighted noise map (blue). 

In the text 
Fig. 7 Restored images of the centre of the Extended Groth Strip, as seen with an eighthour observation using the full LOFAR array. Panel a: shows an image of the field made with good calibration intervals; wellcalibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Used for comparison with the other images. Units of colour bar are Janskys. This image was made with data calibrated following best practice (solution intervals of 8 s, half the bandwidth). rms = 5.87 mJy beam^{−1}. Panel b:shows an image of the field made with poor calibration intervals; poorlycalibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Units of colour bar are Janskys. This image was made with the same data as for panel a, but averaged in time and calibrated using larger gain solution intervals: two minutes and half the bandwidth. rms = 86.4 mJy beam^{−1}. Panel c: shows image made with the same visibilities and imaging parameters, but with the application of the sensitivityoptimal weighting scheme; image made using the same imaging parameters and corrected visibilities as panel b, using sensitivityoptimal weighting. Units of colour bar are Janskys. rms = 9.69 mJy beam^{−1}. Panel d: similarly, differs from panel c only in that artefactoptimal weights, rather than sensitivityoptimal weights, were used; image made using the same imaging parameters and corrected visibilities as panel b, using artefactoptimal weighting. Units of colour bar are Janskys. rms = 15.8 mJy beam^{−1}. The histograms of pixel values in each image have 1000 flux bins ranging from –0.16 Jy to 0.16 Jy. Their ordinates are in log scale. Pixel size is 1.5′′ in all images. 

In the text 
Fig. 8 Comparison between the wellcalibrated image (i.e. the same image as Fig 7a) and antennabased sensitivityoptimal weights. Units of both colour bars are Jansky. Panel a: wellcalibrated, unweighted restored image of the sky near the centre of the Extended Groth Strip. Used for comparison with the other images. Calibration solution intervals used were eight seconds, half the bandwidth. rms = 5.87 mJy beam^{−1}. Panel b: image made using the same imaging parameters and corrected visibilities as Fig. 7b, with the application of antennabased, sensitivityoptimal weighting. Solution interval of two minutes, half the bandwidth. rms = 6.69 mJy beam^{−1}. We see that we recover a very similar image, despite the fact that the data used for the weighted image are averaged by a factor of 8compared to those used for the unweighted image. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.