基于幅度的CLEAN天文图像重建算法

An Amplitude-Based CLEAN Astronomical Image Reconstruction Algorithm

摘要: 射电干涉成像的质量部分依赖于图像重建算法。CLEAN算法是一种应用非常广泛的射电天文数据的重建算法，已经成为射电天文干涉数据处理软件的标准配置。典型的CLEAN算法使用delta函数集合来逼近我们真实的天文图像，然而当前没有工作来优化这种集合，从而达到提高图像重建效果的目的。本文改进了当前的CLEAN算法，提出了一种基于幅度的重建算法。通过通用天文软件程序包(Common Astronomy Software Applications, CASA)模拟美国甚大干涉阵(VLA)进行了观测，对得到测量数据MS (Measurement Sets)使用常规的CLEAN算法和本文改进的CLEAN算法进行了重建。实验显示本文改进的算法的重建图像有更高动态范围。

Abstract: The quality of radio interference imaging partly depends on the image reconstruction algorithm. The CLEAN algorithm is a widely used reconstruction algorithm for radio astronomy data and has become the standard configuration of radio interferometric data processing software. A typical CLEAN algorithm uses a set of delta functions to approximate a real astronomical image. However, currently there is no work to optimize such a set to achieve the purpose of improving the image reconstruction effect. This paper improves the current CLEAN algorithm and proposes an amplitude-based reconstruction algorithm. The VLA observation was simulated by the Common Astronomy Software Applications (CASA), and the conventional CLEAN algorithm and the improved CLEAN algorithm were used to reconstruct the measured data MS (Measurement Sets). Experiments show that the improved algorithm has a higher dynamic range of the reconstructed image.

1. 引言

干涉阵和甚长基线干涉VLBI利用多个望远镜联合测量，通过干涉技术可以获得一个很大的综合孔径(synthesis aperture)，这样大大提高了设备的分辨率，如VLA¹，ALMA²和SKA³。理论上地基的VLBI网络的综合孔径可以接近地球的直径，其分辨率是单望远镜无法达到的。由于望远镜的数量和地理位置限制，干涉阵和VLBI测量在大多数情况下无法满足奈奎斯特采样定律。Van Citter-Zernike理论 [1] 指出，真实天文图像空间与干涉测量的空间存在傅立叶变换关系。数据在空间频率域的不完全测量相当于在空间域引入一个模糊核函数。得到的模糊图像是真实图像与模糊核函数卷积的结果。在射电天文领域，这个模糊核函数常常被称为“脏束”(dirty beam)，模糊图像被称为“脏图”(dirty image)，测量的数据被称为“可见度函数”(visibility)。

脏图是不利于科学研究的。在很多情况下，当脏束的旁瓣水平很高的时候，脏图就基本无法用于一些科学研究。所以重建技术对于干涉测量数据的恢复是非常重要的。由于这种模糊是卷积造成的，所以重建干涉测量数据的方法也常常被称为反卷积。在当前，有三类方法来解这个反卷积：CLEAN算法，最大熵方法，压缩感知重建技术。最大熵方法 [2] 使用显式的规则方法来使得图像以正则项的方式来逼近真实图像，图像随着正则方式不同而改变。压缩感知重建技术 [3] 将压缩感知理论应用到图像重建，它是最近几年发展起来的一项新技术，也取得了一些不错的进展 [4] [5] [6]。CLEAN算法是Högbom在1974年的发明 [7]，它的基本思想是用delta函数集合来逼近真实的天文图像。后来Schwarz引进了数学理论来解释CLEAN算法 [8]；Clark通过快速傅立叶变换和波束块(beam patch)引进了快速算法 [9]；Cotton和Schwab将CLEAN算法扩展到测量数据域 [10]，减少了重建误差；近年来一些基于延展分量函数的算法已经被提出来 [11] [12] [13] [14] [15]。在射电天文图像重建领域里，CLEAN算法是应用最为广泛的反卷积技术。CLEAN算法的思想决定了一个更合理的delta函数集合将能更好地逼近真实图像。在当前仅有的优化方法是将每个delta函数乘以一个相同的循环增益。这个固定的循环增益是使用者设定的，通常是一个经验值，大约在0.01~0.25之间 [1]。这种固定的经验循环增益很大程度上限制了图像的重建质量。本文引进基于幅度的方法来优化这个delta函数集合。

2. 射电综合成图的基本问题

设真实天空图像 $I^{true}$ (或者称为源)的可见度函数为 $V^{true}$，测量的可见度函数为 $V^{meas}$，采样模式为M，空间频率域的噪声为 $n_0$，那么测量的数据可以表示为

$V^{meas}=M\left(V^{ture}+n_0\right)$.

