基于RASCIL的W-projection和W-stacking并行算法實測研究

楊秋萍，朵 琳

(1. 中國科學院云南天文臺，云南 昆明 650216；2. 中國科學院大學，北京 100049；3. 昆明理工大學信息工程與自動化學院，云南 昆明 650500)

在射電干涉陣的成像研究中，隨著觀測視場的擴大，射電干涉陣成像時需要進行大視場成像。射電干涉陣的大視場成像方法包括三維傅里葉變換的方法[1-2]，F(xiàn)aceting方法[3]，W-projection方法[4]和W-stacking方法等[5-6]。

平方公里陣列(Square Kilometre Array, SKA)為射電干涉陣成像帶來了新的挑戰(zhàn)。平方公里陣列的巨大規(guī)模和復雜程度遠遠超出了現(xiàn)有射電天文望遠鏡陣列，全規(guī)模運行產(chǎn)生的海量數(shù)據(jù)需要10億億次/秒的處理能力，是2017年最快的超級計算機神威太湖之光處理能力(0.9億億次/秒)的10倍多[7]。因此，在射電干涉陣的數(shù)據(jù)處理中，多節(jié)點的并行化實現(xiàn)一直是研究的重要內(nèi)容。

在射電干涉陣的成像中，大視場成像方法的并行實現(xiàn)已成為研究的重點。文[8]實現(xiàn)了W-projection的中央處理器(Central Processing Unit, CPU)并行和圖形處理器(Graphics Processing Unit, GPU)并行，并在天河二號(MilkyWay-2)超級計算機上進行了實驗，比較了單精度和雙精度情況下數(shù)據(jù)加載時間和網(wǎng)格化運行時間。文[9]對uv域的Faceting成像方法并行化進行了研究。文[10]在uv域的Faceting成像方法和W-projection成像方法的基礎上提出了一種新的方法w-facets，并通過多核中央處理器和圖形處理器進行了并行實現(xiàn)。文[11]提出新的網(wǎng)格化核函數(shù)，并通過并行實現(xiàn)測試核函數(shù)在網(wǎng)格化時的性能。文[12]對提出的網(wǎng)格化方法進行了并行測試，給出了Strong Scaling下測試結(jié)果。

本文主要基于RASCIL對W-projection和W-stacking方法利用DASK進行中央處理器并行實現(xiàn)，分析在RASCIL中兩種成像方法并行時的策略選擇，以及并行資源的消耗，為后續(xù)基于RASCIL的射電干涉陣數(shù)據(jù)處理提供參考。

1 W-projection和W-stacking的實現(xiàn)框架

