吳懷廣 劉琳琳 石永生 李代祎 謝鵬杰

中圖分類號(hào)：F222.3文獻(xiàn)標(biāo)識(shí)碼：ADOI：10.3969/j.issn.2096-1553.2019.02.011

文章編號(hào)：2096-1553（2019）02-0082-06

關(guān)鍵詞：ARL;并行化算法;Gridding算法;CUDA

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

摘要：針對(duì)海量天文數(shù)據(jù)實(shí)時(shí)性處理效率低的問(wèn)題，通過(guò)對(duì)SKA圖像采集及成像ARL算法庫(kù)中耗時(shí)較長(zhǎng)的Gridding算法進(jìn)行耗時(shí)分析，找出了該算法中調(diào)用頻率高且運(yùn)行時(shí)間長(zhǎng)的兩個(gè)函數(shù)convolutional-grid和convolutional-degrid，利用GPU的多線程并行化處理降低兩個(gè)函數(shù)的循環(huán)迭代，實(shí)現(xiàn)了Gridding算法在GPU和CPU上的協(xié)同運(yùn)行.驗(yàn)證實(shí)驗(yàn)結(jié)果表明，在相同的數(shù)據(jù)量下，改進(jìn)后的Gridding算法運(yùn)行時(shí)間大大縮短，特別是在處理海量數(shù)據(jù)時(shí)，有效提高了ARL的整體運(yùn)行效率.

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 引言

在過(guò)去的幾十年里，射電望遠(yuǎn)鏡的靈敏度和圖像分辨率均有很大的提升.SKA[1-2]作為世界上最大綜合孔徑的射電望遠(yuǎn)鏡，采集數(shù)據(jù)的速率非?？?，數(shù)據(jù)采集量非常大，其設(shè)計(jì)目標(biāo)是要大于12 TB/s.海量數(shù)據(jù)的產(chǎn)生和天文成像本身對(duì)實(shí)時(shí)處理的嚴(yán)格要求，給計(jì)算機(jī)的計(jì)算能力帶來(lái)了巨大的挑戰(zhàn).海量數(shù)據(jù)的科學(xué)處理一般需要百億億次量級(jí)的超級(jí)計(jì)算機(jī)來(lái)完成.目前，通過(guò)提高主頻來(lái)提高CPU處理能力的傳統(tǒng)方式已受到集成電路集成度的制約，因此利用多核CPU和多核加速器（如GPU，Cell/BE等）處理大量時(shí)效性數(shù)據(jù)成為發(fā)展趨勢(shì).在天文成像過(guò)程中，網(wǎng)格化和去網(wǎng)格化是最耗時(shí)的兩個(gè)操作.如果處理的可見度數(shù)據(jù)是EB量級(jí)或以上，則消耗的時(shí)間不能通過(guò)調(diào)整計(jì)算機(jī)性能來(lái)緩和.傳統(tǒng)的網(wǎng)格化一般都是在CPU上運(yùn)行，加速器在每一個(gè)浮點(diǎn)上的帶寬很小，應(yīng)用的時(shí)效性也相應(yīng)較低.所以，并行化算法作為提高計(jì)算速度的一個(gè)重要途徑，在射電天文學(xué)數(shù)據(jù)處理[3]方面越來(lái)越受到關(guān)注.

網(wǎng)格化（Gridding）算法的作用是將落在笛卡爾坐標(biāo)外的數(shù)據(jù)點(diǎn)插值到網(wǎng)格上.W.N.Brouw[4]最早提出了Gridding方法，用于實(shí)現(xiàn)極坐標(biāo)網(wǎng)格采樣的離散傅里葉變換.針對(duì)極坐標(biāo)網(wǎng)格采樣，J.D.OSullivan[5]提出了一種快速的sinc函數(shù)網(wǎng)格算法，并將其應(yīng)用在CT圖像重建上，該算法利用有限范圍的卷積函數(shù)得到sinc函數(shù)插值以減小計(jì)算時(shí)間，但是計(jì)算量太大.C.H.Meyer等[6]將Gridding算法應(yīng)用于磁共振成像（MRI）中的螺旋采樣，采用的卷積核是Kaiser-Bessel窗口函數(shù)，但是Kaiser-Bessel窗口函數(shù)不是最優(yōu)的卷積函數(shù).勞保強(qiáng)等[7]使用卷積函數(shù)網(wǎng)格化不規(guī)則的微波全息數(shù)據(jù)，分析了不同卷積函數(shù)的抗混疊性能，驗(yàn)證了球體函數(shù)是最具抗混淆性能的卷積函數(shù).A.L.Varbanescu等[8]以增加附加計(jì)算為代價(jià)，通過(guò)對(duì)數(shù)據(jù)的排序和搜索來(lái)改善空間局部性，進(jìn)而改善內(nèi)存的性能，但其結(jié)果尚未達(dá)到GPU運(yùn)算峰值的14%.

