射电天文图像重建算法研究综述

A Review of Radio Astronomical Image Reconstruction Algorithms

摘要: 由于射电望远镜采样不完备，导致观测到的天文图像质量下降，影响天文与天体物理研究。为了得到质量更好的天体图像，图像重建算法成为射电天文领域具有研究价值的方向。本文综述了射电天文成像重建领域应用广泛的算法。对Högbom CLEAN算法、Cornwell提出的多尺度CLEAN算法和最大熵算法进行具体的分析研究，并分别将其应用于射电成像的重建；通过对比这些算法在成像中的实现结果，提出领域内亟待解决的问题与挑战。

Abstract: Due to the incomplete sampling of radio telescope, the quality of the observed astronomical image is reduced, which seriously affects the scientific research work. In order to obtain better celestial images, image reconstruction algorithms have become a valuable research field in radio astronomy. In this paper, we research the radio field of astronomical imaging reconstruction algorithms used widely, and then take Högbom CLEAN algorithm, the multi-scale CLEAN algorithm and the maximum entropy algorithm detailed analysis and the imaging reconstruction was applied to the radio, finally compare the algorithm implemented in the imaging results and put forward challenges and problems to be solved within the territory.

1. 引言

目前人类对宇宙的大多认识都来源于望远镜所接收的电磁辐射，射电天文观测已成为人类探索宇宙、研究天体的重要方式。Karl G.Jansky于1932年首次发现来源于银河系中心的射电辐射 [1]，5年后，Grote Reber在160 MHz探测到了银河系的射电辐射并进一步开展了银河系的第一次射电巡天观测 [2]，自此在射电波段对天体和宇宙开展研究的射电天文学开始受到广泛关注。射电辐射能够反映天体的特性和状态，根据探测到的辐射频率可以对源区的磁场和电子密度进行判别 [3]。目前研究射电辐射采用由射电望远镜接收后转换成射电天文图像的方法，至今已取得众多新发现，比如脉冲星和中子星 [4]、宇宙大爆炸的CMB [5] 等。

由Martin Ryle和Antony Hewish发明的综合孔径(aperture synthesis)技术 [6] 对射电干涉测量进行改革。由Van Cittert-Zernike定理和傅里叶变换理论可知，通过测量辐射源的辐射分量，就能够恢复真实源的亮度分布，重建真实天空图像 [7]。在满足所有基线共面和小视场成像的前提下，真实天空的亮度分布 $I_v$ 可表示为：

$\frac{I_v\left(l,m\right)}{\sqrt{1-l^2-m^2}}={\displaystyle \iint V\left(u,v\right)\exp \left[2\pi i\left(ul+vm\right)\right]\text{d}u\text{d}v}$ (1)

其中， $V\left(u,v\right)$ 表示采集到的可见度数据， $l,m$ 分别为方向矢量s对 $u,v$ 轴的投影长度。

然而由于望远镜的天线个数和基线长度都是有限的，不可能得到全部的辐射分量，造成UV覆盖面不完全，并且观测不可避免引入的噪声，这些因素都将造成成像模糊。真实情况下得到的天图被称为脏图 $I^{dirty}$，可表示为：

$\frac{I^{dirty}\left(l,m\right)}{\sqrt{1-l^2-m^2}}={\displaystyle \iint V\left(u,v\right)S\left(u,v\right)\exp \left[-2\pi i\left(ul+vm\right)\right]\text{d}l\text{d}m}$ (2)

其中 $S\left(u,v\right)$ 是干涉阵列的点扩展函数 $B\left(l,m\right)$ 的傅里叶变换。根据卷积定理，上式可表示为：

$I^{dirty}\left(l,m\right)=I_v\left(l,m\right)\ast B\left(l,m\right)+I^{noise}$ (3)

其中 $I^{noise}$ 表示真实观测条件下引入的噪声。脏图包含旁瓣等干扰结构，细节信息无法被表达，不能直接用于天文或天体物理研究，因此利用天文图像重建算法从脏图中恢复真实天空亮度分布对射电天文研究具有重要意义。