W-projection和W-stacking是射電干涉陣大視場成像的經(jīng)典算法，W-projection是CASA(https://casadocs.readthedocs.io/en/stable/)處理大視場成像的方法之一，而W-stacking是WSClean (https://gitlab.com/aroffringa/wsclean/)處理大視場成像的方法之一。

圖1是W-projection成像方法實現(xiàn)的原理框圖，圖中最左側(cè)是射電干涉陣接收的可見度數(shù)據(jù)集，在成像前可見度數(shù)據(jù)集已完成校準，數(shù)據(jù)集加載后，可以獲得觀測的(u,v,w)分布，利用這些信息通過

(1)

計算得到成像時所需要的w平面數(shù)Nw，其中，δA=0.02[9]；FOV為成像的視場；wstep為相鄰w平面的間隔；wmax為w的最大值。然后利用卷積核函數(shù)生成Nw個卷積核，這也是W-projection成像方法的關鍵，具有不同w的可見度數(shù)據(jù)與對應的卷積核進行卷積，實現(xiàn)從(u,v,w)向(u,v)平面的投影，卷積完成的可見度數(shù)據(jù)分布在(u,v)平面，并且進行疊加，最后通過傅里葉逆變換得到成像的臟圖。

圖1 W-projection實現(xiàn)的原理框圖Fig.1 Principles of W-projection implementation

圖2 W-stacking實現(xiàn)的原理框圖Fig.2 Principles of W-stacking implementation

我們通過分析兩種成像方法的實現(xiàn)框圖發(fā)現(xiàn)， W-Projectin方法與W-stacking方法相同之處是都需要根據(jù)觀測校準的可見度數(shù)據(jù)(u,v,w)的信息計算w平面數(shù)Nw。兩種成像方法不同之處是W-projection方法需要根據(jù)w平面數(shù)計算投影的卷積核，卷積核數(shù)量與Nw一致，在計算時需要考慮卷積核的寬度和過采樣點的數(shù)量，這兩個因素影響卷積核的計算時間和占用空間；而W-stacking需要根據(jù)w平面數(shù)對可見度數(shù)據(jù)進行劃分，將對應w值的可見度數(shù)據(jù)劃分到同一個數(shù)據(jù)切片，數(shù)據(jù)切片集合的數(shù)量與w平面數(shù)Nw一致。W-stacking在卷積網(wǎng)格化環(huán)節(jié)中只需要一個卷積核，也就是所有數(shù)據(jù)切片會與同樣的卷積核進行卷積網(wǎng)格化，網(wǎng)格化后再通過傅里葉逆變換得到數(shù)據(jù)切片對應的圖像。

在成像處理過程中，W-projection需要計算與Nw數(shù)量相等的卷積核，因此這一部分的計算消耗超過W-stacking，但是由于W-stacking方法需要首先將可見度數(shù)據(jù)集進行劃分數(shù)據(jù)切片，所以當可見度數(shù)據(jù)的規(guī)模達到一定程度時，需要特別考慮這一操作的計算效率；又因為W-stacking在實現(xiàn)時需要存儲與w平面數(shù)相等個數(shù)的數(shù)據(jù)切片，因此，當可見度數(shù)據(jù)集具有一定規(guī)模時，需要的內(nèi)存規(guī)模遠遠超過W-projection。

2 并行實驗

2.1 實驗數(shù)據(jù)

實驗數(shù)據(jù)來自Karl G. Janskey Very Large Array陣列的D陣型觀測校準數(shù)據(jù)SNR_G55_10s.calib.ms，觀測目標是超新星遺跡G055.7+3.4，觀測的相位中心為19: 21: 40，DEC為21.45.00，觀測日期從2010-08-23-01: 07: 14.00(UTC)至2010-08-23-08: 14: 54.00(UTC)，觀測頻段是L頻段1～2 GHz，包含4個頻譜窗口，每個頻譜窗口有64個頻道，每個頻道2種極化方式，數(shù)據(jù)大小為1.4 G。根據(jù)觀測陣列的u,v,w分布，計算最長基線和最大w值，再根據(jù)成像時每個像素大小為8″，圖像每個方向包含1 280個像素，可以根據(jù)(1)式計算得到在上述成像參數(shù)下，改正w-term的影響所必需的w平面數(shù)為68。

2.2 并行計算

(1)并行環(huán)境

本文使用RASCIL進行成像。RASCIL是一個完全開源的射電干涉陣列數(shù)據(jù)處理包，已廣泛用于平方公里陣列數(shù)據(jù)模擬和處理研究，處理結(jié)果與其他軟件進行比較后證實是可靠的。RASCIL用Python開發(fā)，在實驗中使用的版本為v.0.1.9。

并行處理基于DASK進行，DASK也是一個開源庫，為現(xiàn)有的Python堆棧提供并行性，DASK與Python庫集成。在科學計算中數(shù)據(jù)集和計算規(guī)模的擴展速度遠遠高于處理器和內(nèi)存發(fā)展的速度，科學計算往往需要擴展到多臺計算機進行，DASK提供了可以跨多個核心、處理器和計算機實現(xiàn)并行執(zhí)行。

(2)實驗硬件配置

成像實驗使用相同配置的硬件作為計算節(jié)點，硬件配置如表1。各個節(jié)點通過光纖萬兆交換機和路由器進行連接組成以太網(wǎng)，保證各計算節(jié)點在同一局域網(wǎng)組網(wǎng)，每個節(jié)點有相同的網(wǎng)絡延遲和帶寬。

表1 并行計算節(jié)點的硬件配置Table 1 Hardware configuration of parallel computing nodes

(3)實驗軟件配置

操作系統(tǒng)為Ubuntu 22.04 LTS，開發(fā)語言環(huán)境為Python v3.8.2，并行計算框架為DASK2022.6.1，應用軟件為RASCIL v0.1.9。

2.3 兩種成像方法的并行策略

兩種大視場成像方法的并行實現(xiàn)中，當處理器負載最大，或者數(shù)據(jù)流帶寬達到上限時，計算系統(tǒng)的吞吐量達到上限飽和，計算系統(tǒng)遇到瓶頸，無法提高處理速度。因此，成像算法并行化設計的前提是在保證計算系統(tǒng)吞吐量最大時對計算性能結(jié)果進行分析，即在不同算法和數(shù)據(jù)切片并行調(diào)度方式下，所分配的處理器資源或者帶寬資源應接近100%占用率。

兩種大視場成像算法的實現(xiàn)核心可以分為3部分：(1)數(shù)據(jù)集加載和數(shù)據(jù)集預處理；(2)計算任務和數(shù)據(jù)的分布式調(diào)度執(zhí)行；(3)計算結(jié)果聚合，為下一環(huán)節(jié)管線做準備。

2.3.1 W-projection的并行策略

W-projection能夠分解成一些完全獨立的子任務，同時各個子任務之間的數(shù)據(jù)幾乎沒有依賴，沒有通信。以每通道可見度數(shù)據(jù)分配并行計算任務，數(shù)據(jù)由vis_load_ms task讀入，創(chuàng)建各通道圖像和卷積核(create_wp_gcfcf_from_vis)，完成網(wǎng)格化相關invert和sum_invert task計算，并將各通道圖像通過gather_image將圖像結(jié)果匯集，如圖3，在圖中需要持久化的計算環(huán)節(jié)用 “P” 標注。

2.3.2 W-stacking的并行策略

W-stacking以每通道可見度數(shù)據(jù)并行分配計算任務的方式，數(shù)據(jù)由load_ms task讀入，每個通道做數(shù)據(jù)w_slice切片后進行Scatter-Gather并行執(zhí)行，分片(Scatter)之后每個分片進行gridding_invert task計算，計算結(jié)果進行聚合(Gather)，并將各通道圖像結(jié)果匯集gather_image。 W-stacking按照w進行數(shù)據(jù)分片，有兩種并行策略：

(1)多通道并行，每通道分片并行計算，這樣會消耗非常多的內(nèi)存資源， 任務并行關系如圖4(a)，其中持久化的計算結(jié)果環(huán)節(jié)用 “P” 表示。由于RASCIL中數(shù)據(jù)分片數(shù)量和任務數(shù)量相等，導致在并行處理時系統(tǒng)內(nèi)存開銷過大，因此，在目前的硬件配置下，這種策略無法得到成功的計算結(jié)果。

(2)多通道并行，每通道順序分片并行計算，節(jié)省內(nèi)存資源，任務并行關系如圖4(b)，其中持久化的計算結(jié)果環(huán)節(jié)用 “P” 表示，因此，在W-stacking的并行實現(xiàn)時選擇這種策略。

圖3 W-projection的并行策略Fig.3 Parallel computing implementation of W-projection

在圖3和圖4給出了4條數(shù)據(jù)通道的任務圖，在實驗中可見度數(shù)據(jù)的通道數(shù)最高為64。

3 結(jié)果討論

在成像方法的DASK并行計算實驗中，為減少進程切換損耗和并行調(diào)度復雜度，每個執(zhí)行單元(worker)對應一個進程和線程，執(zhí)行單元總數(shù)量不超過處理器邏輯核心數(shù)量，每個執(zhí)行單元最大內(nèi)存為12 GB。

在并行處理中最小參考基準是單一執(zhí)行單元，單一物理核心，單一進程和線程，對于超線程架構(gòu)和睿頻技術的中央處理器核心，邏輯核心計算性能不等同于物理核心計算性能。因此，在并行計算時按照同一計算場景，同一節(jié)點上56邏輯核心和28物理核心計算時間(分別為208.49 s和190.94 s)進行如下?lián)Q算：

