于昂，勞保強(qiáng)，王俊義，安濤

(1.桂林電子科技大學(xué) 信息與通信學(xué)院，桂林541004；2.中國科學(xué)院 上海天文臺，上海200030)

1 引言

為了更深入地探索宇宙奧秘，射電望遠(yuǎn)鏡朝著更高靈敏度、更高分辨率和更快巡天速度的方向發(fā)展。預(yù)計(jì)于2021年開建的平方公里陣(square kilometre array,SKA)射電望遠(yuǎn)鏡，建成后將成為人類有史以來最大的綜合孔徑射電望遠(yuǎn)鏡。該觀測設(shè)備將為宇宙演化、引力的本質(zhì)、生命起源、地外文明等重大前沿科學(xué)研究提供強(qiáng)有力的數(shù)據(jù)[1]。全數(shù)字化的SKA望遠(yuǎn)鏡，規(guī)模巨大，將產(chǎn)生空前龐大的數(shù)據(jù)。SKA第一階段(SKA1)每年將約有600 PB的數(shù)據(jù)產(chǎn)品傳輸?shù)絊KA區(qū)域中心，這些數(shù)據(jù)將有助于天文學(xué)家進(jìn)行深度處理與分析[2]。最終用于發(fā)表的SKA數(shù)據(jù)產(chǎn)品需要經(jīng)過無數(shù)次的反復(fù)處理與驗(yàn)證，因此SKA數(shù)據(jù)處理軟件與算法的運(yùn)行速度和有效性將大大影響天文學(xué)家的工作效率，甚至影響科學(xué)產(chǎn)出。

大視場成像是SKA1低頻陣列(SKA1-low)數(shù)據(jù)處理的重要環(huán)節(jié)，如SKA1的首要科學(xué)任務(wù)宇宙黎明的探索，將利用SKA1-low進(jìn)行宇宙再電離成像觀測；SKA1十大科學(xué)目標(biāo)之一宇宙磁場，也將利用SKA1-low進(jìn)行大規(guī)模偏振成像觀測。目前，應(yīng)用比較廣泛的大視場成像算法有faceting[3–5],w-projection[6],w-stacking[7],snapshots[8]和w-snapshot[9]。faceting和snapshots主要是利用數(shù)據(jù)分割技術(shù)，使得分割后的數(shù)據(jù)滿足二維傅里葉變換近似成像方法。這種分割技術(shù)的成像方法由于需要循環(huán)處理多份分割數(shù)據(jù)，將消耗較多的計(jì)算時(shí)間和CPU資源。其中，faceting又分為圖像域faceting和uv-faceting，圖像域faceting需要對分割后的數(shù)據(jù)圖像進(jìn)行拼接，將會有圖像邊緣效應(yīng)；snapshots最后也需要進(jìn)行圖像拼接，同樣存在類似的問題。uv-faceting由于分割后的數(shù)據(jù)都在同一個(gè)網(wǎng)格上柵格化，經(jīng)過二維逆傅里葉變換后將獲得整個(gè)視場的圖像，不存在邊緣效應(yīng)。

w-projection,w-stacking和w-snapshot這三種方法主要是通過修正w項(xiàng)來完成大視場成像，相比分割技術(shù)的方法，成像質(zhì)量相對高一些。w-projection主要通過可見度數(shù)據(jù)與w核函數(shù)的卷積，將三維可見度數(shù)據(jù)投影到w=0平面，最后利用二維傅里葉變換進(jìn)行成像。由于數(shù)據(jù)處理時(shí)，卷積核函數(shù)將存放在內(nèi)存中，w-projection的卷積核是三維的w核函數(shù)，因此相比傳統(tǒng)的二維卷積核將要消耗更多的內(nèi)存。w-stacking可以說是w-projection的升級版，該算法在圖像域進(jìn)行w項(xiàng)校正，并使用多線程技術(shù)優(yōu)化算法，減少了卷積核的內(nèi)存使用，但是需要使用大量的內(nèi)存存儲w層的傅里葉變換數(shù)據(jù)。由于使用多線程加速，運(yùn)行速度相比w-projection有很大的改善。w-stacking算法已經(jīng)封裝在wsclean軟件中，主要用于SKA先導(dǎo)望遠(yuǎn)鏡默奇森寬視場陣列(Murchison widefield array,MWA)[10]的數(shù)據(jù)處理。由于w核可以單獨(dú)計(jì)算，并根據(jù)內(nèi)存來選擇任意w平面數(shù)，使得w-projection比較容易與其他算法進(jìn)行結(jié)合使用。w-snapshot算法即是snapshots和w-projection結(jié)合的產(chǎn)物，結(jié)合后的算法改善了snapshots邊緣的畸變問題，相比w-projection有更好的可拓展性。w-snapshot主要封裝在ASKAPsoft軟件中，用于SKA先導(dǎo)望遠(yuǎn)鏡澳大利亞平方公里陣(Australian square kilometre array pathfinder,ASKAP)[11]的數(shù)據(jù)處理。

