混合w-facets成像并行算法研究

于昂，勞保強，王俊義，安濤

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

1 引言

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

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

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

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

2 混合w-facets成像算法

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

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

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

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

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

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

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

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

式(3)末尾的2個平方根項，經過泰勒一級近似展開，有：

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

其中，uj和vj是第j個小視場的基線坐標，計算公式為：

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

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

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

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

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

Algorithm 1基于多核CPU的成像并行算法Input:V可見度數據，(u,v,w)基線坐標，Nwplanes w平面數，網格數組柵格存放柵格化結果全零數組(完整圖像大小)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成像算法的數據處理過程與傳統(tǒng)的射電干涉成像基本相同，大體包括4個部分：數據讀取，卷積核函數的創(chuàng)建，柵格化和逆傅里葉變換。其中，柵格化至少包含5層循環(huán)處理過程，相對其他過程，需要消耗更多的運行時間。

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

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

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

3.2 基于GPU的成像并行算法

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

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

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

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

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

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

其中，T是由所有組基線時間戳數目構成的數組，數組的索引與每組基線的編號一致。

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

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

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

4 實驗與分析

4.1 算法驗證

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

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

表1 上海天文臺SKA數據中心原理樣機的GPU節(jié)點指標參數

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

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

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

表2 測試結果對比

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

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

4.2 性能測試與分析

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

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

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

表3 加速比情況

5 總結與展望

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