在数学上，真实的可见度函数实际上是真实天文图像的傅立叶变换

$V^{ture}=F\left(I^{true}\right)$,

采样模式M是由0和1组成的一个矩阵，0代表未采样和1代表采样。 $n_0$ 是加性噪声，常常被看成是高斯白噪声。它是空间频率的噪声的一个整体表达，这些噪声来自于接收系统，天空和地面等。噪声与源信号一起被采样模式所采样。卷积定理指出，空间频率域的两函数相乘等于在空域各自的傅里叶反变换卷积。由此可知，

$D=F^{-1}\left(V^{meas}\right)=B\ast I^{true}+n$

其中D为脏图，是测量的可见度函数的傅里叶反变换的结果； $F^{-1}$ 为傅立叶反变换；B为脏束，是采样模式的傅里叶反变换； $I^{true}$ 是需要求的未知函数—真实天空图像； $\ast$ 表示卷积运算；n是空间域的噪声 $n=F^{-1}M{n}_0=B\ast F^{-1}{n}_0$，显然空域噪声是被脏束所卷积，在像素之间，它不再是独立的。如图1，脏图等于真实天空图像和噪声与点扩展函数的卷积之和。

Figure 1. Image blurring model in space domain

图1. 空域的图像模糊模型

3. CLEAN反卷积算法的基本思想

为求得这个未知量 $I^{true}$，CLEAN算法本质上是在空间频率域最小化下面的目标函数，

${\chi }^2=\frac{1}{2}{\Vert V^{meas}-F\left(I^{appr}\ast B\right)\Vert }_{l_2}^2$

其中 $I^{appr}$ 是逼近的函数模型。CLEAN算法采用了最速下降法来逼近 [11]，也就是说每次更新的方向是残差中的最大绝对值点。CLEAN算法通过迭代的方式找到函数模型来逼近真实天空图像，

$I^{true}=I^{appr}+\epsilon$

其中， $\epsilon$ 为逼近函数模型与真实天空图像的误差。为优化反卷积过程，CLEAN算法使用循环增益来优化重建过程和结果，那么这个逼近的函数模型被表示成

$I^{appr}={\displaystyle \underset{i}{\sum }g{p}_i{\delta }_i}$

其中g是循环增益， $p_i$ 是第i次迭代的最大绝对幅度值， ${\delta }_i$ 在最大绝对幅度值位置的delta函数。

4. 基于幅度的CLEAN重建算法

正如上面提到，当前循环增益是一个算法使用者设置的经验值，而且在整个反卷积过程中是固定的。也就是对于所有的分量，使用相同的循环增益来优化。CLEAN算法是一种迭代算法，优先重建强信号。随着反卷积的深入，残差中的信号会越来越少。所以在不同的反卷积阶段，残差中所包含的信号比例是不一样的。使用相同的循环增益来优化重建过程显然不是最优的方式。根据上面的分析，循环增益应该和残差中所包含的信号比例成正比。除了考虑信号的幅度，由于测量是含噪声的，所以噪声因素必须考虑进来。我们将这些信号的内部特征考虑进入CLEAN重建过程中，循环增益 $g_0$ 表示为

$g_0=\frac{\left|p_{curr}-n_{level}\right|}{\left|p_{maxi}-n_{level}\right|}$

其中 $p_{curr}$ 代表当前残差的幅度最大绝对值， $p_{maxi}$ 代表脏图的幅度最大绝对值， $n_{level}$ 代表空域噪声水平，使用空域噪声的均方差来表示

$n_{level}=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{{\displaystyle \sum }}}{\left(x_i-\bar{x}\right)}^2}$

其中， $\bar{x}$ 为噪声的均值，N为图像的像素点数。在实际处理过程中， $n_{level}$ 是在非源区域测得的。考虑到延展源情况和重建速度，优化上面表达如下

