基于L-BFGS-B局部极小化的自适应尺度CLEAN算法

张利 肖一凡 米立功 卢梅 赵庆超 王蓓 刘祥 张明 谢泉

摘 要：干涉阵列存在的点扩展函数旁瓣使得观测到的射电源出现不同程度的失真，对重建宇宙真实结构图景和理解宇宙起源造成影响。为解决观测中出现的伪影，本文在现有的CLEAN算法的基础上，提出了基于L-BFGS-B局部极小化的自适应尺度CLEAN算法。首先，基于L-BFGS-B局部极小化算法通过最小化目标函数，寻找最优分量，构建自适应尺度模型;其次，通过CASA实现对测试图像的重建，对比目前广泛使用的Hgbom CLEAN算法的重建图像，评估本文算法性能;最后，对反卷积算法在射电天文图像处理领域的发展做出展望。测试结果表明：相比于传统的算法，本文提出的算法能够构建更加精准的天空亮度分布，为天文图像重建提供一种新的方案。

关键词：L-BFGS-B算法;自适应尺度CLEAN算法;射电天文图像处理

中图分类号：TP391.9

文献标志码：A

通过观测所接收到的电磁辐射来认识宇宙结构，理解宇宙起源及演化，是人类不断追求的目标。在射电波段对天体和宇宙开展研究的天文学分支称为射电天文学。近年来，射电天文学取得了显著突破，包括脉冲星和中子星[1]、星系中的暗物质[2]、由超大质量黑洞驱动的射电星系[3]和类星体[4]等，这对于理解宇宙、恢复宇宙图景具有重要意义。

研究宇宙演化的重要手段之一是研究射电天文图像及其结构变化。由于银河系外天体（例如类星体和星系）的射电辐射到达地球时信号已经非常微弱，只有借助具有较高灵敏度的大型射电干涉阵列才能检测到这些信號。目前已建成的大型干涉阵列包括巨米波射电望远镜（giant metrewave radio telescope， GMRT）、默奇森宽场阵列望远镜（murchison widefield array，MWA）、低频阵列射电望远镜（low frequency array，LOFAR），此外还有若干正在建设的新型干涉阵列，例如平方公里射电阵（square kilometre array，SKA）。

随着越来越多的巡天计划的实施和大型干涉阵列的建成，人们能够得到海量天文数据，但是望远镜数量有限和地球自转综合造成的稀疏采样使射电天文图像变得模糊。为了对源的结构特性做更精准的分析，需要对观测到的射电天文图像进行去模糊处理。Hgbom[5]开创性地提出了CLEAN反卷积算法，通过迭代识别点源并消除点扩展函数（PSF）的影响，有效解决由于稀疏采样引起的图像模糊问题。该算法的提出具有超越时代的意义，并延伸应用到其它领域。CLARK算法[6]和Cotton-Schwab算法[7]对其作出改进，在提高反卷积速度的同时，减少算法误差，对展源的处理更加精准。然而点源分解机制和点扩展函数的旁瓣导致包含延展源和复杂特征的重建图像中出现条纹（stripes）。研究人员提出了一些改进算法解决该问题，例如，尺度敏感的CLEAN算法[8-12]，利用尺度基函数将天文图像参数化，以表达像素之间的相关性，从而进一步消除图像残差，但对于重建图像质量仍有提升空间。本文提出一种自适应尺度CLEAN算法，通过局部极小化目标函数，寻找最优分量，以提高反卷积的性能，从而更加精准地重建射电天文图像。

1 理论基础

1.1 干涉成像原理

天空亮度分布可以分解为奇对称成分和偶对称成分，在射电天文望远镜进行动态范围成像时，使用一对cosine 和 sine 相关器的组合，即复相关器，来输出天空图像的可见度：

V=Rc-iRs=∫I（s^）exp（-2π·s^λ）dΩ。 （1）

其中，Rc和Rs分别表示“cosine”和“sine”相关器的输出，I（s^）表示表面亮度分布，b表示天线之间的基线矢量，s^表示辐射源的方向矢量，λ表示波长。

在（u，v，w）直角坐标系中（如图1），其中，基线矢量为b=（u，v，w），方向矢量s^=（l，m，1-l2-m2），辐射源的立体角dΩ=dldm/1-l2-m2，其中，l和m分别表示s^对u轴和v轴的投影长度，天空可见度可表示为：

V（u，v，w）=Iv （l，m）1-l2-m2exp[-2πi（ul+vm+w1-l2-m2）]dldm。（2）

在满足小视场成像和所有基线矢量共面的情况下，上式可消去w额外项，变为二维傅里叶变换，然后通过逆傅里叶变换，得到天空亮度分布。

受限于干涉阵列的天线数目和基线的长度范围，所观测的uv平面覆盖不完全，因此无法获得天空亮度分布的全部信息，傅里叶逆变换仅能得到脏图Idirty（l，m），

