劉慧慧 聞萌莎 錢(qián)慎一 吳懷廣 張偉偉 李代祎

中圖分類(lèi)號(hào)：TP301文獻(xiàn)標(biāo)識(shí)碼：ADOI：10.3969/j.issn.2096-1553.2019.02.012

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

關(guān)鍵詞：ARL;去卷積算法;CUDA;并行計(jì)算;Clean算法

Key words：ARL; deconvolution algorithm; CUDA; parallel computing;Clean algorithm

摘要：針對(duì)SKA算法參考庫(kù)ARL中的去卷積算法運(yùn)行效率低、無(wú)法滿足海量數(shù)據(jù)實(shí)時(shí)處理的問(wèn)題，提出了CPU和GPU協(xié)同工作模式下的并行化Clean算法.該方法將Clean算法中可以并行計(jì)算的步驟利用多線程在GPU上并行執(zhí)行，將無(wú)法并行計(jì)算的步驟在CPU上串行執(zhí)行.驗(yàn)證實(shí)驗(yàn)結(jié)果表明，在數(shù)據(jù)逐漸增大的過(guò)程中，

并行化Clean算法比在CPU上的串行處理運(yùn)行時(shí)間顯著減少，當(dāng)圖達(dá)到4096像素×4096像素時(shí)，可以有10倍的提速.這說(shuō)明并行化Clean算法在處理海量數(shù)據(jù)時(shí)，能夠顯著提高運(yùn)算效率.

Abstract：The deconvolution algorithm in the ARL of the SKA algorithm reference library is inefficient and cannot meet the needs of real-time processing of massive data. The parallelized Clean algorithm in the cooperative working mode of CPU and GPU was proposed. The steps of parallel computing in Clean algorithm were executed in parallel on GPU using multi-threads， and the steps in the Clean algorithm that couldnt be parallelized were executed serially on the CPU. The results showed that the running time of parallel Clean algorithm under CPU and GPU cooperative mode was significantly shorter than that under CPU. When the image size reached 4096×4096， the parallel Clean algorithm GPU cooperative mode could be speeded up by 10 times， which showed that the parallel Clean algorithm could significantly improve the efficiency of operation when dealing with massive data.

0 引言

射電望遠(yuǎn)鏡是觀測(cè)和研究來(lái)自天體的射電波的基本設(shè)備，但天線的數(shù)量有限，從而導(dǎo)致空間頻率覆蓋不完整[1]，影響最終圖像的構(gòu)建.平方公里陣列SKA（square kilometre array）[2]是國(guó)際上建造的最大綜合孔徑射電望遠(yuǎn)鏡.與現(xiàn)有射電望遠(yuǎn)鏡相比，SKA的靈敏度提高10～100倍，測(cè)量速度提高105倍[3].ARL（algorithm reference library）算法參考庫(kù)是SKA的候選算法庫(kù)，其中去卷積算法[4]是ARL的一個(gè)重要組成，在射電天文成像中起著至關(guān)重要的作用.

去卷積算法可以消除圖像不完整、圖像模糊對(duì)觀測(cè)結(jié)果的影響，引發(fā)很多學(xué)者的關(guān)注和研究.在ARL中，去卷積算法采用的是Clean算法.1974年，J.A.Hgbom [5]首次提出了Clean算法，一種非線性去卷積的迭代方法，用以消除旁瓣干擾.該算法雖然消除了旁瓣的干擾，但需要消耗大量的時(shí)間，運(yùn)算效率較低.2004年，S.Bhatnagar等[6]提出了一種用于無(wú)線電干涉圖像的尺度敏感的反卷積（Asp-Clean）算法，將圖像建模為自適應(yīng)尺度像素的集合，可以更準(zhǔn)確地重建非對(duì)稱結(jié)構(gòu)的天空?qǐng)D像，但是該算法增加了算法復(fù)雜性和計(jì)算成本，其計(jì)算時(shí)間是Clean算法的3倍多.2008年，T.J.Cornwell[7]提出了一種Multiscale Clean算法，

該算法

