一种基于多尺度带通滤波的洁化算法与GPU实现*

梅　盈，邓　辉，张晓丽，王　锋1,,4，陈腾达

(1. 中国科学院云南天文台，云南 昆明　650011；2. 昆明理工大学，云南 昆明　650500； 3. 广州大学天体物理中心，广东 广州　510006；4. 中国科学院大学，北京　100049)

射电观测是研究太阳剧烈活动的重要探测手段。建于中国明安图的螺旋天线阵列——明安图射电频谱日像仪(Mingantu Spectral Radioheliograph, MUSER)在400 MHz～15 GHz的584个频点上对太阳进行高分辨率成像，其观测将填补日本野边山日像仪(2个频点，17 GHz, 34 GHz)、法国南茜日像仪(5个频点，150～450 MHz)、俄罗斯伊尔库茨克射电日像仪(5.7 GHz)等在厘米、分米波段的观测空白，对太阳活动及其对人类影响的研究起到重要作用[1]。

明安图射电频谱日像仪分为低频阵(MUSER-I, 400 MHz～2 GHz)和高频阵(MUSER-II, 2 GHz～15 GHz)，其高时间、高空间和高频率分辨率给数据处理系统提出了较大的挑战。高频阵和低频阵的数据接收机每3 ms分别接收一次包含16通道的观测数据。在后续的数据处理过程中，如何高效地获取高质量的太阳图像成为数据处理系统的关键问题。

射电干涉阵成图的主要过程包含加权(Weighting)、网格化(Gridding)、傅里叶变换、退卷积等主要步骤，其中退卷积是对脏图进行洁化的过程，是成像过程中最耗时的部分。对日像仪成像来说，退卷积过程也是整个管线设计中最令人关注的问题。本文在细致调研相关工作的基础上，进一步对射电干涉阵成像中的多尺度退卷积方法进行了讨论，给出了相应的实现方法。

1　相关研究工作

射电干涉成像的基本原理是天空实际亮度是可见度函数的傅里叶变换。由于在实际观测过程中UV覆盖的不完全性，采样函数是离散的。自1974年Högbom提出洁化算法的思想后，洁化算法得到了极大的发展，并在射电望远镜的图像重建算法中得到广泛的应用。经典Högbom CLEAN算法的关键是为天空实际亮度提供一个模型(δ函数)，迭代寻找亮度的峰值位置，并如(1)式进行叠加，其中参数Ip为迭代过程中的峰值亮度，(xp,yp)为峰值亮度的位置，由残留图像的峰值亮度确定[2-3]。(2)式中的ID为脏图亮度，B为脏束，IC为(1)式计算出的叠加后的亮度值，IR为残图的亮度，n为迭代次数。

(1)

IR(n)=ID-B*IC(n-1).

(2)

Högbom CLEAN算法的一些变种算法，如Clark, Cotton-Schwab算法、最大熵算法(the Maximum Entropy Method, MEM)随后被提出。Steer Clean算法被用于日本野边山日像仪的成像，其洁化效率是Högbom算法的近10倍[4]。然而，基本的Högbom CLEAN算法由于在处理点源和展源上表现的优势，一直被广泛使用。

但随着新一代射电干涉仪的不断发展，在确保基本洁化算法的简单性和信噪比的前提下，多尺度成为提高效率的方法之一。Multi-Resolution CLEAN, Multi-Scale Maximum Entropy, Wavelets CLEAN等多尺度算法有效提高了洁化效率，但存在尺度大小固定和计算代价高等问题。

文[3]于2008年提出的多尺度洁化(Multi-Scale CLEAN)算法被认为算法直接且具有稳定性和收敛性好等特点。多尺度洁化为减少算法的迭代次数，提高算法的稳定性，采用多个不同尺度大小的卷积核与脏图进行卷积，得到多个卷积图像，在多个图像中寻找全局亮度峰值。算法的基本流程如下：

(1)选取不同的尺度大小，与脏图进行卷积，得到卷积后的系列图像；