由于混合算法具備一定的優(yōu)勢，本文主要研究uv-faceting和w-projection的混合算法，助力SKA1科學(xué)數(shù)據(jù)處理。本文的混合算法，一方面，可以改善uv-faceting成像質(zhì)量問題，包括噪聲水平和動態(tài)范圍等；另一方面，可以減少內(nèi)存的使用，提高w-projection的可拓展性。為了提高該混合算法的運(yùn)行速度，本文提出了兩種并行優(yōu)化方法：基于多核CPU的并行算法和基于GPU的并行算法。

2 混合w-facets成像算法

本文的混合成像算法，主要是將uv-faceting算法[4,5]和w-projection算法[12]結(jié)合以完成低頻干涉陣列大視場成像，稱之為混合w-facets成像算法，其中facets表示分面。該算法首先采用uv-faceting中的視場和數(shù)據(jù)分割技術(shù)獲得分面數(shù)據(jù)；然后通過相位旋轉(zhuǎn)技術(shù)獲得每個(gè)分面的可見度數(shù)據(jù)，通過基線旋轉(zhuǎn)技術(shù)獲得每個(gè)分面的u和v坐標(biāo)數(shù)據(jù)。其次，對每個(gè)分面的數(shù)據(jù)進(jìn)行柵格化，柵格化的卷積核函數(shù)不再使用傳統(tǒng)的卷積核函數(shù)，而是w-projection算法中的w核函數(shù)。每個(gè)分面的柵格化處理均在同一個(gè)網(wǎng)格中進(jìn)行，因此柵格化后的數(shù)據(jù)進(jìn)行逆傅里葉變換，只輸出一個(gè)天空圖像，即一個(gè)圖像結(jié)果，也稱為臟圖。相比圖像域faceting算法，本文的方法無需進(jìn)行圖像拼接，避免了圖像拼接引起的圖像邊緣畸變，同時(shí)也減少了運(yùn)行時(shí)間。此外，與w-projection結(jié)合后再與只使用uv-faceting相比，本文算法的成圖質(zhì)量和圖像的動態(tài)范圍將會進(jìn)一步提高。

近年，Wortmann[13,14]提出的w-towers算法整體思想上與本文的w-facets算法類似，但具體設(shè)計(jì)思路不同。主要不同為：(1)w-towers的facets分割方法與圖像域faceting方法相同，而w-facets是基于uvw域分割的。(2)w-facets方法中每個(gè)facet的數(shù)據(jù)在同一個(gè)大網(wǎng)格上進(jìn)行柵格化，子網(wǎng)格實(shí)際上是按照facets數(shù)目進(jìn)行嚴(yán)格均分而獲得；而w-towers需要經(jīng)過更復(fù)雜的計(jì)算，且每個(gè)子網(wǎng)格是由所有facets的數(shù)據(jù)構(gòu)成的。(3)本文采用的是w-projection方法進(jìn)行混合，w-towers采用的是w-stacking方法。下面詳細(xì)介紹本文算法的公式推導(dǎo)過程。

在綜合孔徑成像原理中，經(jīng)過射電干涉儀測量獲得的可見度數(shù)據(jù)V(u,v,w)與天空亮度分布函數(shù)I(l,m)(即天空圖像)存在如下關(guān)系[15]：

其中，u,v,w是射電望遠(yuǎn)鏡的基線在u-v-w坐標(biāo)系下的投影坐標(biāo)，(l,m)為u軸和v軸觀測方向上的單位向量的方向余弦，i是虛數(shù)單位。

為了使用二維傅里葉變換近似方法，可以將視場切割為Nfacets個(gè)相同大小的小視場(也稱為分面)，當(dāng)Nfacets滿足下式條件，每個(gè)小視場將可以采用近似方法重建天空圖像：

其中，Bmax是低頻射電干涉陣列臺站兩兩之間距離(基線)的最長值，max表示取最大值，λ是參考頻率對應(yīng)的波長，D是臺站的直徑范圍。

每個(gè)小視場的可見度數(shù)據(jù)可以通過原始可見度數(shù)據(jù)進(jìn)行相位旋轉(zhuǎn)獲得，令第j個(gè)分面的相位中心坐標(biāo)為(lj,mj,nj)，則該分面的可見度數(shù)據(jù)Vj為[4,5]：