以上研究主要是在CPU環(huán)境下對(duì)Gridding算法的改進(jìn)，鮮見在ARL（SKA的算法參考庫(kù)，用以實(shí)現(xiàn)整體天文數(shù)據(jù)成像流程）中針對(duì)Gridding算法在GPU實(shí)現(xiàn)時(shí)的運(yùn)行效率的研究.鑒于此，本文擬分析CPU環(huán)境下ARL中耗時(shí)最長(zhǎng)的Gridding算法，利用性能分析工具找出Gridding算法中耗時(shí)最長(zhǎng)的函數(shù)模塊，然后對(duì)函數(shù)模塊使用CUDA進(jìn)行GPU并行執(zhí)行，以期利用GPU在并行運(yùn)算方面的優(yōu)勢(shì)，實(shí)現(xiàn)Gridding算法的CPU和GPU協(xié)同處理，從而提高ARL的整體運(yùn)行效率.

1 Gridding算法在ARL網(wǎng)格化中的應(yīng)用

在射電天文學(xué)中，所采集到的信號(hào)通常受到各類因素影響，數(shù)據(jù)呈現(xiàn)出不規(guī)律性.在利用射電望遠(yuǎn)鏡采集到的數(shù)據(jù)重建天空?qǐng)D像時(shí)，需要將不規(guī)則采樣的數(shù)據(jù)映射到標(biāo)準(zhǔn)的二維網(wǎng)格中，這稱為網(wǎng)格化.只有經(jīng)歷這一步，網(wǎng)格才可以通過(guò)快速傅里葉變換來(lái)重建天空?qǐng)D像.在射電望遠(yuǎn)鏡數(shù)據(jù)處理中，網(wǎng)格化是計(jì)算量最大、耗費(fèi)時(shí)間最長(zhǎng)的步驟.

在天文成像中，對(duì)于天空亮度I，初始光束模式A和能見度值V，存在以下傅里葉變換關(guān)系：

其中，S表示（u，v）采樣函數(shù)，V′表示觀測(cè)值.

對(duì)式②作數(shù)值估計(jì)的方法有直接傅里葉變換DFT（direct Fourier transform）和快速傅里葉變換FFT（fast Fourier transform）兩種方法.其中，DFT對(duì)N×N網(wǎng)格的每個(gè)點(diǎn)都做計(jì)算，通過(guò)蠻力計(jì)算求和來(lái)估計(jì)ID（l，m），其表達(dá)式為

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

FFT則通過(guò)利用DFT運(yùn)算中的對(duì)稱性和周期性，減少DFT的運(yùn)算量，從而更快地估計(jì)ID（l，m）.但是利用FFT需要先將數(shù)據(jù)插值到矩形網(wǎng)格點(diǎn)上，然后才可以應(yīng)用FFT.插值的過(guò)程就稱為網(wǎng)格化.ARL中Gridding算法的作用是將落在笛卡爾坐標(biāo)外的數(shù)據(jù)點(diǎn)插值到網(wǎng)格上，以便對(duì)網(wǎng)格化的數(shù)據(jù)進(jìn)行FFT計(jì)算.

天文學(xué)中的Gridding算法步驟如下.

1）密度補(bǔ)償函數(shù)與采樣數(shù)據(jù)相乘.Gridding算法首先要用密度補(bǔ)償函數(shù)與采樣數(shù)據(jù)相乘來(lái)彌補(bǔ)采樣數(shù)據(jù)的不均勻.網(wǎng)格算法中的采樣密度補(bǔ)償函數(shù)和插值函數(shù)（卷積函數(shù)）的選擇對(duì)圖像重建呈現(xiàn)的最后效果影響較大，也是網(wǎng)格化算法研究的核心問(wèn)題.

2）在網(wǎng)格點(diǎn)（uc，vc）處計(jì)算出卷積，即

3）對(duì)數(shù)據(jù)進(jìn)行重采樣，使數(shù)據(jù)落到網(wǎng)格點(diǎn)上.在所有網(wǎng)格點(diǎn)上對(duì)C×VW作采樣處理，即

4）應(yīng)用FFT對(duì)VR作傅里葉變換，得出經(jīng)過(guò)簡(jiǎn)單修正的初步圖像.VR是規(guī)則間隔的δ函數(shù)的線性組合，可以由離散傅里葉變換計(jì)算出其傅里葉變換VR的采樣矩陣，