(2)寻找全局亮度最大点，记录亮度最大点的位置和所在的尺度大小；

(3)计算步骤(2)中的尺度与脏束的卷积，再乘以增益值，存储计算结果及尺度大小；

(4)从(1)中所有的卷积图像中减去(3)的计算结果；

(5)重复以上4个步骤直至到达循环阈值或达到预设循环次数；

(6)将(3)中存储的结果与洁束卷积并加上参与图得到最终的洁图。

为了控制算法的迭代次数，提高图片质量，多尺度洁化算法需要指定不同尺度的卷积图像不同的权值，每次迭代循环结束后，残图需要乘以所在尺度的权值来增加或减少亮度幅值。经过系列实验，(3)式中权值α=0.6被认为是较合适的参数值[3]。

(3)

其中，S(a)为权值；a为当前尺度大小；a(maxscale)为最大尺度大小；α为权重参数。

为了减小迭代算法对增益的依赖，文[3]选择抛物线函数作为多尺度洁化算法的卷积核函数如下：

m(r,α)=ψ(r)[1-(r/a)2]，

(4)

由于高斯卷积核与(4)式给出的卷积核函数仅在较高动态范围下存在较小的差异，一般来说，采用不同尺度的高斯卷积核。以洁束大小为基准，卷积核大小选取常用的等比数列，即依次取洁束大小的1、2、4、8倍。

为进一步提高处理速度，近几年利用图形处理器技术加快退卷积处理成为主要的方法，如基于图形处理器的网格化算法使得w-projection的效率较中央处理器提高近百倍*https://arxiv.org/abs/1403.4209。前期基于图形处理器的Högbom CLEAN算法较中央处理器下的Högbom效率提高了十倍以上。但从多尺度洁化算法来看，其并行实现存在较大的难度。因为搜索最大值过程虽然可以并行，但需要利用找出的最大值要同时在所有尺度的图像中退卷积，这意味着一个尺度的数据处理必须与其它尺度的数据处理相交叉，影响了并行的实现。

2　基于频域的多尺度洁化

2.1　基本思想

参考文[3]的洁化方法，要加快洁化的效率，关键是能够同时进行洁化。实际上，由于Högbom Clean的基本原理，洁化迭代的过程始终是串行的，要并行的唯一可能性，是同时能够独立地进行多个洁化迭代。

法国默东天文台太阳干涉阵与明安图射电频谱日像仪一样，也是一个太阳专用的射电干涉阵，其数据处理采用了基于频率带通滤波的多尺度洁化方法[5]。该方法基于两个简化：(1)没有过多的展源存在；(2)空间尺度中只有相对较窄的谱。对于太阳观测来说，完全可以满足这两个简化需要。该方法与多尺度洁化方法有一定的相似性，但更易于并行实现。

参考法国默东天文台的洁化思想，结合明安图射电频谱日像仪的天线分布和最长基线情况，最终在本文实现的算法基本原理如下：

(1)在频域UV平面上，构建一系列连续的带通滤波器Fk(u,v)，k为滤波器编号，滤波器的半径以2的倍数增长。不同的滤波器对应不同的空间尺度。与法国默东太阳干涉阵不同，在实现时每次将计算最长基线和相应的UV分布，反过来计算相应的滤波器带通频率。

(2)通过滤波器，将原来的稀疏UV分布D(u,v)进一步拆分为若干个不同尺度的独立分布，即：Dk(u,v)=Fk(u,v)×D(u,v)；

(3)对于一个不同尺度的UV分布，通过逆傅里叶变换得到脏图：Ik=FT[Dk(u,v)]，然后进行常规的洁化处理；

(4)最后的洁图等于各尺度洁图的累加。

2.2　尺度的考虑

根据5个不同尺度的大小，将脏图分成不同的频率段，随着k值增大，尺度依次减小，直到覆盖所有的UV值。当k=0时，滤波器为高斯滤波。其尺度为最长UV的20%，从成像来看，这也是最大的尺度。k=1、2、3、4时，滤波器在频率域呈现一个环状，实现一个带通滤波。无论分成多少级，但本质上要确保这些滤波器和所合成的滤波器在频率上连续，覆盖全部UV点，如图1和图2中的最后一张图。在使用高斯滤波的情况下，最终合成的滤波器在边缘逐渐减小的滤波性质，正好减小了由于较大的UV处的稀疏导致的成图不可靠信息。

