15.4 A CLEAN-based algorithm for mosaic deconvolution

 

The dirty mosaic

The dirty maps of each field i are computed with the same phase center (i.e. the same coordinate system) and can therefore be written:

$$F_i = D_i * (B_i \times I) + N_i$$ (15.10)

Note that the dirty beams D_i are a priori different for each pointing center, because the uv coverages, even if similar, are slightly different. The dirty mosaic J can then be constructed according to equation 15.9:

$$J = \frac{\displaystyle \sum\nolimits_i \frac{B_i}{\sigma_i^2... ...+N_i\Big]}{\displaystyle \sum\nolimits_i B_i^2\,\sigma_i^{-2}}$$ (15.11)

This relation is homogeneous to the sky brightness distribution I: the mosaic is corrected for the primary beams attenuation. In practice, a slightly modified mosaic is computed, in order to avoid noise propagation (it makes no sense to add to the center of a field noise coming from the external, attenuated regions of an adjacent field). For that purpose, the primary beams used to construct the mosaic are truncated to some value, typically 10 to 30% of the maximum. The mosaic is thus defined by:

$$J = \frac{\displaystyle \sum\nolimits_i \frac{B_i^t}{\sigma_i... ...\Big]}{\displaystyle \sum\nolimits_i {B_i^t}^2\,\sigma_i^{-2}}$$ (15.12)

where B_i^t denotes the truncated primary beam of the field i. This relation is the measurement equation of a mosaic, connecting the observed quantity J to the sky brightness distribution I (equation 15.1 was the measurement equation of a single-field observation).

Noise distribution

  

Figure: One-dimensional mosaic of 10 half-power overlapping fields, with identical noise level $\sigma$. (Lower panel:) Normalized primary beams, truncated to B $_{\rm min} = 0.1$. (Upper panel:) Resulting noise distribution (Eq.15.14). The noise rms in the mosaic is roughly constant, about 20% lower than the noise of each individual field, but strongly increases at the edges. The two thick vertical lines indicate the truncation of the mosaic done by the algorithm at $\sigma _J=\sigma /\sqrt {\rm B_{min}}$.

Due to the correction for the primary beams attenuation, the noise distribution in a mosaic is not uniform. From Eq.15.12, it can be written:

$$N = \frac{\displaystyle\sum\nolimits_i B_i^t\,\sigma_i^{-2}\,N_i} {\displaystyle\sum\nolimits_i {B_i^t}^2\,\sigma_i^{-2}}$$ (15.13)

Accordingly, the rms depends on the position and is given by:

$$\sigma_J = \frac{\sqrt{\displaystyle\sum\nolimits_i {B_i^t}^2... ...sqrt{\displaystyle\sum\nolimits_i {B_i^t}^2 \, \sigma_i^{-2}}}$$ (15.14)

The noise thus strongly increases at the edges of the mosaic (see Fig. 15.1). The non-uniformity of the noise level with the position makes it impossible to use classical CLEAN methods to deconvolve the mosaic: the risk to identify a noise peak as a CLEAN component would be too important. It is thus necessary to identify the CLEAN components on another distribution. For that purpose, the ``signal-to-noise'' distribution is computed:

 

$$\frac{J}{\sigma_J} = \frac{\displaystyle\sum\nolimits_i B_i^t \, ... ...a_i^{-2}\, F_i}{\sqrt{\displaystyle\sum\nolimits_i {B_i^t}^2 \, \sigma_i^{-2}}}$$ H =

$${\rm i.e.:}\hspace{1cm} H \frac{\displaystyle\sum\nolimits_i B_i^t \, \sigma_i^{-2} \, \Big... ... I) + N_i\Big]}{\sqrt{\displaystyle\sum\nolimits_i {B_i^t}^2 \, \sigma_i^{-2}}}$$ = (15.15)

Deconvolution algorithm

The main idea of the algorithm is to iteratively find the positions of the CLEAN components on H, and then to correct the mosaic J. The initial distributions J₀ and H₀ are computed from the observations and the truncated primary beams, according to equations 15.12 and 15.15. The following operations have then to be performed at each iteration k:

1.

Find the position (x_k,y_k) of the maximum of H.

2.

Find the value j_k of J at the position (x_k,y_k).

3.

Remove from J the contribution of a point-like source of intensity $\gamma j_k$, located at (x_k,y_k) ($\gamma$ is the loop gain, as in the normal CLEAN algorithm):

$$J_{k} = J_{k-1} - \frac {\displaystyle\sum\nolimits_i B_i^t \... ...]\Big]}{\displaystyle\sum\nolimits_i {B_i^t}^2\,\sigma_i^{-2}}$$ (15.16)

$\delta(x_k,y_k)$ denotes a Dirac peak located in (x_k,y_k).

4.

Do the same for H: remove the contribution of a point-like source of intensity $\gamma j_k$, located at (x_k,y_k):