可以更好地處理擴(kuò)展源，提高圖像質(zhì)量，但是比Clean算法需要的運(yùn)行時(shí)間更長(zhǎng).2011年，U.Rau等[8]在Multiscale Clean算法的基礎(chǔ)上，結(jié)合Multi-frequency算法，在更高的靈敏度和采樣頻率的情況下，重建了天空?qǐng)D像，但此算法仍存在運(yùn)算時(shí)間較長(zhǎng)、運(yùn)行效率較低的問(wèn)題.2016年，L.Zhang等[9]提出自適應(yīng)環(huán)路增益自適應(yīng)規(guī)模Clean（Algas-Clean）算法，該算法可以提供更準(zhǔn)確的模型，具有更強(qiáng)的收斂性，比Asp-Clean算法的運(yùn)算速度快了50%.2017年，J.Cheng 等[10]

為解決傳統(tǒng)Clean算法出現(xiàn)的問(wèn)題，

提出了一種小波Clean算法，該算法對(duì)小波濾波器的參數(shù)進(jìn)行了優(yōu)化，提高了圖像質(zhì)量，但沒(méi)有縮短運(yùn)行時(shí)間，未提高算法運(yùn)行效率.目前，大部分學(xué)者都是針對(duì)提高圖像質(zhì)量進(jìn)行研究，在運(yùn)算效率方面研究成果很少.然而，隨著SKA項(xiàng)目的實(shí)施，海量的天文數(shù)據(jù)迅速增加，如何有效、快速地處理數(shù)據(jù)成為當(dāng)今計(jì)算機(jī)領(lǐng)域研究的重點(diǎn)，這使去卷積算法運(yùn)算效率的提高變得刻不容緩.

鑒于此，本文擬通過(guò)對(duì)Clean算法進(jìn)行耗時(shí)分析，利用多線程來(lái)設(shè)計(jì)算法的CUDA核函數(shù)，提出一種基于CPU和GPU協(xié)同工作模式的并行化Clean算法，以期提高海量天文數(shù)據(jù)下的去卷積算法的運(yùn)算效率.1 Clean算法在ARL中的應(yīng)用Clean算法 [11]是一種非線性去卷積方法，專門(mén)用于處理不完整的射電數(shù)據(jù)，將其運(yùn)用于ARL中，可以消除圖像中旁瓣（天線方向圖通常都有兩個(gè)或多個(gè)瓣，其中輻射強(qiáng)度最大的瓣稱為主瓣，其余的瓣稱為旁瓣）的影響.若天空亮度分布用I（x，y）表示，其中x，y分別表示點(diǎn)源的橫縱坐標(biāo)，則與之對(duì)應(yīng)的復(fù)可見(jiàn)度函數(shù)用V（u，v）表示，其中u，v稱為空間頻率.在uv平面（基線投影所在的平面）中，u指向東，v指向北，其基線如圖1所示.

由于SKA中天線數(shù)目是固定的，所以u(píng)v采樣點(diǎn)的數(shù)目有限.將uv采樣點(diǎn)的分布函數(shù)稱為采樣函數(shù)，對(duì)采樣函數(shù)進(jìn)行傅里葉逆變換可得到臟束，將射電望遠(yuǎn)鏡經(jīng)過(guò)預(yù)處理后的數(shù)據(jù)進(jìn)行傅里葉逆變換可得到臟圖，通過(guò)一個(gè)迭代過(guò)程找到臟圖中的最大亮度點(diǎn)及其位置，將上述結(jié)果與擬合光束（高斯光束）進(jìn)行卷積即得到恢復(fù)的圖像.

式③說(shuō)明了臟圖是真實(shí)圖像和臟束的卷積.具有旁瓣的臟束稱為點(diǎn)擴(kuò)散函數(shù)PSF（point spread function）.Clean算法可以從已知的臟圖和PSF中得到圖像點(diǎn)Ip的一系列位置，其主要計(jì)算步驟如下：

1）對(duì)采樣函數(shù)進(jìn)行傅里葉逆變換，得到臟束;

2）將射電望遠(yuǎn)鏡經(jīng)過(guò)預(yù)處理的數(shù)據(jù)進(jìn)行傅里葉逆變換得到臟圖，找到臟圖中的最大值點(diǎn)，并記錄該坐標(biāo);

3）將臟束的中心點(diǎn)移到這個(gè)最大值點(diǎn)的位置上，且使臟束乘以因子γ;

4）從臟圖中減去該臟束;

5）重復(fù)步驟2）—4），直到剩余圖像的最大值小于給定的噪聲水平;

6）對(duì)一擬合光束和δ函數(shù)（被減去的最大值點(diǎn)的總和及其位置稱為δ函數(shù)圖）進(jìn)行卷積;

