ARL中Gridding算法的并行化实现

    吴怀广 刘琳琳 石永生 李代祎 谢鹏杰

    

    

    

    中图分类号：F222.3文献标识码：ADOI：10.3969/j.issn.2096-1553.2019.02.011

    文章编号：2096-1553（2019）02-0082-06

    关键词：ARL;并行化算法;Gridding算法;CUDA

    Key words：ARL;parallelization algorithm;Gridding algorithm;CUDA

    摘要：针对海量天文数据实时性处理效率低的问题，通过对SKA图像采集及成像ARL算法库中耗时较长的Gridding算法进行耗时分析，找出了该算法中调用频率高且运行时间长的两个函数convolutional-grid和convolutional-degrid，利用GPU的多线程并行化处理降低两个函数的循环迭代，实现了Gridding算法在GPU和CPU上的协同运行.验证实验结果表明，在相同的数据量下，改进后的Gridding算法运行时间大大缩短，特别是在处理海量数据时，有效提高了ARL的整体运行效率.

    Abstract：Aiming at the low real-time processing efficiency of massive astronomical data， through time-consuming analysis of gridding algorithm in SKA image acquisition and imaging ARL library， two functions of convolutional-grid and convolutional-degrid with high frequency and long running time were found out in this algorithm. Then， two functions were parallelized on GPU by multi-threading to realize the cooperative operation of gridding algorithm on GPU and CPU. The experimental results showed that under the same amount of data， the running time of the improved gridding algorithm was greatly shortened， especially when dealing with massive data， the overall running efficiency of ARL was effectively improved.

    0 引言

    在过去的几十年里，射电望远镜的灵敏度和图像分辨率均有很大的提升.SKA[1-2]作为世界上最大综合孔径的射电望远镜，采集数据的速率非常快，数据采集量非常大，其设计目标是要大于12 TB/s.海量数据的产生和天文成像本身对实时处理的严格要求，给计算机的计算能力带来了巨大的挑战.海量数据的科学处理一般需要百亿亿次量级的超级计算机来完成.目前，通过提高主頻来提高CPU处理能力的传统方式已受到集成电路集成度的制约，因此利用多核CPU和多核加速器（如GPU，Cell/BE等）处理大量时效性数据成为发展趋势.在天文成像过程中，网格化和去网格化是最耗时的两个操作.如果处理的可见度数据是EB量级或以上，则消耗的时间不能通过调整计算机性能来缓和.传统的网格化一般都是在CPU上运行，加速器在每一个浮点上的带宽很小，应用的时效性也相应较低.所以，并行化算法作为提高计算速度的一个重要途径，在射电天文学数据处理[3]方面越来越受到关注.

    网格化（Gridding）算法的作用是将落在笛卡尔坐标外的数据点插值到网格上.W.N.Brouw[4]最早提出了Gridding方法，用于实现极坐标网格采样的离散傅里叶变换.针对极坐标网格采样，J.D.OSullivan[5]提出了一种快速的sinc函数网格算法，并将其应用在CT图像重建上，该算法利用有限范围的卷积函数得到sinc函数插值以减小计算时间，但是计算量太大.C.H.Meyer等[6]将Gridding算法应用于磁共振成像（MRI）中的螺旋采样，采用的卷积核是Kaiser-Bessel窗口函数，但是Kaiser-Bessel窗口函数不是最优的卷积函数.劳保强等[7]使用卷积函数网格化不规则的微波全息数据，分析了不同卷积函数的抗混叠性能，验证了球体函数是最具抗混淆性能的卷积函数.A.L.Varbanescu等[8]以增加附加计算为代价，通过对数据的排序和搜索来改善空间局部性，进而改善内存的性能，但其结果尚未达到GPU运算峰值的14%.

    以上研究主要是在CPU环境下对Gridding算法的改进，鲜见在ARL（SKA的算法参考库，用以实现整体天文数据成像流程）中针对Gridding算法在GPU实现时的运行效率的研究.鉴于此，本文拟分析CPU环境下ARL中耗时最长的Gridding算法，利用性能分析工具找出Gridding算法中耗时最长的函数模块，然后对函数模块使用CUDA进行GPU并行执行，以期利用GPU在并行运算方面的优势，实现Gridding算法的CPU和GPU协同处理，从而提高ARL的整体运行效率.

    1 Gridding算法在ARL网格化中的应用

    在射电天文学中，所采集到的信号通常受到各类因素影响，数据呈现出不规律性.在利用射电望远镜采集到的数据重建天空图像时，需要将不规则采样的数据映射到标准的二维网格中，这称为网格化.只有经历这一步，网格才可以通过快速傅里叶变换来重建天空图像.在射电望远镜数据处理中，网格化是计算量最大、耗费时间最长的步骤.

    在天文成像中，对于天空亮度I，初始光束模式A和能见度值V，存在以下傅里叶变换关系：

    其中，S表示（u，v）采样函数，V′表示观测值.

    对式②作数值估计的方法有直接傅里叶变换DFT（direct Fourier transform）和快速傅里叶变换FFT（fast Fourier transform）两种方法.其中，DFT对N×N网格的每个点都做计算，通过蛮力计算求和来估计ID（l，m），其表达式为

    ID（l，m）=1M∑Mk=1V′（uk，vk）e2πi（ukl+vkm）

    FFT则通过利用DFT运算中的对称性和周期性，减少DFT的运算量，从而更快地估计ID（l，m）.但是利用FFT需要先将数据插值到矩形网格点上，然后才可以应用FFT.插值的过程就称为网格化.ARL中Gridding算法的作用是将落在笛卡尔坐标外的数据点插值到网格上，以便对网格化的数据进行FFT计算.

    天文学中的Gridding算法步骤如下.

    1）密度补偿函数与采样数据相乘.Gridding算法首先要用密度补偿函数与采样数据相乘来弥补采样数据的不均匀.网格算法中的采样密度补偿函数和插值函数（卷积函数）的选择对图像重建呈现的最后效果影响较大，也是网格化算法研究的核心问题.

    2）在网格点（uc，vc）处计算出卷积，即

    3）对数据进行重采样，使数据落到网格点上.在所有网格点上对C×VW作采样处理，即

    4）应用FFT对VR作傅里叶变换，得出经过简单修正的初步图像.VR是规则间隔的δ函数的线性组合，可以由离散傅里叶变换计算出其傅里叶变换VR的采样矩阵，

    5）为了消除卷积函数的影响，需要对重建后的图像作进一步的修正，即除以卷积函数的傅里叶变换.

    2 Gridding算法的耗时分析

    Gridding中一共含有14个函数，对这些函数进行分析，选出调用频率较高、运行时间最长的4个函数：1）grdsf函数，其功能是获取网格函数及对网格进行修正;2）anti_aliasing_calculate 函数，其功能是计算prolate球体反混叠函数;3）convolutional_degrid 函数，其功能是利用被采样的gcf部分uv坐标值，进行卷积解网格;4）convolutional_grid函数，其功能是由频率和偏振独立的gcf进行卷积网格.Gridding算法中其他函数因其调用次数和单个函数耗时都很小，在整个过程中可以忽略不计.

    选出的4个函数被调用次数如图1所示，执行时间如图2所示.分析中可见性数据量为104条.

    由图1可以看出，grdsf函数的被调用次数最高，达到370次，anti_aliasing_calculate函数被调用185次，convolutional_degrid函数被调用47次，convolutional_grid函数被调用44次.

    由图2可以看出，convolutional_grid和convolutional_degrid两个函数运行时间最长，而grdsf和anti_aliasing_calculate两个函数运行时间较少.

    综合来看，需要对convolutional_grid和convolutional_degrid两个函数进行并行化处理.

    3 Gridding算法的并行化实现

    在整个SKA图像采集和最终成像的过程中，Gridding算法是计算量大、耗时长的算法之一.ARL算法库是基于CPU运行的，利用PyCUDA 进行并行化加速可以使整个ARL的运行时间得到有效缩短.

    通常，GPU并行加速方式两种：一种是在串行代码的基础上通过降低循环进行加速;另一种是重新设计算法，将总模块划分为多个独立的子模块，通过分配线程来并行处理子任务.本文采用第一种优化措施.

    这是因为Gridding算法在处理数据的时候迭代次数很多，通过对Gridding算法的调用次数和耗时的分析可以看出，计算量最多的是convolutional_grid和convolutional_degrid函数，而这两个函数运算时for循环较多，且数据之间不存在依赖关系，所以可以直接利用多線程实现其在GPU下的并行计算，从而达到降低for循环次数，加快运算速度的目的.

    由于GPU与CPU之间不能直接通信，所以进行并行化处理之前需将数据从主机端（CPU）拷贝到设备端（GPU），然后才能执行相关的核函数.同样的，对于GPU并行化处理后的结果也需要从GPU拷贝到CPU.

    本文利用全局存储器进行并行化处理.实验所用显卡的静态存储器的大小是11 440 MB，因此所有的线程模块中的线程不会产生访问冲突.利用CUDA[9]中的cudaMemcpy（）将数据载入到GPU的静态存储器中，以便所有的线程共享.根据线程索引获取可见性数据进行求和.等到线程都完成计算，最后将结果载入到内存.进行GPU并行处理的流程图如图3所示.

    4 验证实验结果与分析

    实验环境：CPU 为Intel Xeon E5-2620 V3，内存为816GB DDR4，内存最大支持为768 GB，操作系统是CentOS 7.0，python版本为3.5;GPU采用Tesla K80，为双核GK210架构，主板型号为MG50-G20.

    实验中选用的GPU一共含有4992个CUDA Cores，每个块最多分配1024个线程，对于不同的数据量可以选择不同的线程和线程块.本实验使用的数据量分别为103条，104条，105条，106条，107条，108条，针对不同的数据量，线程块统一采用（32，32，1）的形式，采用的网格块的分配分别是（1，1，1），（10，1，1），（10，10，1），（100，10，1），（100，100，1），（1000，100，1）.

    ARL中convolutional_grid函数和convolutional_degrid函数在CPU上运行的时间与在GPU加速后的运行时间见表1，其中加速比指在CPU上运行时间和GPU上运行时间的比值.

    从表1可以看出，随着数据量的增加，两个函数的加速比不断增加，但是加速比增加的速度逐渐减小.为了更直观地表示加速比和数据量之间的关系，采用直方图进行表示，结果如图4和图5所示.

    由图4和图5可以看出，当数据量较少（103条）时，两个函数的加速比分别是0.22和0.38，利用GPU并行处理的运行时间反而增加，这是因为数据量相对较少时，GPU器件的启动时间、内存和显存之间数据交互等都需要时间.因此在数据量比较小时，GPU因其内部调度相对来说耗时较长，其并行加速并不占优势.随着数据量的增加，GPU加速的优势逐渐

    显现.当数据量大于104条时，加速比不断增加.因此，GPU并行加速适用于大量数据的运算情况，可在一定程度上提高ARL的整体运行效率.

    4 结语

    本文针对SKA海量天文数据实时性处理效率较低的问题，对ARL算法库中Gridding算法进行了耗时分析与优化：通过分析Gridding算法中每个函数的调用次数和运行时间，找出了需要优化的两个函数，即convolutional_grid和convolutional_degrid函数;然后利用GPU的多线程并行化处理降低两个函数的循环迭代，使得ARL中Gridding算法在GPU和CPU上实现了协同运行.验证实验结果显示，在数据量不影响输入输出数据格式和大小的情况下，对Gridding算法的并行化处理缩短了ARL的整体运算时间，提高了运行的效率，且数据量越大，加速优势越明显.

    本文所做的并行化研究对于天文图像的成像过程具有一定的参考意义，对Gridding算法在多GPU下实现算法加速将是下一步的研究方向.

    参考文献：

    [1] RAZAVI-GHODS N，ACEDO E D L，EL-MAKA-DEMA A，et al.Analysis of sky contributions to system temperature for low frequency SKA aperture array geometries[J].Experimental Astronomy，2012，33（1）：141.

    [2] ZHANG Y，BROWN A K.Bunny ear combline antennas for compact wide-band dual-polarized aperture array[J].IEEE Transactions on Antennas and Propagation，2011，59（8）：3071.

    [3] 彭晓明，郭浩然，庞建民.多核处理器——技术、趋势和挑战[J].计算机科学，2012，39（S3）：320.

    [4] BROUW W N.Aperture synthesis[J].Methods in Computational Physics，1975，14：131.

    [5] OSULLIVAN J D.A fast sinc function gridding algorithm for fourier inversion in computer tomography[J].IEEE Transactions on Medical Imaging，1985，4（4）：200.

    [6] MEYER C H，HU B S，NISHIMURA D G，et al.Fast spiral coronary artery imaging[J].Magnetic Resonance in Medicine，1992，28（2）：202.

    [7] 勞保强，王俊义，王锦清，等.基于卷积核网格化二维近程微波全息[J].微波学报，2014，30（5）：82.

    [8] VARBANESCU A L.On the effective parallel programming of multi-core processors[D].Romania：Universitatea Politehnica Bucuresti，2010.

    [9] CHENG J.CUDA by example：an introduction to general-purpose GPU programming[M].Boston：Addison-Wesley Professional，2010.