首次提出的CLEAN算法(Hg-CLEAN算法) [8] 利用点扩展函数通过简单的循环迭代以识别和去除脏图中的点源伪影。Schwarz在1978年对该算法进行了细致分析并验证了收敛条件 [9]。对于尺寸较小的图像，Hg-CLEAN算法的计算速度上具有优势，而在处理尺寸较大的图像时，运行速度明显下降。为了克服这一缺点，Clark提出了CLEAN算法的一个变种，只取脏束中心部分进行CLEAN迭代，这种方法在保证重建质量的基础上将运行速度提高了一个数量级 [10]。邱耀辉等人结合Hg-CLEAN算法与迭代位移叠加法提出一种在空域重建天文图像的算法，得到更好的重建图像质量 [11]。

上述方法虽然在去除点源附近伪影方面取得了较好的效果，但存在以下缺点：1) 对于展源附近的伪影不能很好的处理，对于包含较多展源结构的图像，重建出的图像仍存在较为严重的失真；2) 算法成本较高，所需存储空间大且运行时间长。Bhatnagar和Cornwell提出尺度敏感的反卷积算法(Asp-CLEAN算法) [12] 能够对非对称结构的天空图像进行更精准的建模，但其算法复杂程度上升，使运行时间延长为Hg-CLEAN算法的三倍。Cornwell提出的多尺度CLEAN算法(MS-CLEAN算法) [13] 对于展源的处理效果有明显提升，然而仍需要比Hg-CLEAN算法更长的运行时间。此后Zhang等人提出具有快速收敛特性的自适应CLEAN算法(Asp-CLEAN算法) [14]，相对于尺度敏感反卷积算法的运行速度提升了50%。

最大熵算法是另一类射电天文图像重建方法。最大熵算法的基本原理是在解空间集中找到能使信息熵H最大的解 [15]，

$H\left(p_1,\cdots ,p_n\right)=-k{\displaystyle {\sum }_{i=1}^n{p}_i\log p_i}$ (4)

其中，p代表不同事件的发生概率。

Frieden等人为研究天文图像，首次将最大熵方法引入图像重建领域，并提出图像熵的概念 [16]。该算法的本质是对非线性方程的求解，但由于计算量过大以及收敛性速度慢的问题一直未被解决，这种方法未能实际应用。此后，Skilling和Bryan提出的剑桥算法(Skilling-Bryan算法) [17] 为射电天文图像重建提供了新方案，其利用以多个搜索方向为基的子空间以寻找最优分量，极大地优化了算法效率。1985年Cornwell和Evans对最大熵算法进行优化，并广泛应用于射电天文图像重建领域 [18]。Starck等人于1996年提出了多尺度最大熵算法 [19]，借助小波基变换对图像进行多尺度分解，并分别在不同尺度上采取规定准则区分噪声结构和真实图像结构，再分别进行归一化处理，进而求解多尺度最大熵。通过大量实验证实，该算法能够在不同尺度上处理优化，从而精确重建图像，但缺乏严密的数学证明，尚有完善空间。Bonavito等人提出直接迭代法 [20]，并应用于重建哈勃望远镜拍摄的天文图像。该算法对旁瓣引起的模糊和过饱和区的处理具有明显效果，但由于未能考虑到噪声特性，将其与真实图像做同样处理，随着迭代次数增加，噪声明显被放大，因此迭代次数对重建图像的质量产生严重影响。Lyon等人对此方法进行改进，提出基于最大似然先验约束的最大熵算法 [21]，对于噪声信号做特征处理，得到质量更好的处理图像。

随着干涉测量技术的发展，具备更高灵敏度的望远镜被建造出来，这对射电天文图像的重建算法提出新的挑战。为了解决高频波段的成像问题，Rau结合Multiscale CLEAN算法与Multi-frequency算法在更高的采样频率下重建天空图像 [22]。为了解决视场较大时存在的边界效应问题，Camera对图像域进行划分，分别校正每个小区域的边界效应，该方法对于恒星系统的模拟观测图像展现出很好的去模糊效果 [23]。

本文综述了有关射电天文领域的图像重建算法。在第2节介绍了经典的射电天文图像重建的基本算法框架，即Högbom CLEAN算法、多尺度CLEAN算法和最大熵算法，并对这些方法进行比较和分析，在第3节利用天文观测数据对这些算法的性能进行评估和对比。在第4节讨论了射电天文图像重建领域面临的挑战和潜在的研究方向。