圖4 W-stacking的并行策略Fig.4 Parallel computing implementation of W-stacking

3.1 并行計算時間

(1)在數(shù)據(jù)處理與計算規(guī)模恒定的情況下，即這時處理的可見度數(shù)據(jù)均為64通道， 對應不同計算集群和中央處理器核心數(shù)量(每個執(zhí)行單元對應一個物理或邏輯核心)配置下， 對兩種成像方法并行計算完成的時間進行比較，如表2和表3，可以看出在增加節(jié)點或節(jié)點的中央處理器核心數(shù)量時， 兩種成像方法的計算時間都減少。

表2 W-projection并行計算時間(單位： s)

表3 W-stacking并行計算時間 (單位： s)

(2)在集群節(jié)點數(shù)量恒定為2，中央處理器物理核心數(shù)量隨著處理數(shù)據(jù)規(guī)模等比例增加的情況下，并行計算完成的時間如表4，由表4可以看到，當數(shù)據(jù)規(guī)模不斷增加時，隨著物理核心數(shù)的增加，計算時間也在增加。

3.2 成像方法并行效率分析

在并行處理中，加速比定義為

(2)

其中，Sp為加速比；Ep為并行效率；Ts為最優(yōu)順序算法的單處理核心執(zhí)行時間；Tp為使用p個處理器核心并行計算所花費的時間。