Idirty（l，m）1-l2-m2=V（u，v）S（u，v）exp[-2πi（ul+vm）]dldm。（3）

其中，S（u，v）表示采样函数。采样函数是干涉阵列的点扩展函数B（l，m）的傅里叶变换。根据卷积定理，上式可表示为：

Idirty（l，m）=I（l，m）*B（l，m）。 （4）

脏图存在明显的旁瓣等干扰，无法反映真实的天空亮度分布，因此利用反卷积算法，例如CLEAN算法、最大熵算法[14]，从脏图中尽可能地对真实射电天文图像进行重建。

1.2 Hgbom CLEAN算法

Hgbom提出的CLEAN算法将天空亮度分布假设为一系列点源集合，在图像域中通过迭代洁化的过程，消除脏图中的旁瓣效应。对Hgbom CLEAN（Hg-CLEAN）算法的具体实现过程如下：

步骤1 寻找脏图中最大绝对亮源，取其位置（xi，yi）和幅度ai作为初始参数;

步骤2 更新模型分量，Imodeli=Imodeli-1+g×aiδ（x-xi）（y-yi）。 （5）

其中，g表示循环增益;

步骤3 更新第i次残差图像，

Iiresidual=Idirty-B（l，m）*Iimodel。（6）

其中，*表示卷积运算;

步骤4 重复步骤1—3，直到迭代次数超过预设阈值或残差图像中的最大值小于规定的噪声水平时停止。将洁束卷积模型图像，再加上残差图像，得到恢复出来的图像。

这种算法对于点源集合的处理具有很好的效果，然而点源无法表示延展源中像素间的相关性，因此不能很好地重建展源。

2 基于L-BFGS-B局部极小化的自适应尺度CLEAN算法

受CLEAN算法[5，12]和文献[15]的启发，本文提出一种L-BFGS-B局部极小化的自适应尺度CLEAN算法，和其它CLEAN算法的变种一样，使用相同的框架[10]。本文引入的L-BFGS-B算法由L-BFGS算法[16]改进得到，是一种基于梯度投影的非线性优化方法。该方法结合原有算法的Hessian、线搜索算法和信任域方法应用于更新过程。

L-BFGS-B算法的迭代公式为：

xi+1=xi-αiHifi。（7）

其中，αi表示步长调节因子，Hi为Hessian矩阵，fi为多维一阶向量梯度。Hi通过公式在每次迭代中更新参量，

ρi=1yiTsi;（8）

Vi=I-ρiyisiT; （9）

si=xi+1-xi; （10）

yi=fk+1-fk。 （11）

在本文算法的每一次迭代中，L-BFGS-B算法按照以下步骤寻找最优解：

步骤1 结合信任域方法，根据下列公式在点xi附近用二次型模型对目标函数进行局部拟合，计算柯西点的近似值，

m（x）=f（xi ）+giT（x-xi ）+（x-xi）THi2。（12）

其中，（）T表示转置运算，Hi表示在第i次迭代的Hessian矩阵;

步骤2 通过直接法、共轭梯度法或对偶法计算搜索方向dk;

步骤3 依据问题的约束条件，沿搜索方向dk进行线性搜寻，计算步长调节因子αk，更新参数xk+1=xk+αkdk，以寻找函数的最小值;

步骤4 采用L-BFGS-B算法更新Hessian矩阵Hk，同时检查是否收敛。

重复步骤1—4，直到满足L-BFGS-B算法的3个停止准则之一：达到最大迭代次数、目标函数的减少量变小时，投影梯度的范数足够小时，循环终止。从而得到函数最小值，阈值的选择取决于脏图和点扩展函数的实际情况。

然后，利用一系列高斯函数对残差图像进行平滑处理，取结果的最大值作为初始参数，

ωi=ωic2-ωb2。（13）

其中，ωi表示第i次迭代中重新计算得到的高斯模型的尺寸，ωic表示第i次迭代中由最小化算法拟合残差得到的高斯模型尺寸，ωb是点扩展函数的高斯模型尺寸的近似值。优化后的幅值ai可由下式计算得到：

ai=aicωic22πabωb2ωi2。（14）

其中，aic和ab分别表示拟合残差得到的幅值和点扩展函数的模型幅值。接着，计算模型尺度的大小，并根据高斯尺度得到第i次拟合后的模型分量：

Iiguass=∑Ni=1aiexp-12（x-xi）2+（y-yi）2ωi2。 （15）

使用该最优分量对模型图像进行更新，同时利用循环增益g对模型分量进行优化，以建立更加精确的天空模型，

Iimodel=Imodeli-1+gIigauss。 （16）

其中，Iimodel表示第i次迭代的模型图像。接着利用模型图像計算残差图像，更新第i次残差图像，