5）為了消除卷積函數(shù)的影響，需要對(duì)重建后的圖像作進(jìn)一步的修正，即除以卷積函數(shù)的傅里葉變換.

2 Gridding算法的耗時(shí)分析

Gridding中一共含有14個(gè)函數(shù)，對(duì)這些函數(shù)進(jìn)行分析，選出調(diào)用頻率較高、運(yùn)行時(shí)間最長(zhǎng)的4個(gè)函數(shù)：1）grdsf函數(shù)，其功能是獲取網(wǎng)格函數(shù)及對(duì)網(wǎng)格進(jìn)行修正;2）anti_aliasing_calculate 函數(shù)，其功能是計(jì)算prolate球體反混疊函數(shù);3）convolutional_degrid 函數(shù)，其功能是利用被采樣的gcf部分uv坐標(biāo)值，進(jìn)行卷積解網(wǎng)格;4）convolutional_grid函數(shù)，其功能是由頻率和偏振獨(dú)立的gcf進(jìn)行卷積網(wǎng)格.Gridding算法中其他函數(shù)因其調(diào)用次數(shù)和單個(gè)函數(shù)耗時(shí)都很小，在整個(gè)過(guò)程中可以忽略不計(jì).

選出的4個(gè)函數(shù)被調(diào)用次數(shù)如圖1所示，執(zhí)行時(shí)間如圖2所示.分析中可見性數(shù)據(jù)量為104條.

由圖1可以看出，grdsf函數(shù)的被調(diào)用次數(shù)最高，達(dá)到370次，anti_aliasing_calculate函數(shù)被調(diào)用185次，convolutional_degrid函數(shù)被調(diào)用47次，convolutional_grid函數(shù)被調(diào)用44次.

由圖2可以看出，convolutional_grid和convolutional_degrid兩個(gè)函數(shù)運(yùn)行時(shí)間最長(zhǎng)，而grdsf和anti_aliasing_calculate兩個(gè)函數(shù)運(yùn)行時(shí)間較少.

綜合來(lái)看，需要對(duì)convolutional_grid和convolutional_degrid兩個(gè)函數(shù)進(jìn)行并行化處理.

3 Gridding算法的并行化實(shí)現(xiàn)

在整個(gè)SKA圖像采集和最終成像的過(guò)程中，Gridding算法是計(jì)算量大、耗時(shí)長(zhǎng)的算法之一.ARL算法庫(kù)是基于CPU運(yùn)行的，利用PyCUDA 進(jìn)行并行化加速可以使整個(gè)ARL的運(yùn)行時(shí)間得到有效縮短.

通常，GPU并行加速方式兩種：一種是在串行代碼的基礎(chǔ)上通過(guò)降低循環(huán)進(jìn)行加速;另一種是重新設(shè)計(jì)算法，將總模塊劃分為多個(gè)獨(dú)立的子模塊，通過(guò)分配線程來(lái)并行處理子任務(wù).本文采用第一種優(yōu)化措施.

這是因?yàn)镚ridding算法在處理數(shù)據(jù)的時(shí)候迭代次數(shù)很多，通過(guò)對(duì)Gridding算法的調(diào)用次數(shù)和耗時(shí)的分析可以看出，計(jì)算量最多的是convolutional_grid和convolutional_degrid函數(shù)，而這兩個(gè)函數(shù)運(yùn)算時(shí)for循環(huán)較多，且數(shù)據(jù)之間不存在依賴關(guān)系，所以可以直接利用多線程實(shí)現(xiàn)其在GPU下的并行計(jì)算，從而達(dá)到降低for循環(huán)次數(shù)，加快運(yùn)算速度的目的.

由于GPU與CPU之間不能直接通信，所以進(jìn)行并行化處理之前需將數(shù)據(jù)從主機(jī)端（CPU）拷貝到設(shè)備端（GPU），然后才能執(zhí)行相關(guān)的核函數(shù).同樣的，對(duì)于GPU并行化處理后的結(jié)果也需要從GPU拷貝到CPU.

本文利用全局存儲(chǔ)器進(jìn)行并行化處理.實(shí)驗(yàn)所用顯卡的靜態(tài)存儲(chǔ)器的大小是11 440 MB，因此所有的線程模塊中的線程不會(huì)產(chǎn)生訪問(wèn)沖突.利用CUDA[9]中的cudaMemcpy（）將數(shù)據(jù)載入到GPU的靜態(tài)存儲(chǔ)器中，以便所有的線程共享.根據(jù)線程索引獲取可見性數(shù)據(jù)進(jìn)行求和.等到線程都完成計(jì)算，最后將結(jié)果載入到內(nèi)存.進(jìn)行GPU并行處理的流程圖如圖3所示.

