2. 广州大学天体物理中心/物理与电子工程学院, 广东 广州 510006;
3. 中国科学院云南天文台, 云南 昆明 650011
2. Astrophysics Center/Institute of Physics and Electronic Engineering, Guangzhou University, Guangzhou 510006, China;
3. Yunnan Observatories, Chinese Academy of Sciences, Kunming 650011, China
中国新一代厘米-分米波射电日像仪——明安图射电频谱日像仪(MingantU SpEctral Radioheliograph, MUSER)是基于综合孔径成像原理实现在频率0.4~15.0 GHz对太阳进行高分辨率(时间、空间)观测并成像的射电望远镜[1-2]。
明安图射电频谱日像仪主要由天线(天线阵)、接收系统、数据处理系统3部分组成。天线阵由低频阵(MUSER-Ⅰ，共40面天线)和高频阵(MUSER-Ⅱ，共60面天线)两部分组成，天线阵的拓扑结构为螺旋阵[3-4]。在观测过程中，低频阵和高频阵分别每3 ms产生一个数据帧，每个数据帧分别用100 KB和256 KB存储，数据流量分别为32 MB/s和86 MB/s，以每天观测10小时为例，每日原始观测数据量分别约为1.2 TB和3.2 TB[5-6]。
如此海量的天文数据，对于数据处理工作而言是一个巨大的挑战，但挑战与机遇并存。在前期研究工作中，为了满足数据处理的需求，研究人员一方面采用分布式计算技术，设计并实现了一套分布式数据处理执行框架OpenCluster[7]；另一方面采用并行计算技术并借助高性能计算设备，专注于底层算法研究，实现了一套高效的射电干涉阵数据处理管线[8]。特别地，采用消息传递接口(Message Passing Interface, MPI)、统一计算设备架构(Compute Unified Device Architecture, CUDA)以及开放运算语言(Open Computing Language, OpenCL)技术在天文领域有相关应用，并已经取得一定的成果[9-11]，证明这些技术稳定性与性能可以满足天文应用软件开发的要求。
为了便于明安图射电频谱日像仪数据处理系统的开发、部署、应用以及推广，同时满足单机环境下仍旧能够进行高性能的成像数据处理工作，在前期对OpenCL技术研究的基础上，进一步采用OpenCL技术对射电干涉阵成像中的网格化算法进行并行优化，在保证网格化算法执行效率的同时，提升算法对各种硬件平台的适应性，便于后期数据处理软件的应用与推广。
1 基于OpenCL的网格化算法实现 1.1 网格化算法射电干涉阵成像过程是基于综合孔径成像原理对采样的可见度数据进行网格化、快速傅里叶变换(Fast Fourier Transformation, FFT)、去卷积(Deconvolution)等操作，最后产生射电源的观测图像[12]。
为了借助离散傅里叶变换(Discrete Fourier Transform, DFT)的快速方法——快速傅里叶变换的计算效率进行快速成像，射电干涉阵采集的可见度数据必须落在一个均匀划分的网格上(如图 1(a))。实际上，由于UV采样点分布不均匀，干涉仪得到的是非均匀采样的可见度数据(如图 1(b))，将非均匀采样的可见度数据插值到均匀划分的网格点上的过程称为网格化[13-14]。
|
| 图 1 均匀划分的网格(a)与非均匀分布的可见度数据示意图(b) Fig. 1 Regular grid (left) and the schematic diagram of non-uniform visibility data (right) |
将非均匀采样可见度数据转换成均匀采样可见度数据，即对未知像素点的数据进行重采样，主要有两类方法：内插法和卷积核法。一些内插值方法被开发用于重建非均匀采样，如最邻近插值、双线性插值等。文[15]表明，最优的网格化方法是卷积一个无限扩展的sinc函数，即卷积核法，出于实际考虑，有限的卷积函数取代了无限扩展的sinc函数。常用的网格卷积函数(Gridding Convolution Function, GCF)有截断sinc函数、截断指数函数、拟合的球面波函数。文[16]讨论了这些网格卷积函数在不同尺寸大小下的性能。在日像仪数据处理中，选取大小为6的拟合球面波函数作为网格卷积函数。卷积网格化并成像过程如图 2。
|
| 图 2 卷积网格化并成像 Fig. 2 Convolutional gridding and imaging |
根据OpenCL编程规范，OpenCL程序主要由2部分组成：(1)采用OpenCL语言(类C语言)编写的内核函数，执行在OpenCL设备(多核CPU, GPU, Cell, DSP, FPGA等)上；(2)通过调用OpenCL的应用编程接口(Application Programing Interface, API)管理内核函数的主机程序，执行在主机(CPU)上。
本文将网格化过程中的关键步骤，由OpenCL语言编写成内核函数，以便在OpenCL设备上并行执行，并通过在主机上调用这些内核函数，从而实现网格化算法的并行化。主要改写的内核函数有：gridding_kernel，correct_grid_kernel，normalize_kernel，point_symmetric_kernel，shift_kernel等。此外，成像过程的傅里叶变换是通过调用基于OpenCL实现的快速傅里叶变换，包含于python的pyfft包中。网格化并行执行流程如图 3。
|
| 图 3 网格化并行执行流程图 Fig. 3 The flow diagram of parallel gridding |
本文处理的数据以二维数组的形式表示，在采用OpenCL语言编写内核函数过程中，由OpenCL提供的二维数组索引get_global_id(0)和get_global_id(1)，唯一标识二维数组中的元素，为方便计算，需借助传递给内核函数的二维数组的参数(宽度或高度)，将二维数组按行或按列转换成一维数组，这样，二维索引转换成一维索引id，当内核函数在并行设备中执行时，每个线程处理由索引id唯一标识的二维数组元素。OpenCL二维索引转一维索引如图 4。
|
| 图 4 OpenCL二维索引转一维索引示意图 Fig. 4 The schematic diagram of OpenCL two-dimensional index to one-dimensional index |
在网格化执行过程中，由于可见度数据共轭对称，在对可见度数据进行插值时，为减少计算量，先对一半的可见度数据进行插值，再通过中心点对称方式计算得到另一半数据。将中心点对称操作编写成内核函数point_symmetric_kernel，如图 5。
|
| 图 5 中心点对称操作示意图 Fig. 5 The schematic diagram of point symmetric operation |
| point_symmetric_kernel伪代码如下： |
| OpenCL kernel function: point_symmetric_kernel |
| Input parameter: a two dimensional array im, a flag hfac, the two dimensional array height H and width W |
| //Obtain two dimensional array index (ix, iy) |
| ix←get_global_id(0) |
| iy←get_global_id(1) |
| // Stay within the bounds |
| If ix > 0 and ix < H and iy > 0 and iy < W/2 then |
| nix← H-ix |
| niy← W-iy |
| im[ix* W + iy]. x←im[nix* W + niy]. x |
| im[ix* W + iy]. y←hfac*im[nix* W + niy]. y |
| Output: a new two dimensional array im |
在傅里叶变换过程中，为使变换后的图像是一个完整的周期，需要进行平移操作。将平移操作编写成内核函数shift_kernel，如图 6。
|
| 图 6 傅里叶变换平移操作示意图 Fig. 6 The schematic diagram of Fourier transform shift operation |
在网格化执行过程中，由于可见度数据与一个网格卷积函数进行了卷积，为消除网格卷积函数的影响，需进行网格校正操作。将网格校正操作编写成内核函数correct_grid_kernel。进行网格校正后，需要对脏图和脏束进行归一化或标准化，即将脏图和脏束中所有元素除以脏束中的最大值，从而保证脏束的峰值为1以及脏图与脏束具有相同量纲，便于后续进行洁化。
1.3 并行加权由于UV采样点太少，对可见度函数采样值进行傅里叶逆变换得到的图像(脏图)包含一些虚假信息，通过对UV采样点赋予不同的权值改善图像质量的操作称之为加权[17]。常用的加权方式有自然加权、统一加权以及Taper，本文将这3种加权方式改写成OpenCL内核函数weight_natural_kernel，weight_uniform_kernel和weight_taper_kernel。
| 3种加权方法伪代码如下： |
| OpenCL kernel function: weight_natural_kernel、weight_uniform_kernel、weight_taper_kernel |
| Input parameter: a two dimensional array nim, the number of samples falling into cell cnt, the two dimensional array height H and width W//Obtain two dimensional array index (ix, iy) ix←get_global_id(0) iy←get_global_id(1) id←iy + ix* W //Stay within the bounds If ix > 0 and ix < H and iy > 0 and iy < W then if cnt[id]!=0 then |
|
| nim[id].
x←nim[id].
x*wgt nim[id]. y←nim[id]. y*wgt Output: weighted nim |
自然加权赋予相应像素点的权重为1，统一加权赋予相应像素点的权重为区域内采样点频数的倒数，Taper加权赋予相应像素点的权重为对应高斯函数值。通过改变像素点的权重大小，对可见度数据进行加权操作，从而改善图像质量。
2 实验结果 2.1 实验环境OpenCL定义不同的公司或厂商为不同的平台(Platform)，每个公司或厂商提供不同的OpenCL设备(Device)。本文选用的OpenCL设备有NVIDIA的GPU (GeForce GTX TITAN X, Tesla k20m)和Intel的CPU (Intel Xeon E5-2620 V2)，平台及设备参数见表 1。
| Parameter | Intel CPU | GTX TITAN X GPU | Tesla k20m GPU |
| Platform vendor | Intel(R) Corporation | NVIDIA Corporation | NVIDIA Corporation |
| Platform version | OpenCL 1.2 LINUX | OpenCL 1.2 CUDA 7.5.18 | OpenCL 1.2 CUDA 7.5.18 |
| Device name | Intel(R) Xeon(R) CPU E5-2620 V2@2.10 GHz | GeForce GTX TITAN X | Tesla k20m |
| Device type | CPU | GPU | GPU |
| Device max clock speed | 2 100 MHz | 1 215 MHz | 705 MHz |
| Device compute units | 24 | 24 | 13 |
| Device cores | 6 | 2 496 | 2 496 |
| Device max work group size | 8 192 | 1 024 | 1 024 |
实验数据来源于明安图射电频谱日像仪2015年11月1日12时8分49秒354毫秒对太阳的观测数据，原始观测数据经过一系列预处理(坏天线标记、星表查询等)合成为标准天文数据格式(UVFITS文件)，数据文件被命名为20151101-120849_354161240.uvfits。通过数据处理系统执行网格化和快速傅里叶变换等操作，生成相应时刻太阳亮度分布图像(未经过去卷积操作的脏图, 如图 7(a))以及对应的点扩散函数(Point Spread Function, PSF)，点扩散函数也称为脏束(如图 7(b))。
|
| 图 7 MUSER成的脏图(a)和脏束(b)，观测时间是2015年11月1日12: 08: 49: 354，频率是1.712 5 GHz，极化方式是右旋 Fig. 7 The dirty image (a) and dirtybeam (b) of MUSER at 12: 08: 49: 354 on November 11th, 2015 (Frequency: 1.7125GHz, Polarization: right) |
在同一台服务器上，分别在中央处理器(Intel Xeon E5-2620 V2)、图形处理器(GeForce GTX TITAN X, Tesla k20m)并行设备上执行基于OpenCL实现的网格化算法和在图形处理器(GeForce GTX TITAN X, Tesla k20m)上执行基于CUDA实现的网格化算法，并进行快速傅里叶变换成图，分别生成大小为1 024 × 1 024 pixel和2 048 × 2 048 pixel的脏图。网格化并成图(Gridding + IFFT)过程执行时间见表 2。
| 设备 | 图像大小/pixel | 执行时间/s |
| CPU OpenCL | 1 024 × 1 024 | 0.904 |
| 2 048 × 2 048 | 4.407 | |
| GPU OpenCL(TITAN X) | 1 024 × 1 024 | 0.350 |
| 2 048 × 2 048 | 0.713 | |
| GPU OpenCL(k20m) | 1 024 × 1 024 | 0.372 |
| 2 048 × 2 048 | 0.782 | |
| GPU CUDA(TITAN X) | 1 024 × 1 024 | 0.323 |
| 2 048 × 2 048 | 0.610 | |
| GPU CUDA(k20m) | 1 024 × 1 024 | 0.344 |
| 2 048 × 2 048 | 0.760 |
从表 2可以看出，执行基于OpenCL和基于CUDA实现的网格化算法并成图，在图形处理器环境下消耗的时间基本相当，同时，基于OpenCL实现的网格化算法还能在纯中央处理器环境下较快速地执行。
2.3 讨论本文采用并行计算OpenCL技术，基于Python语言导入PyOpenCL与PyFFT扩展包，实现了对日像仪成像过程中的网格化算法以及快速傅里叶变换成图过程。基于OpenCL实现的网格化算法具有如下优点：
(1) 能够在图形处理器环境中执行，并且不局限于NVIDIA GPU (实验条件限制，只选用了NVIDIA GPU和Intel CPU)，执行效率与基于CUDA实现的网格化算法执行效率大致相当，保证了日像仪成图的性能；
(2) 能够在多核中央处理器环境下执行，适用于在单机环境下进行开发与测试，便于软件的推广与应用。
综上所述，采用OpenCL实现的并行网格化算法，在保证算法执行效率的同时有效地提升了算法对硬件平台的适应性。此外，本工作进一步验证了OpenCL在天文软件开发中的可用性，从CUDA到OpenCL，异构系统从NVIDIA GPU + CPU的异构模式转变成并行设备加中央处理器的异构模式。这种异构模式的转变将扩展并行计算在天文高性能软件开发中的应用。
3 下一步工作本文研究的射电干涉阵成像过程局限于小视场，在大视场成像中w项引起的视场扭曲是不可忽略的，传统的二维快速傅里叶变换成像将不再适用。因此，下一步将在此研究工作的基础上，采用OpenCL技术对大视场成像中的关键算法，例如w-projection，faceting以及大视场混合算法，进行并行优化。
致谢: 感谢国家天文台-阿里云天文大数据联合研究中心对本项工作的支持。
| [1] | YAN Y, ZHANG J, WANG W, et al. The Chinese Spectral Radioheliograph-CSRH[J]. Earth Moon & Planets, 2009, 104(1-4): 97–100. |
| [2] | 赖铖, 梅盈, 邓辉, 等. MUSER可见度数据积分方法与实现[J]. 天文研究与技术, 2018, 15(1): 78–86 |
| [3] | 周鑫磊, 王威, 王锋, 等. 基于QT的MUSER观测数据多屏图形化实时显示的设计与实现[J]. 天文研究与技术, 2015, 12(4): 503–509 DOI: 10.3969/j.issn.1672-7673.2015.04.017 |
| [4] | SHI C, WANG F, DENG H, et al. High performance negative database for massive data management system of the Mingantu Spectral Radioheliograph[J]. Publications of the Astronomical Society of the Pacific, 2017, 129(978): 084501(10pp). |
| [5] | DAI H M, MEI Y, WANG W, et al. An auto-flag method of radio visibility data based on support vector machine[J]. Chinese Astronomy & Astrophysics, 2017, 41(1): 125–135. |
| [6] | WANG F, DENG H, WANG W. High performance distributed data processing pipeline for Chinese Spectral RadioHeliograph[C]//Radio Science Conference. 2015. |
| [7] | WEI S L, WANG F, DENG H, et al. OpenCluster:a flexible distributed computing framework for astronomical data processing[J]. Publications of the Astronomical Society of the Pacific, 2016, 129(972): 024001(14pp). |
| [8] | WANG F, MEI Y, DENG H, et al. Distributed data-processing pipeline for Mingantu Ultrawide Spectral Radioheliograph[J]. Publications of the Astronomical Society of the Pacific, 2015, 127(950): 383–396. DOI: 10.1086/680852 |
| [9] | 陈泰燃, 王威, 王锋, 等. 基于MPI的高性能UVFITS数据合成研究与应用[J]. 天文研究与技术, 2016, 13(2): 184–189 DOI: 10.3969/j.issn.1672-7673.2016.02.007 |
| [10] | MEI Y, WANG F, WANG W, et al. GPU-based high-performance imaging for Mingantu Spectral Radio Heliograph[J]. Publications of the Astronomical Society of the Pacific, 2018, 130(983): 014503(11pp). |
| [11] | 冯勇, 陈坤, 邓辉, 等. 基于OpenCL的MUSER CLEAN算法研究与实现[J]. 天文学报, 2017, 58(2): 55–64 |
| [12] | HÖGBOM J A. Aperture synthesis with a non-regular distribution of interferometer baselines[J]. Astronomy & Astrophysics Supplement, 1974, 15(15): 417–426. |
| [13] | THOMPSON A R, BRACEWELL R N. Interpolation and fourier transformation of fringe visibilities[J]. Astronomical Journal, 1973, 79(1): 11–24. |
| [14] | SEDARAT H, NISHIMURA D G. On the optimality of the gridding reconstruction algorithm[J]. IEEE Transactions on Medical Imaging, 2000, 19(4): 306–317. DOI: 10.1109/42.848182 |
| [15] | O'SULLIVAN J D. A fast sinc function gridding algorithm for fourier inversion in computer tomography[J]. IEEE Transactions on Medical Imaging, 1985, 4(4): 200–207. DOI: 10.1109/TMI.1985.4307723 |
| [16] | JACKSON J I, MEYER C H, NISHIMURA D G, et al. Selection of a convolution function for fourier inversion using gridding[J]. IEEE Transactions on Medical Imaging, 1991, 10(3): 473–478. DOI: 10.1109/42.97598 |
| [17] | BOONE F. Weighting interferometric data for direct imaging[J]. Experimental Astronomy, 2013, 36(1-2): 77–104. DOI: 10.1007/s10686-012-9322-1 |