其中，Δlj=l-lj,Δmj=m-mj，(lj,mj,nj)可以利用分面相位中心的赤經(jīng)RA和赤緯DEC，以及原始相位中心的RA和DEC，通過坐標(biāo)轉(zhuǎn)換公式計(jì)算出來[5]。

式(3)末尾的2個(gè)平方根項(xiàng)，經(jīng)過泰勒一級近似展開，有：

將式(4)代入式(3)后，得：

其中，uj和vj是第j個(gè)小視場的基線坐標(biāo)，計(jì)算公式為：

隨后，每個(gè)小視場的可見度數(shù)據(jù)和(uj,vj)數(shù)據(jù)在同一個(gè)網(wǎng)格上進(jìn)行柵格化，再對整個(gè)網(wǎng)格數(shù)據(jù)進(jìn)行二維逆傅里葉變換，即獲得整個(gè)視場的臟圖。其中，柵格化處理時(shí)，采用w-projection的w核函數(shù)與每個(gè)小視場的可見度數(shù)據(jù)進(jìn)行卷積來完成。w核函數(shù)的表達(dá)式為[17]：

其中，C(u,v)是傳統(tǒng)的卷積核函數(shù)，一般選取性能較優(yōu)的球諧函數(shù)；表示卷積運(yùn)算；(u,v,w)是G(l,m,w)的傅里葉變換，G(l,m,w)的表達(dá)式為[6]：

從w核函數(shù)的表達(dá)式可以看出，柵格化的卷積將多一個(gè)維度，也就是這個(gè)w核函數(shù)將由多個(gè)二維卷積核組成。值得注意的是，每個(gè)二維卷積核的大小不固定，將會隨著w值和視場的增大而增大[6]。通常，在w=0時(shí)，二維卷積核的大小為7×7或者9×9[6]。由于一維卷積核函數(shù)是一種窗函數(shù)，所以一維卷積核函數(shù)的長度又稱為全支撐大小。當(dāng)w＞0時(shí)，全支撐大小可以由計(jì)算得到，其中Iwidth是臟圖的最長邊大小或完整圖像的最長邊大小，單位為rad。從上述分析可知，如果w平面數(shù)和圖像過大則消耗的內(nèi)存較多；w平面數(shù)取值過小，則將導(dǎo)致混疊效應(yīng)產(chǎn)生。因此，w平面數(shù)取值應(yīng)合理，我們給出每個(gè)分面較佳的平面數(shù)計(jì)算公式：

其中，wmax是w的最大取值。

uv-faceting與w-projection結(jié)合后，Nfacets在不滿足公式(2)的情況下也能完成大視場成像，因?yàn)楫?dāng)Nfacets=1時(shí)，即進(jìn)行w-projection成像。此外，Nwplanes與需要處理的圖像大小會隨著Nfacets的增加而減小，w核函數(shù)的最大全支撐大小也隨之減小，因此本文的混合算法將大大減少卷積核內(nèi)存的消耗。

Algorithm 1基于多核CPU的成像并行算法Input:V可見度數(shù)據(jù)，(u,v,w)基線坐標(biāo)，Nwplanes w平面數(shù)，網(wǎng)格數(shù)組柵格存放柵格化結(jié)果全零數(shù)組(完整圖像大小)Output:臟圖1:#pragma omp parallel for schedule(static)2:for j in Nfacets do 3: for r in row count do 4: for c in channel count do 5: V,(u,v,w),RA0,DEC0=read data(r,c);6: RAj,DECj=calculate phase(RA0,DEC0,Nfacets)7: Vj=phase facets rotation(V,u,v,w,RAj,DECj);8: unew,vnew=baseline rotation(u,v,w);9: uscale,vscale=scale uv(unew,vnew);10:windex=w/wmax×(Nwplanes-1)11://w核柵格化12: for sv in 0 to cfullsupport[windex]do 13: for su in 0 to cfullsupport[windex]do 14:conv weight=Cw-kernal[windex][sv][su];15: grid[vscale+sv][uscale+su]+=Vj×conv weight;16: end for 17: end for 18: end for 19: end for 20:end for 21:dirty image=IFFT(grid)22:return dirtyimage

3 并行算法

3.1 基于多核CPU的成像并行算法

混合w-facets成像算法的數(shù)據(jù)處理過程與傳統(tǒng)的射電干涉成像基本相同，大體包括4個(gè)部分：數(shù)據(jù)讀取，卷積核函數(shù)的創(chuàng)建，柵格化和逆傅里葉變換。其中，柵格化至少包含5層循環(huán)處理過程，相對其他過程，需要消耗更多的運(yùn)行時(shí)間。