7）將步驟6）的結(jié)果加上剩余圖像，得到最終的清晰圖像.

2 Clean算法的GPU并行化實(shí)現(xiàn)

ARL中的Clean算法是基于CPU的，運(yùn)行時(shí)間較長(zhǎng)，本文根據(jù)GPU并行的特點(diǎn)，利用多線程，設(shè)計(jì)并實(shí)現(xiàn)基于GPU并行化的Clean算法.

2.1 Clean算法的并行化分析

研究表明，去卷積算法在數(shù)據(jù)處理過(guò)程中迭代計(jì)算耗時(shí)較長(zhǎng)，對(duì)計(jì)算機(jī)內(nèi)存和CPU要求較高.GPU[13]是具有數(shù)千個(gè)計(jì)算核心的大規(guī)模并行架構(gòu)，可以在有限時(shí)間內(nèi)執(zhí)行多次浮點(diǎn)運(yùn)算，因此在GPU上并行實(shí)現(xiàn)去卷積算法可以大大提升運(yùn)算速度，節(jié)省時(shí)間.NVIDIA公司提出了一種并行計(jì)算架構(gòu)CUDA（computer unified device architecture）[14]，通過(guò)CUDA C編程在GPU上實(shí)現(xiàn)并行計(jì)算，其運(yùn)算速度比CPU提高了幾倍甚至上百倍[15].

通過(guò)對(duì)ARL中deconvolution模塊的分析可知，在Clean算法的實(shí)現(xiàn)過(guò)程中，需要進(jìn)行多次迭代計(jì)算，而每次迭代都需要使用上一次迭代的結(jié)果，因此無(wú)法將Clean算法的整個(gè)過(guò)程寫(xiě)成內(nèi)核函數(shù)，但可以將其中的部分操作寫(xiě)成內(nèi)核函數(shù)，使其可以在GPU上進(jìn)行并行計(jì)算，以此來(lái)提高計(jì)算效率，而將無(wú)法進(jìn)行并行計(jì)算的操作在CPU上串行執(zhí)行.

并行處理主要分為任務(wù)并行處理和數(shù)據(jù)并行處理.基于CUDA的組織結(jié)構(gòu)，結(jié)合Clean的算法分析，本文將采用數(shù)據(jù)并行處理方式.并行算法的執(zhí)行過(guò)程從CPU開(kāi)始，當(dāng)內(nèi)核函數(shù)被調(diào)用時(shí)，執(zhí)行過(guò)程轉(zhuǎn)移到GPU設(shè)備上.另外，GPU的內(nèi)存分為全局存儲(chǔ)器、常量存儲(chǔ)器、共享存儲(chǔ)器和寄存器[16].在整個(gè)程序中，數(shù)據(jù)需要被Grid網(wǎng)格中的所有線程訪問(wèn)，因此需要將數(shù)據(jù)放置在全局存儲(chǔ)器中;將在臟圖中找到的最大值放入共享存儲(chǔ)器，方便存取;每一個(gè)點(diǎn)源數(shù)據(jù)都需要相同的高斯核，將高斯核放置在常量存儲(chǔ)器中，以此來(lái)加強(qiáng)程序的快速訪問(wèn).

2.2 Clean算法的并行化實(shí)現(xiàn)

Clean算法的并行化實(shí)現(xiàn)步驟如下.

1）初始化潔圖，并將臟束和臟圖作為形參傳入.

2）在GPU上調(diào)用CUDA編寫(xiě)的dirtyFabs-Max核函數(shù)，找到臟圖中的最大值點(diǎn)，并記錄該點(diǎn)的位置坐標(biāo).

3）將臟束的中心點(diǎn)移到這個(gè)最大值點(diǎn)的位置上，且使臟束乘以因子γ.

4）在GPU上調(diào)用CUDA編寫(xiě)的subPsf核函數(shù)，從臟圖中減去該臟束.

5）判斷剩余圖像的最大值是否小于給定的噪聲水平，若未滿足條件，則繼續(xù)進(jìn)行迭代計(jì)算;若滿足條件，則退出迭代，繼續(xù)向下執(zhí)行.

6）對(duì)一擬合光束和δ函數(shù)進(jìn)行卷積.

7）將步驟6得到的結(jié)果加上剩余圖像，得到最終的清晰圖像.

