ARL中Clean算法的并行化研究

劉慧慧 聞萌莎 錢慎一 吳懷廣 張偉偉 李代祎

中圖分類號：TP301文獻標識碼：ADOI：10.3969/j.issn.2096-1553.2019.02.012

文章編號：2096-1553（2019）02-0088-07

關鍵詞：ARL;去卷積算法;CUDA;并行計算;Clean算法

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

摘要：針對SKA算法參考庫ARL中的去卷積算法運行效率低、無法滿足海量數據實時處理的問題，提出了CPU和GPU協同工作模式下的并行化Clean算法.該方法將Clean算法中可以并行計算的步驟利用多線程在GPU上并行執行，將無法并行計算的步驟在CPU上串行執行.驗證實驗結果表明，在數據逐漸增大的過程中，

并行化Clean算法比在CPU上的串行處理運行時間顯著減少，當圖達到4096像素×4096像素時，可以有10倍的提速.這說明并行化Clean算法在處理海量數據時，能夠顯著提高運算效率.

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

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

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

該算法

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

為解決傳統Clean算法出現的問題，

提出了一種小波Clean算法，該算法對小波濾波器的參數進行了優化，提高了圖像質量，但沒有縮短運行時間，未提高算法運行效率.目前，大部分學者都是針對提高圖像質量進行研究，在運算效率方面研究成果很少.然而，隨著SKA項目的實施，海量的天文數據迅速增加，如何有效、快速地處理數據成為當今計算機領域研究的重點，這使去卷積算法運算效率的提高變得刻不容緩.

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

由于SKA中天線數目是固定的，所以uv采樣點的數目有限.將uv采樣點的分布函數稱為采樣函數，對采樣函數進行傅里葉逆變換可得到臟束，將射電望遠鏡經過預處理后的數據進行傅里葉逆變換可得到臟圖，通過一個迭代過程找到臟圖中的最大亮度點及其位置，將上述結果與擬合光束（高斯光束）進行卷積即得到恢復的圖像.

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

1）對采樣函數進行傅里葉逆變換，得到臟束;

2）將射電望遠鏡經過預處理的數據進行傅里葉逆變換得到臟圖，找到臟圖中的最大值點，并記錄該坐標;

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

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

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

6）對一擬合光束和δ函數（被減去的最大值點的總和及其位置稱為δ函數圖）進行卷積;

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

2 Clean算法的GPU并行化實現

ARL中的Clean算法是基于CPU的，運行時間較長，本文根據GPU并行的特點，利用多線程，設計并實現基于GPU并行化的Clean算法.

2.1 Clean算法的并行化分析

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

通過對ARL中deconvolution模塊的分析可知，在Clean算法的實現過程中，需要進行多次迭代計算，而每次迭代都需要使用上一次迭代的結果，因此無法將Clean算法的整個過程寫成內核函數，但可以將其中的部分操作寫成內核函數，使其可以在GPU上進行并行計算，以此來提高計算效率，而將無法進行并行計算的操作在CPU上串行執行.

并行處理主要分為任務并行處理和數據并行處理.基于CUDA的組織結構，結合Clean的算法分析，本文將采用數據并行處理方式.并行算法的執行過程從CPU開始，當內核函數被調用時，執行過程轉移到GPU設備上.另外，GPU的內存分為全局存儲器、常量存儲器、共享存儲器和寄存器[16].在整個程序中，數據需要被Grid網格中的所有線程訪問，因此需要將數據放置在全局存儲器中;將在臟圖中找到的最大值放入共享存儲器，方便存取;每一個點源數據都需要相同的高斯核，將高斯核放置在常量存儲器中，以此來加強程序的快速訪問.

2.2 Clean算法的并行化實現

Clean算法的并行化實現步驟如下.

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

2）在GPU上調用CUDA編寫的dirtyFabs-Max核函數，找到臟圖中的最大值點，并記錄該點的位置坐標.

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

4）在GPU上調用CUDA編寫的subPsf核函數，從臟圖中減去該臟束.

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

6）對一擬合光束和δ函數進行卷積.

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

8）把GPU上經過處理得到的數據傳給CPU并進行圖像繪制.

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

Clean算法中的步驟4是從臟圖中減去臟束，此操作可以進行并行處理.因為臟圖和臟束上的每個點是相互獨立的，可以同時進行減法操作，而GPU可以為每一個元素分配一個線程，因此每一個元素的相減操作就可以并行執行，提高了Clean算法的計算效率.編寫subPsf核函數用來實現從臟圖中減去臟束的操作.因為實驗所用的GPU每個線程塊最多支持1024個線程，而處理大小為1024像素×1024像素的圖像則需要1024×1024個線程，所以需要在一個網格中劃分多個線程塊來得到1024×1024個線程，函數被調用時的執行配置為<<<1024，1024>>>，這樣就可以滿足GPU為每一個數據元素分配一個線程.subPsf核函數的偽代碼如下：

3 驗證實驗結果與分析

驗證實驗所用平臺是高性能計算平臺，配備有型號為NVIDIA Tesla K80的GPU，內存共128 GB，存儲器為512 GB的SSD硬盤和1 TB的SAS高速硬盤，操作系統為CentOS7，使用PyCUDA作為開發語言.

實驗的測試數據來源于ARL提供的天文測試數據.

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

將ARL中的Clean算法用PyCUDA語言在CPU和GPU協同工作條件下實現，用Python語言在CPU條件下實現，

測試迭代次數均為104次，兩個實現過程的計算時間如表1所示，加速比是Clean算法在CPU下的運行時間與在GPU下運行時間的比值.

從表1可以看出，在CPU與GPU協同作用下，用PyCUDA實現的Clean算法對數據的處理速度更快，用時遠遠小于CPU下用Python實現Clean算法的處理數據時間.當圖像大小為512像素×512像素時，數據量較小，只有26×104條左右，加速的效果不是很明顯，只有兩倍

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

從并行化實現的Clean算法在不同圖像大小下的加速比可以明顯看出，當圖像大小增大時，加速比增大，在GPU上運行的加速效果非常顯著.

4 結論

針對SKA算法參考庫ARL中去卷積算法處理海量數據負擔過重，時間太長的問題，對ARL中Clean算法進行并行化研究，提出了CPU和GPU

協同工作模式下的并行化Clean算法.將Clean算法中耗時較長的兩個部分，即在臟圖中找到最大值點和從臟圖減去臟束，采用多線程并行處理，并使用syncthreads函數來保證線程同步，無法并行處理的步驟在CPU上串行執行，從而完成對海量數據的處理.驗證實驗結果表明，處理同樣的天文數據時，本文基于GPU并行化的Clean算法與CPU上的串行Clean算法相比，可以達到10倍的提速.

本文使用的是單GPU與CPU協同工作，沒有實現算法的最大并行化，因此，今后的工作是將單GPU上的算法移植到多GPU上，以便進一步優化ARL中去卷積算法的并行化.

參考文獻：

[1] DABBECH A，FERRARI 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.