因此，本文提出了一種粗粒度并行處理方法，實(shí)現(xiàn)每個(gè)分面數(shù)據(jù)的并行處理。采用開放式多處理(open multi-processing,OpenMP)并行技術(shù)，在最頂層加入并行任務(wù)調(diào)度，即可實(shí)現(xiàn)該循環(huán)的并行。本文主要采用靜態(tài)調(diào)度方式，具體并行過程：當(dāng)分配m個(gè)CPU核處理Nfacets個(gè)分面數(shù)據(jù)時(shí)，每個(gè)核(或線程)將處理Nfacets/m個(gè)分面的數(shù)據(jù)。當(dāng)分配的CPU核數(shù)與Nfacets相等時(shí)，并行效果最佳。這種并行方法避免了并行化過程中可能出現(xiàn)的同步和鎖的影響，是一種簡單實(shí)用的并行方法。

基于多核CPU的混合w-facets成像并行算法的實(shí)現(xiàn)見Algorithm 1偽代碼。算法的輸入數(shù)據(jù)為校準(zhǔn)后的數(shù)據(jù)，一般為寬帶數(shù)據(jù)，有多個(gè)不同頻率通道。最外的三層循環(huán)分別為：分面處理循環(huán)、按行處理循環(huán)和頻率通道處理循環(huán)。內(nèi)層的兩個(gè)循環(huán)，主要完成柵格化處理；最外層，是根據(jù)CPU核數(shù)來進(jìn)行分面處理的并行。最外三層循環(huán)，每一次將從輸入數(shù)據(jù)中讀取可見度、基線坐標(biāo)u-v-w和原始相位中心RA0和DEC0等主要數(shù)據(jù)；接著，根據(jù)分面數(shù)目Nfacets和原始相位中心計(jì)算出當(dāng)前分面的相位中心RAj和DECj。將公式(3)進(jìn)行相位旋轉(zhuǎn)，得到當(dāng)前分面的可見度數(shù)據(jù)Vj；根據(jù)式(6)進(jìn)行基線旋轉(zhuǎn)，獲得當(dāng)前分面的基線坐標(biāo)數(shù)據(jù)unew和vnew，并利用傅里葉變換縮放原理對旋轉(zhuǎn)后的基線數(shù)據(jù)進(jìn)行縮放。然后，通過w的最大值wmax和Nwplanes計(jì)算出當(dāng)前w在w核函數(shù)的索引位置windex。

最后，通過二層循環(huán)完成柵格化，每次循環(huán)將通過windex獲取w核函數(shù)的值，再與當(dāng)前可見度數(shù)據(jù)進(jìn)行乘累加。cfullsupport是w核中一維卷積核函數(shù)全支撐大小數(shù)組，通過windex獲取當(dāng)前w所需的全支撐大小。網(wǎng)格數(shù)組柵格的大小為完整圖像的大小，經(jīng)過逆傅里葉變換后將得到唯一臟圖。

3.2 基于GPU的成像并行算法

GPU是一種可以一次并行運(yùn)行數(shù)千個(gè)線程來實(shí)現(xiàn)快速計(jì)算的處理器，其計(jì)算能力是CPU的百倍。統(tǒng)一計(jì)算設(shè)備架構(gòu)(compute unified device architecture,CUDA)的內(nèi)部主要框架是以GPU中的線程塊為組進(jìn)行劃分，每塊由一個(gè)多核處理器控制，這意味著只需要控制多核處理器就可以同時(shí)控制和處理上萬個(gè)線程。本文將利用CUDA編程實(shí)現(xiàn)基于GPU的混合w-facets成像并行方法。該方法主要將輸入數(shù)據(jù)拷貝到GPU設(shè)備中,然后在GPU中進(jìn)行并行的柵格化處理。該方法讓主機(jī)負(fù)責(zé)整個(gè)系統(tǒng)的內(nèi)存管理、邏輯控制以及圖像形成，而讓GPU負(fù)責(zé)大量的并行計(jì)算，包括卷積核函數(shù)的計(jì)算，柵格化過程的處理等，分工明確。

在GPU中進(jìn)行柵格化處理時(shí)，使用原理為：同一基線連續(xù)測量的任意兩個(gè)數(shù)據(jù)之間的時(shí)間間隔足夠短，相同基線測量的數(shù)據(jù)集經(jīng)過柵格化處理后將落在相同的網(wǎng)格點(diǎn)上。同時(shí)借鑒Romein[18]的研究成果，將柵格化分成與卷積核函數(shù)相同大小的子網(wǎng)格，線程被分配用來處理每個(gè)子網(wǎng)格中的可見度數(shù)據(jù)，并依靠基線軌跡，將數(shù)據(jù)進(jìn)行累加并存儲到GPU寄存器中。當(dāng)基線軌跡移動到下一個(gè)網(wǎng)格點(diǎn)時(shí)，將當(dāng)前網(wǎng)格點(diǎn)內(nèi)的值進(jìn)行累加取均值，并在處理下一個(gè)網(wǎng)格點(diǎn)之前，將當(dāng)前網(wǎng)格點(diǎn)的處理結(jié)果從寄存器寫入全局內(nèi)存中[19]，這樣能減少GPU內(nèi)存的使用。