2. 射电天文图像重建算法

2.1. Högbom CLEAN算法

Högbom CLEAN算法是目前应用最广泛的反卷积方法，其利用先验知识对uv平面未被覆盖的区域进行插值，能够有效减少由于观测数据不完整造成的影响。通过一系列点源集合建立真实天空亮度分布的模型 [8]。该算法在图像域中进行多次迭代，在每次迭代中通过搜索峰值残差及位置信息以移除脏束，消除脏图中的旁瓣效应。根据以上分析，该算法框架如图1所示，其具体实现过程如下：

Figure 1. CLEAN algorithm framework

图1. CLEAN算法框架

步骤1搜寻脏图 $I^d$ 中绝对强度最大的点，标记其幅度 $\left(x_1,y_1\right)$ 和位置 $a_1$。

步骤2在残差图 $I_i^r$ (第一次迭代为 $I^d$ )中找到最大绝对亮度点源的位置 $\left(x_i,y_i\right)$ 和幅值 $a_i$，第i次分量可以定义为：

$I_i^g={\displaystyle {\sum }_1^i{a}_i\delta \left(x-x_i\right)\left(y-y_i\right)}.$ (5)

步骤3根据公式更新模型分量，在模型图像 $I_i^c$ 中添加高斯分量 $I_i^g$ 与循环增益 $\sigma$ 的乘积，

$I_i^m=I_{i-1}^m+\sigma I_i^g$ (6)

循环增益 $\sigma$ 取值范围为0~1，对于尺度敏感的CLEAN算法，通常取值小于0.5。

步骤4计算并更新得到第i次残差图像，

$I_i^r=I^d-B_D\ast I_i^m.$ (7)

根据步骤3~4，循环迭代更新模型和残差图像，直至达到迭代次数阈值n或残差图像中峰值点的幅值达到设定阈值s时，得到迭代结束后的模型图像 $I^m$ 和残差图像 $I^r$。

步骤5利用椭圆高斯模型拟合脏束的主瓣得到洁束 $B_C$，重建后的图像可以表示为：

$I^c=B_C\ast I^m+I^r.$ (8)

算法伪代码如下所示

研究证明算法具有收敛性，文献 [8] 给出了Högbom CLEAN算法的详细分析。这种简单执行的洁化算法有利于图像中点源的重建，然而由于在搜索残差以及移动脏束中消耗大量运行时间，对于较大尺寸的图像计算成本非常高；并且对于图像中展源结构无法进行更加精确的拟合，不能保证重建出的所有天文图像都具备比较好的质量。

2.2. Multi-Scale算法

传统的多尺度CLEAN算法由Cornwell在2008年首次提出 [13]，该算法作为Högbom CLEAN算法的一个变种，在算法的收敛性和稳定性方面有明显的提升；它将真实天空亮度分解为点结构和不同尺寸的扩展源，并据此进行建模。该算法预设尺度，然后针对每个尺度生成一系列残差图像，在每一个图像中进行迭代反卷积。与Högbom CLEAN算法相比，多尺度CLEAN算法根据预设值优先检测最大尺度的亮源，随着迭代依次寻找更加细小的结构，改进了原算法中对忽略细节的缺陷。根据以上分析，其具体实现过程如下：

步骤1将尺度 $m\left({\alpha }_j\right)$ 与残差图像 $I_j^R$ 进行卷积，

$I_{\alpha }^R=m\left({\alpha }_j\right)\ast I_j^R\left(n\right)$ (9)

常常计算尺度偏差 $L\left(\alpha \right)$ 以提高拟合精度，

$L\left(\alpha \right)=1-0.6\alpha /{\alpha }_{\max }.$ (10)

步骤2根据各个尺度 $m\left(\alpha \right)$ 上卷积后的残差图 $I_j^R$ (在第一次循环中为脏图)，寻找绝对强度最大的点，标记全局峰值 ${\beta }_j$ 、位置 $\left(x_j,y_j\right)$ 及对应的尺度 $m\left({\alpha }_j\right)$。将脏束中心位置移动至峰值所在位置，定义模型分量为：