Iiresidual=Idirty-B（l，m）*Iimodel。 （17）

其中，*表示卷积运算。对模型和残差图像进行循环更新，直到满足最大迭代次数或者达到噪声水平时停止，得到最后的模型图像Imodel和残差图像Iresidual，重建后的图像Irestored表示如下公式：

Irestored=Bclean*Imodel+Iresidual。（18）

其中，Bclean为拟合主瓣得到的干净光束。改进后的算法流程如图2所示。

本文采用L-BFGS-B算法对分量进行更准确的拟合。对比Hg-CLEAN算法，该算法通过对最小化目标函数部分的优化，能够对天空亮度分布做出更精准的建模，使重建后的图像更加接近真实天空图像。

3 试验结果分析与算法评估

本文使用天文通用软件包（common astronomy software applications，CASA）模拟JVLA B阵型观测仙女星系M31图像，并用于验证新算法的性能，图像尺寸为256×256像素。图3展示了本算法重建图像前后的对比效果。对比图3（a）中参考的干净天空图像，图3（b）显示了在点扩展函数旁瓣影响下的脏图，其中包含大量的旁瓣，造成图像中天体的细节结构模糊不清。图3（c）是本文算法实现的图像重建效果，从图中能够观察到，脏束的旁瓣效应基本被消除，天空图像中的点源和展源附近的模糊得到较好的处理，说明本文算法能够从脏图中有效重建天空图像。

相对于Hg-CLEAN算法，本文算法呈现出更好的图像重建性能。图4是本文算法与传统的Hg-CLEAN算法对于同一图像的处理结果，展示了对应的模型图像和残差图像。从图中可以看出，本文算法得到的残差图像中的结构远少于Hg-CLEAN算法，说明采用该算法得到的重建图像，能够恢复原始图像中的大部分信息，同时能够体现L-BFGS-B算法的拟合参数效果良好，提高了算法对天空亮度分布建模的精度;观测最终的重建结果，发现Hg-CLEAN算法也能够对原始图像的部分信息进行重建，然而在最终的结果中仍带有大量的伪影，并且在其残差图像中包含明显的图像结构。

另外，從图4（a2）中能够观察到明显的延展信号，图5展示了模型图像中展源部分的处理效果。与Hg-CLEAN算法相比，本文提出的算法在展源处理方面效果更好，原因在于Hg-CLEAN算法将天空图像分解为一系列delta函数，不能很好地表达图像中的展源结构，而本文算法能够依据射电天文图像中的展源结构自适应地选取尺度参数，因此具有处理图像中延展源的性能。

同时，采用客观指标对天空图像重建效果做量化分析，选取峰值信噪比（peak signal-to-noise ratio，PSNR）[17]、结构相似性（structural similarity，SSIM）[18]和均方根误差（root mean square error，RMSE）[19]作为重建图像的质量评价指标，对Hg-CLEAN算法和本文提出的算法进行评价。为了对比不同算法之间的差异，在同样运行环境下对不同算法进行测试，从表1可以看出，本文提出算法在PSNR、SSIM和RMSE指标上均有所提升。

为了验证算法的稳定性，采用三种加权方案对原始图像做处理，即自然加权， 均匀加权和Briggs 加权。图6展示了不同加权方案下图像的结构变化，观察可以发现经自然加权后的图像中的旁瓣效应最为明显，均匀加权后的图像的旁瓣效应最弱。

图7显示了对于不同加权方案的重建结果。对比残差图像，本算法对不同加权下的脏图均有较好的重建效果，尤其对于均匀加权和Briggs加权具有明显优势。

4 结语

CLEAN算法对于射电天文图像处理领域具有重大意义，并能够延伸到天文学相关的其他领域[20]，这系列算法对传统的图像处理算法造成冲击，对图像处理的相关领域也产生了深远影响。本文基于目前已有的CLEAN算法[12]，结合L-BFGS-B算法对目标函数的最小值求解进行优化，并将算法的重建结果与传统的Hg-CLEAN算法[5]做比较，提高了原有算法的性能，使重建得到的射电天文图像更加接近真实天空图像，为测试该算法对不同观测图像的普适性，在以后的工作中，计划扩展到不同阵列和不同科学目标的观测成像中去。

参考文献：

[1] MUSHTUKOV A A，SULEIMANOV V F，TSYGANKOV S S，et al. On the maximum accretion luminosity for magnetized neutron stars and the nature of ultraluminous X-ray sources[J]. Monthly Notices of the Royal Astronomical Society， 2018， 454（3）：2539-2548.

[2] BEHROOZI P S，WECHSLER R H，CONROY C.The average star formation histories of galaxies in dark matter halos from z=0～8[J]. Astrophysical Journal， 2013， 770（1）：2394-2404.