該方法利用基線進(jìn)行并行化，將同一基線不同時(shí)刻的觀測數(shù)據(jù)作為一組，每個(gè)線程負(fù)責(zé)處理一組數(shù)據(jù)。本文的加速方法所需的最小總線程數(shù)為：csmax×csmax×Nbaselines×Nfacets，其中csmax是cfullsupport數(shù)組中的最大值，Nbaselines為基線個(gè)數(shù)。此外，當(dāng)卷積核全支撐大小較小時(shí)，每個(gè)線程將處理1個(gè)分面中1個(gè)基線的1次柵格化循環(huán)(即1個(gè)卷積點(diǎn))，將最大限度地進(jìn)行并行化，循環(huán)次數(shù)減少到只有2層循環(huán)；當(dāng)卷積核全支撐大小較大時(shí)，每個(gè)線程將處理1個(gè)分面中1個(gè)基線的多個(gè)卷積點(diǎn)，因?yàn)槊總€(gè)線程塊的最大線程數(shù)目通常在256～1 024之間。相對文獻(xiàn)[12]，本文增加了分面的并行處理。為避免兩個(gè)或多個(gè)線程發(fā)生競爭現(xiàn)象，我們對全局內(nèi)存的寫入方式設(shè)置為原子性。

由于實(shí)際觀測的數(shù)據(jù)通常是按照時(shí)間順序記錄在表格文件里的，同一時(shí)刻包含不同基線的數(shù)據(jù)，為了進(jìn)行基線并行化，必須對輸入數(shù)據(jù)進(jìn)行重新排序。首先將輸入的數(shù)據(jù)按照不同基線分組并組成一排，然后對相同基線即一組內(nèi)的數(shù)據(jù)集按照觀測時(shí)間先后排序。對每組基線進(jìn)行編號，該編號是0至Nbaselines-1的自然數(shù)。天線在數(shù)據(jù)文件中有唯一編號，因此可以通過下式計(jì)算出基線的編號iB[12]：

其中，iant1和iant2分別是該基線中天線1和天線2的編號，Nant天線個(gè)數(shù)，imin是天線1和天線2的編號最小值。

由于輸入的數(shù)據(jù)一般已經(jīng)進(jìn)行過射頻干擾的去除和校準(zhǔn)處理，處理之后有些數(shù)據(jù)將會被刪除，每組基線內(nèi)的數(shù)據(jù)長度或時(shí)間步數(shù)(時(shí)間戳數(shù)目Nt)將會不同。因此，每組基線主要通過該組數(shù)據(jù)開始索引和時(shí)間戳數(shù)目循環(huán)(遍歷)處理完該組的數(shù)據(jù)[12]?；€的開始索引可以通過下式計(jì)算得到[12]：

其中，T是由所有組基線時(shí)間戳數(shù)目構(gòu)成的數(shù)組，數(shù)組的索引與每組基線的編號一致。

圖1給出基于GPU的混合w-facets成像并行算法流程。圖中主要分為主機(jī)處理、GPU處理和具體模塊處理三個(gè)部分來描述。在單個(gè)GPU節(jié)點(diǎn)中算法優(yōu)化的過程為：首先，主機(jī)讀取數(shù)據(jù)集；然后主機(jī)與GPU設(shè)備進(jìn)行通訊，完成GPU設(shè)備的初始化工作，包括選取性能最好的GPU卡和數(shù)據(jù)處理所需的內(nèi)存的分配等；接著，將存儲在主機(jī)內(nèi)存中的數(shù)據(jù)進(jìn)行重新排序，并將重排后的數(shù)據(jù)發(fā)送到已分配好的GPU內(nèi)存中。其次，在GPU中以基線并行的方式進(jìn)行數(shù)據(jù)處理，計(jì)算并行參數(shù)：根據(jù)線程索引計(jì)算出該線程所需要處理的基線起始索引iB、柵格化的卷積核全支撐索引(su和sv)和分面索引(j)。具體的計(jì)算方法，見式(12)。然后，利用計(jì)算好的索引，按照頻率通道(c)和時(shí)間戳數(shù)目(Nt)兩層循環(huán)，加載該線程所需要處理的數(shù)據(jù)，以及獲取w核函數(shù)，并進(jìn)行乘累加，從而完成基線的柵格化處理。循環(huán)結(jié)束后，即結(jié)束該線程的處理任務(wù)。最后，將柵格化后的數(shù)據(jù)返回到主機(jī)內(nèi)存中，并對該數(shù)據(jù)進(jìn)行逆傅里葉變換(IFFT)從而得到臟圖。