8）把GPU上經(jīng)過(guò)處理得到的數(shù)據(jù)傳給CPU并進(jìn)行圖像繪制.

Clean算法的步驟2主要是為了找到臟圖的最大值點(diǎn)及其位置，該操作可以采用并行方式進(jìn)行處理.初始化最大值點(diǎn)，將圖像中的每個(gè)點(diǎn)與該點(diǎn)進(jìn)行比較，若圖像中的點(diǎn)大于該點(diǎn)的值，則把圖像中的點(diǎn)賦值給最大值點(diǎn).因?yàn)檎业脚K圖中最大值點(diǎn)及其位置的操作都是比較操作，可以同時(shí)進(jìn)行，所以能夠并行計(jì)算.由于ARL中的Clean算法是采用Python語(yǔ)言實(shí)現(xiàn)的，因此本文采用PyCUDA來(lái)實(shí)現(xiàn)Clean算法的并行化.具體實(shí)現(xiàn)方式是內(nèi)核函數(shù)采用CUDA C編程，其余部分仍采用Python編程，通過(guò)PyCUDA使Python程序可以訪問(wèn)Nivida的 CUDA 并行計(jì)算API，從而達(dá)到GPU并行的目的.dirtyFabsMax核函數(shù)將Blocksize和Gridsize作為相關(guān)參數(shù)傳入，指定網(wǎng)格和線程塊的大小，通過(guò)計(jì)算Blocksize和Gridsize的乘積得出網(wǎng)格中線程的數(shù)量，進(jìn)而指定并行處理多線程運(yùn)算GPU上創(chuàng)建線程的個(gè)數(shù).另外，需要syncthreads函數(shù)來(lái)保證線程同步，否則由于線程執(zhí)行具有無(wú)序性和異步性，可能會(huì)出現(xiàn)元素覆蓋的問(wèn)題，導(dǎo)致運(yùn)行結(jié)果不正確.dirtyFabsMax核函數(shù)的偽代碼為

Clean算法中的步驟4是從臟圖中減去臟束，此操作可以進(jìn)行并行處理.因?yàn)榕K圖和臟束上的每個(gè)點(diǎn)是相互獨(dú)立的，可以同時(shí)進(jìn)行減法操作，而GPU可以為每一個(gè)元素分配一個(gè)線程，因此每一個(gè)元素的相減操作就可以并行執(zhí)行，提高了Clean算法的計(jì)算效率.編寫(xiě)subPsf核函數(shù)用來(lái)實(shí)現(xiàn)從臟圖中減去臟束的操作.因?yàn)閷?shí)驗(yàn)所用的GPU每個(gè)線程塊最多支持1024個(gè)線程，而處理大小為1024像素×1024像素的圖像則需要1024×1024個(gè)線程，所以需要在一個(gè)網(wǎng)格中劃分多個(gè)線程塊來(lái)得到1024×1024個(gè)線程，函數(shù)被調(diào)用時(shí)的執(zhí)行配置為<<<1024，1024>>>，這樣就可以滿足GPU為每一個(gè)數(shù)據(jù)元素分配一個(gè)線程.subPsf核函數(shù)的偽代碼如下：

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

驗(yàn)證實(shí)驗(yàn)所用平臺(tái)是高性能計(jì)算平臺(tái)，配備有型號(hào)為NVIDIA Tesla K80的GPU，內(nèi)存共128 GB，存儲(chǔ)器為512 GB的SSD硬盤(pán)和1 TB的SAS高速硬盤(pán)，操作系統(tǒng)為CentOS7，使用PyCUDA作為開(kāi)發(fā)語(yǔ)言.

實(shí)驗(yàn)的測(cè)試數(shù)據(jù)來(lái)源于ARL提供的天文測(cè)試數(shù)據(jù).

首先將ARL中的數(shù)據(jù)進(jìn)行預(yù)處理，生成臟圖和臟束，然后調(diào)用并行化后的Clean算法對(duì)臟圖進(jìn)行處理，生成潔圖，結(jié)果如圖2所示.從圖2可以看出，臟圖中大量的模糊點(diǎn)被清除，潔圖更清晰地反映出天空亮度的分布情況.

將ARL中的Clean算法用PyCUDA語(yǔ)言在CPU和GPU協(xié)同工作條件下實(shí)現(xiàn)，用Python語(yǔ)言在CPU條件下實(shí)現(xiàn)，