[3] KAUFFMANN G，HECKMAN T M，BEST P N.Radio jets in galaxies with actively accreting black holes： new insights from the SDSS[J]. Monthly Notices of the Royal Astronomical Society， 2010， 384（3）：953-971.

[4] POUNDS K A，KING A R，PAGE K L，et al. Evidence of a high velocity ionised outflow in a second narrow line quasar PG0844+349[J]. Monthly Notices of the Royal Astronomical Society， 2018， 346（4）：1025-1030.

[5] HOGBOM J A.Aperture synthesis with a non-regular distribution of interferometer baselines[J]. Astronony Astrophys Suppl， 1974， 15：417-426.

[6] CLARK B G.An efficient implementation of the algorithm 'CLEAN'[J]. Astronomy & Astrophysics， 1980， 89（3）：377-378.

[7] SCHWAB F R.Relaxing the isoplanatism assumption in self-calibration：applications to low-frequency radio interferometry[J]. Astronomical Journal， 1984， 89（7）：1076-1081.

[8] CORNWELL T J.Multiscale CLEAN deconvolution of radio synthesis images[J]. IEEE Journal of Selected Topics in Signal Processing， 2008， 2（5）：793-801.

[9] ZHANG L.Fused CLEAN deconvolution for compact and diffuse emission[J]. Astronomy and Astrophysics， 2018， A117：618-624.

[10]RAU U，CORNWELL T J.A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry[J]. Astronomy & Astrophysics， 2012， A71：532-549.

[11]ZHANG L，MI L G，ZHANG M，et al. Adaptive-scale wide-field reconstruction for radio synthesis imaging[J]. Astronomy and Astrophysics， 2020，A80：1-10.

[12]ZHANG L，BHATNAGAR S，RAU U，et al. Efficient implementation of the adaptive scale pixel decomposition algorithm[J]. Astronomy and Astrophysics， 2016， A128：592-600.

[13]SMITH F G.Book Review：interferometry and synthesis in radio astronomy by a.richard thompson，james m.moran， george w. swenson jr[J]. IEE Proceedings F-Communications， Radar and Signal Processing， 2017， 133（7）：648-650.

[14]NARAYAN R，NITYANANDA R. Maximum entropy image restoration in astronomy[J]. Annual Review of Astronomy and Astrophysics， 2003， 24（1）：127-170.

[15]ZHU C，BYRD R H，LU P，et al. Algorithm 778： L-BFGS-B：fortran subroutines for large-scale bound-constrained optimization[J]. Acm Transactions on Mathematical Software， 1997， 23（4）：550-560.

[16]KOKKOLARAS M C，AUDET W.Hare：derivative-free and blackbox optimization. Springer series in operations research and financial engineering[J]. Optimization and Engineering， 2019， 20（3）：955-957.

[17]PAPPAS T N，SAFRANEK R J，CHEN J.Perceptual criteria for image quality evaluation[J]. Handbook of Image and Video Processing （Second Edition）， 2005：939-959.

[18]WANG Z，BOVIK A C.A universal image quality index[J]. IEEE Signal Processing Letters， 2002， 9（3）：81-84.

[19]LAMBRECHT C V.Special issue on image and video quality metrics[J]. Signal Processing， 1998， 3（3）：153-154.

[20]張利， 徐龙， 米立功，等. 射电天文图像的反卷积算法研究[J]. 天文学报， 2018， 59（6）：117-124.

（责任编辑：于慧梅）

An Adaptive Scale CLEAN Algorithm Based on

L-BFGS-B Local Minimization

ZHANG Li*1， XIAO Yifan1， MI Ligong1， LU Mei1， ZHAO Qingchao1，

WANG Bei1， LIU Xiang2， ZHANG Ming2， XIE Quan1

（1.College of Big Data and Information Engineering， Guizhou University， Guiyang 550025， China; 2.Xinjiang Astronomical Observatory， Chinese Academy of Sciences， Urumqi 830011， China）

Abstract：

The sidelobes of point-spread function in interferometer arrays make the observed radio sources have varying degrees of distortion， which affects the reconstruction of the real structure of the universe and the understanding of the origin of the universe. In order to solve the artifacts in observation， on the basis of the existing CLEAN algorithm， this paper proposes an adaptive scale CLEAN algorithm based on L-BFGS-B local minimization. Firstly， based on the L-BFGS-B local minimization algorithm， an adaptive scale model is constructed by minimizing the objective function and finding the optimal component; secondly， to evaluate the performance of the proposed algorithm; the test image is reconstructed through CASA， which is compared with the widely used Hgbom CLEAN algorithm; finally， the development of deconvolution algorithm in radio astronomy image processing is prospected. The results show that compared with the traditional algorithm， the algorithm proposed can build a more accurate sky brightness distribution， which provide a new scheme for the reconstruction of astronomical images.

Key words：

L-BFGS-B algorithm; adaptive scale CLEAN algorithm; radio astronomical image processing