$g=\{\begin{array}{l}k{g}_0,k{g}_0>0.05;\\ 0.05,k{g}_0\le \mathrm{0.05.}\end{array}$

其中k是一个考虑延展情况的参量，在实验中我们发现k为0.4~0.8时能取得很好的结果。当 $g_0\le 0.05$ 时，循环增益被置为0.05，主要是为重建速度考虑。这个值也参考了当前CLEAN算法的循环增益的建议值，大量经验表明当前的CLEAN算法的循环增益使用0.1及以下的值，性能会更优，故这里选取0.05。

5. 数值模拟

CLEAN算法的发展已经远远超越发明者的预期，在射电天文及其它应用领域获得了广泛的应用 [16]。CLEAN算法不但在二维图像重建中获得了成功，并且在一维信号恢复中也有广泛的应用，如用在光谱分析中的非等间隔采样的信号恢复 [17] [18]。当前Högbom CLEAN算法仍然被广泛使用，在诸如CASA等射电干涉数据处理软件中都有它的实现。为了验证本文改进算法的性能，我们比较了Högbom CLEAN算法和本文的改进算法。

在本文中，选择了对于这类CLEAN算法具有挑战性的源分布—射电源3C31，也就是延展源存在的情况，这也是实际中占大多数的情况。基于CASA的模拟是非常接近真实观测，从反卷积的角度，可以简单认为它与真实观测无异。

本文也从数值的角度来比较二者的性能。我们引入动态范围 [19] 来判断这两种算法的优劣。动态范围是一个常用且重要的天文图像的性能指标，图像重建追求一个高的动态范围。发展图像算法的一个的重要目标是增加重建图像的动态范围。动态范围被定义为重建图像的最大值与非源所在区域的残差水平之比。如果设动态范围为 $D_r$，重建图像为 $I^{resto}$，残差图像为 $I^{resi}$，非源所在区域的残差水平为 $e_{rms}$，残差中非源所在区域的像素数目为N。那么

$D_r=\frac{\max \left(I^{resto}\right)}{e_{rms}}$

其中非源所在区域的残差水平为

$e_{rms}=\sqrt{\frac{{{\displaystyle \sum }}_i{\left(I_{\left(i\right)}^{resi}\right)}^2}{N}}$。

通用天文软件程序包(CASA)的前身是AIPS，是为了应对ALMA和VLA对于数据处理的新挑战，并使用C++语言重新开发的软件。CASA已经成为处理射电干涉测量数据的主流软件。它自带有模拟观测功能，模拟观测的数据格式等都与实际观测的完全相同。可能的不同在于，模拟的数据不会出现实际中因某一望远镜临时出现问题而无法得到数据的问题。另外，模拟噪声可能与实际的不太一样，实际情况，可能因为天气或故障，引入的噪声会有差别。模拟的噪声常常使用的是高斯白噪声。然而，考虑到噪声的随机性，这个问题可以忽略。所以模拟数据相当于真实测量数据校准之后的数据。

Figure 2. From left to right: (a) The radio 3C31 image; (b) The simulated dirty image; (c) The restored image from the Högbom CLEAN algorithm; (d) The restored image from our algorithm

图2. 从左到右：(a) 射电3C31源；(b) 模拟观测的脏图；(c) Högbom CLEAN算法的重建图像；(d) 本文改进算法的重建图像

在这个实验中，我们使用双鱼射电星系3C31源(图2(a))作为参考图像，模拟了VLA B阵型在L波段观测了6小时。像素大小(cellsize)为1.0 arcsec (角秒)，图像尺寸为512 × 512个像素。在成图过程中使用Briggs权重函数对数据进行加权。观测后的数据变换成的脏图显示在图2(b)，显然脏图是模糊的。图2(c)和图2(d)分别显示了Högbom CLEAN算法和本文改进的算法的重建结果。对比于参考图像，改进的方法能够有效地重建出相应的天体图像。我们也测量了残差图像的均方根误差，如表1。改进的算法在均方根误差上也有更好的表现。这从统计上说明，本文改进的算法能够更好地提取信号。换句话说，更好地分离了源与噪声。同时我们计算了Högbom CLEAN算法和改进算法的重建图像的动态范围，后者的动态范围更高(见表2)。同时也测试了本文算法对于射电M31的重建效果。从表2可以看出，与前一个实验一致，本文改进方法的重建图像的动态范围要高于Högbom CLEAN算法。

Table 1. Comparison of RMS error for the results from the two algorithms

表1. 两种方法重建结果的均方根误差比较

Table 2. Comparison of dynamic range for the results from the two algorithms

表2. 两种方法重建结果的动态范围比较

6. 总结

CLEAN算法使用delta函数集合来逼近真实的天文图像。传统的方法使用经验的固定的循环增益来限制delta函数的幅度。本文将幅度和噪声水平之间的关系引进了CLEAN算法来优化这个delta函数集合，使得改进后的算法有更好的性能。结果显示改进的方法能够得到更高动态范围的结果。本文改进的算法是基于CASA和python编程实现的。

基金项目

国家自然科学基金(11963003)，SKA专项资助(2020SKA0110300)，国家重点研发计划(2018YFA0404602)，贵州省教育厅青年科技人才成长项目(黔教合KY字[2018]119)，贵州大学引进人才科研基金(贵大人基合字(2018)60号)，中国科学院太阳活动重点实验室基金(KLSA201805)，贵州省科技计划项目(黔科合平台人才[2017]5788号)。

NOTES

^*通讯作者。

¹http://www.vla.nrao.edu/。

²http://www.almaobservatory.org/。

³http://www.skatelescope.org/。