$I_j^G={\displaystyle \sum {\beta }_{j}m\left(x-x_j,y-y_j,{\alpha }_j\right)}.$ (11)

步骤3根据下式更新模型分量，

$I_j^M\left(n\right)=I_{j-1}^M\left(n\right)+\sigma m\left({\alpha }_j\right)I_j^G$ (12)

在各个尺度上残差图像(在第一次循环中为脏图)中移除脏束与增益因子的乘积值，从而更新第n个残差图像，

$I_j^R\left(n\right)=I_{j-1}^R\left(n\right)-B_D\ast I_j^M\left(n\right).$ (13)

在所有残差图中重复步骤2~4，直到迭代次数超过预设阈值N，或者残差图中最大值低于预设噪声阈值S。将最后产生的模型图像和残差图像分别记为 $I^M\left(n\right)$ 和 $I^R\left(n\right)$。

步骤4分别在不同尺度上拟合脏束主瓣得到洁束 $B_C$，将其与函数阵中的值进行卷积，以得到不同尺度的重建图像分量：

$I^C\left(n\right)=B_C\ast I^M\left(n\right).$ (14)

步骤5将分量 $I^C\left(n\right)$ 和残差图像叠加，得到最终的重建图像 $I^C$，

$I^C={\displaystyle \sum I^C\left(n\right)+I^R\left(n\right)}.$ (15)

算法伪代码如下所示

多尺度CLEAN算法最大的优势在于选取不同尺度的分量，能够对展源结构进行更精确的拟合。然而经研究 [9] 发现对于所选择尺度的数量和分布，会对搜索速度和重建效果具有明显影响。因此在设定该值使需要考虑天文图像中天体结构的差异，以达到更适合的尺度，同时可以利用尺度偏差达到更好的拟合值。

2.3. 最大熵算法

最大熵算法(MEM算法)是一种非线性的图像重建算法 [18]。该算法在处理旁瓣信息和重建真实图像结构之间取得较好的平衡，并且不需要对图像先验知识做出详细的假设，能够获得优于线性图像重建算法的结果，但其数值求解面临计算量较大等问题。算法具体步骤如下所示：

步骤1将脏束 $B_D$ 在频域内与真实图像卷积得到脏图 $I_D$ 以及对应的傅里叶系数 $D_{uv}$，初始化与脏图等同大小的常数矩阵I，并得到其对应的傅里叶系数 $I_{uv}$。

步骤2计算当前脏图 $I_D$ 的傅里叶系数 $D_{uv}$ 和常数矩阵I的傅里叶系数 $I_{uv}$ 之间的统计量C。

步骤3使用参数 $\alpha$ 更换当前傅里叶系数 $I_{uv}$，得到的新的常数矩阵 $I_{xy}$，并计算其导数 $J_{xy}$ 及导数对应的系数 $J_{xy}$。

步骤4使用参数 $\beta$ 更换当前傅里叶系数 $J_{uv}$，得到的新的常数矩阵 $J_{xy}$，并计算其对应的熵 $I_{xy}$ 及熵的傅里叶系数 $J_{uv}$。

步骤5利用当前的 $J_{uv}$ 和 $I_{uv}$ 可计算得到重建图像的傅里叶系数 $I_{uv}$。

循环进行步骤2~5，直到达到迭代次数阈值N或者统计量C小于规定阈值时停止循环，迭代次数N和统计量C的阈值选取对重建图像质量与运行速度造成影响，因此需在两者之间做出权衡，以设置合理的阈值。

3. 实验与分析

本节将Högbom CLEAN算法、多尺度CLEAN算法和最大熵算法在射电天文成像中的实现结果进行对比。实验利用通用天文软件包CASA观测得到G055.7+3.4的Karl G. Jansky超新星残骸阵列的数据，并通过上述算法对图像重建。图2(a)显示观测得到的脏图，可以看出在点源和展源附近出现大量的伪影。

脏图(a) 点扩展函数(b)

Figure 2. G55 observation image

图2. G55观测图像