測(cè)試迭代次數(shù)均為104次，兩個(gè)實(shí)現(xiàn)過(guò)程的計(jì)算時(shí)間如表1所示，加速比是Clean算法在CPU下的運(yùn)行時(shí)間與在GPU下運(yùn)行時(shí)間的比值.

從表1可以看出，在CPU與GPU協(xié)同作用下，用PyCUDA實(shí)現(xiàn)的Clean算法對(duì)數(shù)據(jù)的處理速度更快，用時(shí)遠(yuǎn)遠(yuǎn)小于CPU下用Python實(shí)現(xiàn)Clean算法的處理數(shù)據(jù)時(shí)間.當(dāng)圖像大小為512像素×512像素時(shí)，數(shù)據(jù)量較小，只有26×104條左右，加速的效果不是很明顯，只有兩倍

左右的加速效果.但是當(dāng)圖像大小增大后，特別是大于4096像素×4096像素后，此時(shí)的數(shù)據(jù)量增加到107條以上，加速比急劇增大，加速效果有了明顯提升.當(dāng)圖像大小為8192像素×8192像素時(shí)，Clean算法在CPU上需要運(yùn)行30 min才能得出結(jié)果，但本文實(shí)現(xiàn)的并行算法只需要2 min左右就可以得到結(jié)果，可以達(dá)到10倍的提速，極大地節(jié)省了時(shí)間，提高了算法的運(yùn)行效率.

從并行化實(shí)現(xiàn)的Clean算法在不同圖像大小下的加速比可以明顯看出，當(dāng)圖像大小增大時(shí)，加速比增大，在GPU上運(yùn)行的加速效果非常顯著.

4 結(jié)論

針對(duì)SKA算法參考庫(kù)ARL中去卷積算法處理海量數(shù)據(jù)負(fù)擔(dān)過(guò)重，時(shí)間太長(zhǎng)的問(wèn)題，對(duì)ARL中Clean算法進(jìn)行并行化研究，提出了CPU和GPU

協(xié)同工作模式下的并行化Clean算法.將Clean算法中耗時(shí)較長(zhǎng)的兩個(gè)部分，即在臟圖中找到最大值點(diǎn)和從臟圖減去臟束，采用多線程并行處理，并使用syncthreads函數(shù)來(lái)保證線程同步，無(wú)法并行處理的步驟在CPU上串行執(zhí)行，從而完成對(duì)海量數(shù)據(jù)的處理.驗(yàn)證實(shí)驗(yàn)結(jié)果表明，處理同樣的天文數(shù)據(jù)時(shí)，本文基于GPU并行化的Clean算法與CPU上的串行Clean算法相比，可以達(dá)到10倍的提速.

本文使用的是單GPU與CPU協(xié)同工作，沒(méi)有實(shí)現(xiàn)算法的最大并行化，因此，今后的工作是將單GPU上的算法移植到多GPU上，以便進(jìn)一步優(yōu)化ARL中去卷積算法的并行化.

參考文獻(xiàn)：

[1] DABBECH A，F(xiàn)ERRARI C，MARY D，et al.Moresane：model reconstruction by synthesis-analysis estimators-a sparse deconvolution algorithm for radio interferometric imaging[J].Astronomy & Astrophysics，2015，576：7.

[2] BROEKEMA P C，VAN NIEUWPOORT R V，BAL H E.The square kilometre array science data processor：Preliminary compute platform design[J].Journal of Instrumentation，2015，10（7）：14.

[3] VAN HEERDEN E，KARASTERGIOU A，ROBERTS S J，et al.New approaches for the real-time detection of binary pulsars with the Square Kilometre Array （SKA）[C]∥ General Assembly and Scientific Symposium.Piscataway：IEEE，2014：1-4.

[4] LA CAMERA A，SCHREIBER L，DIOLAITI E，et al.A method for space-variant deblurring with application to adaptive optics imaging in astronomy[J].Astronomy & Astrophysics，2015，579：1.

[5] HGBOM J A.Aperture synthesis with a non-regular distribution of interferometer baselines[J].Astronomy and Astrophysics Supplement Series，1974，15：417.

[6] BHATNAGAR S，CORNWELL T J.Scale sensitive deconvolution of interferometric images-I：adaptive scale pixel （Asp） decomposition[J].Astronomy & Astrophysics，2004，426（2）：747.

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