其中，csw是w值對應(yīng)的卷積核全支撐大小，thid是線程索引。

圖1 基于GPU的混合w-facets成像并行算法流程

4 實(shí)驗(yàn)與分析

4.1 算法驗(yàn)證

本節(jié)采用的數(shù)據(jù)來自SKA先導(dǎo)望遠(yuǎn)鏡MWA的銀河系與河外星系全天巡天(galactic and extragalactic all-sky MWA survey,GLEAM)項(xiàng)目中于2013年8月9日22:13:18開始并歷時(shí)112 s的快照觀測數(shù)據(jù)，該快照觀測的相位中心為(RA=03h13min54.40s,DEC=-55°04′02.82′′)，觀測頻率范圍是200～231 MHz。首先我們對該觀測數(shù)據(jù)進(jìn)行RFI消除、頻率與時(shí)間積分并轉(zhuǎn)換為Measurement Set格式，然后進(jìn)行校準(zhǔn)處理，處理完后獲得的數(shù)據(jù)將用于本文算法進(jìn)行成像。

使用的測試環(huán)境是上海天文臺SKA數(shù)據(jù)中心原型機(jī)的GPU節(jié)點(diǎn)，該節(jié)點(diǎn)具體的硬件指標(biāo)參數(shù)如表1所示。

表1 上海天文臺SKA數(shù)據(jù)中心原理樣機(jī)的GPU節(jié)點(diǎn)指標(biāo)參數(shù)

根據(jù)GLEAM巡天數(shù)據(jù)處理的經(jīng)驗(yàn)，我們設(shè)置需要計(jì)算的臟圖的大小為：4 000×4 000，每個(gè)像素的大小為25′′，因此臟圖的最大寬度Iwidth約為0.48 rad。該快照數(shù)據(jù)的wmax值約為1 035倍波長，根據(jù)式(9)我們計(jì)算出單獨(dú)使用w-projection時(shí)，最佳的Nwplanes約為248。如表2所示，是不同成像算法的測試結(jié)果對比。在成像質(zhì)量對比方面，我們主要給出的是均方根(RMS)、動態(tài)范圍(DR)和射電源識別數(shù)目。通過從這三方面的對比中發(fā)現(xiàn)，二維傅里葉變換方法(2D FFT，保持Nfacets=0和Nwplanes=0)成像質(zhì)量最差，w-stacking方法成像質(zhì)量最好。uv-faceting方法(保持Nwplanes=0)的成像質(zhì)量比2D FFT好，隨著Nfacets的增大成像質(zhì)量越好。w-projection方法(保持Nfacets=0)在Nwplanes≥31時(shí)，成像質(zhì)量均比2D FFT方法好；在較小Nwplanes時(shí)，w-projection比uv-faceting的成像質(zhì)量要好，成像質(zhì)量也隨著Nwplanes的增大質(zhì)量越高。對于w-facets方法的成像結(jié)果，在Nfacets=16和Nwplanes=62時(shí)質(zhì)量最高；在相同的取值時(shí)，比單獨(dú)使用uv-faceting或w-projection的質(zhì)量都要高，并且最接近w-stacking的成像質(zhì)量。此時(shí)，w-facets的圖像動態(tài)范圍比uv-faceting提高了2.34 dB，圖像RMS降低了4 mJy·beam-1。在Nfacets=4和Nwplanes=124時(shí)，w-facets的成像結(jié)果比相同取值的uv-faceting成像質(zhì)量好，比wprojection的差；在Nfacets=64和Nwplanes=31時(shí)，w-facets的成像結(jié)果比w-projection的都要好，與uv-faceting成像結(jié)果相近。成像質(zhì)量的分析表明，uv-faceting與w-projection結(jié)合后，成像質(zhì)量能夠得到進(jìn)一步提高。