并行加速比和兩種成像方法并行效率如圖5，其中(a)為兩種成像方法的并行加速比，(b)為并行效率。由于并行任務的分配不一定和執(zhí)行單元數(shù)量整除，計算時間和系統(tǒng)并行效率并不隨運行中央處理器核心數(shù)量線性增加或減少，所以圖中出現(xiàn)核心數(shù)多反而更慢的起伏點。處理密集運算時，并行效率衰減較快，這是因為一方面并行調(diào)度開銷的增長，另一方面現(xiàn)代中央處理器睿頻和超線程與多中央處理器的架構(gòu)特性為少量核心運算的計算場景提供更好的核心性能。

表4 數(shù)據(jù)規(guī)模和計算核心同步增長時并行計算時間(單位： s)

圖5 并行加速比和并行效率Fig.5 Parallel speedup and system efficiency

3.3 數(shù)據(jù)規(guī)模分析

當處理的可見度數(shù)據(jù)規(guī)模發(fā)生變化時，針對兩種并行計算資源配置的不同場景進行分析。

(1)Strong Scaling：并行處理時，問題規(guī)模保持不變(即可見度數(shù)據(jù)的通道數(shù)量不變)，增加處理器數(shù)量，用于找到解決該問題最合適的處理器數(shù)量，即所用時間盡可能短又不產(chǎn)生太大的開銷。兩種成像方法在Strong Scaling資源配置模式的并行計算時間與加速比如圖6，(a)為W-projection，(b)為W-stacking。從圖6可以看出，在處理的可見度數(shù)據(jù)通道數(shù)量不變時，W-projection在43核心時達到計算速度與每核心使用效率的最佳平衡點，而W-stacking在28核時達到。

圖6 Strong Scaling資源配置Fig.6 Strong Scaling parallel computing configuration

(2)Weak Scaling：并行處理問題規(guī)模隨處理器數(shù)量增加而增加，適合這種并行資源配置策略時并行效率會保持水平穩(wěn)定。測試結(jié)果如圖7，從圖7可以看出，并行效率迅速下降，因此,兩種成像方法的并行實現(xiàn)在數(shù)據(jù)規(guī)模發(fā)生變化時并不適合Weak Scaling并行資源配置策略。

圖7 Weak Scaling資源配置Fig.7 Weak Scaling parallel computing configuration

4 結(jié) 論

本文首先介紹了W-projection和W-stacking兩種大視場成像方法的實現(xiàn)原理，在此基礎上完成基于RASCIL的W-projection和W-stacking方法的DASK并行實驗，采用VLA-D陣列對超新星遺跡G055.7+3.4的校準觀測數(shù)據(jù)集進行兩種成像方法的并行處理，對兩種成像方法設計了并行計算策略，在不同集群節(jié)點和處理器核心數(shù)配置下，分別在數(shù)據(jù)規(guī)模恒定和數(shù)據(jù)規(guī)模隨處理器增長兩種場景做了并行計算時間統(tǒng)計，對加速比和并行實現(xiàn)效率進行分析，并在數(shù)據(jù)規(guī)模增長時對比了兩種成像方法在Strong Scaling和Weak Scaling兩種策略下并行計算的優(yōu)劣，兩種成像方法適合采用Strong Scaling并行資源配置方式。由于在RASCIL v0.1.9版本中，可見度數(shù)據(jù)是利用visibility對象進行存儲和訪問的，在進行并行計算時該對象變得非常龐大，而W-stacking方法需要多次對visibility對象進行讀寫，所以該方法在并行計算時消耗的時間超過W-projection。在未來基于RASCIL的成像并行計算中，可以考慮對可見度數(shù)據(jù)采用更加簡單的數(shù)據(jù)結(jié)構(gòu)，這樣能大幅提升W-stacking的成像效率。