$$H_k = H_{k-1} - \frac {\displaystyle\sum\nolimits_i B_i^t \,\... ...{\sqrt{\displaystyle\sum\nolimits_i {B_i^t}^2\,\sigma_i^{-2}}}$$ (15.17)

Note that in the two last relations, the CLEAN component is multiplied by the true, not truncated primary beam (taken at the (x_k,y_k) position).

After $k_{\rm max}$ iterations, the mosaic J can therefore be written:

$$J = \frac{\displaystyle \sum\nolimits_i B_i^t\,\sigma_i^{-2}\... ...tyle \sum\nolimits_i {B_i^t}^2\sigma_i^{-2}} + J_{k_{\rm max}}$$ (15.18)

Enough iterations have to be performed to ensure that the residual $H_{k_{\rm max}}$is smaller than some user-specified threshold (typically 1 to 3). The comparison between Eqs.15.12 and 15.18 shows that, within the noise, the sum of the CLEAN components can be identified with the sky brightness distribution I. As with the normal CLEAN algorithm, the final clean image is then reconstructed as:

$$M = C * \left[ \sum_{k=1}^{k_{\rm max}} \gamma\,j_k\,\delta(x_k,y_k) \right] + J_{k_{\rm max}}$$ (15.19)

where C is the chosen clean beam. The modified CLEAN algorithms proposed e.g. by [Clark 1980] or [Steer et al 1984] can be similarly adapted to handle mosaics, the main idea being to identify CLEAN components on H and to correct J. Note however that the multi-resolution CLEAN [Wakker & Schwartz 1988] cannot be directly adapted, as it relies on a linear measurement equation, which is not the case for a mosaic.

The MAPPING software

MAPPING is an superset of the GRAPHIC software, which has been developed to allow more sophisticated deconvolutions to be performed. For instance, it allows to choose a support for the deconvolution (clean window) or to monitor the results of the deconvolution after each iteration. Several enhancements of CLEAN (e.g. multi-resolution CLEAN) as well as the WIPE algorithm (see lecture by E. Anterrieu) are also available. The deconvolution of a mosaic has to be done with MAPPING. The implemented algorithm assumes that the noise levels in each field are similar (i.e. $\forall i\,\, \sigma_i = \sigma$), which is a reasonable hypothesis for Plateau de Bure observations. In that case, the equations of the previous paragraph are slightly simplified: J is independent from $\sigma$, and H can be written as the ratio $H'/\sigma$, where H' is independent from $\sigma$ and is actually used to localize the CLEAN components.

We refer to the Mapping Cookbook for a description of the MAPPING software. In short, to deconvolve mosaics:

$\circ$

Create a uv table for each observed field. Then, run the UV_MAP task to compute a dirty map and a dirty beam for each field, with the same phase center (variable UV_SHIFT = YES).

$\circ$

The task MAKE_MOSAIC is used to combine the fields to construct a dirty mosaic. Two parameters have to be supplied: the width and the truncation level $B_{\rm min}$ of the primary beams. Three images are produced: the dirty mosaic^15.3 (yourfile.lmv), all the dirty beams written in the same file (yourfile.beam) and a file describing the positions and sizes of the primary beams (yourfile.lobe). The dirty maps and beams of each individual field are no longer used after this step and can thus be removed if necessary.

$\circ$

The data have to be loaded into the MAPPING buffers. This is done by the READ DIRTY yourfile.lmv, READ BEAM yourfile.beam, and READ PRIMARY yourfile.lobe commands. The latter automatically switches on the mosaic mode of MAPPING (the prompt is now MOSAIC>). From now, the deconvolution commands HOGBOM, CLARK and SDI can be used and will apply the algorithm described above. Use the command MOSAIC to switch on or off the mosaic mode if necessary.

$\circ$

The clean beam of the final image can be specified by the user (variables MAJOR and MINOR). Otherwise, the clean beam computed from the first field will be used. To check if there are differences between the various dirty beams, just use the FIT i command, which will indicates the clean beam computed for the ith field.

$\circ$

The deconvolution will use the same parameters as a usual CLEAN: support, loop gain, maximal number of iterations, maximal value of the final residual, etc.

$\circ$

In addition, two other parameters, SEARCH_W and RESTORE_W, can be supplied. Due to the strong increase of the noise at its edges, the mosaic has to be truncated above some value of $\sigma_J$, and these two variables are used to define this truncation level (in terms of $(\sigma_J/\sigma)^{-2}$). More precisely, SEARCH_W indicates the limit above which CLEAN components have not to be searched, while RESTORE_W indicates the limit above which the clean image is not reconstructed. Default values of these two parameters (both equal to $B_{\rm min}$) are strongly recommended. The corresponding truncation is shown in Fig. 15.1.

Tests of the method

Several tests of the method described in this paragraph have been performed, either with observations (including the comparison of independent mosaics from the same source) or with simulations. They show that very satisfactory results can be achieved with typical Plateau de Bure observations. Interestingly, MEM deconvolution of the same data set (using the task VTESS in AIPS) seems to give worse results: this is most probably related to the limited uv-coverage obtained with the Plateau de Bure interferometer, as compared to typical VLA observations (MEM is known to be vulnerable when there is a relatively small number of visibilities).