在內(nèi)存使用方面，我們記錄了不同成像方法的卷積核內(nèi)存使用情況或傅里葉變換(FFT)內(nèi)存使用情況。w-stacking統(tǒng)計(jì)的是FFT內(nèi)存使用量，其余方法均是統(tǒng)計(jì)卷積核內(nèi)存使用情況，因?yàn)閣-stacking的w項(xiàng)校正在圖像域進(jìn)行，卷積核消耗的內(nèi)存小(小于1 MB)。由表2可知，2D FFT和uv-faceting的卷積核消耗的內(nèi)存均小于1 MB，w-projection方法的卷積核內(nèi)存消耗隨著Nwplanes的增大而增大，且遠(yuǎn)遠(yuǎn)大于1 MB；w-facets方法的卷積核內(nèi)存消耗隨著Nfacets的增大而減少，當(dāng)Nfacets=64時(shí)，內(nèi)存消耗小于1 MB；w-stacking方法由于使用了較大的w層，F(xiàn)FT內(nèi)存消耗較大，高達(dá)48 GB。分析表明，將w-projection與uv-faceting結(jié)合后，內(nèi)存的消耗大大減少。

對于運(yùn)行時(shí)間，我們給出的是使用20個(gè)CPU核的運(yùn)行時(shí)間和使用1個(gè)GPU卡的運(yùn)行時(shí)間，所有運(yùn)行時(shí)間均是測量10次取平均值。由于基于多核CPU的并行是針對分面，因此僅使用w-projection方法時(shí)，程序是串行執(zhí)行的，比uv-faceting方法消耗更多的時(shí)間。由于使用20個(gè)CPU核，uv-faceting方法在Nfacets=16時(shí)能夠充分并行，所以消耗的時(shí)間都比其他的分面數(shù)少。w-facets方法在Nfacets=16和Nwplanes=62時(shí)，耗時(shí)最短，與相同w平面數(shù)的w-projection方法耗時(shí)相近，但耗時(shí)比uv-faceting方法均長。上述分析表明，uv-faceting結(jié)合了w-projection之后，運(yùn)行時(shí)間會變長。使用了GPU加速后，大多數(shù)方法的運(yùn)行時(shí)間整體大大縮減了；uv-faceting和w-facets的運(yùn)行時(shí)間依然隨著Nfacets的增大而變長；w-projection方法隨著Nwplanes增大，運(yùn)行時(shí)間無較大變化，因?yàn)镚PU的加速方法主要是針對柵格化進(jìn)行并行，因此分面數(shù)越大需要的處理時(shí)間越多。w-facets方法中GPU的運(yùn)行速度比20個(gè)CPU核快7～33倍。

表2 測試結(jié)果對比

如圖2和圖3所示，分別為基于GPU的混合算法的成像結(jié)果和使用w-stacking算法的成像結(jié)果。為了便于比較，我們截取出距離原始圖像較遠(yuǎn)位置的射電源，并得到如圖2b)中和圖3b)中的兩個(gè)圖像。圖2b)圖像的峰值為每束光309.5 mJy，RMS值為每束光23.7 mJy。圖3b)圖像的峰值為每束光311.9 mJy，RMS值為每束光21.8 mJy。結(jié)合表2與整個(gè)圖像的對比，表明本文的GPU混合算法能夠?qū)嚯x較遠(yuǎn)的射電源成圖，且成圖結(jié)果與w-stacking獲得的結(jié)果基本一致。基于多核CPU的混合算法結(jié)果與GPU混合算法結(jié)果相同，這里不作展示。

圖2 基于GPU的混合w-facets成像并行算法的成像結(jié)果(Nfacets=16,Nwplanes=62)

4.2 性能測試與分析

性能測試部分采用模擬數(shù)據(jù)進(jìn)行，主要使用牛津大學(xué)開發(fā)的SKA射電望遠(yuǎn)鏡模擬軟件(OSKAR)仿真并生成SKA1-low觀測數(shù)據(jù)。主要的參數(shù)設(shè)置為：天線采用SKA1-low第四階段的完整規(guī)模天線陣列，臺站個(gè)數(shù)為512，基線個(gè)數(shù)為130 816，最長基線為65 km，觀測相位中心為(RA=201.36,DEC=43.02)，偏振個(gè)數(shù)為4，波束個(gè)數(shù)為1，觀測的起始頻率為50 MHz，頻率帶寬為300 MHz，頻率通道數(shù)為64，每個(gè)頻率通道的帶寬為4 687.5 kHz；積分時(shí)間為2 min，起始觀測時(shí)間設(shè)置為2015年01月01日18點(diǎn)整，總觀測時(shí)長為2 min。觀測相位中心指向的天區(qū)包含一個(gè)巨大的射電星系NGC 5128，也稱為半人馬座A(Centaurus A)。它不僅是距離地球最近的射電星系(距離大約為11.09光年)，而且是南半球最大的星系。本模擬觀測使用的天空模型是由包含Centaurus A的MWA GLEAM巡天圖像制作得到[20]。該模型共有23 811個(gè)獨(dú)立成分，包含23 336個(gè)高斯源和2 475個(gè)點(diǎn)源。所有成分的流量強(qiáng)度和譜指數(shù)分布在74～231 MHz頻率范圍內(nèi)。最終模擬觀測產(chǎn)生的數(shù)據(jù)量大小約為36 GB。該數(shù)據(jù)曾經(jīng)作為SKA科學(xué)數(shù)據(jù)處理器工作包中基于執(zhí)行框架的大視場成像流水線的性能測試[21]，因此比較適合用于本文的性能測試。