4 驗(yàn)證實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)環(huán)境：CPU 為Intel Xeon E5-2620 V3，內(nèi)存為816GB DDR4，內(nèi)存最大支持為768 GB，操作系統(tǒng)是CentOS 7.0，python版本為3.5;GPU采用Tesla K80，為雙核GK210架構(gòu)，主板型號(hào)為MG50-G20.

實(shí)驗(yàn)中選用的GPU一共含有4992個(gè)CUDA Cores，每個(gè)塊最多分配1024個(gè)線程，對(duì)于不同的數(shù)據(jù)量可以選擇不同的線程和線程塊.本實(shí)驗(yàn)使用的數(shù)據(jù)量分別為103條，104條，105條，106條，107條，108條，針對(duì)不同的數(shù)據(jù)量，線程塊統(tǒng)一采用（32，32，1）的形式，采用的網(wǎng)格塊的分配分別是（1，1，1），（10，1，1），（10，10，1），（100，10，1），（100，100，1），（1000，100，1）.

ARL中convolutional_grid函數(shù)和convolutional_degrid函數(shù)在CPU上運(yùn)行的時(shí)間與在GPU加速后的運(yùn)行時(shí)間見表1，其中加速比指在CPU上運(yùn)行時(shí)間和GPU上運(yùn)行時(shí)間的比值.

從表1可以看出，隨著數(shù)據(jù)量的增加，兩個(gè)函數(shù)的加速比不斷增加，但是加速比增加的速度逐漸減小.為了更直觀地表示加速比和數(shù)據(jù)量之間的關(guān)系，采用直方圖進(jìn)行表示，結(jié)果如圖4和圖5所示.

由圖4和圖5可以看出，當(dāng)數(shù)據(jù)量較少（103條）時(shí)，兩個(gè)函數(shù)的加速比分別是0.22和0.38，利用GPU并行處理的運(yùn)行時(shí)間反而增加，這是因?yàn)閿?shù)據(jù)量相對(duì)較少時(shí)，GPU器件的啟動(dòng)時(shí)間、內(nèi)存和顯存之間數(shù)據(jù)交互等都需要時(shí)間.因此在數(shù)據(jù)量比較小時(shí)，GPU因其內(nèi)部調(diào)度相對(duì)來(lái)說(shuō)耗時(shí)較長(zhǎng)，其并行加速并不占優(yōu)勢(shì).隨著數(shù)據(jù)量的增加，GPU加速的優(yōu)勢(shì)逐漸

顯現(xiàn).當(dāng)數(shù)據(jù)量大于104條時(shí)，加速比不斷增加.因此，GPU并行加速適用于大量數(shù)據(jù)的運(yùn)算情況，可在一定程度上提高ARL的整體運(yùn)行效率.

4 結(jié)語(yǔ)

本文針對(duì)SKA海量天文數(shù)據(jù)實(shí)時(shí)性處理效率較低的問(wèn)題，對(duì)ARL算法庫(kù)中Gridding算法進(jìn)行了耗時(shí)分析與優(yōu)化：通過(guò)分析Gridding算法中每個(gè)函數(shù)的調(diào)用次數(shù)和運(yùn)行時(shí)間，找出了需要優(yōu)化的兩個(gè)函數(shù)，即convolutional_grid和convolutional_degrid函數(shù);然后利用GPU的多線程并行化處理降低兩個(gè)函數(shù)的循環(huán)迭代，使得ARL中Gridding算法在GPU和CPU上實(shí)現(xiàn)了協(xié)同運(yùn)行.驗(yàn)證實(shí)驗(yàn)結(jié)果顯示，在數(shù)據(jù)量不影響輸入輸出數(shù)據(jù)格式和大小的情況下，對(duì)Gridding算法的并行化處理縮短了ARL的整體運(yùn)算時(shí)間，提高了運(yùn)行的效率，且數(shù)據(jù)量越大，加速優(yōu)勢(shì)越明顯.

本文所做的并行化研究對(duì)于天文圖像的成像過(guò)程具有一定的參考意義，對(duì)Gridding算法在多GPU下實(shí)現(xiàn)算法加速將是下一步的研究方向.

參考文獻(xiàn)：

[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] 彭曉明，郭浩然，龐建民.多核處理器——技術(shù)、趨勢(shì)和挑戰(zhàn)[J].計(jì)算機(jī)科學(xué)，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] 勞保強(qiáng)，王俊義，王錦清，等.基于卷積核網(wǎng)格化二維近程微波全息[J].微波學(xué)報(bào)，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.