Högbom CLEAN算法重建图像如图3(a0)所示，可见反卷积算法对脏图质量具有明显的提高，图中的点源信息已基本恢复，但是很明显没有对图像进行足够全面的清理，图像中展源部分仍然有许多旁瓣，对图3(a1)残差图像的观察也表明重建图像中仍然含有噪声的源通量和结构。这是由于Högbom CLEAN算法假定图像中的结构由一系列点构成，不能对展源做更准确的拟合。如果需要对展源信息进行更准确的复原则要进行大量反复迭代处理，这将极大减慢算法的运行速度，因此Högbom CLEAN算法更适用于处理包含少量或者没有展源构成的射电天文图像。图3(b0)展示了多尺度CLEAN算法的处理效果，可以看出其对点源有较好的复原，能够重建部分展源信息，其恢复效果优于Högbom CLEAN算法，原因在于该算法设定由小到大的多个不同尺度，可以更加精确完整地处理点源和展源信息，对于脏图中不同尺寸的旁瓣都有很好的处理效果，并保证了较快的算法效率。最大熵算法的处理结果如图3(c0)和图3(c1)，可以看出残差图中展源附近的伪影得到很好的处理，但对于一些点源部分，这种算法的表现不如CLEAN算法。

Hg-CLEAN算法的重建图像(a0) Hg-CLEAN算法的残差图像(b0) MS-CLEAN算法的重建图像(c0)MS-CLEAN算法的残差图像(a1) MEM算法的重建图像(b1) MEM算法的残差图像(c1)

Figure 3. Comparison of reconstruction effects of different algorithms

图3. 不同算法的重建效果对比

点源及展源部分的局部特征在图4及图5中可以被更细致地观察。在展源局部结构的重建中，最大熵算法对伪影的处理明显优于CLEAN算法，多尺度CLEAN算法对比Högbom CLEAN算法也有较大改进，提升了图像质量。而CLEAN算法对点源处理表现更好，因此针对不同结构的观测图像，需要选取不同的重建方法。

Hg-CLEAN算法的处理结果(a) MS-CLEAN算法的处理结果(b) MEM算法的处理结果(c)

Figure 4. Zoom images of development source

图4. 展源局部特征图

在场边缘部分的点源细节图中可以发现CLEAN算法重建后的图像中，距离相位中心较远的一些较亮点源附近仍存在明显的弧线和斑点。究其原因，可能是由于天空曲率和非共面基线效应的影响，使在较低频段内成像的动态范围被限制，导致远离相位中心的源周围的伪影。这将在今后的工作中做进一步研究。

Hg-CLEAN算法的处理结果(a) MS-CLEAN算法的处理结果(b) MEM算法的处理结果(c)

Figure 5. Zoom images of point source

图5. 点源局部特征图

利用图像质量评价指标峰值信噪比(PSNR)、均方根值(RMSE)和结构相似性(SSIM)，对算法重建图像与原始图像进行评估，见表1。

Table 1. Reconstructed image quality index comparison

表1. 重建图像质量指标比较

4. 总结与展望

在过去的几年中，随着干涉测量技术的发展，射电图像重建算法也出现许多创新性进展。本文通过介绍综合孔径望远镜的成像原理，分析观测到的射电天文图像模糊的来源，即由于干涉仪在UV平面的采样点有限，通过这些不完备的观测数据无法恢复真实的天文图像 [24]。本文综述了射电天文成像重建领域两类应用广泛的算法：CLEAN反卷积算法和最大熵算法，然后详细调研了Högbom CLEAN算法 [8]、Cornwell提出的多尺度CLEAN算法 [13] 和最大熵算法 [18] 进行具体的分析研究，并分别将其应用于射电成像的重建。CLEAN算法对点源的处理效果明显，最大熵算法在展源部分的处理上更有优势，本文仅适用模拟观测的处理对其效果做对比，在实际观测数据中噪声更高，这些算法效率可能较低，需要在检测误差收敛和迭代次数中做出折中。如何权衡重建图像的质量和运行速度，如何在较高的噪声水平中重建真实图像，这些挑战是当前亟待解决的问题。

致谢

感谢天文通用软件CASA设计开发人员。

基金项目

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