3　算法的图形处理器实现及效率

3.1　图形处理器实现

明安图射电频谱日像仪成图的图形处理器实现采用GPU-CUDA(Compute Unified Device Architecture)架构，使用Python编程语言，PyCUDA，Scikit-CUDA作为CUDA并行编程接口库。实现了数据处理的单机命令行模式及分布式运行模式[6]，其中数据预处理(数据格式转换[7]、相位校准、异常数据标记等)在中央处理器环境下运行，成图及洁化在图形处理器环境下运行，包含8台服务器的高性能计算机群(CPUs + GPUs)已在观测站搭建好。成像流程如图3，其中日像仪多尺度洁化算法描述如表1。以低频阵观测数据 “2015-11-01 04:08:49. 354161240(UTC), 1.712 5 GHz，右旋” 为例，最终洁化图像如图4(a)。在实验过程中，常用天文学软件应用包(the Common Astronomy Software Applications package, CASA)用于验证数据处理各个阶段的正确性。将数据处理系统生成的UVFITS文件通过CASA转换成MS文件并使用多尺度洁化成图可得到基本一致的成图结果。为说明当前处理的正确性，文中给出了日本野边山天文台已经公布的观测处理结果图4(b)作为对比。由于当前的数据处理过程中没有进行日轴方位角改正，图像中存在一些误差，在后期的工作中将进行改正。

图1滤波器设计
Fig.1Filters profiles alongu(-) andv(-) axes

图2UV滤波器设计
Fig.2Filters in the UV-plane

图3基于图形处理器的成像流程
Fig.3GPU-based implementation of imaging

表1　明安图射电频谱日像仪多尺度洁化Table 1　Multi-Scale CLEAN for MSUER

3.2　算法效率

由于由多个滤波器得到的图像同时进行洁化，每个滤波后的图像有不同的迭代次数，且洁化算法的迭代次数动态改变，因此本文仅以具体实例说明多尺度带通滤波的性能。在NVIDIA Corporation GM200测试环境下，对于处理图4中的2015-11-01 1.712 5 GHz的一张1 024 × 1 024的图像，使用基于带通滤波的多尺度洁化时间为1.424 s，前期实验中的Högbom CLEAN算法[8]的洁化时间为4.037 48 s。经过多次洁化处理实验，表明当前的多尺度洁化能有效将成图效率提高近3倍，同时可以减少迭代次数，对数据处理流水线性能的提升有重要意义。

图4(a) 明安图射电频谱日像仪洁化图(Multi-scale)(MUSER-I at 2015-11-01 04:08:49. 354161240 (UTC)，频率: 1.712 5 GHz，右旋)； (b) 野边山天文台洁化图(2015-11-01 at 04: 10: 00 (UTC), 频率: 17 GHz (R + L)

Fig.4(a) Clean image of MUSER on 2015-11-01 at 04:08:49. 354161240 (UTC), frequency: 1.7125GHz, polarization: right); (b) Clean image from Nobeyama on 2015-11-01 at 04:10:00 (UTC), frequency: 17GHz (R + L)

4　总　结

本文在分析多尺度洁化算法原理的基础上，实现了基于带通滤波的多尺度洁化算法，较文[3]的多尺度洁化算法在全局不同尺度的图像中的迭代而言，完全实现了不同尺度的并行。实验选取了针对明安图射电频谱日像仪成像的滤波器、尺度大小及不同尺度的取值。同时，基于图形处理器的实现表明，基于带通滤波的多尺度洁化算法比同等条件下的Högbom CLEAN性能提高了近3倍，并且可有效减少迭代次数。本文的研究结果将有效提高整个数据处理系统的性能，并为射电干涉成像提供了有价值的参考。