實(shí)驗(yàn)設(shè)置的成像參數(shù)是：臟圖大小為1 024×1 024，每個(gè)像素大小為14.06′，Nfacets=36。該模擬數(shù)據(jù)的wmax=61 026倍波長，根據(jù)臟圖的大小和分面數(shù)，我們計(jì)算出最佳的Nwpalanes=355。該部分的實(shí)驗(yàn)采用的測試環(huán)境與驗(yàn)證性實(shí)驗(yàn)相同。如圖4所示，是性能測試的結(jié)果，主要使用了不同CPU核數(shù)(線程數(shù))測試基于多核CPU方法的可拓展性，使用一個(gè)GPU卡測試基于GPU方法的加速性能。圖中每個(gè)測試點(diǎn)的運(yùn)行時(shí)間均是連續(xù)測試10次取平均值。我們根據(jù)圖4計(jì)算出加速比情況，如表3所示。根據(jù)多核CPU的方法，當(dāng)使用的線程數(shù)為36時(shí)，線程數(shù)與Nfacets值相等，因此消耗的運(yùn)行時(shí)間最短；使用1個(gè)線程時(shí)，相當(dāng)于串行程序，運(yùn)行時(shí)間最長；線程數(shù)為1～12，消耗的運(yùn)行時(shí)間近似線性減少；線程數(shù)大于12后，沒有明顯的線性關(guān)系，在線程數(shù)為24和36時(shí)，運(yùn)行時(shí)間出現(xiàn)下降。使用V100 GPU時(shí)，單精度加速比為201.8，運(yùn)行速度是36個(gè)線程數(shù)的8.9倍左右；雙精度加速比為173.5，運(yùn)行速度是36個(gè)線程數(shù)的7.7倍左右。上述分析表明，基于多核CPU的方法在少量線程數(shù)內(nèi)有一定的可拓展性，當(dāng)線程數(shù)與Nfacets相等時(shí)，運(yùn)行速度最快；基于GPU的方法在運(yùn)行速度上比多核CPU方法有明顯的優(yōu)勢。

圖4 混合w-facets成像并行算法性能測試結(jié)果

表3 加速比情況

5 總結(jié)與展望

由于低頻射電望遠(yuǎn)鏡陣列對成像具有高動態(tài)范圍和高分辨率的要求，因此面對海量天文觀測數(shù)據(jù)時(shí)，既須保證大視場成像的成圖質(zhì)量也必須具有較快的運(yùn)行速度。故本文研究了uv-faceting與w-projection相混合的算法，稱為w-facets算法，為將來開展低頻射電望遠(yuǎn)鏡陣列巡天觀測的數(shù)據(jù)處理提供技術(shù)手段。與uv-faceting算法相比，w-facets算法能夠進(jìn)一步減少噪聲水平，提高成圖的動態(tài)范圍，改善成像質(zhì)量；與w-projection算法相比，w-facets算法大大減少了內(nèi)存的消耗，提高了可拓展性。w-facets在uv-faceting的基礎(chǔ)上使用w核進(jìn)行柵格化，運(yùn)行時(shí)間比uv-faceting慢很多，但與w-projection的運(yùn)行時(shí)間相近。為此，本文提出了兩種并行優(yōu)化方法：基于多核CPU的混合w-facets成像并行方法和基于GPU的混合w-facets成像并行方法。通過圖像驗(yàn)證性實(shí)驗(yàn)和性能測試實(shí)驗(yàn)表明，本文提出的兩種方法均能夠正確實(shí)現(xiàn)混合w-facets算法成像，并提升混合算法的成像速度?；诙嗪薈PU的并行加速方法運(yùn)行速度最大提升約22.7倍，基于GPU的并行加速方法運(yùn)行速度最大可以提升約201.8倍；這表明基于GPU的并行加速方法比基于多核CPU的并行加速方法更高效。但是通過本文測試分析可知，基于多核CPU的并行加速方法效果略差，因?yàn)椴捎玫氖谴至６鹊牟⑿蟹椒ǎ罄m(xù)工作將考慮通過更細(xì)粒度的CPU并行方案進(jìn)行優(yōu)化。